Removing sync generator is now working, bus plot implemented

6 files changed, 157 insertions, 80 deletions

diff --git a/Project/Bus.cpp b/Project/Bus.cpp
index 69bca6d..3ccb55d 100644
--- a/Project/Bus.cpp
+++ b/Project/Bus.cpp
@@ -3,12 +3,8 @@
 #include "DegreesAndRadians.h"
 #endif
 
-Bus::Bus()
-    : PowerElement()
-{
-}
-Bus::Bus(wxPoint2DDouble position)
-    : PowerElement()
+Bus::Bus() : PowerElement() {}
+Bus::Bus(wxPoint2DDouble position) : PowerElement()
 {
     m_width = 100.0;
     m_height = 5.0;
@@ -41,14 +37,14 @@ void Bus::Draw(wxPoint2DDouble translation, double scale) const
 
         glColor4d(0.0, 0.5, 1.0, 0.5);
 
-        wxPoint2DDouble pts[4] = { WorldToScreen(translation, scale, -(m_width / 2.0), -(m_height / 2.0)) -
-                wxPoint2DDouble(m_borderSize, m_borderSize),
-            WorldToScreen(translation, scale, -(m_width / 2.0), (m_height / 2.0)) -
-                wxPoint2DDouble(m_borderSize, -m_borderSize),
-            WorldToScreen(translation, scale, (m_width / 2.0), (m_height / 2.0)) -
-                wxPoint2DDouble(-m_borderSize, -m_borderSize),
-            WorldToScreen(translation, scale, (m_width / 2.0), -(m_height / 2.0)) -
-                wxPoint2DDouble(-m_borderSize, m_borderSize) };
+        wxPoint2DDouble pts[4] = {WorldToScreen(translation, scale, -(m_width / 2.0), -(m_height / 2.0)) -
+                                      wxPoint2DDouble(m_borderSize, m_borderSize),
+                                  WorldToScreen(translation, scale, -(m_width / 2.0), (m_height / 2.0)) -
+                                      wxPoint2DDouble(m_borderSize, -m_borderSize),
+                                  WorldToScreen(translation, scale, (m_width / 2.0), (m_height / 2.0)) -
+                                      wxPoint2DDouble(-m_borderSize, -m_borderSize),
+                                  WorldToScreen(translation, scale, (m_width / 2.0), -(m_height / 2.0)) -
+                                      wxPoint2DDouble(-m_borderSize, m_borderSize)};
         DrawRectangle(pts);
         glPopMatrix();
     }
@@ -75,8 +71,8 @@ void Bus::Draw(wxPoint2DDouble translation, double scale) const
         glRotated(m_angle, 0.0, 0.0, 1.0);
         glTranslated(-screenPt.m_x, -screenPt.m_y, 0.0);
 
-        wxPoint2DDouble pbPosition[2] = { WorldToScreen(translation, scale, m_width / 2.0),
-            WorldToScreen(translation, scale, -m_width / 2.0) };
+        wxPoint2DDouble pbPosition[2] = {WorldToScreen(translation, scale, m_width / 2.0),
+                                         WorldToScreen(translation, scale, -m_width / 2.0)};
         DrawPickbox(pbPosition[0]);
         DrawPickbox(pbPosition[1]);
 
@@ -204,12 +200,12 @@ wxString Bus::GetTipText() const
     tipText += wxString::Format(" (%d)", m_electricalData.number + 1);
     tipText += "\n";
     tipText += StringFromDouble(m_electricalData.nominalVoltage, 1) +
-        (m_electricalData.nominalVoltageUnit == UNIT_V ? _(" V") : _(" kV"));
+               (m_electricalData.nominalVoltageUnit == UNIT_V ? _(" V") : _(" kV"));
     tipText += "\n";
     tipText += _("\nV = ") + wxString::FromDouble(std::abs(m_electricalData.voltage), 5) + _(" p.u.");
     tipText += "\n";
     tipText += wxString(L'\u03B8') + " = " + wxString::FromDouble(wxRadToDeg(std::arg(m_electricalData.voltage)), 5) +
