diff --git a/dpsim-models/include/dpsim-models/Components.h b/dpsim-models/include/dpsim-models/Components.h index d92170fe28..14305532ed 100644 --- a/dpsim-models/include/dpsim-models/Components.h +++ b/dpsim-models/include/dpsim-models/Components.h @@ -82,6 +82,7 @@ #endif #include #include +#include #include #include #include diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Diode.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Diode.h new file mode 100644 index 0000000000..49a592d02f --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_Diode.h @@ -0,0 +1,82 @@ +/* Copyright 2017-2021 Institute for Automation of Complex Power Systems, + * EONERC, RWTH Aachen University + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + *********************************************************************************/ + +#pragma once + +#include +#include +//#include +#include + +namespace CPS { +namespace EMT { +namespace Ph1 { +/// EMT Diode +class Diode : public MNASimPowerComp, + public SharedFactory, + public MNANonlinearVariableCompInterface { +protected: + ///TODO: mI_S and mV_T as const CPS::Attribute::Ptr, also change in + // pybind EMTComponents + Real mI_S = 0.000001; + Real mV_T = 0.027; + Matrix Jacobian = Matrix::Zero(1, 1); + Real itVoltage = 0.; + Real itCurrent = 0.; + +public: + /// Defines UID, name, component parameters and logging level + Diode(String uid, String name, Logger::Level logLevel = Logger::Level::off); + /// Defines name, component parameters and logging level + Diode(String name, Logger::Level logLevel = Logger::Level::off) + : Diode(name, name, logLevel) {} + + // #### General #### + /// + SimPowerComp::Ptr clone(String name) override; + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + + void setParameters(Real, Real); + + // #### MNA section #### + /// Initializes internal variables of the component + void mnaCompInitialize(Real omega, Real timeStep, + Attribute::Ptr leftSideVector) override; + /// Stamps system matrix + void mnaCompApplySystemMatrixStamp(SparseMatrixRow &systemMatrix) override; + /// Stamps right side (source) vector + void mnaCompApplyRightSideVectorStamp(Matrix &rightVector) override { + } //No right side vector stamps + /// Update interface voltage from MNA system result + void mnaCompUpdateVoltage(const Matrix &leftVector) override; + /// Update interface current from MNA system result + void mnaCompUpdateCurrent(const Matrix &leftVector) override; + /// MNA pre and post step operations + void mnaCompPreStep(Real time, + Int timeStepCount) override; //No right side vector stamps + void mnaCompPostStep(Real time, Int timeStepCount, + 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) override; + + virtual void iterationUpdate(const Matrix &leftVector) override; + + virtual bool hasParameterChanged() override { return true; } + + void calculateNonlinearFunctionResult() override; + + void updateJacobian(); +}; +} // namespace Ph1 +} // namespace EMT +} // namespace CPS diff --git a/dpsim-models/include/dpsim-models/Solver/MNANonlinearVariableCompInterface.h b/dpsim-models/include/dpsim-models/Solver/MNANonlinearVariableCompInterface.h new file mode 100644 index 0000000000..be466000ef --- /dev/null +++ b/dpsim-models/include/dpsim-models/Solver/MNANonlinearVariableCompInterface.h @@ -0,0 +1,36 @@ +/* Copyright 2017-2021 Institute for Automation of Complex Power Systems, + * EONERC, RWTH Aachen University + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + *********************************************************************************/ + +#pragma once + +#include +#include +#include + +namespace CPS { +/// MNA interface to be used by nonlinear elements that require recomputing +// of the system matrix +class MNANonlinearVariableCompInterface + : virtual public MNAVariableCompInterface { +public: + typedef std::shared_ptr NonlinearPtr; + typedef std::vector NonlinearList; + + const Attribute::Ptr mNonlinearFunctionStamp; + std::vector> mNonlinearVariableSystemMatrixEntries; + + /// Returns value of the component's defining voltage/current equation + // based on current circuit quantities + virtual void calculateNonlinearFunctionResult() = 0; + virtual void iterationUpdate(const Matrix &leftVector) = 0; + +protected: + MNANonlinearVariableCompInterface() + : mNonlinearFunctionStamp(AttributeStatic::make()){}; +}; +} // namespace CPS diff --git a/dpsim-models/src/CMakeLists.txt b/dpsim-models/src/CMakeLists.txt index 795c44cd97..79b1f6a047 100644 --- a/dpsim-models/src/CMakeLists.txt +++ b/dpsim-models/src/CMakeLists.txt @@ -70,6 +70,7 @@ list(APPEND MODELS_SOURCES EMT/EMT_Ph1_VoltageSourceRamp.cpp EMT/EMT_Ph1_VoltageSourceNorton.cpp EMT/EMT_Ph1_Switch.cpp + EMT/EMT_Ph1_Diode.cpp EMT/EMT_Ph3_CurrentSource.cpp EMT/EMT_Ph3_VoltageSource.cpp diff --git a/dpsim-models/src/EMT/EMT_Ph1_Diode.cpp b/dpsim-models/src/EMT/EMT_Ph1_Diode.cpp new file mode 100644 index 0000000000..c568e3a977 --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph1_Diode.cpp @@ -0,0 +1,166 @@ +/* Copyright 2017-2021 Institute for Automation of Complex Power Systems, + * EONERC, RWTH Aachen University + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + ******************************************************************************/ + +#include + +using namespace CPS; + +EMT::Ph1::Diode::Diode(String uid, String name, Logger::Level logLevel) + : MNASimPowerComp(uid, name, false, false, logLevel) { + mPhaseType = PhaseType::Single; + setTerminalNumber(2); + **mIntfVoltage = Matrix::Zero(1, 1); + **mIntfCurrent = Matrix::Zero(1, 1); +} + +void EMT::Ph1::Diode::setParameters(Real I_S, Real V_T) { + mI_S = I_S; + mV_T = V_T; +} + +SimPowerComp::Ptr EMT::Ph1::Diode::clone(String name) { + auto copy = Diode::make(name, mLogLevel); + copy->setParameters(mI_S, mV_T); + return copy; +} + +void EMT::Ph1::Diode::initializeFromNodesAndTerminals(Real frequency) { + + // IntfVoltage initialization for each phase + MatrixComp vInitABC = Matrix::Zero(1, 1); + vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(1) - + RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); + (**mIntfVoltage)(0, 0) = vInitABC(0, 0).real(); + + ///FIXME: Initialization should include solving the system once to obtain + // the actual values solving the + // system for the 0th time step. As of now abnormal current + // values for the 0th time step are to be expected. + + (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.); + + SPDLOG_LOGGER_INFO(mSLog, + "\n--- Initialization from powerflow ---" + "\nVoltage across: {:f}" + "\nCurrent: {:f}" + "\nTerminal 0 voltage: {:f}" + "\nTerminal 1 voltage: {:f}" + "\n--- Initialization from powerflow finished ---", + (**mIntfVoltage)(0, 0), (**mIntfCurrent)(0, 0), + initialSingleVoltage(0).real(), + initialSingleVoltage(1).real()); +} + +void EMT::Ph1::Diode::mnaCompInitialize(Real omega, Real timeStep, + Attribute::Ptr leftVector) { + + updateMatrixNodeIndices(); + + **mRightVector = Matrix::Zero(leftVector->get().rows(), 1); + **mNonlinearFunctionStamp = Matrix::Zero(leftVector->get().rows(), 1); + calculateNonlinearFunctionResult(); + updateJacobian(); + + mMnaTasks.push_back(std::make_shared(*this, leftVector)); +} + +void EMT::Ph1::Diode::mnaCompApplySystemMatrixStamp( + SparseMatrixRow &systemMatrix) { + // Set diagonal entries + if (terminalNotGrounded(0)) { + // set upper left block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), + matrixNodeIndex(0, 0), Jacobian(0, 0)); + } + if (terminalNotGrounded(1)) { + // set buttom right block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), + matrixNodeIndex(1, 0), Jacobian(0, 0)); + } + // Set off diagonal blocks, 2x3x3 entries + if (terminalNotGrounded(0) && terminalNotGrounded(1)) { + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), + matrixNodeIndex(1, 0), -Jacobian(0, 0)); + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), + matrixNodeIndex(0, 0), -Jacobian(0, 0)); + } +} + +void EMT::Ph1::Diode::mnaCompAddPostStepDependencies( + AttributeBase::List &prevStepDependencies, + AttributeBase::List &attributeDependencies, + AttributeBase::List &modifiedAttributes, + Attribute::Ptr &leftVector) { + attributeDependencies.push_back(leftVector); + modifiedAttributes.push_back(attribute("v_intf")); + modifiedAttributes.push_back(attribute("i_intf")); +} + +void EMT::Ph1::Diode::mnaCompPreStep(Real time, Int timeStepCount) { + mnaCompApplyRightSideVectorStamp(**mRightVector); +} + +void EMT::Ph1::Diode::mnaCompPostStep(Real time, Int timeStepCount, + Attribute::Ptr &leftVector) { + mnaUpdateVoltage(**leftVector); + mnaUpdateCurrent(**leftVector); +} + +void EMT::Ph1::Diode::mnaCompUpdateVoltage(const Matrix &leftVector) { + // v1 - v0 + **mIntfVoltage = Matrix::Zero(3, 1); + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = + (**mIntfVoltage)(0, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + } +} + +void EMT::Ph1::Diode::mnaCompUpdateCurrent(const Matrix &leftVector) { + (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.); +} + +void EMT::Ph1::Diode::iterationUpdate(const Matrix &leftVector) { + //Update phase voltages + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = + Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = + (**mIntfVoltage)(0, 0) - + Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + } + + calculateNonlinearFunctionResult(); + + updateJacobian(); +} + +void EMT::Ph1::Diode::calculateNonlinearFunctionResult() { + + (**mIntfCurrent)(0, 0) = mI_S * (expf((**mIntfVoltage)(0, 0) / mV_T) - 1.); + + if (terminalNotGrounded(1)) { + Math::setVectorElement(**mNonlinearFunctionStamp, matrixNodeIndex(1, 0), + (**mIntfCurrent)(0, 0)); + } + if (terminalNotGrounded(0)) { + Math::setVectorElement(**mNonlinearFunctionStamp, matrixNodeIndex(0, 0), + -(**mIntfCurrent)(0, 0)); + } +} + +void EMT::Ph1::Diode::updateJacobian() { + Jacobian(0, 0) = (mI_S / mV_T) * expf((**mIntfVoltage)(0, 0) / mV_T); +} \ No newline at end of file diff --git a/dpsim/examples/cxx/CMakeLists.txt b/dpsim/examples/cxx/CMakeLists.txt index dba9ec1117..f5ec131767 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_Ph1_Diode_test.cpp # EMT examples with PF initialization Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp diff --git a/dpsim/examples/cxx/Circuits/EMT_Ph1_Diode_test.cpp b/dpsim/examples/cxx/Circuits/EMT_Ph1_Diode_test.cpp new file mode 100644 index 0000000000..d078c1e45a --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_Ph1_Diode_test.cpp @@ -0,0 +1,59 @@ +#include "../Examples.h" +#include + +using namespace DPsim; +using namespace CPS::EMT; + +void EMT_Ph1_Diode() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simName = "EMT_Ph1_Diode_test"; + Logger::setLogDir("logs/" + simName); + + // Nodes + auto n1 = SimNode::make("n1", PhaseType::Single); + auto n2 = SimNode::make("n1", PhaseType::Single); + + // Components + + auto vs0 = Ph1::VoltageSource::make("vs0"); + vs0->setParameters(CPS::Complex(1., 0.), 50.0); + + auto load = CPS::EMT::Ph1::Resistor::make("Load"); + load->setParameters(10.); + + auto diode = Ph1::Diode::make("Diode"); + + // Topology + load->connect(SimNode::List{n1, n2}); + + diode->connect(SimNode::List{n2, SimNode::GND}); //SimNode::GND, n2 + + vs0->connect(SimNode::List{SimNode::GND, n1}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2}, + SystemComponentList{load, vs0, diode}); + + // Logging + auto logger = DataLogger::make(simName); + logger->logAttribute("I_Diode", diode->attribute("i_intf")); + logger->logAttribute("V_Diode", diode->attribute("v_intf")); + + Simulation sim(simName, Logger::Level::info); + sim.doInitFromNodesAndTerminals(true); + sim.setSystem(sys); + sim.addLogger(logger); + sim.setDomain(Domain::EMT); + sim.setSolverType(Solver::Type::ITERATIVEMNA); + sim.setDirectLinearSolverConfiguration(DirectLinearSolverConfiguration()); + sim.setTimeStep(timeStep); + sim.setFinalTime(finalTime); + sim.run(); +} + +int main() { + EMT_Ph1_Diode(); + return 0; +} \ No newline at end of file diff --git a/dpsim/include/dpsim/IterativeMnaSolverDirect.h b/dpsim/include/dpsim/IterativeMnaSolverDirect.h new file mode 100644 index 0000000000..1bde5db094 --- /dev/null +++ b/dpsim/include/dpsim/IterativeMnaSolverDirect.h @@ -0,0 +1,101 @@ +/* Copyright 2017-2021 Institute for Automation of Complex Power Systems, + * EONERC, RWTH Aachen University + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + *********************************************************************************/ + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#ifdef WITH_KLU +#include +#endif +#include +#ifdef WITH_CUDA +#include +#ifdef WITH_CUDA_SPARSE +#include +#endif +#ifdef WITH_MAGMA +#include +#endif +#endif +#include +#include +#include +#include +#include +#include +#include + +namespace DPsim { + +template +class IterativeMnaSolverDirect : public MnaSolverDirect { + +public: + IterativeMnaSolverDirect( + String name, CPS::Domain domain = CPS::Domain::EMT, + CPS::Logger::Level logLevel = CPS::Logger::Level::info) + : MnaSolverDirect(name, domain, logLevel) { + MnaSolverDirect::mImplementationInUse = + DirectLinearSolverImpl::KLU; + } + + virtual ~IterativeMnaSolverDirect() = default; + + virtual void initialize() override; + +protected: + // Delta of iterative Solutions + Matrix mLeftStep; + // Number of Iterations per step + unsigned iterations = 0; + // If mRightSideVector deviates less than Epsilon per element from the + // result of the system defining node equations, mesh equations + // and auxhiliary equations (calculationError), + // the solution is good enough + bool isConverged = true; + Real Epsilon = 0.0001; + Matrix calculationError; + Real calculationErrorElement; + + // List of all nonlinear component function contributions + std::vector mNonlinearFunctionStamps; + // Sum of all nonlinear component function contributions + Matrix mNonlinearFunctionResult = Matrix::Zero(1, 1); + + CPS::MNANonlinearVariableCompInterface::NonlinearList + mMNANonlinearVariableComponents; + + using MnaSolverDirect::mVariableSystemMatrix; + using MnaSolverDirect::mBaseSystemMatrix; + using MnaSolverDirect::mMNAIntfSwitches; + using MnaSolverDirect::mDirectLinearSolverVariableSystemMatrix; + using MnaSolverDirect::mRecomputationTimes; + using MnaSolverDirect::mNumRecomputations; + using MnaSolverDirect::mMNAIntfVariableComps; + using MnaSolverDirect::mRightSideVector; + + using MnaSolver::mNumMatrixNodeIndices; + + void solveWithSystemMatrixRecomputation(Real time, + Int timeStepCount) override; + virtual void recomputeSystemMatrix(Real time) override; + std::shared_ptr createSolveTaskRecomp() override; +}; +} // namespace DPsim \ No newline at end of file diff --git a/dpsim/include/dpsim/Simulation.h b/dpsim/include/dpsim/Simulation.h index 4d58db8dd0..b7e2ba2c0f 100644 --- a/dpsim/include/dpsim/Simulation.h +++ b/dpsim/include/dpsim/Simulation.h @@ -141,7 +141,9 @@ class Simulation : public CPS::AttributeList { template void createSolvers(); /// Subroutine for MNA only because there are many MNA options template void createMNASolver(); - /// Prepare schedule for simulation + + template void createIterativeMNASolver(); + // Prepare schedule for simulation void prepSchedule(); /// ### SynGen Interface ### diff --git a/dpsim/include/dpsim/Solver.h b/dpsim/include/dpsim/Solver.h index 72f12ddc68..56823b1b4e 100644 --- a/dpsim/include/dpsim/Solver.h +++ b/dpsim/include/dpsim/Solver.h @@ -75,7 +75,7 @@ class Solver { // #### Solver settings #### /// Solver types: /// Modified Nodal Analysis, Differential Algebraic, Newton Raphson - enum class Type { MNA, DAE, NRP }; + enum class Type { MNA, ITERATIVEMNA, DAE, NRP }; /// void setTimeStep(Real timeStep) { mTimeStep = timeStep; } /// @@ -124,4 +124,4 @@ class Solver { mMaxIterations = maxIterations; } }; -} // namespace DPsim +} // namespace DPsim \ No newline at end of file diff --git a/dpsim/src/CMakeLists.txt b/dpsim/src/CMakeLists.txt index 0f0460a280..f74732fc29 100644 --- a/dpsim/src/CMakeLists.txt +++ b/dpsim/src/CMakeLists.txt @@ -3,6 +3,7 @@ set(DPSIM_SOURCES RealTimeSimulation.cpp MNASolver.cpp MNASolverDirect.cpp + IterativeMnaSolverDirect.cpp DenseLUAdapter.cpp SparseLUAdapter.cpp DirectLinearSolverConfiguration.cpp diff --git a/dpsim/src/IterativeMnaSolverDirect.cpp b/dpsim/src/IterativeMnaSolverDirect.cpp new file mode 100644 index 0000000000..1386f5e275 --- /dev/null +++ b/dpsim/src/IterativeMnaSolverDirect.cpp @@ -0,0 +1,195 @@ +/* Copyright 2017-2021 Institute for Automation of Complex Power Systems, + * EONERC, RWTH Aachen University + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + ******************************************************************************/ + +#include +#include + +using namespace DPsim; +using namespace CPS; + +template +std::shared_ptr +IterativeMnaSolverDirect::createSolveTaskRecomp() { + return std::make_shared< + typename IterativeMnaSolverDirect::SolveTaskRecomp>(*this); +} + +template +void IterativeMnaSolverDirect::recomputeSystemMatrix(Real time) { + // Start from base matrix + mVariableSystemMatrix = mBaseSystemMatrix; + + // Now stamp switches into matrix + for (auto sw : mMNAIntfSwitches) + sw->mnaApplySystemMatrixStamp(mVariableSystemMatrix); + + // Now stamp variable elements into matrix + for (auto comp : mMNAIntfVariableComps) + comp->mnaApplySystemMatrixStamp(mVariableSystemMatrix); + + auto start = std::chrono::steady_clock::now(); + mDirectLinearSolverVariableSystemMatrix->factorize(mVariableSystemMatrix); + auto end = std::chrono::steady_clock::now(); + std::chrono::duration diff = end - start; + mRecomputationTimes.push_back(diff.count()); + + ++mNumRecomputations; +} + +template +void IterativeMnaSolverDirect::solveWithSystemMatrixRecomputation( + Real time, Int timeStepCount) { + // Reset source vector + this->mRightSideVector.setZero(); + // Add together the right side vector (computed by the components' + // pre-step tasks) + for (auto stamp : this->mRightVectorStamps) + this->mRightSideVector += *stamp; + + //Number of Iterations per step + unsigned iterations = 0; + do { + // Get switch and variable comp status and update system matrix and lu + //factorization accordingly + if (this->hasVariableComponentChanged()) { + this->recomputeSystemMatrix(time); + } + + // Calculate new delta solution vector: + // systemMatrix*leftStep = mRightSideVector-mNonlinearSSNfunctionResult + // Corresponds to Newton-Raphson: + // + // f(x2) =~ f(x1) + Df(x1)*(x2-x1) + // + // f(x) is the function vector containing the system describing + // equations WITHOUT injections and SSN history terms + // These are the right side source vector and are set equal to the + // function: + // f(x2) = mRightSideVector =~ f(x1) + Df(x1)*(x2-x1) + // Subtracting the past function value leaves: + // systemMatrix*leftStep = mRightSideVector-(mBaseSystemMatrix*(**mLeftSideVector)+ mNonlinearSSNfunctionResult) + // (mBaseSystemMatrix*(**mLeftSideVector)+ mNonlinearSSNfunctionResult) = f(x1), + // that is the old mRightSideVector of the previous step + // mRightSideVector is stamped by mna-prestep tasks and only updated each + // step, not iteration. + std::cout << mVariableSystemMatrix << std::endl << std::endl; + std::cout << **(this->mLeftSideVector) << std::endl << std::endl; + // mDirectLinearSolvers.solve takes a Matrix reference, + // thus we need a temporary variable as an argument + Matrix temp = (this->mRightSideVector) - + ((this->mBaseSystemMatrix) * (**(this->mLeftSideVector)) + + mNonlinearFunctionResult); + mLeftStep = mDirectLinearSolverVariableSystemMatrix->solve(temp); + + // x2 = x1 + (x2-x1) + **(this->mLeftSideVector) += mLeftStep; + + // *Update all CURRENT nonexplicit states that have a defining equation + // in the system matrix since they could not be expressed as a function + // of system inputs due to nonlinearity + + // *Update all CURRENT states + + // *Update system Jacobian with new system solution + // (including new nonexplicit states) + + // *Calculate the actual system function result with + // the new solution vector + + for (auto comp : mMNANonlinearVariableComponents) + comp->iterationUpdate(**(this->mLeftSideVector)); + + // Collect all System equation contributions from + // nonlinear SSN components + mNonlinearFunctionResult -= mNonlinearFunctionResult; + for (auto stamp : mNonlinearFunctionStamps) + mNonlinearFunctionResult += *stamp; + + // Check Convergence: + // Is the deviation of the system function result + // with the new solution vector + // + // (mBaseSystemMatrix * **mLeftSideVector + mNonlinearSSNfunctionResult) + // + // with mBaseSystemMatrix * **mLeftSideVector the contribution of + // static, linear components such as capacitors, inductors, resistors + // (Jacobian times vector is not only approximation but + // actual function of those components) + // + // and mNonlinearSSNfunctionResult the non-approximated, + // actual contribution of nonlinear SSN components + // using the new Solution vector + // + // from the system excitation in the source vector + // mRightSideVector small enough? If yes: Newton-Raphson has converged. + + calculationError = this->mRightSideVector - + (this->mBaseSystemMatrix * (**(this->mLeftSideVector)) + + mNonlinearFunctionResult); + + isConverged = true; + for (int i = 0; i < calculationError.rows(); i++) { + calculationErrorElement = calculationError(i, 0); + //std::cout << calculationErrorElement << std::endl; + if (abs(calculationErrorElement) >= Epsilon) { + isConverged = false; + break; + } + } + iterations++; + } while (!isConverged && iterations < 100); + //if(iterations > 1) std::cout << iterations << std::endl; + //std::cout << mVariableSystemMatrix << std::endl; + /// TODO: split into separate task? + //(dependent on x, updating all v attributes) + for (UInt nodeIdx = 0; nodeIdx < this->mNumNetNodes; ++nodeIdx) + this->mNodes[nodeIdx]->mnaUpdateVoltage(**(this->mLeftSideVector)); +} + +template +void IterativeMnaSolverDirect::initialize() { + this->mFrequencyParallel = false; + this->mSystemMatrixRecomputation = true; + MnaSolver::initialize(); + + mMNANonlinearVariableComponents.clear(); + mMNANonlinearVariableComponents.shrink_to_fit(); + + mNonlinearFunctionStamps.clear(); + mNonlinearFunctionStamps.shrink_to_fit(); + + for (auto comp : MnaSolver::mSystem.mComponents) { + auto MNANonlinearComp = + std::dynamic_pointer_cast(comp); + if (MNANonlinearComp) + mMNANonlinearVariableComponents.push_back(MNANonlinearComp); + } + + for (auto comp : this->mMNANonlinearVariableComponents) { + const Matrix &stamp = comp->mNonlinearFunctionStamp.get()->get(); + if (stamp.size() != 0) { + this->mNonlinearFunctionStamps.push_back(&stamp); + } + } + + mNonlinearFunctionResult = Matrix::Zero(this->mRightSideVector.rows(), 1); + for (auto stamp : mNonlinearFunctionStamps) + mNonlinearFunctionResult += *stamp; + + // Delta of iterative Solutions + mLeftStep = Matrix::Zero(this->mLeftSideVector->get().rows(), 1); + + //If mRightSideVector deviates less than Epsilon per element + //from the result of the system defining node equations, mesh equations + //and auxhiliary equations (calculationError), the solution is good enough + calculationError = Matrix::Zero(this->mRightSideVector.rows(), 1); + calculationErrorElement = 0.; +} + +template class DPsim::IterativeMnaSolverDirect; +template class DPsim::IterativeMnaSolverDirect; \ No newline at end of file diff --git a/dpsim/src/Simulation.cpp b/dpsim/src/Simulation.cpp index 1d79b4f8f4..668e22345a 100644 --- a/dpsim/src/Simulation.cpp +++ b/dpsim/src/Simulation.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -91,6 +92,9 @@ template void Simulation::createSolvers() { case Solver::Type::MNA: createMNASolver(); break; + case Solver::Type::ITERATIVEMNA: + createIterativeMNASolver(); + break; #ifdef WITH_SUNDIALS case Solver::Type::DAE: solver = std::make_shared(**mName, mSystem, **mTimeStep, 0.0); @@ -170,6 +174,31 @@ template void Simulation::createMNASolver() { } } +template void Simulation::createIterativeMNASolver() { + Solver::Ptr solver; + + std::shared_ptr> Itsolver = + std::make_shared>(**mName, mDomain, + mLogLevel); + Itsolver->setDirectLinearSolverImplementation( + DirectLinearSolverImpl::SparseLU); + + Itsolver->setTimeStep(**mTimeStep); + Itsolver->doSteadyStateInit(**mSteadyStateInit); + Itsolver->doFrequencyParallelization(false); + Itsolver->setSteadStIniTimeLimit(mSteadStIniTimeLimit); + Itsolver->setSteadStIniAccLimit(mSteadStIniAccLimit); + Itsolver->setSystem(mSystem); + Itsolver->setSolverAndComponentBehaviour(mSolverBehaviour); + Itsolver->doInitFromNodesAndTerminals(mInitFromNodesAndTerminals); + Itsolver->doSystemMatrixRecomputation(true); + Itsolver->setDirectLinearSolverConfiguration( + mDirectLinearSolverConfiguration); + Itsolver->initialize(); + + mSolvers.push_back(Itsolver); +} + void Simulation::sync() const { SPDLOG_LOGGER_INFO(mLog, "Start synchronization with remotes on interfaces"); diff --git a/dpsim/src/pybind/EMTComponents.cpp b/dpsim/src/pybind/EMTComponents.cpp index 92afccbe00..17441ab2ea 100644 --- a/dpsim/src/pybind/EMTComponents.cpp +++ b/dpsim/src/pybind/EMTComponents.cpp @@ -102,13 +102,27 @@ void addEMTPh1Components(py::module_ mEMTPh1) { .def("connect", &CPS::EMT::Ph1::Inductor::connect) .def_property("L", createAttributeGetter("L"), createAttributeSetter("L")); - - py::class_, CPS::SimPowerComp, CPS::Base::Ph1::Switch>(mEMTPh1, "Switch", py::multiple_inheritance()) - .def(py::init(), "name"_a, "loglevel"_a = CPS::Logger::Level::off) - .def("set_parameters", &CPS::EMT::Ph1::Switch::setParameters, "open_resistance"_a, "closed_resistance"_a, "closed"_a = false) // cppcheck-suppress assignBoolToPointer + + py::class_, + CPS::SimPowerComp, CPS::Base::Ph1::Switch>( + mEMTPh1, "Switch", py::multiple_inheritance()) + .def(py::init(), "name"_a, + "loglevel"_a = CPS::Logger::Level::off) + .def("set_parameters", &CPS::EMT::Ph1::Switch::setParameters, + "open_resistance"_a, "closed_resistance"_a, + "closed"_a = false) // cppcheck-suppress assignBoolToPointer .def("open", &CPS::EMT::Ph1::Switch::open) .def("close", &CPS::EMT::Ph1::Switch::close) - .def("connect", &CPS::EMT::Ph1::Switch::connect); + .def("connect", &CPS::EMT::Ph1::Switch::connect); + + py::class_, + CPS::SimPowerComp>(mEMTPh1, "Diode", + py::multiple_inheritance()) + .def(py::init()) + .def(py::init()) + .def("set_parameters", &CPS::EMT::Ph1::Diode::setParameters, "I_S"_a, + "V_T"_a) + .def("connect", &CPS::EMT::Ph1::Diode::connect); } void addEMTPh3Components(py::module_ mEMTPh3) { @@ -128,7 +142,6 @@ void addEMTPh3Components(py::module_ mEMTPh3) { .def_property("f_src", createAttributeGetter("f_src"), createAttributeSetter("f_src")); - py::class_, CPS::SimPowerComp>(mEMTPh3, "Resistor", py::multiple_inheritance()) @@ -136,7 +149,6 @@ void addEMTPh3Components(py::module_ mEMTPh3) { .def(py::init()) .def("set_parameters", &CPS::EMT::Ph3::Resistor::setParameters, "R"_a) .def("connect", &CPS::EMT::Ph3::Resistor::connect); - py::class_, @@ -194,7 +206,7 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def("set_parameters", &CPS::EMT::Ph3::Switch::setParameters, "open_resistance"_a, "closed_resistance"_a, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "closed"_a = false) .def("open", &CPS::EMT::Ph3::Switch::openSwitch) .def("close", &CPS::EMT::Ph3::Switch::closeSwitch) @@ -369,7 +381,7 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def(py::init(), "uid"_a, "name"_a, "loglevel"_a = CPS::Logger::Level::off, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "with_trafo"_a = false) .def("set_parameters", &CPS::EMT::Ph3::AvVoltageSourceInverterDQ::setParameters, @@ -402,9 +414,8 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def(py::init(), "uid"_a, "name"_a, "loglevel"_a = CPS::Logger::Level::off, - // cppcheck-suppress assignBoolToPointer - "with_resistive_losses"_a = - false) + // cppcheck-suppress assignBoolToPointer + "with_resistive_losses"_a = false) .def("set_parameters", &CPS::EMT::Ph3::Transformer::setParameters, "nom_voltage_end_1"_a, "nom_voltage_end_2"_a, "rated_power"_a, "ratio_abs"_a, "ratio_phase"_a, "resistance"_a, "inductance"_a) @@ -428,7 +439,7 @@ void addEMTPh3Components(py::module_ mEMTPh3) { "loglevel"_a = CPS::Logger::Level::off) .def("set_parameters", &CPS::EMT::Ph3::SeriesSwitch::setParameters, "open_resistance"_a, "closed_resistance"_a, - // cppcheck-suppress assignBoolToPointer + // cppcheck-suppress assignBoolToPointer "closed"_a = false) .def("open", &CPS::EMT::Ph3::SeriesSwitch::open) .def("close", &CPS::EMT::Ph3::SeriesSwitch::close) diff --git a/dpsim/src/pybind/main.cpp b/dpsim/src/pybind/main.cpp index a1392119e1..9905c9a42b 100644 --- a/dpsim/src/pybind/main.cpp +++ b/dpsim/src/pybind/main.cpp @@ -99,6 +99,7 @@ PYBIND11_MODULE(dpsimpy, m) { py::enum_(m, "Solver") .value("MNA", DPsim::Solver::Type::MNA) + .value("ITERATIVEMNA", DPsim::Solver::Type::ITERATIVEMNA) .value("DAE", DPsim::Solver::Type::DAE) .value("NRP", DPsim::Solver::Type::NRP); diff --git a/examples/Notebooks/Circuits/EMT_Ph1_Diode.ipynb b/examples/Notebooks/Circuits/EMT_Ph1_Diode.ipynb new file mode 100644 index 0000000000..a857925926 --- /dev/null +++ b/examples/Notebooks/Circuits/EMT_Ph1_Diode.ipynb @@ -0,0 +1,125 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# EMT Ph3 Linear Circuit SSN" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run C++ examples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import subprocess\n", + "\n", + "#%matplotlib widget\n", + "\n", + "name = 'EMT_Ph1_Diode_test'\n", + "\n", + "dpsim_path = subprocess.Popen(['git', 'rev-parse', '--show-toplevel'], stdout=subprocess.PIPE).communicate()[0].rstrip().decode('utf-8')\n", + "\n", + "path_exec = dpsim_path + '/build/dpsim/examples/cxx/'\n", + "sim = subprocess.Popen([path_exec + name], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", + "print(sim.communicate()[0].decode())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from villas.dataprocessing.readtools import *\n", + "from villas.dataprocessing.timeseries import *\n", + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import matplotlib.pyplot as plt\n", + "import re\n", + "import numpy as np\n", + "import math\n", + "\n", + "work_dir = os.getcwd() + \"/logs/\"\n", + "path_logfile = work_dir + name + '/' + name + '.csv'\n", + "ts_dpsim_EMT = read_timeseries_dpsim(path_logfile)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.close('all')\n", + "fig1 = plt.figure()\n", + "\n", + "plt.plot(ts_dpsim_EMT['V_Diode'].time, ts_dpsim_EMT['V_Diode'].values, \"r-\", label='V_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "#plt.legend(bbox_to_anchor=(1,1))\n", + "\n", + "plt.title('Diode Voltage')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Voltage [V]')\n", + "\n", + "fig2 = plt.figure()\n", + "\n", + "plt.plot(ts_dpsim_EMT['I_Diode'].time, ts_dpsim_EMT['I_Diode'].values, \"r-\", label='I_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "#plt.legend(bbox_to_anchor=(1,1))\n", + "\n", + "plt.title('Diode Current')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Current [A]')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + }, + "tests": { + "skip": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/Notebooks/Circuits/EMT_Ph1_VS_R1_Diode.ipynb b/examples/Notebooks/Circuits/EMT_Ph1_VS_R1_Diode.ipynb new file mode 100644 index 0000000000..11b1eb6406 --- /dev/null +++ b/examples/Notebooks/Circuits/EMT_Ph1_VS_R1_Diode.ipynb @@ -0,0 +1,139 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from villas.dataprocessing.readtools import *\n", + "from villas.dataprocessing.timeseries import *\n", + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import dpsimpy\n", + "import re\n", + "import os\n", + "\n", + "%matplotlib widget\n", + "\n", + "#Define simulation scenario\n", + "time_step = 0.0001\n", + "final_time = 0.1 #0.006\n", + "simName = \"EMT_Ph1_VS_R1_Diode\"\n", + "\n", + "dpsimpy.Logger.set_log_dir('logs/' + simName)\n", + "\n", + "#Nodes\n", + "gnd = dpsimpy.emt.SimNode.gnd\n", + "n1 = dpsimpy.emt.SimNode(\"n1\", dpsimpy.PhaseType.Single)\n", + "n2 = dpsimpy.emt.SimNode(\"n2\", dpsimpy.PhaseType.Single)\n", + "\n", + "#Components\n", + "vs = dpsimpy.emt.ph1.VoltageSource(\"vs\")\n", + "vs.set_parameters(complex(1.,0.), 50.0);\n", + "\n", + "load = dpsimpy.emt.ph1.Resistor(\"r1\")\n", + "load.set_parameters(10.)\n", + "\n", + "diode = dpsimpy.emt.ph1.Diode(\"Diode\");\n", + "\n", + "#Topology\n", + "vs.connect([gnd, n1]);\n", + "load.connect([n2, n1]);\n", + "diode.connect([gnd, n2]);\n", + "\n", + "sys = dpsimpy.SystemTopology(50, [n1, n2], [vs, load, diode])\n", + "\n", + "#Logging\n", + "logger = dpsimpy.Logger(simName)\n", + "logger.log_attribute(\"I_Diode\", \"i_intf\", diode)\n", + "logger.log_attribute(\"V_Diode\", \"v_intf\", diode)\n", + "logger.log_attribute(\"V_Resistor\", \"v_intf\", load)\n", + "logger.log_attribute(\"V_VS\", \"v_intf\", vs)\n", + "\n", + "# Simulation\n", + "sim = dpsimpy.Simulation(simName)\n", + "sim.set_system(sys)\n", + "sim.set_time_step(time_step)\n", + "sim.set_final_time(final_time)\n", + "sim.set_domain(dpsimpy.Domain.EMT)\n", + "sim.set_solver(dpsimpy.Solver.ITERATIVEMNA)\n", + "sim.add_logger(logger)\n", + "sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_logfile = os.getcwd() + \"/logs/\" + simName + '/' + simName + '.csv'\n", + "ts_dpsim_EMT = read_timeseries_dpsim(path_logfile)\n", + "\n", + "plt.close('all')\n", + "fig1 = plt.figure()\n", + "\n", + "plt.plot(ts_dpsim_EMT['V_Diode'].time, ts_dpsim_EMT['V_Diode'].values, \"b-\", label='V_D')\n", + "plt.plot(ts_dpsim_EMT['V_Resistor'].time, ts_dpsim_EMT['V_Resistor'].values, \"g-\", label='V_R')\n", + "#plt.plot(ts_dpsim_EMT['V_VS'].time, ts_dpsim_EMT['V_VS'].values, \"r-\", label='V_VS')\n", + "\n", + "#difference = ts_dpsim_EMT['V_VS'].values - (ts_dpsim_EMT['V_Resistor'].values+ts_dpsim_EMT['V_Diode'].values)\n", + "#plt.plot(ts_dpsim_EMT['V_VS'].time, difference, \"b-\", label='difference')\n", + "\n", + "plt.legend(loc = 4)\n", + "#plt.legend(bbox_to_anchor=(1,1))\n", + "\n", + "plt.title('Voltages')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Voltage [V]')\n", + "\n", + "fig2 = plt.figure()\n", + "\n", + "plt.plot(ts_dpsim_EMT['I_Diode'].time, ts_dpsim_EMT['I_Diode'].values, \"r-\", label='I_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "\n", + "plt.title('Diode Current')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Current [A]')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}