Merge pull request #51 from Thales1330/wip/induction-motor

Newton bug fixed

1 files changed, 550 insertions, 82 deletions

diff --git a/Project/PowerFlow.cpp b/Project/PowerFlow.cpp
index bafbbfb..1b2ae13 100644
--- a/Project/PowerFlow.cpp
+++ b/Project/PowerFlow.cpp
@@ -20,11 +20,13 @@
 PowerFlow::PowerFlow() : ElectricCalculation() {}
 PowerFlow::PowerFlow(std::vector<Element*> elementList) : ElectricCalculation() { GetElementsFromList(elementList); }
 PowerFlow::~PowerFlow() {}
-bool PowerFlow::RunGaussSeidel(double systemPowerBase,
-                               int maxIteration,
-                               double error,
-                               double initAngle,
-                               double accFactor)
+bool PowerFlow::InitPowerFlow(std::vector<BusType>& busType,
+                              std::vector<std::complex<double> >& voltage,
+                              std::vector<std::complex<double> >& power,
+                              std::vector<std::complex<double> >& loadPower,
+                              std::vector<ReactiveLimits>& reactiveLimit,
+                              double systemPowerBase,
+                              double initAngle)
 {
     double radInitAngle = wxDegToRad(initAngle);
     // Calculate the Ybus.
@@ -33,16 +35,16 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
         return false;
     }
 
-    // Number of buses on the system.
-    int numberOfBuses = static_cast<int>(m_busList.size());
+    // Number of buses in the system.
+    m_numberOfBuses = static_cast<int>(m_busList.size());
 
-    std::vector<BusType> busType;                  // Bus type
-    std::vector<std::complex<double> > voltage;    // Voltage of buses
-    std::vector<std::complex<double> > power;      // Injected power
-    std::vector<std::complex<double> > loadPower;  // Only the load power
-    std::vector<ReactiveLimits> reactiveLimit;     // Limit of reactive power on PV buses
+    busType.clear();
+    voltage.clear();
+    power.clear();
+    loadPower.clear();
+    reactiveLimit.clear();
 