-        " " + wxString(L'\u00B0');
+               " " + wxString(L'\u00B0');
 
     tipText += _("\n\nFault info:");
     tipText += _("\nVa = ") + wxString::FromDouble(std::abs(m_electricalData.faultVoltage[0]), 5) + _(" p.u.");
@@ -220,8 +216,24 @@ wxString Bus::GetTipText() const
         tipText += _("\nIb = ") + wxString::FromDouble(std::abs(m_electricalData.faultCurrent[1]), 5) + _(" p.u.");
         tipText += _("\nIc = ") + wxString::FromDouble(std::abs(m_electricalData.faultCurrent[2]), 5) + _(" p.u.");
     }
-    
+
     tipText += _("\n\nSsc = ") + wxString::FromDouble(std::abs(m_electricalData.scPower), 5) + _(" p.u.");
 
     return tipText;
 }
+
+bool Bus::GetPlotData(ElementPlotData& plotData)
+{
+    if(!m_electricalData.plotBus) return false;
+    plotData.SetName(m_electricalData.name);
+    plotData.SetCurveType(ElementPlotData::CT_BUS);
+
+    std::vector<double> absVoltage, argVoltage;
+    for(unsigned int i = 0; i < m_electricalData.stabVoltageVector.size(); ++i) {
+        absVoltage.push_back(std::abs(m_electricalData.stabVoltageVector[i]));
+        argVoltage.push_back(wxRadToDeg(std::arg(m_electricalData.stabVoltageVector[i])));
+    }
+    plotData.AddData(absVoltage, _("Voltage"));
+    plotData.AddData(argVoltage, _("Angle"));
+    return true;
+}
diff --git a/Project/Bus.h b/Project/Bus.h
index 3f802df..cbc8287 100644
--- a/Project/Bus.h
+++ b/Project/Bus.h
@@ -11,7 +11,7 @@ struct BusElectricalData {
     ElectricalUnit nominalVoltageUnit = UNIT_kV;
     bool isVoltageControlled = false;
     double controlledVoltage = 1.0;
-    int controlledVoltageUnitChoice = 0; // 0 = p.u., 1 = same as nominalVoltageUnit (UNIT_V or UNIT_kV).
+    int controlledVoltageUnitChoice = 0;  // 0 = p.u., 1 = same as nominalVoltageUnit (UNIT_V or UNIT_kV).
     bool slackBus = false;
 
     // Power flow (p.u.)
@@ -24,10 +24,10 @@ struct BusElectricalData {
     // p.u. fault data
     double faultResistance = 0.0;
     double faultReactance = 0.0;
-    std::complex<double> faultCurrent[3] = { std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0),
-        std::complex<double>(0.0, 0.0) };
-    std::complex<double> faultVoltage[3] = { std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0),
-        std::complex<double>(0.0, 0.0) };
+    std::complex<double> faultCurrent[3] = {std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0),
+                                            std::complex<double>(0.0, 0.0)};
+    std::complex<double> faultVoltage[3] = {std::complex<double>(0.0, 0.0), std::complex<double>(0.0, 0.0),
+                                            std::complex<double>(0.0, 0.0)};
     double scPower = 0.0;
 
     // Stability
@@ -37,11 +37,12 @@ struct BusElectricalData {
     double stabFaultLength = 0.0;
     double stabFaultResistance = 0.0;
     double stabFaultReactance = 0.0;
+    std::vector<std::complex<double> > stabVoltageVector;
 };
 
 class Bus : public PowerElement
 {
-public:
+   public:
     Bus();
     Bus(wxPoint2DDouble position);
     Bus(wxPoint2DDouble position, wxString name);
@@ -60,9 +61,10 @@ public:
     virtual BusElectricalData GetElectricalData() const { return m_electricalData; }
     virtual void SetElectricalData(BusElectricalData electricalData) { m_electricalData = electricalData; }
     virtual bool ShowForm(wxWindow* parent, Element* element);
+    virtual bool GetPlotData(ElementPlotData& plotData);
 
-protected:
+   protected:
     BusElectricalData m_electricalData;
 };
 
