Skip to content

Commit

Permalink
Make MNAEigenvalueExtractor a template class
Browse files Browse the repository at this point in the history
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 <georgii.tishenin@eonerc.rwth-aachen.de>
  • Loading branch information
georgii-tishenin committed Jan 16, 2024
1 parent 410da27 commit ffe67b6
Show file tree
Hide file tree
Showing 6 changed files with 156 additions and 43 deletions.
24 changes: 16 additions & 8 deletions dpsim/include/dpsim/MNAEigenvalueExtractor.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,37 +3,45 @@
#include <dpsim/Definitions.h>
#include <dpsim/DataLogger.h>
#include <dpsim-models/SystemTopology.h>
#include <dpsim-models/Solver/EigenvalueCompInterface.h>
#include <dpsim-models/Solver/EigenvalueDynamicCompInterface.h>

namespace DPsim
{
/// Extracts eigenvalues using MNA power system conductance matrix
template <typename MatrixType>
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<MatrixType>::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();
};
}
7 changes: 5 additions & 2 deletions dpsim/include/dpsim/MNASolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
namespace DPsim {
/// Solver class using Modified Nodal Analysis (MNA).
template <typename VarType>
class MnaSolver : public Solver {
class MnaSolver : public Solver {
protected:
// #### General simulation settings ####
/// Simulation domain, which can be dynamic phasor (DP) or EMT
Expand Down Expand Up @@ -123,7 +123,10 @@ namespace DPsim {
std::vector<Real> 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<std::is_same<VarType, Real>::value, Matrix, MatrixComp>::type;
MNAEigenvalueExtractor<MatrixType> mMNAEigenvalueExtractor;

/// Constructor should not be called by users but by Simulation
MnaSolver(String name,
CPS::Domain domain = CPS::Domain::DP,
Expand Down
140 changes: 112 additions & 28 deletions dpsim/src/MNAEigenvalueExtractor.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,30 +2,42 @@

namespace DPsim
{
MNAEigenvalueExtractor::MNAEigenvalueExtractor() : mSLog(CPS::Logger::get("MNAEigenvalueExtractor", CPS::Logger::Level::info, CPS::Logger::Level::info)) {}
template <typename MatrixType>
MNAEigenvalueExtractor<MatrixType>::MNAEigenvalueExtractor() : mSLog(CPS::Logger::get("MNAEigenvalueExtractor", CPS::Logger::Level::info, CPS::Logger::Level::info)) {}

MNAEigenvalueExtractor::MNAEigenvalueExtractor(const CPS::SystemTopology &topology, UInt numMatrixNodeIndices)
template <typename MatrixType>
MNAEigenvalueExtractor<MatrixType>::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 <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::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 <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::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 <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::identifyEigenvalueComponents(const CPS::IdentifiedObject::List &components)
{
// TODO: [Georgii] throw exception if topology contains components that do not implement EigenvalueCompInterface
for (auto comp : components)
Expand All @@ -34,11 +46,17 @@ namespace DPsim
if (eigenvalueComponent)
{
mEigenvalueComponents.push_back(eigenvalueComponent);
auto eigenvalueDynamicComponent = std::dynamic_pointer_cast<CPS::EigenvalueDynamicCompInterface<MatrixType>>(comp);
if (eigenvalueDynamicComponent)
{
mEigenvalueDynamicComponents.push_back(eigenvalueDynamicComponent);
}
}
}
}

void MNAEigenvalueExtractor::setBranchIndices()
template <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::setBranchIndices()
{
int size = mEigenvalueComponents.size();
for (int i = 0; i < size; i++)
Expand All @@ -47,51 +65,117 @@ namespace DPsim
}
}

void MNAEigenvalueExtractor::createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices)
template <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::createEmptyEigenvalueMatrices(UInt numMatrixNodeIndices)
{
int nBranches = mEigenvalueComponents.size();
mSignMatrix = Matrix(nBranches, nBranches);
mDiscretizationMatrix = Matrix(nBranches, nBranches);
mBranchNodeIncidenceMatrix = Matrix(nBranches, numMatrixNodeIndices);
}

void MNAEigenvalueExtractor::stampEigenvalueMatrices()
template <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::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<Matrix>::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<MatrixComp>::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 <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::extractEigenvalues(const Matrix &powerSystemMatrix)
{
calculateStateMatrix(powerSystemMatrix);
computeDiscreteEigenvalues();
recoverEigenvalues();
logExtraction();
}

template <>
void MNAEigenvalueExtractor<Matrix>::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<MatrixComp>::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 <typename MatrixType>
void MNAEigenvalueExtractor<MatrixType>::computeDiscreteEigenvalues()
{
auto discreteEigenvaluesIncludingZeros = mStateMatrix.eigenvalues();
mDiscreteEigenvalues = CPS::Math::returnNonZeroElements(discreteEigenvaluesIncludingZeros);
}

void MNAEigenvalueExtractor::recoverEigenvalues(Real timeStep)
template <>
void MNAEigenvalueExtractor<Matrix>::recoverEigenvalues()
{
mEigenvalues = 2.0 / mTimeStep * (mDiscreteEigenvalues.array() - 1.0) / (mDiscreteEigenvalues.array() + 1.0);
}

template <>
void MNAEigenvalueExtractor<MatrixComp>::recoverEigenvalues()
{
mEigenvalues = 2.0 / mTimeStep * (mDiscreteEigenvalues.array() - 1.0) / (mDiscreteEigenvalues.array() + 1.0) + 1.0j * mSystemOmega;
}

template <>
void MNAEigenvalueExtractor<Matrix>::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<MatrixComp>::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();
}
}

template class MNAEigenvalueExtractor<Matrix>;
template class MNAEigenvalueExtractor<MatrixComp>;
}
2 changes: 1 addition & 1 deletion dpsim/src/MNASolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ void MnaSolver<VarType>::identifyTopologyObjects() {
template <typename VarType>
void MnaSolver<VarType>::extractEigenvalues()
{
mMNAEigenvalueExtractor.initialize(mSystem, mNumMatrixNodeIndices);
mMNAEigenvalueExtractor.initialize(mSystem, mNumMatrixNodeIndices, Solver::mTimeStep);
}

template <typename VarType>
Expand Down
2 changes: 1 addition & 1 deletion dpsim/src/MNASolverDirect.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -307,7 +307,7 @@ template <typename VarType>
void MnaSolverDirect<VarType>::extractEigenvalues()
{
MnaSolver<VarType>::extractEigenvalues();
MnaSolver<VarType>::mMNAEigenvalueExtractor.extractEigenvalues(((Matrix)mSwitchedMatrices[mCurrentSwitchStatus][0]), Solver::mTimeStep);
MnaSolver<VarType>::mMNAEigenvalueExtractor.extractEigenvalues(((Matrix)mSwitchedMatrices[mCurrentSwitchStatus][0]));
}

template <typename VarType>
Expand Down
24 changes: 21 additions & 3 deletions dpsim/src/Simulation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<DPsim::MnaSolver<Real>>(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<DPsim::MnaSolver<Complex>>(mSolvers[0]);
if (mnaComplexSolver)
{
mnaComplexSolver->extractEigenvalues();
}
break;
}
case Domain::EMT:
{
auto mnaRealSolver = std::dynamic_pointer_cast<DPsim::MnaSolver<Real>>(mSolvers[0]);
if (mnaRealSolver)
{
mnaRealSolver->extractEigenvalues();
}
break;
}
}
}
}
Expand Down

0 comments on commit ffe67b6

Please sign in to comment.