Electromechanical class and several methods implemented

11 files changed, 434 insertions, 106 deletions

diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp
index d024364..86653d5 100644
--- a/Project/ElectricCalculation.cpp
+++ b/Project/ElectricCalculation.cpp
@@ -7,6 +7,7 @@ ElectricCalculation::ElectricCalculation() {}
 ElectricCalculation::~ElectricCalculation() {}
 void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
 {
+    m_powerElementList.clear();
     m_busList.clear();
     m_capacitorList.clear();
     m_indMotorList.clear();
@@ -19,6 +20,8 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
     // 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))
@@ -41,9 +44,10 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
 }
 
 bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus,
-    double systemPowerBase,
-    YBusSequence sequence,
-    bool includeSyncMachines)
+                                  double systemPowerBase,
+                                  YBusSequence sequence,
+                                  bool includeSyncMachines,
+                                  bool allLoadsAsImpedances)
 {
     if(m_busList.size() == 0) return false;
 
@@ -73,7 +77,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
         if(load->IsOnline()) {
             int n = static_cast<Bus*>(load->GetParentList()[0])->GetEletricalData().number;
             LoadElectricalData data = load->GetPUElectricalData(systemPowerBase);
-            if(data.loadType == CONST_IMPEDANCE)
+            if(data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances)
                 yBus[n][n] += std::complex<double>(data.activePower, -data.reactivePower);
         }
     }
@@ -160,8 +164,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
             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));
+                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);
@@ -181,10 +185,11 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
                 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));
+                            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);
@@ -195,14 +200,14 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
                     case DELTA_GWYE: {
                         std::complex<double> y =
                             1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.secondaryGrndResistance),
-                                      data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance));
+                                                       data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance));
                         yBus[n2][n2] += y;
                         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));
+                                                       data.zeroIndReactance + 3.0 * (data.primaryGrndReactance));
                         yBus[n1][n1] += y;
                         break;
                     }
@@ -258,10 +263,10 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
 }
 
 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)
+                                                  std::vector<std::complex<double> > power,
+                                                  std::vector<BusType> busType,
+                                                  std::vector<ReactiveLimits> reactiveLimit,
+                                                  double systemPowerBase)
 {
     for(int i = 0; i < (int)reactiveLimit.size(); ++i) {
         if(reactiveLimit[i].maxLimit > -1e-5 && reactiveLimit[i].maxLimit < 1e-5) reactiveLimit[i].maxLimit = 1e-5;
@@ -287,9 +292,9 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl
             std::complex<double> v2 = voltage[n2];
 
             data.current[0] = (v1 - v2) / std::complex<double>(data.resistance, data.indReactance) +
-                v1 * std::complex<double>(0.0, data.capSusceptance / 2.0);
+                              v1 * std::complex<double>(0.0, data.capSusceptance / 2.0);
             data.current[1] = (v2 - v1) / std::complex<double>(data.resistance, data.indReactance) +
-                v2 * std::complex<double>(0.0, data.capSusceptance / 2.0);
+                              v2 * std::complex<double>(0.0, data.capSusceptance / 2.0);
 
             data.powerFlow[0] = v1 * std::conj(data.current[0]);
             data.powerFlow[1] = v2 * std::conj(data.current[1]);
@@ -310,8 +315,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl
             TransformerElectricalData data = transformer->GetElectricalData();
             int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetEletricalData().number;
             int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetEletricalData().number;
-            std::complex<double> v1 = voltage[n1]; // Primary voltage
-            std::complex<double> v2 = voltage[n2]; // Secondary voltage
+            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>(data.resistance, data.indReactance);
@@ -491,8 +496,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl
                         reactivePower = childData_PU.minReactive;
                         reachedMachineLimit = true;
                     } else if((!childData_PU.haveMaxReactive && reactiveLimit[i].limitReached == RL_MAX_REACHED) ||
-                        (!childData_PU.haveMinReactive && reactiveLimit[i].limitReached == RL_MIN_REACHED) ||
-                        (!childData_PU.haveMaxReactive && !childData_PU.haveMaxReactive)) {
+                              (!childData_PU.haveMinReactive && reactiveLimit[i].limitReached == RL_MIN_REACHED) ||
+                              (!childData_PU.haveMaxReactive && !childData_PU.haveMaxReactive)) {
                         reactivePower += exceededReactive;
                         exceededReactive = 0.0;
                     }
@@ -532,7 +537,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl
 }
 
 bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix,
-    std::vector<std::vector<std::complex<double> > >& inverse)
+                                       std::vector<std::vector<std::complex<double> > >& inverse)
 {
     int order = static_cast<int>(matrix.size());
 
@@ -599,3 +604,193 @@ bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<doub
 
     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] http://pt.wikipedia.org/wiki/Elimina%C3%A7%C3%A3o_de_Gauss
+
+    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(unsigned int i = matrix.size() - 2; i >= 0; i--) {
+        for(unsigned int j = matrix.size() - 1; j >= i + 1; j--) {
+            solution[i] -= triangMatrix[i][j] * solution[j];
+        }
+    }
+
+    return solution;
+}
+
+SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* generator)
+{
+    auto data = generator->GetElectricalData();
+    if(data.transTd0 != 0.0) {
+        if(data.transTq0 != 0.0) {
+            if(data.subTd0 != 0.0) {
+                if(data.subTq0 != 0.0)
+                    return SM_MODEL_5;
+                else
+                    return SM_MODEL_4;
+            } else
+                return SM_MODEL_3;
+        } else
+            return SM_MODEL_2;
+    }
+
+    return 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];
+    }
+
+    // Upper 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;
+}
diff --git a/Project/ElectricCalculation.h b/Project/ElectricCalculation.h
index 3eea5dd..be9aa89 100644
--- a/Project/ElectricCalculation.h
+++ b/Project/ElectricCalculation.h
@@ -5,6 +5,7 @@
 #include <complex>
 
 #include "Element.h"
+#include "PowerElement.h"
 #include "Bus.h"
 #include "Capacitor.h"
 #include "IndMotor.h"
