ElectricCalculation.cpp

/*
 *  Copyright (C) 2017  Thales Lima Oliveira <thales@ufu.br>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#include "ElectricCalculation.h"
#ifdef USING_WX_3_0_X
#include "../utils/DegreesAndRadians.h"
#endif
 
#include <wx/clipbrd.h>
 
ElectricCalculation::ElectricCalculation() {}
ElectricCalculation::~ElectricCalculation() {}
void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
{
    m_powerElementList.clear();
    m_busList.clear();
    m_capacitorList.clear();
    m_indMotorList.clear();
    m_inductorList.clear();
    m_lineList.clear();
    m_loadList.clear();
    m_syncGeneratorList.clear();
    m_syncMotorList.clear();
    m_transformerList.clear();
    m_harmCurrentList.clear();
    m_emtElementList.clear();
    // TODO: Bad design?
    for (auto it = elementList.begin(); it != elementList.end(); it++) {
        Element* element = *it;
        m_powerElementList.push_back(static_cast<PowerElement*>(element));
 
        if (Bus* bus = dynamic_cast<Bus*>(element))
            m_busList.push_back(bus);
        else if (Capacitor* capacitor = dynamic_cast<Capacitor*>(element))
            m_capacitorList.push_back(capacitor);
        else if (IndMotor* indMotor = dynamic_cast<IndMotor*>(element))
            m_indMotorList.push_back(indMotor);
        else if (Inductor* inductor = dynamic_cast<Inductor*>(element))
            m_inductorList.push_back(inductor);
        else if (Line* line = dynamic_cast<Line*>(element))
            m_lineList.push_back(line);
        else if (Load* load = dynamic_cast<Load*>(element))
            m_loadList.push_back(load);
        else if (SyncGenerator* syncGenerator = dynamic_cast<SyncGenerator*>(element))
            m_syncGeneratorList.push_back(syncGenerator);
        else if (SyncMotor* syncMotor = dynamic_cast<SyncMotor*>(element))
            m_syncMotorList.push_back(syncMotor);
        else if (Transformer* transformer = dynamic_cast<Transformer*>(element))
            m_transformerList.push_back(transformer);
        else if (HarmCurrent* harmCurrent = dynamic_cast<HarmCurrent*>(element))
            m_harmCurrentList.push_back(harmCurrent);
        else if (EMTElement* emtElement = dynamic_cast<EMTElement*>(element))
            m_emtElementList.push_back(emtElement);
    }
 
    // Set buses numbers
    int busNumber = 0;
    for (auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
        Bus* bus = *itb;
        BusElectricalData data = bus->GetElectricalData();
        data.number = busNumber;
        bus->SetElectricalData(data);
        busNumber++;
    }
}
 
bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus,
    double systemPowerBase,
    YBusSequence sequence,
    bool includeSyncMachines,
    bool allLoadsAsImpedances,
    bool usePowerFlowVoltagesOnImpedances)
{
    if (m_busList.size() == 0) return false;
 
    // Clear and fill with zeros the Ybus
    yBus.clear();
    for (int i = 0; i < (int)m_busList.size(); i++) {
        std::vector<std::complex<double> > line;
        for (int j = 0; j < (int)m_busList.size(); j++) { line.push_back(std::complex<double>(0.0, 0.0)); }
        yBus.push_back(line);
    }
 
    // Set buses numbers and reset connection status
    int busNumber = 0;
    for (auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
        Bus* bus = *itb;
        BusElectricalData data = bus->GetElectricalData();
        data.number = busNumber;
        data.isConnected = true;
        //data.voltage = std::complex<double>(1.0, 0.0);
        bus->SetElectricalData(data);
        busNumber++;
    }
 
    // Load
    for (auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) {
        Load* load = *it;
        if (load->IsOnline()) {
            int n = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().number;
            LoadElectricalData data = load->GetPUElectricalData(systemPowerBase);
            if (data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) {
                std::complex<double> yLoad = std::complex<double>(data.activePower, -data.reactivePower);
                if (allLoadsAsImpedances && usePowerFlowVoltagesOnImpedances) {
                    std::complex<double> v = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().voltage;
                    yLoad /= (std::abs(v) * std::abs(v));
                }
                yBus[n][n] += yLoad;
            }
        }
    }
 
    // Capacitor
    for (auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) {
        Capacitor* capacitor = *it;
        if (capacitor->IsOnline()) {
            int n = static_cast<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number;
            CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase);
            yBus[n][n] += std::complex<double>(0.0, data.reactivePower);
        }
    }
 
    // Inductor
    for (auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) {
        Inductor* inductor = *it;
        if (inductor->IsOnline()) {
            int n = static_cast<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number;
            InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase);
            yBus[n][n] += std::complex<double>(0.0, -data.reactivePower);
        }
    }
 
    //EMT Elements
    for (EMTElement* emtElement : m_emtElementList) {
        if (emtElement->IsOnline()) {
            auto data = emtElement->GetEMTElementData();
            int n = static_cast<Bus*>(emtElement->GetParentList()[0])->GetElectricalData().number;
            yBus[n][n] += data.y0;
            //wxMessageBox(wxString::Format("%f +j %f", data.y0.real(), data.y0.imag()));
        }
    }
 
    // Power line
    for (auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) {
        Line* line = *it;
        if (line->IsOnline()) {
            LineElectricalData data = line->GetPUElectricalData(systemPowerBase);
 
            int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number;
            int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number;
 
            switch (sequence) {
            case POSITIVE_SEQ:
            case NEGATIVE_SEQ: {
                yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
                yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
                yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
                yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
                yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0);
                yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0);
            } break;
            case ZERO_SEQ: {
                yBus[n1][n2] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
                yBus[n2][n1] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
 
                yBus[n1][n1] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
                yBus[n2][n2] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
 
                yBus[n1][n1] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
                yBus[n2][n2] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
            }
            }
        }
    }
 
    // Transformer
    for (auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) {
        Transformer* transformer = *it;
 
        if (transformer->IsOnline()) {
            TransformerElectricalData data = transformer->GetPUElectricalData(systemPowerBase);
 
            int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number;
            int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().number;
 
            // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as series
            // impedance.
            if (data.turnsRatio == 1.0 && data.phaseShift == 0.0 && sequence != ZERO_SEQ) {
                yBus[n1][n2] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
                yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
 
                yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
                yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
            }
            // If the transformer have no-nominal turn ratio and/or phase shifting, it will be modelled in a different
            // way (see references).
            //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232]
            //[Ref. 2: http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf]
            //[Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf]
            else if (sequence != ZERO_SEQ) {
                // Complex turns ratio
                double radPhaseShift = wxDegToRad(data.phaseShift);
                std::complex<double> a = std::complex<double>(data.turnsRatio * std::cos(radPhaseShift),
                    -data.turnsRatio * std::sin(radPhaseShift));
 
                // Transformer admitance
                std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
                if (sequence == POSITIVE_SEQ) {
                    yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
                    yBus[n1][n2] += -(y / std::conj(a));
                    yBus[n2][n1] += -(y / a);
                    yBus[n2][n2] += y;
                }
                else if (sequence == NEGATIVE_SEQ) {
                    yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
                    yBus[n1][n2] += -(y / a);
                    yBus[n2][n1] += -(y / std::conj(a));
                    yBus[n2][n2] += y;
                }
            }
            else if (sequence == ZERO_SEQ) {
                switch (data.connection) {
                case GWYE_GWYE: {
                    std::complex<double> y =
                        1.0 /
                        std::complex<double>(
                            data.zeroResistance + 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance),
                            data.zeroIndReactance +
                            3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance));
                    std::complex<double> a = std::complex<double>(data.turnsRatio, 0.0);
 
                    yBus[n1][n1] += y / (a * a);
                    yBus[n1][n2] += -(y / a);
                    yBus[n2][n1] += -(y / a);
                    yBus[n2][n2] += y;
                } break;
                case DELTA_GWYE: {
                    std::complex<double> y =
                        1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.secondaryGrndResistance),
                            data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance));
                    yBus[n2][n2] += y;
                    yBus[n1][n1] += 1e-4;
                    yBus[n1][n2] += 1e-4;
                    yBus[n2][n1] += 1e-4;
                    break;
                }
                case GWYE_DELTA: {
                    std::complex<double> y =
                        1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.primaryGrndResistance),
                            data.zeroIndReactance + 3.0 * (data.primaryGrndReactance));
                    yBus[n1][n1] += y;
                    yBus[n2][n2] += 1e-4;
                    yBus[n1][n2] += 1e-4;
                    yBus[n2][n1] += 1e-4;
                    break;
                }
                default:
                    break;
                }
            }
        }
    }
 
    if (includeSyncMachines) {
        // Synchronous generator
        for (auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
            SyncGenerator* syncGenerator = *it;
            if (syncGenerator->IsOnline()) {
                int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
                SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase);
                switch (sequence) {
                case POSITIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
                } break;
                case NEGATIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
                } break;
                case ZERO_SEQ: {
                    if (data.groundNeutral) {
                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance + 3.0 * data.groundResistance,
                            data.zeroReactance + 3.0 * data.groundReactance);
                    }
                } break;
                }
            }
        }
        // Synchronous motor
        for (auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) {
            SyncMotor* syncMotor = *it;
            if (syncMotor->IsOnline()) {
                int n = static_cast<Bus*>(syncMotor->GetParentList()[0])->GetElectricalData().number;
                SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase);
                switch (sequence) {
                case POSITIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
                } break;
                case NEGATIVE_SEQ: {
                    yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
                } break;
                case ZERO_SEQ: {
                    if (data.groundNeutral) {
                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance + 3.0 * data.groundResistance,
                            data.zeroReactance + 3.0 * data.groundReactance);
                    }
                } break;
                }
            }
        }
    }
 
    // Identify buses to slack bus (power flow)
    //if (sequence != ZERO_SEQ) {
    std::vector<bool> connectedToSlack(m_busList.size(), false);
    busNumber = 0;
    int slackBus = 0;
    std::vector<int> checkBusVector;
    for (auto* bus : m_busList)
    {
        if (bus->GetElectricalData().slackBus)
        {
            connectedToSlack[busNumber] = true;
            slackBus = busNumber;
            break;
        }
        busNumber++;
    }
    GetNextConnection(slackBus, yBus, connectedToSlack);
 
    bool hasIsolatedBus = false;
    for (auto isConnectedToSlack : connectedToSlack)
    {
        if (!isConnectedToSlack) hasIsolatedBus = true;
    }
 
    if (hasIsolatedBus) {
        // Update isolated and connected buses status and numbers
        busNumber = 0;
        for (unsigned int i = 0; i < m_busList.size(); ++i)
        {
            auto data = m_busList[i]->GetElectricalData();
            data.isConnected = connectedToSlack[i];
            if (connectedToSlack[i])
            {
                data.number = busNumber;
                busNumber++;
            }
            else {
                data.number = -1;
            }
            m_busList[i]->SetElectricalData(data);
        }
 
        // Create new Ybus without isolated buses (remove rows and columns)
        std::vector< std::vector< std::complex<double> > > newYBus;
        for (unsigned int i = 0; i < yBus.size(); ++i) {
            if (connectedToSlack[i]) {
                std::vector< std::complex<double> > newLine;
                for (unsigned int j = 0; j < yBus.size(); ++j) {
                    if (connectedToSlack[j]) {
                        newLine.push_back(yBus[i][j]);
                    }
                }
                newYBus.push_back(newLine);
            }
        }
 
        yBus.clear();
        yBus = newYBus;
    }
    //}
 
    return true;
}
 
void ElectricCalculation::GetNextConnection(const unsigned int& checkBusNumber,
    const std::vector< std::vector< std::complex<double> > >& yBus,
    std::vector<bool>& connToSlack)
{
    for (unsigned int i = 0; i < yBus.size(); i++) {
 
        if (i != checkBusNumber && !connToSlack[i]) {
            if (std::abs(yBus[i][checkBusNumber]) > 1e-6)
            {
                connToSlack[i] = true;
                GetNextConnection(i, yBus, connToSlack);
            }
        }
    }
}
 
void ElectricCalculation::DistributeReactivePower(std::vector<ReactiveMachine>& machines, double qTotal)
{
    int freeMachines = machines.size();
    double qRemaining = qTotal;
 
    bool changed = true;
 
    while (changed && freeMachines > 0) {
        changed = false;
 
        double qShare = qRemaining / freeMachines;
 
        for (auto& m : machines) {
            if (m.fixed)
                continue;
 
            if (m.hasMax && qShare > m.qMax) {
                m.q = m.qMax;
                m.fixed = true;
                qRemaining -= m.qMax;
                freeMachines--;
                changed = true;
            }
            else if (m.hasMin && qShare < m.qMin) {
                m.q = m.qMin;
                m.fixed = true;
                qRemaining -= m.qMin;
                freeMachines--;
                changed = true;
            }
            else {
                m.q = qShare;
            }
        }
    }
 
    if (freeMachines > 0) {
        double qShare = qRemaining / freeMachines;
        for (auto& m : machines) {
            if (!m.fixed)
                m.q = qShare;
        }
    }
}
 
void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<double> > voltage,
    std::vector<std::complex<double> > power,
    std::vector<BusType> busType,
    std::vector<ReactiveLimits> reactiveLimit,
    double systemPowerBase)
{
    double zeroLimit = 1e-4;
    for (unsigned int i = 0; i < reactiveLimit.size(); ++i) {
        if (reactiveLimit[i].maxLimit > -zeroLimit && reactiveLimit[i].maxLimit < zeroLimit)
            reactiveLimit[i].maxLimit = zeroLimit;
        if (reactiveLimit[i].minLimit > -zeroLimit && reactiveLimit[i].minLimit < zeroLimit)
            reactiveLimit[i].minLimit = zeroLimit;
    }
    for (unsigned int i = 0; i < power.size(); ++i) {
        if (std::real(power[i]) > -zeroLimit && std::real(power[i]) < zeroLimit)
            power[i] = std::complex<double>(0.0, std::imag(power[i]));
        if (std::imag(power[i]) > -zeroLimit && std::imag(power[i]) < zeroLimit)
            power[i] = std::complex<double>(std::real(power[i]), 0.0);
    }
    // Buses
    for (int i = 0; i < (int)m_busList.size(); i++) {
        Bus* bus = m_busList[i];
        BusElectricalData data = bus->GetElectricalData();
        if (data.isConnected) {
            data.voltage = voltage[data.number];
            data.power = power[data.number];
            data.busType = busType[data.number];
        }
        else {
            data.voltage = std::complex<double>(0.0, 0.0);
            data.power = std::complex<double>(0.0, 0.0);
        }
        bus->SetElectricalData(data);
    }
 
    // Power line
    for (int i = 0; i < (int)m_lineList.size(); i++) {
        Line* line = m_lineList[i];
        if (line->IsOnline()) {
            auto dataBus1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData();
            auto dataBus2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData();
 
            if (dataBus1.isConnected && dataBus2.isConnected) {
                int n1 = dataBus1.number;
                int n2 = dataBus2.number;
 
                LineElectricalData data = line->GetElectricalData();
                LineElectricalData dataPU = line->GetPUElectricalData(systemPowerBase);
                std::complex<double> v1 = voltage[n1];
                std::complex<double> v2 = voltage[n2];
 
                data.current[0] = (v1 - v2) / std::complex<double>(dataPU.resistance, dataPU.indReactance) +
                    v1 * std::complex<double>(0.0, dataPU.capSusceptance / 2.0);
                data.current[1] = (v2 - v1) / std::complex<double>(dataPU.resistance, dataPU.indReactance) +
                    v2 * std::complex<double>(0.0, dataPU.capSusceptance / 2.0);
 
                data.powerFlow[0] = v1 * std::conj(data.current[0]);
                data.powerFlow[1] = v2 * std::conj(data.current[1]);
 
                if (data.powerFlow[0].real() > data.powerFlow[1].real())
                    line->SetPowerFlowDirection(PowerFlowDirection::PF_BUS1_TO_BUS2);
                else
                    line->SetPowerFlowDirection(PowerFlowDirection::PF_BUS2_TO_BUS1);
 
                line->SetElectricalData(data);
            }
        }
    }
 
    // Transformer
    for (int i = 0; i < (int)m_transformerList.size(); i++) {
        Transformer* transformer = m_transformerList[i];
        if (transformer->IsOnline()) {
            auto dataBus1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData();
            auto dataBus2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData();
 
            if (dataBus1.isConnected && dataBus2.isConnected) {
                TransformerElectricalData data = transformer->GetElectricalData();
                TransformerElectricalData dataPU = transformer->GetPUElectricalData(systemPowerBase);
 
                int n1 = dataBus1.number;
                int n2 = dataBus2.number;
 
                std::complex<double> v1 = voltage[n1];  // Primary voltage
                std::complex<double> v2 = voltage[n2];  // Secondary voltage
 
                // Transformer admitance
                std::complex<double> y = 1.0 / std::complex<double>(dataPU.resistance, dataPU.indReactance);
 
                if (dataPU.turnsRatio == 1.0 && dataPU.phaseShift == 0.0) {
                    data.current[0] = (v1 - v2) * y;
                    data.current[1] = (v2 - v1) * y;
                }
                else {
                    double radPS = wxDegToRad(dataPU.phaseShift);
                    std::complex<double> a =
                        std::complex<double>(dataPU.turnsRatio * std::cos(radPS), -dataPU.turnsRatio * std::sin(radPS));
 
                    data.current[0] = v1 * (y / std::pow(std::abs(a), 2)) - v2 * (y / std::conj(a));
                    data.current[1] = -v1 * (y / a) + v2 * y;
                }
 
                data.powerFlow[0] = v1 * std::conj(data.current[0]);
                data.powerFlow[1] = v2 * std::conj(data.current[1]);
 
                if (data.powerFlow[0].real() > data.powerFlow[1].real())
                    transformer->SetPowerFlowDirection(PowerFlowDirection::PF_BUS1_TO_BUS2);
                else
                    transformer->SetPowerFlowDirection(PowerFlowDirection::PF_BUS2_TO_BUS1);
 
                transformer->SetElectricaData(data);
            }
        }
    }
 
    // Ind motor
    for (unsigned int i = 0; i < m_indMotorList.size(); i++) {
        IndMotor* motor = m_indMotorList[i];
        auto data = motor->GetElectricalData();
        if (motor->IsOnline() && data.calcQInPowerFlow) {
            double reactivePower = data.qValue * systemPowerBase;
 
            switch (data.reactivePowerUnit) {
            case ElectricalUnit::UNIT_PU: {
                reactivePower /= systemPowerBase;
            } break;
            case ElectricalUnit::UNIT_kvar: {
                reactivePower /= 1e3;
            } break;
            case ElectricalUnit::UNIT_Mvar: {
                reactivePower /= 1e6;
            } break;
            default:
                break;
            }
 
            data.reactivePower = reactivePower;
 
            motor->SetElectricalData(data);
        }
    }
 
    //EMT Elements
    for (EMTElement* emtElement : m_emtElementList) {
 
        if (!emtElement->IsOnline()) {
            auto data = emtElement->GetEMTElementData();
            data.power = std::complex<double>(0.0, 0.0);
            data.currHarmonics[1] = std::complex<double>(0.0, 0.0);
            emtElement->SetEMTElementData(data);
        }
        else {
            emtElement->UpdateData();
            auto data = emtElement->GetEMTElementData();
            double baseCurrent = systemPowerBase / (sqrt(3.0) * data.baseVoltage);
            std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
            std::complex<double> power = data.puVoltage * std::conj(current);
            data.power = power;
            emtElement->SetEMTElementData(data);
        }
    }
 
    // Synchronous machines
    for (auto* bus : m_busList) {
        BusElectricalData data = bus->GetElectricalData();
        int i = data.number;
 
        if (data.isConnected) {
            // Get the synchronous machines connected and calculate the load power on the bus.
            std::vector<SyncGenerator*> syncGeneratorsOnBus;
            std::vector<SyncMotor*> syncMotorsOnBus;
            std::complex<double> loadPower(0.0, 0.0);
 
            for (auto itsg = m_syncGeneratorList.begin(); itsg != m_syncGeneratorList.end(); itsg++) {
                SyncGenerator* syncGenerator = *itsg;
                if (bus == syncGenerator->GetParentList()[0] && syncGenerator->IsOnline())
                    syncGeneratorsOnBus.push_back(syncGenerator);
            }
            for (auto itsm = m_syncMotorList.begin(); itsm != m_syncMotorList.end(); itsm++) {
                SyncMotor* syncMotor = *itsm;
                if (bus == syncMotor->GetParentList()[0] && syncMotor->IsOnline()) {
                    syncMotorsOnBus.push_back(syncMotor);
                    SyncMotorElectricalData childData = syncMotor->GetPUElectricalData(systemPowerBase);
                    loadPower += std::complex<double>(childData.activePower, 0.0);
                }
            }
            for (auto itlo = m_loadList.begin(); itlo != m_loadList.end(); itlo++) {
                Load* load = *itlo;
                if (bus == load->GetParentList()[0] && load->IsOnline()) {
                    LoadElectricalData childData = load->GetPUElectricalData(systemPowerBase);
                    if (childData.loadType == CONST_POWER)
                        loadPower += std::complex<double>(childData.activePower, childData.reactivePower);
 
                    if (childData.activePower >= 0.0)
                        load->SetPowerFlowDirection(PowerFlowDirection::PF_TO_ELEMENT);
                    else
                        load->SetPowerFlowDirection(PowerFlowDirection::PF_TO_BUS);
                }
            }
            for (auto itim = m_indMotorList.begin(); itim != m_indMotorList.end(); itim++) {
                IndMotor* indMotor = *itim;
                if (bus == indMotor->GetParentList()[0] && indMotor->IsOnline()) {
                    IndMotorElectricalData childData = indMotor->GetPUElectricalData(systemPowerBase);
                    loadPower += std::complex<double>(childData.activePower, childData.reactivePower);
 
                    if (childData.activePower >= 0.0)
                        indMotor->SetPowerFlowDirection(PowerFlowDirection::PF_TO_ELEMENT);
                    else
                        indMotor->SetPowerFlowDirection(PowerFlowDirection::PF_TO_BUS);
                }
            }
 
            // Setup the reactive machines vector for the distribution of reactive power
            std::vector<ReactiveMachine> machines;
            for (auto* gen : syncGeneratorsOnBus) {
 
                auto data = gen->GetPUElectricalData(systemPowerBase);
 
                ReactiveMachine m;
                m.qMax = data.maxReactive;
                m.qMin = data.minReactive;
 
                m.hasMax = data.haveMaxReactive;
                m.hasMin = data.haveMinReactive;
 
                m.machine = gen;
                m.isGenerator = true;
 
                machines.push_back(m);
            }
            for (auto* mot : syncMotorsOnBus) {
 
                auto data = mot->GetPUElectricalData(systemPowerBase);
 
                ReactiveMachine m;
                m.qMax = data.maxReactive;
                m.qMin = data.minReactive;
 
                m.hasMax = data.haveMaxReactive;
                m.hasMin = data.haveMinReactive;
 
                m.machine = mot;
                m.isGenerator = false;
 
                machines.push_back(m);
            }
 
            // Distribute the reactive power among the synchronous machines and set their reactive power
            double qTotal = power[i].imag() + loadPower.imag();
            DistributeReactivePower(machines, qTotal);
 
            // Set the reactive power on the machines
            for (auto& m : machines) {
 
                if (m.isGenerator) {
 
                    SyncGenerator* gen = static_cast<SyncGenerator*>(m.machine);
                    auto data = gen->GetElectricalData();
 
                    double reactivePower = m.q * systemPowerBase;
 
                    switch (data.reactivePowerUnit) {
                    case ElectricalUnit::UNIT_PU:
                        reactivePower /= systemPowerBase;
                        break;
                    case ElectricalUnit::UNIT_kvar:
                        reactivePower /= 1e3;
                        break;
                    case ElectricalUnit::UNIT_Mvar:
                        reactivePower /= 1e6;
                        break;
                    default:
                        break;
                    }
 
                    data.reactivePower = reactivePower;
                    gen->SetElectricalData(data);
                }
                else {
 
                    SyncMotor* mot = static_cast<SyncMotor*>(m.machine);
                    auto data = mot->GetElectricalData();
 
                    double reactivePower = m.q * systemPowerBase;
 
                    switch (data.reactivePowerUnit) {
                    case ElectricalUnit::UNIT_PU:
                        reactivePower /= systemPowerBase;
                        break;
                    case ElectricalUnit::UNIT_kvar:
                        reactivePower /= 1e3;
                        break;
                    case ElectricalUnit::UNIT_Mvar:
                        reactivePower /= 1e6;
                        break;
                    default:
                        break;
                    }
 
                    data.reactivePower = reactivePower;
                    mot->SetElectricalData(data);
                }
            }
        }
    }
}
 
bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix,
    std::vector<std::vector<std::complex<double> > >& inverse)
{
    int order = static_cast<int>(matrix.size());
 
    inverse.clear();
    // Fill the inverse matrix with identity.
    for (int i = 0; i < order; ++i) {
        std::vector<std::complex<double> > line;
        for (int j = 0; j < order; ++j) {
            line.push_back(i == j ? std::complex<double>(1.0, 0.0) : std::complex<double>(0.0, 0.0));
        }
        inverse.push_back(line);
    }
 
    // Check if a main diagonal value of the matrix is zero, if one is zero, try a linear combination to remove it.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i == j && matrix[i][j] == std::complex<double>(0.0, 0.0)) {
                int row = 0;
                while (row < order) {
                    if (matrix[row][j] != std::complex<double>(0.0, 0.0)) {
                        for (int k = 0; k < order; ++k) {
                            matrix[i][k] += matrix[row][k];
                            inverse[i][k] += inverse[row][k];
                        }
                        break;
                    }
                    row++;
                }
                // If all line values are zero, the matrix is singular and the solution is impossible.
                if (row == order) return false;
            }
        }
    }
 
    // Linear combinations are made in both matrices, the goal is the input matrix become the identity. The final result
    // have two matrices: the identity and the inverse of the input.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i != j) {
                if (matrix[i][i] == std::complex<double>(0.0, 0.0)) return false;
 
                std::complex<double> factor = matrix[j][i] / matrix[i][i];
                for (int k = 0; k < order; ++k) {
                    matrix[j][k] -= factor * matrix[i][k];
                    inverse[j][k] -= factor * inverse[i][k];
                }
            }
        }
    }
    // Main diagonal calculation.
    for (int i = 0; i < order; ++i) {
        for (int j = 0; j < order; ++j) {
            if (i == j) {
                if (matrix[i][j] == std::complex<double>(0.0, 0.0)) return false;
 
                std::complex<double> factor = (matrix[i][j] - std::complex<double>(1.0, 0.0)) / matrix[i][j];
                for (int k = 0; k < order; ++k) {
                    matrix[j][k] -= factor * matrix[i][k];
                    inverse[j][k] -= factor * inverse[i][k];
                }
            }
        }
    }
 
    return true;
}
 
void ElectricCalculation::ABCtoDQ0(std::complex<double> complexValue, double angle, double& dValue, double& qValue)
{
    dValue = -std::real(complexValue) * std::sin(angle) + std::imag(complexValue) * std::cos(angle);
    qValue = std::real(complexValue) * std::cos(angle) + std::imag(complexValue) * std::sin(angle);
}
 
void ElectricCalculation::DQ0toABC(double dValue, double qValue, double angle, std::complex<double>& complexValue)
{
    double real = qValue * std::cos(angle) - dValue * std::sin(angle);
    double imag = qValue * std::sin(angle) + dValue * std::cos(angle);
    complexValue = std::complex<double>(real, imag);
}
 
std::vector<std::complex<double> > ElectricCalculation::GaussianElimination(
    std::vector<std::vector<std::complex<double> > > matrix,
    std::vector<std::complex<double> > array)
{
    //[Ref] https://en.wikipedia.org/wiki/Gaussian_elimination
 
    std::vector<std::complex<double> > solution;
 
    std::vector<std::vector<std::complex<double> > > triangMatrix;
    triangMatrix.resize(matrix.size());
    for (unsigned int i = 0; i < matrix.size(); i++) { triangMatrix[i].resize(matrix.size()); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) { solution.push_back(array[i]); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) {
        for (unsigned int j = 0; j < matrix.size(); j++) { triangMatrix[i][j] = matrix[i][j]; }
    }
 
    for (unsigned int k = 0; k < matrix.size(); k++) {
        unsigned int k1 = k + 1;
        for (unsigned int i = k; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != std::complex<double>(0.0, 0.0)) {
                for (unsigned int j = k1; j < matrix.size(); j++) {
                    triangMatrix[i][j] = triangMatrix[i][j] / triangMatrix[i][k];
                }
                solution[i] = solution[i] / triangMatrix[i][k];
            }
        }
        for (unsigned int i = k1; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != std::complex<double>(0.0, 0.0)) {
                for (unsigned int j = k1; j < matrix.size(); j++) { triangMatrix[i][j] -= triangMatrix[k][j]; }
                solution[i] -= solution[k];
            }
        }
    }
    for (int i = static_cast<int>(matrix.size()) - 2; i >= 0; i--) {
        for (int j = static_cast<int>(matrix.size()) - 1; j >= i + 1; j--) {
            solution[i] -= triangMatrix[i][j] * solution[j];
        }
    }
 
    return solution;
}
 
std::vector<double> ElectricCalculation::GaussianElimination(std::vector<std::vector<double> > matrix,
    std::vector<double> array)
{
    //[Ref] https://en.wikipedia.org/wiki/Gaussian_elimination
 
    std::vector<double> solution;
 
    std::vector<std::vector<double> > triangMatrix;
    triangMatrix.resize(matrix.size());
    for (unsigned int i = 0; i < matrix.size(); i++) { triangMatrix[i].resize(matrix.size()); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) { solution.push_back(array[i]); }
 
    for (unsigned int i = 0; i < matrix.size(); i++) {
        for (unsigned int j = 0; j < matrix.size(); j++) { triangMatrix[i][j] = matrix[i][j]; }
    }
 
    for (unsigned int k = 0; k < matrix.size(); k++) {
        unsigned int k1 = k + 1;
        for (unsigned int i = k; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != 0.0) {
                for (unsigned int j = k1; j < matrix.size(); j++) {
                    triangMatrix[i][j] = triangMatrix[i][j] / triangMatrix[i][k];
                }
                solution[i] = solution[i] / triangMatrix[i][k];
            }
        }
        for (unsigned int i = k1; i < matrix.size(); i++) {
            if (triangMatrix[i][k] != 0.0) {
                for (unsigned int j = k1; j < matrix.size(); j++) { triangMatrix[i][j] -= triangMatrix[k][j]; }
                solution[i] -= solution[k];
            }
        }
    }
    for (int i = static_cast<int>(matrix.size()) - 2; i >= 0; i--) {
        for (int j = static_cast<int>(matrix.size()) - 1; j >= i + 1; j--) {
            solution[i] -= triangMatrix[i][j] * solution[j];
        }
    }
 
    return solution;
}
 
Machines::SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* generator)
{
    auto data = generator->GetElectricalData();
    if (data.transTd0 != 0.0) {
        if (data.transTq0 != 0.0) {
            if (data.subTd0 != 0.0 || data.subTq0 != 0.0) { return Machines::SM_MODEL_5; }
            return Machines::SM_MODEL_3;
        }
        else {
            if (data.subTd0 != 0.0 || data.subTq0 != 0.0) { return Machines::SM_MODEL_4; }
            return Machines::SM_MODEL_2;
        }
    }
 
    return Machines::SM_MODEL_1;
}
 
std::vector<std::complex<double> > ElectricCalculation::ComplexMatrixTimesVector(
    std::vector<std::vector<std::complex<double> > > matrix,
    std::vector<std::complex<double> > vector)
{
    std::vector<std::complex<double> > solution;
    for (unsigned int i = 0; i < matrix.size(); i++) {
        solution.push_back(std::complex<double>(0.0, 0.0));
 
        for (unsigned int j = 0; j < matrix.size(); j++) { solution[i] += matrix[i][j] * vector[j]; }
    }
 
    return solution;
}
 
void ElectricCalculation::GetLUDecomposition(std::vector<std::vector<std::complex<double> > > matrix,
    std::vector<std::vector<std::complex<double> > >& matrixL,
    std::vector<std::vector<std::complex<double> > >& matrixU)
{
    // Doolittle method
    // [Ref] http://www3.nd.edu/~zxu2/acms40390F11/Alg-LU-Crout.pdf
    // [Ref] http://www.engr.colostate.edu/~thompson/hPage/CourseMat/Tutorials/CompMethods/doolittle.pdf
 
    int size = static_cast<int>(matrix.size());  // Decomposed matrix size.
 
    // Set upper and lower matrices sizes.
    matrixL.resize(size);
    matrixU.resize(size);
    for (int i = 0; i < size; i++) {
        matrixL[i].resize(size);
        matrixU[i].resize(size);
    }
 
    // First row of upper matrix and first column of lower matrix.
    for (int i = 0; i < size; i++) {
        matrixU[0][i] = matrix[0][i];
        matrixL[i][0] = matrix[i][0] / matrixU[0][0];
    }
 
    // Lower matrix main diagonal.
    for (int i = 1; i < size; i++) { matrixL[i][i] = std::complex<double>(1.0, 0.0); }
 
    for (int i = 1; i < size - 1; i++) {
        // Upper matrix main diagonal.
        matrixU[i][i] = matrix[i][i];
        for (int k = 0; k < i; k++) { matrixU[i][i] -= matrixL[i][k] * matrixU[k][i]; }
 
        // Others elements of upper matrix
        for (int j = i + 1; j < size; j++) {
            matrixU[i][j] = matrix[i][j];
            for (int k = 0; k < i; k++) { matrixU[i][j] -= matrixL[i][k] * matrixU[k][j]; }
        }
 
        // Lower matrix elements
        for (int j = i + 1; j < size; j++) {
            matrixL[j][i] = matrix[j][i];
            for (int k = 0; k < i; k++) { matrixL[j][i] -= matrixL[j][k] * matrixU[k][i]; }
            matrixL[j][i] = matrixL[j][i] / matrixU[i][i];
        }
    }
 
    // Last element of upper matrix.
    matrixU[size - 1][size - 1] = matrix[size - 1][size - 1];
    for (int k = 0; k < size - 1; k++) { matrixU[size - 1][size - 1] -= matrixL[size - 1][k] * matrixU[k][size - 1]; }
}
 
std::vector<std::complex<double> > ElectricCalculation::LUEvaluate(std::vector<std::vector<std::complex<double> > > u,
    std::vector<std::vector<std::complex<double> > > l,
    std::vector<std::complex<double> > b)
{
    int size = static_cast<int>(b.size());
    std::vector<std::complex<double> > x;
    std::vector<std::complex<double> > y;
    x.resize(size);
    y.resize(size);
 
    // Forward
    for (int i = 0; i < size; i++) {
        y[i] = b[i];
        for (int j = 0; j < i; j++) { y[i] -= l[i][j] * y[j]; }
        y[i] /= l[i][i];
    }
    // Backward
    for (int i = size - 1; i >= 0; i--) {
        x[i] = y[i];
        for (int j = i + 1; j < size; j++) { x[i] -= u[i][j] * x[j]; }
        x[i] /= u[i][i];
    }
    return x;
}
 
bool ElectricCalculation::GetParentBus(Element* childElement, Bus*& parentBus)
{
    auto parentList = childElement->GetParentList();
    if (parentList.size() < 1) return false;
    parentBus = static_cast<Bus*>(parentList[0]);
    if (parentBus == nullptr) return false;
    if (parentBus->GetElectricalData().number < 0) return false;
    return true;
}
 
bool ElectricCalculation::GetParentBus(Element* childElement, Bus*& parentBus1, Bus*& parentBus2)
{
    auto parentList = childElement->GetParentList();
    if (parentList.size() < 2) return false;
    parentBus1 = static_cast<Bus*>(parentList[0]);
    parentBus2 = static_cast<Bus*>(parentList[1]);
    if (parentBus1 == nullptr || parentBus2 == nullptr) return false;
    if (parentBus1->GetElectricalData().number < 0 || parentBus2->GetElectricalData().number < 0) return false;
    return true;
}
 
bool ElectricCalculation::CalculateEMTElementsAdmittance(const double& basePower, wxString& errorMsg)
{
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        emtElement->UpdateData();
        if (!emtElement->CalculateCurrent(errorMsg)) return false;
        auto data = emtElement->GetEMTElementData();
        double baseCurrent = basePower / (sqrt(3.0) * data.baseVoltage);
        std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
        std::complex<double> y0 = current / data.puVoltage;
        data.y0 = y0;
        data.power = std::complex<double>(0.0, 0.0);
        emtElement->SetEMTElementData(data);
    }
    return true;
}
 
bool ElectricCalculation::CalculateEMTElementsPower(const double& basePower, wxString& errorMsg, bool updateCurrent)
{
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        emtElement->UpdateData();
        if (updateCurrent) {
            if (!emtElement->CalculateCurrent(errorMsg)) return false;
        }
        auto data = emtElement->GetEMTElementData();
        double baseCurrent = basePower / (sqrt(3.0) * data.baseVoltage);
        //std::complex<double> current = data.currHarmonics[1];
        std::complex<double> current = data.currHarmonics[1] / baseCurrent; // Fundamental frequency current
        
        std::complex<double> i0 = data.y0 * data.puVoltage; // Current already injected via admittance matrix
        //std::complex<double> voltage = data.baseVoltage * data.puVoltage;
        //std::complex<double> power = (sqrt(3.0) * voltage * std::conj(current)) / basePower;
        std::complex<double> power = data.puVoltage * std::conj(current - i0);
        data.powerDiff = power - data.power;
        data.power = power;
        emtElement->SetEMTElementData(data);
 
        //wxString currentStr = data.name + ":\n";
        //for (auto const& [harm, current] : data.currHarmonics)
        //{
        //  currentStr += wxString::Format("I(%dh) = %f p.u.\n", harm, abs(current) / baseCurrent);
        //}
        //currentStr += wxString::Format("S = %e + j%e p.u.\n", data.power.real(), data.power.imag());
        //currentStr += wxString::Format("SDiff = %e p.u.\n", abs(data.powerDiff));
        //wxMessageBox(currentStr);
    }
    return true;
}
 
double ElectricCalculation::CalculateEMTPowerError(const std::vector<std::complex<double>>& voltage, std::vector<std::complex<double>>& power, const double& basePower, wxString& errorMsg)
{
    // Update buses voltages
    for (Bus* bus : m_busList) {
        auto data = bus->GetElectricalData();
        if (data.isConnected) {
            data.voltage = voltage[data.number];
            bus->SetElectricalData(data);
        }
    }
 
    CalculateEMTElementsPower(basePower, errorMsg);
    // Update power array and calculate error
    double error = 0.0;
    for (EMTElement* emtElement : m_emtElementList) {
        if (!emtElement->IsOnline()) continue;
 
        if (!emtElement->GetParentList().empty()) {
            if (emtElement->GetParentList()[0] != nullptr) {
                int numBus = static_cast<Bus*>(emtElement->GetParentList()[0])->GetElectricalData().number;
                auto data = emtElement->GetEMTElementData();
                power[numBus] += data.powerDiff;
                if (error < std::abs(data.powerDiff))
                    error = std::abs(data.powerDiff);
            }
        }
    }
 
    return error;
}