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_ExponentialDiode.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_ExponentialDiode.h

@@ -0,0 +1,84 @@
+/* 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 ExponentialDiode using the Shockley ideal diode equation
+class ExponentialDiode : public MNASimPowerComp<Real>,
+                         public SharedFactory<ExponentialDiode>,
+                         public MNANonlinearVariableCompInterface {
+protected:
+  Matrix Jacobian = Matrix::Zero(1, 1);
+
+public:
+  //Reverse-bias saturation current. Default: mI_S = 0.000001
+  const CPS::Attribute<Real>::Ptr mI_S;
+  //Thermal voltage. Default: mV_T = 0.027
+  const CPS::Attribute<Real>::Ptr mV_T;
+
+  /// Defines UID, name, component parameters and logging level
+  ExponentialDiode(String uid, String name,
+                   Logger::Level logLevel = Logger::Level::off);
+  /// Defines name, component parameters and logging level
+  ExponentialDiode(String name, Logger::Level logLevel = Logger::Level::off)
+      : ExponentialDiode(name, name, logLevel) {}
+
+  // #### General ####
+  ///
+  SimPowerComp<Real>::Ptr clone(String name) override;
+  /// Initializes component from power flow data
+  void initializeFromNodesAndTerminals(Real frequency) override;
+  //Sets reverse-bias saturation current I_S and thermal voltage V_T.
+  //If this method is not called, the diode will have the following
+  //values by default: I_S = 0.000001 (A) and V_T = 0.027 (V).
+  void setParameters(Real I_S, Real V_T);
+
+  // #### 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_ExponentialDiode.cpp
 
 	EMT/EMT_Ph3_CurrentSource.cpp
 	EMT/EMT_Ph3_VoltageSource.cpp

dpsim-models/src/EMT/EMT_Ph1_ExponentialDiode.cpp

@@ -0,0 +1,178 @@
+/* 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_ExponentialDiode.h>
+
+using namespace CPS;
+
+EMT::Ph1::ExponentialDiode::ExponentialDiode(String uid, String name,
+                                             Logger::Level logLevel)
+    : MNASimPowerComp<Real>(uid, name, false, false, logLevel),
+      mI_S(mAttributes->create<Real>("I_S")),
+      mV_T(mAttributes->create<Real>("V_T")) {
+  mPhaseType = PhaseType::Single;
+  setTerminalNumber(2);
+  **mI_S = 0.000001;
+  **mV_T = 0.027;
+  **mIntfVoltage = Matrix::Zero(1, 1);
+  **mIntfCurrent = Matrix::Zero(1, 1);
+}
+
+void EMT::Ph1::ExponentialDiode::setParameters(Real I_S, Real V_T) {
+  **mI_S = I_S;
+  **mV_T = V_T;
+}
+
+SimPowerComp<Real>::Ptr EMT::Ph1::ExponentialDiode::clone(String name) {
+  auto copy = ExponentialDiode::make(name, mLogLevel);
+  copy->setParameters(**mI_S, **mV_T);
+  return copy;
+}
+
+void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::mnaCompPreStep(Real time, Int timeStepCount) {
+  mnaCompApplyRightSideVectorStamp(**mRightVector);
+}
+
+void EMT::Ph1::ExponentialDiode::mnaCompPostStep(
+    Real time, Int timeStepCount, Attribute<Matrix>::Ptr &leftVector) {
+  mnaUpdateVoltage(**leftVector);
+  mnaUpdateCurrent(**leftVector);
+}
+
+void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::mnaCompUpdateCurrent(
+    const Matrix &leftVector) {
+  (**mIntfCurrent)(0, 0) =
+      (**mI_S) * (expf((**mIntfVoltage)(0, 0) / (**mV_T)) - 1.);
+}
+
+void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::updateJacobian() {
+  Jacobian(0, 0) =
+      (**mI_S / (**mV_T)) * expf((**mIntfVoltage)(0, 0) / (**mV_T));
+}

dpsim/examples/cxx/CMakeLists.txt

@@ -35,6 +35,7 @@ set(CIRCUIT_SOURCES
 	Circuits/EMT_PiLine.cpp
 	Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp
 	Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp
+	Circuits/EMT_Ph1_ExponentialDiode_test.cpp
 
 	# EMT examples with PF initialization
 	Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp
@@ -136,6 +137,7 @@ list(APPEND TEST_SOURCES
 	Circuits/DP_SMIB_ReducedOrderSGIterative_LoadStep.cpp
 	Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp
 	Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp
+	Circuits/EMT_Ph1_ExponentialDiode_test.cpp
 )
 
 if(WITH_SUNDIALS)