diff --git a/docs/hugo/content/en/docs/Development/Guidelines.md b/docs/hugo/content/en/docs/Development/Guidelines.md index a2d3a20e1d..ecdcf238d5 100644 --- a/docs/hugo/content/en/docs/Development/Guidelines.md +++ b/docs/hugo/content/en/docs/Development/Guidelines.md @@ -5,6 +5,20 @@ linkTitle: "Guidelines" This is a summary of general guidelines for the development of DPsim. +## Scaling of Voltages and Currents + +Voltage quantities are expressed either as phase-to-phase RMS values (denominated as `RMS3PH`) or as phase-to-ground peak values (denominated as `PEAK1PH`): + +- Initialisation quantities (e.g. `initialSingleVoltage` of `SimPowerComp`) as `RMS3PH` values +- Simulation quantities in both `SP` and `DP` domain (e.g. `mIntfVoltage` of `DP::Ph1::PiLine`) as `RMS3PH values` +- Simulation quantities in the `EMT` domain (e.g. `mIntfVoltage` of `EMT::Ph3::Transformer`) as `PEAK1PH` values + +Current quantities are expressed either as `RMS` or as `PEAK` values: + +- Simulation quantities in both `SP` and `DP` domain (e.g. `mIntfCurrent` of `DP::Ph1::PiLine`) as `RMS` values +- Simulation quantities in the `EMT` domain (e.g. `mIntfCurrent` of `EMT::Ph3::Transformer`) as `PEAK` values + + ## Logging Debug or trace should be the default log level for information that might be nice to have but not necessary for every simulation case. @@ -23,4 +37,4 @@ A new version of DPsim has to be indicated as follows: - Update sonar-project.properties Due to the creation of a new tag, a new PyPi package will be deployed automatically. -To release an updated Docker image, the container workflow needs to be triggered manually. \ No newline at end of file +To release an updated Docker image, the container workflow needs to be triggered manually. diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h index 6e9a990f17..0f663c5f40 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph1_VoltageSource.h @@ -25,6 +25,8 @@ namespace Ph1 { class VoltageSource : public MNASimPowerComp, public SharedFactory { + private: + Real mTimeStep; protected: void updateVoltage(Real time); public: diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorVBR.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorVBR.h index 0295639477..855273c8e6 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorVBR.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorVBR.h @@ -241,12 +241,13 @@ namespace Ph3 { void CalculateAuxiliarConstants(Real dt); void CalculateAuxiliarVariables(); + /// Getters //Matrix& rotorFluxes() { return mRotorFlux; } - Matrix& dqStatorCurrents(); - Real electricalTorque() const; - Real rotationalSpeed() const; - Real rotorPosition() const; - Matrix& statorCurrents(); + Matrix& dqStatorCurrents() { return mDqStatorCurrents; } + Real electricalTorque() const { return **mElecTorque * mBase_T; } + Real rotationalSpeed() const { return **mOmMech * mBase_OmMech; } + Real rotorPosition() const { return mThetaMech; } + Matrix& statorCurrents() { return mIabc; } // #### MNA section #### /// Stamps system matrix @@ -262,7 +263,7 @@ namespace Ph3 { /// Add MNA post step dependencies void mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector); /// Mark that parameter changes so that system matrix is updated - Bool hasParameterChanged() override; + Bool hasParameterChanged() override { return true; } }; } } diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_VoltageSource.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_VoltageSource.h index b5bae13e3a..76e903941d 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_VoltageSource.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_VoltageSource.h @@ -31,6 +31,8 @@ namespace CPS { private: /// CPS::Signal::SignalGenerator::Ptr mSrcSig; + /// + Real mTimeStep; protected: // Updates voltage according to reference phasor and frequency void updateVoltage(Real time); diff --git a/dpsim-models/include/dpsim-models/SimNode.h b/dpsim-models/include/dpsim-models/SimNode.h index edacdeca03..bfcbae2072 100644 --- a/dpsim-models/include/dpsim-models/SimNode.h +++ b/dpsim-models/include/dpsim-models/SimNode.h @@ -58,6 +58,8 @@ namespace CPS { SimNode(UInt matrixNodeIndex, PhaseType phaseType = PhaseType::Single) : SimNode("N" + std::to_string(matrixNodeIndex), "N" + std::to_string(matrixNodeIndex), matrixNodeIndex, phaseType) { } + /// Initialize mVoltage according to mInitialVoltage + void initialize(); /// Initialize state matrices with size according to phase type and frequency number void initialize(Matrix frequencies); /// Returns matrix index for specified phase diff --git a/dpsim-models/src/Base/Base_ReducedOrderSynchronGenerator.cpp b/dpsim-models/src/Base/Base_ReducedOrderSynchronGenerator.cpp index 507366797b..87f2f2b716 100644 --- a/dpsim-models/src/Base/Base_ReducedOrderSynchronGenerator.cpp +++ b/dpsim-models/src/Base/Base_ReducedOrderSynchronGenerator.cpp @@ -349,6 +349,16 @@ void Base::ReducedOrderSynchronGenerator::initializeFromNodesAndTerminals( // initialize theta and calculate transform matrix **mThetaMech = **mDelta - PI / 2.; + // set initial interface current + (**mIntfCurrent)(0,0) = (mInitCurrent * mBase_I).real(); + (**mIntfCurrent)(1,0) = (mInitCurrent * mBase_I * SHIFT_TO_PHASE_B).real(); + (**mIntfCurrent)(2,0) = (mInitCurrent * mBase_I * SHIFT_TO_PHASE_C).real(); + + // set initial interface voltage + (**mIntfVoltage)(0,0) = (mInitVoltage * mBase_V).real(); + (**mIntfVoltage)(1,0) = (mInitVoltage * mBase_V * SHIFT_TO_PHASE_B).real(); + (**mIntfVoltage)(2,0) = (mInitVoltage * mBase_V * SHIFT_TO_PHASE_C).real(); + SPDLOG_LOGGER_DEBUG(this->mSLog, "\n--- Initialization from power flow ---" "\nInitial Vd (per unit): {:f}" @@ -423,6 +433,12 @@ void Base::ReducedOrderSynchronGenerator::initializeFromNodesAndTermina // initialize theta and calculate transform matrix **mThetaMech = **mDelta - PI / 2.; + // set initial value of current + (**mIntfCurrent)(0,0) = mInitCurrent * mBase_I_RMS; + + // set initial interface voltage + (**mIntfVoltage)(0,0) = mInitVoltage * mBase_V_RMS; + SPDLOG_LOGGER_DEBUG(this->mSLog, "\n--- Initialization from power flow ---" "\nInitial Vd (per unit): {:f}" @@ -476,16 +492,11 @@ void Base::ReducedOrderSynchronGenerator::mnaCompPreStep(Real time, Int **mMechTorque = mTurbineGovernor->step(**mOmMech, mTimeStep); } - // predict mechanical vars for all reduced-order models in the same manner - if (mSimTime > 0.0) { - // predict omega at t=k+1 (forward euler) - **mElecTorque = (**mVdq)(0,0) * (**mIdq)(0,0) + (**mVdq)(1,0) * (**mIdq)(1,0); - **mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque)); - - // predict theta and delta at t=k+1 (backward euler) - **mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech); - **mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech; - } + // calculate mechanical variables at t=k+1 with forward euler + **mElecTorque = (**mVdq)(0,0) * (**mIdq)(0,0) + (**mVdq)(1,0) * (**mIdq)(1,0); + **mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque)); + **mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech); + **mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech; // model specific calculation of electrical vars stepInPerUnit(); @@ -509,16 +520,11 @@ void Base::ReducedOrderSynchronGenerator::mnaCompPreStep(Real time, Int ti **mMechTorque = mTurbineGovernor->step(**mOmMech, mTimeStep); } - // predict mechanical vars for all reduced-order models in the same manner - if (mSimTime > 0.0) { - // predict omega at t=k+1 (forward euler) - **mElecTorque = (**mVdq0)(0,0) * (**mIdq0)(0,0) + (**mVdq0)(1,0) * (**mIdq0)(1,0); - **mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque)); - - // predict theta and delta at t=k+1 (backward euler) - **mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech); - **mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech; - } + // calculate mechanical variables at t=k+1 with forward euler + **mElecTorque = ((**mVdq0)(0,0) * (**mIdq0)(0,0) + (**mVdq0)(1,0) * (**mIdq0)(1,0)); + **mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque)); + **mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech); + **mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech; // model specific calculation of electrical vars stepInPerUnit(); diff --git a/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp b/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp index aa56bfa4a6..0a37d61dcc 100644 --- a/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp +++ b/dpsim-models/src/DP/DP_Ph1_Capacitor.cpp @@ -55,6 +55,7 @@ void DP::Ph1::Capacitor::initializeFromNodesAndTerminals(Real frequency) { Logger::phasorToString((**mIntfCurrent)(0,0)), Logger::phasorToString(initialSingleVoltage(0)), Logger::phasorToString(initialSingleVoltage(1))); + mSLog->flush(); } void DP::Ph1::Capacitor::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { @@ -82,6 +83,7 @@ void DP::Ph1::Capacitor::mnaCompInitialize(Real omega, Real timeStep, Attribute< Logger::phasorToString((**mIntfVoltage)(0,0)), Logger::phasorToString((**mIntfCurrent)(0,0)), Logger::complexToString(mEquivCurrent(0,0))); + mSLog->flush(); } void DP::Ph1::Capacitor::mnaCompInitializeHarm(Real omega, Real timeStep, std::vector::Ptr> leftVectors) { diff --git a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator4OrderTPM.cpp b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator4OrderTPM.cpp index f98181aa3a..16705622b8 100644 --- a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator4OrderTPM.cpp +++ b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator4OrderTPM.cpp @@ -74,6 +74,11 @@ void DP::Ph1::SynchronGenerator4OrderTPM::specificInitialization() { // initial emf in the dq reference frame (**mEdq_t)(0,0) = (**mVdq)(0,0) - (**mIdq)(1,0) * mLq_t; (**mEdq_t)(1,0) = (**mVdq)(1,0) + (**mIdq)(0,0) * mLd_t; + + // set previous values of stator current at simulation start + //(**mIntfCurrent)(0,0) = std::conj(mInitElecPower / (mInitVoltage * mBase_V_RMS)); + mIdpTwoPrevStep = **mIntfCurrent; + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nInitial Ed_t (per unit): {:f}" @@ -147,13 +152,6 @@ void DP::Ph1::SynchronGenerator4OrderTPM::stepInPerUnit() { mResistanceMatrixVarying(1,1) = (mA + mB) / 2.0 * sin(2*DeltaTheta); SPDLOG_LOGGER_DEBUG(mSLog, "\nR_var [pu] (t={:f}): {:s}", mSimTime, Logger::matrixToString(mResistanceMatrixVarying)); - // predict electrical vars - // set previous values of stator current at simulation start - if (mSimTime == 0.0) { - (**mIntfCurrent)(0,0) = std::conj(mInitElecPower / (mInitVoltage * mBase_V_RMS)); - mIdpTwoPrevStep = **mIntfCurrent; - } - // predict stator current (linear extrapolation) Matrix IdpPrediction = Matrix::Zero(2,1); IdpPrediction(0,0) = 2 * (**mIntfCurrent)(0,0).real() - mIdpTwoPrevStep(0,0).real(); diff --git a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6aOrderVBR.cpp b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6aOrderVBR.cpp index 432eace345..e01fd9b03d 100644 --- a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6aOrderVBR.cpp +++ b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6aOrderVBR.cpp @@ -41,6 +41,10 @@ void DP::Ph1::SynchronGenerator6aOrderVBR::specificInitialization() { (**mEdq_s)(0,0) = (**mVdq)(0,0) - mLq_s * (**mIdq)(1,0); (**mEdq_s)(1,0) = (**mVdq)(1,0) + mLd_s * (**mIdq)(0,0); + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq)(1,0) + mBd_t * (**mEdq_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq)(0,0) + mBq_t * (**mEdq_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type a (Marconato's model)" @@ -63,15 +67,13 @@ void DP::Ph1::SynchronGenerator6aOrderVBR::stepInPerUnit() { mDomainInterface.updateDQToDPTransform(**mThetaMech, mSimTime); mDomainInterface.updateDPToDQTransform(**mThetaMech, mSimTime); - if (mSimTime>0.0){ - // calculate Edq_t at t=k - (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); - (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); + // calculate Edq_t at t=k + (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); + (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); - // calculate Edq_s at t=k - (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); - (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); - } + // calculate Edq_s at t=k + (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); + (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); // Update time-varying reactance matrix calculateConductanceMatrix(); diff --git a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6bOrderVBR.cpp b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6bOrderVBR.cpp index fdabacaeee..77d6655990 100644 --- a/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6bOrderVBR.cpp +++ b/dpsim-models/src/DP/DP_Ph1_SynchronGenerator6bOrderVBR.cpp @@ -41,7 +41,11 @@ void DP::Ph1::SynchronGenerator6bOrderVBR::specificInitialization() { (**mEdq_s)(0,0) = (**mVdq)(0,0) - mLq_s * (**mIdq)(1,0); (**mEdq_s)(1,0) = (**mVdq)(1,0) + mLd_s * (**mIdq)(0,0); - SPDLOG_LOGGER_INFO(mSLog, + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq)(1,0) + mBd_t * (**mEdq_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq)(0,0) + mBq_t * (**mEdq_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type b (Anderson - Fouad's model)" "\nInitial Ed_t (per unit): {:f}" @@ -64,15 +68,13 @@ void DP::Ph1::SynchronGenerator6bOrderVBR::stepInPerUnit() { mDomainInterface.updateDQToDPTransform(**mThetaMech, mSimTime); mDomainInterface.updateDPToDQTransform(**mThetaMech, mSimTime); - if (mSimTime>0.0){ - // calculate Edq_t at t=k - (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); - (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); + // calculate Edq_t at t=k + (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); + (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); - // calculate Edq_s at t=k - (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); - (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); - } + // calculate Edq_s at t=k + (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); + (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); // Update time-varying reactance matrix calculateConductanceMatrix(); diff --git a/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp index af095071b6..3cb33fdbdd 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Capacitor.cpp @@ -28,8 +28,9 @@ void EMT::Ph1::Capacitor::initializeFromNodesAndTerminals(Real frequency) { Real omega = 2 * PI * frequency; Complex impedance = { 0, - 1. / (omega * **mCapacitance) }; - (**mIntfVoltage)(0,0) = (initialSingleVoltage(1) - initialSingleVoltage(0)).real(); - (**mIntfCurrent)(0,0) = ((initialSingleVoltage(1) - initialSingleVoltage(0)) / impedance).real(); + Complex voltage = RMS3PH_TO_PEAK1PH * (initialSingleVoltage(1) - initialSingleVoltage(0)); + (**mIntfVoltage)(0,0) = voltage.real(); + (**mIntfCurrent)(0,0) = (voltage / impedance).real(); SPDLOG_LOGGER_INFO(mSLog, "\n--- Initialization from powerflow ---" @@ -40,8 +41,9 @@ void EMT::Ph1::Capacitor::initializeFromNodesAndTerminals(Real frequency) { "\n--- Initialization from powerflow finished ---", (**mIntfVoltage)(0,0), (**mIntfCurrent)(0,0), - initialSingleVoltage(0).real(), - initialSingleVoltage(1).real()); + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)).real(), + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(1)).real()); + mSLog->flush(); } void EMT::Ph1::Capacitor::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { diff --git a/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp index f10766a7d6..efa8b73557 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Inductor.cpp @@ -28,8 +28,9 @@ void EMT::Ph1::Inductor::initializeFromNodesAndTerminals(Real frequency) { Real omega = 2 * PI * frequency; Complex impedance = { 0, omega * **mInductance }; - (**mIntfVoltage)(0,0) = (initialSingleVoltage(1) - initialSingleVoltage(0)).real(); - (**mIntfCurrent)(0,0) = ((initialSingleVoltage(1) - initialSingleVoltage(0)) / impedance).real(); + Complex voltage = RMS3PH_TO_PEAK1PH * (initialSingleVoltage(1) - initialSingleVoltage(0)); + (**mIntfVoltage)(0,0) = voltage.real(); + (**mIntfCurrent)(0,0) = (voltage / impedance).real(); SPDLOG_LOGGER_INFO(mSLog, "\n--- Initialization from powerflow ---" @@ -40,8 +41,8 @@ void EMT::Ph1::Inductor::initializeFromNodesAndTerminals(Real frequency) { "\n--- Initialization from powerflow finished ---", (**mIntfVoltage)(0,0), (**mIntfCurrent)(0,0), - initialSingleVoltage(0).real(), - initialSingleVoltage(1).real()); + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)).real(), + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(1)).real()); } void EMT::Ph1::Inductor::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { diff --git a/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp b/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp index b7b2456201..8fb45d6ad6 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_Resistor.cpp @@ -24,8 +24,8 @@ SimPowerComp::Ptr EMT::Ph1::Resistor::clone(String name) { } void EMT::Ph1::Resistor::initializeFromNodesAndTerminals(Real frequency) { - - (**mIntfVoltage)(0,0) = (initialSingleVoltage(1) - initialSingleVoltage(0)).real(); + Complex voltage = RMS3PH_TO_PEAK1PH * (initialSingleVoltage(1) - initialSingleVoltage(0)); + (**mIntfVoltage)(0,0) = voltage.real(); (**mIntfCurrent)(0,0) = (**mIntfVoltage)(0,0) / **mResistance; SPDLOG_LOGGER_INFO(mSLog, @@ -37,8 +37,9 @@ void EMT::Ph1::Resistor::initializeFromNodesAndTerminals(Real frequency) { "\n--- Initialization from powerflow finished ---", (**mIntfVoltage)(0,0), (**mIntfCurrent)(0,0), - initialSingleVoltage(0).real(), - initialSingleVoltage(1).real()); + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(0)).real(), + (RMS3PH_TO_PEAK1PH * initialSingleVoltage(1)).real()); + mSLog->flush(); } void EMT::Ph1::Resistor::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { @@ -67,6 +68,7 @@ void EMT::Ph1::Resistor::mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMa SPDLOG_LOGGER_INFO(mSLog, "Add {:f} to system at ({:d},{:d})", -conductance, matrixNodeIndex(0), matrixNodeIndex(1)); SPDLOG_LOGGER_INFO(mSLog, "Add {:f} to system at ({:d},{:d})", -conductance, matrixNodeIndex(1), matrixNodeIndex(0)); } + mSLog->flush(); } void EMT::Ph1::Resistor::mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) { diff --git a/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp b/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp index 617e3cf7a6..cc4afb5cba 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_VoltageSource.cpp @@ -34,8 +34,10 @@ SimPowerComp::Ptr EMT::Ph1::VoltageSource::clone(String name) { } void EMT::Ph1::VoltageSource::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { - updateMatrixNodeIndices(); + updateMatrixNodeIndices(); (**mIntfVoltage)(0,0) = Math::abs(**mVoltageRef) * cos(Math::phase(**mVoltageRef)); + + mTimeStep = timeStep; } void EMT::Ph1::VoltageSource::mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMatrix) { @@ -66,7 +68,7 @@ void EMT::Ph1::VoltageSource::updateVoltage(Real time) { Complex voltageRef = mVoltageRef->get(); Real srcFreq = mSrcFreq->get(); if (srcFreq > 0) - (**mIntfVoltage)(0,0) = Math::abs(voltageRef) * cos(time * 2.*PI*srcFreq + Math::phase(voltageRef)); + (**mIntfVoltage)(0,0) = Math::abs(voltageRef) * cos((time) * 2.*PI*srcFreq + Math::phase(voltageRef)); else (**mIntfVoltage)(0,0) = voltageRef.real(); } diff --git a/dpsim-models/src/EMT/EMT_Ph1_VoltageSourceRamp.cpp b/dpsim-models/src/EMT/EMT_Ph1_VoltageSourceRamp.cpp index 90b8625ead..e6675fa073 100644 --- a/dpsim-models/src/EMT/EMT_Ph1_VoltageSourceRamp.cpp +++ b/dpsim-models/src/EMT/EMT_Ph1_VoltageSourceRamp.cpp @@ -43,8 +43,8 @@ void EMT::Ph1::VoltageSourceRamp::setParameters(Complex voltage, Complex addVolt void EMT::Ph1::VoltageSourceRamp::initialize(Matrix frequencies) { SimPowerComp::initialize(frequencies); - if (**mVoltageRef == Complex(0, 0)) - **mVoltageRef = initialSingleVoltage(1) - initialSingleVoltage(0); + if (**mVoltageRef == Complex(0, 0)) + **mVoltageRef = RMS3PH_TO_PEAK1PH * (initialSingleVoltage(1) - initialSingleVoltage(0)); mSubVoltageSource = VoltageSource::make(**mName + "_src", mLogLevel); mSubVoltageSource->setParameters(**mVoltageRef, 0); diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6aOrderVBR.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6aOrderVBR.cpp index d27774b519..8351dfc752 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6aOrderVBR.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6aOrderVBR.cpp @@ -32,6 +32,7 @@ EMT::Ph3::SynchronGenerator6aOrderVBR::SynchronGenerator6aOrderVBR } void EMT::Ph3::SynchronGenerator6aOrderVBR::specificInitialization() { + // initial voltage behind the transient reactance in the dq reference frame (**mEdq0_t)(0,0) = (mLq - mLq_t - mYq) * (**mIdq0)(1,0); (**mEdq0_t)(1,0) = (1 - mTaa / mTd0_t) * **mEf - (mLd - mLd_t - mYd) * (**mIdq0)(0,0); @@ -40,6 +41,11 @@ void EMT::Ph3::SynchronGenerator6aOrderVBR::specificInitialization() { (**mEdq0_s)(0,0) = (**mVdq0)(0,0) - mLq_s * (**mIdq0)(1,0); (**mEdq0_s)(1,0) = (**mVdq0)(1,0) + mLd_s * (**mIdq0)(0,0); + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq0)(1,0) + mBd_t * (**mEdq0_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq0)(0,0) + mBq_t * (**mEdq0_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + mEh_t(2,0) = 0.0; + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type a (Marconato's model)" @@ -59,16 +65,14 @@ void EMT::Ph3::SynchronGenerator6aOrderVBR::specificInitialization() { void EMT::Ph3::SynchronGenerator6aOrderVBR::stepInPerUnit() { - if (mSimTime>0.0) { - // calculate Edq_t at t=k - (**mEdq0_t)(0,0) = mAd_t * (**mIdq0)(1,0) + mEh_t(0,0); - (**mEdq0_t)(1,0) = mAq_t * (**mIdq0)(0,0) + mEh_t(1,0); - (**mEdq0_t)(2,0) = 0.0; + // calculate Edq_t at t=k + (**mEdq0_t)(0,0) = mAd_t * (**mIdq0)(1,0) + mEh_t(0,0); + (**mEdq0_t)(1,0) = mAq_t * (**mIdq0)(0,0) + mEh_t(1,0); + (**mEdq0_t)(2,0) = 0.0; - // calculate Edq_s at t=k - (**mEdq0_s)(0,0) = -(**mIdq0)(1,0) * mLq_s + (**mVdq0)(0,0); - (**mEdq0_s)(1,0) = (**mIdq0)(0,0) * mLd_s + (**mVdq0)(1,0); - } + // calculate Edq_s at t=k + (**mEdq0_s)(0,0) = -(**mIdq0)(1,0) * mLq_s + (**mVdq0)(0,0); + (**mEdq0_s)(1,0) = (**mIdq0)(0,0) * mLd_s + (**mVdq0)(1,0); // get transformation matrix mAbcToDq0 = get_parkTransformMatrix(); diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6bOrderVBR.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6bOrderVBR.cpp index 1f31cc2e80..4183718fc2 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6bOrderVBR.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGenerator6bOrderVBR.cpp @@ -32,6 +32,7 @@ EMT::Ph3::SynchronGenerator6bOrderVBR::SynchronGenerator6bOrderVBR } void EMT::Ph3::SynchronGenerator6bOrderVBR::specificInitialization() { + // initial voltage behind the transient reactance in the dq reference frame (**mEdq0_t)(0,0) = (mLq - mLq_t) * (**mIdq0)(1,0); (**mEdq0_t)(1,0) = **mEf - (mLd - mLd_t) * (**mIdq0)(0,0); @@ -40,7 +41,12 @@ void EMT::Ph3::SynchronGenerator6bOrderVBR::specificInitialization() { (**mEdq0_s)(0,0) = (**mVdq0)(0,0) - mLq_s * (**mIdq0)(1,0); (**mEdq0_s)(1,0) = (**mVdq0)(1,0) + mLd_s * (**mIdq0)(0,0); - SPDLOG_LOGGER_INFO(mSLog, + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq0)(1,0) + mBd_t * (**mEdq0_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq0)(0,0) + mBq_t * (**mEdq0_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + mEh_t(2,0) = 0.0; + + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type b (Anderson - Fouad's model)" "\nInitial Ed_t (per unit): {:f}" @@ -59,16 +65,14 @@ void EMT::Ph3::SynchronGenerator6bOrderVBR::specificInitialization() { void EMT::Ph3::SynchronGenerator6bOrderVBR::stepInPerUnit() { - if (mSimTime>0.0) { - // calculate Edq_t at t=k - (**mEdq0_t)(0,0) = mAd_t * (**mIdq0)(1,0) + mEh_t(0,0); - (**mEdq0_t)(1,0) = mAq_t * (**mIdq0)(0,0) + mEh_t(1,0); - (**mEdq0_t)(2,0) = 0.0; + // calculate Edq_t at t=k + (**mEdq0_t)(0,0) = mAd_t * (**mIdq0)(1,0) + mEh_t(0,0); + (**mEdq0_t)(1,0) = mAq_t * (**mIdq0)(0,0) + mEh_t(1,0); + (**mEdq0_t)(2,0) = 0.0; - // calculate Edq_s at t=k - (**mEdq0_s)(0,0) = -(**mIdq0)(1,0) * mLq_s + (**mVdq0)(0,0); - (**mEdq0_s)(1,0) = (**mIdq0)(0,0) * mLd_s + (**mVdq0)(1,0); - } + // calculate Edq_s at t=k + (**mEdq0_s)(0,0) = -(**mIdq0)(1,0) * mLq_s + (**mVdq0)(0,0); + (**mEdq0_s)(1,0) = (**mIdq0)(0,0) * mLd_s + (**mVdq0)(1,0); // get transformation matrix mAbcToDq0 = get_parkTransformMatrix(); diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQ.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQ.cpp index 43e93bb09a..f7315bfdab 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQ.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQ.cpp @@ -147,11 +147,10 @@ void EMT::Ph3::SynchronGeneratorDQ::initializeFromNodesAndTerminals(Real frequen "\n--- Initialization from powerflow finished ---", Logger::phasorToString(initialSingleVoltage(0)), Logger::complexToString(terminal(0)->singlePower())); - mSLog->flush(); } else { SPDLOG_LOGGER_INFO(mSLog, "Initial values already set, skipping initializeFromNodesAndTerminals."); - mSLog->flush(); } + mSLog->flush(); } void EMT::Ph3::SynchronGeneratorDQ::initialize(Matrix frequencies) { diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQTrapez.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQTrapez.cpp index e24470765e..15a94e8a36 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQTrapez.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorDQTrapez.cpp @@ -24,7 +24,7 @@ EMT::Ph3::SynchronGeneratorDQTrapez::SynchronGeneratorDQTrapez(String name, Logg } void EMT::Ph3::SynchronGeneratorDQTrapez::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { - updateMatrixNodeIndices(); + updateMatrixNodeIndices(); mTimeStep = timeStep; SynchronGeneratorDQ::initializeMatrixAndStates(); diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorVBR.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorVBR.cpp index bcd1b11466..98684775b7 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorVBR.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorVBR.cpp @@ -29,26 +29,15 @@ void EMT::Ph3::SynchronGeneratorVBR::initialize(Matrix frequencies) { SimPowerComp::initialize(frequencies); } -Matrix& EMT::Ph3::SynchronGeneratorVBR::dqStatorCurrents() { return mDqStatorCurrents; } - -Real EMT::Ph3::SynchronGeneratorVBR::electricalTorque() const { return **mElecTorque * mBase_T; } - -Real EMT::Ph3::SynchronGeneratorVBR::rotationalSpeed() const { return **mOmMech * mBase_OmMech; } - -Real EMT::Ph3::SynchronGeneratorVBR::rotorPosition() const { return mThetaMech; } - -Matrix& EMT::Ph3::SynchronGeneratorVBR::statorCurrents() { return mIabc; } - -Bool EMT::Ph3::SynchronGeneratorVBR::hasParameterChanged() { return true; } - void EMT::Ph3::SynchronGeneratorVBR::setBaseAndOperationalPerUnitParameters( Real nomPower, Real nomVolt, Real nomFreq, Int poleNumber, Real nomFieldCur, Real Rs, Real Ld, Real Lq, Real Ld_t, Real Lq_t, Real Ld_s, Real Lq_s, Real Ll, Real Td0_t, Real Tq0_t, Real Td0_s, Real Tq0_s, Real inertia) { - Base::SynchronGenerator::setBaseAndOperationalPerUnitParameters(nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur, - Rs, Ld, Lq, Ld_t, Lq_t, Ld_s, Lq_s, - Ll, Td0_t, Tq0_t, Td0_s, Tq0_s, inertia); + Base::SynchronGenerator::setBaseAndOperationalPerUnitParameters( + nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur, + Rs, Ld, Lq, Ld_t, Lq_t, Ld_s, Lq_s, + Ll, Td0_t, Tq0_t, Td0_s, Tq0_s, inertia); SPDLOG_LOGGER_INFO(mSLog, "Set base parameters: \n" "nomPower: {:e}\nnomVolt: {:e}\nnomFreq: {:e}\n nomFieldCur: {:e}\n", @@ -79,12 +68,13 @@ void EMT::Ph3::SynchronGeneratorVBR::setBaseAndFundamentalPerUnitParameters( nomPower, nomVolt, nomFreq, nomFieldCur, poleNumber, Rs, Ll, Lmd, Lmq, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, inertia); - SPDLOG_LOGGER_INFO(mSLog, "Set base and fundamental parameters in per unit: \n" - "nomPower: {:e}\nnomVolt: {:e}\nnomFreq: {:e}\npoleNumber: {:d}\nnomFieldCur: {:e}\n" - "Rs: {:e}\nLl: {:e}\nLmd: {:e}\nLmq: {:e}\nRfd: {:e}\nLlfd: {:e}\nRkd: {:e}\n" - "Llkd: {:e}\nRkq1: {:e}\nLlkq1: {:e}\nRkq2: {:e}\nLlkq2: {:e}\ninertia: {:e}", - nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur, - Rs, Ll, Lmd, Lmq, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, inertia); + SPDLOG_LOGGER_INFO(mSLog, + "Set base and fundamental parameters in per unit: \n" + "nomPower: {:e}\nnomVolt: {:e}\nnomFreq: {:e}\npoleNumber: {:d}\nnomFieldCur: {:e}\n" + "Rs: {:e}\nLl: {:e}\nLmd: {:e}\nLmq: {:e}\nRfd: {:e}\nLlfd: {:e}\nRkd: {:e}\n" + "Llkd: {:e}\nRkq1: {:e}\nLlkq1: {:e}\nRkq2: {:e}\nLlkq2: {:e}\ninertia: {:e}", + nomPower, nomVolt, nomFreq, poleNumber, nomFieldCur, + Rs, Ll, Lmd, Lmq, Rfd, Llfd, Rkd, Llkd, Rkq1, Llkq1, Rkq2, Llkq2, inertia); } void EMT::Ph3::SynchronGeneratorVBR::setInitialValues(Real initActivePower, Real initReactivePower, @@ -93,11 +83,12 @@ void EMT::Ph3::SynchronGeneratorVBR::setInitialValues(Real initActivePower, Real Base::SynchronGenerator::setInitialValues(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, initMechPower); - SPDLOG_LOGGER_INFO(mSLog, "Set initial values: \n" - "initActivePower: {:e}\ninitReactivePower: {:e}\ninitTerminalVolt: {:e}\n" - "initVoltAngle: {:e} \ninitMechPower: {:e}", - initActivePower, initReactivePower, initTerminalVolt, - initVoltAngle, initMechPower); + SPDLOG_LOGGER_INFO(mSLog, + "Set initial values: \n" + "initActivePower: {:e}\ninitReactivePower: {:e}\ninitTerminalVolt: {:e}\n" + "initVoltAngle: {:e} \ninitMechPower: {:e}", + initActivePower, initReactivePower, initTerminalVolt, + initVoltAngle, initMechPower); } void EMT::Ph3::SynchronGeneratorVBR::initializeFromNodesAndTerminals(Real frequency) { @@ -118,11 +109,10 @@ void EMT::Ph3::SynchronGeneratorVBR::initializeFromNodesAndTerminals(Real freque "\n--- Initialization from powerflow finished ---", Logger::phasorToString(initialSingleVoltage(0)), Logger::complexToString(terminal(0)->singlePower())); - mSLog->flush(); } else { SPDLOG_LOGGER_INFO(mSLog, "Initial values already set, skipping initializeFromNodesAndTerminals."); - mSLog->flush(); } + mSLog->flush(); } void EMT::Ph3::SynchronGeneratorVBR::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { @@ -136,7 +126,6 @@ void EMT::Ph3::SynchronGeneratorVBR::mnaCompInitialize(Real omega, Real timeStep for (auto indexPair : mVariableSystemMatrixEntries) SPDLOG_LOGGER_INFO(mSLog, "({}, {})", indexPair.first, indexPair.second); - mSystemOmega = omega; mTimeStep = timeStep; @@ -241,6 +230,14 @@ void EMT::Ph3::SynchronGeneratorVBR::mnaCompInitialize(Real omega, Real timeStep SPDLOG_LOGGER_INFO(mSLog, "Initialize right side vector of size {}", leftVector->get().rows()); SPDLOG_LOGGER_INFO(mSLog, "Component affects right side vector entries {}, {} and {}", matrixNodeIndex(0,0), matrixNodeIndex(0,1), matrixNodeIndex(0,2)); + + // set initial interface current + mVabc << mVa, mVb, mVc; + **mIntfVoltage = mVabc * mBase_V; + + // set initial interface voltage + mIabc << mIa, mIb, mIc; + **mIntfCurrent = mIabc * mBase_I; } void EMT::Ph3::SynchronGeneratorVBR::mnaCompPreStep(Real time, Int timeStepCount) { @@ -336,10 +333,7 @@ void EMT::Ph3::SynchronGeneratorVBR::mnaCompPostStep(Real time, Int timeStepCoun } // ################ Update machine stator and rotor variables ############################ - mVabc << - mVa, - mVb, - mVc; + mVabc << mVa, mVb, mVc; mVq = parkTransform(mThetaMech, mVa, mVb, mVc)(0); mVd = parkTransform(mThetaMech, mVa, mVb, mVc)(1); diff --git a/dpsim-models/src/EMT/EMT_Ph3_Transformer.cpp b/dpsim-models/src/EMT/EMT_Ph3_Transformer.cpp index 6aa249bd01..686eaedc52 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_Transformer.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_Transformer.cpp @@ -79,7 +79,7 @@ void EMT::Ph3::Transformer::initializeFromNodesAndTerminals(Real frequency) { SPDLOG_LOGGER_INFO(mSLog, "Reactance (referred to higher voltage side) = {} [Ohm]", Logger::matrixToString(omega * mInductance)); MatrixComp vInitABC = MatrixComp::Zero(3, 1); - vInitABC(0, 0) = mVirtualNodes[0]->initialSingleVoltage() - RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); + vInitABC(0, 0) = RMS3PH_TO_PEAK1PH * (mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(0)); vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B; vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C; diff --git a/dpsim-models/src/EMT/EMT_Ph3_VoltageSource.cpp b/dpsim-models/src/EMT/EMT_Ph3_VoltageSource.cpp index a860e418aa..1ad0c938cc 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_VoltageSource.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_VoltageSource.cpp @@ -99,6 +99,8 @@ SimPowerComp::Ptr EMT::Ph3::VoltageSource::clone(String name) { void EMT::Ph3::VoltageSource::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { updateMatrixNodeIndices(); + + mTimeStep = timeStep; } void EMT::Ph3::VoltageSource::mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMatrix) { diff --git a/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6aOrderVBR.cpp b/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6aOrderVBR.cpp index 79000c1812..6f01a9ff61 100644 --- a/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6aOrderVBR.cpp +++ b/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6aOrderVBR.cpp @@ -41,6 +41,10 @@ void SP::Ph1::SynchronGenerator6aOrderVBR::specificInitialization() { (**mEdq_s)(0,0) = (**mVdq)(0,0) - mLq_s * (**mIdq)(1,0); (**mEdq_s)(1,0) = (**mVdq)(1,0) + mLd_s * (**mIdq)(0,0); + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq)(1,0) + mBd_t * (**mEdq_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq)(0,0) + mBq_t * (**mEdq_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type a (Marconato's model)" @@ -59,15 +63,14 @@ void SP::Ph1::SynchronGenerator6aOrderVBR::specificInitialization() { } void SP::Ph1::SynchronGenerator6aOrderVBR::stepInPerUnit() { - if (mSimTime>0.0) { - // calculate Edq_t at t=k - (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); - (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); - - // calculate Edq_s at t=k - (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); - (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); - } + + // calculate Edq_t at t=k + (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); + (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); + + // calculate Edq_s at t=k + (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); + (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); mDqToComplexA = get_DqToComplexATransformMatrix(); mComplexAToDq = mDqToComplexA.transpose(); diff --git a/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6bOrderVBR.cpp b/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6bOrderVBR.cpp index 3c724ee413..5d53bb0ca9 100644 --- a/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6bOrderVBR.cpp +++ b/dpsim-models/src/SP/SP_Ph1_SynchronGenerator6bOrderVBR.cpp @@ -41,6 +41,10 @@ void SP::Ph1::SynchronGenerator6bOrderVBR::specificInitialization() { (**mEdq_s)(0,0) = (**mVdq)(0,0) - mLq_s * (**mIdq)(1,0); (**mEdq_s)(1,0) = (**mVdq)(1,0) + mLd_s * (**mIdq)(0,0); + // initialize history term behind the transient reactance + mEh_t(0,0) = mAd_t * (**mIdq)(1,0) + mBd_t * (**mEdq_t)(0,0); + mEh_t(1,0) = mAq_t * (**mIdq)(0,0) + mBq_t * (**mEdq_t)(1,0) + mDq_t * (**mEf) + mDq_t * mEf_prev; + SPDLOG_LOGGER_INFO(mSLog, "\n--- Model specific initialization ---" "\nSG model: 6th order type b (Anderson-Fouad's model)" @@ -59,15 +63,14 @@ void SP::Ph1::SynchronGenerator6bOrderVBR::specificInitialization() { } void SP::Ph1::SynchronGenerator6bOrderVBR::stepInPerUnit() { - if (mSimTime>0.0) { - // calculate Edq_t at t=k - (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); - (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); - - // calculate Edq_s at t=k - (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); - (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); - } + + // calculate Edq_t at t=k + (**mEdq_t)(0,0) = mAd_t * (**mIdq)(1,0) + mEh_t(0,0); + (**mEdq_t)(1,0) = mAq_t * (**mIdq)(0,0) + mEh_t(1,0); + + // calculate Edq_s at t=k + (**mEdq_s)(0,0) = -(**mIdq)(1,0) * mLq_s + (**mVdq)(0,0); + (**mEdq_s)(1,0) = (**mIdq)(0,0) * mLd_s + (**mVdq)(1,0); mDqToComplexA = get_DqToComplexATransformMatrix(); mComplexAToDq = mDqToComplexA.transpose(); diff --git a/dpsim-models/src/SimNode.cpp b/dpsim-models/src/SimNode.cpp index 557c1a19ec..71115ff8bb 100644 --- a/dpsim-models/src/SimNode.cpp +++ b/dpsim-models/src/SimNode.cpp @@ -37,6 +37,23 @@ SimNode::SimNode(PhaseType phaseType) **mVoltage = MatrixVar::Zero(3, 1); } +template <> +void SimNode::initialize() { + if (phaseType() == PhaseType::Single) + (**mVoltage)(0,0) = (RMS3PH_TO_PEAK1PH * (**mInitialVoltage)(0,0)).real(); + else + **mVoltage = (RMS3PH_TO_PEAK1PH * **mInitialVoltage).real(); +} + +template <> +void SimNode::initialize() { + (**mVoltage)(0,0) = (**mInitialVoltage)(0,0); + if (phaseType() == PhaseType::ABC) { + (**mVoltage)(1,0) = (**mInitialVoltage)(1,0); + (**mVoltage)(2,0) = (**mInitialVoltage)(2,0); + } +} + template void SimNode::initialize(Matrix frequencies) { mFrequencies = frequencies; diff --git a/dpsim-models/src/TopologicalNode.cpp b/dpsim-models/src/TopologicalNode.cpp index ebae180f37..400064b645 100644 --- a/dpsim-models/src/TopologicalNode.cpp +++ b/dpsim-models/src/TopologicalNode.cpp @@ -16,7 +16,15 @@ MatrixComp TopologicalNode::initialVoltage() const { return **mInitialVoltage; } void TopologicalNode::setInitialVoltage(MatrixComp voltage) const { **mInitialVoltage = voltage; } -void TopologicalNode::setInitialVoltage(Complex voltage) const { (**mInitialVoltage)(0,0) = voltage; } +void TopologicalNode::setInitialVoltage(Complex voltage) const { + if (mPhaseType == PhaseType::Single) { + (**mInitialVoltage)(0,0) = voltage; + } else { + (**mInitialVoltage)(0,0) = voltage; + (**mInitialVoltage)(1,0) = SHIFT_TO_PHASE_B * voltage; + (**mInitialVoltage)(2,0) = SHIFT_TO_PHASE_C * voltage; + } +} void TopologicalNode::setInitialVoltage(Complex voltage, Int phaseIndex) const { (**mInitialVoltage)(phaseIndex, 0) = voltage; diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenDQ7odTrapez_SMIB_Fault.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenDQ7odTrapez_SMIB_Fault.cpp index 18b1c52454..8cc125fdd1 100644 --- a/dpsim/examples/cxx/Circuits/EMT_SynGenDQ7odTrapez_SMIB_Fault.cpp +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenDQ7odTrapez_SMIB_Fault.cpp @@ -128,9 +128,9 @@ int main(int argc, char* argv[]) { auto gen = CPS::EMT::Ph3::SynchronGeneratorDQTrapez::make("SynGen", Logger::Level::debug); gen->setParametersFundamentalPerUnit( syngenKundur.nomPower, syngenKundur.nomVoltage, syngenKundur.nomFreq, syngenKundur.poleNum, syngenKundur.nomFieldCurr, - syngenKundur.Rs, syngenKundur.Ll, syngenKundur.Lmd, syngenKundur.Lmq, syngenKundur.Rfd, syngenKundur.Llfd, syngenKundur.Rkd, syngenKundur.Llkd, syngenKundur.Rkq1, syngenKundur.Llkq1, syngenKundur.Rkq2, syngenKundur.Llkq2, syngenKundur.H, - initActivePower, initReactivePower, initTerminalVolt, - initVoltAngle, initMechPower); + syngenKundur.Rs, syngenKundur.Ll, syngenKundur.Lmd, syngenKundur.Lmq, syngenKundur.Rfd, syngenKundur.Llfd, syngenKundur.Rkd, + syngenKundur.Llkd, syngenKundur.Rkq1, syngenKundur.Llkq1, syngenKundur.Rkq2, syngenKundur.Llkq2, syngenKundur.H, + initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, initMechPower); //Grid bus as Slack auto extnet = EMT::Ph3::NetworkInjection::make("Slack", Logger::Level::debug); diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenVBR_SMIB_Fault.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenVBR_SMIB_Fault.cpp index c11a6bcd10..5d24224bbd 100644 --- a/dpsim/examples/cxx/Circuits/EMT_SynGenVBR_SMIB_Fault.cpp +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenVBR_SMIB_Fault.cpp @@ -106,6 +106,13 @@ int main(int argc, char* argv[]) { // ----- Dynamic simulation ------ Logger::setLogDir("logs/"+simName); + // Extract relevant powerflow results + Real initTerminalVolt = std::abs(n1PF->singleVoltage())*RMS3PH_TO_PEAK1PH; + Real initVoltAngle = Math::phase(n1PF->singleVoltage()); // angle in rad + Real initActivePower = genPF->getApparentPower().real(); + Real initReactivePower = genPF->getApparentPower().imag(); + Real initMechPower = initActivePower; + // Nodes auto n1 = SimNode::make("n1", PhaseType::ABC); auto n2 = SimNode::make("n2", PhaseType::ABC); @@ -117,6 +124,7 @@ int main(int argc, char* argv[]) { syngenKundur.nomPower, syngenKundur.nomVoltage, syngenKundur.nomFreq, syngenKundur.poleNum, syngenKundur.nomFieldCurr, syngenKundur.Rs, syngenKundur.Ll, syngenKundur.Lmd, syngenKundur.Lmq, syngenKundur.Rfd, syngenKundur.Llfd, syngenKundur.Rkd, syngenKundur.Llkd, syngenKundur.Rkq1, syngenKundur.Llkq1, syngenKundur.Rkq2, syngenKundur.Llkq2, syngenKundur.H); + gen->setInitialValues(initActivePower, initReactivePower, initTerminalVolt, initVoltAngle, initMechPower); //Grid bus as Slack auto extnet = EMT::Ph3::NetworkInjection::make("Slack", Logger::Level::debug); @@ -145,7 +153,6 @@ int main(int argc, char* argv[]) { // Initialization of dynamic topology system.initWithPowerflow(systemPF); - gen->terminal(0)->setPower(-genPF->getApparentPower()); // Logging auto logger = DataLogger::make(simName); diff --git a/dpsim/include/dpsim/Definitions.h b/dpsim/include/dpsim/Definitions.h index dbb7db2136..d3d2e80e8c 100644 --- a/dpsim/include/dpsim/Definitions.h +++ b/dpsim/include/dpsim/Definitions.h @@ -10,6 +10,9 @@ #include +// This macro defines the tolerance used to compare double numbers +#define DOUBLE_EPSILON 1E-12 + namespace DPsim { // #### Types #### using Real = CPS::Real; diff --git a/dpsim/src/MNASolver.cpp b/dpsim/src/MNASolver.cpp index 22f4ff6f1d..5551b638ef 100644 --- a/dpsim/src/MNASolver.cpp +++ b/dpsim/src/MNASolver.cpp @@ -128,6 +128,10 @@ void MnaSolver::initializeComponents() { for (auto comp : mMNAIntfSwitches) comp->mnaInitialize(mSystem.mSystemOmega, mTimeStep, mLeftSideVector); + + // Initialize nodes + for (UInt nodeIdx = 0; nodeIdx < mNodes.size(); ++nodeIdx) + mNodes[nodeIdx]->initialize(); } template <> @@ -177,6 +181,10 @@ void MnaSolver::initializeComponents() { for (auto comp : mMNAIntfSwitches) comp->mnaInitialize(mSystem.mSystemOmega, mTimeStep, mLeftSideVector); + + // Initialize nodes + for (UInt nodeIdx = 0; nodeIdx < mNodes.size(); ++nodeIdx) + mNodes[nodeIdx]->initialize(); } } diff --git a/dpsim/src/Simulation.cpp b/dpsim/src/Simulation.cpp index a1c94d00c0..3f2770f859 100644 --- a/dpsim/src/Simulation.cpp +++ b/dpsim/src/Simulation.cpp @@ -347,6 +347,16 @@ void Simulation::start() { SPDLOG_LOGGER_INFO(mLog, "Time step: {:e}", **mTimeStep); SPDLOG_LOGGER_INFO(mLog, "Final time: {:e}", **mFinalTime); + // In PF we dont log the initial conditions of the componentes because they are not calculated + // In dynamic simulations log initial values of attributes (t=0) + if (mSolverType != Solver::Type::NRP) { + if (mLoggers.size() > 0) + mLoggers[0]->log(0, 0); + + // In dynamic simulations increase simulation time to calculate first results at t=timestep + mTime += **mTimeStep; + } + mSimulationStartTimePoint = std::chrono::steady_clock::now(); } @@ -369,7 +379,7 @@ void Simulation::stop() { } Real Simulation::next() { - if (mTime < **mFinalTime) + if (mTime < **mFinalTime + DOUBLE_EPSILON) step(); else stop(); @@ -381,7 +391,7 @@ Real Simulation::next() { void Simulation::run() { start(); - while (mTime < **mFinalTime) { + while (mTime < **mFinalTime + DOUBLE_EPSILON) { step(); } @@ -390,8 +400,8 @@ void Simulation::run() { Real Simulation::step() { auto start = std::chrono::steady_clock::now(); + mEvents.handleEvents(mTime); - mScheduler->step(mTime, mTimeStepCount); mTime += **mTimeStep; diff --git a/examples/Notebooks/Circuits/CS_R2CL.ipynb b/examples/Notebooks/Circuits/CS_R2CL.ipynb index 7496d77402..70fc1173f3 100644 --- a/examples/Notebooks/Circuits/CS_R2CL.ipynb +++ b/examples/Notebooks/Circuits/CS_R2CL.ipynb @@ -15,7 +15,11 @@ "source": [ "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.plottools as pt\n", - "from villas.dataprocessing.timeseries import TimeSeries as ts" + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import math \n", + "import dpsimpy\n", + "\n", + "#%matplotlib widget\n" ] }, { @@ -32,8 +36,6 @@ "outputs": [], "source": [ "# DPsim EMT simulation\n", - "import dpsimpy\n", - "\n", "name = 'EMT_CS_R2CL'\n", "\n", "# Nodes\n", @@ -41,25 +43,29 @@ "n1 = dpsimpy.emt.SimNode('n1')\n", "n2 = dpsimpy.emt.SimNode('n2')\n", "\n", + "# initialize node voltages as in modelica\n", + "n1.set_initial_voltage(complex(6.26676, -2.13813) * math.sqrt(3))\n", + "n2.set_initial_voltage(complex(-0.539123, 0.42205) * math.sqrt(3))\n", + "\n", "# Components\n", "cs = dpsimpy.emt.ph1.CurrentSource('cs')\n", "cs.I_ref = complex(10, 0)\n", "cs.f_src = 50\n", - "r1 = dpsimpy.emt.ph1.Resistor('r_1')\n", + "r1 = dpsimpy.emt.ph1.Resistor('r_1', dpsimpy.LogLevel.off)\n", "r1.R = 1\n", - "c1 = dpsimpy.emt.ph1.Capacitor('c_1')\n", + "c1 = dpsimpy.emt.ph1.Capacitor('c_1', dpsimpy.LogLevel.off)\n", "c1.C = 0.001\n", - "l1 = dpsimpy.emt.ph1.Inductor('l_1')\n", + "l1 = dpsimpy.emt.ph1.Inductor('l_1', dpsimpy.LogLevel.off)\n", "l1.L = 0.001\n", "r2 = dpsimpy.emt.ph1.Resistor('r_2')\n", "r2.R = 1\n", "\n", "# Connections\n", "cs.connect([gnd, n1])\n", - "r1.connect([n1, gnd])\n", - "c1.connect([n1, n2]);\n", - "l1.connect([n2, gnd]);\n", - "r2.connect([n2, gnd]);\n", + "r1.connect([gnd, n1])\n", + "c1.connect([n2, n1])\n", + "l1.connect([gnd, n2])\n", + "r2.connect([gnd, n2])\n", "\n", "# Define system topology\n", "system = dpsimpy.SystemTopology(50, [gnd, n1, n2], [cs, r1, c1, l1, r2])\n", @@ -70,6 +76,7 @@ "logger.log_attribute('n2.v', 'v', n2)\n", "logger.log_attribute('cs.i_intf', 'i_intf', cs)\n", "logger.log_attribute('c_1.i_intf', 'i_intf', c1)\n", + "logger.log_attribute('l_1.i_intf', 'i_intf', l1)\n", "\n", "sim = dpsimpy.Simulation(name)\n", "sim.set_domain(dpsimpy.Domain.EMT)\n", @@ -118,17 +125,17 @@ "source": [ "v1_emt = 'n1.v'\n", "v2_emt = 'n2.v'\n", - "i01_emt = 'cs.i_intf'\n", - "i12_emt = 'c_1.i_intf'\n", + "ic1_emt = 'c_1.i_intf'\n", + "il1_emt = 'l_1.i_intf'\n", "\n", "ts_dpsim_emt[v1_emt].label = 'v1 EMT'\n", "ts_dpsim_emt[v2_emt].label = 'v2 EMT'\n", - "ts_dpsim_emt[i01_emt].label = 'i01 EMT'\n", - "ts_dpsim_emt[i12_emt].label = 'i12 EMT'\n", + "ts_dpsim_emt[ic1_emt].label = 'ic1 EMT'\n", + "ts_dpsim_emt[il1_emt].label = 'il1 EMT'\n", "pt.plot_timeseries(1, ts_dpsim_emt[v1_emt])\n", "pt.plot_timeseries(1, ts_dpsim_emt[v2_emt])\n", - "pt.plot_timeseries(2, ts_dpsim_emt[i01_emt])\n", - "pt.plot_timeseries(2, ts_dpsim_emt[i12_emt])" + "pt.plot_timeseries(2, ts_dpsim_emt[ic1_emt])\n", + "pt.plot_timeseries(2, ts_dpsim_emt[il1_emt])" ] }, { @@ -145,8 +152,6 @@ "outputs": [], "source": [ "# DPsim DP simulation\n", - "import dpsimpy\n", - "\n", "name = 'DP_CS_R2CL'\n", "\n", "# Nodes\n", @@ -154,27 +159,31 @@ "n1 = dpsimpy.dp.SimNode('n1')\n", "n2 = dpsimpy.dp.SimNode('n2')\n", "\n", + "# initialize node voltages as in simulunk\n", + "n1.set_initial_voltage(complex(6.26676, -2.13813) * math.sqrt(2))\n", + "n2.set_initial_voltage(complex(-0.539123, 0.42205) * math.sqrt(2))\n", + "\n", "# Components\n", "cs = dpsimpy.dp.ph1.CurrentSource('cs')\n", "cs.I_ref = complex(10,0)\n", - "r1 = dpsimpy.dp.ph1.Resistor('r_1');\n", + "r1 = dpsimpy.dp.ph1.Resistor('r_1')\n", "r1.R = 1\n", - "c1 = dpsimpy.dp.ph1.Capacitor('c_1');\n", + "c1 = dpsimpy.dp.ph1.Capacitor('c_1')\n", "c1.C = 0.001\n", - "l1 = dpsimpy.dp.ph1.Inductor('l_1');\n", + "l1 = dpsimpy.dp.ph1.Inductor('l_1')\n", "l1.L = 0.001\n", - "r2 = dpsimpy.dp.ph1.Resistor('r_2');\n", + "r2 = dpsimpy.dp.ph1.Resistor('r_2')\n", "r2.R = 1\n", "\n", "# Connections\n", "cs.connect([gnd, n1])\n", - "r1.connect([n1, gnd])\n", - "c1.connect([n1, n2]);\n", - "l1.connect([n2, gnd]);\n", - "r2.connect([n2, gnd]);\n", + "r1.connect([gnd, n1])\n", + "c1.connect([n2, n1])\n", + "l1.connect([gnd, n2])\n", + "r2.connect([gnd, n2])\n", "\n", "# Define system topology\n", - "system = dpsimpy.SystemTopology(50, [gnd, n1, n2], [cs, r1, c1, l1, r2]);\n", + "system = dpsimpy.SystemTopology(50, [gnd, n1, n2], [cs, r1, c1, l1, r2])\n", "\n", "# Logging\n", "logger = dpsimpy.Logger(name)\n", @@ -182,6 +191,7 @@ "logger.log_attribute('n2.v', 'v', n2)\n", "logger.log_attribute('cs.i_intf', 'i_intf', cs)\n", "logger.log_attribute('c_1.i_intf', 'i_intf', c1)\n", + "logger.log_attribute('l_1.i_intf', 'i_intf', l1)\n", "\n", "sim = dpsimpy.Simulation(name)\n", "sim.set_system(system)\n", @@ -233,24 +243,25 @@ "source": [ "v1_dp = 'n1.v_shift'\n", "v2_dp = 'n2.v_shift'\n", - "i01_dp = 'cs.i_intf_shift'\n", - "i12_dp = 'c_1.i_intf_shift'\n", + "ic1_dp = 'c_1.i_intf_shift'\n", + "il1_dp = 'l_1.i_intf_shift'\n", "\n", "ts_dpsim_dp_emt[v1_dp].label = 'v1 DP'\n", "ts_dpsim_dp_emt[v2_dp].label = 'v2 DP'\n", - "ts_dpsim_dp_emt[i01_dp].label = 'i01 DP'\n", - "ts_dpsim_dp_emt[i12_dp].label = 'i12 DP'\n", + "ts_dpsim_dp_emt[ic1_dp].label = 'ic1 DP'\n", + "ts_dpsim_dp_emt[il1_dp].label = 'il1 DP'\n", "pt.plot_timeseries(1, ts_dpsim_dp_emt[v1_dp])\n", "pt.plot_timeseries(1, ts_dpsim_dp_emt[v2_dp])\n", - "pt.plot_timeseries(2, ts_dpsim_dp_emt[i01_dp])\n", - "pt.plot_timeseries(2, ts_dpsim_dp_emt[i12_dp])" + "pt.plot_timeseries(2, ts_dpsim_dp_emt[ic1_dp])\n", + "pt.plot_timeseries(2, ts_dpsim_dp_emt[il1_dp])" ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "## Simulink reference results" + "## Modelica reference results" ] }, { @@ -266,11 +277,11 @@ "if not os.path.exists('reference-results'):\n", " os.mkdir('reference-results')\n", "\n", - "url = 'https://raw.githubusercontent.com/dpsim-simulator/reference-results/master/Simulink/Circuits/SL_CS_R2CL.csv'\n", - "local_file = 'reference-results/SL_CS_R2CL.csv'\n", + "url = 'https://raw.githubusercontent.com/dpsim-simulator/reference-results/62d774a31d40663368487e2751f93ed3ee53139d/Modelica/BasicGrids/CS_R2CL_Modelica.csv'\n", + "local_file = 'reference-results/CS_R2CL_Modelica.csv'\n", "urllib.request.urlretrieve(url, local_file) \n", "\n", - "ts_sl = rt.read_timeseries_simulink(local_file)" + "ts_mod = rt.read_timeseries_simulink(local_file)" ] }, { @@ -279,19 +290,19 @@ "metadata": {}, "outputs": [], "source": [ - "v1_sl = 'v1'\n", - "v2_sl = 'v2'\n", - "i01_sl = 'i12'\n", - "i12_sl = 'i34'\n", + "v1_mod = 'capacitor.v1'\n", + "v2_mod = 'capacitor.v2'\n", + "ic1_mod = 'capacitor.i'\n", + "il1_mod = 'inductor.i'\n", "\n", - "ts_sl[v1_sl].label = 'v1 SL'\n", - "ts_sl[v2_sl].label = 'v2 SL'\n", - "ts_sl[i01_sl].label = 'i01 SL'\n", - "ts_sl[i12_sl].label = 'i12 SL'\n", - "pt.plot_timeseries(1, ts_sl[v1_sl])\n", - "pt.plot_timeseries(1, ts_sl[v2_sl])\n", - "pt.plot_timeseries(2, ts_sl[i01_sl])\n", - "pt.plot_timeseries(2, ts_sl[i12_sl])" + "ts_mod[v1_mod].label = 'v1 Mod'\n", + "ts_mod[v2_mod].label = 'v2 Mod'\n", + "ts_mod[ic1_mod].label = 'ic1 Mod'\n", + "ts_mod[il1_mod].label = 'il1 Mod'\n", + "pt.plot_timeseries(1, ts_mod[v1_mod])\n", + "pt.plot_timeseries(1, ts_mod[v2_mod])\n", + "pt.plot_timeseries(2, ts_mod[ic1_mod])\n", + "pt.plot_timeseries(2, ts_mod[il1_mod])" ] }, { @@ -307,25 +318,26 @@ "metadata": {}, "outputs": [], "source": [ + "v1_emt = 'n1.v'\n", + "v2_emt = 'n2.v'\n", + "\n", "# plot v1\n", "pt.plot_timeseries(1, ts_dpsim_emt[v1_emt])\n", "pt.plot_timeseries(1, ts_dpsim_dp_emt[v1_dp])\n", "pt.plot_timeseries(1, ts_dpsim_dp['n1.v'].abs())\n", - "pt.plot_timeseries(1, ts_sl[v1_sl])\n", + "pt.plot_timeseries(1, ts_mod[v1_mod])\n", "# plot v2\n", "pt.plot_timeseries(2, ts_dpsim_emt[v2_emt])\n", "pt.plot_timeseries(2, ts_dpsim_dp_emt[v2_dp])\n", - "pt.plot_timeseries(2, ts_sl[v2_sl])\n", - "# plot i01\n", - "pt.plot_timeseries(5, ts_dpsim_emt[i01_emt])\n", - "pt.plot_timeseries(5, ts_dpsim_dp_emt[i01_dp])\n", - "pt.plot_timeseries(5, ts_sl[i01_sl])\n", - "# plot i12\n", - "pt.plot_timeseries(6, ts_dpsim_emt[i12_emt])\n", - "pt.plot_timeseries(6, ts_dpsim_dp_emt[i12_dp])\n", - "ts_i12_sl = ts_sl[i12_sl].scale(-1)\n", - "ts_i12_sl.label = '-i12 SL'\n", - "pt.plot_timeseries(6, ts_i12_sl)" + "pt.plot_timeseries(2, ts_mod[v2_mod])\n", + "# plot ic1\n", + "pt.plot_timeseries(5, ts_dpsim_emt[ic1_emt])\n", + "pt.plot_timeseries(5, ts_dpsim_dp_emt[ic1_dp])\n", + "pt.plot_timeseries(5, ts_mod[ic1_mod])\n", + "# plot il1\n", + "pt.plot_timeseries(6, ts_dpsim_emt[il1_emt])\n", + "pt.plot_timeseries(6, ts_dpsim_dp_emt[il1_dp])\n", + "pt.plot_timeseries(6, ts_mod[il1_mod])" ] }, { @@ -334,19 +346,19 @@ "metadata": {}, "outputs": [], "source": [ - "# calculate the RMSE between Simulink (ts_sl) and EMT (ts_dpsim_emt)\n", - "err_sl_emt = 0\n", - "err_sl_emt += ts.rmse(ts_sl[v1_sl], ts_dpsim_emt[v1_emt])\n", - "err_sl_emt += ts.rmse(ts_sl[v2_sl], ts_dpsim_emt[v2_emt])\n", - "err_sl_emt = err_sl_emt / 2\n", - "print(\"Total RMSE of Simulink reference and DPsim EMT: %g\" % (err_sl_emt))\n", + "# calculate the RMSE between Modelica (ts_sl) and EMT (ts_dpsim_emt)\n", + "err_mod_emt = 0\n", + "err_mod_emt += ts.rmse(ts_mod[v1_mod], ts_dpsim_emt[v1_emt])\n", + "err_mod_emt += ts.rmse(ts_mod[v2_mod], ts_dpsim_emt[v2_emt])\n", + "err_mod_emt = err_mod_emt / 2\n", + "print(\"Total RMSE of Simulink reference and DPsim EMT: %g\" % (err_mod_emt))\n", "\n", - "# calculate the RMSE between Simulink (ts_sl) and DP (ts_dpsim_dp_emt)\n", - "err_sl_dp = 0\n", - "err_sl_dp += ts.rmse(ts_sl[v1_sl], ts_dpsim_dp_emt[v1_dp])\n", - "err_sl_dp += ts.rmse(ts_sl[v2_sl], ts_dpsim_dp_emt[v2_dp])\n", - "err_sl_dp = err_sl_dp / 2\n", - "print(\"Total RMSE of Simulink reference and DPsim DP: %g\" % (err_sl_dp))" + "# calculate the RMSE between Modelica (ts_sl) and DP (ts_dpsim_dp_emt)\n", + "err_mod_dp = 0\n", + "err_mod_dp += ts.rmse(ts_mod[v1_mod], ts_dpsim_dp_emt[v1_dp])\n", + "err_mod_dp += ts.rmse(ts_mod[v2_mod], ts_dpsim_dp_emt[v2_dp])\n", + "err_mod_dp = err_mod_dp / 2\n", + "print(\"Total RMSE of Simulink reference and DPsim DP: %g\" % (err_mod_dp))" ] }, { @@ -355,16 +367,9 @@ "metadata": {}, "outputs": [], "source": [ - "assert err_sl_emt < 0.1\n", - "assert err_sl_dp < 0.1" + "assert err_mod_emt < 0.00017\n", + "assert err_mod_dp < 2.26e-06" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -383,7 +388,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.9" + "version": "3.9.13" }, "tests": { "skip": false diff --git a/examples/Notebooks/Circuits/Compare_DP_SMIB_ReducedOrderSG_VBR_TPM_LoadStep.ipynb b/examples/Notebooks/Circuits/Compare_DP_SMIB_ReducedOrderSG_VBR_TPM_LoadStep.ipynb index c1205244e5..8b698f4583 100644 --- a/examples/Notebooks/Circuits/Compare_DP_SMIB_ReducedOrderSG_VBR_TPM_LoadStep.ipynb +++ b/examples/Notebooks/Circuits/Compare_DP_SMIB_ReducedOrderSG_VBR_TPM_LoadStep.ipynb @@ -1,6 +1,7 @@ { "cells": [ { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -8,6 +9,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -34,6 +36,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -82,6 +85,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -119,6 +123,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -165,6 +170,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -172,6 +178,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -196,6 +203,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -220,6 +228,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -240,6 +249,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -259,13 +269,6 @@ "assert(rmse_list[4]<0.003)\n", "assert(rmse_list[5]<2.4e-5)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/examples/Notebooks/Circuits/Compare_EMT_SynGenDQ7odTrapez_EMT_SynGenVBR_SMIB_Fault.ipynb b/examples/Notebooks/Circuits/Compare_EMT_SynGenDQ7odTrapez_EMT_SynGenVBR_SMIB_Fault.ipynb index 32432e084c..2efd159990 100644 --- a/examples/Notebooks/Circuits/Compare_EMT_SynGenDQ7odTrapez_EMT_SynGenVBR_SMIB_Fault.ipynb +++ b/examples/Notebooks/Circuits/Compare_EMT_SynGenDQ7odTrapez_EMT_SynGenVBR_SMIB_Fault.ipynb @@ -328,19 +328,12 @@ "for name in ['v_gen_0', 'v_gen_1', 'v_gen_2']:\n", " diff[name] = ts_dcim[name].rmse(ts_dcim[name], ts_vbr[name])/np.max(ts_dcim[name].values)\n", " print(name + ': ' + str(diff[name]))\n", - " assert(diff[name]) < 1e-5\n", + " assert(diff[name]) < 1.86e-6\n", "for name in ['i_gen_0', 'i_gen_1', 'i_gen_2']:\n", " diff[name] = ts_dcim[name].rmse(ts_dcim[name], ts_vbr[name])/np.max(ts_dcim[name].values)\n", " print(name + ': ' + str(diff[name]))\n", - " assert(diff[name]) < 1e-3" + " assert(diff[name]) < 5e-4" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -359,7 +352,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.9.13" }, "tests": { "skip": false diff --git a/examples/Notebooks/Circuits/SP_Validation_ReducedOrderSG_VBR_SMIB_Fault_withPSAT.ipynb b/examples/Notebooks/Circuits/SP_Validation_ReducedOrderSG_VBR_SMIB_Fault_withPSAT.ipynb index ff0fe7a733..b238975bc9 100644 --- a/examples/Notebooks/Circuits/SP_Validation_ReducedOrderSG_VBR_SMIB_Fault_withPSAT.ipynb +++ b/examples/Notebooks/Circuits/SP_Validation_ReducedOrderSG_VBR_SMIB_Fault_withPSAT.ipynb @@ -517,7 +517,7 @@ "metadata": {}, "outputs": [], "source": [ - "def sp_reducedOrderSG_SMIB_fault(sim_pf, gen_pf, gen_model, final_time=20, time_step=0.001):\n", + "def sp_reducedOrderSG_SMIB_fault(sim_pf, gen_pf, gen_model, final_time=10, time_step=0.001):\n", " \n", " # ## Extract relevant powerflow results\n", " init_electrical_power = gen_pf.get_apparent_power()\n", diff --git a/examples/Notebooks/Circuits/VS_CS_R4.ipynb b/examples/Notebooks/Circuits/VS_CS_R4.ipynb index a6eaf718bd..a9b0e8893d 100644 --- a/examples/Notebooks/Circuits/VS_CS_R4.ipynb +++ b/examples/Notebooks/Circuits/VS_CS_R4.ipynb @@ -15,7 +15,11 @@ "source": [ "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.plottools as pt\n", - "from villas.dataprocessing.timeseries import TimeSeries as ts" + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import math\n", + "import dpsimpy\n", + "\n", + "#%matplotlib widget" ] }, { @@ -32,28 +36,30 @@ "outputs": [], "source": [ "# DPsim EMT simulation\n", - "import dpsimpy\n", - "\n", "name = 'EMT_VS_CS_R4_AC'\n", "\n", "# Nodes\n", + "# Obs.: in DPsim are initial values 3Ph, RMS quantities\n", "gnd = dpsimpy.emt.SimNode.gnd\n", "n1 = dpsimpy.emt.SimNode('n1')\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", "n2 = dpsimpy.emt.SimNode('n2')\n", + "n2.set_initial_voltage(complex(3.53553, 0) * math.sqrt(3))\n", "n3 = dpsimpy.emt.SimNode('n3')\n", + "n3.set_initial_voltage(complex(3.53553, 0) * math.sqrt(3))\n", "\n", "# Components\n", - "vs = dpsimpy.emt.ph1.VoltageSource('vs')\n", + "vs = dpsimpy.emt.ph1.VoltageSource('vs', dpsimpy.LogLevel.debug)\n", "vs.set_parameters(V_ref=complex(10, 0), f_src=50)\n", - "r1 = dpsimpy.emt.ph1.Resistor('r1')\n", + "r1 = dpsimpy.emt.ph1.Resistor('r1', dpsimpy.LogLevel.debug)\n", "r1.set_parameters(R=1)\n", - "r2 = dpsimpy.emt.ph1.Resistor('r2')\n", + "r2 = dpsimpy.emt.ph1.Resistor('r2', dpsimpy.LogLevel.debug)\n", "r2.set_parameters(R=1)\n", - "r3 = dpsimpy.emt.ph1.Resistor('r3')\n", + "r3 = dpsimpy.emt.ph1.Resistor('r3', dpsimpy.LogLevel.debug)\n", "r3.set_parameters(R=10)\n", - "r4 = dpsimpy.emt.ph1.Resistor('r4')\n", + "r4 = dpsimpy.emt.ph1.Resistor('r4', dpsimpy.LogLevel.debug)\n", "r4.set_parameters(R=5)\n", - "cs = dpsimpy.emt.ph1.CurrentSource('cs')\n", + "cs = dpsimpy.emt.ph1.CurrentSource('cs', dpsimpy.LogLevel.debug)\n", "cs.set_parameters(I_ref=complex(1,0), f_src=50)\n", "\n", "vs.connect([gnd, n1])\n", @@ -66,11 +72,11 @@ "system = dpsimpy.SystemTopology(50, [gnd, n1, n2, n3], [vs, r1, r2, r3, r4, cs])\n", "\n", "logger = dpsimpy.Logger(name)\n", - "logger.log_attribute('n1.v', 'v', n1);\n", - "logger.log_attribute('n2.v', 'v', n2);\n", - "logger.log_attribute('n3.v', 'v', n3);\n", - "logger.log_attribute('r1.i_intf', 'i_intf', r1);\n", - "logger.log_attribute('r3.i_intf', 'i_intf', r3);\n", + "logger.log_attribute('n1.v', 'v', n1)\n", + "logger.log_attribute('n2.v', 'v', n2)\n", + "logger.log_attribute('n3.v', 'v', n3)\n", + "logger.log_attribute('r1.i_intf', 'i_intf', r1)\n", + "logger.log_attribute('r3.i_intf', 'i_intf', r3)\n", "\n", "sim = dpsimpy.Simulation(name)\n", "sim.set_domain(dpsimpy.Domain.EMT)\n", @@ -128,19 +134,21 @@ "outputs": [], "source": [ "# DPsim DP simulation\n", - "import dpsimpy\n", - "\n", "name = 'DP_VS_CS_R4_AC'\n", "\n", "# Nodes\n", + "# Obs.: in DPsim are initial values 3Ph, RMS quantities\n", "gnd = dpsimpy.dp.SimNode.gnd\n", "n1 = dpsimpy.dp.SimNode('n1')\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", "n2 = dpsimpy.dp.SimNode('n2')\n", + "n2.set_initial_voltage(complex(3.53553, 0) * math.sqrt(3))\n", "n3 = dpsimpy.dp.SimNode('n3')\n", + "n3.set_initial_voltage(complex(3.53553, 0) * math.sqrt(3))\n", "\n", "# Components\n", "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", - "vs.set_parameters(V_ref=complex(10, 0))\n", + "vs.set_parameters(V_ref=complex(10, 0) * math.sqrt(3/2))\n", "r1 = dpsimpy.dp.ph1.Resistor('r1')\n", "r1.set_parameters(R=1)\n", "r2 = dpsimpy.dp.ph1.Resistor('r2')\n", @@ -150,7 +158,7 @@ "r4 = dpsimpy.dp.ph1.Resistor('r4')\n", "r4.set_parameters(R=5)\n", "cs = dpsimpy.dp.ph1.CurrentSource('cs')\n", - "cs.set_parameters(I_ref=complex(1,0))\n", + "cs.set_parameters(I_ref=complex(1,0) * math.sqrt(3/2))\n", "\n", "vs.connect([gnd, n1])\n", "r1.connect([n1, n2])\n", @@ -162,11 +170,11 @@ "system = dpsimpy.SystemTopology(50, [gnd, n1, n2, n3], [vs, r1, r2, r3, r4, cs])\n", "\n", "logger = dpsimpy.Logger(name)\n", - "logger.log_attribute('n1.v', 'v', n1);\n", - "logger.log_attribute('n2.v', 'v', n2);\n", - "logger.log_attribute('n3.v', 'v', n3);\n", - "logger.log_attribute('r1.i_intf', 'i_intf', r1);\n", - "logger.log_attribute('r3.i_intf', 'i_intf', r3);\n", + "logger.log_attribute('n1.v', 'v', n1)\n", + "logger.log_attribute('n2.v', 'v', n2)\n", + "logger.log_attribute('n3.v', 'v', n3)\n", + "logger.log_attribute('r1.i_intf', 'i_intf', r1)\n", + "logger.log_attribute('r3.i_intf', 'i_intf', r3)\n", "\n", "sim = dpsimpy.Simulation(name)\n", "sim.set_domain(dpsimpy.Domain.DP)\n", @@ -199,7 +207,11 @@ "outputs": [], "source": [ "#convert to emt\n", - "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)" + "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)\n", + "\n", + "# 3Ph, RMS to 1Ph \n", + "for name, ts in ts_dpsim_dp_emt.items():\n", + " ts.values = ts.values * math.sqrt(2/3) " ] }, { @@ -287,8 +299,8 @@ "pt.plot_timeseries(2, ts_dpsim_dp_emt['n2.v_shift'])\n", "pt.plot_timeseries(2, ts_sl['v2'])\n", "# plot v3\n", - "pt.plot_timeseries(3, ts_dpsim_emt['n1.v'])\n", - "pt.plot_timeseries(3, ts_dpsim_dp_emt['n1.v_shift'])\n", + "pt.plot_timeseries(3, ts_dpsim_emt['n3.v'])\n", + "pt.plot_timeseries(3, ts_dpsim_dp_emt['n3.v_shift'])\n", "pt.plot_timeseries(3, ts_sl['v3'])\n", "# plot i12\n", "pt.plot_timeseries(4, ts_dpsim_emt['r1.i_intf'])\n", @@ -331,16 +343,9 @@ "metadata": {}, "outputs": [], "source": [ - "assert err_sl_emt < 0.1\n", - "assert err_sl_dp < 0.1" + "assert err_sl_emt < 3.4e-07\n", + "assert err_sl_dp < 2.31e-07" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -359,7 +364,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.9" + "version": "3.9.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { diff --git a/examples/Notebooks/Circuits/VS_R2L3.ipynb b/examples/Notebooks/Circuits/VS_R2L3.ipynb index d800cf4546..40e61a5843 100644 --- a/examples/Notebooks/Circuits/VS_R2L3.ipynb +++ b/examples/Notebooks/Circuits/VS_R2L3.ipynb @@ -15,7 +15,11 @@ "source": [ "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.plottools as pt\n", - "from villas.dataprocessing.timeseries import TimeSeries as ts" + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import math\n", + "import dpsimpy\n", + "\n", + "#%matplotlib widget" ] }, { @@ -32,8 +36,6 @@ "outputs": [], "source": [ "# DPsim EMT simulation\n", - "import dpsimpy\n", - "\n", "name = 'EMT_VS_R2L3'\n", "\n", "# Nodes\n", @@ -43,6 +45,13 @@ "n3 = dpsimpy.emt.SimNode('n3')\n", "n4 = dpsimpy.emt.SimNode('n4')\n", "\n", + "# initialize node voltages as in simulink\n", + "# Obs.: in DPsim are initial values 3Ph, RMS quantities\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", + "n2.set_initial_voltage(complex(7.02433, 0.415825) * math.sqrt(3))\n", + "n3.set_initial_voltage(complex(4.41163, 0.122191) * math.sqrt(3))\n", + "n4.set_initial_voltage(complex(0.0856874, -0.550796) * math.sqrt(3))\n", + "\n", "# Components\n", "vs = dpsimpy.emt.ph1.VoltageSource('vs')\n", "vs.set_parameters(V_ref=complex(10,0),f_src=50)\n", @@ -59,11 +68,11 @@ "\n", "# Connections\n", "vs.connect([gnd, n1])\n", - "r1.connect([n1, n2])\n", - "l1.connect([n2, n3])\n", - "l2.connect([n3, gnd])\n", - "l3.connect([n3, n4])\n", - "r2.connect([n4, gnd]);\n", + "r1.connect([n2, n1])\n", + "l1.connect([n3, n2])\n", + "l2.connect([gnd, n3])\n", + "l3.connect([n4, n3])\n", + "r2.connect([gnd, n4,])\n", "\n", "# Define system topology\n", "system = dpsimpy.SystemTopology(50, [gnd, n1, n2, n3, n4], [vs, r1, l1, l2, l3, r2]);\n", @@ -142,8 +151,6 @@ "outputs": [], "source": [ "# DPsim DP simulation\n", - "import dpsimpy\n", - "\n", "name = 'DP_VS_R2L3'\n", "\n", "# Nodes\n", @@ -153,9 +160,15 @@ "n3 = dpsimpy.dp.SimNode('n3')\n", "n4 = dpsimpy.dp.SimNode('n4')\n", "\n", + "# initialize node voltages as in simulink\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", + "n2.set_initial_voltage(complex(7.02433, 0.415825) * math.sqrt(3))\n", + "n3.set_initial_voltage(complex(4.41163, 0.122191) * math.sqrt(3))\n", + "n4.set_initial_voltage(complex(0.0856874, -0.550796) * math.sqrt(3))\n", + "\n", "# Components\n", "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", - "vs.set_parameters(V_ref=complex(10,0))\n", + "vs.set_parameters(V_ref=complex(10,0) * math.sqrt(3/2))\n", "r1 = dpsimpy.dp.ph1.Resistor('r_1')\n", "r1.set_parameters(R=1)\n", "l1 = dpsimpy.dp.ph1.Inductor('l_1')\n", @@ -218,7 +231,11 @@ "outputs": [], "source": [ "#convert to emt\n", - "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)" + "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)\n", + "\n", + "# 3Ph, RMS to 1Ph \n", + "for name, ts in ts_dpsim_dp_emt.items():\n", + " ts.values = ts.values * math.sqrt(2/3)" ] }, { @@ -362,8 +379,8 @@ "metadata": {}, "outputs": [], "source": [ - "assert err_sl_emt < 0.1\n", - "assert err_sl_dp < 0.1" + "assert err_sl_emt < 0.0172\n", + "assert err_sl_dp < 0.0172" ] } ], @@ -383,7 +400,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.9" + "version": "3.9.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { diff --git a/examples/Notebooks/Circuits/VS_RC1.ipynb b/examples/Notebooks/Circuits/VS_RC1.ipynb index 4475942c2d..c2c1415b65 100644 --- a/examples/Notebooks/Circuits/VS_RC1.ipynb +++ b/examples/Notebooks/Circuits/VS_RC1.ipynb @@ -15,7 +15,11 @@ "source": [ "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.plottools as pt\n", - "from villas.dataprocessing.timeseries import TimeSeries as ts" + "from villas.dataprocessing.timeseries import TimeSeries as ts\n", + "import math\n", + "import dpsimpy\n", + "\n", + "#%matplotlib widget" ] }, { @@ -32,8 +36,6 @@ "outputs": [], "source": [ "# DPsim EMT simulation\n", - "import dpsimpy\n", - "\n", "name = 'EMT_VS_RC1'\n", "\n", "# Nodes\n", @@ -41,6 +43,11 @@ "n1 = dpsimpy.emt.SimNode('n1')\n", "n2 = dpsimpy.emt.SimNode('n2')\n", "\n", + "# initialize node voltages as in modelica\n", + "# Obs.: in DPsim are initial values 3Ph, RMS quantities\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", + "n2.set_initial_voltage(complex(6.43587, 2.02189) * math.sqrt(3))\n", + "\n", "# Components\n", "vs = dpsimpy.emt.ph1.VoltageSource('vs')\n", "vs.set_parameters(V_ref=complex(10,0),f_src=50)\n", @@ -115,8 +122,6 @@ "outputs": [], "source": [ "# DPsim DP simulation\n", - "import dpsimpy\n", - "\n", "name = 'DP_VS_RC1'\n", "\n", "# Nodes\n", @@ -124,9 +129,13 @@ "n1 = dpsimpy.dp.SimNode('n1')\n", "n2 = dpsimpy.dp.SimNode('n2')\n", "\n", + "# initialize node voltages as in simulink\n", + "n1.set_initial_voltage(complex(10, 0) * math.sqrt(3/2))\n", + "n2.set_initial_voltage(complex(6.43587, 2.02189) * math.sqrt(3))\n", + "\n", "# Components\n", "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", - "vs.set_parameters(V_ref=complex(10,0))\n", + "vs.set_parameters(V_ref=complex(10,0) * math.sqrt(3/2))\n", "r1 = dpsimpy.dp.ph1.Resistor('r_1')\n", "r1.set_parameters(R=1)\n", "c1 = dpsimpy.dp.ph1.Capacitor('c_1')\n", @@ -177,7 +186,11 @@ "outputs": [], "source": [ "# convert to emt\n", - "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)" + "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)\n", + "\n", + "# 3Ph, RMS to 1Ph \n", + "for name, ts in ts_dpsim_dp_emt.items():\n", + " ts.values = ts.values * math.sqrt(2/3) " ] }, { @@ -214,11 +227,11 @@ "if not os.path.exists('reference-results'):\n", " os.mkdir('reference-results')\n", "\n", - "url = 'https://raw.githubusercontent.com/dpsim-simulator/reference-results/master/Simulink/Circuits/SL_VS_RC1.csv'\n", - "local_file = 'reference-results/SL_VS_RC1.csv'\n", + "url = 'https://raw.githubusercontent.com/dpsim-simulator/reference-results/master/Modelica/BasicGrids/VS_RC1_Modelica.csv'\n", + "local_file = 'reference-results/VS_RC1_Modelica.csv'\n", "urllib.request.urlretrieve(url, local_file) \n", "\n", - "ts_sl = rt.read_timeseries_simulink(local_file)" + "ts_mod = rt.read_timeseries_simulink(local_file)" ] }, { @@ -227,12 +240,12 @@ "metadata": {}, "outputs": [], "source": [ - "pt.set_timeseries_labels(ts_sl['v1'], 'v1 SL')\n", - "pt.set_timeseries_labels(ts_sl['v2'], 'v2 SL')\n", - "pt.set_timeseries_labels(ts_sl['i12'], 'i12 SL')\n", - "pt.plot_timeseries(1, ts_sl['v1'])\n", - "pt.plot_timeseries(1, ts_sl['v2'])\n", - "pt.plot_timeseries(2, ts_sl['i12'])" + "pt.set_timeseries_labels(ts_mod['v1'], 'v1 Mod')\n", + "pt.set_timeseries_labels(ts_mod['v2'], 'v2 Mod')\n", + "pt.set_timeseries_labels(ts_mod['i'], 'i12 Mod')\n", + "pt.plot_timeseries(1, ts_mod['v1'])\n", + "pt.plot_timeseries(1, ts_mod['v2'])\n", + "pt.plot_timeseries(2, ts_mod['i'])" ] }, { @@ -251,16 +264,16 @@ "# plot v1\n", "pt.plot_timeseries(1, ts_dpsim_emt['n1.v'])\n", "pt.plot_timeseries(1, ts_dpsim_dp_emt['n1.v_shift'])\n", - "pt.plot_timeseries(1, ts_sl['v1'])\n", + "pt.plot_timeseries(1, ts_mod['v1'])\n", "# plot v2\n", "pt.plot_timeseries(2, ts_dpsim_emt['n2.v'])\n", "pt.plot_timeseries(2, ts_dpsim_dp_emt['n2.v_shift'])\n", - "pt.plot_timeseries(2, ts_sl['v2'])\n", + "pt.plot_timeseries(2, ts_mod['v2'])\n", "# plot i12\n", "pt.plot_timeseries(3, ts_dpsim_emt['r_1.i_intf'])\n", "pt.plot_timeseries(3, ts_dpsim_dp_emt['r_1.i_intf_shift'])\n", - "ts_sl_scale = ts_sl['i12'].scale(-1)\n", - "ts_sl_scale.label = '-i12 SL'\n", + "ts_sl_scale = ts_mod['i'].scale(-1)\n", + "ts_sl_scale.label = '-i Mod'\n", "pt.plot_timeseries(3, ts_sl_scale)" ] }, @@ -272,17 +285,17 @@ "source": [ "# calculate the RMSE between Simulink (ts_sl) and EMT (ts_dpsim_emt)\n", "err_sl_emt = 0\n", - "err_sl_emt += ts.rmse(ts_sl['v1'], ts_dpsim_emt['n1.v'])\n", - "err_sl_emt += ts.rmse(ts_sl['v2'], ts_dpsim_emt['n2.v'])\n", + "err_sl_emt += ts.rmse(ts_mod['v1'], ts_dpsim_emt['n1.v'])\n", + "err_sl_emt += ts.rmse(ts_mod['v2'], ts_dpsim_emt['n2.v'])\n", "err_sl_emt = err_sl_emt / 2\n", - "print(\"Total RMSE of Simulink reference and DPsim EMT: %g\" % (err_sl_emt))\n", + "print(\"Total RMSE of Modelica reference and DPsim EMT: %g\" % (err_sl_emt))\n", "\n", "# calculate the RMSE between Simulink (ts_sl) and DP (ts_dpsim_dp_emt)\n", "err_sl_dp = 0\n", - "err_sl_dp += ts.rmse(ts_sl['v1'], ts_dpsim_dp_emt['n1.v_shift'])\n", - "err_sl_dp += ts.rmse(ts_sl['v2'], ts_dpsim_dp_emt['n2.v_shift'])\n", + "err_sl_dp += ts.rmse(ts_mod['v1'], ts_dpsim_dp_emt['n1.v_shift'])\n", + "err_sl_dp += ts.rmse(ts_mod['v2'], ts_dpsim_dp_emt['n2.v_shift'])\n", "err_sl_dp = err_sl_dp / 2\n", - "print(\"Total RMSE of Simulink reference and DPsim DP: %g\" % (err_sl_dp))" + "print(\"Total RMSE of Modelica reference and DPsim DP: %g\" % (err_sl_dp))" ] }, { @@ -291,8 +304,8 @@ "metadata": {}, "outputs": [], "source": [ - "assert err_sl_emt < 0.1\n", - "assert err_sl_dp < 0.1" + "assert err_sl_emt < 0.062\n", + "assert err_sl_dp < 0.062" ] } ], @@ -312,7 +325,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.9" + "version": "3.9.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { diff --git a/examples/Notebooks/Circuits/VS_RL1.ipynb b/examples/Notebooks/Circuits/VS_RL1.ipynb index d41f6a1114..5f46d005a7 100644 --- a/examples/Notebooks/Circuits/VS_RL1.ipynb +++ b/examples/Notebooks/Circuits/VS_RL1.ipynb @@ -13,6 +13,7 @@ "metadata": {}, "outputs": [], "source": [ + "import math\n", "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.plottools as pt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", @@ -42,7 +43,7 @@ "\n", "# Components\n", "vs_pf = dpsimpy.sp.ph1.VoltageSource('vs')\n", - "vs_pf.set_parameters(V_ref=complex(10,0), f_src=50)\n", + "vs_pf.set_parameters(V_ref=complex(10,0) * math.sqrt(3/2), f_src=50)\n", "r1_pf = dpsimpy.sp.ph1.Resistor('r1_pf')\n", "r1_pf.set_parameters(R=5)\n", "l1_pf = dpsimpy.sp.ph1.Inductor('l1_pf')\n", @@ -132,7 +133,6 @@ "outputs": [], "source": [ "# read EMT results\n", - "#work_dir = '../../dpsim/Logs/'\n", "work_dir = 'logs/EMT_VS_RL1/'\n", "log_name = 'EMT_VS_RL1'\n", "print(work_dir + log_name + '.csv')\n", @@ -183,7 +183,7 @@ "\n", "# Components\n", "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n", - "vs.set_parameters(V_ref=complex(10,0))\n", + "vs.set_parameters(V_ref=complex(10,0) * math.sqrt(3/2))\n", "r1 = dpsimpy.dp.ph1.Resistor('r1')\n", "r1.set_parameters(R=5)\n", "l1 = dpsimpy.dp.ph1.Inductor('l1')\n", @@ -231,7 +231,11 @@ "outputs": [], "source": [ "# convert to emt\n", - "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)" + "ts_dpsim_dp_emt = ts.frequency_shift_list(ts_dpsim_dp, 50)\n", + "\n", + "# 3Ph, RMS to 1Ph \n", + "for name, ts in ts_dpsim_dp_emt.items():\n", + " ts.values = ts.values * math.sqrt(2/3)" ] }, { @@ -302,9 +306,6 @@ "metadata": {}, "outputs": [], "source": [ - "#v1 = 'v1'\n", - "#v2 = 'v2'\n", - "#i12 = 'i12'\n", "v1 = 'n1.v'\n", "v2 = 'n2.v'\n", "i12 = 'r1.i_intf'\n", @@ -352,8 +353,8 @@ "metadata": {}, "outputs": [], "source": [ - "assert err_mod_emt < 0.009\n", - "assert err_mod_dp < 0.002" + "assert err_mod_emt < 0.00012\n", + "assert err_mod_dp < 0.0000042" ] } ], diff --git a/examples/Notebooks/Components/Inverter_Grid_Test.ipynb b/examples/Notebooks/Components/Inverter_Grid_Test.ipynb index 2f448862b6..0c2d392be4 100644 --- a/examples/Notebooks/Components/Inverter_Grid_Test.ipynb +++ b/examples/Notebooks/Components/Inverter_Grid_Test.ipynb @@ -235,6 +235,7 @@ "source": [ "import numpy as np\n", "diff = ts_ml_inverter['v_c'].values[:5000] - np.interp(ts_ml_inverter['v_c'].time[:5000], ts_dpsim_emt_v3.time, ts_dpsim_emt_v3.values)\n", + "print(max(diff))\n", "assert np.all(diff) < 10" ] }, @@ -334,7 +335,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.9.13" }, "orig_nbformat": 3 }, diff --git a/examples/Notebooks/Components/SynGenDq7od.ipynb b/examples/Notebooks/Components/SynGenDq7od.ipynb index 76fb4ca3a4..6eeb698fd3 100644 --- a/examples/Notebooks/Components/SynGenDq7od.ipynb +++ b/examples/Notebooks/Components/SynGenDq7od.ipynb @@ -1,6 +1,7 @@ { "cells": [ { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -17,11 +18,12 @@ "import villas.dataprocessing.plottools as pt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", "import matplotlib.pyplot as plt\n", - "import dpsimpy\n", - "import numpy as np" + "import numpy as np\n", + "import dpsimpy" ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -156,6 +158,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -195,6 +198,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -217,6 +221,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -224,6 +229,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -275,6 +281,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -434,6 +441,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -493,10 +501,12 @@ "source": [ "import numpy as np\n", "diff = ts_sl_tpf_i1.values - ts_dpsim_ode_tpf_emt['i_gen_0_shift'].values\n", - "assert np.max(diff[:4000]) < 200" + "print(max(diff[:4000]))\n", + "assert np.max(diff[:4000]) < 28.5" ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -623,6 +633,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -656,7 +667,8 @@ "source": [ "import numpy as np\n", "diff = ts_sl_tpf_i1.values - ts_dpsim_trapez_tpf_emt['i_gen_0_shift'].values\n", - "assert np.max(diff[:4000]) < 200" + "print(np.max(diff[:4000]))\n", + "assert np.max(diff[:4000]) < 30.6" ] } ], @@ -675,7 +687,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.9.13" }, "orig_nbformat": 3 }, diff --git a/examples/Notebooks/Components/SynGenDq7od_ThreePhFault_DP_EMT.ipynb b/examples/Notebooks/Components/SynGenDq7od_ThreePhFault_DP_EMT.ipynb index b435dc8e71..780edc6f62 100644 --- a/examples/Notebooks/Components/SynGenDq7od_ThreePhFault_DP_EMT.ipynb +++ b/examples/Notebooks/Components/SynGenDq7od_ThreePhFault_DP_EMT.ipynb @@ -1,6 +1,7 @@ { "cells": [ { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -8,6 +9,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -23,11 +25,12 @@ "import villas.dataprocessing.readtools as rt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", "import matplotlib.pyplot as plt\n", - "import dpsimpy\n", - "import numpy as np" + "import numpy as np\n", + "import dpsimpy" ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -131,6 +134,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -234,6 +238,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -335,6 +340,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -437,6 +443,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -507,6 +514,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -554,6 +562,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -601,6 +610,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -655,10 +665,12 @@ "source": [ "import numpy as np\n", "diff = ts_sl_tpf_i1.values - ts_dpsim_ode_tpf_emt['i_gen_0_shift'].values\n", - "assert np.max(diff[:4000]) < 200" + "print(np.max(diff[:4000]))\n", + "assert np.max(diff[:4000]) < 28.5" ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -693,7 +705,8 @@ "source": [ "import numpy as np\n", "diff = ts_sl_tpf_i1.values - ts_dpsim_trapez_tpf_emt['i_gen_0_shift'].values\n", - "assert np.max(diff[:4000]) < 200" + "print(np.max(diff[:4000]))\n", + "assert np.max(diff[:4000]) < 30.6" ] } ], @@ -713,7 +726,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.9.13" }, "tests": { "skip": false diff --git a/examples/Notebooks/Components/Trafo.ipynb b/examples/Notebooks/Components/Trafo.ipynb index 8ca3d82bfa..e54a49a58c 100644 --- a/examples/Notebooks/Components/Trafo.ipynb +++ b/examples/Notebooks/Components/Trafo.ipynb @@ -957,7 +957,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.9.13" } }, "nbformat": 4, diff --git a/examples/Notebooks/Grids/DP_WSCC9bus_SGTrStab.ipynb b/examples/Notebooks/Grids/DP_WSCC9bus_SGTrStab.ipynb index 22df4fdce6..204e721545 100644 --- a/examples/Notebooks/Grids/DP_WSCC9bus_SGTrStab.ipynb +++ b/examples/Notebooks/Grids/DP_WSCC9bus_SGTrStab.ipynb @@ -41,8 +41,8 @@ "import villas.dataprocessing.readtools as rt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", "import matplotlib.pyplot as plt\n", - "import dpsimpy\n", - "import numpy as np\n" + "import numpy as np\n", + "import dpsimpy" ] }, { @@ -268,7 +268,7 @@ "source": [ "for node, phasor in phasors.items():\n", " abs_diff = phasor['abs'].values[-1] - phasor['abs'].values[0]\n", - " assert abs_diff < 20\n", + " assert abs_diff < 27\n", " phase_diff = phasor['phase'].values[-1] - phasor['phase'].values[0]\n", " assert phase_diff < 0.001\n", " print(node + ': ' + str(abs_diff) + '<' + str(phase_diff))" diff --git a/examples/Notebooks/Grids/IEEE_LV_powerflow.ipynb b/examples/Notebooks/Grids/IEEE_LV_powerflow.ipynb index 487c42067a..556dfcf315 100644 --- a/examples/Notebooks/Grids/IEEE_LV_powerflow.ipynb +++ b/examples/Notebooks/Grids/IEEE_LV_powerflow.ipynb @@ -184,7 +184,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.12" + "version": "3.9.13" } }, "nbformat": 4,