-#endif // BUS_H
+#endif  // BUS_H
diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp
index 6cf5d9e..a81912f 100644
--- a/Project/Electromechanical.cpp
+++ b/Project/Electromechanical.cpp
@@ -60,10 +60,8 @@ bool Electromechanical::RunStabilityCalculation()
 
         if(currentPbdTime > pbdTime) {
             if(!pbd.Update((currentTime / simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) {
+                m_errorMsg = wxString::Format(_("Simulation cancelled at %.2fs."), currentTime);
                 pbd.Update(100);
-                currentTime = simTime;
-                
-                m_errorMsg = _("Simulation cancelled.");
                 return false;
             }
             currentPbdTime = 0.0;
@@ -417,6 +415,13 @@ std::complex<double> Electromechanical::GetSyncMachineAdmittance(SyncGenerator*
 
 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;
@@ -670,8 +675,16 @@ bool Electromechanical::SolveSynchronousMachines()
             // Calculate integration constants.
             CalculateIntegrationConstants(syncGenerator, id, iq);
         }
+        else {
+            CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f);
+        }
     }
 
+    m_wError = 0;
+    m_deltaError = 0;
+    m_transEdError = 0;
+    m_transEqError = 0;
+
     double error = 1.0;
     int iterations = 0;
     while(error > m_tolerance) {
@@ -688,64 +701,79 @@ bool Electromechanical::SolveSynchronousMachines()
             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;
-                }
-                int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
+            double k = 1.0;  // Power base change factor.
+            if(data.useMachineBase) {
+                double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit);
+                k = m_powerSystemBase / oldBase;
+            }
+            int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number;
 
+            if(syncGenerator->IsOnline()) {
                 data.terminalVoltage = m_vBus[n];
+            }
 
-                // Mechanical differential equations.
-                double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe);
-                error = std::max(error, std::abs(data.speed - w) / w0);
+            // Mechanical differential equations.
+            double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe);
+            error = std::max(error, std::abs(data.speed - w) / w0);
 
-                double delta = data.icDelta.c + data.icDelta.m * w;
-                error = std::max(error, std::abs(data.delta - delta));
+            m_wError += std::abs(data.speed - w) / w0;
 
-                data.speed = w;
-                data.delta = delta;
+            double delta = data.icDelta.c + data.icDelta.m * w;
+            error = std::max(error, std::abs(data.delta - delta));
 
-                // Electric power.
-                double id, iq, vd, vq;
-                std::complex<double> iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]);
+            data.speed = w;
+            data.delta = delta;
+
+            m_deltaError += std::abs(data.delta - delta);
 
+            // Electric power.
+            double id, iq, vd, vq, pe;
+            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);
-                ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq);
-
-                double pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k;
-                // data.pe = (2 * pe - data.pe);  // Extrapolating Pe.
-                data.pe = pe;  // Don't extrapolating Pe
-
-                // Electrical differential equations
-                switch(GetMachineModel(syncGenerator)) {
-                    case SM_MODEL_1: {
-                        // There is no differential equations.
-                    } break;
-                    case SM_MODEL_2: {
-                    } break;
-                    case SM_MODEL_3: {
-                        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));
-
-                        double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq;
-                        error = std::max(error, std::abs(data.tranEd - tranEd));
-
-                        data.tranEq = tranEq;
-                        data.tranEd = tranEd;
-
-                    } break;
-                    case SM_MODEL_4: {
-                    } break;
-                    case SM_MODEL_5: {
-                    } break;
-                }
+
+                pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k;
+                // pe = (2 * pe - data.pe);  // Extrapolating Pe.
             } else {
-                // Set values to open circuit machine.
+                pe = id = iq = 0.0f;
+            }
+
+            data.pe = pe;
+
+            // Electrical differential equations
+            switch(GetMachineModel(syncGenerator)) {
+                case SM_MODEL_1: {
+                    // There is no differential equations.
+                } break;
+                case SM_MODEL_2: {
+                } break;
+                case SM_MODEL_3: {
+                    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));
+
+                    m_transEqError += std::abs(data.tranEq - tranEq);
+
+                    double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq;
+                    error = std::max(error, std::abs(data.tranEd - tranEd));
+
+                    m_transEdError += 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 SM_MODEL_4: {
+                } break;
+                case SM_MODEL_5: {
+                } break;
             }
 
             syncGenerator->SetElectricalData(data);
