Merge pull request #33 from Thales1330/wip/electromechanical-calc

Wip electromechanical calc

1 files changed, 1094 insertions, 0 deletions

diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp
new file mode 100644
index 0000000..57a275a
--- /dev/null
+++ b/Project/Electromechanical.cpp
@@ -0,0 +1,1094 @@
+#include "Electromechanical.h"
+#include "ControlElementSolver.h"
+
+Electromechanical::Electromechanical(wxWindow* parent, std::vector<Element*> elementList, SimulationData data)
+{
+    m_parent = parent;
+    GetElementsFromList(elementList);
+    SetEventTimeList();
+
+    m_powerSystemBase = GetPowerValue(data.basePower, data.basePowerUnit);
+    m_systemFreq = data.stabilityFrequency;
+    m_simTime = data.stabilitySimulationTime;
+    m_timeStep = data.timeStep;
+    m_tolerance = data.stabilityTolerance;
+    m_maxIterations = data.stabilityMaxIterations;
+
+    m_ctrlTimeStepMultiplier = 1.0 / static_cast<double>(data.controlTimeStepRatio);
+
+    m_plotTime = data.plotTime;
+    m_useCOI = data.useCOI;
+}
+
+Electromechanical::~Electromechanical() {}
+bool Electromechanical::RunStabilityCalculation()
+{
+    wxProgressDialog pbd(_("Running simulation"), _("Initializing..."), 100, m_parent,
+                         wxPD_APP_MODAL | wxPD_AUTO_HIDE | wxPD_CAN_ABORT | wxPD_SMOOTH);
+
+    SetSyncMachinesModel();
+
+    // Calculate the admittance matrix with the synchronous machines.
+    if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) {
+        m_errorMsg = _("It was not possible to build the admittance matrix.");
+        return false;
+    }
+    InsertSyncMachinesOnYBus();
+    GetLUDecomposition(m_yBus, m_yBusL, m_yBusU);
+
+    // Get buses voltages.
+    m_vBus.clear();
+    m_vBus.resize(m_busList.size());
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetElectricalData();
+        m_vBus[data.number] = data.voltage;
+    }
+
+    // Calculate injected currents
+    m_iBus = ComplexMatrixTimesVector(m_yBus, m_vBus);
+    for(unsigned int i = 0; i < m_iBus.size(); ++i) {
+        if(std::abs(m_iBus[i]) < 1e-5) m_iBus[i] = std::complex<double>(0.0, 0.0);
+    }
+
+    if(!InitializeDynamicElements()) return false;
+
+    double pbdTime = m_plotTime;
+    double currentTime = 0.0;
+    double currentPlotTime = 0.0;
+    double currentPbdTime = 0.0;
+    while(currentTime < m_simTime) {
+        if(HasEvent(currentTime)) {
+            SetEvent(currentTime);
+            GetLUDecomposition(m_yBus, m_yBusL, m_yBusU);
+        }
+
+        if(currentPlotTime >= m_plotTime || currentTime == 0.0) {
+            m_timeVector.push_back(currentTime);
+            SaveData();
+            currentPlotTime = 0.0;
+        }
+
+        if(currentPbdTime > pbdTime) {
+            if(!pbd.Update((currentTime / m_simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) {
+                m_errorMsg = wxString::Format(_("Simulation cancelled at %.2fs."), currentTime);
+                pbd.Update(100);
+                return false;
+            }
+            currentPbdTime = 0.0;
+        }
+
+        if(!SolveSynchronousMachines()) return false;
+
+        currentTime += m_timeStep;
+        currentPlotTime += m_timeStep;
+        currentPbdTime += m_timeStep;
+    }
+    return true;
+}
+
+void Electromechanical::SetEventTimeList()
+{
+    // Fault
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetElectricalData();
+        if(data.stabHasFault) {
+            m_eventTimeList.push_back(data.stabFaultTime);
+            m_eventOccurrenceList.push_back(false);
+            m_eventTimeList.push_back(data.stabFaultTime + data.stabFaultLength);
+            m_eventOccurrenceList.push_back(false);
+        }
+    }
+    // Switching
+    for(auto it = m_powerElementList.begin(), itEnd = m_powerElementList.end(); it != itEnd; ++it) {
+        PowerElement* element = *it;
+        SwitchingData swData = element->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swTime.size(); ++i) {
+            m_eventTimeList.push_back(swData.swTime[i]);
+            m_eventOccurrenceList.push_back(false);
+        }
+    }
+}
+
+bool Electromechanical::HasEvent(double currentTime)
+{
+    for(unsigned int i = 0; i < m_eventTimeList.size(); ++i) {
+        if(!m_eventOccurrenceList[i]) {
+            if(EventTrigger(m_eventTimeList[i], currentTime)) {
+                m_eventOccurrenceList[i] = true;
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+void Electromechanical::SetEvent(double currentTime)
+{
+    // Fault
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetElectricalData();
+        if(data.stabHasFault) {
+            int n = data.number;
+
+            // Insert fault
+            if(EventTrigger(data.stabFaultTime, currentTime)) {
+                double r, x;
+                r = data.stabFaultResistance;
+                x = data.stabFaultReactance;
+                if(x < 1e-5) x = 1e-5;
+                m_yBus[n][n] += std::complex<double>(1.0, 0.0) / std::complex<double>(r, x);
+            }
+
+            // Remove fault
+            else if(EventTrigger(data.stabFaultTime + data.stabFaultLength, currentTime)) {
+                double r, x;
+                r = data.stabFaultResistance;
+                x = data.stabFaultReactance;
+                if(x < 1e-5) x = 1e-5;
+                m_yBus[n][n] -= std::complex<double>(1.0, 0.0) / std::complex<double>(r, x);
+            }
+        }
+    }
+
+    // SyncGenerator switching
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* generator = *it;
+        auto swData = generator->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove machine (only connected machines)
+                if(swData.swType[i] == SW_REMOVE && generator->IsOnline()) {
+                    generator->SetOnline(false);
+                    int n = static_cast<Bus*>(generator->GetParentList()[0])->GetElectricalData().number;
+                    m_yBus[n][n] -= GetSyncMachineAdmittance(generator);
+                }
+
+                // Insert machine (only disconnected machines)
+                if(swData.swType[i] == SW_INSERT && !generator->IsOnline() && generator->GetParentList().size() == 1) {
+                    if(generator->SetOnline(true)) {
+                        int n = static_cast<Bus*>(generator->GetParentList()[0])->GetElectricalData().number;
+                        m_yBus[n][n] += GetSyncMachineAdmittance(generator);
+                    }
+                }
+            }
+        }
+    }
+
+    // Load switching
+    for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) {
+        Load* load = *it;
+        auto swData = load->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove load (only connected loads)
+                if(swData.swType[i] == SW_REMOVE && load->IsOnline()) {
+                    load->SetOnline(false);
+                    auto data = load->GetPUElectricalData(m_powerSystemBase);
+                    Bus* parentBus = static_cast<Bus*>(load->GetParentList()[0]);
+                    int n = parentBus->GetElectricalData().number;
+                    std::complex<double> v = parentBus->GetElectricalData().voltage;
+                    m_yBus[n][n] -= std::complex<double>(data.activePower, -data.reactivePower) / (v * v);
+                }
+
+                // Insert load (only disconnected load)
+                if(swData.swType[i] == SW_INSERT && !load->IsOnline() && load->GetParentList().size() == 1) {
+                    if(load->SetOnline(true)) {
+                        auto data = load->GetPUElectricalData(m_powerSystemBase);
+                        Bus* parentBus = static_cast<Bus*>(load->GetParentList()[0]);
+                        int n = parentBus->GetElectricalData().number;
+                        std::complex<double> v = parentBus->GetElectricalData().voltage;
+                        m_yBus[n][n] += std::complex<double>(data.activePower, -data.reactivePower) / (v * v);
+                    }
+                }
+            }
+        }
+    }
+
+    // Line switching
+    for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) {
+        Line* line = *it;
+        auto swData = line->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove line (only connected lines)
+                if(swData.swType[i] == SW_REMOVE && line->IsOnline()) {
+                    line->SetOnline(false);
+                    auto data = line->GetElectricalData();
+
+                    int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number;
+                    int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number;
+
+                    m_yBus[n1][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    m_yBus[n2][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                    m_yBus[n1][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    m_yBus[n2][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                    m_yBus[n1][n1] -= std::complex<double>(0.0, data.capSusceptance / 2.0);
+                    m_yBus[n2][n2] -= std::complex<double>(0.0, data.capSusceptance / 2.0);
+                }
+
+                // Insert line (only disconnected lines)
+                if(swData.swType[i] == SW_INSERT && !line->IsOnline() && line->GetParentList().size() == 2) {
+                    if(line->SetOnline(true)) {
+                        auto data = line->GetElectricalData();
+
+                        int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number;
+                        int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number;
+
+                        m_yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                        m_yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                        m_yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                        m_yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                        m_yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+                        m_yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+                    }
+                }
+            }
+        }
+    }
+
+    // Transformer switching
+    for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) {
+        Transformer* transformer = *it;
+        auto swData = transformer->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove transformer (only connected transformers)
+                if(swData.swType[i] == SW_REMOVE && transformer->IsOnline()) {
+                    transformer->SetOnline(false);
+                    auto data = transformer->GetElectricalData();
+
+                    int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number;
+                    int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().number;
+
+                    if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) {
+                        m_yBus[n1][n2] -= -1.0 / std::complex<double>(data.resistance, data.indReactance);
+                        m_yBus[n2][n1] -= -1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                        m_yBus[n1][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                        m_yBus[n2][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    } else {
+                        // 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);
+                        m_yBus[n1][n1] -= y / std::pow(std::abs(a), 2.0);
+                        m_yBus[n1][n2] -= -(y / std::conj(a));
+                        m_yBus[n2][n1] -= -(y / a);
+                        m_yBus[n2][n2] -= y;
+                    }
+                }
+
+                // Insert transformer (only disconnected transformers)
+                if(swData.swType[i] == SW_INSERT && !transformer->IsOnline() &&
+                   transformer->GetParentList().size() == 2) {
+                    if(transformer->SetOnline(true)) {
+                        auto data = transformer->GetElectricalData();
+
+                        int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number;
+                        int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().number;
+
+                        if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) {
+                            m_yBus[n1][n2] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
+                            m_yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
+
+                            m_yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                            m_yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                        } else {
+                            // 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);
+                            m_yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
+                            m_yBus[n1][n2] += -(y / std::conj(a));
+                            m_yBus[n2][n1] += -(y / a);
+                            m_yBus[n2][n2] += y;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // Capacitor switching
+    for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) {
+        Capacitor* capacitor = *it;
+        auto swData = capacitor->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove capacitor (only connected capacitors)
+                if(swData.swType[i] == SW_REMOVE && capacitor->IsOnline()) {
+                    capacitor->SetOnline(false);
+                    auto data = capacitor->GetPUElectricalData(m_powerSystemBase);
+                    int n = static_cast<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number;
+                    m_yBus[n][n] -= std::complex<double>(0.0, data.reactivePower);
+                }
+
+                // Insert capacitor (only disconnected capacitors)
+                if(swData.swType[i] == SW_INSERT && !capacitor->IsOnline() && capacitor->GetParentList().size() == 1) {
+                    if(capacitor->SetOnline(true)) {
+                        auto data = capacitor->GetPUElectricalData(m_powerSystemBase);
+                        int n = static_cast<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number;
+                        m_yBus[n][n] += std::complex<double>(0.0, data.reactivePower);
+                    }
+                }
+            }
+        }
+    }
+
+    // Inductor switching
+    for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) {
+        Inductor* inductor = *it;
+        auto swData = inductor->GetSwitchingData();
+        for(unsigned int i = 0; i < swData.swType.size(); ++i) {
+            if(EventTrigger(swData.swTime[i], currentTime)) {
+                // Remove inductor (only connected inductors)
+                if(swData.swType[i] == SW_REMOVE && inductor->IsOnline()) {
+                    inductor->SetOnline(false);
+                    auto data = inductor->GetPUElectricalData(m_powerSystemBase);
+                    int n = static_cast<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number;
+                    m_yBus[n][n] -= std::complex<double>(0.0, -data.reactivePower);
+                }
+
+                // Insert inductor (only disconnected inductors)
+                if(swData.swType[i] == SW_INSERT && !inductor->IsOnline() && inductor->GetParentList().size() == 1) {
+                    if(inductor->SetOnline(true)) {
+                        auto data = inductor->GetPUElectricalData(m_powerSystemBase);
+                        int n = static_cast<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number;
+                        m_yBus[n][n] += std::complex<double>(0.0, -data.reactivePower);
+                    }
+                }
+            }
+        }
+    }
+}
+
+void Electromechanical::InsertSyncMachinesOnYBus()
+{
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* generator = *it;
+        if(generator->IsOnline()) {
+            auto data = generator->GetElectricalData();
+            int n = static_cast<Bus*>(generator->GetParentList()[0])->GetElectricalData().number;
+            m_yBus[n][n] += GetSyncMachineAdmittance(generator);
+        }
+    }
+}
+
+bool Electromechanical::EventTrigger(double eventTime, double currentTime)
+{
+    return (((eventTime - m_timeStep) < currentTime) && (eventTime >= currentTime));
+}
+
+std::complex<double> Electromechanical::GetSyncMachineAdmittance(SyncGenerator* generator)
+{
+    auto data = generator->GetElectricalData();
+    double k = 1.0;  // Power base change factor.
+    if(data.useMachineBase) {
+        double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+        k = m_powerSystemBase / oldBase;
+    }
+
+    double xd = 0.0;
+    double xq = 0.0;
+    double ra = data.armResistance * k;
+
+    switch(data.model) {
+        case Machines::SM_MODEL_1: {
+            xq = data.transXd * k;
+            xd = xq;
+        } break;
+        case Machines::SM_MODEL_2: {
+            xd = data.transXd * k;
+            xq = data.transXq * k;
+            if(xq == 0.0) {
+                xq = data.syncXq * k;
+                if(xq == 0.0) {
+                    xq = data.syncXd * k;
+                }
+            }
+        } break;
+        case Machines::SM_MODEL_3: {
+            xd = data.transXd * k;
+            xq = data.transXq * k;
+            if(xq == 0.0) xq = xd;
+        } break;
+        case Machines::SM_MODEL_4:
+        case Machines::SM_MODEL_5: {
+            xd = data.subXd * k;
+            xq = data.subXq * k;
+            if(xd == 0.0) xd = xq;
+            if(xq == 0.0) xq = xd;
+        } break;
+    }
+    double xdq = 0.5 * (xd + xq);
+    return (std::complex<double>(ra, -xdq) / std::complex<double>(ra * ra + xd * xq, 0.0));
+}
+
+bool Electromechanical::InitializeDynamicElements()
+{
+    // Buses
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetElectricalData();
+        data.stabVoltageVector.clear();
+        bus->SetElectricalData(data);
+    }
+    // Synchronous generators
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto dataPU = syncGenerator->GetPUElectricalData(m_powerSystemBase);
+        auto data = syncGenerator->GetElectricalData();
+        if(syncGenerator->IsOnline()) {
+            double k = 1.0;  // Power base change factor.
+            if(data.useMachineBase) {
+                double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                k = m_powerSystemBase / oldBase;
+            }
+            data.terminalVoltage = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().voltage;
+
+            std::complex<double> conjS(dataPU.activePower, -dataPU.reactivePower);
+            std::complex<double> conjV = std::conj(data.terminalVoltage);
+            std::complex<double> ia = conjS / conjV;
+
+            double xd = data.syncXd * k;
+            double xq = data.syncXq * k;
+            double ra = data.armResistance * k;
+
+            if(data.model == Machines::SM_MODEL_1) {
+                xq = data.transXd * k;
+                xd = xq;
+            } else if(data.syncXq == 0.0)
+                xq = data.syncXd * k;
+
+            // Initialize state variables
+            std::complex<double> eq0 = data.terminalVoltage + std::complex<double>(ra, xq) * ia;
+            data.delta = std::arg(eq0);
+
+            double fi0 = std::arg(ia);
+            double id0, iq0;
+            // ABCtoDQ0(ia, data.delta - fi0, id0, iq0);
+            iq0 = std::abs(ia) * std::cos(data.delta - fi0);
+            id0 = -std::abs(ia) * std::sin(data.delta - fi0);
+
+            data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0;
+            data.fieldVoltage = data.initialFieldVoltage;
+            data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra));
+            data.speed = 2.0 * M_PI * m_systemFreq;
+
+            data.pe = data.pm;
+            data.electricalPower = std::complex<double>(dataPU.activePower, dataPU.reactivePower);
+
+            // Variables to extrapolate.
+            data.oldIq = iq0;
+            data.oldId = id0;
+            data.oldPe = data.pe;
+
+            switch(data.model) {
+                case Machines::SM_MODEL_1: {
+                    // double tranXd = data.transXd * k;
+
+                    // data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0;
+                    data.tranEq = std::abs(eq0);
+
+                    data.tranEd = 0.0;
+                    data.subEq = 0.0;
+                    data.subEd = 0.0;
+                } break;
+                case Machines::SM_MODEL_2: {
+                    double tranXd = data.transXd * k;
+
+                    data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0;
+                    data.tranEd = 0.0;
+                    data.subEd = 0.0;
+                    data.subEq = 0.0;
+                } break;
+                case Machines::SM_MODEL_3: {
+                    double tranXd = data.transXd * k;
+                    double tranXq = data.transXq * k;
+                    if(tranXq == 0.0) tranXq = tranXd;
+
+                    data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0;
+                    data.tranEd = -(xq - tranXq) * iq0;
+
+                    data.subEd = 0.0;
+                    data.subEq = 0.0;
+                } break;
+                case Machines::SM_MODEL_4: {
+                    double tranXd = data.transXd * k;
+                    double subXd = data.subXd * k;
+                    double subXq = data.subXq * k;
+                    if(subXd == 0.0) subXd = subXq;
+                    if(subXq == 0.0) subXq = subXd;
+
+                    data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0;
+                    data.tranEd = 0.0;
+                    data.subEq = data.tranEq + (tranXd - subXd) * id0;
+                    data.subEd = -(xq - subXq) * iq0;
+                } break;
+                case Machines::SM_MODEL_5: {
+                    double tranXd = data.transXd * k;
+                    double tranXq = data.transXq * k;
+                    double subXd = data.subXd * k;
+                    double subXq = data.subXq * k;
+                    if(subXd == 0.0) subXd = subXq;
+                    if(subXq == 0.0) subXq = subXd;
+
+                    data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0;
+                    data.tranEd = -(xq - tranXq) * iq0;
+                    data.subEq = data.tranEq + (tranXd - subXd) * id0;
+                    data.subEd = data.tranEd - (tranXq - subXq) * iq0;
+                } break;
+                default:
+                    break;
+            }
+
+            // Initialize controllers
+            if(data.useAVR) {
+                if(data.avrSolver) delete data.avrSolver;
+                data.avrSolver = new ControlElementSolver(data.avr, m_timeStep * m_ctrlTimeStepMultiplier, m_tolerance,
+                                                          false, std::abs(data.terminalVoltage), m_parent);
+                if(!data.avrSolver->IsOK()) {
+                    m_errorMsg = _("Error on initializate the AVR of \"") + data.name + _("\".");
+                    syncGenerator->SetElectricalData(data);
+                    return false;
+                }
+            }
+            if(data.useSpeedGovernor) {
+                if(data.speedGovSolver) delete data.speedGovSolver;
+                data.speedGovSolver = new ControlElementSolver(data.speedGov, m_timeStep * m_ctrlTimeStepMultiplier,
+                                                               m_tolerance, false, data.speed, m_parent);
+                if(!data.speedGovSolver->IsOK()) {
+                    m_errorMsg = _("Error on initializate the speed governor of \"") + data.name + _("\".");
+                    syncGenerator->SetElectricalData(data);
+                    return false;
+                }
+            }
+        } else {
+            // Initialize open circuit machine.
+        }
+        // Reset plot data
+        data.terminalVoltageVector.clear();
+        data.electricalPowerVector.clear();
+        data.mechanicalPowerVector.clear();
+        data.freqVector.clear();
+        data.fieldVoltageVector.clear();
+        data.deltaVector.clear();
+
+        syncGenerator->SetElectricalData(data);
+    }
+    CalculateReferenceSpeed();
+    return true;
+}
+
+void Electromechanical::CalculateMachinesCurrents()
+{
+    // Reset injected currents vector
+    for(unsigned int i = 0; i < m_iBus.size(); ++i) m_iBus[i] = std::complex<double>(0.0, 0.0);
+
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto data = syncGenerator->GetElectricalData();
+        if(syncGenerator->IsOnline()) {
+            double k = 1.0;  // Power base change factor.
+            if(data.useMachineBase) {
+                double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                k = m_powerSystemBase / oldBase;
+            }
+
+            double xd = 0.0;
+            double xq = 0.0;
+            double ra = data.armResistance * k;
+            int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
+            std::complex<double> e = std::complex<double>(0.0, 0.0);
+            std::complex<double> v = m_vBus[n];
+            std::complex<double> iInj = m_iBus[n];
+            double xdq = 0.0;
+
+            switch(data.model) {
+                case Machines::SM_MODEL_1: {
+                    DQ0toABC(data.tranEd, data.tranEq, data.delta, e);
+                    xq = data.transXd * k;
+                    xd = xq;
+                } break;
+                case Machines::SM_MODEL_2: {
+                    DQ0toABC(data.tranEd, data.tranEq, data.delta, e);
+                    xd = data.transXd * k;
+                    xq = data.transXq * k;
+                    if(xq == 0.0) {
+                        xq = data.syncXq * k;
+                        if(xq == 0.0) {
+                            xq = data.syncXd * k;
+                        }
+                    }
+                } break;
+                case Machines::SM_MODEL_3: {
+                    DQ0toABC(data.tranEd, data.tranEq, data.delta, e);
+                    xd = data.transXd * k;
+                    xq = data.transXq * k;
+                    if(xq == 0.0) xq = xd;
+                } break;
+                case Machines::SM_MODEL_4:
+                case Machines::SM_MODEL_5: {
+                    DQ0toABC(data.subEd, data.subEq, data.delta, e);
+                    xd = data.subXd * k;
+                    xq = data.subXq * k;
+                    if(xd == 0.0) xd = xq;
+                    if(xq == 0.0) xq = xd;
+                } break;
+            }
+            xdq = 0.5 * (xd + xq);
+
+            std::complex<double> y0 = std::complex<double>(ra, -xdq) / std::complex<double>(ra * ra + xd * xq, 0.0);
+            std::complex<double> iUnadj = y0 * e;
+
+            std::complex<double> iAdj =
+                std::complex<double>(0.0, -((0.5 * (xq - xd)) / (ra * ra + xd * xq))) * (std::conj(e) - std::conj(v));
+            iAdj = iAdj * std::cos(2.0 * data.delta) + iAdj * std::complex<double>(0.0, std::sin(2.0 * data.delta));
+
+            iInj = iUnadj + iAdj;
+
+            m_iBus[n] += iInj;
+
+            std::complex<double> iMachine = iInj - y0 * v;
+
+            data.electricalPower = v * std::conj(iMachine);
+        } else {
+            data.electricalPower = std::complex<double>(0.0, 0.0);
+        }
+
+        syncGenerator->SetElectricalData(data);
+    }
+}
+
+void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq, double k)
+{
+    CalculateReferenceSpeed();
+    auto data = syncGenerator->GetElectricalData();
+
+    double syncXd, syncXq, transXd, transXq, subXd, subXq;
+    syncXd = data.syncXd * k;
+    syncXq = data.syncXq * k;
+    transXd = data.transXd * k;
+    transXq = data.transXq * k;
+    subXd = data.subXd * k;
+    subXq = data.subXq * k;
+
+    if(syncXq == 0.0) syncXq = syncXd;
+    if(transXq == 0.0) transXq = transXd;
+    if(subXd == 0.0) subXd = subXq;
+    if(subXq == 0.0) subXq = subXd;
+
+    double transTd0, transTq0, subTd0, subTq0;
+    transTd0 = data.transTd0;
+    transTq0 = data.transTq0;
+    subTd0 = data.subTd0;
+    subTq0 = data.subTq0;
+
+    if(subTd0 == 0.0) subTd0 = subTq0;
+    if(subTq0 == 0.0) subTq0 = subTd0;
+
+    // Speed
+    data.icSpeed.m = m_timeStep / ((4.0f * data.inertia / m_refSpeed) / k + m_timeStep * data.damping * k);
+    data.icSpeed.c = (1.0f - 2.0f * data.icSpeed.m * data.damping * k) * data.speed +
+                     data.icSpeed.m * (data.pm - data.pe + 2.0f * m_refSpeed * data.damping * k);
+
+    // Delta
+    data.icDelta.m = 0.5f * m_timeStep;
+    data.icDelta.c = data.delta + data.icDelta.m * (data.speed - 2.0f * m_refSpeed);
+
+    // Eq'
+    if(data.model == Machines::SM_MODEL_2 || data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 ||
+       data.model == Machines::SM_MODEL_5) {
+        data.icTranEq.m = m_timeStep / (2.0f * transTd0 + m_timeStep);
+        data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq +
+                          data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id);
+    }
+
+    // Ed'
+    if(data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) {
+        data.icTranEd.m = m_timeStep / (2.0f * transTq0 + m_timeStep);
+        data.icTranEd.c = (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (syncXq - transXq) * iq;
+    }
+
+    // Eq''
+    if(data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) {
+        data.icSubEq.m = m_timeStep / (2.0f * subTd0 + m_timeStep);
+        data.icSubEq.c =
+            (1.0f - 2.0f * data.icSubEq.m) * data.subEq + data.icSubEq.m * (data.tranEq + (transXd - subXd) * id);
+    }
+    // Ed''
+    if(data.model == Machines::SM_MODEL_4) {
+        data.icSubEd.m = m_timeStep / (2.0f * subTq0 + m_timeStep);
+        data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd - data.icSubEd.m * (syncXq - subXq) * iq;
+    }
+    if(data.model == Machines::SM_MODEL_5) {
+        data.icSubEd.m = m_timeStep / (2.0f * subTq0 + m_timeStep);
+        data.icSubEd.c =
+            (1.0f - 2.0f * data.icSubEd.m) * data.subEd + data.icSubEd.m * (data.tranEd - (transXq - subXq) * iq);
+    }
+
+    syncGenerator->SetElectricalData(data);
+}
+
+bool Electromechanical::SolveSynchronousMachines()
+{
+    // CalculateMachinesCurrents();
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto data = syncGenerator->GetElectricalData();
+
+        if(syncGenerator->IsOnline()) {
+            int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
+            double id, iq, pe;
+
+            pe = data.pe;
+
+            double k = 1.0;  // Power base change factor.
+            if(data.useMachineBase) {
+                double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                k = m_powerSystemBase / oldBase;
+            }
+
+            std::complex<double> iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]);
+
+            ABCtoDQ0(iMachine, data.delta, id, iq);
+
+            // Calculate integration constants.
+            CalculateIntegrationConstants(syncGenerator, id, iq, k);
+
+            CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k);
+            // Extrapolate nonintegrable variables.
+            id = 2.0 * id - data.oldId;
+            iq = 2.0 * iq - data.oldIq;
+            pe = 2.0 * pe - data.oldPe;
+
+            CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k);
+        } else {
+            CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f);
+        }
+    }
+
+    m_wError = 0;
+
+    double error = 1.0;
+    int iterations = 0;
+    while(error > m_tolerance) {
+        error = 0.0;
+
+        // Calculate the injected currents.
+        CalculateMachinesCurrents();
+
+        // Calculate the buses voltages.
+        m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus);
+
+        // Solve machine equations.
+        for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+            SyncGenerator* syncGenerator = *it;
+
+            auto data = syncGenerator->GetElectricalData();
+
+            double id, iq, pe;
+            double k = 1.0;  // Power base change factor.
+            if(data.useMachineBase) {
+                double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                k = m_powerSystemBase / oldBase;
+            }
+
+            CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k);
+
+            double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k);
+
+            if(genError > error) error = genError;
+        }
+
+        ++iterations;
+
+        if(iterations > m_maxIterations) {
+            m_errorMsg = _("Impossible to solve the synchronous generators.\nCheck the system parameters and/or "
+                           "decrease the time step.");
+            return false;
+        }
+    }
+    m_numIt = iterations;
+
+    // Solve controllers.
+    int ctrlRatio = static_cast<int>(1 / m_ctrlTimeStepMultiplier);
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto data = syncGenerator->GetElectricalData();
+        if(data.useAVR && data.avrSolver) {
+            for(int i = 0; i < ctrlRatio; ++i) data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage));
+            data.fieldVoltage = data.initialFieldVoltage + data.avrSolver->GetLastSolution();
+        }
+
+        if(data.useSpeedGovernor && data.speedGovSolver) {
+            for(int i = 0; i < ctrlRatio; ++i) data.speedGovSolver->SolveNextStep(data.speed);
+            data.pm = data.speedGovSolver->GetLastSolution();
+        }
+        syncGenerator->SetElectricalData(data);
+    }
+
+    return true;
+}
+
+double Electromechanical::GetPowerValue(double value, ElectricalUnit unit)
+{
+    switch(unit) {
+        case UNIT_PU: {
+            return value;
+        } break;
+        case UNIT_VA: {
+            return value;
+        } break;
+        case UNIT_kVA: {
+            return value * 1e3;
+        } break;
+        case UNIT_MVA: {
+            return value * 1e6;
+        } break;
+        case UNIT_W: {
+            return value;
+        } break;
+        case UNIT_kW: {
+            return value * 1e3;
+        } break;
+        case UNIT_MW: {
+            return value * 1e6;
+        } break;
+        case UNIT_VAr: {
+            return value;
+        } break;
+        case UNIT_kVAr: {
+            return value * 1e3;
+        } break;
+        case UNIT_MVAr: {
+            return value * 1e6;
+        } break;
+        default: {
+            return 0.0;
+        } break;
+    }
+    return 0.0;
+}
+
+void Electromechanical::SaveData()
+{
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto data = syncGenerator->GetElectricalData();
+        if(data.plotSyncMachine) {
+            data.terminalVoltageVector.push_back(data.terminalVoltage);
+            data.electricalPowerVector.push_back(data.electricalPower);
+            data.mechanicalPowerVector.push_back(data.pm);
+            data.freqVector.push_back(data.speed / (2.0f * M_PI));
+            data.fieldVoltageVector.push_back(data.fieldVoltage);
+            data.deltaVector.push_back(wxRadToDeg(data.delta));
+            syncGenerator->SetElectricalData(data);
+        }
+    }
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetElectricalData();
+        if(data.plotBus) {
+            data.stabVoltageVector.push_back(m_vBus[data.number]);
+            bus->SetElectricalData(data);
+        }
+    }
+
+    m_wErrorVector.push_back(m_wError);
+    m_numItVector.push_back(m_numIt);
+}
+
+void Electromechanical::SetSyncMachinesModel()
+{
+    for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+        SyncGenerator* syncGenerator = *it;
+        auto data = syncGenerator->GetElectricalData();
+        data.model = GetMachineModel(syncGenerator);
+        syncGenerator->SetElectricalData(data);
+    }
+}
+
+void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator,
+                                                            double& id,
+                                                            double& iq,
+                                                            double& pe,
+                                                            double k)
+{
+    auto data = syncGenerator->GetElectricalData();
+    int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
+
+    if(syncGenerator->IsOnline()) {
+        data.terminalVoltage = m_vBus[n];
+    }
+
+    double vd, vq;
+    ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq);
+
+    if(syncGenerator->IsOnline()) {
+        std::complex<double> iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]);
+        ABCtoDQ0(iMachine, data.delta, id, iq);
+
+        pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k;
+    } else {
+        pe = id = iq = 0.0f;
+    }
+    data.pe = pe;
+    data.oldId = id;
+    data.oldIq = iq;
+    syncGenerator->SetElectricalData(data);
+}
+
+double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGenerator,
+                                                           double id,
+                                                           double iq,
+                                                           double pe,
+                                                           double k)
+{
+    double error = 0.0;
+    auto data = syncGenerator->GetElectricalData();
+
+    // Mechanical differential equations.
+    double w = data.icSpeed.c + data.icSpeed.m * (data.pm - pe);
+    error = std::max(error, std::abs(data.speed - w) / m_refSpeed);
+
+    m_wError += std::abs(data.speed - w) / m_refSpeed;
+
+    double delta = data.icDelta.c + data.icDelta.m * w;
+    error = std::max(error, std::abs(data.delta - delta));
+
+    data.speed = w;
+    data.delta = delta;
+
+    // Electrical differential equations
+    switch(data.model) {
+        case Machines::SM_MODEL_1: {
+            // There is no differential equations.
+        } break;
+        case Machines::SM_MODEL_2: {
+            double tranEq =
+                data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id);
+            error = std::max(error, std::abs(data.tranEq - tranEq));
+
+            data.tranEq = tranEq;
+        } break;
+        case Machines::SM_MODEL_3: {
+            double syncXd, syncXq, transXd, transXq;
+            syncXd = data.syncXd * k;
+            syncXq = data.syncXq * k;
+            transXd = data.transXd * k;
+            transXq = data.transXq * k;
+            if(syncXq == 0.0) syncXq = syncXd;
+            if(transXq == 0.0) transXq = transXd;
+
+            double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id);
+            error = std::max(error, std::abs(data.tranEq - tranEq));
+
+            double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq;
+            error = std::max(error, std::abs(data.tranEd - tranEd));
+
+            data.tranEq = tranEq;
+            data.tranEd = tranEd;
+
+            if(!syncGenerator->IsOnline()) {
+                std::complex<double> e;
+                DQ0toABC(data.tranEd, data.tranEq, data.delta, e);
+                data.terminalVoltage = e;
+            }
+        } break;
+        case Machines::SM_MODEL_4: {
+            double syncXd, syncXq, transXd, subXd, subXq;
+            syncXd = data.syncXd * k;
+            syncXq = data.syncXq * k;
+            transXd = data.transXd * k;
+            subXd = data.subXd * k;
+            subXq = data.subXq * k;
+            if(syncXq == 0.0) syncXq = syncXd;
+            if(subXd == 0.0) subXd = subXq;
+            if(subXq == 0.0) subXq = subXd;
+
+            double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id);
+            error = std::max(error, std::abs(data.tranEq - tranEq));
+
+            double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id);
+            error = std::max(error, std::abs(data.subEq - subEq));
+
+            double subEd = data.icSubEd.c - data.icSubEd.m * (syncXq - subXq) * iq;
+            error = std::max(error, std::abs(data.subEd - subEd));
+
+            data.tranEq = tranEq;
+            data.subEq = subEq;
+            data.subEd = subEd;
+        } break;
+        case Machines::SM_MODEL_5: {
+            double syncXd, syncXq, transXd, transXq, subXd, subXq;
+            syncXd = data.syncXd * k;
+            syncXq = data.syncXq * k;
+            transXd = data.transXd * k;
+            transXq = data.transXq * k;
+            subXd = data.subXd * k;
+            subXq = data.subXq * k;
+            if(syncXq == 0.0) syncXq = syncXd;
+            if(transXq == 0.0) transXq = transXd;
+            if(subXd == 0.0) subXd = subXq;
+            if(subXq == 0.0) subXq = subXd;
+
+            double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id);
+            error = std::max(error, std::abs(data.tranEq - tranEq));
+
+            double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq;
+            error = std::max(error, std::abs(data.tranEd - tranEd));
+
+            double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id);
+            error = std::max(error, std::abs(data.subEq - subEq));
+
+            double subEd = data.icSubEd.c + data.icSubEd.m * (tranEd - (transXq - subXq) * iq);
+            error = std::max(error, std::abs(data.subEd - subEd));
+
+            data.tranEq = tranEq;
+            data.tranEd = tranEd;
+            data.subEq = subEq;
+            data.subEd = subEd;
+        } break;
+    }
+
+    syncGenerator->SetElectricalData(data);
+    return error;
+}
+
+void Electromechanical::CalculateReferenceSpeed()
+{
+    if(m_useCOI) {
+        double sumHW = 0.0;
+        double sumH = 0.0;
+        for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+            SyncGenerator* syncGenerator = *it;
+            if(syncGenerator->IsOnline()) {
+                auto data = syncGenerator->GetElectricalData();
+                double k = 1.0;  // Power base change factor.
+                if(data.useMachineBase) {
+                    double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                    k = m_powerSystemBase / oldBase;
+                }
+                sumH += data.inertia / k;
+                sumHW += data.inertia * data.speed / k;
+            }
+        }
+        m_refSpeed = sumHW / sumH;
+    } else {
+        m_refSpeed = 2.0 * M_PI * m_systemFreq;
+    }
+}