diff --git a/CMakeLists.txt b/CMakeLists.txt index 55bbcc8466..48252b92d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -59,6 +59,8 @@ if(MSVC) # Silence Visual Studio error C2220 add_definitions(-D_SILENCE_STDEXT_ARR_ITERS_DEPRECATION_WARNING) + + add_compile_options(/bigobj) # Fixes error C1128 for compiling MNAEigenvalueExtractor with MSVC compiler. # Set exception handling for portability set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /EHsc") diff --git a/docs/hugo/content/en/docs/Overview/eigenvalues.md b/docs/hugo/content/en/docs/Overview/eigenvalues.md new file mode 100644 index 0000000000..6699598435 --- /dev/null +++ b/docs/hugo/content/en/docs/Overview/eigenvalues.md @@ -0,0 +1,43 @@ +--- +title: "Eigenvalues" +linkTitle: "Eigenvalues" +date: 2024-03-20 +--- + +In parallel to simulating a power system, DPsim allows to extract eigenvalues of the power system state-space model's state matrix for each simulation time step. Eigenvalue extraction can be enabled for a ``Simulation``. + + + +## Equations +``discreteEigenvalues`` $\mathbf{z}$ are computed from discretized state matrix $A_{d}$, that in EMT domain is calculated from: + +$$ +\mathbf{A}_{d}=\mathbf{S}-\mathbf{G}_{b}\mathbf{A}_{bn}\mathbf{G}^{-1} \mathbf{A}_{nb} +$$ + +where $\mathbf{S}$ is a sign matrix, $\mathbf{G}_{b}$ is a discretization matrix, $\mathbf{G}$ is a Modified Nodal Analysis (MNA) power system conductance matrix, $\mathbf{A}_{bn}$ and $\mathbf{A}_{nb}$ are branch x node and node x branch incidence matrices respectively. + +The MNA power system conductance matrix $\mathbf{G}$ is available from power system MNA model. To prepare the rest of the matrices, each power system component needs to be stamped into $\mathbf{A}_{bn}$ and $\mathbf{A}_{nb}$, while dynamic components also need to be stamped into $\mathbf{S}$ and $\mathbf{G}_{b}$ matrices. + + +``eigenvalues`` $\mathbf{\lambda}$ of the time-continuous state-space model matrix $A$ can then be calculated from ``discreteEigenvalues`` $\mathbf{z}$. Assuming the Trapezoidal rule of discretization in EMT domain, the equation is: + +$$ +\mathbf{\lambda}=\frac{2}{\Delta t} \frac{\mathbf{z} - 1}{\mathbf{z} + 1} +$$ + + +## Implementation + +The ``MNAEigenvalueExtractor`` class is a template class responsible for extracting eigenvalues. +The ``EigenvalueCompInterface`` is interface to be implemented by all components participating in eigenvalue extraction. The ``EigenvalueDynamicCompInterface`` is interface to be implemented by all _dynamic_ components participating in eigenvalue extraction. These interfaces provide the signatures of the functions that are used to stamp a component into the matrices $\mathbf{A}_{bn}$, $\mathbf{A}_{nb}$, $\mathbf{S}$ and $\mathbf{G}_{b}$. + +The table below provides an overview of the components, that support eigenvalue extraction in EMT and in DP domains: + +| Component |Is dynamic?| EMT domain, Ph1 | DP domain, Ph1 | +| -------- |-------- | -------- | -------- | +| ``Resistor`` |—| ✔ | ✔ | +| ``Inductor`` |✔| ✔ | ✔ | +| ``Capacitor`` |✔| ✔ | ✔ | +| ``Switch`` |—| — | ✔ | +| ``VoltageSource`` |—| ✔ | ✔ | diff --git a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Capacitor.h b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Capacitor.h index d98fff531d..0028a8c100 100644 --- a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Capacitor.h +++ b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Capacitor.h @@ -10,6 +10,7 @@ #include #include +#include namespace CPS { namespace DP { @@ -23,7 +24,8 @@ namespace Ph1 { /// frequency and the current source changes for each iteration. class Capacitor : public MNASimPowerComp, public Base::Ph1::Capacitor, - public SharedFactory { + public SharedFactory, + public EigenvalueDynamicCompInterface { protected: /// DC equivalent current source for harmonics [A] MatrixComp mEquivCurrent; @@ -119,6 +121,15 @@ class Capacitor : public MNASimPowerComp, Capacitor &mCapacitor; std::vector::Ptr> mLeftVectors; }; + + // #### Implementation of eigenvalue dynamic component interface #### + void stampSignMatrix(UInt branchIdx, MatrixVar &signMatrix, + Complex coeffDP) final; + void stampDiscretizationMatrix(UInt branchIdx, + MatrixVar &discretizationMatrix, + Complex coeffDP) final; + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace DP diff --git a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Inductor.h b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Inductor.h index 131a56b1b7..5e713209ab 100644 --- a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Inductor.h +++ b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Inductor.h @@ -10,6 +10,7 @@ #include #include +#include #include namespace CPS { @@ -25,7 +26,8 @@ namespace Ph1 { class Inductor : public MNASimPowerComp, public Base::Ph1::Inductor, public MNATearInterface, - public SharedFactory { + public SharedFactory, + public EigenvalueDynamicCompInterface { protected: /// DC equivalent current source for harmonics [A] MatrixComp mEquivCurrent; @@ -127,6 +129,15 @@ class Inductor : public MNASimPowerComp, Inductor &mInductor; std::vector::Ptr> mLeftVectors; }; + + // #### Implementation of eigenvalue dynamic component interface #### + void stampSignMatrix(UInt branchIdx, MatrixVar &signMatrix, + Complex coeffDP) final; + void stampDiscretizationMatrix(UInt branchIdx, + MatrixVar &discretizationMatrix, + Complex coeffDP) final; + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace DP diff --git a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Resistor.h b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Resistor.h index a057c5bcf5..68726c830c 100644 --- a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Resistor.h +++ b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Resistor.h @@ -11,6 +11,7 @@ #include #include #include +#include #include namespace CPS { @@ -21,7 +22,8 @@ class Resistor : public MNASimPowerComp, public Base::Ph1::Resistor, public MNATearInterface, public DAEInterface, - public SharedFactory { + public SharedFactory, + public EigenvalueCompInterface { public: /// Defines UID, name and logging level Resistor(String uid, String name, @@ -30,37 +32,38 @@ class Resistor : public MNASimPowerComp, Resistor(String name, Logger::Level logLevel = Logger::Level::off) : Resistor(name, name, logLevel) {} - SimPowerComp::Ptr clone(String name); + SimPowerComp::Ptr clone(String name) override; // #### General #### /// Initializes component from power flow data - void initializeFromNodesAndTerminals(Real frequency); + void initializeFromNodesAndTerminals(Real frequency) override; // #### MNA section #### void mnaCompInitialize(Real omega, Real timeStep, - Attribute::Ptr leftVector); - void mnaCompInitializeHarm(Real omega, Real timeStep, - std::vector::Ptr> leftVector); + Attribute::Ptr leftVector) override; + void mnaCompInitializeHarm( + Real omega, Real timeStep, + std::vector::Ptr> leftVector) override; /// Stamps system matrix - void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix); + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; /// Stamps system matrix considering the frequency index void mnaCompApplySystemMatrixStampHarm(SparseMatrixRow &systemMatrix, - Int freqIdx); + Int freqIdx) override; /// Update interface voltage from MNA system result - void mnaCompUpdateVoltage(const Matrix &leftVector); + void mnaCompUpdateVoltage(const Matrix &leftVector) override; void mnaCompUpdateVoltageHarm(const Matrix &leftVector, Int freqIdx); /// Update interface current from MNA system result - void mnaCompUpdateCurrent(const Matrix &leftVector); + void mnaCompUpdateCurrent(const Matrix &leftVector) override; void mnaCompUpdateCurrentHarm(); /// MNA pre and post step operations void mnaCompPostStep(Real time, Int timeStepCount, - Attribute::Ptr &leftVector); + Attribute::Ptr &leftVector) override; /// add MNA pre and post step dependencies void mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, - Attribute::Ptr &leftVector); + Attribute::Ptr &leftVector) override; class MnaPostStepHarm : public Task { public: @@ -86,9 +89,13 @@ class Resistor : public MNASimPowerComp, // #### DAE Section #### ///Residual Function for DAE Solver void daeResidual(double ttime, const double state[], const double dstate_dt[], - double resid[], std::vector &off); + double resid[], std::vector &off) override; ///Voltage Getter - Complex daeInitialize(); + Complex daeInitialize() override; + + // #### Implementation of eigenvalue component interface #### + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace DP diff --git a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Switch.h b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Switch.h index 268d8798d4..47bff17aa9 100644 --- a/dpsim-models/include/dpsim-models/DP/DP_Ph1_Switch.h +++ b/dpsim-models/include/dpsim-models/DP/DP_Ph1_Switch.h @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -25,7 +26,8 @@ namespace Ph1 { class Switch : public MNASimPowerComp, public Base::Ph1::Switch, public SharedFactory, - public MNASwitchInterface { + public MNASwitchInterface, + public EigenvalueCompInterface { protected: public: /// Defines UID, name, component parameters and logging level @@ -34,40 +36,44 @@ class Switch : public MNASimPowerComp, Switch(String name, Logger::Level logLevel = Logger::Level::off) : Switch(name, name, logLevel) {} - SimPowerComp::Ptr clone(String name); + SimPowerComp::Ptr clone(String name) override; // #### General #### /// Initializes component from power flow data - void initializeFromNodesAndTerminals(Real frequency); + void initializeFromNodesAndTerminals(Real frequency) override; // #### General MNA section #### void mnaCompInitialize(Real omega, Real timeStep, - Attribute::Ptr leftVector); + Attribute::Ptr leftVector) override; /// Stamps system matrix - void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix); + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; /// Stamps right side (source) vector - void mnaCompApplyRightSideVectorStamp(Matrix &rightVector); + void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override; /// Update interface voltage from MNA system result - void mnaCompUpdateVoltage(const Matrix &leftVector); + void mnaCompUpdateVoltage(const Matrix &leftVector) override; /// Update interface current from MNA system result - void mnaCompUpdateCurrent(const Matrix &leftVector); + void mnaCompUpdateCurrent(const Matrix &leftVector) override; /// MNA post step operations void mnaCompPostStep(Real time, Int timeStepCount, - Attribute::Ptr &leftVector); + Attribute::Ptr &leftVector) override; /// Add MNA post step dependencies void mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, - Attribute::Ptr &leftVector); + Attribute::Ptr &leftVector) override; // #### MNA section for switch #### /// Check if switch is closed - Bool mnaIsClosed(); + Bool mnaIsClosed() override; /// Stamps system matrix considering the defined switch position void mnaCompApplySwitchSystemMatrixStamp(Bool closed, SparseMatrixRow &systemMatrix, - Int freqIdx); + Int freqIdx) override; + + // #### Implementation of eigenvalue component interface #### + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace DP diff --git a/dpsim-models/include/dpsim-models/DP/DP_Ph1_VoltageSource.h b/dpsim-models/include/dpsim-models/DP/DP_Ph1_VoltageSource.h index 84f3dc204e..72bbc048e5 100644 --- a/dpsim-models/include/dpsim-models/DP/DP_Ph1_VoltageSource.h +++ b/dpsim-models/include/dpsim-models/DP/DP_Ph1_VoltageSource.h @@ -14,6 +14,7 @@ #include #include #include +#include #include namespace CPS { @@ -38,7 +39,8 @@ namespace Ph1 { /// a new equation ej - ek = V is added to the problem. class VoltageSource : public MNASimPowerComp, public DAEInterface, - public SharedFactory { + public SharedFactory, + public EigenvalueCompInterface { private: /// void updateVoltage(Real time); @@ -151,6 +153,10 @@ class VoltageSource : public MNASimPowerComp, double resid[], std::vector &off) override; ///Voltage Getter Complex daeInitialize() override; + + // #### Implementation of eigenvalue component interface #### + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace DP diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Capacitor.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Capacitor.h index 51bb01bf5c..e7cc873577 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Capacitor.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Capacitor.h @@ -10,6 +10,7 @@ #include #include +#include #include namespace CPS { @@ -24,7 +25,8 @@ namespace Ph1 { ///frequency and the current source changes for each iteration. class Capacitor : public MNASimPowerComp, public Base::Ph1::Capacitor, - public SharedFactory { + public SharedFactory, + public EigenvalueDynamicCompInterface { protected: /// DC equivalent current source [A] Real mEquivCurrent; @@ -73,6 +75,15 @@ class Capacitor : public MNASimPowerComp, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + + // #### Implementation of eigenvalue dynamic component interface #### + void stampSignMatrix(UInt branchIdx, MatrixVar &signMatrix, + Complex coeffDP) final; + void stampDiscretizationMatrix(UInt branchIdx, + MatrixVar &discretizationMatrix, + Complex coeffDP) final; + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace EMT diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Inductor.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Inductor.h index d87769454c..f9d5744fe7 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Inductor.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Inductor.h @@ -10,6 +10,7 @@ #include #include +#include #include namespace CPS { @@ -24,7 +25,8 @@ namespace Ph1 { /// frequency and the current source changes for each iteration. class Inductor : public MNASimPowerComp, public Base::Ph1::Inductor, - public SharedFactory { + public SharedFactory, + public EigenvalueDynamicCompInterface { protected: /// DC equivalent current source [A] Real mEquivCurrent; @@ -74,6 +76,15 @@ class Inductor : public MNASimPowerComp, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + + // #### Implementation of eigenvalue dynamic component interface #### + void stampSignMatrix(UInt branchIdx, MatrixVar &signMatrix, + Complex coeffDP) final; + void stampDiscretizationMatrix(UInt branchIdx, + MatrixVar &discretizationMatrix, + Complex coeffDP) final; + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace EMT diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Resistor.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Resistor.h index bb7f2947ab..23ac90cfb9 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Resistor.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Resistor.h @@ -12,6 +12,7 @@ #include #include +#include #include namespace CPS { @@ -20,7 +21,8 @@ namespace Ph1 { /// EMT Resistor class Resistor : public MNASimPowerComp, public Base::Ph1::Resistor, - public SharedFactory { + public SharedFactory, + public EigenvalueCompInterface { protected: public: /// Defines UID, name, component parameters and logging level @@ -56,6 +58,10 @@ class Resistor : public MNASimPowerComp, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + + // #### Implementation of eigenvalue component interface #### + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace EMT diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h index c687af334d..57b69bcc2c 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h @@ -9,6 +9,7 @@ #pragma once #include +#include #include namespace CPS { @@ -23,7 +24,8 @@ namespace Ph1 { /// of node k as negative. Moreover /// a new equation ej - ek = V is added to the problem. class VoltageSource : public MNASimPowerComp, - public SharedFactory { + public SharedFactory, + public EigenvalueCompInterface { private: Real mTimeStep; @@ -74,6 +76,9 @@ class VoltageSource : public MNASimPowerComp, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + // #### Implementation of eigenvalue component interface #### + void stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) final; }; } // namespace Ph1 } // namespace EMT diff --git a/dpsim-models/include/dpsim-models/Logger.h b/dpsim-models/include/dpsim-models/Logger.h index ac42949ee3..5e866d2cd3 100644 --- a/dpsim-models/include/dpsim-models/Logger.h +++ b/dpsim-models/include/dpsim-models/Logger.h @@ -54,6 +54,8 @@ class Logger { std::string pattern); // #### to string methods #### + template + static String matrixVarToString(const MatrixVar &mat); static String matrixToString(const Matrix &mat); static String matrixCompToString(const MatrixComp &mat); static String sparseMatrixToString(const SparseMatrix &mat); diff --git a/dpsim-models/include/dpsim-models/MathUtils.h b/dpsim-models/include/dpsim-models/MathUtils.h index 6841df8ebc..0ffdb38e6a 100644 --- a/dpsim-models/include/dpsim-models/MathUtils.h +++ b/dpsim-models/include/dpsim-models/MathUtils.h @@ -56,12 +56,13 @@ class Math { static Real realFromVectorElement(const Matrix &mat, Matrix::Index row); - // #### Matric Operations #### + // #### Matrix Operations #### // // | Re-Re(row,col)_harm1 | Im-Re(row,col)_harm1 | Interharmonics harm1-harm2 // | Re-Im(row,col)_harm1 | Im-Im(row,col)_harm1 | Interharmonics harm1-harm2 // | Interharmonics harm1-harm2 | Re(row,col)_harm2 | Re(row,col)_harm2 | // | Interharmonics harm1-harm2 | Im(row,col)_harm2 | Im(row,col)_harm2 | + static void setMatrixElement(SparseMatrixRow &mat, Matrix::Index row, Matrix::Index column, Complex value, Int maxFreq = 1, Int freqIdx = 0); diff --git a/dpsim-models/include/dpsim-models/Solver/EigenvalueCompInterface.h b/dpsim-models/include/dpsim-models/Solver/EigenvalueCompInterface.h new file mode 100644 index 0000000000..815cb5369c --- /dev/null +++ b/dpsim-models/include/dpsim-models/Solver/EigenvalueCompInterface.h @@ -0,0 +1,14 @@ +#pragma once +#include + +namespace CPS { +/// Interface to be implemented by all components taking part in eigenvalue extraction. +class EigenvalueCompInterface { +public: + using Ptr = std::shared_ptr; + /// Stamp component into branch<->node incidence matrix used for eigenvalue extraction. + virtual void + stampBranchNodeIncidenceMatrix(UInt branchIdx, + Matrix &branchNodeIncidenceMatrix) = 0; +}; +} // namespace CPS \ No newline at end of file diff --git a/dpsim-models/include/dpsim-models/Solver/EigenvalueDynamicCompInterface.h b/dpsim-models/include/dpsim-models/Solver/EigenvalueDynamicCompInterface.h new file mode 100644 index 0000000000..6c4e6e8246 --- /dev/null +++ b/dpsim-models/include/dpsim-models/Solver/EigenvalueDynamicCompInterface.h @@ -0,0 +1,21 @@ +#pragma once +#include + +namespace CPS { +template +/// Interface to be implemented by all dynamic components taking part in eigenvalue extraction. +class EigenvalueDynamicCompInterface : public EigenvalueCompInterface { +public: + using Ptr = std::shared_ptr; + using List = std::vector; + + /// Stamp component into sign matrix used for eigenvalue extraction. + virtual void stampSignMatrix(UInt branchIdx, MatrixVar &signMatrix, + Complex coeffDP) = 0; + /// Stamp component into discretization matrix used for eigenvalue extraction. + virtual void + stampDiscretizationMatrix(UInt branchIdx, + MatrixVar &discretizationMatrix, + Complex coeffDP) = 0; +}; +} // namespace CPS \ No newline at end of file diff --git a/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp b/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp index 08f0195da4..605eae0c4e 100644 --- a/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp +++ b/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp @@ -295,3 +295,25 @@ void DP::Ph1::Capacitor::mnaCompUpdateCurrentHarm() { Logger::phasorToString((**mIntfCurrent)(0, freq))); } } + +void DP::Ph1::Capacitor::stampSignMatrix(UInt branchIdx, + MatrixVar &signMatrix, + Complex coeffDP) { + signMatrix(branchIdx, branchIdx) = -1.0; +} + +void DP::Ph1::Capacitor::stampDiscretizationMatrix( + UInt branchIdx, MatrixVar &discretizationMatrix, Complex coeffDP) { + discretizationMatrix(branchIdx, branchIdx) = + -mEquivCond(0, 0) * (1.0 + coeffDP); +} + +void DP::Ph1::Capacitor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/DP/DP_Ph1_Inductor.cpp b/dpsim-models/src/DP/DP_Ph1_Inductor.cpp index 48fa904970..5772502b15 100644 --- a/dpsim-models/src/DP/DP_Ph1_Inductor.cpp +++ b/dpsim-models/src/DP/DP_Ph1_Inductor.cpp @@ -284,3 +284,25 @@ void DP::Ph1::Inductor::mnaTearPostStep(Complex voltage, Complex current) { (**mIntfVoltage)(0, 0) = voltage; (**mIntfCurrent)(0, 0) = mEquivCond(0, 0) * voltage + mEquivCurrent(0, 0); } + +void DP::Ph1::Inductor::stampSignMatrix(UInt branchIdx, + MatrixVar &signMatrix, + Complex coeffDP) { + signMatrix(branchIdx, branchIdx) = coeffDP; +} + +void DP::Ph1::Inductor::stampDiscretizationMatrix( + UInt branchIdx, MatrixVar &discretizationMatrix, Complex coeffDP) { + discretizationMatrix(branchIdx, branchIdx) = + mEquivCond(0, 0) * (1.0 + coeffDP); +} + +void DP::Ph1::Inductor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/DP/DP_Ph1_Resistor.cpp b/dpsim-models/src/DP/DP_Ph1_Resistor.cpp index 3ef1f4fb28..bd2ef0eb90 100644 --- a/dpsim-models/src/DP/DP_Ph1_Resistor.cpp +++ b/dpsim-models/src/DP/DP_Ph1_Resistor.cpp @@ -205,3 +205,13 @@ void DP::Ph1::Resistor::daeResidual(double ttime, const double state[], } Complex DP::Ph1::Resistor::daeInitialize() { return (**mIntfVoltage)(0, 0); } + +void DP::Ph1::Resistor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/DP/DP_Ph1_Switch.cpp b/dpsim-models/src/DP/DP_Ph1_Switch.cpp index 6f7f24d35d..b113d3a108 100644 --- a/dpsim-models/src/DP/DP_Ph1_Switch.cpp +++ b/dpsim-models/src/DP/DP_Ph1_Switch.cpp @@ -107,3 +107,13 @@ void DP::Ph1::Switch::mnaCompPostStep(Real time, Int timeStepCount, mnaCompUpdateVoltage(**leftVector); mnaCompUpdateCurrent(**leftVector); } + +void DP::Ph1::Switch::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/DP/DP_Ph1_VoltageSource.cpp b/dpsim-models/src/DP/DP_Ph1_VoltageSource.cpp index 0e06e04737..804817c5b5 100644 --- a/dpsim-models/src/DP/DP_Ph1_VoltageSource.cpp +++ b/dpsim-models/src/DP/DP_Ph1_VoltageSource.cpp @@ -307,3 +307,13 @@ Complex DP::Ph1::VoltageSource::daeInitialize() { (**mIntfVoltage)(0, 0) = mSrcSig->getSignal(); return mSrcSig->getSignal(); } + +void DP::Ph1::VoltageSource::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp index 5fcbba1961..2e1d6546ef 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp @@ -119,3 +119,24 @@ void EMT::Ph1::Capacitor::mnaCompUpdateVoltage(const Matrix &leftVector) { void EMT::Ph1::Capacitor::mnaCompUpdateCurrent(const Matrix &leftVector) { (**mIntfCurrent)(0, 0) = mEquivCond * (**mIntfVoltage)(0, 0) + mEquivCurrent; } + +void EMT::Ph1::Capacitor::stampSignMatrix(UInt branchIdx, + MatrixVar &signMatrix, + Complex coeffDP) { + signMatrix(branchIdx, branchIdx) = -1.0; +} + +void EMT::Ph1::Capacitor::stampDiscretizationMatrix( + UInt branchIdx, MatrixVar &discretizationMatrix, Complex coeffDP) { + discretizationMatrix(branchIdx, branchIdx) = -2 * mEquivCond; +} + +void EMT::Ph1::Capacitor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp index 63e8a275da..de20b732bc 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp @@ -116,3 +116,24 @@ void EMT::Ph1::Inductor::mnaCompUpdateVoltage(const Matrix &leftVector) { void EMT::Ph1::Inductor::mnaCompUpdateCurrent(const Matrix &leftVector) { (**mIntfCurrent)(0, 0) = mEquivCond * (**mIntfVoltage)(0, 0) + mEquivCurrent; } + +void EMT::Ph1::Inductor::stampSignMatrix(UInt branchIdx, + MatrixVar &signMatrix, + Complex coeffDP) { + signMatrix(branchIdx, branchIdx) = 1.0; +} + +void EMT::Ph1::Inductor::stampDiscretizationMatrix( + UInt branchIdx, MatrixVar &discretizationMatrix, Complex coeffDP) { + discretizationMatrix(branchIdx, branchIdx) = 2 * mEquivCond; +} + +void EMT::Ph1::Inductor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp index a96950151c..ec016b9d95 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp @@ -88,3 +88,13 @@ void EMT::Ph1::Resistor::mnaCompUpdateVoltage(const Matrix &leftVector) { void EMT::Ph1::Resistor::mnaCompUpdateCurrent(const Matrix &leftVector) { (**mIntfCurrent)(0, 0) = (**mIntfVoltage)(0, 0) / **mResistance; } + +void EMT::Ph1::Resistor::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp b/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp index 41688a8137..09b0efabff 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp @@ -84,7 +84,7 @@ void EMT::Ph1::VoltageSource::updateVoltage(Real time) { if (srcFreq > 0) (**mIntfVoltage)(0, 0) = Math::abs(voltageRef) * - cos((time)*2. * PI * srcFreq + Math::phase(voltageRef)); + cos((time) * 2. * PI * srcFreq + Math::phase(voltageRef)); else (**mIntfVoltage)(0, 0) = voltageRef.real(); } @@ -121,3 +121,13 @@ void EMT::Ph1::VoltageSource::mnaCompUpdateCurrent(const Matrix &leftVector) { (**mIntfCurrent)(0, 0) = Math::realFromVectorElement( leftVector, mVirtualNodes[0]->matrixNodeIndex()); } + +void EMT::Ph1::VoltageSource::stampBranchNodeIncidenceMatrix( + UInt branchIdx, Matrix &branchNodeIncidenceMatrix) { + if (terminalNotGrounded(0)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(0)) = 1.0; + } + if (terminalNotGrounded(1)) { + branchNodeIncidenceMatrix(branchIdx, matrixNodeIndex(1)) = -1.0; + } +} diff --git a/dpsim-models/src/Logger.cpp b/dpsim-models/src/Logger.cpp index 62fe3df163..ba53e5b9bb 100644 --- a/dpsim-models/src/Logger.cpp +++ b/dpsim-models/src/Logger.cpp @@ -28,6 +28,16 @@ void Logger::setLogPattern(std::shared_ptr logger, } // #### to string methods #### +template +String Logger::matrixVarToString(const MatrixVar &mat) { + std::stringstream ss; + ss << std::scientific << "\n" << mat; + return ss.str(); +} +template String Logger::matrixVarToString(const MatrixVar &mat); +template String +Logger::matrixVarToString(const MatrixVar &mat); + String Logger::matrixToString(const Matrix &mat) { std::stringstream ss; ss << std::scientific << "\n" << mat; diff --git a/dpsim/examples/cxx/CMakeLists.txt b/dpsim/examples/cxx/CMakeLists.txt index dba9ec1117..6f9eaa1eaa 100644 --- a/dpsim/examples/cxx/CMakeLists.txt +++ b/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_DP_SinglePhaseRLC.cpp # EMT examples with PF initialization Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp diff --git a/dpsim/examples/cxx/Circuits/EMT_DP_SinglePhaseRLC.cpp b/dpsim/examples/cxx/Circuits/EMT_DP_SinglePhaseRLC.cpp new file mode 100644 index 0000000000..4929810865 --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_DP_SinglePhaseRLC.cpp @@ -0,0 +1,127 @@ +#include + +using namespace DPsim; +using namespace CPS; + +void runEMT(Real resistance) { + Real timeStep = 1e-4; + Real finalTime = 1e-3; + String simName = + "EMT_SinglePhaseRLC_resistance_" + std::to_string(resistance) + "_Ohm"; + Logger::setLogDir("logs/" + simName); + + // Nodes + auto n1 = SimNode::make("n1"); + auto n2 = SimNode::make("n2"); + auto n3 = SimNode::make("n3"); + + // Components + auto vs = EMT::Ph1::VoltageSource::make("vs"); + vs->setParameters(Complex(100, 0)); + auto r = EMT::Ph1::Resistor::make("r", Logger::Level::info); + r->setParameters(resistance); + auto l = EMT::Ph1::Inductor::make("l", Logger::Level::info); + l->setParameters(5); + auto c = EMT::Ph1::Capacitor::make("c", Logger::Level::info); + c->setParameters(250e-6); + + // Connections + vs->connect(SimNode::List{SimNode::GND, n3}); + r->connect(SimNode::List{n3, n1}); + l->connect(SimNode::List{n1, n2}); + c->connect(SimNode::List{n2, SimNode::GND}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2, n3}, + SystemComponentList{vs, r, l, c}); + + // Logger + auto logger = DataLogger::make(simName); + logger->logAttribute("v1", n1->attribute("v")); + logger->logAttribute("v2", n2->attribute("v")); + logger->logAttribute("i", r->attribute("i_intf")); + + Simulation sim(simName); + sim.setSystem(sys); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.setDomain(Domain::EMT); + sim.doEigenvalueExtraction(true); + sim.addLogger(logger); + + sim.run(); +} + +void runDP() { + Real timeStep = 1e-4; + Real finalTime = 1e-3; + String simName = "DP_SinglePhaseRLC"; + Logger::setLogDir("logs/" + simName); + + // Nodes + auto n1 = SimNode::make("n1"); + auto n2 = SimNode::make("n2"); + auto n3 = SimNode::make("n3"); + auto n4 = SimNode::make("n2"); + auto n5 = SimNode::make("n3"); + + // Components + auto vs = DP::Ph1::VoltageSource::make("vs"); + vs->setParameters(Complex(100, 0)); + auto r = DP::Ph1::Resistor::make("r", Logger::Level::info); + r->setParameters(100); + auto r2 = DP::Ph1::Resistor::make("r2", Logger::Level::info); + r2->setParameters(1e-3); + auto l = DP::Ph1::Inductor::make("l", Logger::Level::info); + l->setParameters(5); + auto c = DP::Ph1::Capacitor::make("c", Logger::Level::info); + c->setParameters(250e-6); + auto s = DP::Ph1::Switch::make("s", Logger::Level::info); + s->setParameters(1e8, 1e-4, true); + auto s2 = DP::Ph1::Switch::make("s2", Logger::Level::info); + s2->setParameters(1e8, 1e-4, false); + + // Switch events + auto sEvent = SwitchEvent::make(5e-4, s, false); + auto s2Event = SwitchEvent::make(5e-4, s2, true); + + // Connections + vs->connect(SimNode::List{SimNode::GND, n1}); + r->connect(SimNode::List{n4, SimNode::GND}); + r2->connect(SimNode::List{n5, SimNode::GND}); + l->connect(SimNode::List{n1, n2}); + c->connect(SimNode::List{n2, n3}); + s->connect(SimNode::List{n3, n4}); + s2->connect(SimNode::List{n3, n5}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2, n3, n4, n5}, + SystemComponentList{vs, r, r2, l, c, s, s2}); + + // Logger + auto logger = DataLogger::make(simName); + logger->logAttribute("v1", n1->attribute("v")); + logger->logAttribute("v2", n2->attribute("v")); + logger->logAttribute("ir", r->attribute("i_intf")); + logger->logAttribute("ir1", r2->attribute("i_intf")); + + Simulation sim(simName); + sim.setSystem(sys); + sim.addEvent(sEvent); + sim.addEvent(s2Event); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.setDomain(Domain::DP); + sim.doEigenvalueExtraction(true); + sim.addLogger(logger); + + sim.run(); +} + +int main(int argc, char *argv[]) { + runEMT(100); + runEMT(1.1e-3); + runDP(); + + return 0; +} diff --git a/dpsim/include/dpsim/EigenvalueUtils.h b/dpsim/include/dpsim/EigenvalueUtils.h new file mode 100644 index 0000000000..7032e2798a --- /dev/null +++ b/dpsim/include/dpsim/EigenvalueUtils.h @@ -0,0 +1,14 @@ +#pragma once + +#include + +namespace DPsim { +class EigenvalueUtils { +public: + static MatrixComp returnNonZeroElements(const MatrixComp &mat); + + static MatrixComp + convertRealEquivalentToComplexMatrix(const Matrix &realEquivalentMatrix); +}; + +} // namespace DPsim \ No newline at end of file diff --git a/dpsim/include/dpsim/MNAEigenvalueExtractor.h b/dpsim/include/dpsim/MNAEigenvalueExtractor.h new file mode 100644 index 0000000000..236f051559 --- /dev/null +++ b/dpsim/include/dpsim/MNAEigenvalueExtractor.h @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include +#include + +namespace DPsim { +/// Extracts eigenvalues using MNA power system conductance matrix +template class MNAEigenvalueExtractor { +public: + MNAEigenvalueExtractor(CPS::Logger::Level logLevel); + + void initialize(const CPS::SystemTopology &topology, + UInt numMatrixNodeIndices, Real timeStep); + void extractEigenvalues(const Matrix &powerSystemMatrix, Real time, + Int timeStepCount); + void closeLogger(); + +private: + std::map + mEigenvalueComponentToBranchIdx; + typename CPS::EigenvalueDynamicCompInterface::List + mEigenvalueDynamicComponents; + Real mTimeStep; + Real mSystemOmega; + Complex mCoeffDP; + MatrixVar mSignMatrix; + MatrixVar mDiscretizationMatrix; + Matrix mBranchNodeIncidenceMatrix; + Matrix mNodeBranchIncidenceMatrix; + MatrixVar mStateMatrix; + CPS::AttributeStatic::Ptr mDiscreteEigenvalues; + CPS::AttributeStatic::Ptr mEigenvalues; + MNAEigenvalueLogger mLogger; + + void setParameters(const CPS::SystemTopology &topology, Real timeStep); + void + identifyEigenvalueComponents(const CPS::IdentifiedObject::List &components); + void createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices); + void stampEigenvalueMatrices(); + void calculateStateMatrix(const Matrix &powerSystemMatrix); + void computeDiscreteEigenvalues(); + void recoverEigenvalues(); +}; +} // namespace DPsim diff --git a/dpsim/include/dpsim/MNAEigenvalueLogger.h b/dpsim/include/dpsim/MNAEigenvalueLogger.h new file mode 100644 index 0000000000..c8d6108cec --- /dev/null +++ b/dpsim/include/dpsim/MNAEigenvalueLogger.h @@ -0,0 +1,29 @@ +#pragma once + +#include + +namespace DPsim { +/// Logs eigenvalues during simulation +class MNAEigenvalueLogger { +public: + MNAEigenvalueLogger(String name, CPS::Logger::Level logLevel); + + void + setLogAttributes(CPS::AttributeStatic::Ptr eigenvalues, + CPS::AttributeStatic::Ptr discreteEigenvalues); + template + void logInitialization(const MatrixVar &signMatrix, + const MatrixVar &discretizationMatrix, + const Matrix &branchNodeIncidenceMatrix, + const Matrix &nodeBranchIncidenceMatrix); + template + void logExtraction(Real time, Int timeStepCount, + const MatrixVar &stateMatrix); + void close(); + +private: + CPS::Logger::Log mSLog; + DataLogger mEigenvaluesLogger; + DataLogger mDiscreteEigenvaluesLogger; +}; +} // namespace DPsim diff --git a/dpsim/include/dpsim/MNASolver.h b/dpsim/include/dpsim/MNASolver.h index 9967523098..a0ffdeec9b 100644 --- a/dpsim/include/dpsim/MNASolver.h +++ b/dpsim/include/dpsim/MNASolver.h @@ -23,6 +23,7 @@ #include #include #include +#include #include /* std::size_t is the largest data type. No container can store @@ -120,6 +121,9 @@ template class MnaSolver : public Solver { /// LU refactorization measurements std::vector mRecomputationTimes; + // #### Eigenvalue extraction #### + MNAEigenvalueExtractor mMNAEigenvalueExtractor; + /// Constructor should not be called by users but by Simulation MnaSolver(String name, CPS::Domain domain = CPS::Domain::DP, CPS::Logger::Level logLevel = CPS::Logger::Level::info); @@ -181,6 +185,8 @@ template class MnaSolver : public Solver { virtual std::shared_ptr createLogTask() = 0; /// Create a solve task for this solver implementation virtual std::shared_ptr createSolveTaskHarm(UInt freqIdx) = 0; + /// Create a task for eigenvalues extraction + virtual std::shared_ptr createExtractEigenvaluesTask() = 0; // #### Scheduler Task Methods #### /// Solves system for single frequency diff --git a/dpsim/include/dpsim/MNASolverDirect.h b/dpsim/include/dpsim/MNASolverDirect.h index cc1be7b62c..9cabe816d7 100644 --- a/dpsim/include/dpsim/MNASolverDirect.h +++ b/dpsim/include/dpsim/MNASolverDirect.h @@ -136,6 +136,8 @@ template class MnaSolverDirect : public MnaSolver { std::shared_ptr createLogTask() override; /// Create a solve task for this solver implementation std::shared_ptr createSolveTaskHarm(UInt freqIdx) override; + /// Create a task for eigenvalue extraction + std::shared_ptr createExtractEigenvaluesTask() override; /// Logging of system matrices and source vector void logSystemMatrices() override; /// Solves system for single frequency @@ -163,6 +165,12 @@ template class MnaSolverDirect : public MnaSolver { /// Destructor virtual ~MnaSolverDirect() = default; + // ### Eigenvalue extraction ### + /// Extract eigenvalues from power system conductance matrix + void extractEigenvalues(Real time, Int timeStepCount) override; + /// + void closeEigenvalueLogger() override; + /// Sets the linear solver to "implementation" and creates an object void setDirectLinearSolverImplementation(DirectLinearSolverImpl implementation); @@ -258,6 +266,23 @@ template class MnaSolverDirect : public MnaSolver { MnaSolverDirect &mSolver; }; + /// + class ExtractEigenvaluesTask : public CPS::Task { + public: + ExtractEigenvaluesTask(MnaSolverDirect &solver) + : Task(solver.mName + ".ExtractEigenvalues"), mSolver(solver) { + mAttributeDependencies.push_back(solver.mLeftSideVector); + mModifiedAttributes.push_back(Scheduler::external); + } + + void execute(Real time, Int timeStepCount) { + mSolver.extractEigenvalues(time, timeStepCount); + } + + private: + MnaSolverDirect &mSolver; + }; + /// class LogTask : public CPS::Task { public: diff --git a/dpsim/include/dpsim/Simulation.h b/dpsim/include/dpsim/Simulation.h index 4d58db8dd0..364058f9c9 100644 --- a/dpsim/include/dpsim/Simulation.h +++ b/dpsim/include/dpsim/Simulation.h @@ -95,6 +95,8 @@ class Simulation : public CPS::AttributeList { Bool mInitFromNodesAndTerminals = true; /// Enable recomputation of system matrix during simulation Bool mSystemMatrixRecomputation = false; + /// Enable eigenvalue extraction during simulation + Bool mIsEigenvalueExtractionEnabled = false; /// If tearing components exist, the Diakoptics /// solver is selected automatically. @@ -174,6 +176,11 @@ class Simulation : public CPS::AttributeList { void setSolverType(Solver::Type solverType = Solver::Type::MNA) { mSolverType = solverType; } + + /// Enable eigenvalue extraction during simulation. + void doEigenvalueExtraction(Bool isEigenvalueExtractionEnabled = true) { + mIsEigenvalueExtractionEnabled = isEigenvalueExtractionEnabled; + } /// set solver and component to initialization or simulation behaviour void setSolverAndComponentBehaviour(Solver::Behaviour behaviour) { mSolverBehaviour = behaviour; diff --git a/dpsim/include/dpsim/Solver.h b/dpsim/include/dpsim/Solver.h index 72f12ddc68..4b5c2e2ad5 100644 --- a/dpsim/include/dpsim/Solver.h +++ b/dpsim/include/dpsim/Solver.h @@ -64,6 +64,9 @@ class Solver { /// Solver behaviour initialization or simulation Behaviour mBehaviour = Solver::Behaviour::Simulation; + /// Enables extraction of eigenvalues + Bool mIsEigenvalueExtractionEnabled = false; + public: Solver(String name, CPS::Logger::Level logLevel) : mName(name), mLogLevel(logLevel), @@ -123,5 +126,18 @@ class Solver { void setMaxNumberOfIterations(int maxIterations) { mMaxIterations = maxIterations; } + + /// ### Eigenvalue Extraction ### + void doEigenvalueExtraction(Bool isEigenvalueExtractionEnabled) { + mIsEigenvalueExtractionEnabled = isEigenvalueExtractionEnabled; + } + /// + virtual void extractEigenvalues(Real time, Int timeStepCount){ + // no default implementation for all types of solvers + }; + /// + virtual void closeEigenvalueLogger(){ + // no default implementation for all types of solvers + }; }; } // namespace DPsim diff --git a/dpsim/src/CMakeLists.txt b/dpsim/src/CMakeLists.txt index 0f0460a280..7ca7a0a91d 100644 --- a/dpsim/src/CMakeLists.txt +++ b/dpsim/src/CMakeLists.txt @@ -3,6 +3,9 @@ set(DPSIM_SOURCES RealTimeSimulation.cpp MNASolver.cpp MNASolverDirect.cpp + MNAEigenvalueExtractor.cpp + MNAEigenvalueLogger.cpp + EigenvalueUtils.cpp DenseLUAdapter.cpp SparseLUAdapter.cpp DirectLinearSolverConfiguration.cpp diff --git a/dpsim/src/EigenvalueUtils.cpp b/dpsim/src/EigenvalueUtils.cpp new file mode 100644 index 0000000000..cfe4e274fb --- /dev/null +++ b/dpsim/src/EigenvalueUtils.cpp @@ -0,0 +1,45 @@ +#include + +using namespace DPsim; + +// TODO: consider resizing existing matrix instead of creating a new one (performance, memory) +MatrixComp EigenvalueUtils::returnNonZeroElements(const MatrixComp &mat) { + Eigen::Matrix mask = + (mat.array().abs() > 1e-14); + int nonZeroCount = mask.cast().sum(); + + Eigen::MatrixXcd nonZeroMatrix(nonZeroCount, 1); + int index = 0; + for (int i = 0; i < mask.rows(); ++i) { + for (int j = 0; j < mask.cols(); ++j) { + if (mask(i, j)) { + nonZeroMatrix(index++) = mat(i, j); + } + } + } + return nonZeroMatrix; +} + +MatrixComp EigenvalueUtils::convertRealEquivalentToComplexMatrix( + const Matrix &realEquivalentMatrix) { + // The size of the complex matrix is half the size of the real matrix + int size = realEquivalentMatrix.rows() / 2; + + // Create a complex matrix of the appropriate size + MatrixComp complexMatrix(size, size); + + // Iterate over the complex matrix + for (int i = 0; i < size; ++i) { + for (int j = 0; j < size; ++j) { + // The real part is in the upper left quadrant of the real matrix + double realPart = realEquivalentMatrix(i, j); + + // The imaginary part is in the lower left quadrant of the real matrix + double imagPart = realEquivalentMatrix(i + size, j); + + // Assign the complex number to the complex matrix + complexMatrix(i, j) = std::complex(realPart, imagPart); + } + } + return complexMatrix; +} diff --git a/dpsim/src/MNAEigenvalueExtractor.cpp b/dpsim/src/MNAEigenvalueExtractor.cpp new file mode 100644 index 0000000000..4561482a6e --- /dev/null +++ b/dpsim/src/MNAEigenvalueExtractor.cpp @@ -0,0 +1,147 @@ +#include + +namespace DPsim { +template +MNAEigenvalueExtractor::MNAEigenvalueExtractor( + CPS::Logger::Level logLevel) + : mLogger("MNAEigenvalueExtractor", logLevel) { + mEigenvalues = CPS::AttributeStatic::make(); + mDiscreteEigenvalues = CPS::AttributeStatic::make(); +} + +template +void MNAEigenvalueExtractor::initialize( + const CPS::SystemTopology &topology, UInt numMatrixNodeIndices, + Real timeStep) { + setParameters(topology, timeStep); + identifyEigenvalueComponents(topology.mComponents); + createEmptyEigenvalueMatrices(numMatrixNodeIndices); + stampEigenvalueMatrices(); + mLogger.setLogAttributes(mEigenvalues, mDiscreteEigenvalues); + mLogger.logInitialization(mSignMatrix, mDiscretizationMatrix, + mBranchNodeIncidenceMatrix, + mNodeBranchIncidenceMatrix); +} + +template +void MNAEigenvalueExtractor::setParameters( + const CPS::SystemTopology &topology, Real timeStep) { + mTimeStep = timeStep; + // Relevant only for DP + mSystemOmega = topology.mSystemOmega; + + Real denominator = mTimeStep * mTimeStep * mSystemOmega * mSystemOmega + 4.0; + Real realPart = + (4.0 - mTimeStep * mTimeStep * mSystemOmega * mSystemOmega) / denominator; + Real imagPart = (-4.0 * mTimeStep * mSystemOmega) / denominator; + mCoeffDP = Complex(realPart, imagPart); +} + +template +void MNAEigenvalueExtractor::identifyEigenvalueComponents( + const CPS::IdentifiedObject::List &components) { + // TODO: throw exception if topology contains components that do not implement EigenvalueCompInterface + UInt branchIdx = 0; + for (auto component : components) { + auto eigenvalueComponent = + std::dynamic_pointer_cast(component); + if (eigenvalueComponent) { + mEigenvalueComponentToBranchIdx[eigenvalueComponent] = branchIdx; + auto eigenvalueDynamicComponent = std::dynamic_pointer_cast< + CPS::EigenvalueDynamicCompInterface>(component); + if (eigenvalueDynamicComponent) { + mEigenvalueDynamicComponents.push_back(eigenvalueDynamicComponent); + } + branchIdx++; + } + } +} + +template +void MNAEigenvalueExtractor::createEmptyEigenvalueMatrices( + UInt numMatrixNodeIndices) { + int nBranches = mEigenvalueComponentToBranchIdx.size(); + mSignMatrix = MatrixVar::Zero(nBranches, nBranches); + mDiscretizationMatrix = MatrixVar::Zero(nBranches, nBranches); + mBranchNodeIncidenceMatrix = Matrix::Zero(nBranches, numMatrixNodeIndices); + **mEigenvalues = MatrixComp::Zero(mEigenvalueDynamicComponents.size(), 1); + **mDiscreteEigenvalues = + MatrixComp::Zero(mEigenvalueDynamicComponents.size(), 1); +} + +template +void MNAEigenvalueExtractor::stampEigenvalueMatrices() { + for (const auto &compToBranchIdx : mEigenvalueComponentToBranchIdx) { + auto comp = compToBranchIdx.first; + UInt branchIdx = compToBranchIdx.second; + comp->stampBranchNodeIncidenceMatrix(branchIdx, mBranchNodeIncidenceMatrix); + } + mNodeBranchIncidenceMatrix = mBranchNodeIncidenceMatrix.transpose(); + for (const auto &dynamicComp : mEigenvalueDynamicComponents) { + UInt branchIdx = mEigenvalueComponentToBranchIdx[dynamicComp]; + dynamicComp->stampSignMatrix(branchIdx, mSignMatrix, mCoeffDP); + dynamicComp->stampDiscretizationMatrix(branchIdx, mDiscretizationMatrix, + mCoeffDP); + } +} + +template +void MNAEigenvalueExtractor::extractEigenvalues( + const Matrix &powerSystemMatrix, Real time, Int timeStepCount) { + calculateStateMatrix(powerSystemMatrix); + computeDiscreteEigenvalues(); + recoverEigenvalues(); + mLogger.logExtraction(time, timeStepCount, mStateMatrix); +} + +template <> +void MNAEigenvalueExtractor::calculateStateMatrix( + const Matrix &powerSystemMatrix) { + // TODO: use right hand side solving of factorized power system matrix instead of inversion (performance). + Matrix intermediateResult = + powerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix; + mStateMatrix = mSignMatrix - mDiscretizationMatrix * + mBranchNodeIncidenceMatrix * + intermediateResult; +} + +template <> +void MNAEigenvalueExtractor::calculateStateMatrix( + const Matrix &powerSystemMatrix) { + // TODO: use right hand side solving of factorized power system matrix instead of inversion (performance). + MatrixComp compPowerSystemMatrix = + EigenvalueUtils::convertRealEquivalentToComplexMatrix(powerSystemMatrix); + MatrixComp intermediateResult = + compPowerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix; + mStateMatrix = mSignMatrix - mDiscretizationMatrix * + mBranchNodeIncidenceMatrix * + intermediateResult; +} + +template +void MNAEigenvalueExtractor::computeDiscreteEigenvalues() { + auto discreteEigenvaluesIncludingZeros = mStateMatrix.eigenvalues(); + **mDiscreteEigenvalues = + EigenvalueUtils::returnNonZeroElements(discreteEigenvaluesIncludingZeros); + // TODO: filter out eigenvalues = -1 + 0i to avoid division by zero in recoverEigenvalues(). +} + +template <> void MNAEigenvalueExtractor::recoverEigenvalues() { + **mEigenvalues = 2.0 / mTimeStep * ((**mDiscreteEigenvalues).array() - 1.0) / + ((**mDiscreteEigenvalues).array() + 1.0); +} + +template <> void MNAEigenvalueExtractor::recoverEigenvalues() { + **mEigenvalues = 2.0 / mTimeStep * ((**mDiscreteEigenvalues).array() - 1.0) / + ((**mDiscreteEigenvalues).array() + 1.0) + + Complex(0.0, 1.0) * mSystemOmega; +} + +template +void MNAEigenvalueExtractor::closeLogger() { + mLogger.close(); +} + +template class MNAEigenvalueExtractor; +template class MNAEigenvalueExtractor; +} // namespace DPsim diff --git a/dpsim/src/MNAEigenvalueLogger.cpp b/dpsim/src/MNAEigenvalueLogger.cpp new file mode 100644 index 0000000000..4f32003dd5 --- /dev/null +++ b/dpsim/src/MNAEigenvalueLogger.cpp @@ -0,0 +1,67 @@ +#include + +namespace DPsim { +MNAEigenvalueLogger::MNAEigenvalueLogger(String name, + CPS::Logger::Level logLevel) + : mSLog(CPS::Logger::get(name, logLevel, logLevel)), + mEigenvaluesLogger("eigenvalues", true, 1), + mDiscreteEigenvaluesLogger("discreteEigenvalues", true, 1) {} + +void MNAEigenvalueLogger::setLogAttributes( + CPS::AttributeStatic::Ptr eigenvalues, + CPS::AttributeStatic::Ptr discreteEigenvalues) { + mEigenvaluesLogger.logAttribute("eigenvalues", eigenvalues); + mDiscreteEigenvaluesLogger.logAttribute("discreteEigenvalues", + discreteEigenvalues); +} + +template +void MNAEigenvalueLogger::logInitialization( + const MatrixVar &signMatrix, + const MatrixVar &discretizationMatrix, + const Matrix &branchNodeIncidenceMatrix, + const Matrix &nodeBranchIncidenceMatrix) { + SPDLOG_LOGGER_DEBUG(mSLog, "---- Initialize ----"); + SPDLOG_LOGGER_DEBUG(mSLog, "sign matrix: {}", + CPS::Logger::matrixVarToString(signMatrix)); + SPDLOG_LOGGER_DEBUG(mSLog, "discretization matrix: {}", + CPS::Logger::matrixVarToString(discretizationMatrix)); + SPDLOG_LOGGER_DEBUG(mSLog, "branch <-> node incidence matrix: {}", + CPS::Logger::matrixToString(branchNodeIncidenceMatrix)); + SPDLOG_LOGGER_DEBUG(mSLog, "node <-> branch incidence matrix: {}", + CPS::Logger::matrixToString(nodeBranchIncidenceMatrix)); +} +template void MNAEigenvalueLogger::logInitialization( + const MatrixVar &signMatrix, + const MatrixVar &discretizationMatrix, + const Matrix &branchNodeIncidenceMatrix, + const Matrix &nodeBranchIncidenceMatrix); +template void MNAEigenvalueLogger::logInitialization( + const MatrixVar &signMatrix, + const MatrixVar &discretizationMatrix, + const Matrix &branchNodeIncidenceMatrix, + const Matrix &nodeBranchIncidenceMatrix); + +template +void MNAEigenvalueLogger::logExtraction(Real time, Int timeStepCount, + const MatrixVar &stateMatrix) { + SPDLOG_LOGGER_DEBUG(mSLog, "---- Extract eigenvalues ----"); + SPDLOG_LOGGER_DEBUG(mSLog, "time: {}", CPS::Logger::realToString(time)); + SPDLOG_LOGGER_DEBUG(mSLog, "discretized state matrix: {}", + CPS::Logger::matrixVarToString(stateMatrix)); + + mEigenvaluesLogger.log(time, timeStepCount); + mDiscreteEigenvaluesLogger.log(time, timeStepCount); +} +template void +MNAEigenvalueLogger::logExtraction(Real time, Int timeStepCount, + const MatrixVar &stateMatrix); +template void MNAEigenvalueLogger::logExtraction( + Real time, Int timeStepCount, const MatrixVar &stateMatrix); + +void MNAEigenvalueLogger::close() { + mEigenvaluesLogger.close(); + mDiscreteEigenvaluesLogger.close(); + mSLog->flush(); +} +} // namespace DPsim diff --git a/dpsim/src/MNASolver.cpp b/dpsim/src/MNASolver.cpp index 3bb1291fd0..629dbfdfe4 100644 --- a/dpsim/src/MNASolver.cpp +++ b/dpsim/src/MNASolver.cpp @@ -18,7 +18,8 @@ namespace DPsim { template MnaSolver::MnaSolver(String name, CPS::Domain domain, CPS::Logger::Level logLevel) - : Solver(name, logLevel), mDomain(domain) { + : Solver(name, logLevel), mDomain(domain), + mMNAEigenvalueExtractor(logLevel) { // Raw source and solution vector logging mLeftVectorLog = std::make_shared( @@ -96,6 +97,12 @@ template void MnaSolver::initialize() { SPDLOG_LOGGER_INFO(mSLog, "--- Initial system matrices and vectors ---"); logSystemMatrices(); + if (Solver::mIsEigenvalueExtractionEnabled) { + SPDLOG_LOGGER_INFO(mSLog, "--- Initialize eigenvalue extractor ---"); + mMNAEigenvalueExtractor.initialize(mSystem, mNumMatrixNodeIndices, + Solver::mTimeStep); + } + mSLog->flush(); } @@ -615,6 +622,9 @@ template Task::List MnaSolver::getTasks() { l.push_back(createSolveTaskRecomp()); } else { l.push_back(createSolveTask()); + if (mIsEigenvalueExtractionEnabled) { + l.push_back(createExtractEigenvaluesTask()); + } l.push_back(createLogTask()); } return l; diff --git a/dpsim/src/MNASolverDirect.cpp b/dpsim/src/MNASolverDirect.cpp index 58404822fc..82e12cdd57 100644 --- a/dpsim/src/MNASolverDirect.cpp +++ b/dpsim/src/MNASolverDirect.cpp @@ -225,6 +225,13 @@ MnaSolverDirect::createSolveTaskHarm(UInt freqIdx) { freqIdx); } +template +std::shared_ptr +MnaSolverDirect::createExtractEigenvaluesTask() { + return std::make_shared::ExtractEigenvaluesTask>( + *this); +} + template std::shared_ptr MnaSolverDirect::createLogTask() { return std::make_shared::LogTask>(*this); @@ -325,6 +332,19 @@ void MnaSolverDirect::solveWithHarmonics(Real time, Int timeStepCount, mRightSideVectorHarm[freqIdx]); } +template +void MnaSolverDirect::extractEigenvalues(Real time, + Int timeStepCount) { + MnaSolver::mMNAEigenvalueExtractor.extractEigenvalues( + ((Matrix)mSwitchedMatrices[mCurrentSwitchStatus][0]), time, + timeStepCount); +} + +template +void MnaSolverDirect::closeEigenvalueLogger() { + MnaSolver::mMNAEigenvalueExtractor.closeLogger(); +} + template void MnaSolverDirect::logSystemMatrices() { if (mFrequencyParallel) { for (UInt i = 0; i < mSwitchedMatrices[std::bitset(0)].size(); diff --git a/dpsim/src/Simulation.cpp b/dpsim/src/Simulation.cpp index 1d79b4f8f4..662803c504 100644 --- a/dpsim/src/Simulation.cpp +++ b/dpsim/src/Simulation.cpp @@ -163,6 +163,7 @@ template void Simulation::createMNASolver() { solver->doSystemMatrixRecomputation(mSystemMatrixRecomputation); solver->setDirectLinearSolverConfiguration( mDirectLinearSolverConfiguration); + solver->doEigenvalueExtraction(mIsEigenvalueExtractionEnabled); solver->initialize(); solver->setMaxNumberOfIterations(mMaxIterations); } @@ -379,6 +380,9 @@ void Simulation::stop() { for (auto lg : mLoggers) lg->close(); + for (auto solver : mSolvers) + solver->closeEigenvalueLogger(); + SPDLOG_LOGGER_INFO(mLog, "Simulation finished."); mLog->flush(); } diff --git a/dpsim/src/pybind/main.cpp b/dpsim/src/pybind/main.cpp index a1392119e1..7e5d87e3fd 100644 --- a/dpsim/src/pybind/main.cpp +++ b/dpsim/src/pybind/main.cpp @@ -208,7 +208,9 @@ PYBIND11_MODULE(dpsimpy, m) { &DPsim::Simulation::setDirectLinearSolverImplementation) .def("set_direct_linear_solver_configuration", &DPsim::Simulation::setDirectLinearSolverConfiguration) - .def("log_lu_times", &DPsim::Simulation::logLUTimes); + .def("log_lu_times", &DPsim::Simulation::logLUTimes) + .def("do_eigenvalue_extraction", + &DPsim::Simulation::doEigenvalueExtraction); py::class_(m, "RealTimeSimulation") diff --git a/examples/Notebooks/Eigenvalues/DP_EMT_SinglePhaseRLC_Tests.ipynb b/examples/Notebooks/Eigenvalues/DP_EMT_SinglePhaseRLC_Tests.ipynb new file mode 100644 index 0000000000..ae25ddba23 --- /dev/null +++ b/examples/Notebooks/Eigenvalues/DP_EMT_SinglePhaseRLC_Tests.ipynb @@ -0,0 +1,226 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Eigenvalue extraction tests in a single-phase RLC circuit in DP & EMT domains" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import dpsimpy\n", + "import math\n", + "import villas.dataprocessing.readtools as rt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# function to assert that two complex numbers are close\n", + "def assert_complex_isclose(expected, actual, rel_tol=1e-6):\n", + " assert math.isclose(expected.real, actual.real, rel_tol=rel_tol), \"Real parts not close: {} vs {}\".format(expected.real, actual.real)\n", + " assert math.isclose(expected.imag, actual.imag, rel_tol=rel_tol), \"Imaginary parts not close: {} vs {}\".format(expected.imag, actual.imag)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Circuit parameters\n", + "v_volt = 8.5\n", + "r_ohm = 100\n", + "l_henry = 5\n", + "c_farad = 250e-6\n", + "deltaT = 1e-4\n", + "final_time = 1e-3" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Analytical solution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eigenvalue0_expected = ((c_farad**2*r_ohm**2 - 4*l_henry*c_farad)**(1/2) - c_farad*r_ohm)/(2*c_farad*l_henry)\n", + "eigenvalue1_expected = -((c_farad**2*r_ohm**2 - 4*l_henry*c_farad)**(1/2) + c_farad*r_ohm)/(2*c_farad*l_henry)\n", + "\n", + "print('Expected eigenvalue 0: ' + str(eigenvalue0_expected))\n", + "print('Expected eigenvalue 1: ' + str(eigenvalue1_expected))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Test 1: Eigenvalues extracted from simulation in DP domain match analytical solution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# DPsim DP simulation\n", + "name = 'DP_SinglePhaseRLC_Test'\n", + "log_dir = \"logs/\" + name\n", + "dpsimpy.Logger.set_log_dir(log_dir)\n", + "\n", + "# Create nodes\n", + "gnd = dpsimpy.dp.SimNode.gnd\n", + "n1 = dpsimpy.dp.SimNode('n1')\n", + "n2 = dpsimpy.dp.SimNode('n2')\n", + "n3 = dpsimpy.dp.SimNode('n3')\n", + "\n", + "# Create components\n", + "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", + "vs.set_parameters(V_ref=complex(v_volt,0))\n", + "r = dpsimpy.dp.ph1.Resistor('r')\n", + "r.set_parameters(R=r_ohm)\n", + "l = dpsimpy.dp.ph1.Inductor('l')\n", + "l.set_parameters(L=l_henry)\n", + "c = dpsimpy.dp.ph1.Capacitor('c')\n", + "c.set_parameters(C=c_farad)\n", + "\n", + "# Set connections\n", + "vs.connect([gnd, n1])\n", + "l.connect([n1, n2])\n", + "c.connect([n2, n3])\n", + "r.connect([n3, gnd])\n", + "\n", + "# Create topology\n", + "system = dpsimpy.SystemTopology(50, [n1, n2, n3], [vs, r, l, c])\n", + "\n", + "# Configure and run simulation\n", + "sim = dpsimpy.Simulation(name)\n", + "sim.set_domain(dpsimpy.Domain.DP)\n", + "sim.set_system(system)\n", + "sim.do_eigenvalue_extraction(True)\n", + "sim.set_time_step(deltaT)\n", + "sim.set_final_time(final_time)\n", + "sim.run()\n", + "\n", + "# Read log file\n", + "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n", + "eigenvalue0 = eigenvalues_timeseries['eigenvalues_0'].values\n", + "eigenvalue1 = eigenvalues_timeseries['eigenvalues_1'].values\n", + "\n", + "# Assert\n", + "nSamples = int(final_time/deltaT)\n", + "\n", + "assert len(eigenvalue0) == nSamples\n", + "assert len(eigenvalue1) == nSamples\n", + "\n", + "for i in range(nSamples):\n", + " assert_complex_isclose(eigenvalue0_expected, eigenvalue0[i], 1e-8)\n", + " assert_complex_isclose(eigenvalue1_expected, eigenvalue1[i], 1e-8)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Test 2: Eigenvalues extracted from simulation in EMT domain match analytical solution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# DPsim EMT simulation\n", + "name = 'EMT_SinglePhaseRLC_Test'\n", + "log_dir = \"logs/\" + name\n", + "dpsimpy.Logger.set_log_dir(log_dir)\n", + "\n", + "# Create nodes\n", + "gnd = dpsimpy.emt.SimNode.gnd\n", + "n1 = dpsimpy.emt.SimNode('n1')\n", + "n2 = dpsimpy.emt.SimNode('n2')\n", + "n3 = dpsimpy.emt.SimNode('n3')\n", + "\n", + "# Create components\n", + "vs = dpsimpy.emt.ph1.VoltageSource('vs')\n", + "vs.set_parameters(V_ref=complex(v_volt,0))\n", + "r = dpsimpy.emt.ph1.Resistor('r')\n", + "r.set_parameters(R=r_ohm)\n", + "l = dpsimpy.emt.ph1.Inductor('l')\n", + "l.set_parameters(L=l_henry)\n", + "c = dpsimpy.emt.ph1.Capacitor('c')\n", + "c.set_parameters(C=c_farad)\n", + "\n", + "# Set connections\n", + "vs.connect([gnd, n1])\n", + "l.connect([n1, n2])\n", + "c.connect([n2, n3])\n", + "r.connect([n3, gnd])\n", + "\n", + "# Create topology\n", + "system = dpsimpy.SystemTopology(50, [n1, n2, n3], [vs, r, l, c])\n", + "\n", + "# Configure and run simulation\n", + "sim = dpsimpy.Simulation(name)\n", + "sim.set_domain(dpsimpy.Domain.EMT)\n", + "sim.set_system(system)\n", + "sim.do_eigenvalue_extraction(True)\n", + "sim.set_time_step(deltaT)\n", + "sim.set_final_time(final_time)\n", + "sim.run()\n", + "\n", + "# Read log file\n", + "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n", + "eigenvalue0 = eigenvalues_timeseries['eigenvalues_0'].values\n", + "eigenvalue1 = eigenvalues_timeseries['eigenvalues_1'].values\n", + "\n", + "# Assert\n", + "nSamples = int(final_time/deltaT)\n", + "\n", + "assert len(eigenvalue0) == nSamples\n", + "assert len(eigenvalue1) == nSamples\n", + "\n", + "for i in range(nSamples):\n", + " assert_complex_isclose(eigenvalue0_expected, eigenvalue0[i], 1e-8)\n", + " assert_complex_isclose(eigenvalue1_expected, eigenvalue1[i], 1e-8)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/Notebooks/Eigenvalues/DP_SinglePhaseRLC_Demo.ipynb b/examples/Notebooks/Eigenvalues/DP_SinglePhaseRLC_Demo.ipynb new file mode 100644 index 0000000000..909d0c8df2 --- /dev/null +++ b/examples/Notebooks/Eigenvalues/DP_SinglePhaseRLC_Demo.ipynb @@ -0,0 +1,164 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Eigenvalue extraction during simulation of a single-phase RLC circuit.\n", + "\n", + "- Simulate RLC circuit in DP domain.\n", + "- Simulation includes switch events changing resistance value of the circuit.\n", + "- Extract eigenvalues at each time step of simulation.\n", + "- Visualize eigenvalue trajectory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import dpsimpy\n", + "import villas.dataprocessing.readtools as rt\n", + "import villas.dataprocessing.plottools as pt\n", + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# DPsim DP simulation\n", + "name = 'DP_SinglePhaseRLC_Demo'\n", + "log_dir = \"logs/\" + name\n", + "dpsimpy.Logger.set_log_dir(log_dir)\n", + "\n", + "# Create nodes\n", + "gnd = dpsimpy.dp.SimNode.gnd\n", + "n1 = dpsimpy.dp.SimNode('n1')\n", + "n2 = dpsimpy.dp.SimNode('n2')\n", + "n3 = dpsimpy.dp.SimNode('n3')\n", + "n4 = dpsimpy.dp.SimNode('n4')\n", + "n5 = dpsimpy.dp.SimNode('n5')\n", + "\n", + "# Create components\n", + "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", + "vs.set_parameters(V_ref=complex(100,0))\n", + "r1 = dpsimpy.dp.ph1.Resistor('r_1')\n", + "r1.set_parameters(R=10)\n", + "r2 = dpsimpy.dp.ph1.Resistor('r_2')\n", + "r2.set_parameters(R=0.1)\n", + "l = dpsimpy.dp.ph1.Inductor('l')\n", + "l.set_parameters(L=5)\n", + "c = dpsimpy.dp.ph1.Capacitor('c')\n", + "c.set_parameters(C=250e-6)\n", + "s1 = dpsimpy.dp.ph1.Switch('s1')\n", + "s1.set_parameters(1e8, 1e-4, True)\n", + "s2 = dpsimpy.dp.ph1.Switch('s2')\n", + "s2.set_parameters(1e8, 1e-4, False)\n", + "\n", + "# Create switch events\n", + "s1Event = dpsimpy.event.SwitchEvent(0.05, s1, False)\n", + "s2Event = dpsimpy.event.SwitchEvent(0.05, s2, True)\n", + "\n", + "# Set connections\n", + "vs.connect([gnd, n1])\n", + "r1.connect([n4, gnd])\n", + "r2.connect([n5, gnd])\n", + "l.connect([n1, n2])\n", + "c.connect([n2, n3])\n", + "s1.connect([n3, n4])\n", + "s2.connect([n3, n5])\n", + "\n", + "# Create topology\n", + "system = dpsimpy.SystemTopology(50, [n1, n2, n3, n4, n5], [vs, r1, r2, l, c, s1, s2])\n", + "\n", + "# Configure logging\n", + "logger = dpsimpy.Logger(name)\n", + "logger.log_attribute('current', 'i_intf', l)\n", + "logger.log_attribute('voltage', 'v', n3)\n", + "\n", + "# Configure and run simulation\n", + "sim = dpsimpy.Simulation(name)\n", + "sim.set_domain(dpsimpy.Domain.DP)\n", + "sim.set_system(system)\n", + "sim.add_event(s1Event)\n", + "sim.add_event(s2Event)\n", + "sim.do_eigenvalue_extraction(True)\n", + "sim.set_time_step(0.0001)\n", + "sim.set_final_time(0.1)\n", + "sim.add_logger(logger)\n", + "sim.run()\n", + "\n", + "# Read log file\n", + "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n", + "time = eigenvalues_timeseries['eigenvalues_0'].time\n", + "eigenvalue1_real = eigenvalues_timeseries['eigenvalues_0'].values.real\n", + "eigenvalue1_imag = eigenvalues_timeseries['eigenvalues_0'].values.imag\n", + "eigenvalue2_real = eigenvalues_timeseries['eigenvalues_1'].values.real\n", + "eigenvalue2_imag = eigenvalues_timeseries['eigenvalues_1'].values.imag\n", + "\n", + "simulation_timeseries = rt.read_timeseries_dpsim(log_dir + '/' + name + '.csv')\n", + "simulation_timeseries_time_domain = ts.frequency_shift_list(simulation_timeseries, 50)\n", + "print(simulation_timeseries_time_domain)\n", + "current_timeseries = simulation_timeseries_time_domain['current_shift']\n", + "voltage_timeseries = simulation_timeseries_time_domain['voltage_shift']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.title('Current')\n", + "plt.plot(current_timeseries.time, current_timeseries.values)\n", + "plt.show()\n", + "\n", + "plt.figure()\n", + "plt.title('Voltage across resistor')\n", + "plt.plot(voltage_timeseries.time, voltage_timeseries.values)\n", + "plt.show()\n", + "\n", + "plt.figure()\n", + "plt.title('First eigenvalue')\n", + "plt.plot(time, eigenvalue1_real, label='Real part')\n", + "plt.plot(time, eigenvalue1_imag, label='Imaginary part')\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "plt.figure()\n", + "plt.title('Second eigenvalue')\n", + "plt.plot(time, eigenvalue2_real, label='Real part')\n", + "plt.plot(time, eigenvalue2_imag, label='Imaginary part')\n", + "plt.legend()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}