-    reactiveLimit.resize(numberOfBuses);
+    reactiveLimit.resize(m_numberOfBuses);
 
     int busNumber = 0;
     for(auto itb = m_busList.begin(); itb != m_busList.end(); itb++) {
@@ -165,7 +167,7 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
     // Check if have slack bus and if have generation on the slack bus
     bool haveSlackBus = false;
     bool slackBusHaveGeneration = false;
-    for(int i = 0; i < (int)busType.size(); i++) {
+    for(unsigned int i = 0; i < busType.size(); i++) {
         if(busType[i] == BUS_SLACK) {
             auto itb = m_busList.begin();
             std::advance(itb, i);
@@ -187,6 +189,23 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
         return false;
     }
 
+    return true;
+}
+
+bool PowerFlow::RunGaussSeidel(double systemPowerBase,
+                               int maxIteration,
+                               double error,
+                               double initAngle,
+                               double accFactor)
+{
+    std::vector<BusType> busType;                  // Bus type
+    std::vector<std::complex<double> > voltage;    // Voltage of buses
+    std::vector<std::complex<double> > power;      // Injected power
+    std::vector<std::complex<double> > loadPower;  // Only the load power
+    std::vector<ReactiveLimits> reactiveLimit;     // Limit of reactive power on PV buses
+
+    if(!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
+
     // Gauss-Seidel method
     std::vector<std::complex<double> > oldVoltage;  // Old voltage array.
     oldVoltage.resize(voltage.size());
@@ -203,98 +222,235 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
         }
 
         // Update the old voltage array to current iteration values.
-        for(int i = 0; i < numberOfBuses; i++) oldVoltage[i] = voltage[i];
+        for(int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
 
-        double iterationError = 0.0;
+        double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
 
-        for(int i = 0; i < numberOfBuses; i++) {
-            if(busType[i] == BUS_PQ) {
-                std::complex<double> yeSum(0.0, 0.0);
-                for(int k = 0; k < numberOfBuses; k++) {
-                    if(i != k) {
-                        // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
-                        yeSum += m_yBus[i][k] * voltage[k];
-                    }
-                }
+        if(iterationError < error) {
+            if(!CheckReactiveLimits(busType, reactiveLimit, power, loadPower)) break;
+        }
 
-                // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
-                std::complex<double> newVolt =
-                    (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
+        iteration++;
+    }
 
-                // Apply the acceleration factor.
-                newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
-                                               accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+    // Adjust the power array.
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        std::complex<double> sBus = std::complex<double>(0.0, 0.0);
+        for(int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
+        power[i] = sBus;
+    }
 
-                voltage[i] = newVolt;
-            }
-            if(busType[i] == BUS_PV) {
-                std::complex<double> yeSum(0.0, 0.0);
-                for(int k = 0; k < numberOfBuses; k++) {
-                    if(i != k) {
-                        // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
-                        yeSum += m_yBus[i][k] * voltage[k];
-                    }
+    UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
+
+    return true;
+}
+
+bool PowerFlow::RunNewtonRaphson(double systemPowerBase, int maxIteration, double error, double initAngle)
+{
+    std::vector<BusType> busType;                  // Bus type
+    std::vector<std::complex<double> > voltage;    // Voltage of buses
+    std::vector<std::complex<double> > power;      // Injected power
+    std::vector<std::complex<double> > loadPower;  // Only the load power
+    std::vector<ReactiveLimits> reactiveLimit;     // Limit of reactive power on PV buses
+
+    if(!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
+    auto oldBusType = busType;
+
+    // Newton-Raphson method
+    int numPQ = 0;  // Number of PQ buses
+    int numPV = 0;  // Number of PV buses
+    GetNumPVPQ(busType, numPQ, numPV);
+
+    // DeltaP and DeltaQ array
+    std::vector<double> dPdQ;
+    dPdQ.resize(numPV + 2 * numPQ, 0.0);
+
+    int iteration = 0;  // Current iteration number.
+    while(true) {
+        // Reach the max number of iterations.
+        if(iteration >= maxIteration) {
+            m_errorMsg = _("The maximum number of iterations was reached.");
+            return false;
+        }
+
+        // Calculate dPdQ array
+
+        // Fill it with zeros
+        std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
+
+        int indexDP = 0;
+        int indexDQ = numPQ + numPV;
+        for(int i = 0; i < m_numberOfBuses; i++) {
+            if(busType[i] != BUS_SLACK) {
+                for(int j = 0; j < m_numberOfBuses; j++) {
+                    // PV ou PQ bus
+                    std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
+                    dPdQ[indexDP] += sInj.real();
+
+                    // PQ bus
+                    if(busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
                 }
-                std::complex<double> yeSumT = yeSum + (m_yBus[i][i] * voltage[i]);
 
-                // Q[i] = - Im( E*[i] * Sum { Y[i,k] * E[k] } )
-                std::complex<double> qCalc = std::conj(voltage[i]) * yeSumT;
-                power[i] = std::complex<double>(power[i].real(), -qCalc.imag());
+                // PQ or PV bus
+                dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
+                indexDP++;
 
-                // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
-                std::complex<double> newVolt =
-                    (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
+                // PQ bus
+                if(busType[i] == BUS_PQ) {
+                    dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
+                    indexDQ++;
+                }
+            }
+        }
 
-                // Apply the acceleration factor.
-                newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
-                                               accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+        // Calculate the iteration error
+        double iterationError = 0.0;
+        for(unsigned int i = 0; i < dPdQ.size(); ++i) {
+            if(iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
+        }
 
-                // Keep the same voltage magnitude.
-                voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)),
-                                                  std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
+        // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
+        // If any reactive limit was reached, change the bus type.
+        if(iterationError < error) {
+            if(!CheckReactiveLimits(busType, reactiveLimit, power, loadPower))
+                break;
+            else {
+                GetNumPVPQ(busType, numPQ, numPV);
+                dPdQ.clear();
+                dPdQ.resize(numPV + 2 * numPQ, 0.0);
             }
+        }
+
+        NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ);
+
+        iteration++;
+    }
+
+    // Adjust the power array.
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        std::complex<double> sBus = std::complex<double>(0.0, 0.0);
+        for(int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
+        power[i] = sBus;
+    }
+
+    UpdateElementsPowerFlow(voltage, power, oldBusType, reactiveLimit, systemPowerBase);
+
+    return true;
+}
+
+bool PowerFlow::RunGaussNewton(double systemPowerBase,
+                               int maxIteration,
+                               double error,
+                               double initAngle,
+                               double accFactor,
+                               double gaussTol)
+{
+    std::vector<BusType> busType;                  // Bus type
+    std::vector<std::complex<double> > voltage;    // Voltage of buses
+    std::vector<std::complex<double> > power;      // Injected power
+    std::vector<std::complex<double> > loadPower;  // Only the load power
+    std::vector<ReactiveLimits> reactiveLimit;     // Limit of reactive power on PV buses
+
+    if(!InitPowerFlow(busType, voltage, power, loadPower, reactiveLimit, systemPowerBase, initAngle)) return false;
+
+    // Gauss-Seidel method
+    std::vector<std::complex<double> > oldVoltage;  // Old voltage array.
+    oldVoltage.resize(voltage.size());
+
+    auto oldBusType = busType;
+
+    int iteration = 0;  // Current itaration number.
+
+    while(true) {
+        // Reach the max number of iterations.
+        if(iteration >= maxIteration) {
+            m_errorMsg = _("The maximum number of iterations was reached.");
+            return false;
+        }
+
+        // Update the old voltage array to current iteration values.
+        for(int i = 0; i < m_numberOfBuses; i++) oldVoltage[i] = voltage[i];
+
+        double iterationError = GaussSeidel(busType, voltage, oldVoltage, power, accFactor);
 
-            double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()),
-                                       std::abs(voltage[i].imag() - oldVoltage[i].imag()));
+        if(iterationError < gaussTol) break;
+
+        iteration++;
+    }
 
-            if(busError > iterationError) iterationError = busError;
+    // Newton-Raphson method
+    int numPQ = 0;  // Number of PQ buses
+    int numPV = 0;  // Number of PV buses
+    GetNumPVPQ(busType, numPQ, numPV);
+
+    // DeltaP and DeltaQ array
+    std::vector<double> dPdQ;
+    dPdQ.resize(numPV + 2 * numPQ, 0.0);
+
+    while(true) {
+        // Reach the max number of iterations.
+        if(iteration >= maxIteration) {
+            m_errorMsg = _("The maximum number of iterations was reached.");
+            return false;
         }
 
-        if(iterationError < error) {
-            bool limitReach = false;
-            for(int i = 0; i < numberOfBuses; i++) {
-                if(busType[i] == BUS_PV) {
-                    if(reactiveLimit[i].maxLimitType == RL_LIMITED) {
-                        if(power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) {
-                            power[i] =
-                                std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
-                            busType[i] = BUS_PQ;
-                            reactiveLimit[i].limitReached = RL_MAX_REACHED;
-                            limitReach = true;
-                        }
-                    }
-                    if(reactiveLimit[i].minLimitType == RL_LIMITED) {
-                        if(power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) {
-                            power[i] =
-                                std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
-                            busType[i] = BUS_PQ;
-                            reactiveLimit[i].limitReached = RL_MIN_REACHED;
-                            limitReach = true;
-                        }
-                    }
+        // Calculate dPdQ array
+
+        // Fill it with zeros
+        std::fill(dPdQ.begin(), dPdQ.end(), 0.0);
+
+        int indexDP = 0;
+        int indexDQ = numPQ + numPV;
+        for(int i = 0; i < m_numberOfBuses; i++) {
+            if(busType[i] != BUS_SLACK) {
+                for(int j = 0; j < m_numberOfBuses; j++) {
+                    // PV ou PQ bus
+                    std::complex<double> sInj = std::conj(m_yBus[i][j]) * voltage[i] * std::conj(voltage[j]);
+                    dPdQ[indexDP] += sInj.real();
+
+                    // PQ bus
+                    if(busType[i] == BUS_PQ) dPdQ[indexDQ] += sInj.imag();
+                }
+
+                // PQ or PV bus
+                dPdQ[indexDP] = power[i].real() - dPdQ[indexDP];
+                indexDP++;
+
+                // PQ bus
+                if(busType[i] == BUS_PQ) {
+                    dPdQ[indexDQ] = power[i].imag() - dPdQ[indexDQ];
+                    indexDQ++;
                 }
             }
-            if(!limitReach) break;
         }
 
+        // Calculate the iteration error
+        double iterationError = 0.0;
+        for(unsigned int i = 0; i < dPdQ.size(); ++i) {
+            if(iterationError < std::abs(dPdQ[i])) iterationError = std::abs(dPdQ[i]);
+        }
+
+        // Check if the iteration error is less than tolerance, also check if any reactive limit was reached.
+        // If any reactive limit was reached, change the bus type.
+        if(iterationError < error) {
+            if(!CheckReactiveLimits(busType, reactiveLimit, power, loadPower))
+                break;
+            else {
+                GetNumPVPQ(busType, numPQ, numPV);
+                dPdQ.clear();
+                dPdQ.resize(numPV + 2 * numPQ, 0.0);
+            }
+        }
+
+        NewtonRaphson(busType, voltage, power, numPV, numPQ, dPdQ);
+
         iteration++;
     }
 
     // Adjust the power array.
-    // TODO: Only the slack bus??
-    for(int i = 0; i < numberOfBuses; i++) {
+    for(int i = 0; i < m_numberOfBuses; i++) {
         std::complex<double> sBus = std::complex<double>(0.0, 0.0);
-        for(int j = 0; j < numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
+        for(int j = 0; j < m_numberOfBuses; j++) sBus += voltage[i] * std::conj(voltage[j]) * std::conj(m_yBus[i][j]);
         power[i] = sBus;
     }
 
@@ -302,3 +458,315 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
 
     return true;
 }
+
+void PowerFlow::GetNumPVPQ(std::vector<BusType> busType, int& numPQ, int& numPV)
+{
+    numPQ = 0;
+    numPV = 0;
+    for(auto it = busType.begin(), itEnd = busType.end(); it != itEnd; ++it) {
+        if(*it == BUS_PQ)
+            numPQ++;
+        else if(*it == BUS_PV)
+            numPV++;
+    }
+}
+
+double PowerFlow::GaussSeidel(std::vector<BusType> busType,
+                              std::vector<std::complex<double> >& voltage,
+                              std::vector<std::complex<double> > oldVoltage,
+                              std::vector<std::complex<double> >& power,
+                              double accFactor)
+{
+    double error = 0.0;
+
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] == BUS_PQ) {
+            std::complex<double> yeSum(0.0, 0.0);
+            for(int k = 0; k < m_numberOfBuses; k++) {
+                if(i != k) {
+                    // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
+                    yeSum += m_yBus[i][k] * voltage[k];
+                }
+            }
+
+            // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
+            std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
+
+            // Apply the acceleration factor.
+            newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
+                                           accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+
+            voltage[i] = newVolt;
+        }
+        if(busType[i] == BUS_PV) {
+            std::complex<double> yeSum(0.0, 0.0);
+            for(int k = 0; k < m_numberOfBuses; k++) {
+                if(i != k) {
+                    // Sum { Y[i,k] * E[k] } | k = 1->n; k diff i
+                    yeSum += m_yBus[i][k] * voltage[k];
+                }
+            }
+            std::complex<double> yeSumT = yeSum + (m_yBus[i][i] * voltage[i]);
+
+            // Q[i] = - Im( E*[i] * Sum { Y[i,k] * E[k] } )
+            std::complex<double> qCalc = std::conj(voltage[i]) * yeSumT;
+            power[i] = std::complex<double>(power[i].real(), -qCalc.imag());
+
+            // E[i] = (1/Y[i,i])*((P[i]-jQ[i])/E*[i] - Sum { Y[i,k] * E[k] (k diff i) })
+            std::complex<double> newVolt = (1.0 / m_yBus[i][i]) * (std::conj(power[i]) / std::conj(voltage[i]) - yeSum);
+
+            // Apply the acceleration factor.
+            newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
+                                           accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+
+            // Keep the same voltage magnitude.
+            voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)),
+                                              std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
+        }
+
+        double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()),
+                                   std::abs(voltage[i].imag() - oldVoltage[i].imag()));
+
+        if(busError > error) error = busError;
+    }
+    return error;
+}
+
+void PowerFlow::NewtonRaphson(std::vector<BusType> busType,
+                              std::vector<std::complex<double> >& voltage,
+                              std::vector<std::complex<double> > power,
+                              int numPV,
+                              int numPQ,
+                              std::vector<double> dPdQ)
+{
+    // Jacobian matrix
+    std::vector<std::vector<double> > jacobMatrix = CalculateJacobianMatrix(voltage, busType, numPV, numPQ);
+
+    // Calculate DeltaTheta DeltaV array
+    std::vector<double> dTdV = GaussianElimination(jacobMatrix, dPdQ);
+
+    // Update voltage array
+    int indexDT = 0;
+    int indexDV = numPQ + numPV;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] != BUS_SLACK) {
+            if(busType[i] == BUS_PV) {
+                double newV = std::abs(voltage[i]);
+                double newT = std::arg(voltage[i]) + dTdV[indexDT];
+                voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
+                indexDT++;
+            } else {
+                double newV = std::abs(voltage[i]) + dTdV[indexDV];
+                double newT = std::arg(voltage[i]) + dTdV[indexDT];
+                voltage[i] = std::complex<double>(newV * std::cos(newT), newV * std::sin(newT));
+                indexDV++;
+                indexDT++;
+            }
+        }
+    }
+}
+
+std::vector<std::vector<double> > PowerFlow::CalculateJacobianMatrix(std::vector<std::complex<double> > voltage,
+                                                                     std::vector<BusType> busType,
+                                                                     int numPV,
+                                                                     int numPQ)
+{
+    // Jacobian matrix
+    std::vector<std::vector<double> > jacobMatrix;
+    jacobMatrix.resize(2 * numPQ + numPV);  // Jabobian matrix size
+    for(int i = 0; i < (2 * numPQ + numPV); i++) jacobMatrix[i].resize(2 * numPQ + numPV);
+
+    // Calculate the submatrices
+    //{ H
+    std::vector<std::vector<double> > sMatrixH;
+    sMatrixH.resize(numPQ + numPV);  // submatrix H size
+    for(int i = 0; i < (numPQ + numPV); i++) sMatrixH[i].resize(numPQ + numPV);
+
+    int busNumI = 0;
+    int busNumJ = 0;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] != BUS_SLACK) {
+            for(int j = 0; j < m_numberOfBuses; j++) {
+                if(busType[j] != BUS_SLACK) {
+                    if(busNumI != busNumJ) {
+                        double vi = std::abs(voltage[i]);
+                        double ti = std::arg(voltage[i]);
+                        double vj = std::abs(voltage[j]);
+                        double tj = std::arg(voltage[j]);
+                        sMatrixH[busNumI][busNumJ] =
+                            vi * vj *
+                            (m_yBus[i][j].real() * std::sin(ti - tj) - m_yBus[i][j].imag() * std::cos(ti - tj));
+                    } else  // busNumI == busNumJ
+                    {
+                        std::complex<double> sj = std::complex<double>(0.0, 0.0);
+                        for(int k = 0; k < m_numberOfBuses; k++)
+                            sj += voltage[j] * std::conj(m_yBus[j][k]) * std::conj(voltage[k]);
+                        sMatrixH[busNumI][busNumJ] =
+                            -sj.imag() - m_yBus[j][j].imag() * std::pow(std::abs(voltage[j]), 2.0);
+                    }
+                    busNumJ++;
+                }
+            }
+            busNumI++;
+            busNumJ = 0;
+        }
+    }
+    //}
+    //{ N
+    std::vector<std::vector<double> > sMatrixN;
+    sMatrixN.resize(numPQ + numPV);  // submatrix N size
+    for(int i = 0; i < (numPQ + numPV); i++) sMatrixN[i].resize(numPQ);
+
+    busNumI = 0;
+    busNumJ = 0;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] != BUS_SLACK) {
+            for(int j = 0; j < m_numberOfBuses; j++) {
+                if(busType[j] == BUS_PQ)  // Only PQ buses
+                {
+                    if(busNumI != busNumJ) {
+                        double vi = std::abs(voltage[i]);
+                        double ti = std::arg(voltage[i]);
+                        double vj = std::abs(voltage[j]);
+                        double tj = std::arg(voltage[j]);
+                        sMatrixN[busNumI][busNumJ] =
+                            vi * vj *
+                            (m_yBus[i][j].real() * std::cos(ti - tj) + m_yBus[i][j].imag() * std::sin(ti - tj));
+                    } else  // busNumI == busNumJ
+                    {
+                        std::complex<double> sj = std::complex<double>(0.0, 0.0);
+                        for(int k = 0; k < m_numberOfBuses; k++)
+                            sj += voltage[j] * std::conj(m_yBus[j][k]) * std::conj(voltage[k]);
+                        sMatrixN[busNumI][busNumJ] =
+                            sj.real() + m_yBus[j][j].real() * std::pow(std::abs(voltage[j]), 2.0);
+                    }
+                    busNumJ++;
+                }
+            }
+            busNumI++;
+            busNumJ = 0;
+        }
+    }
+    //}
+    //{ J
+    std::vector<std::vector<double> > sMatrixJ;
+    sMatrixJ.resize(numPQ);  // submatrix J size
+    for(int i = 0; i < (numPQ); i++) sMatrixJ[i].resize(numPQ + numPV);
+
+    busNumI = 0;
+    busNumJ = 0;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] == BUS_PQ)  // Only PQ
+        {
+            for(int j = 0; j < m_numberOfBuses; j++) {
+                if(busType[j] == BUS_SLACK) {
+                    if(busNumI != busNumJ) {
+                        double vi = std::abs(voltage[i]);
+                        double ti = std::arg(voltage[i]);
+                        double vj = std::abs(voltage[j]);
+                        double tj = std::arg(voltage[j]);
+                        sMatrixJ[busNumI][busNumJ] =
+                            -vi * vj *
+                            (m_yBus[i][j].real() * std::cos(ti - tj) + m_yBus[i][j].imag() * std::sin(ti - tj));
+                    } else  // busNumI == busNumJ
+                    {
+                        std::complex<double> sj = std::complex<double>(0.0, 0.0);
+                        for(int k = 0; k < m_numberOfBuses; k++)
+                            sj += voltage[j] * std::conj(m_yBus[j][k]) * std::conj(voltage[k]);
+                        sMatrixJ[busNumI][busNumJ] =
+                            sj.real() - m_yBus[j][j].real() * std::pow(std::abs(voltage[j]), 2.0);
+                    }
+                    busNumJ++;
+                }
+            }
+            busNumI++;
+            busNumJ = 0;
+        }
+    }
+    //}
+    //{ L
+    std::vector<std::vector<double> > sMatrixL;
+    sMatrixL.resize(numPQ);  // submatrix L size
+    for(int i = 0; i < (numPQ); i++) sMatrixL[i].resize(numPQ);
+
+    busNumI = 0;
+    busNumJ = 0;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] == BUS_PQ)  // nćo deve ser referźncia e nćo deve ser PV, ou seja é PQ
+        {
+            for(int j = 0; j < m_numberOfBuses; j++) {
+                if(busType[j] == BUS_PQ)  // nćo deve ser referźncia e nćo deve ser PV, ou seja é PQ
+                {
+                    if(busNumI != busNumJ) {
+                        double vi = std::abs(voltage[i]);
+                        double ti = std::arg(voltage[i]);
+                        double vj = std::abs(voltage[j]);
+                        double tj = std::arg(voltage[j]);
+                        sMatrixL[busNumI][busNumJ] =
+                            vi * vj *
+                            (m_yBus[i][j].real() * std::sin(ti - tj) - m_yBus[i][j].imag() * std::cos(ti - tj));
+                    } else  // busNumI == busNumJ
+                    {
+                        std::complex<double> sj = std::complex<double>(0.0, 0.0);
+                        for(int k = 0; k < m_numberOfBuses; k++)
+                            sj += voltage[j] * std::conj(m_yBus[j][k]) * std::conj(voltage[k]);
+                        sMatrixL[busNumI][busNumJ] =
+                            -sj.imag() - m_yBus[j][j].imag() * std::pow(std::abs(voltage[j]), 2.0);
+                    }
+                    busNumJ++;
+                }
+            }
+            busNumI++;
+            busNumJ = 0;
+        }
+    }
+    //}
+
+    // Fill Jacobian matrix
+    for(int i = 0; i < (2 * numPQ + numPV); i++) {
+        for(int j = 0; j < (2 * numPQ + numPV); j++) {
+            if(i < (numPQ + numPV)) {
+                if(j < (numPQ + numPV))
+                    jacobMatrix[i][j] = sMatrixH[i][j];
+                else  // j >= numPQ + numPV
+                    jacobMatrix[i][j] = sMatrixN[i][j - (numPQ + numPV)];
+            } else  // i >= numPQ + numPV
+            {
+                if(j < (numPQ + numPV))
+                    jacobMatrix[i][j] = sMatrixJ[i - (numPQ + numPV)][j];
+                else  // j >= numPQ + numPV
+                    jacobMatrix[i][j] = sMatrixL[i - (numPQ + numPV)][j - (numPQ + numPV)];
+            }
+        }
+    }
+    return jacobMatrix;
+}
+
+bool PowerFlow::CheckReactiveLimits(std::vector<BusType>& busType,
+                                    std::vector<ReactiveLimits>& reactiveLimit,
+                                    std::vector<std::complex<double> > power,
+                                    std::vector<std::complex<double> > loadPower)
+{
+    bool limitReach = false;
+    for(int i = 0; i < m_numberOfBuses; i++) {
+        if(busType[i] == BUS_PV) {
+            if(reactiveLimit[i].maxLimitType == RL_LIMITED) {
+                if(power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) {
+                    power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
+                    busType[i] = BUS_PQ;
+                    reactiveLimit[i].limitReached = RL_MAX_REACHED;
+                    limitReach = true;
+                }
+            }
+            if(reactiveLimit[i].minLimitType == RL_LIMITED) {
+                if(power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) {
+                    power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
+                    busType[i] = BUS_PQ;
+                    reactiveLimit[i].limitReached = RL_MIN_REACHED;
+                    limitReach = true;
+                }
+            }
+        }
+    }
+    return limitReach;
+}