Iterative MNA solver

dpsim-models/include/dpsim-models/Components.h

@@ -82,6 +82,7 @@
 #endif
 #include <dpsim-models/EMT/EMT_Ph1_Capacitor.h>
 #include <dpsim-models/EMT/EMT_Ph1_CurrentSource.h>
+#include <dpsim-models/EMT/EMT_Ph1_Diode.h>
 #include <dpsim-models/EMT/EMT_Ph1_Inductor.h>
 #include <dpsim-models/EMT/EMT_Ph1_Resistor.h>
 #include <dpsim-models/EMT/EMT_Ph1_VoltageSource.h>

dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Diode.h

@@ -0,0 +1,82 @@
+/* Copyright 2017-2021 Institute for Automation of Complex Power Systems,
+ *                     EONERC, RWTH Aachen University
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at https://mozilla.org/MPL/2.0/.
+ *********************************************************************************/
+
+#pragma once
+
+#include <dpsim-models/MNASimPowerComp.h>
+#include <dpsim-models/Solver/MNAInterface.h>
+//#include <dpsim-models/Solver/MNAVariableCompInterface.h>
+#include <dpsim-models/Solver/MNANonlinearVariableCompInterface.h>
+
+namespace CPS {
+namespace EMT {
+namespace Ph1 {
+/// EMT Diode
+class Diode : public MNASimPowerComp<Real>,
+              public SharedFactory<Diode>,
+              public MNANonlinearVariableCompInterface {
+protected:
+  ///TODO: mI_S and mV_T as const CPS::Attribute<Real>::Ptr, also change in
+  // pybind EMTComponents
+  Real mI_S = 0.000001;
+  Real mV_T = 0.027;
+  Matrix Jacobian = Matrix::Zero(1, 1);
+  Real itVoltage = 0.;
+  Real itCurrent = 0.;
+
+public:
+  /// Defines UID, name, component parameters and logging level
+  Diode(String uid, String name, Logger::Level logLevel = Logger::Level::off);
+  /// Defines name, component parameters and logging level
+  Diode(String name, Logger::Level logLevel = Logger::Level::off)
+      : Diode(name, name, logLevel) {}
+
+  // #### General ####
+  ///
+  SimPowerComp<Real>::Ptr clone(String name) override;
+  /// Initializes component from power flow data
+  void initializeFromNodesAndTerminals(Real frequency) override;
+
+  void setParameters(Real, Real);
+
+  // #### MNA section ####
+  /// Initializes internal variables of the component
+  void mnaCompInitialize(Real omega, Real timeStep,
+                         Attribute<Matrix>::Ptr leftSideVector) override;
+  /// Stamps system matrix
+  void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override;
+  /// Stamps right side (source) vector
+  void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override {
+  } //No right side vector stamps
+  /// Update interface voltage from MNA system result
+  void mnaCompUpdateVoltage(const Matrix &leftVector) override;
+  /// Update interface current from MNA system result
+  void mnaCompUpdateCurrent(const Matrix &leftVector) override;
+  /// MNA pre and post step operations
+  void mnaCompPreStep(Real time,
+                      Int timeStepCount) override; //No right side vector stamps
+  void mnaCompPostStep(Real time, Int timeStepCount,
+                       Attribute<Matrix>::Ptr &leftVector) override;
+  /// add MNA pre and post step dependencies
+  void
+  mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies,
+                                 AttributeBase::List &attributeDependencies,
+                                 AttributeBase::List &modifiedAttributes,
+                                 Attribute<Matrix>::Ptr &leftVector) override;
+
+  virtual void iterationUpdate(const Matrix &leftVector) override;
+
+  virtual bool hasParameterChanged() override { return true; }
+
+  void calculateNonlinearFunctionResult() override;
+
+  void updateJacobian();
+};
+} // namespace Ph1
+} // namespace EMT
+} // namespace CPS

