From fa77a6917a513308b823896182f82f369f526640 Mon Sep 17 00:00:00 2001 From: Georgii Tishenin Date: Mon, 24 Jun 2024 17:50:12 +0200 Subject: [PATCH] Add a separate EigenvalueUtils class Signed-off-by: Georgii Tishenin --- dpsim-models/include/dpsim-models/MathUtils.h | 4 -- dpsim-models/src/MathUtils.cpp | 41 ----------------- dpsim/include/dpsim/EigenvalueUtils.h | 14 ++++++ dpsim/include/dpsim/MNAEigenvalueExtractor.h | 1 + dpsim/src/CMakeLists.txt | 1 + dpsim/src/EigenvalueUtils.cpp | 45 +++++++++++++++++++ dpsim/src/MNAEigenvalueExtractor.cpp | 4 +- 7 files changed, 63 insertions(+), 47 deletions(-) create mode 100644 dpsim/include/dpsim/EigenvalueUtils.h create mode 100644 dpsim/src/EigenvalueUtils.cpp diff --git a/dpsim-models/include/dpsim-models/MathUtils.h b/dpsim-models/include/dpsim-models/MathUtils.h index ddfb920734..0ffdb38e6a 100644 --- a/dpsim-models/include/dpsim-models/MathUtils.h +++ b/dpsim-models/include/dpsim-models/MathUtils.h @@ -63,10 +63,6 @@ class Math { // | Interharmonics harm1-harm2 | Re(row,col)_harm2 | Re(row,col)_harm2 | // | Interharmonics harm1-harm2 | Im(row,col)_harm2 | Im(row,col)_harm2 | - static MatrixComp returnNonZeroElements(const MatrixComp &mat); - - static MatrixComp convertRealEquivalentToComplexMatrix(const Matrix &realEquivalentMatrix); - static void setMatrixElement(SparseMatrixRow &mat, Matrix::Index row, Matrix::Index column, Complex value, Int maxFreq = 1, Int freqIdx = 0); diff --git a/dpsim-models/src/MathUtils.cpp b/dpsim-models/src/MathUtils.cpp index 5ff0f62992..4439464f3d 100644 --- a/dpsim-models/src/MathUtils.cpp +++ b/dpsim-models/src/MathUtils.cpp @@ -97,47 +97,6 @@ Real Math::realFromVectorElement(const Matrix &mat, Matrix::Index row) { return mat(row, 0); } -// TODO: consider resizing existing matrix instead of creating a new one (performance, memory) -MatrixComp Math::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 Math::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; -} - void Math::setMatrixElement(SparseMatrixRow &mat, Matrix::Index row, Matrix::Index column, Complex value, Int maxFreq, Int freqIdx) { 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 index 21ecad20ee..236f051559 100644 --- a/dpsim/include/dpsim/MNAEigenvalueExtractor.h +++ b/dpsim/include/dpsim/MNAEigenvalueExtractor.h @@ -2,6 +2,7 @@ #include #include +#include #include namespace DPsim { diff --git a/dpsim/src/CMakeLists.txt b/dpsim/src/CMakeLists.txt index 999d85a82d..7ca7a0a91d 100644 --- a/dpsim/src/CMakeLists.txt +++ b/dpsim/src/CMakeLists.txt @@ -5,6 +5,7 @@ set(DPSIM_SOURCES 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 index c5aa5f8ab8..4561482a6e 100644 --- a/dpsim/src/MNAEigenvalueExtractor.cpp +++ b/dpsim/src/MNAEigenvalueExtractor.cpp @@ -110,7 +110,7 @@ void MNAEigenvalueExtractor::calculateStateMatrix( const Matrix &powerSystemMatrix) { // TODO: use right hand side solving of factorized power system matrix instead of inversion (performance). MatrixComp compPowerSystemMatrix = - CPS::Math::convertRealEquivalentToComplexMatrix(powerSystemMatrix); + EigenvalueUtils::convertRealEquivalentToComplexMatrix(powerSystemMatrix); MatrixComp intermediateResult = compPowerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix; mStateMatrix = mSignMatrix - mDiscretizationMatrix * @@ -122,7 +122,7 @@ template void MNAEigenvalueExtractor::computeDiscreteEigenvalues() { auto discreteEigenvaluesIncludingZeros = mStateMatrix.eigenvalues(); **mDiscreteEigenvalues = - CPS::Math::returnNonZeroElements(discreteEigenvaluesIncludingZeros); + EigenvalueUtils::returnNonZeroElements(discreteEigenvaluesIncludingZeros); // TODO: filter out eigenvalues = -1 + 0i to avoid division by zero in recoverEigenvalues(). }