Commit 346b2aa2 authored by Markus Mirz's avatar Markus Mirz
Browse files

Merge branch 'merge-powerflow' into 'development'

Merge powerflow

See merge request acs/public/simulation/dpsim!72

Former-commit-id: 1cb74c54
parents 7ca7b50b 0de4558e
Subproject commit 2a1de9a9dbd92c97f86c368187e9f361ce264b3b
Subproject commit 3197577719d0cbde8faa769f7f37be5ed5b67dde
Note: CIM data have been created with NEPLAN export option
Original file: https://git.rwth-aachen.de/acs/core/simulation/reference-examples/blob/development/Neplan/IEEE_EU_LV/IEEE_EU_LV_reduced.nepprj
Export options in NEPLAN: Use Name as RID
24e317fa301a73ebb622ce56facc47e9f568c3c2
\ No newline at end of file
963b1473266ab72656cd1ae0d059d8c8b8518563
\ No newline at end of file
This diff is collapsed.
5079c72a813811b36bb1f3a41fb2d83f7644957e
\ No newline at end of file
......@@ -20,7 +20,7 @@
#include "cps/CIM/Reader.h"
#include <DPsim.h>
#include "cps/loadProfileReader.h"
#include "cps/LoadProfileReader.h"
using namespace std;
using namespace DPsim;
......
......@@ -20,7 +20,7 @@
#include "cps/CIM/Reader.h"
#include <DPsim.h>
#include "cps/loadProfileReader.h"
#include "cps/LoadProfileReader.h"
using namespace std;
using namespace DPsim;
......@@ -44,6 +44,26 @@ int main(int argc, char** argv){
#elif defined(__linux__) || defined(__APPLE__)
String loadProfilePath("Examples/CSV/CIGRE_MV_NoTap/");
#endif
std::map<String,String> assignList = {
// {load mRID, file name}
{"LOAD-H-1", "Load_H_1"},
{"LOAD-H-3", "Load_H_3"},
{"LOAD-H-4", "Load_H_4"},
{"LOAD-H-5", "Load_H_5"},
{"LOAD-H-6", "Load_H_6"},
{"LOAD-H-8", "Load_H_8"},
{"LOAD-H-10", "Load_H_10"},
{"LOAD-H-11", "Load_H_11"},
{"LOAD-H-12", "Load_H_12"},
{"LOAD-H-14", "Load_H_14"},
{"LOAD-I-1", "Load_I_1"},
{"LOAD-I-3", "Load_I_3"},
{"LOAD-I-7", "Load_I_7"},
{"LOAD-I-9", "Load_I_9"},
{"LOAD-I-10", "Load_I_10"},
{"LOAD-I-12", "Load_I_12"},
{"LOAD-I-13", "Load_I_13"},
{"LOAD-I-14", "Load_I_14"}};
std::list<string> filenames = {
path + "Rootnet_FULL_NE_06J16h_DI.xml",
......@@ -56,15 +76,69 @@ int main(int argc, char** argv){
CIM::Reader reader(simName, Logger::Level::INFO, Logger::Level::NONE);
SystemTopology system = reader.loadCIM(system_freq, filenames, CPS::Domain::Static);
loadProfileAssigner(simName, loadProfilePath, system, Logger::Level::INFO);
//load profile lpreader
LoadProfileReader lpreader(simName, loadProfilePath, assignList, Logger::Level::INFO);
lpreader.assign(system, 0, 1.5, 15, LoadProfileReader::Mode::MANUAL);
auto logger = DPsim::DataLogger::make(simName);
for (auto node : system.mNodes)
{
logger->addAttribute(node->name(), node->attribute("v"));
std::list<std::shared_ptr<CPS::Static::Ph1::PiLine>> lines;
for (auto comp : system.mComponentsAtNode[node]) {
if (std::shared_ptr<CPS::Static::Ph1::PiLine> line =
std::dynamic_pointer_cast<CPS::Static::Ph1::PiLine>(comp))
{
String current=(node->name() == line->node(0)->name()) ? ("current") : ("current_1");
String p_branch = (node->name() == line->node(0)->name()) ? ("p_branch") : ("p_branch_1");
String q_branch = (node->name() == line->node(0)->name()) ? ("q_branch") : ("q_branch_1");
logger->addAttribute(line->name() + "." + node->name() + ".I", line->attribute<Complex>(current));
logger->addAttribute(line->name() + "." + node->name() + ".P", line->attribute<Real>(p_branch));
logger->addAttribute(line->name() + "." + node->name() + ".Q", line->attribute<Real>(q_branch));
lines.push_back(line);
}
else if (std::shared_ptr<CPS::Static::Ph1::externalGridInjection> extnet =
std::dynamic_pointer_cast<CPS::Static::Ph1::externalGridInjection>(comp))
{
logger->addAttribute(node->name() + ".Pinj", extnet->attribute<Real>("p_inj"));
logger->addAttribute(node->name() + ".Qinj", extnet->attribute<Real>("q_inj"));
}
else if (std::shared_ptr<CPS::Static::Ph1::Transformer> trafo =
std::dynamic_pointer_cast<CPS::Static::Ph1::Transformer>(comp))
{
String current = (node->name() == trafo->node(0)->name()) ? ("current") : ("current_1");
String p_branch = (node->name() == trafo->node(0)->name()) ? ("p_branch") : ("p_branch_1");
String q_branch = (node->name() == trafo->node(0)->name()) ? ("q_branch") : ("q_branch_1");
logger->addAttribute(trafo->name() + "." + node->name() + ".I", trafo->attribute<Complex>(current));
logger->addAttribute(trafo->name() + "." + node->name() + ".P", trafo->attribute<Real>(p_branch));
logger->addAttribute(trafo->name() + "." + node->name() + ".Q", trafo->attribute<Real>(q_branch));
}
}
// get nodal injection from specific line or transformer
// (the first line obj connected to the node or, if none, the first trafo)
if (!lines.empty()) {
logger->addAttribute(node->name() + ".Pinj", lines.front()->attribute<Real>("p_inj"));
logger->addAttribute(node->name() + ".Qinj", lines.front()->attribute<Real>("q_inj"));
}
else
{
for (auto comp : system.mComponentsAtNode[node]) {
if (std::shared_ptr<CPS::Static::Ph1::Transformer> trafo =
std::dynamic_pointer_cast<CPS::Static::Ph1::Transformer>(comp))
{
logger->addAttribute(node->name() + ".Pinj", trafo->attribute<Real>("p_inj"));
logger->addAttribute(node->name() + ".Qinj", trafo->attribute<Real>("q_inj"));
break;
}
}
}
}
Simulation sim(simName, system, 1, 120, Domain::Static, Solver::Type::NRP, Logger::Level::NONE, true);
Simulation sim(simName, system, 1.5, 15, Domain::Static, Solver::Type::NRP, Logger::Level::INFO, true);
sim.addLogger(logger);
sim.run();
......
/**
* @author Jan Dinkelbach <jdinkelbach@eonerc.rwth-aachen.de>
* @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
*
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*********************************************************************************/
#include "cps/CIM/Reader.h"
#include <DPsim.h>
#include "cps/LoadProfileReader.h"
using namespace std;
using namespace DPsim;
using namespace CPS;
using namespace CPS::CIM;
/*
* This example runs the powerflow for the CIGRE MV benchmark system (neglecting the tap changers of the transformers)
*/
int main(int argc, char** argv){
#ifdef _WIN32
String path("..\\..\\..\\..\\dpsim\\Examples\\CIM\\IEEE_EU_LV_reduced\\");
#elif defined(__linux__) || defined(__APPLE__)
String path("Examples/CIM/IEEE_EU_LV_reduced/");
#endif
std::list<string> filenames = {
path + "Rootnet_FULL_NE_13J16h_DI.xml",
path + "Rootnet_FULL_NE_13J16h_EQ.xml",
path + "Rootnet_FULL_NE_13J16h_SV.xml",
path + "Rootnet_FULL_NE_13J16h_TP.xml"
};
String simName = "CIGRE-LV-reduced";
CPS::Real system_freq = 50;
CIM::Reader reader(simName, Logger::Level::INFO, Logger::Level::NONE);
SystemTopology system = reader.loadCIM(system_freq, filenames, CPS::Domain::Static);
auto logger = DPsim::DataLogger::make(simName);
for (auto node : system.mNodes)
{
logger->addAttribute(node->name(), node->attribute("v"));
}
Simulation sim(simName, system, 1, 5, Domain::Static, Solver::Type::NRP, Logger::Level::INFO, true);
sim.addLogger(logger);
sim.run();
return 0;
}
/**
* @author Jan Dinkelbach <jdinkelbach@eonerc.rwth-aachen.de>
* @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
*
* DPsim
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*********************************************************************************/
#include "cps/CIM/Reader.h"
#include <DPsim.h>
#include "cps/LoadProfileReader.h"
using namespace std;
using namespace DPsim;
using namespace CPS;
using namespace CPS::CIM;
/*
* This example runs the powerflow for the IEEE LV EU benchmark system
*/
int main(int argc, char** argv){
#ifdef _WIN32
String path("..\\..\\..\\..\\dpsim\\Examples\\CIM\\IEEE_EU_LV_reduced\\");
#elif defined(__linux__) || defined(__APPLE__)
String path("Examples/CIM/IEEE_EU_LV_reduced/");
#endif
#ifdef _WIN32
String loadProfilePath("..\\..\\..\\..\\SimulationData\\Load\\IEEE_European_LV_TestFeeder_v2_Profiles\\csv\\");
#elif defined(__linux__) || defined(__APPLE__)
String loadProfilePath("../SimulationData/Load/IEEE_European_LV_TestFeeder_v2_Profiles/csv/");
#endif
std::list<string> filenames = {
path + "Rootnet_FULL_NE_13J16h_DI.xml",
path + "Rootnet_FULL_NE_13J16h_EQ.xml",
path + "Rootnet_FULL_NE_13J16h_SV.xml",
path + "Rootnet_FULL_NE_13J16h_TP.xml"
};
String simName = "IEEE-LV-reduced";
CPS::Real system_freq = 50;
/*
create assign list
*/
std::map<String, String> assignList;
for (const auto & entry : std::experimental::filesystem::directory_iterator(loadProfilePath))
{
string filename = entry.path().filename().string();
string load_number;
for (auto &c : filename) {
if (isdigit(c)) {
load_number += c;
}
}
assignList.insert(std::pair<string,string>("PQ"+load_number,"Load_profile_"+load_number));
}
CIM::Reader reader(simName, Logger::Level::INFO, Logger::Level::NONE);
SystemTopology system = reader.loadCIM(system_freq, filenames, CPS::Domain::Static);
//load profile lpreader
LoadProfileReader lpreader(simName, loadProfilePath, assignList, Logger::Level::INFO);
lpreader.assign(system, 1, 1, 30, LoadProfileReader::Mode::MANUAL, LoadProfileReader::DataFormat::HHMMSS);
auto logger = DPsim::DataLogger::make(simName);
for (auto node : system.mNodes)
{
logger->addAttribute(node->name(), node->attribute("v"));
}
Simulation sim(simName, system, 1, 30, Domain::Static, Solver::Type::NRP, Logger::Level::INFO, true);
sim.addLogger(logger);
sim.run();
return 0;
}
......@@ -97,6 +97,10 @@ if(WITH_CIM)
# powerflow example
CIM/CIGRE_MV_PowerFlowTest.cpp
CIM/CIGRE_MV_PowerFlowTest_LoadProfiles.cpp
CIM/IEEE_LV_PowerFlowTest.cpp
CIM/IEEE_LV_PowerFlowTest_LP.cpp
)
......
......@@ -59,7 +59,8 @@ int main(int argc, char** argv){
CIM::Reader reader(simName, Logger::Level::DEBUG, Logger::Level::NONE);
SystemTopology sys = reader.loadCIM(system_freq, filenames, CPS::Domain::Static);
loadProfileAssigner(simName, loadProfilePath, sys, Logger::Level::INFO);
loadProfileReader lpreader(simName, loadProfilePath, Logger::Level::INFO);
lpreader.assign(sys, 1, 1, 60, loadProfileReader::Mode::AUTO);
RealTimeSimulation sim(simName, sys, args.timeStep, args.duration, args.solver.domain, args.solver.type, args.logLevel);
......
......@@ -87,7 +87,6 @@ namespace DPsim {
///Imaginary part of Y
CPS::Real B(int i, int j);
// solver settings
double tolerance=10e-9;
/// Maximum number of iterations
......@@ -128,9 +127,9 @@ namespace DPsim {
*
* \sa Solver_State
*/
virtual void NRP_initialize();
virtual void NRP_initialize(CPS::Real time);
virtual void powerFlow();
virtual Bool powerFlow(Bool with_iwamoto = false);
void determinePowerFlowBusType();
......@@ -158,17 +157,17 @@ namespace DPsim {
std::vector<int> LastPV;
//******* vectors of node indices for pf calculation ********
// #### vectors of node indices for pf calculation ####
std::vector<CPS::UInt> PQBusIndices;
std::vector<CPS::UInt> PVBusIndices;
std::vector<CPS::UInt> slackBusIndex;
//******* vectors of nodes for pf calculation ********
// #### vectors of nodes for pf calculation ####
CPS::TopologicalNode::List PQBuses;
CPS::TopologicalNode::List PVBuses;
CPS::TopologicalNode::List slackBus;
/****** vectors of components ******/
// #### vectors of components ####
std::vector<std::shared_ptr<CPS::Static::Ph1::Transformer>> Transformers;
std::vector<std::shared_ptr<CPS::Static::Ph1::SynchronGenerator>> SynchronGenerators;
std::vector<std::shared_ptr<CPS::Static::Ph1::Load>> Loads;
......@@ -194,7 +193,7 @@ namespace DPsim {
void update_solution(CPS::Vector X, CPS::UInt npq, CPS::UInt npv);
void get_increments(CPS::Vector X, CPS::Vector &incV, CPS::Vector &incD, CPS::UInt npq, CPS::UInt npv);
void get_increments(const CPS::Vector& X, CPS::Vector &incV, CPS::Vector &incD, CPS::UInt npq, CPS::UInt npv);
void calculate_slack_power(); //calculate the slack bus power
......@@ -205,9 +204,9 @@ namespace DPsim {
void compile();
void generate_initial_solution(bool keep_last_solution=false);
void generate_initial_solution(CPS::Real time, bool keep_last_solution = false);
void set_solution();
void set_solution(CPS::Bool didConverge);
void calculate_flows();
/* TODO: Z matrix composition to be fixed after moving from Circuit to Solver_NRpolar
......@@ -222,6 +221,10 @@ namespace DPsim {
void correct_PVbuses_violating_Q(CPS::UInt &npq, CPS::UInt &npv, CPS::Matrix &J, CPS::Vector &K, CPS::Vector &X);
void calculate_branch_flow();
void calculate_nodal_injection();
// TODO proper tasking system integration for parallelism
class SolveTask : public CPS::Task {
public:
......@@ -235,6 +238,7 @@ namespace DPsim {
private:
NRpolarSolver& mSolver;
};
};
......
......@@ -97,14 +97,16 @@ NRpolarSolver::NRpolarSolver(CPS::String simName, CPS::SystemTopology & sysTopol
determinePowerFlowBusType();
}
void NRpolarSolver::NRP_initialize(){
void NRpolarSolver::NRP_initialize(Real time){
generate_initial_solution();
generate_initial_solution(time);
mLog.info() << "#### NEWTON-RAPHSON POLAR SOLVER " << std::endl;
/* print addmittance matrix for debug only
mLog.info() << "#### Admittance Matrix: " <<std::endl
<< Eigen::Matrix<CPS::Complex,Eigen::Dynamic,Eigen::Dynamic>(Y) << std::endl;
*/
mLog.info() << "#### Create index vectors for power flow solver:" << std::endl;
BUSES.reserve(
......@@ -347,7 +349,7 @@ void NRpolarSolver::determinePowerFlowBusType() {
// this could be integrated into the function that determines node type (PV,PQ)
void NRpolarSolver::generate_initial_solution(bool keep_last_solution) {
void NRpolarSolver::generate_initial_solution(Real time, bool keep_last_solution) {
resize_sol(SysTopology.mNodes.size());
resize_complex_sol(SysTopology.mNodes.size());
......@@ -361,7 +363,7 @@ void NRpolarSolver::generate_initial_solution(bool keep_last_solution) {
if (std::shared_ptr<CPS::Static::Ph1::Load> load = std::dynamic_pointer_cast<CPS::Static::Ph1::Load>(comp)) {
if (load->use_profile) {
load->updatePQ();
load->updatePQ(time);
}
sol_P(pq->simNode()) -= load->mPQ->attribute<CPS::Real>("P")->get() / Sbase;
sol_Q(pq->simNode()) -= load->mPQ->attribute<CPS::Real>("Q")->get() / Sbase;
......@@ -411,11 +413,13 @@ void NRpolarSolver::generate_initial_solution(bool keep_last_solution) {
sol_initialized = true;
sol_complex_initialized = true;
/*
mLog.info() << "#### Initial solution: " << std::endl;
mLog.info() << "P\t\tQ\t\tV\t\tD" << std::endl;
for (UInt i = 0; i < sol_length; i++) {
mLog.info() << sol_P[i] << "\t" << sol_Q[i] << "\t" << sol_V[i] << "\t" << sol_D[i] << std::endl;
}
*/
}
......@@ -454,7 +458,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
for (UInt a = 0; a < npqpv; a++) { //rows
k = PQPV[a];
//diagonal
J(a, a) = -Q(k) - B(k, k) * V.coeff(k) * V.coeff(k);
J.coeffRef(a, a) = -Q(k) - B(k, k) * V.coeff(k) * V.coeff(k);
//non diagonal elements
for (UInt b = 0; b < npqpv; b++) {
......@@ -464,7 +468,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
*(G(k, j) * sin(D.coeff(k) - D.coeff(j))
- B(k, j) * cos(D.coeff(k) - D.coeff(j)));
//if (val != 0.0)
J(a, b) = val;
J.coeffRef(a, b) = val;
}
}
}
......@@ -477,7 +481,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
//diagonal
//std::cout << "J2D:" << (a + da) << "," << (a + db) << std::endl;
if (a < npq)
J(a + da, a + db) = P(k) + G(k, k) * V.coeff(k) * V.coeff(k);
J.coeffRef(a + da, a + db) = P(k) + G(k, k) * V.coeff(k) * V.coeff(k);
//non diagonal elements
for (UInt b = 0; b < npq; b++) {
......@@ -488,7 +492,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
+ B(k, j) * sin(sol_D.coeff(k) - sol_D.coeff(j)));
//if (val != 0.0)
//std::cout << "J2ij:" << (a + da) << "," << (b + db) << std::endl;
J(a + da, b + db) = val;
J.coeffRef(a + da, b + db) = val;
}
}
}
......@@ -501,7 +505,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
k = PQPV[a];
//diagonal
//std::cout << "J3:" << (a + da) << "," << (a + db) << std::endl;
J(a + da, a + db) = P(k) - G(k, k) * V.coeff(k) * V.coeff(k);
J.coeffRef(a + da, a + db) = P(k) - G(k, k) * V.coeff(k) * V.coeff(k);
//non diagonal elements
for (UInt b = 0; b < npqpv; b++) {
......@@ -512,7 +516,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
+ B(k, j) * sin(D.coeff(k) - D.coeff(j)));
//if (val != 0.0)
//std::cout << "J3:" << (a + da) << "," << (b + db) << std::endl;
J(a + da, b + db) = -val;
J.coeffRef(a + da, b + db) = -val;
}
}
}
......@@ -524,7 +528,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
k = PQPV[a];
//diagonal
//std::cout << "J4:" << (a + da) << "," << (a + db) << std::endl;
J(a + da, a + db) = Q(k) - B(k, k) * V.coeff(k) * V.coeff(k);
J.coeffRef(a + da, a + db) = Q(k) - B(k, k) * V.coeff(k) * V.coeff(k);
//non diagonal elements
for (UInt b = 0; b < npq; b++) {
......@@ -535,7 +539,7 @@ void NRpolarSolver::Jacobian(Matrix &J, Vector &V, Vector &D, UInt npq, UInt npv
- B(k, j) * cos(D.coeff(k) - D.coeff(j)));
if (val != 0.0) {
//std::cout << "J4:" << (a + da) << "," << (b + db) << std::endl;
J(a + da, b + db) = val;
J.coeffRef(a + da, b + db) = val;
}
}
}
......@@ -627,8 +631,10 @@ double NRpolarSolver::mu(Matrix &J, Matrix &J2, Vector &F, Vector &dV, Vector &d
Vector b = J * (dx);
Vector c(2*npq+npv); //= dx. * b * 0.5;
for (UInt i=0;i<(2*npq+npv); i++) //this loop is because EIGEN does not want to perform this simple element wise vector multiplication...
c(i) = dx.coeff(i) * b.coeff(i) * 0.5;
//for (UInt i=0;i<(2*npq+npv); i++) //this loop is because EIGEN does not want to perform this simple element wise vector multiplication...
// c(i) = dx.coeff(i) * b.coeff(i) * 0.5;
Vector temp_c = J2 * (dx);
c = dx.cwiseProduct(temp_c)*0.5;
double g0 = -1* a.dot(b);
double g1 = b.dot(b) + 2 * a.dot(c);
......@@ -672,7 +678,7 @@ bool NRpolarSolver::converged(const Vector& PQinc, UInt npqpvpq) const
}
void NRpolarSolver::get_increments(Vector X, Vector &incV, Vector &incD, UInt npq, UInt npv){
void NRpolarSolver::get_increments(const Vector& X, Vector &incV, Vector &incD, UInt npq, UInt npv){
UInt npqpv = npq + npv;
UInt k;
......@@ -711,8 +717,7 @@ void NRpolarSolver::update_solution(Vector X, UInt npq, UInt npv) {
}
}
void NRpolarSolver::powerFlow()
Bool NRpolarSolver::powerFlow(Bool with_iwamoto)
{
UInt npq = PQBusIndices.size();
......@@ -730,20 +735,26 @@ void NRpolarSolver::powerFlow()
// First shot: Perhaps the model already converged?
get_power_inc(K, npq, npv);
auto didConverge = converged(K, npqpvpq);
Bool didConverge = converged(K, npqpvpq);
Iterations = 0;
for (unsigned i = 0; i < maxIterations && ! didConverge; ++i) {
for (unsigned i = 1; i < maxIterations && ! didConverge; ++i) {
Jacobian(J, sol_V, sol_D, npq, npv);
Eigen::FullPivLU<Matrix>lu(J); //Full pivot LU
auto sparseJ = J.sparseView();
Eigen::SparseLU<SparseMatrix>lu(sparseJ); // sparse matrix LU decomposition
X = lu.solve(K);
get_increments(X, incV, incD, npq, npv);
auto mu_ = mu(J, J2, K, incV, incD, X, npq, npv);
//upgrade the solution
update_solution(X * mu_, npq, npv);
//update the solution
if (with_iwamoto)
{
auto mu_ = mu(J, J2, K, incV, incD, X, npq, npv);
update_solution(X * mu_, npq, npv);
}
else
{
update_solution(X, npq, npv);
}
//Calculate the increment of power for the new iteration
get_power_inc(K, npq, npv);
......@@ -751,23 +762,22 @@ void NRpolarSolver::powerFlow()
didConverge = converged(K, npqpvpq);
Iterations = i;
}
if (didConverge)
{