dpsim-models/include/dpsim-models/Solver/MNANonlinearVariableCompInterface.h

@@ -0,0 +1,36 @@
+/* Copyright 2017-2021 Institute for Automation of Complex Power Systems,
+ *                     EONERC, RWTH Aachen University
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at https://mozilla.org/MPL/2.0/.
+ *********************************************************************************/
+
+#pragma once
+
+#include <dpsim-models/Config.h>
+#include <dpsim-models/Definitions.h>
+#include <dpsim-models/Solver/MNAVariableCompInterface.h>
+
+namespace CPS {
+/// MNA interface to be used by nonlinear elements that require recomputing
+//	of the system matrix
+class MNANonlinearVariableCompInterface
+    : virtual public MNAVariableCompInterface {
+public:
+  typedef std::shared_ptr<MNANonlinearVariableCompInterface> NonlinearPtr;
+  typedef std::vector<NonlinearPtr> NonlinearList;
+
+  const Attribute<Matrix>::Ptr mNonlinearFunctionStamp;
+  std::vector<std::pair<UInt, UInt>> mNonlinearVariableSystemMatrixEntries;
+
+  /// Returns value of the component's defining voltage/current equation
+  //	based on current circuit quantities
+  virtual void calculateNonlinearFunctionResult() = 0;
+  virtual void iterationUpdate(const Matrix &leftVector) = 0;
+
+protected:
+  MNANonlinearVariableCompInterface()
+      : mNonlinearFunctionStamp(AttributeStatic<Matrix>::make()){};
+};
+} // namespace CPS

dpsim-models/src/CMakeLists.txt