@@ -758,6 +786,7 @@ bool Electromechanical::SolveSynchronousMachines()
             return false;
         }
     }
+    m_numIt = iterations;
 
     // Solve controllers.
     for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
@@ -833,4 +862,18 @@ void Electromechanical::SaveData()
             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_deltaErrorVector.push_back(m_deltaError);
+    m_transEdErrorVector.push_back(m_transEdError);
+    m_transEqErrorVector.push_back(m_transEqError);
+    m_numItVector.push_back(m_numIt);
 }
diff --git a/Project/Electromechanical.h b/Project/Electromechanical.h
index 808b936..8dabe05 100644
--- a/Project/Electromechanical.h
+++ b/Project/Electromechanical.h
@@ -18,6 +18,12 @@ class Electromechanical : public ElectricCalculation
     
     std::vector<double> GetTimeVector() const { return m_timeVector; }
     
+    std::vector<double> m_wErrorVector;
+    std::vector<double> m_deltaErrorVector;
+    std::vector<double> m_transEdErrorVector;
+    std::vector<double> m_transEqErrorVector;
+    std::vector<double> m_numItVector;
+    
    protected:
     void SetEventTimeList();
     bool HasEvent(double currentTime);
@@ -57,6 +63,13 @@ class Electromechanical : public ElectricCalculation
     std::vector<bool> m_eventOccurrenceList;
     
     std::vector<double> m_timeVector;
+    
+    //tests
+    double m_wError = 0.0;
+    double m_deltaError = 0.0;
+    double m_transEdError = 0.0;
+    double m_transEqError = 0.0;
+    double m_numIt = 0;
 };
 
 #endif  // ELECTROMECHANICAL_H
diff --git a/Project/Project.mk b/Project/Project.mk
index de32785..a5c96c3 100644
--- a/Project/Project.mk
+++ b/Project/Project.mk
@@ -13,7 +13,7 @@ CurrentFileName        :=
 CurrentFilePath        :=
 CurrentFileFullPath    :=
 User                   :=Thales
-Date                   :=28/05/2017
+Date                   :=29/05/2017
 CodeLitePath           :="C:/Program Files (x86)/CodeLite"
 LinkerName             :=C:/TDM-GCC-64/bin/g++.exe
 SharedObjectLinkerName :=C:/TDM-GCC-64/bin/g++.exe -shared -fPIC
diff --git a/Project/Workspace.cpp b/Project/Workspace.cpp
index 31d6849..f2caee8 100644
--- a/Project/Workspace.cpp
+++ b/Project/Workspace.cpp
@@ -1416,6 +1416,13 @@ bool Workspace::RunStability()
         ElementPlotData plotData;
         if(element->GetPlotData(plotData)) plotDataList.push_back(plotData);
     }
+    ElementPlotData tests(_("Error"), ElementPlotData::CT_TEST);
+    tests.AddData(stability.m_wErrorVector, _("Speed error"));
+    tests.AddData(stability.m_deltaErrorVector, _("Delta error"));
+    tests.AddData(stability.m_transEqErrorVector, _("Eq error"));
+    tests.AddData(stability.m_transEdErrorVector, _("Ed error"));
+    tests.AddData(stability.m_numItVector, _("Number iterations"));
+    plotDataList.push_back(tests);
 
     ChartView* cView = new ChartView(this, plotDataList, stability.GetTimeVector());
     cView->Show();