From ffe67b68dc8d79ebc9918e30b4e1161812266c90 Mon Sep 17 00:00:00 2001 From: Georgii Tishenin Date: Tue, 16 Jan 2024 12:06:04 +0100 Subject: [PATCH] Make MNAEigenvalueExtractor a template class with specialization for Matrix and MatrixComp types. Matrix type is used for EMT, MatrixComp type is used for DP domain. Signed-off-by: Georgii Tishenin --- dpsim/include/dpsim/MNAEigenvalueExtractor.h | 24 ++-- dpsim/include/dpsim/MNASolver.h | 7 +- dpsim/src/MNAEigenvalueExtractor.cpp | 140 +++++++++++++++---- dpsim/src/MNASolver.cpp | 2 +- dpsim/src/MNASolverDirect.cpp | 2 +- dpsim/src/Simulation.cpp | 24 +++- 6 files changed, 156 insertions(+), 43 deletions(-) diff --git a/dpsim/include/dpsim/MNAEigenvalueExtractor.h b/dpsim/include/dpsim/MNAEigenvalueExtractor.h index 5d0af92b05..1e797fc814 100644 --- a/dpsim/include/dpsim/MNAEigenvalueExtractor.h +++ b/dpsim/include/dpsim/MNAEigenvalueExtractor.h @@ -3,37 +3,45 @@ #include #include #include -#include +#include namespace DPsim { /// Extracts eigenvalues using MNA power system conductance matrix + template class MNAEigenvalueExtractor { public: MNAEigenvalueExtractor(); - MNAEigenvalueExtractor(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices); + MNAEigenvalueExtractor(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices, Real timeStep); - void initialize(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices); - void extractEigenvalues(const Matrix &powerSystemMatrix, Real timeStep); + void initialize(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices, Real timeStep); + void extractEigenvalues(const Matrix &powerSystemMatrix); private: CPS::EigenvalueCompInterface::List mEigenvalueComponents; - Matrix mSignMatrix; - Matrix mDiscretizationMatrix; + typename CPS::EigenvalueDynamicCompInterface::List mEigenvalueDynamicComponents; + Real mTimeStep; + Real mSystemOmega; + Complex mCoeffDP; + MatrixType mSignMatrix; + MatrixType mDiscretizationMatrix; Matrix mBranchNodeIncidenceMatrix; Matrix mNodeBranchIncidenceMatrix; - Matrix mStateMatrix; + MatrixType mStateMatrix; MatrixComp mDiscreteEigenvalues; MatrixComp mEigenvalues; CPS::Logger::Log mSLog; + void setParameters(const CPS::SystemTopology &topology, Real timeStep); void identifyEigenvalueComponents(const CPS::IdentifiedObject::List &components); void setBranchIndices(); void createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices); void stampEigenvalueMatrices(); + void logInitialization(); void calculateStateMatrix(const Matrix &powerSystemMatrix); void computeDiscreteEigenvalues(); - void recoverEigenvalues(Real timeStep); + void recoverEigenvalues(); + void logExtraction(); }; } diff --git a/dpsim/include/dpsim/MNASolver.h b/dpsim/include/dpsim/MNASolver.h index 742a897243..9c7136e40b 100644 --- a/dpsim/include/dpsim/MNASolver.h +++ b/dpsim/include/dpsim/MNASolver.h @@ -37,7 +37,7 @@ namespace DPsim { /// Solver class using Modified Nodal Analysis (MNA). template - class MnaSolver : public Solver { + class MnaSolver : public Solver { protected: // #### General simulation settings #### /// Simulation domain, which can be dynamic phasor (DP) or EMT @@ -123,7 +123,10 @@ namespace DPsim { std::vector mRecomputationTimes; // #### Eigenvalue extraction #### - MNAEigenvalueExtractor mMNAEigenvalueExtractor; + // Define MatrixType as an alias for Matrix if VarType is Real, otherwise define it as an alias for MatrixComp + using MatrixType = typename std::conditional::value, Matrix, MatrixComp>::type; + MNAEigenvalueExtractor mMNAEigenvalueExtractor; + /// Constructor should not be called by users but by Simulation MnaSolver(String name, CPS::Domain domain = CPS::Domain::DP, diff --git a/dpsim/src/MNAEigenvalueExtractor.cpp b/dpsim/src/MNAEigenvalueExtractor.cpp index 7fd3466a5c..eefc2beb60 100644 --- a/dpsim/src/MNAEigenvalueExtractor.cpp +++ b/dpsim/src/MNAEigenvalueExtractor.cpp @@ -2,30 +2,42 @@ namespace DPsim { - MNAEigenvalueExtractor::MNAEigenvalueExtractor() : mSLog(CPS::Logger::get("MNAEigenvalueExtractor", CPS::Logger::Level::info, CPS::Logger::Level::info)) {} + template + MNAEigenvalueExtractor::MNAEigenvalueExtractor() : mSLog(CPS::Logger::get("MNAEigenvalueExtractor", CPS::Logger::Level::info, CPS::Logger::Level::info)) {} - MNAEigenvalueExtractor::MNAEigenvalueExtractor(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices) + template + MNAEigenvalueExtractor::MNAEigenvalueExtractor(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices, Real timeStep) : mSLog(CPS::Logger::get("MNAEigenvalueExtractor", CPS::Logger::Level::info, CPS::Logger::Level::info)) { - initialize(topology, numMatrixNodeIndices); + initialize(topology, numMatrixNodeIndices, timeStep); } - void MNAEigenvalueExtractor::initialize(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices) + template + void MNAEigenvalueExtractor::initialize(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices, Real timeStep) { + setParameters(topology, timeStep); identifyEigenvalueComponents(topology.mComponents); setBranchIndices(); createEmptyEigenvalueMatrices(numMatrixNodeIndices); stampEigenvalueMatrices(); + logInitialization(); + } - SPDLOG_LOGGER_INFO(mSLog, "---- Initialize ----"); - SPDLOG_LOGGER_INFO(mSLog, "sign matrix: {}", CPS::Logger::matrixToString(mSignMatrix)); - SPDLOG_LOGGER_INFO(mSLog, "discretization matrix: {}", CPS::Logger::matrixToString(mDiscretizationMatrix)); - SPDLOG_LOGGER_INFO(mSLog, "branch <-> node incidence matrix: {}", CPS::Logger::matrixToString(mBranchNodeIncidenceMatrix)); - SPDLOG_LOGGER_INFO(mSLog, "node <-> branch incidence matrix: {}", CPS::Logger::matrixToString(mNodeBranchIncidenceMatrix)); - mSLog->flush(); + 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); } - void MNAEigenvalueExtractor::identifyEigenvalueComponents(const CPS::IdentifiedObject::List &components) + template + void MNAEigenvalueExtractor::identifyEigenvalueComponents(const CPS::IdentifiedObject::List &components) { // TODO: [Georgii] throw exception if topology contains components that do not implement EigenvalueCompInterface for (auto comp : components) @@ -34,11 +46,17 @@ namespace DPsim if (eigenvalueComponent) { mEigenvalueComponents.push_back(eigenvalueComponent); + auto eigenvalueDynamicComponent = std::dynamic_pointer_cast>(comp); + if (eigenvalueDynamicComponent) + { + mEigenvalueDynamicComponents.push_back(eigenvalueDynamicComponent); + } } } } - void MNAEigenvalueExtractor::setBranchIndices() + template + void MNAEigenvalueExtractor::setBranchIndices() { int size = mEigenvalueComponents.size(); for (int i = 0; i < size; i++) @@ -47,7 +65,8 @@ namespace DPsim } } - void MNAEigenvalueExtractor::createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices) + template + void MNAEigenvalueExtractor::createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices) { int nBranches = mEigenvalueComponents.size(); mSignMatrix = Matrix(nBranches, nBranches); @@ -55,43 +74,108 @@ namespace DPsim mBranchNodeIncidenceMatrix = Matrix(nBranches, numMatrixNodeIndices); } - void MNAEigenvalueExtractor::stampEigenvalueMatrices() + template + void MNAEigenvalueExtractor::stampEigenvalueMatrices() { for (auto comp : mEigenvalueComponents) { - comp->stampEigenvalueMatrices(mSignMatrix, mDiscretizationMatrix, mBranchNodeIncidenceMatrix); + comp->stampBranchNodeIncidenceMatrix(mBranchNodeIncidenceMatrix); } mNodeBranchIncidenceMatrix = mBranchNodeIncidenceMatrix.transpose(); + for (auto dynamicComp : mEigenvalueDynamicComponents) + { + dynamicComp->stampSignMatrix(mSignMatrix, mCoeffDP); + dynamicComp->stampDiscretizationMatrix(mDiscretizationMatrix, mCoeffDP); + } } - void MNAEigenvalueExtractor::extractEigenvalues(const Matrix &powerSystemMatrix, Real timeStep) + template <> + void MNAEigenvalueExtractor::logInitialization() { - calculateStateMatrix(powerSystemMatrix); - computeDiscreteEigenvalues(); - recoverEigenvalues(timeStep); + SPDLOG_LOGGER_INFO(mSLog, "---- Initialize ----"); + SPDLOG_LOGGER_INFO(mSLog, "sign matrix: {}", CPS::Logger::matrixToString(mSignMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "discretization matrix: {}", CPS::Logger::matrixToString(mDiscretizationMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "branch <-> node incidence matrix: {}", CPS::Logger::matrixToString(mBranchNodeIncidenceMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "node <-> branch incidence matrix: {}", CPS::Logger::matrixToString(mNodeBranchIncidenceMatrix)); + mSLog->flush(); + } - SPDLOG_LOGGER_INFO(mSLog, "---- Extract eigenvalues ----"); - SPDLOG_LOGGER_INFO(mSLog, "discretized state matrix: {}", CPS::Logger::matrixToString(mStateMatrix)); - SPDLOG_LOGGER_INFO(mSLog, "discrete eigenvalues: {}", CPS::Logger::matrixCompToString(mDiscreteEigenvalues)); - SPDLOG_LOGGER_INFO(mSLog, "eigenvalues: {}", CPS::Logger::matrixCompToString(mEigenvalues)); + template <> + void MNAEigenvalueExtractor::logInitialization() + { + SPDLOG_LOGGER_INFO(mSLog, "---- Initialize ----"); + SPDLOG_LOGGER_INFO(mSLog, "sign matrix: {}", CPS::Logger::matrixCompToString(mSignMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "discretization matrix: {}", CPS::Logger::matrixCompToString(mDiscretizationMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "branch <-> node incidence matrix: {}", CPS::Logger::matrixToString(mBranchNodeIncidenceMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "node <-> branch incidence matrix: {}", CPS::Logger::matrixToString(mNodeBranchIncidenceMatrix)); mSLog->flush(); } - void MNAEigenvalueExtractor::calculateStateMatrix(const Matrix &powerSystemMatrix) + template + void MNAEigenvalueExtractor::extractEigenvalues(const Matrix &powerSystemMatrix) + { + calculateStateMatrix(powerSystemMatrix); + computeDiscreteEigenvalues(); + recoverEigenvalues(); + logExtraction(); + } + + template <> + void MNAEigenvalueExtractor::calculateStateMatrix(const Matrix &powerSystemMatrix) { // TODO: [Georgii] use back substitution of factorized power system matrix instead of inversion (performance) Matrix intermediateResult = powerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix; mStateMatrix = mSignMatrix + mDiscretizationMatrix * mBranchNodeIncidenceMatrix * intermediateResult; } - void MNAEigenvalueExtractor::computeDiscreteEigenvalues() + template <> + void MNAEigenvalueExtractor::calculateStateMatrix(const Matrix &powerSystemMatrix) + { + // TODO: [Georgii] use back substitution of factorized power system matrix instead of inversion (performance) + MatrixComp compPowerSystemMatrix = CPS::Math::convertToComplex(powerSystemMatrix); + MatrixComp intermediateResult = compPowerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix; + mStateMatrix = mSignMatrix + mDiscretizationMatrix * mBranchNodeIncidenceMatrix * intermediateResult; + } + + template + void MNAEigenvalueExtractor::computeDiscreteEigenvalues() { auto discreteEigenvaluesIncludingZeros = mStateMatrix.eigenvalues(); mDiscreteEigenvalues = CPS::Math::returnNonZeroElements(discreteEigenvaluesIncludingZeros); } - void MNAEigenvalueExtractor::recoverEigenvalues(Real timeStep) + 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) + 1.0j * mSystemOmega; + } + + template <> + void MNAEigenvalueExtractor::logExtraction() + { + SPDLOG_LOGGER_INFO(mSLog, "---- Extract eigenvalues ----"); + SPDLOG_LOGGER_INFO(mSLog, "discretized state matrix: {}", CPS::Logger::matrixToString(mStateMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "discrete eigenvalues: {}", CPS::Logger::matrixCompToString(mDiscreteEigenvalues)); + SPDLOG_LOGGER_INFO(mSLog, "eigenvalues: {}", CPS::Logger::matrixCompToString(mEigenvalues)); + mSLog->flush(); + } + + template <> + void MNAEigenvalueExtractor::logExtraction() { - mEigenvalues = 2.0 / timeStep * (mDiscreteEigenvalues.array() - 1.0) / (mDiscreteEigenvalues.array() + 1.0); + SPDLOG_LOGGER_INFO(mSLog, "---- Extract eigenvalues ----"); + SPDLOG_LOGGER_INFO(mSLog, "discretized state matrix: {}", CPS::Logger::matrixCompToString(mStateMatrix)); + SPDLOG_LOGGER_INFO(mSLog, "discrete eigenvalues: {}", CPS::Logger::matrixCompToString(mDiscreteEigenvalues)); + SPDLOG_LOGGER_INFO(mSLog, "eigenvalues: {}", CPS::Logger::matrixCompToString(mEigenvalues)); + mSLog->flush(); } -} \ No newline at end of file + + template class MNAEigenvalueExtractor; + template class MNAEigenvalueExtractor; +} diff --git a/dpsim/src/MNASolver.cpp b/dpsim/src/MNASolver.cpp index e6c4257de7..7fb280f9fb 100644 --- a/dpsim/src/MNASolver.cpp +++ b/dpsim/src/MNASolver.cpp @@ -351,7 +351,7 @@ void MnaSolver::identifyTopologyObjects() { template void MnaSolver::extractEigenvalues() { - mMNAEigenvalueExtractor.initialize(mSystem, mNumMatrixNodeIndices); + mMNAEigenvalueExtractor.initialize(mSystem, mNumMatrixNodeIndices, Solver::mTimeStep); } template diff --git a/dpsim/src/MNASolverDirect.cpp b/dpsim/src/MNASolverDirect.cpp index 591f4ec054..2e9c3f7bb9 100644 --- a/dpsim/src/MNASolverDirect.cpp +++ b/dpsim/src/MNASolverDirect.cpp @@ -307,7 +307,7 @@ template void MnaSolverDirect::extractEigenvalues() { MnaSolver::extractEigenvalues(); - MnaSolver::mMNAEigenvalueExtractor.extractEigenvalues(((Matrix)mSwitchedMatrices[mCurrentSwitchStatus][0]), Solver::mTimeStep); + MnaSolver::mMNAEigenvalueExtractor.extractEigenvalues(((Matrix)mSwitchedMatrices[mCurrentSwitchStatus][0])); } template diff --git a/dpsim/src/Simulation.cpp b/dpsim/src/Simulation.cpp index 15f95b2841..b48b4d80bd 100644 --- a/dpsim/src/Simulation.cpp +++ b/dpsim/src/Simulation.cpp @@ -422,10 +422,28 @@ void Simulation::extractEigenvalues() // TODO: [Georgii] throw exceptions for multiple solvers and for non-supported solver types if (mSolvers.size() == 1) { - auto mnaSolver = std::dynamic_pointer_cast>(mSolvers[0]); - if (mnaSolver) + switch (mDomain) { - mnaSolver->extractEigenvalues(); + case Domain::SP: + break; // SP domain is not supported + case Domain::DP: + { + auto mnaComplexSolver = std::dynamic_pointer_cast>(mSolvers[0]); + if (mnaComplexSolver) + { + mnaComplexSolver->extractEigenvalues(); + } + break; + } + case Domain::EMT: + { + auto mnaRealSolver = std::dynamic_pointer_cast>(mSolvers[0]); + if (mnaRealSolver) + { + mnaRealSolver->extractEigenvalues(); + } + break; + } } } }