@@ -70,6 +70,7 @@ list(APPEND MODELS_SOURCES
 	EMT/EMT_Ph1_VoltageSourceRamp.cpp
 	EMT/EMT_Ph1_VoltageSourceNorton.cpp
 	EMT/EMT_Ph1_Switch.cpp
+	EMT/EMT_Ph1_Diode.cpp
 
 	EMT/EMT_Ph3_CurrentSource.cpp
 	EMT/EMT_Ph3_VoltageSource.cpp

dpsim-models/src/EMT/EMT_Ph1_Diode.cpp

@@ -0,0 +1,166 @@
+/* Copyright 2017-2021 Institute for Automation of Complex Power Systems,
+ *                     EONERC, RWTH Aachen University
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at https://mozilla.org/MPL/2.0/.
+ ******************************************************************************/
+
+#include <dpsim-models/EMT/EMT_Ph1_Diode.h>
+
+using namespace CPS;
+
+EMT::Ph1::Diode::Diode(String uid, String name, Logger::Level logLevel)
+    : MNASimPowerComp<Real>(uid, name, false, false, logLevel) {
+  mPhaseType = PhaseType::Single;
+  setTerminalNumber(2);
+  **mIntfVoltage = Matrix::Zero(1, 1);
+  **mIntfCurrent = Matrix::Zero(1, 1);
+}
+
+void EMT::Ph1::Diode::setParameters(Real I_S, Real V_T) {
+  mI_S = I_S;
+  mV_T = V_T;
+}
+
+SimPowerComp<Real>::Ptr EMT::Ph1::Diode::clone(String name) {
+  auto copy = Diode::make(name, mLogLevel);
+  copy->setParameters(mI_S, mV_T);
+  return copy;
+}
+
+void EMT::Ph1::Diode::initializeFromNodesAndTerminals(Real frequency) {
+
+  // IntfVoltage initialization for each phase
+  MatrixComp vInitABC = Matrix::Zero(1, 1);
+  vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) -
+                   RMS3PH_TO_PEAK1PH * initialSingleVoltage(0);
+  (**mIntfVoltage)(0, 0) = vInitABC(0, 0).real();
+
+  ///FIXME: 	Initialization should include solving the system once to obtain
+  //			the actual values solving the
+  //			system for the 0th time step. As of now abnormal current
+  //			values for the 0th time step are to be expected.
+
+  (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.);
+
+  SPDLOG_LOGGER_INFO(mSLog,
+                     "\n--- Initialization from powerflow ---"
+                     "\nVoltage across: {:f}"
+                     "\nCurrent: {:f}"
+                     "\nTerminal 0 voltage: {:f}"
+                     "\nTerminal 1 voltage: {:f}"
+                     "\n--- Initialization from powerflow finished ---",
+                     (**mIntfVoltage)(0, 0), (**mIntfCurrent)(0, 0),
+                     initialSingleVoltage(0).real(),
+                     initialSingleVoltage(1).real());
+}
+
+void EMT::Ph1::Diode::mnaCompInitialize(Real omega, Real timeStep,
+                                        Attribute<Matrix>::Ptr leftVector) {
+
+  updateMatrixNodeIndices();
+
+  **mRightVector = Matrix::Zero(leftVector->get().rows(), 1);
+  **mNonlinearFunctionStamp = Matrix::Zero(leftVector->get().rows(), 1);
+  calculateNonlinearFunctionResult();
+  updateJacobian();
+
+  mMnaTasks.push_back(std::make_shared<MnaPostStep>(*this, leftVector));
+}
+
+void EMT::Ph1::Diode::mnaCompApplySystemMatrixStamp(
+    SparseMatrixRow &systemMatrix) {
+  // Set diagonal entries
+  if (terminalNotGrounded(0)) {
+    // set upper left block, 3x3 entries
+    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0),
+                             matrixNodeIndex(0, 0), Jacobian(0, 0));
+  }
+  if (terminalNotGrounded(1)) {
+    // set buttom right block, 3x3 entries
+    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0),
+                             matrixNodeIndex(1, 0), Jacobian(0, 0));
+  }
+  // Set off diagonal blocks, 2x3x3 entries
+  if (terminalNotGrounded(0) && terminalNotGrounded(1)) {
+    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0),
+                             matrixNodeIndex(1, 0), -Jacobian(0, 0));
+
+    Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0),
+                             matrixNodeIndex(0, 0), -Jacobian(0, 0));
+  }
+}
+
+void EMT::Ph1::Diode::mnaCompAddPostStepDependencies(
+    AttributeBase::List &prevStepDependencies,
+    AttributeBase::List &attributeDependencies,
+    AttributeBase::List &modifiedAttributes,
+    Attribute<Matrix>::Ptr &leftVector) {
+  attributeDependencies.push_back(leftVector);
+  modifiedAttributes.push_back(attribute("v_intf"));
+  modifiedAttributes.push_back(attribute("i_intf"));
+}
+
+void EMT::Ph1::Diode::mnaCompPreStep(Real time, Int timeStepCount) {
+  mnaCompApplyRightSideVectorStamp(**mRightVector);
+}
+
+void EMT::Ph1::Diode::mnaCompPostStep(Real time, Int timeStepCount,
+                                      Attribute<Matrix>::Ptr &leftVector) {
+  mnaUpdateVoltage(**leftVector);
+  mnaUpdateCurrent(**leftVector);
+}
+
+void EMT::Ph1::Diode::mnaCompUpdateVoltage(const Matrix &leftVector) {
+  // v1 - v0
+  **mIntfVoltage = Matrix::Zero(3, 1);
+  if (terminalNotGrounded(1)) {
+    (**mIntfVoltage)(0, 0) =
+        Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0));
+  }
+  if (terminalNotGrounded(0)) {
+    (**mIntfVoltage)(0, 0) =
+        (**mIntfVoltage)(0, 0) -
+        Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0));
+  }
+}
+
+void EMT::Ph1::Diode::mnaCompUpdateCurrent(const Matrix &leftVector) {
+  (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.);
+}
+
+void EMT::Ph1::Diode::iterationUpdate(const Matrix &leftVector) {
+  //Update phase voltages
+  if (terminalNotGrounded(1)) {
+    (**mIntfVoltage)(0, 0) =
+        Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0));
+  }
+  if (terminalNotGrounded(0)) {
+    (**mIntfVoltage)(0, 0) =
+        (**mIntfVoltage)(0, 0) -
+        Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0));
+  }
+
+  calculateNonlinearFunctionResult();
+
+  updateJacobian();
+}
+
+void EMT::Ph1::Diode::calculateNonlinearFunctionResult() {
+
+  (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.);
+
+  if (terminalNotGrounded(1)) {
+    Math::setVectorElement(**mNonlinearFunctionStamp, matrixNodeIndex(1, 0),
+                           (**mIntfCurrent)(0, 0));
+  }
+  if (terminalNotGrounded(0)) {
+    Math::setVectorElement(**mNonlinearFunctionStamp, matrixNodeIndex(0, 0),
+                           -(**mIntfCurrent)(0, 0));
+  }
+}
+
+void EMT::Ph1::Diode::updateJacobian() {
+  Jacobian(0, 0) = (mI_S / mV_T) * expf((**mIntfVoltage)(0, 0) / mV_T);
+}

