Merge branch 'merge-powerflow' into 'development'

Merge powerflow

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

libcps @ 31975777

-Subproject commit 2a1de9a9dbd92c97f86c368187e9f361ce264b3b
+Subproject commit 3197577719d0cbde8faa769f7f37be5ed5b67dde

Examples/CIM/IEEE_EU_LV_reduced/README.md

+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

Examples/Cxx/CIM/CIGRE_MV_PowerFlowTest.cpp

@@ -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;

Examples/Cxx/CIM/CIGRE_MV_PowerFlowTest_LoadProfiles.cpp

@@ -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();

Examples/Cxx/CIM/IEEE_LV_PowerFlowTest.cpp

+/**
+ * @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;
+}
+

Examples/Cxx/CIM/IEEE_LV_PowerFlowTest_LP.cpp

+/**
+ * @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;
+}
+

Examples/Cxx/CMakeLists.txt

@@ -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
+
 	)
 
 	

Examples/Cxx/Shmem/Shmem_CIGRE_MV_PowerFlowTest_LoadProfiles.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);
 

Include/dpsim/NRPSolver.h

@@ -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;
 		};
+
     };
 
 

Source/NRPSolver.cpp

@@ -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)
+	{
+		//Calculate the reactive power for the PV buses:
+		calculate_Q(npq, npv);
+	}