diff --git a/dpsim-models/include/dpsim-models/Components.h b/dpsim-models/include/dpsim-models/Components.h index 29a4f0bd58..38a5e5fb1c 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_ExponentialDiode.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_ExponentialDiode.h new file mode 100644 index 0000000000..73408fa65a --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_ExponentialDiode.h @@ -0,0 +1,84 @@ +/* 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 ExponentialDiode using the Shockley ideal diode equation +class ExponentialDiode : public MNASimPowerComp, + public SharedFactory, + public MNANonlinearVariableCompInterface { +protected: + Matrix Jacobian = Matrix::Zero(1, 1); + +public: + //Reverse-bias saturation current. Default: mI_S = 0.000001 + const CPS::Attribute::Ptr mI_S; + //Thermal voltage. Default: mV_T = 0.027 + const CPS::Attribute::Ptr mV_T; + + /// Defines UID, name, component parameters and logging level + ExponentialDiode(String uid, String name, + Logger::Level logLevel = Logger::Level::off); + /// Defines name, component parameters and logging level + ExponentialDiode(String name, Logger::Level logLevel = Logger::Level::off) + : ExponentialDiode(name, name, logLevel) {} + + // #### General #### + /// + SimPowerComp::Ptr clone(String name) override; + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + //Sets reverse-bias saturation current I_S and thermal voltage V_T. + //If this method is not called, the diode will have the following + //values by default: I_S = 0.000001 (A) and V_T = 0.027 (V). + void setParameters(Real I_S, Real V_T); + + // #### 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 cb46ea72cd..45a235f5ed 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_ExponentialDiode.cpp EMT/EMT_Ph3_CurrentSource.cpp EMT/EMT_Ph3_VoltageSource.cpp diff --git a/dpsim-models/src/EMT/EMT_Ph1_ExponentialDiode.cpp b/dpsim-models/src/EMT/EMT_Ph1_ExponentialDiode.cpp new file mode 100644 index 0000000000..a867668ae4 --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph1_ExponentialDiode.cpp @@ -0,0 +1,178 @@ +/* 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::ExponentialDiode::ExponentialDiode(String uid, String name, + Logger::Level logLevel) + : MNASimPowerComp(uid, name, false, false, logLevel), + mI_S(mAttributes->create("I_S")), + mV_T(mAttributes->create("V_T")) { + mPhaseType = PhaseType::Single; + setTerminalNumber(2); + **mI_S = 0.000001; + **mV_T = 0.027; + **mIntfVoltage = Matrix::Zero(1, 1); + **mIntfCurrent = Matrix::Zero(1, 1); +} + +void EMT::Ph1::ExponentialDiode::setParameters(Real I_S, Real V_T) { + **mI_S = I_S; + **mV_T = V_T; +} + +SimPowerComp::Ptr EMT::Ph1::ExponentialDiode::clone(String name) { + auto copy = ExponentialDiode::make(name, mLogLevel); + copy->setParameters(**mI_S, **mV_T); + return copy; +} + +void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::mnaCompPreStep(Real time, Int timeStepCount) { + mnaCompApplyRightSideVectorStamp(**mRightVector); +} + +void EMT::Ph1::ExponentialDiode::mnaCompPostStep( + Real time, Int timeStepCount, Attribute::Ptr &leftVector) { + mnaUpdateVoltage(**leftVector); + mnaUpdateCurrent(**leftVector); +} + +void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::mnaCompUpdateCurrent( + const Matrix &leftVector) { + (**mIntfCurrent)(0, 0) = + (**mI_S) * (expf((**mIntfVoltage)(0, 0) / (**mV_T)) - 1.); +} + +void EMT::Ph1::ExponentialDiode::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::ExponentialDiode::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::ExponentialDiode::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 0edc927961..2fac7bde7c 100644 --- a/dpsim/examples/cxx/CMakeLists.txt +++ b/dpsim/examples/cxx/CMakeLists.txt @@ -35,6 +35,7 @@ set(CIRCUIT_SOURCES Circuits/EMT_PiLine.cpp Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp + Circuits/EMT_Ph1_ExponentialDiode_test.cpp # EMT examples with PF initialization Circuits/EMT_Slack_PiLine_PQLoad_with_PF_Init.cpp @@ -136,6 +137,7 @@ list(APPEND TEST_SOURCES Circuits/DP_SMIB_ReducedOrderSGIterative_LoadStep.cpp Circuits/EMT_Ph3_R3C1L1CS1_RC_vs_SSN.cpp Circuits/EMT_Ph3_RLC1VS1_RC_vs_SSN.cpp + Circuits/EMT_Ph1_ExponentialDiode_test.cpp ) if(WITH_SUNDIALS) diff --git a/dpsim/examples/cxx/Circuits/EMT_Ph1_ExponentialDiode_test.cpp b/dpsim/examples/cxx/Circuits/EMT_Ph1_ExponentialDiode_test.cpp new file mode 100644 index 0000000000..70bc6adcf9 --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_Ph1_ExponentialDiode_test.cpp @@ -0,0 +1,68 @@ +#include "../Examples.h" +#include + +using namespace DPsim; +using namespace CPS::EMT; + +void EMT_Ph1_ExponentialDiode() { + // Define simulation scenario + Real timeStep = 0.0001; + Real finalTime = 0.1; + String simName = "EMT_Ph1_ExponentialDiode_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 expDiode = Ph1::ExponentialDiode::make("ExponentialDiode"); + expDiode->setParameters( + 0.000001, 0.027); //Calling this is optional. If this method call + //is omitted, the diode will get the following + //values by default: + //I_S = 0.000001 (A) and V_T = 0.027 (V). + + // Topology + load->connect(SimNode::List{n1, n2}); + + expDiode->connect(SimNode::List{ + SimNode::GND, + n2}); //Looking at the diode's current flow direction, always + //connect by listing the output/low voltage side first + //and the input/high voltage side second + + vs0->connect(SimNode::List{SimNode::GND, n1}); + + // Define system topology + auto sys = SystemTopology(50, SystemNodeList{n1, n2}, + SystemComponentList{load, vs0, expDiode}); + + // Logging + auto logger = DataLogger::make(simName); + logger->logAttribute("I_ExponentialDiode", expDiode->attribute("i_intf")); + logger->logAttribute("V_ExponentialDiode", expDiode->attribute("v_intf")); + + Simulation sim(simName, Logger::Level::info); + sim.doInitFromNodesAndTerminals(false); + 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_ExponentialDiode(); + 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..b14a76e675 --- /dev/null +++ b/dpsim/src/IterativeMnaSolverDirect.cpp @@ -0,0 +1,196 @@ +/* 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. + + // 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 + + Matrix ResetFunctionResult = + mNonlinearFunctionResult; //use automatic adjustment of dimensions + mNonlinearFunctionResult -= + ResetFunctionResult; //to reset function result vector + + 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); + if (abs(calculationErrorElement) >= Epsilon) { + isConverged = false; + break; + } + } + iterations++; + } while (!isConverged && iterations < 100); + /// 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 edb5b06fae..ffd360e9b5 100644 --- a/dpsim/src/pybind/EMTComponents.cpp +++ b/dpsim/src/pybind/EMTComponents.cpp @@ -114,6 +114,16 @@ void addEMTPh1Components(py::module_ mEMTPh1) { .def("open", &CPS::EMT::Ph1::Switch::open) .def("close", &CPS::EMT::Ph1::Switch::close) .def("connect", &CPS::EMT::Ph1::Switch::connect); + + py::class_, + CPS::SimPowerComp>(mEMTPh1, "ExponentialDiode", + py::multiple_inheritance()) + .def(py::init()) + .def(py::init()) + .def("set_parameters", &CPS::EMT::Ph1::ExponentialDiode::setParameters, + "I_S"_a, "V_T"_a) + .def("connect", &CPS::EMT::Ph1::ExponentialDiode::connect); } void addEMTPh3Components(py::module_ mEMTPh3) { 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_ExponentialDiode_test.ipynb b/examples/Notebooks/Circuits/EMT_Ph1_ExponentialDiode_test.ipynb new file mode 100644 index 0000000000..79baf3d346 --- /dev/null +++ b/examples/Notebooks/Circuits/EMT_Ph1_ExponentialDiode_test.ipynb @@ -0,0 +1,125 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# EMT Ph1 Exponential Diode Test Circuit" + ] + }, + { + "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_ExponentialDiode_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_ExponentialDiode'].time, ts_dpsim_EMT['V_ExponentialDiode'].values, \"r-\", label='V_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "#plt.legend(bbox_to_anchor=(1,1))\n", + "\n", + "plt.title('ExponentialDiode Voltage')\n", + "plt.xlabel('t [s]')\n", + "plt.ylabel('Voltage [V]')\n", + "\n", + "fig2 = plt.figure()\n", + "\n", + "plt.plot(ts_dpsim_EMT['I_ExponentialDiode'].time, ts_dpsim_EMT['I_ExponentialDiode'].values, \"r-\", label='I_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "#plt.legend(bbox_to_anchor=(1,1))\n", + "\n", + "plt.title('ExponentialDiode 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_ExponentialDiode.ipynb b/examples/Notebooks/Circuits/EMT_Ph1_VS_R1_ExponentialDiode.ipynb new file mode 100644 index 0000000000..a3eb65383e --- /dev/null +++ b/examples/Notebooks/Circuits/EMT_Ph1_VS_R1_ExponentialDiode.ipynb @@ -0,0 +1,148 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Test of Exponential Diode: VS-R-Diode in series" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup Simulation" + ] + }, + { + "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_ExponentialDiode\"\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", + "expDiode = dpsimpy.emt.ph1.ExponentialDiode(\"ExponentialDiode\");\n", + "\n", + "#Topology\n", + "vs.connect([gnd, n1]);\n", + "load.connect([n2, n1]);\n", + "expDiode.connect([gnd, n2]);\n", + "\n", + "sys = dpsimpy.SystemTopology(50, [n1, n2], [vs, load, expDiode])\n", + "\n", + "#Logging\n", + "logger = dpsimpy.Logger(simName)\n", + "logger.log_attribute(\"I_ExponentialDiode\", \"i_intf\", expDiode)\n", + "logger.log_attribute(\"V_ExponentialDiode\", \"v_intf\", expDiode)\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": "markdown", + "metadata": {}, + "source": [ + "## Read and visualize results" + ] + }, + { + "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_ExponentialDiode'].time, ts_dpsim_EMT['V_ExponentialDiode'].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", + "\n", + "plt.legend(loc = 4)\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_ExponentialDiode'].time, ts_dpsim_EMT['I_ExponentialDiode'].values, \"r-\", label='I_D')\n", + "\n", + "plt.legend(loc = 4)\n", + "\n", + "plt.title('ExponentialDiode 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 +}