dpsim/examples/cxx/CMakeLists.txt

@@ -33,6 +33,7 @@ set(CIRCUIT_SOURCES
 	#Circuits/EMT_ResVS_RL_Switch.cpp
 	Circuits/EMT_VSI.cpp
 	Circuits/EMT_PiLine.cpp
+	Circuits/EMT_Ph1_Diode_test.cpp
 
 	# EMT examples with PF initialization
 	Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp

dpsim/examples/cxx/Circuits/EMT_Ph1_Diode_test.cpp

@@ -0,0 +1,59 @@
+#include "../Examples.h"
+#include <DPsim.h>
+
+using namespace DPsim;
+using namespace CPS::EMT;
+
+void EMT_Ph1_Diode() {
+  // Define simulation scenario
+  Real timeStep = 0.0001;
+  Real finalTime = 0.1;
+  String simName = "EMT_Ph1_Diode_test";
+  Logger::setLogDir("logs/" + simName);
+
+  // Nodes
+  auto n1 = SimNode::make("n1", PhaseType::Single);
+  auto n2 = SimNode::make("n1", PhaseType::Single);
+
+  // Components
+
+  auto vs0 = Ph1::VoltageSource::make("vs0");
+  vs0->setParameters(CPS::Complex(1., 0.), 50.0);
+
+  auto load = CPS::EMT::Ph1::Resistor::make("Load");
+  load->setParameters(10.);
+
+  auto diode = Ph1::Diode::make("Diode");
+
+  // Topology
+  load->connect(SimNode::List{n1, n2});
+
+  diode->connect(SimNode::List{n2, SimNode::GND}); //SimNode::GND, n2
+
+  vs0->connect(SimNode::List{SimNode::GND, n1});
+
+  // Define system topology
+  auto sys = SystemTopology(50, SystemNodeList{n1, n2},
+                            SystemComponentList{load, vs0, diode});
+
+  // Logging
+  auto logger = DataLogger::make(simName);
+  logger->logAttribute("I_Diode", diode->attribute("i_intf"));
+  logger->logAttribute("V_Diode", diode->attribute("v_intf"));
+
+  Simulation sim(simName, Logger::Level::info);
+  sim.doInitFromNodesAndTerminals(true);
+  sim.setSystem(sys);
+  sim.addLogger(logger);
+  sim.setDomain(Domain::EMT);
+  sim.setSolverType(Solver::Type::ITERATIVEMNA);
+  sim.setDirectLinearSolverConfiguration(DirectLinearSolverConfiguration());
+  sim.setTimeStep(timeStep);
+  sim.setFinalTime(finalTime);
+  sim.run();
+}
+
+int main() {
+  EMT_Ph1_Diode();
+  return 0;
+}