From 73a48d92000c3eabb5ebc22e1bf87e64db5caf6d Mon Sep 17 00:00:00 2001 From: Marvin Tollnitsch Date: Mon, 30 Sep 2024 16:25:57 +0200 Subject: [PATCH 1/4] Add EMT three-phase SSN component models for L, C and RLC + add "EMT_Ph3_SSN_Inductor.h" and "EMT_Ph3_SSN_Inductor.cpp". For the inductor, the resulting SSN model is identical to the resistive companion model. + add "EMT_Ph3_SSN_Capacitor.h" and "EMT_Ph3_SSN_Capacitor.cpp". For the capacitor, the SSN model results in an additional mesh equation similar to a voltage source, expanding the dimensions of the system matrix by 1x1. The resistive companion model should be used instead. + add "EMT_Ph3_SSN_Full_Serial_RLC.h" and "EMT_Ph3_SSN_Full_Serial_RLC.cpp". The implemented RLC circuit consists of R, L and C , completely wired in series. Instead of four required nodes for those three components, the combined SSN model only requires two, reducing the dimensions of the system matrix. All three implemented SSN components can be used together with RC component models and with the already implemented MNA Solver classes. * modify "Components.h" to include "EMT_Ph3_SSN_Capacitor.h", "EMT_Ph3_SSN_Inductor.h" and "EMT_Ph3_SSN_Full_Serial_RLC.h" * modify "dpsim-models/src/CMakeLists.txt" to include "EMT_Ph3_SSN_Capacitor.cpp", "EMT_Ph3_SSN_Inductor.cpp" and "EMT_Ph3_SSN_Full_Serial_RLC.cpp" Signed-off-by: Marvin Tollnitsch --- .../include/dpsim-models/Components.h | 3 + .../dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h | 86 +++++ .../EMT/EMT_Ph3_SSN_Full_Serial_RLC.h | 97 +++++ .../dpsim-models/EMT/EMT_Ph3_SSN_Inductor.h | 86 +++++ dpsim-models/src/CMakeLists.txt | 3 + .../src/EMT/EMT_Ph3_SSN_Capacitor.cpp | 238 ++++++++++++ .../src/EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp | 362 ++++++++++++++++++ dpsim-models/src/EMT/EMT_Ph3_SSN_Inductor.cpp | 193 ++++++++++ 8 files changed, 1068 insertions(+) create mode 100644 dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h create mode 100644 dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Full_Serial_RLC.h create mode 100644 dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Inductor.h create mode 100644 dpsim-models/src/EMT/EMT_Ph3_SSN_Capacitor.cpp create mode 100644 dpsim-models/src/EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp create mode 100644 dpsim-models/src/EMT/EMT_Ph3_SSN_Inductor.cpp diff --git a/dpsim-models/include/dpsim-models/Components.h b/dpsim-models/include/dpsim-models/Components.h index d92170fe28..29a4f0bd58 100644 --- a/dpsim-models/include/dpsim-models/Components.h +++ b/dpsim-models/include/dpsim-models/Components.h @@ -108,6 +108,9 @@ #include #include #include +#include +#include +#include #ifdef WITH_SUNDIALS #include #endif diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h new file mode 100644 index 0000000000..dbe14f4085 --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Capacitor.h @@ -0,0 +1,86 @@ +/* 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 +#include +#include + +namespace CPS { +namespace EMT { +namespace Ph3 { +namespace SSN { +/// \brief Capacitor model +/// +/// The SSN-capacitor is represented by a DC equivalent circuit which +/// corresponds to one iteration of the trapezoidal integration method. +/// The equivalent DC circuit is a resistance in series with a voltage +/// source (real voltage source). The resistance is constant for a defined +/// time step and system frequency and the voltage source changes for each +/// iteration. This component is meant to show an I-type SSN model +/// implementation based on a simple example element. The RC model should be +/// used over this if focussing on SSN model implementation is not the +/// circuit's goal since this model increases the system matrix dimensions +/// by 1x1 (added virtual node, real voltage source MNA scheme). +class Capacitor final : public MNASimPowerComp, + public Base::Ph3::Capacitor, + public SharedFactory { +public: + /// Defines UID, name and logging level + Capacitor(String uid, String name, + Logger::Level logLevel = Logger::Level::off); + /// Defines name and logging level + Capacitor(String name, Logger::Level logLevel = Logger::Level::off) + : Capacitor(name, name, logLevel) {} + + SimPowerComp::Ptr clone(String name) override; + + // #### General #### + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + + // #### MNA section #### + /// Initializes internal variables of the component + void mnaCompInitialize(Real omega, Real timeStep, + Attribute::Ptr leftVector) override; + /// Stamps system matrix + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; + /// Stamps right side (source) vector + void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override; + /// 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 step operations + void mnaCompPreStep(Real time, Int timeStepCount) override; + /// MNA post step operations + void mnaCompPostStep(Real time, Int timeStepCount, + Attribute::Ptr &leftVector) override; + /// Add MNA pre step dependencies + void mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) override; + /// Add MNA post step dependencies + void + mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) override; + +private: + Matrix mDufourBKHat = Matrix::Zero(3, 3); + Matrix mDufourWKN = Matrix::Zero(3, 3); + //rightsideVector history term + Matrix mHistoricVoltage = Matrix::Zero(3, 1); +}; +} // namespace SSN +} // namespace Ph3 +} // namespace EMT +} // namespace CPS diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Full_Serial_RLC.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Full_Serial_RLC.h new file mode 100644 index 0000000000..145109cbe0 --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Full_Serial_RLC.h @@ -0,0 +1,97 @@ +/* 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 +#include +#include +#include +#include + +namespace CPS { +namespace EMT { +namespace Ph3 { +namespace SSN { +/// \brief Full_Serial_RLC +/// +/// This element represents an one port circuit consisting of a resistor, +/// an inductor and a capacitor connected in series. The terminals are at +/// the beginning and the end of the component chain. +/// The states are the capacitor voltage and the inductor current, the output +/// is the latter of those states (inductor current). The input is the voltage +/// across the whole circuit. States and past inputs are updated after each +/// time step and are used to calculate the current (input) voltage, +/// represented as MNA node voltages. +class Full_Serial_RLC final : public MNASimPowerComp, + public SharedFactory, + public Base::Ph3::Resistor, + public Base::Ph3::Inductor, + public Base::Ph3::Capacitor { +public: + /// Defines UID, name, component parameters and logging level + Full_Serial_RLC(String uid, String name, + Logger::Level logLevel = Logger::Level::off); + /// Defines name and logging level + Full_Serial_RLC(String name, Logger::Level logLevel = Logger::Level::off) + : Full_Serial_RLC(name, name, logLevel) {} + + SimPowerComp::Ptr clone(String name) override; + void setParameters(Matrix resistance, Matrix inductance, Matrix capacitance); + + // #### General #### + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + + // #### MNA section #### + /// Initializes internal variables of the component + void mnaCompInitialize(Real omega, Real timeStep, + Attribute::Ptr leftVector) override; + /// Stamps system matrix + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; + /// Stamps right side (source) vector + void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override; + /// 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 step operations + void mnaCompPreStep(Real time, Int timeStepCount) override; + /// MNA post step operations + void mnaCompPostStep(Real time, Int timeStepCount, + Attribute::Ptr &leftVector) override; + /// Add MNA pre step dependencies + void mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) override; + /// Add MNA post step dependencies + void + mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) override; + +private: + Matrix mState = Matrix::Zero(6, 1); + Matrix mYHistory = Matrix::Zero(3, 1); + + Matrix mDufourUNT = Matrix::Zero(3, 1); + + Matrix mDufourAKHat = Matrix::Zero(6, 6); + Matrix mDufourBKHat = Matrix::Zero(6, 3); + Matrix mDufourBKNHat = Matrix::Zero(6, 3); + Matrix mDufourWKN = Matrix::Zero(3, 3); + Matrix mDufourCKN = Matrix::Zero(3, 6); + + void ssnUpdateState(); +}; +} // namespace SSN +} // namespace Ph3 +} // namespace EMT +} // namespace CPS \ No newline at end of file diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Inductor.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Inductor.h new file mode 100644 index 0000000000..a1e1871496 --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SSN_Inductor.h @@ -0,0 +1,86 @@ +/* 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 +#include +#include + +namespace CPS { +namespace EMT { +namespace Ph3 { +namespace SSN { +/// \brief Inductor +/// +/// The SSN-inductor is represented by a DC equivalent circuit which +/// corresponds to one iteration of the trapezoidal integration method. +/// The equivalent DC circuit is a resistance in parallel with a current +/// source. The resistance is constant for a defined time step and system +/// frequency and the current source changes for each iteration. +/// This means that for the inductor the SSN and Resistive Companion models +/// are conceptionally the same and only variable names, structures and +/// state update timings differ. This component is meant to show a V-type SSN +/// model implementation based on a simple example element. Due to additional +/// calculation steps, the RC model should be used over this otherwise. +class Inductor final : public MNASimPowerComp, + public Base::Ph3::Inductor, + public SharedFactory { +public: + /// Defines UID, name, component parameters and logging level + Inductor(String uid, String name, + Logger::Level logLevel = Logger::Level::off); + /// Defines name and logging level + Inductor(String name, Logger::Level logLevel = Logger::Level::off) + : Inductor(name, name, logLevel) {} + + SimPowerComp::Ptr clone(String name) override; + + // #### General #### + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + + // #### MNA section #### + /// Initializes internal variables of the component + void mnaCompInitialize(Real omega, Real timeStep, + Attribute::Ptr leftVector) override; + /// Stamps system matrix + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; + /// Stamps right side (source) vector + void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override; + /// 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 step operations + void mnaCompPreStep(Real time, Int timeStepCount) override; + /// MNA post step operations + void mnaCompPostStep(Real time, Int timeStepCount, + Attribute::Ptr &leftVector) override; + /// Add MNA pre step dependencies + void mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) override; + /// Add MNA post step dependencies + void + mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) override; + +private: + //rightsideVector history term + Matrix mHistoricCurrent = Matrix::Zero(3, 1); + //dependency on latest Voltage, represented by Conductance in system matrix + Matrix mDufourBKHat = Matrix::Zero(3, 3); + Matrix mDufourWKN = Matrix::Zero(3, 3); +}; +} // namespace SSN +} // namespace Ph3 +} // namespace EMT +} // namespace CPS \ No newline at end of file diff --git a/dpsim-models/src/CMakeLists.txt b/dpsim-models/src/CMakeLists.txt index 795c44cd97..cb46ea72cd 100644 --- a/dpsim-models/src/CMakeLists.txt +++ b/dpsim-models/src/CMakeLists.txt @@ -102,6 +102,9 @@ list(APPEND MODELS_SOURCES EMT/EMT_Ph3_SynchronGeneratorIdeal.cpp # EMT/EMT_Ph3_SynchronGeneratorVBRStandalone.cpp EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp + EMT/EMT_Ph3_SSN_Capacitor.cpp + EMT/EMT_Ph3_SSN_Inductor.cpp + EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp SP/SP_Ph1_VoltageSource.cpp SP/SP_Ph1_Capacitor.cpp diff --git a/dpsim-models/src/EMT/EMT_Ph3_SSN_Capacitor.cpp b/dpsim-models/src/EMT/EMT_Ph3_SSN_Capacitor.cpp new file mode 100644 index 0000000000..948d9d9467 --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph3_SSN_Capacitor.cpp @@ -0,0 +1,238 @@ +/* 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 + +using namespace CPS; + +EMT::Ph3::SSN::Capacitor::Capacitor(String uid, String name, + Logger::Level logLevel) + : MNASimPowerComp(uid, name, true, true, logLevel), + Base::Ph3::Capacitor(mAttributes) { + mPhaseType = PhaseType::ABC; + setVirtualNodeNumber(1); + setTerminalNumber(2); + mHistoricVoltage = Matrix::Zero(3, 1); + **mIntfVoltage = Matrix::Zero(3, 1); + **mIntfCurrent = Matrix::Zero(3, 1); +} + +SimPowerComp::Ptr EMT::Ph3::SSN::Capacitor::clone(String name) { + auto copy = Capacitor::make(name, mLogLevel); + copy->setParameters(**mCapacitance); + return copy; +} +void EMT::Ph3::SSN::Capacitor::initializeFromNodesAndTerminals(Real frequency) { + + Real omega = 2 * PI * frequency; + MatrixComp admittance = MatrixComp::Zero(3, 3); + admittance << Complex(0, omega * (**mCapacitance)(0, 0)), + Complex(0, omega * (**mCapacitance)(0, 1)), + Complex(0, omega * (**mCapacitance)(0, 2)), + Complex(0, omega * (**mCapacitance)(1, 0)), + Complex(0, omega * (**mCapacitance)(1, 1)), + Complex(0, omega * (**mCapacitance)(1, 2)), + Complex(0, omega * (**mCapacitance)(2, 0)), + Complex(0, omega * (**mCapacitance)(2, 1)), + Complex(0, omega * (**mCapacitance)(2, 2)); + + MatrixComp vInitABC = Matrix::Zero(3, 1); + vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) - + RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); + vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B; + vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C; + **mIntfVoltage = vInitABC.real(); + **mIntfCurrent = (admittance * vInitABC).real(); + + SPDLOG_LOGGER_INFO(mSLog, + "\nCapacitance [F]: {:s}" + "\nAdmittance [S]: {:s}", + Logger::matrixToString(**mCapacitance), + Logger::matrixCompToString(admittance)); + SPDLOG_LOGGER_INFO( + mSLog, + "\n--- Initialization from powerflow ---" + "\nVoltage across: {:s}" + "\nCurrent: {:s}" + "\nTerminal 0 voltage: {:s}" + "\nTerminal 1 voltage: {:s}" + "\n--- Initialization from powerflow finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(1))); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompInitialize( + Real omega, Real timeStep, Attribute::Ptr leftVector) { + updateMatrixNodeIndices(); + + mHistoricVoltage = mDufourBKHat * **mIntfCurrent + **mIntfVoltage; + mDufourBKHat = timeStep * (2.0 * **mCapacitance).inverse(); + mDufourWKN = mDufourBKHat; + + **mRightVector = Matrix::Zero(leftVector->get().rows(), 1); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompApplySystemMatrixStamp( + SparseMatrixRow &systemMatrix) { + if (terminalNotGrounded(0)) { + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + -1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + matrixNodeIndex(0, 0), -1); + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + -1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + matrixNodeIndex(0, 1), -1); + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + -1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + matrixNodeIndex(0, 2), -1); + } + if (terminalNotGrounded(1)) { + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + 1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + matrixNodeIndex(1, 0), 1); + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + 1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + matrixNodeIndex(1, 1), 1); + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + 1); + Math::addToMatrixElement(systemMatrix, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + matrixNodeIndex(1, 2), 1); + } + //mesh equations are independent from grounded terminals + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(0, 0)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(0, 1)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(0, 2)); + + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(1, 0)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(1, 1)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(1, 2)); + + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), -mDufourWKN(2, 0)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), -mDufourWKN(2, 1)); + Math::addToMatrixElement( + systemMatrix, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), -mDufourWKN(2, 2)); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompApplyRightSideVectorStamp( + Matrix &rightVector) { + mHistoricVoltage = mDufourBKHat * **mIntfCurrent + **mIntfVoltage; + Math::setVectorElement(rightVector, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::A), + mHistoricVoltage(0, 0)); + Math::setVectorElement(rightVector, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::B), + mHistoricVoltage(1, 0)); + Math::setVectorElement(rightVector, + mVirtualNodes[0]->matrixNodeIndex(PhaseType::C), + mHistoricVoltage(2, 0)); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) { + // actually depends on C, but then we'd have to modify the system matrix anyway + prevStepDependencies.push_back(mIntfCurrent); + prevStepDependencies.push_back(mIntfVoltage); + modifiedAttributes.push_back(mRightVector); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompPreStep(Real time, Int timeStepCount) { + mnaCompApplyRightSideVectorStamp(**mRightVector); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompAddPostStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) { + attributeDependencies.push_back(leftVector); + modifiedAttributes.push_back(mIntfVoltage); + modifiedAttributes.push_back(mIntfCurrent); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompPostStep( + Real time, Int timeStepCount, Attribute::Ptr &leftVector) { + mnaCompUpdateVoltage(**leftVector); + mnaCompUpdateCurrent(**leftVector); +} + +void EMT::Ph3::SSN::Capacitor::mnaCompUpdateVoltage(const Matrix &leftVector) { + // v1 - v0 + **mIntfVoltage = Matrix::Zero(3, 1); + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + (**mIntfVoltage)(1, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 1)); + (**mIntfVoltage)(2, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 2)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = + (**mIntfVoltage)(0, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + (**mIntfVoltage)(1, 0) = + (**mIntfVoltage)(1, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1)); + (**mIntfVoltage)(2, 0) = + (**mIntfVoltage)(2, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); + } +} + +void EMT::Ph3::SSN::Capacitor::mnaCompUpdateCurrent(const Matrix &leftVector) { + (**mIntfCurrent)(0, 0) = Math::realFromVectorElement( + leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::A)); + (**mIntfCurrent)(1, 0) = Math::realFromVectorElement( + leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::B)); + (**mIntfCurrent)(2, 0) = Math::realFromVectorElement( + leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C)); + + SPDLOG_LOGGER_DEBUG(mSLog, "\nCurrent: {:s}", + Logger::matrixToString(**mIntfCurrent)); +} diff --git a/dpsim-models/src/EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp b/dpsim-models/src/EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp new file mode 100644 index 0000000000..1b53bb2a77 --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph3_SSN_Full_Serial_RLC.cpp @@ -0,0 +1,362 @@ +/* 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 + +using namespace CPS; + +EMT::Ph3::SSN::Full_Serial_RLC::Full_Serial_RLC(String uid, String name, + Logger::Level logLevel) + : MNASimPowerComp(uid, name, true, true, logLevel), + Base::Ph3::Resistor(mAttributes), Base::Ph3::Inductor(mAttributes), + Base::Ph3::Capacitor(mAttributes) { + + mPhaseType = PhaseType::ABC; + + **mIntfVoltage = Matrix::Zero(3, 1); + **mIntfCurrent = Matrix::Zero(3, 1); + setTerminalNumber(2); +} + +SimPowerComp::Ptr EMT::Ph3::SSN::Full_Serial_RLC::clone(String name) { + auto copy = Full_Serial_RLC::make(name, mLogLevel); + copy->setParameters(**mResistance, **mInductance, **mCapacitance); + return copy; +} + +void EMT::Ph3::SSN::Full_Serial_RLC::initializeFromNodesAndTerminals( + Real frequency) { + + Real omega = 2 * PI * frequency; + + MatrixComp impedance = MatrixComp::Zero(3, 3); + + impedance << Complex((**mResistance)(0, 0), + omega * (**mInductance)(0, 0) - + 1. / (omega * (**mCapacitance)(0, 0))), + Complex((**mResistance)(0, 1), omega * (**mInductance)(0, 1) - + 1. / (omega * (**mCapacitance)(0, 1))), + Complex((**mResistance)(0, 2), omega * (**mInductance)(0, 2) - + 1. / (omega * (**mCapacitance)(0, 2))), + Complex((**mResistance)(1, 0), omega * (**mInductance)(1, 0) - + 1. / (omega * (**mCapacitance)(1, 0))), + Complex((**mResistance)(1, 1), omega * (**mInductance)(1, 1) - + 1. / (omega * (**mCapacitance)(1, 1))), + Complex((**mResistance)(1, 2), omega * (**mInductance)(1, 2) - + 1. / (omega * (**mCapacitance)(1, 2))), + Complex((**mResistance)(2, 0), omega * (**mInductance)(2, 0) - + 1. / (omega * (**mCapacitance)(2, 0))), + Complex((**mResistance)(2, 1), omega * (**mInductance)(2, 1) - + 1. / (omega * (**mCapacitance)(2, 1))), + Complex((**mResistance)(2, 2), omega * (**mInductance)(2, 2) - + 1. / (omega * (**mCapacitance)(2, 2))); + + MatrixComp vInitABC = Matrix::Zero(3, 1); + vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) - + RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); + vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B; + vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C; + **mIntfVoltage = vInitABC.real(); + MatrixComp admittance = impedance.inverse(); + **mIntfCurrent = (admittance * vInitABC).real(); + + ///FIXME: mIntfCurrent gets initialized to a vector of NaNs! + **mIntfCurrent = Matrix::Zero(3, 3); + + SPDLOG_LOGGER_INFO(mSLog, "\nImpedance [Ohm]: {:s}", + Logger::matrixCompToString(impedance)); + SPDLOG_LOGGER_INFO( + mSLog, + "\n--- Initialization from powerflow ---" + "\nVoltage across: {:s}" + "\nCurrent: {:s}" + "\nTerminal 0 voltage: {:s}" + "\nTerminal 1 voltage: {:s}" + "\n--- Initialization from powerflow finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(1))); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompInitialize( + Real omega, Real timeStep, Attribute::Ptr leftVector) { + + updateMatrixNodeIndices(); + + mState = Matrix::Zero(6, 1); + mYHistory = Matrix::Zero(3, 1); + + //Fill mDufourAKHat + //State Equation one, phases A,B,C: top left submatrix + mDufourAKHat(0, 0) = + 1. - ((2. * (timeStep * timeStep)) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mCapacitance)(0, 0) * (**mResistance)(0, 0) + + timeStep * timeStep)); + mDufourAKHat(1, 1) = + 1. - ((2. * (timeStep * timeStep)) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mCapacitance)(1, 1) * (**mResistance)(1, 1) + + timeStep * timeStep)); + mDufourAKHat(2, 2) = + 1. - ((2. * (timeStep * timeStep)) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mCapacitance)(2, 2) * (**mResistance)(2, 2) + + timeStep * timeStep)); + //State Equation one, phases A,B,C: top right submatrix + mDufourAKHat(0, 3) = + (timeStep / (2. * (**mCapacitance)(0, 0))) * + (1. + ((4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) - + 2. * timeStep * (**mResistance)(0, 0) * (**mCapacitance)(0, 0) - + (timeStep * timeStep)) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mResistance)(0, 0) * (**mCapacitance)(0, 0) + + (timeStep * timeStep)))); + mDufourAKHat(1, 4) = + (timeStep / (2. * (**mCapacitance)(1, 1))) * + (1. + ((4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) - + 2. * timeStep * (**mResistance)(1, 1) * (**mCapacitance)(1, 1) - + (timeStep * timeStep)) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mResistance)(1, 1) * (**mCapacitance)(1, 1) + + (timeStep * timeStep)))); + mDufourAKHat(2, 5) = + (timeStep / (2. * (**mCapacitance)(2, 2))) * + (1. + ((4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) - + 2. * timeStep * (**mResistance)(2, 2) * (**mCapacitance)(2, 2) - + (timeStep * timeStep)) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mResistance)(2, 2) * (**mCapacitance)(2, 2) + + (timeStep * timeStep)))); + //State Equation two, phases A,B,C: bottom left submatrix + mDufourAKHat(3, 0) = + -((4. * (**mCapacitance)(0, 0) * timeStep) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mCapacitance)(0, 0) * (**mResistance)(0, 0) + + (timeStep * timeStep))); + mDufourAKHat(4, 1) = + -((4. * (**mCapacitance)(1, 1) * timeStep) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mCapacitance)(1, 1) * (**mResistance)(1, 1) + + (timeStep * timeStep))); + mDufourAKHat(5, 2) = + -((4. * (**mCapacitance)(2, 2) * timeStep) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mCapacitance)(2, 2) * (**mResistance)(2, 2) + + (timeStep * timeStep))); + //State Equation two, phases A,B,C: bottom right submatrix + mDufourAKHat(3, 3) = + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) - + 2. * timeStep * (**mResistance)(0, 0) * (**mCapacitance)(0, 0) - + (timeStep * timeStep)) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mResistance)(0, 0) * (**mCapacitance)(0, 0) + + (timeStep * timeStep)); + mDufourAKHat(4, 4) = + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) - + 2. * timeStep * (**mResistance)(1, 1) * (**mCapacitance)(1, 1) - + (timeStep * timeStep)) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mResistance)(1, 1) * (**mCapacitance)(1, 1) + + (timeStep * timeStep)); + mDufourAKHat(5, 5) = + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) - + 2. * timeStep * (**mResistance)(2, 2) * (**mCapacitance)(2, 2) - + (timeStep * timeStep)) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mResistance)(2, 2) * (**mCapacitance)(2, 2) + + (timeStep * timeStep)); + + ///Fill mDufourBKHat + //State Equation one, phases A,B,C: top submatrix + mDufourBKHat(0, 0) = + (timeStep * timeStep) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mCapacitance)(0, 0) * (**mResistance)(0, 0) + + (timeStep * timeStep)); + mDufourBKHat(1, 1) = + (timeStep * timeStep) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mCapacitance)(1, 1) * (**mResistance)(1, 1) + + (timeStep * timeStep)); + mDufourBKHat(2, 2) = + (timeStep * timeStep) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mCapacitance)(2, 2) * (**mResistance)(2, 2) + + (timeStep * timeStep)); + + //State Equation two, phases A,B,C: bottom submatrix + mDufourBKHat(3, 0) = + (timeStep * 2. * (**mCapacitance)(0, 0)) / + (4. * (**mInductance)(0, 0) * (**mCapacitance)(0, 0) + + 2. * timeStep * (**mCapacitance)(0, 0) * (**mResistance)(0, 0) + + (timeStep * timeStep)); + mDufourBKHat(4, 1) = + (timeStep * 2. * (**mCapacitance)(1, 1)) / + (4. * (**mInductance)(1, 1) * (**mCapacitance)(1, 1) + + 2. * timeStep * (**mCapacitance)(1, 1) * (**mResistance)(1, 1) + + (timeStep * timeStep)); + mDufourBKHat(5, 2) = + (timeStep * 2. * (**mCapacitance)(2, 2)) / + (4. * (**mInductance)(2, 2) * (**mCapacitance)(2, 2) + + 2. * timeStep * (**mCapacitance)(2, 2) * (**mResistance)(2, 2) + + (timeStep * timeStep)); + + mDufourBKNHat = mDufourBKHat; + + mDufourCKN(0, 3) = 1.; + mDufourCKN(1, 4) = 1.; + mDufourCKN(2, 5) = 1.; + + mDufourWKN = mDufourCKN * mDufourBKHat; + + ///FIXME: mIntfCurrent is state 2 and is potentially directly initialized by other initialization methodes (e.g. FromNodesAndTerminals). + /// State 1, which is Voltage over the capacitor, is not directly initialized and has to be calculated from the states. This is why + /// the current must be reset as it would be altered as well. However, the old value of state one (time step "-1") is unknown or "zero" + /// in this case, so calculation of state 1 would always assume zero as the past value of state 1 and also takes mIntfCurrent + /// for the calculation ->State 1 is ahead of state 2 by one step, but also always (wrongly?) assumes past state 1 to be zero. + /// How to handle properly? + mState(3, 0) = (**mIntfCurrent)(0, 0); + mState(4, 0) = (**mIntfCurrent)(1, 0); + mState(5, 0) = (**mIntfCurrent)(2, 0); + ssnUpdateState(); + mState(3, 0) = (**mIntfCurrent)(0, 0); + mState(4, 0) = (**mIntfCurrent)(1, 0); + mState(5, 0) = (**mIntfCurrent)(2, 0); + + **mRightVector = Matrix::Zero(leftVector->get().rows(), 1); + + SPDLOG_LOGGER_INFO(mSLog, + "\n--- MNA initialization ---" + "\nInitial voltage {:s}" + "\nInitial current {:s}" + "\n--- MNA initialization finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompApplySystemMatrixStamp( + SparseMatrixRow &systemMatrix) { + MNAStampUtils::stampConductanceMatrix( + mDufourWKN, systemMatrix, matrixNodeIndex(0), matrixNodeIndex(1), + terminalNotGrounded(0), terminalNotGrounded(1), mSLog); + + SPDLOG_LOGGER_INFO(mSLog, "\nConductance matrix: {:s}", + Logger::matrixToString(mDufourWKN)); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompApplyRightSideVectorStamp( + Matrix &rightVector) { + // Update history term + mYHistory = + mDufourCKN * (mDufourAKHat * mState + mDufourBKHat * **mIntfVoltage); + + if (terminalNotGrounded(0)) { + Math::setVectorElement(rightVector, matrixNodeIndex(0, 0), mYHistory(0, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(0, 1), mYHistory(1, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(0, 2), mYHistory(2, 0)); + } + if (terminalNotGrounded(1)) { + Math::setVectorElement(rightVector, matrixNodeIndex(1, 0), + -mYHistory(0, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(1, 1), + -mYHistory(1, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(1, 2), + -mYHistory(2, 0)); + } + SPDLOG_LOGGER_DEBUG( + mSLog, "\nHistory current term (mnaCompApplyRightSideVectorStamp): {:s}", + Logger::matrixToString(mYHistory)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) { + // actually depends on L,C, but then we'd have to modify the system matrix anyway + prevStepDependencies.push_back(mIntfVoltage); + prevStepDependencies.push_back(mIntfCurrent); + modifiedAttributes.push_back(mRightVector); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompPreStep(Real time, + Int timeStepCount) { + mnaCompApplyRightSideVectorStamp(**mRightVector); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompAddPostStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) { + attributeDependencies.push_back(leftVector); + modifiedAttributes.push_back(mIntfVoltage); + modifiedAttributes.push_back(mIntfCurrent); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompPostStep( + Real time, Int timeStepCount, Attribute::Ptr &leftVector) { + mnaCompUpdateVoltage(**leftVector); + mnaCompUpdateCurrent(**leftVector); + ssnUpdateState(); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompUpdateVoltage( + const Matrix &leftVector) { + // v1 - v0 + mDufourUNT = **mIntfVoltage; + **mIntfVoltage = Matrix::Zero(3, 1); + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + (**mIntfVoltage)(1, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 1)); + (**mIntfVoltage)(2, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 2)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = + (**mIntfVoltage)(0, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + (**mIntfVoltage)(1, 0) = + (**mIntfVoltage)(1, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1)); + (**mIntfVoltage)(2, 0) = + (**mIntfVoltage)(2, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); + } + SPDLOG_LOGGER_DEBUG(mSLog, "\nUpdate Voltage: {:s}", + Logger::matrixToString(**mIntfVoltage)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::mnaCompUpdateCurrent( + const Matrix &leftVector) { + **mIntfCurrent = mYHistory + mDufourWKN * **mIntfVoltage; + + SPDLOG_LOGGER_DEBUG(mSLog, "\nUpdate Current: {:s}", + Logger::matrixToString(**mIntfCurrent)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Full_Serial_RLC::setParameters(Matrix resistance, + Matrix inductance, + Matrix capacitance) { + **mInductance = inductance; + **mCapacitance = capacitance; + **mResistance = resistance; + mParametersSet = true; +} + +void EMT::Ph3::SSN::Full_Serial_RLC::ssnUpdateState() { + mState = mDufourAKHat * mState + mDufourBKHat * mDufourUNT + + mDufourBKNHat * **mIntfVoltage; +} \ No newline at end of file diff --git a/dpsim-models/src/EMT/EMT_Ph3_SSN_Inductor.cpp b/dpsim-models/src/EMT/EMT_Ph3_SSN_Inductor.cpp new file mode 100644 index 0000000000..8d4d36ac8a --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph3_SSN_Inductor.cpp @@ -0,0 +1,193 @@ +/* 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 + +using namespace CPS; + +EMT::Ph3::SSN::Inductor::Inductor(String uid, String name, + Logger::Level logLevel) + : MNASimPowerComp(uid, name, true, true, logLevel), + Base::Ph3::Inductor(mAttributes) { + mPhaseType = PhaseType::ABC; + setTerminalNumber(2); + **mIntfVoltage = Matrix::Zero(3, 1); + **mIntfCurrent = Matrix::Zero(3, 1); + mHistoricCurrent = Matrix::Zero(3, 1); +} + +SimPowerComp::Ptr EMT::Ph3::SSN::Inductor::clone(String name) { + auto copy = Inductor::make(name, mLogLevel); + copy->setParameters(**mInductance); + return copy; +} + +void EMT::Ph3::SSN::Inductor::initializeFromNodesAndTerminals(Real frequency) { + + Real omega = 2 * PI * frequency; + MatrixComp impedance = MatrixComp::Zero(3, 3); + impedance << Complex(0, omega * (**mInductance)(0, 0)), + Complex(0, omega * (**mInductance)(0, 1)), + Complex(0, omega * (**mInductance)(0, 2)), + Complex(0, omega * (**mInductance)(1, 0)), + Complex(0, omega * (**mInductance)(1, 1)), + Complex(0, omega * (**mInductance)(1, 2)), + Complex(0, omega * (**mInductance)(2, 0)), + Complex(0, omega * (**mInductance)(2, 1)), + Complex(0, omega * (**mInductance)(2, 2)); + + MatrixComp vInitABC = Matrix::Zero(3, 1); + vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) - + RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); + vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B; + vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C; + **mIntfVoltage = vInitABC.real(); + MatrixComp admittance = impedance.inverse(); + **mIntfCurrent = (admittance * vInitABC).real(); + + SPDLOG_LOGGER_INFO(mSLog, + "\nInductance [H]: {:s}" + "\nImpedance [Ohm]: {:s}", + Logger::matrixToString(**mInductance), + Logger::matrixCompToString(impedance)); + SPDLOG_LOGGER_INFO( + mSLog, + "\n--- Initialization from powerflow ---" + "\nVoltage across: {:s}" + "\nCurrent: {:s}" + "\nTerminal 0 voltage: {:s}" + "\nTerminal 1 voltage: {:s}" + "\n--- Initialization from powerflow finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)), + Logger::phasorToString(RMS3PH_TO_PEAK1PH * initialSingleVoltage(1))); +} + +void EMT::Ph3::SSN::Inductor::mnaCompInitialize( + Real omega, Real timeStep, Attribute::Ptr leftVector) { + + updateMatrixNodeIndices(); + //update history term + mDufourBKHat = (timeStep / 2. * (**mInductance).inverse()); + mDufourWKN = mDufourBKHat; + mHistoricCurrent = mDufourBKHat * **mIntfVoltage + **mIntfCurrent; + + **mRightVector = Matrix::Zero(leftVector->get().rows(), 1); + + SPDLOG_LOGGER_INFO(mSLog, + "\n--- MNA initialization ---" + "\nInitial voltage {:s}" + "\nInitial current {:s}" + "\nhistoric current {:s}" + "\n--- MNA initialization finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), + Logger::matrixToString(mHistoricCurrent)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Inductor::mnaCompApplySystemMatrixStamp( + SparseMatrixRow &systemMatrix) { + + MNAStampUtils::stampConductanceMatrix( + mDufourWKN, systemMatrix, matrixNodeIndex(0), matrixNodeIndex(1), + terminalNotGrounded(0), terminalNotGrounded(1), mSLog); + + SPDLOG_LOGGER_INFO(mSLog, "\nConductance matrix: {:s}", + Logger::matrixToString(mDufourWKN)); +} + +void EMT::Ph3::SSN::Inductor::mnaCompApplyRightSideVectorStamp( + Matrix &rightVector) { + // Update internal state + mHistoricCurrent = mDufourBKHat * **mIntfVoltage + **mIntfCurrent; + if (terminalNotGrounded(0)) { + Math::setVectorElement(rightVector, matrixNodeIndex(0, 0), + mHistoricCurrent(0, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(0, 1), + mHistoricCurrent(1, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(0, 2), + mHistoricCurrent(2, 0)); + } + if (terminalNotGrounded(1)) { + Math::setVectorElement(rightVector, matrixNodeIndex(1, 0), + -mHistoricCurrent(0, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(1, 1), + -mHistoricCurrent(1, 0)); + Math::setVectorElement(rightVector, matrixNodeIndex(1, 2), + -mHistoricCurrent(2, 0)); + } + SPDLOG_LOGGER_DEBUG( + mSLog, "\nHistory current term (mnaCompApplyRightSideVectorStamp): {:s}", + Logger::matrixToString(mHistoricCurrent)); + mSLog->flush(); +} + +void EMT::Ph3::SSN::Inductor::mnaCompAddPreStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes) { + // actually depends on L, but then we'd have to modify the system matrix anyway + prevStepDependencies.push_back(mIntfVoltage); + prevStepDependencies.push_back(mIntfCurrent); + modifiedAttributes.push_back(mRightVector); +} + +void EMT::Ph3::SSN::Inductor::mnaCompPreStep(Real time, Int timeStepCount) { + mnaCompApplyRightSideVectorStamp(**mRightVector); +} + +void EMT::Ph3::SSN::Inductor::mnaCompAddPostStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) { + attributeDependencies.push_back(leftVector); + modifiedAttributes.push_back(mIntfVoltage); + modifiedAttributes.push_back(mIntfCurrent); +} + +void EMT::Ph3::SSN::Inductor::mnaCompPostStep( + Real time, Int timeStepCount, Attribute::Ptr &leftVector) { + mnaCompUpdateVoltage(**leftVector); + mnaCompUpdateCurrent(**leftVector); +} + +void EMT::Ph3::SSN::Inductor::mnaCompUpdateVoltage(const Matrix &leftVector) { + // v1 - v0 + **mIntfVoltage = Matrix::Zero(3, 1); + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + (**mIntfVoltage)(1, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 1)); + (**mIntfVoltage)(2, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 2)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = + (**mIntfVoltage)(0, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + (**mIntfVoltage)(1, 0) = + (**mIntfVoltage)(1, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1)); + (**mIntfVoltage)(2, 0) = + (**mIntfVoltage)(2, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); + } + SPDLOG_LOGGER_DEBUG(mSLog, "\nUpdate Voltage: {:s}", + Logger::matrixToString(**mIntfVoltage)); +} + +void EMT::Ph3::SSN::Inductor::mnaCompUpdateCurrent(const Matrix &leftVector) { + **mIntfCurrent = mHistoricCurrent + mDufourWKN * **mIntfVoltage; + SPDLOG_LOGGER_DEBUG(mSLog, "\nUpdate Current: {:s}", + Logger::matrixToString(**mIntfCurrent)); + mSLog->flush(); +} From a8615ff6ca2eb58628de50bbdad877626cb5fcf3 Mon Sep 17 00:00:00 2001 From: Marvin Tollnitsch Date: Mon, 30 Sep 2024 16:27:46 +0200 Subject: [PATCH 2/4] Implementation and declaration adjustments for EMT Ph3 CurrentSource + add implementation of "setParameters(MatrixComp currentRef, Real srcFreq = 50.0)" method in "EMT_Ph3_CurrentSource.h". Arguments are reference current (complex static phasor, amplitude and initial angle) and frequency * rename "voltageRef" parameters to "currentRef" in all "setParameters()" method declarations in "EMT_Ph3_CurrentSource.h" * update "clone()" method to use the newly implemented "setParameters()" method in "EMT_Ph3_CurrentSource.cpp" Signed-off-by: Marvin Tollnitsch --- .../dpsim-models/EMT/EMT_Ph3_CurrentSource.h | 8 ++--- .../src/EMT/EMT_Ph3_CurrentSource.cpp | 33 ++++++++++++++----- 2 files changed, 29 insertions(+), 12 deletions(-) diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_CurrentSource.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_CurrentSource.h index 58a439c987..e4c8274d70 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_CurrentSource.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_CurrentSource.h @@ -47,13 +47,13 @@ class CurrentSource : public MNASimPowerComp, // #### General #### /// Initializes component from power flow data void initializeFromNodesAndTerminals(Real frequency) override; - /// Setter for reference voltage - void setParameters(MatrixComp voltageRef, Real srcFreq = 50.0); + /// Setter for reference current + void setParameters(MatrixComp currentRef, Real srcFreq = 50.0); /// Setter for reference signal of type frequency ramp - void setParameters(MatrixComp voltageRef, Real freqStart, Real rocof, + void setParameters(MatrixComp currentRef, Real freqStart, Real rocof, Real timeStart, Real duration, bool smoothRamp = true); /// Setter for reference signal of type cosine frequency modulation - void setParameters(MatrixComp voltageRef, Real modulationFrequency, + void setParameters(MatrixComp currentRef, Real modulationFrequency, Real modulationAmplitude, Real baseFrequency = 50.0, bool zigzag = false); diff --git a/dpsim-models/src/EMT/EMT_Ph3_CurrentSource.cpp b/dpsim-models/src/EMT/EMT_Ph3_CurrentSource.cpp index cc16714b4f..a2aab0319a 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_CurrentSource.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_CurrentSource.cpp @@ -22,6 +22,30 @@ EMT::Ph3::CurrentSource::CurrentSource(String uid, String name, **mIntfVoltage = Matrix::Zero(3, 1); **mIntfCurrent = Matrix::Zero(3, 1); } +SimPowerComp::Ptr EMT::Ph3::CurrentSource::clone(String name) { + auto copy = CurrentSource::make(name, mLogLevel); + copy->setParameters(attributeTyped("I_ref")->get(), + attributeTyped("f_src")->get()); + return copy; +} + +void EMT::Ph3::CurrentSource::setParameters(MatrixComp currentRef, + Real srcFreq) { + auto srcSigSine = Signal::SineWaveGenerator::make(**mName + "_sw"); + // Complex(1,0) is used as initialPhasor for signal generator as only phase is used + srcSigSine->setParameters(Complex(1, 0), srcFreq); + mSrcSig = srcSigSine; + + **mCurrentRef = currentRef; + mSrcFreq->setReference(mSrcSig->mFreq); + + mSLog->info("\nCurrent reference phasor [I]: {:s}" + "\nFrequency [Hz]: {:s}", + Logger::matrixCompToString(currentRef), + Logger::realToString(srcFreq)); + + mParametersSet = true; +} void EMT::Ph3::CurrentSource::initializeFromNodesAndTerminals(Real frequency) { SPDLOG_LOGGER_INFO( @@ -70,13 +94,6 @@ void EMT::Ph3::CurrentSource::initializeFromNodesAndTerminals(Real frequency) { mSLog->flush(); } -SimPowerComp::Ptr EMT::Ph3::CurrentSource::clone(String name) { - auto copy = CurrentSource::make(name, mLogLevel); - // TODO: implement setParameters - // copy->setParameters(attributeTyped("I_ref")->get(), attributeTyped("f_src")->get()); - return copy; -} - void EMT::Ph3::CurrentSource::mnaCompInitialize( Real omega, Real timeStep, Attribute::Ptr leftVector) { updateMatrixNodeIndices(); @@ -167,4 +184,4 @@ void EMT::Ph3::CurrentSource::mnaCompUpdateVoltage(const Matrix &leftVector) { (**mIntfVoltage)(2, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); } -} +} \ No newline at end of file From c0881ed91a282ef23b9348e325b8b8a12c02a758 Mon Sep 17 00:00:00 2001 From: Marvin Tollnitsch Date: Mon, 30 Sep 2024 16:30:22 +0200 Subject: [PATCH 3/4] Make EMT::Ph3 CurrentSource and SSN Full_Serial_RLC usable in pybind + add EMT::Ph3::CurrentSource and EMT::Ph3::SSN::Full_Serial_RLC as classes in "EMTComponents.cpp" * apply clang-format to "EMTComponents.cpp" Signed-off-by: Marvin Tollnitsch --- dpsim/src/pybind/EMTComponents.cpp | 59 ++++++++++++++++++++++++------ 1 file changed, 47 insertions(+), 12 deletions(-) diff --git a/dpsim/src/pybind/EMTComponents.cpp b/dpsim/src/pybind/EMTComponents.cpp index 92afccbe00..edb5b06fae 100644 --- a/dpsim/src/pybind/EMTComponents.cpp +++ b/dpsim/src/pybind/EMTComponents.cpp @@ -102,13 +102,18 @@ void addEMTPh1Components(py::module_ mEMTPh1) { .def("connect", &CPS::EMT::Ph1::Inductor::connect) .def_property("L", createAttributeGetter("L"), createAttributeSetter("L")); - - py::class_, CPS::SimPowerComp, CPS::Base::Ph1::Switch>(mEMTPh1, "Switch", py::multiple_inheritance()) - .def(py::init(), "name"_a, "loglevel"_a = CPS::Logger::Level::off) - .def("set_parameters", &CPS::EMT::Ph1::Switch::setParameters, "open_resistance"_a, "closed_resistance"_a, "closed"_a = false) // cppcheck-suppress assignBoolToPointer + + py::class_, + CPS::SimPowerComp, CPS::Base::Ph1::Switch>( + mEMTPh1, "Switch", py::multiple_inheritance()) + .def(py::init(), "name"_a, + "loglevel"_a = CPS::Logger::Level::off) + .def("set_parameters", &CPS::EMT::Ph1::Switch::setParameters, + "open_resistance"_a, "closed_resistance"_a, + "closed"_a = false) // cppcheck-suppress assignBoolToPointer .def("open", &CPS::EMT::Ph1::Switch::open) .def("close", &CPS::EMT::Ph1::Switch::close) - .def("connect", &CPS::EMT::Ph1::Switch::connect); + .def("connect", &CPS::EMT::Ph1::Switch::connect); } void addEMTPh3Components(py::module_ mEMTPh3) { @@ -128,6 +133,21 @@ void addEMTPh3Components(py::module_ mEMTPh3) { .def_property("f_src", createAttributeGetter("f_src"), createAttributeSetter("f_src")); + py::class_, + CPS::SimPowerComp>(mEMTPh3, "CurrentSource", + py::multiple_inheritance()) + .def(py::init()) + .def(py::init()) + .def("set_parameters", + py::overload_cast( + &CPS::EMT::Ph3::CurrentSource::setParameters), + "I_ref"_a, "f_src"_a = 50) + .def("connect", &CPS::EMT::Ph3::VoltageSource::connect) + .def_property("I_ref", createAttributeGetter("I_ref"), + createAttributeSetter("V_ref")) + .def_property("f_src", createAttributeGetter("f_src"), + createAttributeSetter("f_src")); py::class_, CPS::SimPowerComp>(mEMTPh3, "Resistor", @@ -136,7 +156,6 @@ void addEMTPh3Components(py::module_ mEMTPh3) { .def(py::init()) .def("set_parameters", &CPS::EMT::Ph3::Resistor::setParameters, "R"_a) .def("connect", &CPS::EMT::Ph3::Resistor::connect); - py::class_, @@ -194,7 +213,7 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def("set_parameters", &CPS::EMT::Ph3::Switch::setParameters, "open_resistance"_a, "closed_resistance"_a, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "closed"_a = false) .def("open", &CPS::EMT::Ph3::Switch::openSwitch) .def("close", &CPS::EMT::Ph3::Switch::closeSwitch) @@ -369,7 +388,7 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def(py::init(), "uid"_a, "name"_a, "loglevel"_a = CPS::Logger::Level::off, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "with_trafo"_a = false) .def("set_parameters", &CPS::EMT::Ph3::AvVoltageSourceInverterDQ::setParameters, @@ -402,9 +421,8 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def(py::init(), "uid"_a, "name"_a, "loglevel"_a = CPS::Logger::Level::off, - // cppcheck-suppress assignBoolToPointer - "with_resistive_losses"_a = - false) + // cppcheck-suppress assignBoolToPointer + "with_resistive_losses"_a = false) .def("set_parameters", &CPS::EMT::Ph3::Transformer::setParameters, "nom_voltage_end_1"_a, "nom_voltage_end_2"_a, "rated_power"_a, "ratio_abs"_a, "ratio_phase"_a, "resistance"_a, "inductance"_a) @@ -428,9 +446,26 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def("set_parameters", &CPS::EMT::Ph3::SeriesSwitch::setParameters, "open_resistance"_a, "closed_resistance"_a, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "closed"_a = false) .def("open", &CPS::EMT::Ph3::SeriesSwitch::open) .def("close", &CPS::EMT::Ph3::SeriesSwitch::close) .def("connect", &CPS::EMT::Ph3::SeriesSwitch::connect); + + py::class_, + CPS::SimPowerComp>(mEMTPh3, "Full_Serial_RLC", + py::multiple_inheritance()) + .def(py::init()) + .def(py::init()) + .def("set_parameters", + &CPS::EMT::Ph3::SSN::Full_Serial_RLC::setParameters, "R"_a, "L"_a, + "C"_a) + .def("connect", &CPS::EMT::Ph3::SSN::Full_Serial_RLC::connect) + .def_property("R", createAttributeGetter("R"), + createAttributeSetter("R")) + .def_property("L", createAttributeGetter("L"), + createAttributeSetter("L")) + .def_property("C", createAttributeGetter("C"), + createAttributeSetter("C")); } From dbbc8d08f5282c936bb939eb115d04a8ce3823b3 Mon Sep 17 00:00:00 2001 From: Marvin Tollnitsch Date: Mon, 30 Sep 2024 16:31:43 +0200 Subject: [PATCH 4/4] Add test circuits for new Ph3 SSN models of L, C and RLC + add "EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp" which contains a linear circuit implemented using resistive companion component models in one simulation and the new SSN component models in another simulation for comparison + add "EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp" which contains a linear RLC circuit implemented using resistive compantion component models in one simulation and the new SSN component model in another simulation for comparison + add "EMT_Ph3_compare_RC_SSN.ipynb". This notebook simulates the circuits in the files mentioned above to show, compare and assert the results (SSN against RC) * adjust "dpsim/examples/cxx/CMakeLists.txt" to include "EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp" and "EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp" Signed-off-by: Marvin Tollnitsch --- dpsim/examples/cxx/CMakeLists.txt | 4 + .../Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp | 143 ++++++++++ .../Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp | 114 ++++++++ .../Circuits/EMT_Ph3_compare_RC_SSN.ipynb | 245 ++++++++++++++++++ 4 files changed, 506 insertions(+) create mode 100644 dpsim/examples/cxx/Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp create mode 100644 dpsim/examples/cxx/Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp create mode 100644 examples/Notebooks/Circuits/EMT_Ph3_compare_RC_SSN.ipynb diff --git a/dpsim/examples/cxx/CMakeLists.txt b/dpsim/examples/cxx/CMakeLists.txt index dba9ec1117..0edc927961 100644 --- a/dpsim/examples/cxx/CMakeLists.txt +++ b/dpsim/examples/cxx/CMakeLists.txt @@ -33,6 +33,8 @@ set(CIRCUIT_SOURCES #Circuits/EMT_ResVS_RL_Switch.cpp Circuits/EMT_VSI.cpp Circuits/EMT_PiLine.cpp + Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp + Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp # EMT examples with PF initialization Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp @@ -132,6 +134,8 @@ list(APPEND TEST_SOURCES Circuits/EMT_DP_SP_Slack_PiLine_PQLoad_FrequencyRamp_CosineFM.cpp Circuits/DP_SMIB_ReducedOrderSG_LoadStep.cpp Circuits/DP_SMIB_ReducedOrderSGIterative_LoadStep.cpp + Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp + Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp ) if(WITH_SUNDIALS) diff --git a/dpsim/examples/cxx/Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp b/dpsim/examples/cxx/Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp new file mode 100644 index 0000000000..465e9149ba --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp @@ -0,0 +1,143 @@ +#include + +using namespace DPsim; +using namespace CPS::EMT; + +void EMT_PH3_R3_C1_L1_CS1() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simNameRC = "EMT_Ph3_R3C1L1CS1_RC_vs_SSN_RC"; + + // Nodes + auto n1 = SimNode::make("n1", PhaseType::ABC); + auto n2 = SimNode::make("n2", PhaseType::ABC); + + // Components + + Matrix param = Matrix::Zero(3, 3); + param << 1., 0, 0, 0, 1., 0, 0, 0, 1.; + + auto cs0 = Ph3::CurrentSource::make("CS0"); + cs0->setParameters( + CPS::Math::singlePhaseVariableToThreePhase(CPS::Math::polar(1.0, 0.0)), + 50.0); + + auto r1 = Ph3::Resistor::make("R1"); + r1->setParameters(10 * param); + + auto r2 = Ph3::Resistor::make("R2"); + r2->setParameters(param); + + auto r3 = Ph3::Resistor::make("R3"); + r3->setParameters(5 * param); + + auto l1 = Ph3::Inductor::make("L1"); + l1->setParameters(0.02 * param); + + auto c1 = Ph3::Capacitor::make("C1"); + c1->setParameters(0.001 * param); + + // Topology + cs0->connect(SimNode::List{n1, SimNode::GND}); + + r1->connect(SimNode::List{n2, n1}); + r2->connect(SimNode::List{n2, SimNode::GND}); + r3->connect(SimNode::List{n2, SimNode::GND}); + + l1->connect(SimNode::List{n2, SimNode::GND}); + + c1->connect(SimNode::List{n1, n2}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2}, + SystemComponentList{cs0, r1, r2, r3, l1, c1}); + + // Logging + Logger::setLogDir("logs/" + simNameRC); + auto logger = DataLogger::make(simNameRC); + logger->logAttribute("V1", n1->attribute("v")); + logger->logAttribute("V2", n2->attribute("v")); + logger->logAttribute("V_C1", c1->attribute("v_intf")); + logger->logAttribute("I_L1", l1->attribute("i_intf")); + + Simulation sim(simNameRC, Logger::Level::info); + sim.setSystem(sys); + sim.addLogger(logger); + sim.setDomain(Domain::EMT); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.run(); +} + +void EMT_PH3_SSN_R3_C1_L1_CS1() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simNameSSN = "EMT_Ph3_R3C1L1CS1_RC_vs_SSN_SSN"; + Logger::setLogDir("logs/" + simNameSSN); + + // Nodes + auto n1 = CPS::EMT::SimNode::make("n1", PhaseType::ABC); + auto n2 = CPS::EMT::SimNode::make("n2", PhaseType::ABC); + + // Components + + Matrix param = Matrix::Zero(3, 3); + param << 1., 0, 0, 0, 1., 0, 0, 0, 1.; + + auto cs0 = Ph3::CurrentSource::make("CS0"); + cs0->setParameters( + CPS::Math::singlePhaseVariableToThreePhase(CPS::Math::polar(1.0, 0.0)), + 50.0); + + auto r1 = Ph3::Resistor::make("R1", Logger::Level::debug); + r1->setParameters(10 * param); + + auto r2 = Ph3::Resistor::make("R2", Logger::Level::debug); + r2->setParameters(param); + + auto r3 = Ph3::Resistor::make("R3", Logger::Level::debug); + r3->setParameters(5 * param); + + auto l1 = Ph3::SSN::Inductor::make("L1"); + l1->setParameters(0.02 * param); + + auto c1 = Ph3::SSN::Capacitor::make("C1"); + c1->setParameters(0.001 * param); + + // Topology + cs0->connect(CPS::EMT::SimNode::List{n1, CPS::EMT::SimNode::GND}); + + r1->connect(CPS::EMT::SimNode::List{n2, n1}); + r2->connect(CPS::EMT::SimNode::List{n2, CPS::EMT::SimNode::GND}); + r3->connect(CPS::EMT::SimNode::List{n2, CPS::EMT::SimNode::GND}); + + l1->connect(CPS::EMT::SimNode::List{n2, CPS::EMT::SimNode::GND}); + + c1->connect(CPS::EMT::SimNode::List{n1, n2}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2}, + SystemComponentList{cs0, r1, r2, r3, l1, c1}); + + // Logging + auto logger = DataLogger::make(simNameSSN); + logger->logAttribute("V1_SSN", n1->attribute("v")); + logger->logAttribute("V2_SSN", n2->attribute("v")); + logger->logAttribute("V_C1_SSN", c1->attribute("v_intf")); + logger->logAttribute("I_L1_SSN", l1->attribute("i_intf")); + + Simulation sim(simNameSSN, Logger::Level::info); + sim.setSystem(sys); + sim.addLogger(logger); + sim.setDomain(Domain::EMT); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.run(); +} + +int main(int argc, char *argv[]) { + EMT_PH3_R3_C1_L1_CS1(); + EMT_PH3_SSN_R3_C1_L1_CS1(); +} diff --git a/dpsim/examples/cxx/Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp b/dpsim/examples/cxx/Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp new file mode 100644 index 0000000000..60918fd37f --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp @@ -0,0 +1,114 @@ +#include + +using namespace DPsim; +using namespace CPS::EMT; + +void EMT_Ph3_RLC1_VS1() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simNameRC = "EMT_Ph3_RLC1VS1_RC_vs_SSN_RC"; + + // Nodes + auto n1 = SimNode::make("n1", PhaseType::ABC); + auto n2 = SimNode::make("n2", PhaseType::ABC); + auto n3 = SimNode::make("n3", PhaseType::ABC); + + // Components + + Matrix param = Matrix::Zero(3, 3); + param << 1., 0, 0, 0, 1., 0, 0, 0, 1.; + + auto vs0 = Ph3::VoltageSource::make("VS0"); + vs0->setParameters( + CPS::Math::singlePhaseVariableToThreePhase(CPS::Math::polar(1.0, 0.0)), + 50.0); + + auto r1 = Ph3::Resistor::make("R1"); + r1->setParameters(1. * param); + + auto l1 = Ph3::Inductor::make("L1"); + l1->setParameters(0.05 * param); + + auto c1 = Ph3::Capacitor::make("C1"); + c1->setParameters(0.01 * param); + + // Topology + vs0->connect(SimNode::List{n1, SimNode::GND}); + + r1->connect(SimNode::List{n1, n2}); + + l1->connect(SimNode::List{n2, n3}); + + c1->connect(SimNode::List{n3, SimNode::GND}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2, n3}, + SystemComponentList{vs0, r1, l1, c1}); + + // Logging + Logger::setLogDir("logs/" + simNameRC); + auto logger = DataLogger::make(simNameRC); + logger->logAttribute("I_R", r1->attribute("i_intf")); + logger->logAttribute("V1_RC", n1->attribute("v")); + + Simulation sim(simNameRC, Logger::Level::info); + sim.setSystem(sys); + sim.addLogger(logger); + sim.setDomain(Domain::EMT); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.run(); +} + +void EMT_Ph3_SSN_RLC1_VS1() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simNameSSN = "EMT_Ph3_RLC1VS1_RC_vs_SSN_SSN"; + + // Nodes + auto n1 = CPS::EMT::SimNode::make("n1", PhaseType::ABC); + + // Components + + Matrix param = Matrix::Zero(3, 3); + param << 1., 0, 0, 0, 1., 0, 0, 0, 1.; + + auto vs0 = Ph3::VoltageSource::make("VS0"); + vs0->setParameters( + CPS::Math::singlePhaseVariableToThreePhase(CPS::Math::polar(1.0, 0.0)), + 50.0); + + auto rlc = Ph3::SSN::Full_Serial_RLC::make("RLC"); + rlc->setParameters(1. * param, 0.05 * param, 0.01 * param); + + // Topology + vs0->connect(CPS::EMT::SimNode::List{n1, CPS::EMT::SimNode::GND}); + + rlc->connect(CPS::EMT::SimNode::List{n1, CPS::EMT::SimNode::GND}); + + // Define system topology + auto sys = + SystemTopology(50, SystemNodeList{n1}, SystemComponentList{vs0, rlc}); + + // Logging + Logger::setLogDir("logs/" + simNameSSN); + auto logger = DataLogger::make(simNameSSN); + logger->logAttribute("I_RLC_SSN", rlc->attribute("i_intf")); + logger->logAttribute("V1_SSN", n1->attribute("v")); + + Simulation sim(simNameSSN, Logger::Level::info); + sim.setSystem(sys); + sim.addLogger(logger); + sim.setDomain(Domain::EMT); + sim.setSolverType(Solver::Type::MNA); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.run(); +} + +int main(int argc, char *argv[]) { + EMT_Ph3_RLC1_VS1(); + EMT_Ph3_SSN_RLC1_VS1(); +} diff --git a/examples/Notebooks/Circuits/EMT_Ph3_compare_RC_SSN.ipynb b/examples/Notebooks/Circuits/EMT_Ph3_compare_RC_SSN.ipynb new file mode 100644 index 0000000000..4c762e70d5 --- /dev/null +++ b/examples/Notebooks/Circuits/EMT_Ph3_compare_RC_SSN.ipynb @@ -0,0 +1,245 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Comparing RC and SSN linear circuit simulations in EMT Ph3 domain\n", + "### Comparing EMT domain simulations of Ph3 linear circuits built from Resistive Companion (RC) component models against the same linear circuits build from SSN (State Space Nodal) component models." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run C++ examples: R3_C1_L1_CS1 circuit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import subprocess\n", + "\n", + "#%matplotlib widget\n", + "\n", + "name = 'EMT_Ph3_R3C1L1CS1_RC_vs_SSN'\n", + "\n", + "dpsim_path = subprocess.Popen(['git', 'rev-parse', '--show-toplevel'], stdout=subprocess.PIPE).communicate()[0].rstrip().decode('utf-8')\n", + "\n", + "path_exec = dpsim_path + '/build/dpsim/examples/cxx/'\n", + "sim = subprocess.Popen([path_exec + name], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", + "print(sim.communicate()[0].decode())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from villas.dataprocessing.readtools import *\n", + "from villas.dataprocessing.timeseries import *\n", + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import matplotlib.pyplot as plt\n", + "import re\n", + "import numpy as np\n", + "import math\n", + "\n", + "work_dir = os.getcwd() + \"/logs/\"\n", + "path_logfile_RC = work_dir + name + '_RC/' + name + '_RC' + '.csv'\n", + "ts_EMT_Ph3_R3C1L1CS1_RC = read_timeseries_dpsim(path_logfile_RC)\n", + "path_logfile_SSN = work_dir + name + '_SSN/' + name + '_SSN' + '.csv'\n", + "ts_EMT_Ph3_R3C1L1CS1_SSN = read_timeseries_dpsim(path_logfile_SSN)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run C++ examples: RLC_VS circuit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "name = 'EMT_Ph3_RLC1VS1_RC_vs_SSN'\n", + "\n", + "dpsim_path = subprocess.Popen(['git', 'rev-parse', '--show-toplevel'], stdout=subprocess.PIPE).communicate()[0].rstrip().decode('utf-8')\n", + "\n", + "path_exec = dpsim_path + '/build/dpsim/examples/cxx/'\n", + "sim = subprocess.Popen([path_exec + name], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", + "print(sim.communicate()[0].decode())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "work_dir = os.getcwd() + \"/logs/\"\n", + "path_logfile_RC = work_dir + name + '_RC/' + name + '_RC' + '.csv'\n", + "ts_EMT_Ph3_RLC1VS1_RC = read_timeseries_dpsim(path_logfile_RC)\n", + "path_logfile_SSN = work_dir + name + '_SSN/' + name + '_SSN' + '.csv'\n", + "ts_EMT_Ph3_RLC1VS1_SSN = read_timeseries_dpsim(path_logfile_SSN)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot R3_C1_L1_CS1 circuit results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.close('all')\n", + "fig1 = plt.figure()\n", + "\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_0'].time, ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_0'].values, \"r-\", label='I_L1_a')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_1'].time, ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_1'].values, \"g-\", label='I_L1_b')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_2'].time, ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_2'].values, \"b-\", label='I_L1_c')\n", + "\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_0'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_0'].values, \"rx\", markevery=10, label='I_L1_SSN_a')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_1'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_1'].values, \"gx\", markevery=10, label='I_L1_SSN_b')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_2'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_2'].values, \"bx\", markevery=10, label='I_L1_SSN_c')\n", + "\n", + "plt.legend(loc = 4)\n", + "\n", + "plt.title('Comparison of resistive companion and SSN simulation: Inductor current')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Phase current [A]')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig2 = plt.figure()\n", + "\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_0'].time, ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_0'].values, \"r-\", label='V_C1_a')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_1'].time, ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_1'].values, \"g-\", label='V_C1_b')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_2'].time, ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_2'].values, \"b-\", label='V_C1_c')\n", + "\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_0'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_0'].values, \"rx\", markevery=10, label='V_C1_SSN_a')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_1'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_1'].values, \"gx\", markevery=10, label='V_C1_SSN_b')\n", + "plt.plot(ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_2'].time, ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_2'].values, \"bx\", markevery=10, label='V_C1_SSN_c')\n", + "\n", + "plt.legend(loc = 4)\n", + "\n", + "plt.title('Comparison of resistive companion and SSN simulation: Capacitor voltage')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Phase voltage [V]')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot RLC_VS circuit results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig3 = plt.figure()\n", + "\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_RC['I_R_0'].time, ts_EMT_Ph3_RLC1VS1_RC['I_R_0'].values, \"r-\", label='I_R_a')\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_RC['I_R_1'].time, ts_EMT_Ph3_RLC1VS1_RC['I_R_1'].values, \"g-\", label='I_R_b')\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_RC['I_R_2'].time, ts_EMT_Ph3_RLC1VS1_RC['I_R_2'].values, \"b-\", label='I_R_c')\n", + "\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_0'].time, ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_0'].values, \"rx\", markevery=10, label='I_RLC_SSN_a')\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_1'].time, ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_1'].values, \"gx\", markevery=10, label='I_RLC_SSN_b')\n", + "plt.plot(ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_2'].time, ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_2'].values, \"bx\", markevery=10, label='I_RLC_SSN_c')\n", + "\n", + "plt.legend(loc = 4)\n", + "\n", + "plt.title('Comparison of resistive companion and SSN simulation results: RLC current')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Phase currents [A]')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Assert" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "epsilon = 1e-100\n", + "\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_0'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_0'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_1'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_1'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['I_L1_2'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['I_L1_SSN_2'].values) < epsilon)\n", + "\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_0'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_0'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_1'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_1'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_R3C1L1CS1_RC['V_C1_2'].values-ts_EMT_Ph3_R3C1L1CS1_SSN['V_C1_SSN_2'].values) < epsilon)\n", + "\n", + "assert(np.max(ts_EMT_Ph3_RLC1VS1_RC['I_R_0'].values-ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_0'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_RLC1VS1_RC['I_R_1'].values-ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_1'].values) < epsilon)\n", + "assert(np.max(ts_EMT_Ph3_RLC1VS1_RC['I_R_2'].values-ts_EMT_Ph3_RLC1VS1_SSN['I_RLC_SSN_2'].values) < epsilon)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + }, + "tests": { + "skip": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}