@@ -18,16 +19,18 @@
 enum BusType { BUS_SLACK = 0, BUS_PV, BUS_PQ };
 
 enum ReactiveLimitsType {
-    RL_UNLIMITED = 0,    // The bus can generate any ammount of reactive power.
-    RL_LIMITED,          // The bus reactive power generation is limited.
-    RL_UNLIMITED_SOURCE, // The bus have at least one source of infinite reative power.
-    RL_MAX_REACHED,      // Max limit reached
-    RL_MIN_REACHED,      // Min limit reached
-    RL_NONE_REACHED      // No limits reached
+    RL_UNLIMITED = 0,     // The bus can generate any ammount of reactive power.
+    RL_LIMITED,           // The bus reactive power generation is limited.
+    RL_UNLIMITED_SOURCE,  // The bus have at least one source of infinite reative power.
+    RL_MAX_REACHED,       // Max limit reached
+    RL_MIN_REACHED,       // Min limit reached
+    RL_NONE_REACHED       // No limits reached
 };
 
 enum YBusSequence { POSITIVE_SEQ = 0, NEGATIVE_SEQ, ZERO_SEQ };
 
+enum SyncMachineModel { SM_MODEL_1 = 0, SM_MODEL_2, SM_MODEL_3, SM_MODEL_4, SM_MODEL_5 };
+
 struct ReactiveLimits {
     double maxLimit = 0.0;
     double minLimit = 0.0;
@@ -45,7 +48,7 @@ struct ReactiveLimits {
  */
 class ElectricCalculation
 {
-public:
+   public:
     /**
      * @brief Constructor.
      */
@@ -57,7 +60,7 @@ public:
     ~ElectricCalculation();
 
     /**
-     * @brief Separate the elements from a generic list.
+     * @brief Separate the power elements from a generic list.
      * @param elementList List of generic elements.
      */
     virtual void GetElementsFromList(std::vector<Element*> elementList);
@@ -71,9 +74,10 @@ public:
      * @return Return true if was possible to build the admittance matrix.
      */
     virtual bool GetYBus(std::vector<std::vector<std::complex<double> > >& yBus,
-        double systemPowerBase,
-        YBusSequence sequence = POSITIVE_SEQ,
-        bool includeSyncMachines = false);
+                         double systemPowerBase,
+                         YBusSequence sequence = POSITIVE_SEQ,
+                         bool includeSyncMachines = false,
+                         bool allLoadsAsImpedances = false);
 
     /**
      * @brief Invert a matrix.
@@ -82,7 +86,7 @@ public:
      * @return Return true if was possible to invert the matrix.
      */
     virtual bool InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix,
-        std::vector<std::vector<std::complex<double> > >& inverse);
+                              std::vector<std::vector<std::complex<double> > >& inverse);
 
     /**
      * @brief Update the elements after the power flow calculation.
@@ -93,66 +97,77 @@ public:
      * @param systemPowerBase Base power of the system.
      */
     virtual void UpdateElementsPowerFlow(std::vector<std::complex<double> > voltage,
-        std::vector<std::complex<double> > power,
-        std::vector<BusType> busType,
-        std::vector<ReactiveLimits> reactiveLimit,
-        double systemPowerBase);
+                                         std::vector<std::complex<double> > power,
+                                         std::vector<BusType> busType,
+                                         std::vector<ReactiveLimits> reactiveLimit,
+                                         double systemPowerBase);
+
+    void ABCtoDQ0(std::complex<double> complexValue, double angle, double& dValue, double& qValue);
+    void DQ0toABC(double dValue, double qValue, double angle, std::complex<double>& complexValue);
+
+    std::vector<std::complex<double> > GaussianElimination(std::vector<std::vector<std::complex<double> > > matrix,
+                                                           std::vector<std::complex<double> > array);
+
+    SyncMachineModel GetMachineModel(SyncGenerator* generator);
+
+    std::vector<std::complex<double> > ComplexMatrixTimesVector(std::vector<std::vector<std::complex<double> > > matrix,
+                                                                std::vector<std::complex<double> > vector);
+
+    void 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);
+
+    std::vector<std::complex<double> > LUEvaluate(std::vector<std::vector<std::complex<double> > > u,
+                                                  std::vector<std::vector<std::complex<double> > > l,
+                                                  std::vector<std::complex<double> > b);
 
     /**
      * @brief Get the buses of the system (use GetElementsFromList first).
      * @return A list of bus elements.
      */
     const std::vector<Bus*> GetBusList() const { return m_busList; }
-
     /**
      * @brief Get the capacitors of the system (use GetElementsFromList first).
      * @return A list of capacitor elements.
      */
     const std::vector<Capacitor*> GetCapacitorList() const { return m_capacitorList; }
-    
     /**
      * @brief Get the induction motors of the system (use GetElementsFromList first).
      * @return A list of induction motor elements.
      */
     const std::vector<IndMotor*> GetIndMotorList() const { return m_indMotorList; }
-    
     /**
      * @brief Get the inductors of the system (use GetElementsFromList first).
      * @return A list of inductor elements.
      */
     const std::vector<Inductor*> GetInductorList() const { return m_inductorList; }
-    
     /**
      * @brief Get the lines of the system (use GetElementsFromList first).
      * @return A list of line elements.
      */
     const std::vector<Line*> GetLineList() const { return m_lineList; }
-    
     /**
      * @brief Get the loads of the system (use GetElementsFromList first).
      * @return A list of load elements.
      */
     const std::vector<Load*> GetLoadList() const { return m_loadList; }
-    
     /**
      * @brief Get the synchronous generators of the system (use GetElementsFromList first).
      * @return A list of synchronous generator elements.
      */
     const std::vector<SyncGenerator*> GetSyncGeneratorList() const { return m_syncGeneratorList; }
-    
     /**
      * @brief Get the synchronous motors of the system (use GetElementsFromList first).
      * @return A list of synchronous motor elements.
      */
     const std::vector<SyncMotor*> GetSyncMotorList() const { return m_syncMotorList; }
-    
     /**
      * @brief Get the transformers of the system (use GetElementsFromList first).
      * @return A list of transformer elements.
      */
     const std::vector<Transformer*> GetTransformerList() const { return m_transformerList; }
-
-protected:
+   protected:
+    std::vector<PowerElement*> m_powerElementList;
     std::vector<Bus*> m_busList;
     std::vector<Capacitor*> m_capacitorList;
     std::vector<IndMotor*> m_indMotorList;
@@ -164,4 +179,4 @@ protected:
     std::vector<Transformer*> m_transformerList;
 };
 
-#endif // ELECTRICCALCULATION_H
+#endif  // ELECTRICCALCULATION_H
diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp
new file mode 100644
index 0000000..85e10b4
--- /dev/null
+++ b/Project/Electromechanical.cpp
@@ -0,0 +1,59 @@
+#include "Electromechanical.h"
+
+Electromechanical::Electromechanical(std::vector<Element*> elementList)
+{
+    GetElementsFromList(elementList);
+    SetEventTimeList();
+}
+
+Electromechanical::~Electromechanical() {}
+void Electromechanical::SetEventTimeList()
+{
+    // Fault
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetEletricalData();
+        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((m_eventTimeList[i] - m_timeStep) < currentTime && m_eventTimeList[i] >= currentTime) {
+                m_eventOccurrenceList[i] = true;
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+void Electromechanical::SetEvent(double currentTime) {}
+bool Electromechanical::RunStabilityCalculation()
+{
+    // test
+    double simTime = 10.0;
+    double currentTime = 0.0;
+    while(currentTime <= simTime) {
+        if(HasEvent(currentTime)) {
+        }
+        currentTime += m_timeStep;
+    }
+    return true;
+}
diff --git a/Project/Electromechanical.h b/Project/Electromechanical.h
new file mode 100644
index 0000000..c8555b7
--- /dev/null
+++ b/Project/Electromechanical.h
@@ -0,0 +1,32 @@
+#ifndef ELECTROMECHANICAL_H
+#define ELECTROMECHANICAL_H
+
+#include "ElectricCalculation.h"
+
+class Electromechanical : public ElectricCalculation
+{
+public:
+    Electromechanical(std::vector<Element*> elementList);
+    ~Electromechanical();
+    
+    bool RunStabilityCalculation();
+    wxString GetErrorMessage() const { return m_errorMsg; }
+    
+protected:
+    void SetEventTimeList();
+    bool HasEvent(double currentTime);
+    void SetEvent(double currentTime);
+    
+    void Insert
+    
+    wxString m_errorMsg = _("Unknown error");
+    
+    std::vector<std::vector<std::complex<double> > > m_yBus;
+    
+    double m_timeStep = 1e-3;
+    
+    std::vector<double> m_eventTimeList;
+    std::vector<bool> m_eventOccurrenceList;
+};
+
+#endif // ELECTROMECHANICAL_H
diff --git a/Project/Machines.h b/Project/Machines.h
index 8efcab6..615d82e 100644
--- a/Project/Machines.h
+++ b/Project/Machines.h
@@ -5,10 +5,10 @@
 
 class Machines : public PowerElement
 {
-public:
-	Machines();
-	~Machines();
-	
+   public:
+    Machines();
+    ~Machines();
+
     virtual bool AddParent(Element* parent, wxPoint2DDouble position);
     virtual bool Contains(wxPoint2DDouble position) const { return m_rect.Contains(position); }
     virtual void Draw(wxPoint2DDouble translation, double scale) const;
@@ -22,14 +22,13 @@ public:
     virtual bool SetNodeParent(Element* parent);
     virtual void UpdateNodes();
     virtual void Rotate(bool clockwise = true);
-	virtual void DrawSymbol() const {}
+    virtual void DrawSymbol() const {}
     virtual void SetPowerFlowDirection(PowerFlowDirection pfDirection);
 
    protected:
     void UpdateSwitchesPosition();
     void UpdatePowerFlowArrowsPosition();
     bool m_inserted = false;
-
 };
 
-#endif // MACHINES_H
+#endif  // MACHINES_H
diff --git a/Project/MainFrame.cpp b/Project/MainFrame.cpp
index 87c547e..df9fe52 100644
--- a/Project/MainFrame.cpp
+++ b/Project/MainFrame.cpp
@@ -247,7 +247,14 @@ void MainFrame::OnPowerFlowClick(wxRibbonButtonBarEvent& event)
 void MainFrame::OnProjectSettingsClick(wxRibbonButtonBarEvent& event) {}
 void MainFrame::OnRedoClick(wxRibbonButtonBarEvent& event) {}
 void MainFrame::OnResetVoltagesClick(wxRibbonButtonBarEvent& event) {}
-void MainFrame::OnRunStabilityClick(wxRibbonButtonBarEvent& event) {}
+void MainFrame::OnRunStabilityClick(wxRibbonButtonBarEvent& event)
+{
+    Workspace* workspace = static_cast<Workspace*>(m_auiNotebook->GetCurrentPage());
+    if(workspace) {
+        workspace->RunStability();
+    }
+}
+
 void MainFrame::OnSCPowerClick(wxRibbonButtonBarEvent& event)
 {
     Workspace* workspace = static_cast<Workspace*>(m_auiNotebook->GetCurrentPage());
diff --git a/Project/Project.mk b/Project/Project.mk
index cea75af..7f113a8 100644
--- a/Project/Project.mk
+++ b/Project/Project.mk
@@ -13,7 +13,7 @@ CurrentFileName        :=
 CurrentFilePath        :=
 CurrentFileFullPath    :=
 User                   :=NDSE-69
-Date                   :=18/05/2017
+Date                   :=20/05/2017
 CodeLitePath           :="C:/Program Files/CodeLite"
 LinkerName             :=C:/TDM-GCC-64/bin/g++.exe
 SharedObjectLinkerName :=C:/TDM-GCC-64/bin/g++.exe -shared -fPIC
@@ -68,10 +68,10 @@ Objects0=$(IntermediateDirectory)/main.cpp$(ObjectSuffix) $(IntermediateDirector
 	$(IntermediateDirectory)/ControlEditor.cpp$(ObjectSuffix) $(IntermediateDirectory)/Camera.cpp$(ObjectSuffix) $(IntermediateDirectory)/ChartView.cpp$(ObjectSuffix) $(IntermediateDirectory)/MainFrameBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/WorkspaceBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/BusFormBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/ElementFormBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlEditorBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/ChartViewBitmaps.cpp$(ObjectSuffix) $(IntermediateDirectory)/MainFrameBase.cpp$(ObjectSuffix) \
 	$(IntermediateDirectory)/WorkspaceBase.cpp$(ObjectSuffix) $(IntermediateDirectory)/ElementForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlEditorBase.cpp$(ObjectSuffix) $(IntermediateDirectory)/ChartViewBase.cpp$(ObjectSuffix) $(IntermediateDirectory)/Bus.cpp$(ObjectSuffix) $(IntermediateDirectory)/Line.cpp$(ObjectSuffix) $(IntermediateDirectory)/Transformer.cpp$(ObjectSuffix) $(IntermediateDirectory)/Machines.cpp$(ObjectSuffix) $(IntermediateDirectory)/SyncGenerator.cpp$(ObjectSuffix) $(IntermediateDirectory)/IndMotor.cpp$(ObjectSuffix) \
 	$(IntermediateDirectory)/Branch.cpp$(ObjectSuffix) $(IntermediateDirectory)/SyncMotor.cpp$(ObjectSuffix) $(IntermediateDirectory)/Shunt.cpp$(ObjectSuffix) $(IntermediateDirectory)/Load.cpp$(ObjectSuffix) $(IntermediateDirectory)/Inductor.cpp$(ObjectSuffix) $(IntermediateDirectory)/Capacitor.cpp$(ObjectSuffix) $(IntermediateDirectory)/PowerElement.cpp$(ObjectSuffix) $(IntermediateDirectory)/ElectricCalculation.cpp$(ObjectSuffix) $(IntermediateDirectory)/PowerFlow.cpp$(ObjectSuffix) $(IntermediateDirectory)/Fault.cpp$(ObjectSuffix) \
-	$(IntermediateDirectory)/Text.cpp$(ObjectSuffix) $(IntermediateDirectory)/GraphicalElement.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElement.cpp$(ObjectSuffix) $(IntermediateDirectory)/TransferFunction.cpp$(ObjectSuffix) $(IntermediateDirectory)/ConnectionLine.cpp$(ObjectSuffix) $(IntermediateDirectory)/Sum.cpp$(ObjectSuffix) $(IntermediateDirectory)/Multiplier.cpp$(ObjectSuffix) $(IntermediateDirectory)/Limiter.cpp$(ObjectSuffix) $(IntermediateDirectory)/RateLimiter.cpp$(ObjectSuffix) $(IntermediateDirectory)/Exponential.cpp$(ObjectSuffix) \
-	$(IntermediateDirectory)/Constant.cpp$(ObjectSuffix) $(IntermediateDirectory)/Gain.cpp$(ObjectSuffix) $(IntermediateDirectory)/IOControl.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElementContainer.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElementSolver.cpp$(ObjectSuffix) $(IntermediateDirectory)/wxMathPlot_mathplot.cpp$(ObjectSuffix) $(IntermediateDirectory)/BusForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/GeneratorStabForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LineForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/SwitchingForm.cpp$(ObjectSuffix) \
-	$(IntermediateDirectory)/TransformerForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LoadForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ReactiveShuntElementForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/IndMotorForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/SyncMachineForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/TextForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/TransferFunctionForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/SumForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LimiterForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/RateLimiterForm.cpp$(ObjectSuffix) \
-	$(IntermediateDirectory)/ExponentialForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ConstantForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/GainForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/IOControlForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlSystemTest.cpp$(ObjectSuffix) 
+	$(IntermediateDirectory)/Electromechanical.cpp$(ObjectSuffix) $(IntermediateDirectory)/Text.cpp$(ObjectSuffix) $(IntermediateDirectory)/GraphicalElement.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElement.cpp$(ObjectSuffix) $(IntermediateDirectory)/TransferFunction.cpp$(ObjectSuffix) $(IntermediateDirectory)/ConnectionLine.cpp$(ObjectSuffix) $(IntermediateDirectory)/Sum.cpp$(ObjectSuffix) $(IntermediateDirectory)/Multiplier.cpp$(ObjectSuffix) $(IntermediateDirectory)/Limiter.cpp$(ObjectSuffix) $(IntermediateDirectory)/RateLimiter.cpp$(ObjectSuffix) \
+	$(IntermediateDirectory)/Exponential.cpp$(ObjectSuffix) $(IntermediateDirectory)/Constant.cpp$(ObjectSuffix) $(IntermediateDirectory)/Gain.cpp$(ObjectSuffix) $(IntermediateDirectory)/IOControl.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElementContainer.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlElementSolver.cpp$(ObjectSuffix) $(IntermediateDirectory)/wxMathPlot_mathplot.cpp$(ObjectSuffix) $(IntermediateDirectory)/BusForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/GeneratorStabForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LineForm.cpp$(ObjectSuffix) \
+	$(IntermediateDirectory)/SwitchingForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/TransformerForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LoadForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ReactiveShuntElementForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/IndMotorForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/SyncMachineForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/TextForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/TransferFunctionForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/SumForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/LimiterForm.cpp$(ObjectSuffix) \
+	$(IntermediateDirectory)/RateLimiterForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ExponentialForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ConstantForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/GainForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/IOControlForm.cpp$(ObjectSuffix) $(IntermediateDirectory)/ControlSystemTest.cpp$(ObjectSuffix) 
 
 
 
@@ -416,6 +416,14 @@ $(IntermediateDirectory)/Fault.cpp$(DependSuffix): Fault.cpp
 $(IntermediateDirectory)/Fault.cpp$(PreprocessSuffix): Fault.cpp
 	$(CXX) $(CXXFLAGS) $(IncludePCH) $(IncludePath) $(PreprocessOnlySwitch) $(OutputSwitch) $(IntermediateDirectory)/Fault.cpp$(PreprocessSuffix) Fault.cpp
 
+$(IntermediateDirectory)/Electromechanical.cpp$(ObjectSuffix): Electromechanical.cpp $(IntermediateDirectory)/Electromechanical.cpp$(DependSuffix)
+	$(CXX) $(IncludePCH) $(SourceSwitch) "C:/Users/NDSE-69/Documents/GitHub/PSP/Project/Electromechanical.cpp" $(CXXFLAGS) $(ObjectSwitch)$(IntermediateDirectory)/Electromechanical.cpp$(ObjectSuffix) $(IncludePath)
+$(IntermediateDirectory)/Electromechanical.cpp$(DependSuffix): Electromechanical.cpp
+	@$(CXX) $(CXXFLAGS) $(IncludePCH) $(IncludePath) -MG -MP -MT$(IntermediateDirectory)/Electromechanical.cpp$(ObjectSuffix) -MF$(IntermediateDirectory)/Electromechanical.cpp$(DependSuffix) -MM Electromechanical.cpp
+
+$(IntermediateDirectory)/Electromechanical.cpp$(PreprocessSuffix): Electromechanical.cpp
+	$(CXX) $(CXXFLAGS) $(IncludePCH) $(IncludePath) $(PreprocessOnlySwitch) $(OutputSwitch) $(IntermediateDirectory)/Electromechanical.cpp$(PreprocessSuffix) Electromechanical.cpp
+
 $(IntermediateDirectory)/Text.cpp$(ObjectSuffix): Text.cpp $(IntermediateDirectory)/Text.cpp$(DependSuffix)
 	$(CXX) $(IncludePCH) $(SourceSwitch) "C:/Users/NDSE-69/Documents/GitHub/PSP/Project/Text.cpp" $(CXXFLAGS) $(ObjectSwitch)$(IntermediateDirectory)/Text.cpp$(ObjectSuffix) $(IncludePath)
 $(IntermediateDirectory)/Text.cpp$(DependSuffix): Text.cpp
diff --git a/Project/Project.project b/Project/Project.project
index 5be2acf..2fe870c 100644
--- a/Project/Project.project
+++ b/Project/Project.project
@@ -29,6 +29,7 @@
         <File Name="ElectricCalculation.cpp"/>
         <File Name="PowerFlow.cpp"/>
         <File Name="Fault.cpp"/>
+        <File Name="Electromechanical.cpp"/>
       </VirtualDirectory>
       <File Name="ElementDataObject.cpp"/>
       <VirtualDirectory Name="graphical element">
@@ -116,6 +117,7 @@
         <File Name="ElectricCalculation.h"/>
         <File Name="PowerFlow.h"/>
         <File Name="Fault.h"/>
+        <File Name="Electromechanical.h"/>
       </VirtualDirectory>
       <VirtualDirectory Name="rapidXML">
         <File Name="rapidXML/rapidxml.hpp"/>
diff --git a/Project/Project.txt b/Project/Project.txt
index 7fdcdf6..eeefc35 100644
--- a/Project/Project.txt
+++ b/Project/Project.txt
@@ -1 +1 @@
-./Release/main.cpp.o ./Release/win_resources.rc.o ./Release/ElementDataObject.cpp.o ./Release/Element.cpp.o ./Release/ElementPlotData.cpp.o ./Release/ArtMetro.cpp.o ./Release/wxGLString.cpp.o ./Release/MainFrame.cpp.o ./Release/Workspace.cpp.o ./Release/FileHanding.cpp.o ./Release/ControlEditor.cpp.o ./Release/Camera.cpp.o ./Release/ChartView.cpp.o ./Release/MainFrameBitmaps.cpp.o ./Release/WorkspaceBitmaps.cpp.o ./Release/BusFormBitmaps.cpp.o ./Release/ElementFormBitmaps.cpp.o ./Release/ControlEditorBitmaps.cpp.o ./Release/ChartViewBitmaps.cpp.o ./Release/MainFrameBase.cpp.o ./Release/WorkspaceBase.cpp.o ./Release/ElementForm.cpp.o ./Release/ControlEditorBase.cpp.o ./Release/ChartViewBase.cpp.o ./Release/Bus.cpp.o ./Release/Line.cpp.o ./Release/Transformer.cpp.o ./Release/Machines.cpp.o ./Release/SyncGenerator.cpp.o ./Release/IndMotor.cpp.o ./Release/Branch.cpp.o ./Release/SyncMotor.cpp.o ./Release/Shunt.cpp.o ./Release/Load.cpp.o ./Release/Inductor.cpp.o ./Release/Capacitor.cpp.o ./Release/PowerElement.cpp.o ./Release/ElectricCalculation.cpp.o ./Release/PowerFlow.cpp.o ./Release/Fault.cpp.o ./Release/Text.cpp.o ./Release/GraphicalElement.cpp.o ./Release/ControlElement.cpp.o ./Release/TransferFunction.cpp.o ./Release/ConnectionLine.cpp.o ./Release/Sum.cpp.o ./Release/Multiplier.cpp.o ./Release/Limiter.cpp.o ./Release/RateLimiter.cpp.o ./Release/Exponential.cpp.o ./Release/Constant.cpp.o ./Release/Gain.cpp.o ./Release/IOControl.cpp.o ./Release/ControlElementContainer.cpp.o ./Release/ControlElementSolver.cpp.o ./Release/wxMathPlot_mathplot.cpp.o ./Release/BusForm.cpp.o ./Release/GeneratorStabForm.cpp.o ./Release/LineForm.cpp.o ./Release/SwitchingForm.cpp.o ./Release/TransformerForm.cpp.o ./Release/LoadForm.cpp.o ./Release/ReactiveShuntElementForm.cpp.o ./Release/IndMotorForm.cpp.o ./Release/SyncMachineForm.cpp.o ./Release/TextForm.cpp.o ./Release/TransferFunctionForm.cpp.o ./Release/SumForm.cpp.o ./Release/LimiterForm.cpp.o ./Release/RateLimiterForm.cpp.o ./Release/ExponentialForm.cpp.o ./Release/ConstantForm.cpp.o ./Release/GainForm.cpp.o ./Release/IOControlForm.cpp.o ./Release/ControlSystemTest.cpp.o   
+./Release/main.cpp.o ./Release/win_resources.rc.o ./Release/ElementDataObject.cpp.o ./Release/Element.cpp.o ./Release/ElementPlotData.cpp.o ./Release/ArtMetro.cpp.o ./Release/wxGLString.cpp.o ./Release/MainFrame.cpp.o ./Release/Workspace.cpp.o ./Release/FileHanding.cpp.o ./Release/ControlEditor.cpp.o ./Release/Camera.cpp.o ./Release/ChartView.cpp.o ./Release/MainFrameBitmaps.cpp.o ./Release/WorkspaceBitmaps.cpp.o ./Release/BusFormBitmaps.cpp.o ./Release/ElementFormBitmaps.cpp.o ./Release/ControlEditorBitmaps.cpp.o ./Release/ChartViewBitmaps.cpp.o ./Release/MainFrameBase.cpp.o ./Release/WorkspaceBase.cpp.o ./Release/ElementForm.cpp.o ./Release/ControlEditorBase.cpp.o ./Release/ChartViewBase.cpp.o ./Release/Bus.cpp.o ./Release/Line.cpp.o ./Release/Transformer.cpp.o ./Release/Machines.cpp.o ./Release/SyncGenerator.cpp.o ./Release/IndMotor.cpp.o ./Release/Branch.cpp.o ./Release/SyncMotor.cpp.o ./Release/Shunt.cpp.o ./Release/Load.cpp.o ./Release/Inductor.cpp.o ./Release/Capacitor.cpp.o ./Release/PowerElement.cpp.o ./Release/ElectricCalculation.cpp.o ./Release/PowerFlow.cpp.o ./Release/Fault.cpp.o ./Release/Electromechanical.cpp.o ./Release/Text.cpp.o ./Release/GraphicalElement.cpp.o ./Release/ControlElement.cpp.o ./Release/TransferFunction.cpp.o ./Release/ConnectionLine.cpp.o ./Release/Sum.cpp.o ./Release/Multiplier.cpp.o ./Release/Limiter.cpp.o ./Release/RateLimiter.cpp.o ./Release/Exponential.cpp.o ./Release/Constant.cpp.o ./Release/Gain.cpp.o ./Release/IOControl.cpp.o ./Release/ControlElementContainer.cpp.o ./Release/ControlElementSolver.cpp.o ./Release/wxMathPlot_mathplot.cpp.o ./Release/BusForm.cpp.o ./Release/GeneratorStabForm.cpp.o ./Release/LineForm.cpp.o ./Release/SwitchingForm.cpp.o ./Release/TransformerForm.cpp.o ./Release/LoadForm.cpp.o ./Release/ReactiveShuntElementForm.cpp.o ./Release/IndMotorForm.cpp.o ./Release/SyncMachineForm.cpp.o ./Release/TextForm.cpp.o ./Release/TransferFunctionForm.cpp.o ./Release/SumForm.cpp.o ./Release/LimiterForm.cpp.o ./Release/RateLimiterForm.cpp.o ./Release/ExponentialForm.cpp.o ./Release/ConstantForm.cpp.o ./Release/GainForm.cpp.o ./Release/IOControlForm.cpp.o ./Release/ControlSystemTest.cpp.o   
diff --git a/Project/Workspace.cpp b/Project/Workspace.cpp
index 5575764..4b2cf29 100644
--- a/Project/Workspace.cpp
+++ b/Project/Workspace.cpp
@@ -16,14 +16,11 @@
 
 #include "PowerFlow.h"
 #include "Fault.h"
+#include "Electromechanical.h"
 
 // Workspace
-Workspace::Workspace()
-    : WorkspaceBase(NULL)
-{
-}
-Workspace::Workspace(wxWindow* parent, wxString name, wxStatusBar* statusBar)
-    : WorkspaceBase(parent)
+Workspace::Workspace() : WorkspaceBase(NULL) {}
+Workspace::Workspace(wxWindow* parent, wxString name, wxStatusBar* statusBar) : WorkspaceBase(parent)
 {
     m_timer->Start();
     m_name = name;
@@ -36,7 +33,7 @@ Workspace::Workspace(wxWindow* parent, wxString name, wxStatusBar* statusBar)
         m_elementNumber[i] = 1;
     }
 
-    const int widths[4] = { -3, -1, 100, 100 };
+    const int widths[4] = {-3, -1, 100, 100};
     m_statusBar->SetStatusWidths(4, widths);
 }
 
@@ -60,8 +57,8 @@ void Workspace::OnPaint(wxPaintEvent& event)
     SetViewport();
 
     // Set GLCanvas scale and translation.
-    glScaled(m_camera->GetScale(), m_camera->GetScale(), 0.0);                         // Scale
-    glTranslated(m_camera->GetTranslation().m_x, m_camera->GetTranslation().m_y, 0.0); // Translation
+    glScaled(m_camera->GetScale(), m_camera->GetScale(), 0.0);                          // Scale
+    glTranslated(m_camera->GetTranslation().m_x, m_camera->GetTranslation().m_y, 0.0);  // Translation
 
     // Draw
 
@@ -94,14 +91,14 @@ void Workspace::OnPaint(wxPaintEvent& event)
     glVertex2d(m_selectionRect.m_x + m_selectionRect.m_width, m_selectionRect.m_y);
     glEnd();
 
-    glFlush(); // Sends all pending information directly to the GPU.
+    glFlush();  // Sends all pending information directly to the GPU.
     m_glCanvas->SwapBuffers();
     event.Skip();
 }
 
 void Workspace::SetViewport()
 {
-    glClearColor(1.0, 1.0, 1.0, 1.0); // White background.
+    glClearColor(1.0, 1.0, 1.0, 1.0);  // White background.
     glClear(GL_COLOR_BUFFER_BIT);
     glDisable(GL_DEPTH_TEST);
     glDisable(GL_TEXTURE_2D);
@@ -143,7 +140,7 @@ void Workspace::OnLeftClickDown(wxMouseEvent& event)
                 if(typeid(*element) == typeid(Bus)) {
                     // Select the bus.
                     element->SetSelected();
-                    foundElement = true; // Element found.
+                    foundElement = true;  // Element found.
                     // Add the new element's parent. If the element being inserted returns true, back to
                     // edit mode.
                     if(newElement->AddParent(element, m_camera->ScreenToWorld(clickPoint))) {
@@ -166,7 +163,7 @@ void Workspace::OnLeftClickDown(wxMouseEvent& event)
         bool clickPickbox = false;
         for(auto it = m_elementList.begin(), itEnd = m_elementList.end(); it != itEnd; ++it) {
             PowerElement* element = *it;
-            element->ResetPickboxes(); // Reset pickbox state.
+            element->ResetPickboxes();  // Reset pickbox state.
 
             // Set movement initial position (not necessarily will be moved).
             element->StartMove(m_camera->ScreenToWorld(clickPoint));
@@ -260,14 +257,14 @@ void Workspace::OnLeftDoubleClick(wxMouseEvent& event)
             if(elementIsBus) {
                 // The voltage was changed
                 if(oldBus.GetEletricalData().nominalVoltage != currentBus->GetEletricalData().nominalVoltage ||
-                    oldBus.GetEletricalData().nominalVoltageUnit != currentBus->GetEletricalData().nominalVoltageUnit) {
+                   oldBus.GetEletricalData().nominalVoltageUnit != currentBus->GetEletricalData().nominalVoltageUnit) {
                     // Check if the bus has line as child.
                     std::vector<Element*> childList = element->GetChildList();
                     for(auto itc = childList.begin(), itcEnd = childList.end(); itc != itcEnd; ++itc) {
                         Element* child = *itc;
                         if(typeid(*child) == typeid(Line)) {
                             wxMessageDialog msgDialog(this, _("Do you want to change the rated voltage of the path?"),
-                                _("Warning"), wxYES_NO | wxCENTRE | wxICON_WARNING);
+                                                      _("Warning"), wxYES_NO | wxCENTRE | wxICON_WARNING);
                             if(msgDialog.ShowModal() == wxID_YES)
                                 ValidateBusesVoltages(element);
                             else {
@@ -445,7 +442,7 @@ void Workspace::OnMouseMotion(wxMouseEvent& event)
     bool redraw = false;
     switch(m_mode) {
         case MODE_INSERT: {
-            Element* newElement = *(m_elementList.end() - 1); // Get the last element in the list.
+            Element* newElement = *(m_elementList.end() - 1);  // Get the last element in the list.
             newElement->SetPosition(m_camera->ScreenToWorld(event.GetPosition()));
             redraw = true;
         } break;
@@ -563,7 +560,7 @@ void Workspace::OnMouseMotion(wxMouseEvent& event)
     if(redraw) Redraw();
     m_camera->UpdateMousePosition(event.GetPosition());
     UpdateStatusBar();
-    m_timer->Start(); // Restart the timer.
+    m_timer->Start();  // Restart the timer.
     event.Skip();
 }
 
@@ -631,10 +628,10 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
     char key = event.GetUnicodeKey();
     if(key != WXK_NONE) {
         switch(key) {
-            case WXK_ESCAPE: // Cancel operations.
+            case WXK_ESCAPE:  // Cancel operations.
             {
                 if(m_mode == MODE_INSERT) {
-                    m_elementList.pop_back(); // Removes the last element being inserted.
+                    m_elementList.pop_back();  // Removes the last element being inserted.
                     m_mode = MODE_EDIT;
                     Redraw();
                 } else if(m_mode == MODE_INSERT_TEXT) {
@@ -643,7 +640,7 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                     Redraw();
                 }
             } break;
-            case WXK_DELETE: // Delete selected elements
+            case WXK_DELETE:  // Delete selected elements
             {
                 DeleteSelectedElements();
             } break;
@@ -661,15 +658,15 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                     Fit();
                 }
             } break;
-            case 'R': // Rotate the selected elements.
+            case 'R':  // Rotate the selected elements.
             {
                 RotateSelectedElements(event.GetModifiers() != wxMOD_SHIFT);
             } break;
-            case 'B': // Insert a bus.
+            case 'B':  // Insert a bus.
             {
                 if(!insertingElement) {
                     Bus* newBus = new Bus(m_camera->ScreenToWorld(event.GetPosition()),
-                        wxString::Format(_("Bus %d"), GetElementNumber(ID_BUS)));
+                                          wxString::Format(_("Bus %d"), GetElementNumber(ID_BUS)));
                     IncrementElementNumber(ID_BUS);
                     m_elementList.push_back(newBus);
                     m_mode = MODE_INSERT;
@@ -679,13 +676,13 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
             } break;
             case 'L': {
                 if(!insertingElement) {
-                    if(!event.ControlDown() && event.ShiftDown()) { // Insert a load.
+                    if(!event.ControlDown() && event.ShiftDown()) {  // Insert a load.
                         Load* newLoad = new Load(wxString::Format(_("Load %d"), GetElementNumber(ID_LOAD)));
                         IncrementElementNumber(ID_LOAD);
                         m_elementList.push_back(newLoad);
                         m_mode = MODE_INSERT;
                         m_statusBar->SetStatusText(_("Insert Load: Click on a buses, ESC to cancel."));
-                    } else if(!event.ControlDown() && !event.ShiftDown()) { // Insert a power line.
+                    } else if(!event.ControlDown() && !event.ShiftDown()) {  // Insert a power line.
                         Line* newLine = new Line(wxString::Format(_("Line %d"), GetElementNumber(ID_LINE)));
                         IncrementElementNumber(ID_LINE);
                         m_elementList.push_back(newLine);
@@ -696,11 +693,12 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                 }
                 // Tests - Ctrl + Shift + L
                 if(event.ControlDown() && event.ShiftDown()) {
-                    ControlEditor* ce = new ControlEditor(this, IOControl::IN_TERMINAL_VOLTAGE | IOControl::OUT_FIELD_VOLTAGE);
+                    ControlEditor* ce =
+                        new ControlEditor(this, IOControl::IN_TERMINAL_VOLTAGE | IOControl::OUT_FIELD_VOLTAGE);
                     ce->Show();
                 }
             } break;
-            case 'T': // Insert a transformer.
+            case 'T':  // Insert a transformer.
             {
                 if(!insertingElement) {
                     Transformer* newTransformer =
@@ -712,7 +710,7 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                     Redraw();
                 }
             } break;
-            case 'G': // Insert a generator.
+            case 'G':  // Insert a generator.
             {
                 if(!insertingElement) {
                     SyncGenerator* newGenerator =
@@ -726,14 +724,14 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
             } break;
             case 'I': {
                 if(!insertingElement) {
-                    if(event.GetModifiers() == wxMOD_SHIFT) { // Insert an inductor.
+                    if(event.GetModifiers() == wxMOD_SHIFT) {  // Insert an inductor.
                         Inductor* newInductor =
                             new Inductor(wxString::Format(_("Inductor %d"), GetElementNumber(ID_INDUCTOR)));
                         IncrementElementNumber(ID_INDUCTOR);
                         m_elementList.push_back(newInductor);
                         m_mode = MODE_INSERT;
                         m_statusBar->SetStatusText(_("Insert Inductor: Click on a buses, ESC to cancel."));
-                    } else // Insert an induction motor.
+                    } else  // Insert an induction motor.
                     {
                         IndMotor* newIndMotor =
                             new IndMotor(wxString::Format(_("Induction motor %d"), GetElementNumber(ID_INDMOTOR)));
@@ -745,7 +743,7 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                     Redraw();
                 }
             } break;
-            case 'K': // Insert a synchronous condenser.
+            case 'K':  // Insert a synchronous condenser.
             {
                 if(!insertingElement) {
                     SyncMotor* newSyncCondenser =
@@ -759,7 +757,7 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
             } break;
             case 'C': {
                 if(!insertingElement) {
-                    if(event.GetModifiers() == wxMOD_SHIFT) { // Insert a capacitor.
+                    if(event.GetModifiers() == wxMOD_SHIFT) {  // Insert a capacitor.
                         Capacitor* newCapacitor =
                             new Capacitor(wxString::Format(_("Capacitor %d"), GetElementNumber(ID_CAPACITOR)));
                         IncrementElementNumber(ID_CAPACITOR);
@@ -767,7 +765,7 @@ void Workspace::OnKeyDown(wxKeyEvent& event)
                         m_mode = MODE_INSERT;
                         m_statusBar->SetStatusText(_("Insert Capacitor: Click on a buses, ESC to cancel."));
                         Redraw();
-                    } else if(event.GetModifiers() == wxMOD_CONTROL) { // Copy.
+                    } else if(event.GetModifiers() == wxMOD_CONTROL) {  // Copy.
                         CopySelection();
                     }
                 }
@@ -916,7 +914,7 @@ void Workspace::RotateSelectedElements(bool clockwise)
         // Parent's element rotating...
         for(int i = 0; i < (int)element->GetParentList().size(); i++) {
             Element* parent = element->GetParentList()[i];
-            if(parent) { // Check if parent is not null
+            if(parent) {  // Check if parent is not null
                 if(parent->IsSelected()) {
                     element->RotateNode(parent, clockwise);
                     // Update the positions used on motion action, the element will not be necessarily
@@ -949,7 +947,6 @@ void Workspace::DeleteSelectedElements()
         Element* element = *it;
 
         if(element->IsSelected()) {
-
             // Remove child/parent.
             std::vector<Element*> childList = element->GetChildList();
             for(auto itc = childList.begin(), itEnd = childList.end(); itc != itEnd; ++itc) {
@@ -992,8 +989,8 @@ void Workspace::DeleteSelectedElements()
 }
 
 bool Workspace::GetElementsCorners(wxPoint2DDouble& leftUpCorner,
-    wxPoint2DDouble& rightDownCorner,
-    std::vector<Element*> elementList)
+                                   wxPoint2DDouble& rightDownCorner,
+                                   std::vector<Element*> elementList)
 {
     if(elementList.size() == 0) return false;
 
@@ -1054,7 +1051,7 @@ void Workspace::ValidateBusesVoltages(Element* initialBus)
                 BusElectricalData data2 = static_cast<Bus*>(child->GetParentList()[1])->GetEletricalData();
 
                 if(data1.nominalVoltage != data2.nominalVoltage ||
-                    data1.nominalVoltageUnit != data2.nominalVoltageUnit) {
+                   data1.nominalVoltageUnit != data2.nominalVoltageUnit) {
                     data1.nominalVoltage = nominalVoltage;
                     data2.nominalVoltage = nominalVoltage;
                     data1.nominalVoltageUnit = nominalVoltageUnit;
@@ -1063,7 +1060,7 @@ void Workspace::ValidateBusesVoltages(Element* initialBus)
                     static_cast<Bus*>(child->GetParentList()[0])->SetElectricalData(data1);
                     static_cast<Bus*>(child->GetParentList()[1])->SetElectricalData(data2);
 
-                    it = m_elementList.begin(); // Restart search.
+                    it = m_elementList.begin();  // Restart search.
                 }
             }
         }
@@ -1151,7 +1148,7 @@ bool Workspace::Paste()
         if(wxTheClipboard->IsSupported(wxDataFormat("PSPCopy"))) {
             if(!wxTheClipboard->GetData(dataObject)) {
                 wxMessageDialog dialog(this, _("It was not possible to paste from clipboard."), _("Error"),
-                    wxOK | wxCENTER | wxICON_ERROR, wxDefaultPosition);
+                                       wxOK | wxCENTER | wxICON_ERROR, wxDefaultPosition);
                 dialog.ShowModal();
                 wxTheClipboard->Close();
                 return false;
@@ -1169,7 +1166,7 @@ bool Workspace::Paste()
 
         // Paste buses (parents).
         auto parentList = elementsLists->parentList;
-        std::vector<Bus*> pastedBusList; // To set new parents;
+        std::vector<Bus*> pastedBusList;  // To set new parents;
         for(auto it = parentList.begin(), itEnd = parentList.end(); it != itEnd; ++it) {
             Element* copy = (*it)->GetCopy();
             if(copy) {
@@ -1259,7 +1256,7 @@ bool Workspace::Paste()
         }
     } else {
         wxMessageDialog dialog(this, _("It was not possible to paste from clipboard."), _("Error"),
-            wxOK | wxCENTER | wxICON_ERROR, wxDefaultPosition);
+                               wxOK | wxCENTER | wxICON_ERROR, wxDefaultPosition);
         dialog.ShowModal();
         return false;
     }
@@ -1393,3 +1390,15 @@ bool Workspace::RunSCPower()
 
     return result;
 }
+
+bool Workspace::RunStability()
+{
+    Electromechanical stability(GetElementList());
+    bool result = stability.RunStabilityCalculation();
+    if(!result) {
+        wxMessageDialog msgDialog(this, stability.GetErrorMessage(), _("Error"), wxOK | wxCENTRE | wxICON_ERROR);
+        msgDialog.ShowModal();
+    }
+
+    return result;
+}
diff --git a/Project/Workspace.h b/Project/Workspace.h
index bf03a9e..f007b59 100644
--- a/Project/Workspace.h
+++ b/Project/Workspace.h
@@ -31,6 +31,7 @@ class Text;
 
 class PowerFlow;
 class Fault;
+class Electromechanical;
 
 enum ElementID {
     ID_BUS = 0,
@@ -110,6 +111,7 @@ public:
     bool RunPowerFlow();
     bool RunFault();
     bool RunSCPower();
+    bool RunStability();
 
 protected:
     virtual void OnIdle(wxIdleEvent& event);