diff --git a/dpsim-models/include/dpsim-models/Base/Base_Ph3_Switch.h b/dpsim-models/include/dpsim-models/Base/Base_Ph3_Switch.h index 5f44ebe279..83610edf01 100644 --- a/dpsim-models/include/dpsim-models/Base/Base_Ph3_Switch.h +++ b/dpsim-models/include/dpsim-models/Base/Base_Ph3_Switch.h @@ -36,6 +36,9 @@ namespace Ph3 { } void closeSwitch() { **mSwitchClosed = true; } void openSwitch() { **mSwitchClosed = false; } + + /// Check if switch is closed + Bool isClosed() { return **mSwitchClosed; } }; } } diff --git a/dpsim-models/include/dpsim-models/Components.h b/dpsim-models/include/dpsim-models/Components.h index 3326608c27..2aeb3f7907 100644 --- a/dpsim-models/include/dpsim-models/Components.h +++ b/dpsim-models/include/dpsim-models/Components.h @@ -115,6 +115,7 @@ #include #include #include +#include #include #include diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorTrStab.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorTrStab.h index 85f4042da1..1874278f7d 100644 --- a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorTrStab.h +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_SynchronGeneratorTrStab.h @@ -42,7 +42,11 @@ namespace Ph3 { Matrix mStates; /// Nominal system angle Real mThetaN = 0; - + /// Flag for usage of attribute of w_ref (otherwise mNomOmega is used) + Bool mUseOmegaRef = false; + /// Flag for usage of actual mechanical speed for torque conversion (otherwise mNomOmega is used) + Bool mConvertWithOmegaMech = true; + public: // #### Model specific variables #### /// emf behind transient reactance @@ -73,6 +77,8 @@ namespace Ph3 { Matrix getParkTransformMatrixPowerInvariant(Real theta); // #### General Functions #### + /// Flags to modify model behavior + void setModelFlags(Bool convertWithOmegaMech); /// void setInitialValues(Complex elecPower, Real mechPower); /// \brief Initializes the machine parameters @@ -105,6 +111,8 @@ namespace Ph3 { void mnaParentAddPreStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes) override; void mnaParentAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + void setReferenceOmega(Attribute::Ptr refOmegaPtr, Attribute::Ptr refDeltaPtr); + class AddBStep : public Task { public: AddBStep(SynchronGeneratorTrStab& generator) : diff --git a/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_varResSwitch.h b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_varResSwitch.h new file mode 100644 index 0000000000..a53d1f88df --- /dev/null +++ b/dpsim-models/include/dpsim-models/EMT/EMT_Ph3_varResSwitch.h @@ -0,0 +1,89 @@ +/* Copyright 2017-2020 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 + +namespace CPS { +namespace EMT { +namespace Ph3 { + /// \brief + /// + /// Switch with variable resistance to avoid numerical oscillations, when an inductive current is suddenly interrupted. + /// It is useful as a fault switch especially on a faulted generator or transformer winding. + + /// The switch resistance changes at a defined fixed rate by multiplying previous resistance value with a factor for the rate of change + /// The MNA variable component interface is used to recompute the system Matrix after each change. + class varResSwitch : + public MNASimPowerComp, + public Base::Ph3::Switch, + public MNAVariableCompInterface, + public MNASwitchInterface, + public SharedFactory { + + protected: + + Bool mPrevState=false; + Real mDeltaResClosed = 0; + Real mDeltaResOpen = 1.5; + Matrix mPrevRes; // previous resistance value to multiply with rate of change + // because we change the base value mClosedResistance & mOpenResistance to recompute the system Matrix + // we need to save the initialisation values to use them as target values in the transition + Matrix mInitClosedRes; + Matrix mInitOpenRes; + + + + public: + /// Defines UID, name, component parameters and logging level + varResSwitch(String uid, String name, Logger::Level logLevel = Logger::Level::off); + /// Defines name, component parameters and logging level + varResSwitch(String name, Logger::Level logLevel = Logger::Level::off) + : varResSwitch(name, name, logLevel) { } + + SimPowerComp::Ptr clone(String name); + + // #### General #### + void setInitParameters(Real timestep); + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency); + + // #### General MNA section #### + void mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector); + /// Stamps system matrix + void mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMatrix); + /// Stamps right side (source) vector + void mnaCompApplyRightSideVectorStamp(Matrix& rightVector); + /// Update interface voltage from MNA system result + void mnaCompUpdateVoltage(const Matrix& leftVector); + /// Update interface current from MNA system result + void mnaCompUpdateCurrent(const Matrix& leftVector); + + // #### MNA section for switches #### + /// Check if switch is closed + Bool mnaIsClosed() { return isClosed(); } + /// Stamps system matrix considering the defined switch position + void mnaCompApplySwitchSystemMatrixStamp(Bool closed, SparseMatrixRow& systemMatrix, Int freqIdx); + + void mnaCompPostStep(Real time, Int timeStepCount, Attribute::Ptr &leftVector) override; + + /// Add MNA post step dependencies + void mnaCompAddPostStepDependencies(AttributeBase::List &prevStepDependencies, AttributeBase::List &attributeDependencies, AttributeBase::List &modifiedAttributes, Attribute::Ptr &leftVector) override; + + Bool hasParameterChanged(); + }; +} +} +} diff --git a/dpsim-models/src/CMakeLists.txt b/dpsim-models/src/CMakeLists.txt index 79334c602c..cb11f1307c 100644 --- a/dpsim-models/src/CMakeLists.txt +++ b/dpsim-models/src/CMakeLists.txt @@ -98,6 +98,7 @@ list(APPEND MODELS_SOURCES EMT/EMT_Ph3_SynchronGeneratorIdeal.cpp # EMT/EMT_Ph3_SynchronGeneratorVBRStandalone.cpp EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp + EMT/EMT_Ph3_varResSwitch.cpp SP/SP_Ph1_VoltageSource.cpp SP/SP_Ph1_Capacitor.cpp diff --git a/dpsim-models/src/DP/DP_Ph1_SynchronGeneratorTrStab.cpp b/dpsim-models/src/DP/DP_Ph1_SynchronGeneratorTrStab.cpp index fa2b664ed4..87cdc15805 100644 --- a/dpsim-models/src/DP/DP_Ph1_SynchronGeneratorTrStab.cpp +++ b/dpsim-models/src/DP/DP_Ph1_SynchronGeneratorTrStab.cpp @@ -176,7 +176,7 @@ void DP::Ph1::SynchronGeneratorTrStab::initializeFromNodesAndTerminals(Real freq mSubVoltageSource->setVirtualNodeAt(mVirtualNodes[1], 0); mSubVoltageSource->initialize(mFrequencies); mSubVoltageSource->initializeFromNodesAndTerminals(frequency); - addMNASubComponent(mSubVoltageSource, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true); + addMNASubComponent(mSubVoltageSource, MNA_SUBCOMP_TASK_ORDER::TASK_AFTER_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true); // Create sub inductor as Xpd mSubInductor = DP::Ph1::Inductor::make(**mName + "_ind", mLogLevel); @@ -184,7 +184,7 @@ void DP::Ph1::SynchronGeneratorTrStab::initializeFromNodesAndTerminals(Real freq mSubInductor->connect({mVirtualNodes[0], terminal(0)->node()}); mSubInductor->initialize(mFrequencies); mSubInductor->initializeFromNodesAndTerminals(frequency); - addMNASubComponent(mSubInductor, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true); + addMNASubComponent(mSubInductor, MNA_SUBCOMP_TASK_ORDER::TASK_AFTER_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true); SPDLOG_LOGGER_INFO(mSLog, "\n--- Initialize according to powerflow ---" "\nTerminal 0 voltage: {:e}<{:e}" diff --git a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp index d32cd5435c..2d591a4fc8 100644 --- a/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp +++ b/dpsim-models/src/EMT/EMT_Ph3_SynchronGeneratorTrStab.cpp @@ -9,27 +9,6 @@ #include using namespace CPS; -Matrix EMT::Ph3::SynchronGeneratorTrStab::parkTransformPowerInvariant(Real theta, const Matrix &fabc) { - // Calculates fdq = Tdq * fabc - // Assumes that d-axis starts aligned with phase a - Matrix Tdq = getParkTransformMatrixPowerInvariant(theta); - Matrix dqvector = Tdq * fabc; - return dqvector; -} - - -Matrix EMT::Ph3::SynchronGeneratorTrStab::getParkTransformMatrixPowerInvariant(Real theta) { - // Return park matrix for theta - // Assumes that d-axis starts aligned with phase a - Matrix Tdq = Matrix::Zero(2, 3); - Real k = sqrt(2. / 3.); - Tdq << - k * cos(theta), k * cos(theta - 2. * M_PI / 3.), k * cos(theta + 2. * M_PI / 3.), - -k * sin(theta), -k * sin(theta - 2. * M_PI / 3.), -k * sin(theta + 2. * M_PI / 3.); - return Tdq; -} - - EMT::Ph3::SynchronGeneratorTrStab::SynchronGeneratorTrStab(String uid, String name, Logger::Level logLevel) : Base::SynchronGenerator(mAttributes), CompositePowerComp(uid, name, true, true, logLevel), mEp(mAttributes->create("Ep")), @@ -38,6 +17,7 @@ EMT::Ph3::SynchronGeneratorTrStab::SynchronGeneratorTrStab(String uid, String na mDelta_p(mAttributes->create("delta_r")), mRefOmega(mAttributes->createDynamic("w_ref")), mRefDelta(mAttributes->createDynamic("delta_ref")) { + mPhaseType = PhaseType::ABC; setVirtualNodeNumber(2); setTerminalNumber(1); **mIntfVoltage = Matrix::Zero(3,1); @@ -52,9 +32,34 @@ SimPowerComp::Ptr EMT::Ph3::SynchronGeneratorTrStab::clone(String name) { return copy; } +Matrix EMT::Ph3::SynchronGeneratorTrStab::parkTransformPowerInvariant(Real theta, const Matrix &fabc) { + // Calculates fdq = Tdq * fabc + // Assumes that d-axis starts aligned with phase a + Matrix Tdq = getParkTransformMatrixPowerInvariant(theta); + Matrix dqvector = Tdq * fabc; + return dqvector; +} + + +Matrix EMT::Ph3::SynchronGeneratorTrStab::getParkTransformMatrixPowerInvariant(Real theta) { + // Return park matrix for theta + // Assumes that d-axis starts aligned with phase a + Matrix Tdq = Matrix::Zero(2, 3); + // Real k = sqrt(2. / 3.); + Real k = 2. / 3.; // choice of k affects the DQ power eq. (see Eremia pp 26) + Tdq << + k * cos(theta), k * cos(theta - 2. * M_PI / 3.), k * cos(theta + 2. * M_PI / 3.), + -k * sin(theta), -k * sin(theta - 2. * M_PI / 3.), -k * sin(theta + 2. * M_PI / 3.); + return Tdq; +} + void EMT::Ph3::SynchronGeneratorTrStab::setFundamentalParametersPU(Real nomPower, Real nomVolt, Real nomFreq, Real Ll, Real Lmd, Real Llfd, Real inertia, Real D) { setBaseParameters(nomPower, nomVolt, nomFreq); + mSLog->info("\n--- Base Parameters ---" + "\nnomPower: {:f}" + "\nnomVolt: {:f}" + "\nnomFreq: {:f}", mNomPower, mNomVolt, mNomFreq); // Input is in per unit but all values are converted to absolute values. mParameterType = ParameterType::statorReferred; @@ -71,14 +76,24 @@ void EMT::Ph3::SynchronGeneratorTrStab::setFundamentalParametersPU(Real nomPower mXpd = mNomOmega * (**mLd - mLmd*mLmd / mLfd) * mBase_L; mLpd = (**mLd - mLmd*mLmd / mLfd) * mBase_L; - SPDLOG_LOGGER_INFO(mSLog, "\n--- Parameters ---" + //The units of D are per unit power divided by per unit speed deviation. + // D is transformed to an absolute value to obtain Kd, which will be used in the swing equation + mKd= D*mNomPower/mNomOmega; + + mSLog->info("\n--- Parameters ---" "\nimpedance: {:f}" - "\ninductance: {:f}", mXpd, mLpd); + "\ninductance: {:f}" + "\ninertia: {:f}" + "\ndamping: {:f}", mXpd, mLpd, **mInertia, mKd); } void EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersSI(Real nomPower, Real nomVolt, Real nomFreq, Int polePairNumber, Real Rs, Real Lpd, Real inertiaJ, Real Kd) { setBaseParameters(nomPower, nomVolt, nomFreq); + mSLog->info("\n--- Base Parameters ---" + "\nnomPower: {:f}" + "\nnomVolt: {:f}" + "\nnomFreq: {:f}", mNomPower, mNomVolt, mNomFreq); mParameterType = ParameterType::statorReferred; mStateType = StateType::statorReferred; @@ -90,14 +105,23 @@ void EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersSI(Real nomPower, R mXpd = mNomOmega * Lpd; mLpd = Lpd; - SPDLOG_LOGGER_INFO(mSLog, "\n--- Parameters ---" - "\nimpedance: {:f}" - "\ninductance: {:f}", mXpd, mLpd); + mSLog->info("\n--- Parameters ---" + "\nimpedance: {:f}" + "\ninductance: {:f}" + "\ninertia: {:f}" + "\ndamping: {:f}", mXpd, mLpd, **mInertia, mKd); } -void EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersPU(Real nomPower, Real nomVolt, Real nomFreq, Real Xpd, - Real inertia, Real Rs, Real D) { +void EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersPU(Real nomPower, Real nomVolt, Real nomFreq, + Real Xpd, Real inertia, Real Rs, Real D) { setBaseParameters(nomPower, nomVolt, nomFreq); + mSLog->info("\n--- Base Parameters ---" + "\nnomPower: {:f}" + "\nnomVolt: {:f}" + "\nnomFreq: {:f}", mNomPower, mNomVolt, mNomFreq); + + mSLog->info("\n--- Parameters Per-Unit ---" + "\n Xpd: {:f} [p.u.]", Xpd); // Input is in per unit but all values are converted to absolute values. mParameterType = ParameterType::statorReferred; @@ -114,9 +138,18 @@ void EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersPU(Real nomPower, R // D is transformed to an absolute value to obtain Kd, which will be used in the swing equation mKd= D*mNomPower/mNomOmega; - SPDLOG_LOGGER_INFO(mSLog, "\n--- Parameters ---" - "\nimpedance: {:f}" - "\ninductance: {:f}", mXpd, mLpd); + mSLog->info("\n--- Parameters ---" + "\nimpedance: {:f}" + "\ninductance: {:f}" + "\ninertia: {:f}" + "\ndamping: {:f}", mXpd, mLpd, **mInertia, mKd); +} + +void EMT::Ph3::SynchronGeneratorTrStab::setModelFlags(Bool convertWithOmegaMech) { + mConvertWithOmegaMech = convertWithOmegaMech; + + mSLog->info("\n--- Model flags ---" + "\nconvertWithOmegaMech: {:s}", std::to_string(mConvertWithOmegaMech)); } void EMT::Ph3::SynchronGeneratorTrStab::setInitialValues(Complex elecPower, Real mechPower) { @@ -140,8 +173,7 @@ void EMT::Ph3::SynchronGeneratorTrStab::initializeFromNodesAndTerminals(Real fre MatrixComp intfVoltageComplex = MatrixComp::Zero(3, 1); MatrixComp intfCurrentComplex = MatrixComp::Zero(3, 1); // // derive complex threephase initialization from single phase initial values (only valid for balanced systems) - // intfVoltageComplex(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); - intfVoltageComplex(0, 0) = initialSingleVoltage(0); + intfVoltageComplex(0, 0) = RMS3PH_TO_PEAK1PH * initialSingleVoltage(0); intfVoltageComplex(1, 0) = intfVoltageComplex(0, 0) * SHIFT_TO_PHASE_B; intfVoltageComplex(2, 0) = intfVoltageComplex(0, 0) * SHIFT_TO_PHASE_C; intfCurrentComplex(0, 0) = std::conj(-2./3.*mInitElecPower / intfVoltageComplex(0, 0)); @@ -162,25 +194,30 @@ void EMT::Ph3::SynchronGeneratorTrStab::initializeFromNodesAndTerminals(Real fre // Delta_p is the angular position of mEp with respect to the synchronously rotating reference **mDelta_p= Math::phase(**mEp); - // // Update active electrical power that is compared with the mechanical power - **mElecActivePower = ( 3./2. * intfVoltageComplex(0,0) * std::conj( - intfCurrentComplex(0,0)) ).real(); - // mElecActivePower = ( (mEp - (**mIntfVoltage)(0,0)) / mImpedance * (**mIntfVoltage)(0,0) ).real(); - // For infinite power bus - // mElecActivePower = (Math::abs(mEp) * Math::abs((**mIntfVoltage)(0,0)) / mXpd) * sin(mDelta_p); + mThetaN = **mDelta_p - PI / 2.; + + // // without DQ transformation + // **mElecActivePower = 3./2. *( intfVoltageComplex(0, 0) * std::conj( -intfCurrentComplex(0, 0)) ).real(); + // **mElecReactivePower = 3./2. *( intfVoltageComplex(0, 0) * std::conj( -intfCurrentComplex(0, 0)) ).imag(); + + // Transform interface quantities to synchronously rotating DQ reference frame + Matrix intfVoltageDQ = parkTransformPowerInvariant(mThetaN, **mIntfVoltage); + Matrix intfCurrentDQ = parkTransformPowerInvariant(mThetaN, **mIntfCurrent); + + // Electric active power from DQ quantities + **mElecActivePower = -3./2. * (intfVoltageDQ(0, 0)*intfCurrentDQ(0, 0) + intfVoltageDQ(1, 0)*intfCurrentDQ(1, 0)); //using DQ with k= 2. / 3. + // **mElecActivePower = -1. * (intfVoltageDQ(0, 0)*intfCurrentDQ(0, 0) + intfVoltageDQ(1, 0)*intfCurrentDQ(1, 0)); //using DQ with k=sqrt(2. / 3.) // Start in steady state so that electrical and mech. power are the same // because of the initial condition mOmMech = mNomOmega the damping factor is not considered at the initialisation **mMechPower = **mElecActivePower - mKd*(**mOmMech - mNomOmega); // Initialize node between X'd and Ep - mVirtualNodes[0]->setInitialVoltage(PEAK1PH_TO_RMS3PH* **mEp); - - MatrixComp vref = MatrixComp::Zero(3,1); - vref= CPS::Math::singlePhaseVariableToThreePhase(**mEp); + mVirtualNodes[0]->setInitialVoltage(PEAK1PH_TO_RMS3PH*(**mEp)); // Create sub voltage source for emf mSubVoltageSource = EMT::Ph3::VoltageSource::make(**mName + "_src", mLogLevel); - mSubVoltageSource->setParameters(vref,frequency); + mSubVoltageSource->setParameters(CPS::Math::singlePhaseVariableToThreePhase(PEAK1PH_TO_RMS3PH*(**mEp)), mNomFreq); mSubVoltageSource->connect({SimNode::GND, mVirtualNodes[0]}); mSubVoltageSource->setVirtualNodeAt(mVirtualNodes[1], 0); mSubVoltageSource->initialize(mFrequencies); @@ -190,55 +227,66 @@ void EMT::Ph3::SynchronGeneratorTrStab::initializeFromNodesAndTerminals(Real fre // Create sub inductor as Xpd mSubInductor = EMT::Ph3::Inductor::make(**mName + "_ind", mLogLevel); mSubInductor->setParameters(CPS::Math::singlePhaseParameterToThreePhase(mLpd)); - mSubInductor->connect({mVirtualNodes[0],terminal(0)->node()}); + mSubInductor->connect({mVirtualNodes[0], mTerminals[0]->node()}); mSubInductor->initialize(mFrequencies); mSubInductor->initializeFromNodesAndTerminals(frequency); addMNASubComponent(mSubInductor, MNA_SUBCOMP_TASK_ORDER::TASK_AFTER_PARENT, MNA_SUBCOMP_TASK_ORDER::TASK_BEFORE_PARENT, true); - SPDLOG_LOGGER_INFO(mSLog, "\n--- Initialize according to powerflow ---" - "\nTerminal 0 voltage: {:e}<{:e}" - "\nVoltage behind reactance: {:e}<{:e}" + mSLog->info("\n--- Initialize according to powerflow ---" + "\nTerminal Voltage (ABC): {:s}" + "\nTerminal Current (ABC): {:s}" + "\nVoltage behind reactance (Phasor): {:e}<{:e}" "\ninitial electrical power: {:e}+j{:e}" "\nactive electrical power: {:e}" "\nmechanical power: {:e}" "\n--- End of powerflow initialization ---", - Math::abs((**mIntfVoltage)(0,0)), Math::phaseDeg((**mIntfVoltage)(0,0)), + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), Math::abs(**mEp), Math::phaseDeg(**mEp), mInitElecPower.real(), mInitElecPower.imag(), **mElecActivePower, **mMechPower); } void EMT::Ph3::SynchronGeneratorTrStab::step(Real time) { - // #### Calculations on input of time step k #### + // // #### Calculations on input of time step k #### // Transform interface quantities to synchronously rotating DQ reference frame Matrix intfVoltageDQ = parkTransformPowerInvariant(mThetaN, **mIntfVoltage); Matrix intfCurrentDQ = parkTransformPowerInvariant(mThetaN, **mIntfCurrent); // Update electrical power (minus sign to calculate generated power from consumed current) - **mElecActivePower = - 1. * (intfVoltageDQ(0, 0)*intfCurrentDQ(0, 0) + intfVoltageDQ(1, 0)*intfCurrentDQ(1, 0)); - - // The damping factor Kd is adjusted to obtain a damping ratio of 0.3 - // Real MaxElecActivePower= Math::abs(mEp) * Math::abs((**mIntfVoltage)(0,0)) / mXpd; - // mKd=4*0.3*sqrt(mNomOmega*mInertia*MaxElecActivePower*0.5); - mKd=1*mNomPower; - - // #### Calculate state for time step k+1 #### - // semi-implicit Euler or symplectic Euler method for mechanical equations - Real dOmMech = mNomOmega / (2. * **mInertia * mNomPower) * (**mMechPower - **mElecActivePower-mKd*(**mOmMech - mNomOmega)); + **mElecActivePower = -3./2. * (intfVoltageDQ(0, 0)*intfCurrentDQ(0, 0) + intfVoltageDQ(1, 0)*intfCurrentDQ(1, 0)); //using DQ with k= 2. / 3. + // **mElecActivePower = -1. * (intfVoltageDQ(0, 0)*intfCurrentDQ(0, 0) + intfVoltageDQ(1, 0)*intfCurrentDQ(1, 0)); //using DQ with k=sqrt(2. / 3.) + + // Mechanical speed derivative at time step k + // convert torque to power with actual rotor angular velocity or nominal omega + Real dOmMech; + if (mConvertWithOmegaMech) + dOmMech = mNomOmega*mNomOmega / (2.* (**mInertia) * mNomPower * (**mOmMech)) * (**mMechPower - **mElecActivePower - mKd*(**mOmMech - mNomOmega)); + else + dOmMech = mNomOmega / (2. * **mInertia * mNomPower) * (**mMechPower - **mElecActivePower - mKd*(**mOmMech - mNomOmega)); + + // #### Calculate states for time step k+1 applying semi-implicit Euler #### + // Mechanical speed at time step k+1 applying Euler forward if (mBehaviour == Behaviour::MNASimulation) **mOmMech = **mOmMech + mTimeStep * dOmMech; - Real dDelta_p = **mOmMech - mNomOmega; - if (mBehaviour == Behaviour::MNASimulation) - **mDelta_p = **mDelta_p + mTimeStep * dDelta_p; + + // Derivative of rotor angle at time step k + 1 + // if reference omega is set, calculate delta with respect to reference + Real dDelta_p = **mOmMech - (mUseOmegaRef ? **mRefOmega : mNomOmega); + + // Rotor angle at time step k + 1 applying Euler backward // Update emf - only phase changes - if (mBehaviour == Behaviour::MNASimulation) + if (mBehaviour == Behaviour::MNASimulation) { + **mDelta_p = **mDelta_p + mTimeStep * dDelta_p; **mEp = Complex(**mEp_abs * cos(**mDelta_p), **mEp_abs * sin(**mDelta_p)); + } // Update nominal system angle - mThetaN = mThetaN + mTimeStep * mNomOmega; + // mThetaN = mThetaN + mTimeStep * mNomOmega; + mThetaN = mThetaN + mTimeStep * **mOmMech; - // mStates << Math::abs(mEp), Math::phaseDeg(mEp), mElecActivePower, mMechPower, - // mDelta_p, mOmMech, dOmMech, dDelta_p, (**mIntfVoltage)(0,0).real(), (**mIntfVoltage)(0,0).imag(); - // SPDLOG_LOGGER_DEBUG(mSLog, "\nStates, time {:f}: \n{:s}", time, Logger::matrixToString(mStates)); + mStates << Math::abs(**mEp), Math::phaseDeg(**mEp), **mElecActivePower, **mMechPower, + **mDelta_p, **mOmMech, dOmMech, dDelta_p, intfVoltageDQ(0, 0), intfVoltageDQ(1, 0) ; + SPDLOG_LOGGER_DEBUG(mSLog, "\nStates, time {:f}: \n{:s}", time, Logger::matrixToString(mStates)); } void EMT::Ph3::SynchronGeneratorTrStab::mnaParentInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { @@ -265,8 +313,8 @@ void EMT::Ph3::SynchronGeneratorTrStab::mnaParentPreStep(Real time, Int timeStep void EMT::Ph3::SynchronGeneratorTrStab::AddBStep::execute(Real time, Int timeStepCount) { **mGenerator.mRightVector = - **mGenerator.mSubInductor->mRightVector - + **mGenerator.mSubVoltageSource->mRightVector; + mGenerator.mSubInductor->mRightVector->get() + + mGenerator.mSubVoltageSource->mRightVector->get(); } void EMT::Ph3::SynchronGeneratorTrStab::mnaParentPostStep(Real time, Int timeStepCount, Attribute::Ptr &leftVector) { @@ -276,6 +324,11 @@ void EMT::Ph3::SynchronGeneratorTrStab::mnaParentPostStep(Real time, Int timeSte void EMT::Ph3::SynchronGeneratorTrStab::mnaCompUpdateVoltage(const Matrix& leftVector) { SPDLOG_LOGGER_DEBUG(mSLog, "Read voltage from {:d}", matrixNodeIndex(0)); + for (auto virtualNode : mVirtualNodes) + virtualNode->mnaUpdateVoltage(leftVector); + + **mIntfVoltage = Matrix::Zero(3,1); + (**mIntfVoltage)(0, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); (**mIntfVoltage)(1, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1)); (**mIntfVoltage)(2, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); @@ -284,5 +337,13 @@ void EMT::Ph3::SynchronGeneratorTrStab::mnaCompUpdateVoltage(const Matrix& leftV void EMT::Ph3::SynchronGeneratorTrStab::mnaCompUpdateCurrent(const Matrix& leftVector) { SPDLOG_LOGGER_DEBUG(mSLog, "Read current from {:d}", matrixNodeIndex(0)); - **mIntfCurrent = **mSubInductor->mIntfCurrent; + **mIntfCurrent = mSubInductor->mIntfCurrent->get(); +} + +void EMT::Ph3::SynchronGeneratorTrStab::setReferenceOmega(Attribute::Ptr refOmegaPtr, Attribute::Ptr refDeltaPtr) { + mRefOmega->setReference(refOmegaPtr); + mRefDelta->setReference(refDeltaPtr); + mUseOmegaRef=true; + + mSLog->info("Use of reference omega."); } diff --git a/dpsim-models/src/EMT/EMT_Ph3_varResSwitch.cpp b/dpsim-models/src/EMT/EMT_Ph3_varResSwitch.cpp new file mode 100644 index 0000000000..e705155841 --- /dev/null +++ b/dpsim-models/src/EMT/EMT_Ph3_varResSwitch.cpp @@ -0,0 +1,255 @@ +/* Copyright 2017-2020 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; + +// !!! TODO: Adaptions to use in EMT_Ph3 models phase-to-ground peak variables +// !!! with initialization from phase-to-phase RMS variables + +EMT::Ph3::varResSwitch::varResSwitch(String uid, String name, Logger::Level logLevel) + : MNASimPowerComp(uid, name, false, true, logLevel), Base::Ph3::Switch(mAttributes) { + setTerminalNumber(2); + **mIntfVoltage = Matrix::Zero(1,1); + **mIntfCurrent = Matrix::Zero(1,1); +} + +SimPowerComp::Ptr EMT::Ph3::varResSwitch::clone(String name) { + auto copy = varResSwitch::make(name, mLogLevel); + copy->setParameters(**mOpenResistance, **mClosedResistance, **mSwitchClosed); + return copy; +} + +void EMT::Ph3::varResSwitch::initializeFromNodesAndTerminals(Real frequency) { + + Matrix impedance = (**mSwitchClosed) ? **mClosedResistance : **mOpenResistance; + MatrixComp vInitABC = MatrixComp::Zero(3, 1); + vInitABC(0,0) = initialSingleVoltage(1) - initialSingleVoltage(0); + vInitABC(1, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_B; + vInitABC(2, 0) = vInitABC(0, 0) * SHIFT_TO_PHASE_C; + **mIntfVoltage = vInitABC.real(); + **mIntfCurrent = (impedance.inverse() * vInitABC).real(); + + SPDLOG_LOGGER_INFO(mSLog, + "\n--- Initialization from powerflow ---" + "\nVoltage across: {:s}" + "\nCurrent: {:s}" + "\nTerminal 0 voltage: {:s}" + "\nTerminal 1 voltage: {:s}" + "\n--- Initialization from powerflow finished ---", + Logger::matrixToString(**mIntfVoltage), + Logger::matrixToString(**mIntfCurrent), + Logger::phasorToString(initialSingleVoltage(0)), + Logger::phasorToString(initialSingleVoltage(1))); +} + +void EMT::Ph3::varResSwitch::mnaCompInitialize(Real omega, Real timeStep, Attribute::Ptr leftVector) { + + updateMatrixNodeIndices(); + **mRightVector = Matrix::Zero(0, 0); +} + +void EMT::Ph3::varResSwitch::mnaCompApplySystemMatrixStamp(SparseMatrixRow& systemMatrix) { + Matrix conductance = (**mSwitchClosed) ? + (**mClosedResistance).inverse() : (**mOpenResistance).inverse(); + + // Set diagonal entries + if (terminalNotGrounded(0)) { + // set upper left block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 0), conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 1), conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 2), conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 0), conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 1), conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 2), conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 0), conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 1), conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 2), conductance(2, 2)); + } + if (terminalNotGrounded(1)) { + // set buttom right block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 0), conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 1), conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 2), conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 0), conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 1), conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 2), conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 0), conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 1), conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 2), conductance(2, 2)); + } + // Set off diagonal blocks, 2x3x3 entries + if (terminalNotGrounded(0) && terminalNotGrounded(1)) { + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 0), -conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 1), -conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 2), -conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 0), -conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 1), -conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 2), -conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 0), -conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 1), -conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 2), -conductance(2, 2)); + + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 0), -conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 1), -conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 2), -conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 0), -conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 1), -conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 2), -conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 0), -conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 1), -conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 2), -conductance(2, 2)); + } + SPDLOG_LOGGER_TRACE(mSLog, + "\nConductance matrix: {:s}", + Logger::matrixToString(conductance)); +} + +void EMT::Ph3::varResSwitch::mnaCompApplySwitchSystemMatrixStamp(Bool closed, SparseMatrixRow& systemMatrix, Int freqIdx) { + MatrixFixedSize<3, 3> conductance = (closed) ? + (**mClosedResistance).inverse() : (**mOpenResistance).inverse(); + + // Set diagonal entries + if (terminalNotGrounded(0)) { + // set upper left block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 0), conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 1), conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(0, 2), conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 0), conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 1), conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(0, 2), conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 0), conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 1), conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(0, 2), conductance(2, 2)); + } + if (terminalNotGrounded(1)) { + // set buttom right block, 3x3 entries + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 0), conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 1), conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(1, 2), conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 0), conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 1), conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(1, 2), conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 0), conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 1), conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(1, 2), conductance(2, 2)); + } + // Set off diagonal blocks, 2x3x3 entries + if (terminalNotGrounded(0) && terminalNotGrounded(1)) { + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 0), -conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 1), -conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 0), matrixNodeIndex(1, 2), -conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 0), -conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 1), -conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 1), matrixNodeIndex(1, 2), -conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 0), -conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 1), -conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(0, 2), matrixNodeIndex(1, 2), -conductance(2, 2)); + + + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 0), -conductance(0, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 1), -conductance(0, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 0), matrixNodeIndex(0, 2), -conductance(0, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 0), -conductance(1, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 1), -conductance(1, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 1), matrixNodeIndex(0, 2), -conductance(1, 2)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 0), -conductance(2, 0)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 1), -conductance(2, 1)); + Math::addToMatrixElement(systemMatrix, matrixNodeIndex(1, 2), matrixNodeIndex(0, 2), -conductance(2, 2)); + } + + SPDLOG_LOGGER_TRACE(mSLog, + "\nConductance matrix: {:s}", + Logger::matrixToString(conductance)); +} + +void EMT::Ph3::varResSwitch::mnaCompApplyRightSideVectorStamp(Matrix& rightVector) { } + +void EMT::Ph3::varResSwitch::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::Ph3::varResSwitch::mnaCompPostStep(Real time, Int timeStepCount, Attribute::Ptr &leftVector) { + mnaUpdateVoltage(*leftVector); + mnaUpdateCurrent(*leftVector); +} + +void EMT::Ph3::varResSwitch::mnaCompUpdateVoltage(const Matrix& leftVector) { + // Voltage across component is defined as V1 - V0 + **mIntfVoltage = Matrix::Zero(3, 1); + if (terminalNotGrounded(1)) { + (**mIntfVoltage)(0, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 0)); + (**mIntfVoltage)(1, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 1)); + (**mIntfVoltage)(2, 0) = Math::realFromVectorElement(leftVector, matrixNodeIndex(1, 2)); + } + if (terminalNotGrounded(0)) { + (**mIntfVoltage)(0, 0) = (**mIntfVoltage)(0, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 0)); + (**mIntfVoltage)(1, 0) = (**mIntfVoltage)(1, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 1)); + (**mIntfVoltage)(2, 0) = (**mIntfVoltage)(2, 0) - Math::realFromVectorElement(leftVector, matrixNodeIndex(0, 2)); + } +} + +void EMT::Ph3::varResSwitch::mnaCompUpdateCurrent(const Matrix& leftVector) { + **mIntfCurrent = (**mSwitchClosed) ? + (**mClosedResistance).inverse() * **mIntfVoltage: + (**mOpenResistance).inverse() * **mIntfVoltage; +} + +Bool EMT::Ph3::varResSwitch::hasParameterChanged() { +//Get present state +Bool presentState=this->mnaIsClosed(); + +// Check if state of switch changed from open to closed +if (!(mPrevState == presentState)) { + // Switch is closed : change with 1/mDeltaRes + if (this->mnaIsClosed()==true) { + // **mClosedResistance= 1./mDeltaRes*mPrevRes; + **mClosedResistance= mDeltaResClosed*mPrevRes; + mPrevRes= **mClosedResistance; + // check if target value is reached + if ((**mClosedResistance)(0,0) < mInitClosedRes(0,0)) { + **mClosedResistance= mInitClosedRes; + mPrevRes= **mClosedResistance; + mPrevState= this->mnaIsClosed(); + } + } + // Switch is opened : change with mDeltaRes + else if (this->mnaIsClosed()==false) { + **mOpenResistance= mDeltaResOpen*mPrevRes; + mPrevRes= **mOpenResistance; + // check if target value is reached + if ( (**mOpenResistance)(0,0) > mInitOpenRes(0,0)) { + **mOpenResistance= mInitOpenRes; + mPrevRes= **mOpenResistance; + mPrevState= this->mnaIsClosed(); + } + } + return 1; //recompute system matrix +} +else{ + return 0; // do not recompute system matrix + } +} + +void EMT::Ph3::varResSwitch::setInitParameters(Real timestep) { + //Define variables for the transition + mDeltaResClosed= 0; + // mDeltaResOpen = 1.5; // assumption for 1ms step size + mDeltaResOpen= 0.5*timestep/0.001 + 1; + mPrevState= **mSwitchClosed; + mPrevRes= (**mSwitchClosed) ? **mClosedResistance : **mOpenResistance; + mInitClosedRes=**mClosedResistance; + mInitOpenRes=**mOpenResistance; +} \ No newline at end of file diff --git a/dpsim/examples/cxx/CMakeLists.txt b/dpsim/examples/cxx/CMakeLists.txt index dba9ec1117..1e2df519a1 100644 --- a/dpsim/examples/cxx/CMakeLists.txt +++ b/dpsim/examples/cxx/CMakeLists.txt @@ -86,8 +86,10 @@ set(CIRCUIT_SOURCES #3Bus System Circuits/SP_SynGenTrStab_3Bus_SteadyState.cpp Circuits/DP_SynGenTrStab_3Bus_SteadyState.cpp + Circuits/EMT_SynGenTrStab_3Bus_SteadyState.cpp Circuits/DP_SynGenTrStab_3Bus_Fault.cpp Circuits/SP_SynGenTrStab_3Bus_Fault.cpp + Circuits/EMT_SynGenTrStab_3Bus_Fault.cpp ) diff --git a/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_Fault.cpp b/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_Fault.cpp index 0cb3fcbc10..b2d7fe369c 100644 --- a/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_Fault.cpp +++ b/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_Fault.cpp @@ -16,10 +16,10 @@ using namespace CIM::Examples::Grids::ThreeBus; ScenarioConfig ThreeBus; //Switch to trigger fault at generator terminal -Real SwitchOpen = 1e12; -Real SwitchClosed = 0.1; +Real SwitchOpen = 1e6; +Real SwitchClosed = 1; -void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, bool startFaultEvent, bool endFaultEvent, Real startTimeFault, Real endTimeFault, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { +void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, Bool startFaultEvent, Bool endFaultEvent, Real startTimeFault, Real endTimeFault, Bool useVarResSwitch, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { // ----- POWERFLOW FOR INITIALIZATION ----- Real timeStepPF = finalTime; Real finalTimePF = finalTime+timeStepPF; @@ -32,13 +32,13 @@ void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b auto n3PF = SimNode::make("n3", PhaseType::Single); //Synchronous generator 1 - auto gen1PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); gen1PF->setBaseVoltage(ThreeBus.Vnom); //Synchronous generator 2 - auto gen2PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); gen2PF->setBaseVoltage(ThreeBus.Vnom); @@ -122,7 +122,6 @@ void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b // Get actual active and reactive power of generator's Terminal from Powerflow solution Complex initApparentPower_G2= gen2PF->getApparentPower(); gen2DP->setInitialValues(initApparentPower_G2, ThreeBus.initMechPower_G2); - gen2DP->setModelFlags(true); gen2DP->setReferenceOmega(gen1DP->attributeTyped("w_r"), gen1DP->attributeTyped("delta_r")); @@ -139,18 +138,18 @@ void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b //Line23 auto line23DP = DP::Ph1::PiLine::make("PiLine23", Logger::Level::debug); line23DP->setParameters(ThreeBus.lineResistance23, ThreeBus.lineInductance23, ThreeBus.lineCapacitance23, ThreeBus.lineConductance23); + + //two state Switch + auto faultDP = DP::Ph1::Switch::make("Br_fault", Logger::Level::debug); + faultDP->setParameters(SwitchOpen, SwitchClosed); + faultDP->open(); - // //Switch - // auto faultDP = DP::Ph1::Switch::make("Br_fault", Logger::Level::debug); + // //Variable resistance Switch + // auto faultDP = DP::Ph1::varResSwitch::make("Br_fault", Logger::Level::debug); // faultDP->setParameters(SwitchOpen, SwitchClosed); + // faultDP->setInitParameters(timeStep); // faultDP->open(); - //Variable resistance Switch - auto faultDP = DP::Ph1::varResSwitch::make("Br_fault", Logger::Level::debug); - faultDP->setParameters(SwitchOpen, SwitchClosed); - faultDP->setInitParameters(timeStep); - faultDP->open(); - // Topology gen1DP->connect({ n1DP }); gen2DP->connect({ n2DP }); @@ -204,9 +203,11 @@ void DP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b simDP.setFinalTime(finalTime); simDP.setDomain(Domain::DP); simDP.addLogger(loggerDP); - simDP.doSystemMatrixRecomputation(true); simDP.setDirectLinearSolverImplementation(DPsim::DirectLinearSolverImpl::SparseLU); + if (useVarResSwitch == true) { + simDP.doSystemMatrixRecomputation(true); + } // Events if (startFaultEvent){ auto sw1 = SwitchEvent::make(startTimeFault, faultDP, true); @@ -233,8 +234,9 @@ int main(int argc, char* argv[]) { Real timeStep = 0.001; Bool startFaultEvent=true; Bool endFaultEvent=true; + Bool useVarResSwitch=false; Real startTimeFault=10; - Real endTimeFault=10.2; + Real endTimeFault=10.1; Real cmdInertia_G1= 1.0; Real cmdInertia_G2= 1.0; Real cmdDamping_G1=1.0; @@ -258,7 +260,11 @@ int main(int argc, char* argv[]) { startTimeFault = args.getOptionReal("STARTTIMEFAULT"); if (args.options.find("ENDTIMEFAULT") != args.options.end()) endTimeFault = args.getOptionReal("ENDTIMEFAULT"); + // if (args.options.find("USEVARRESSWITCH") != args.options.end()) + // useVarResSwitch = args.options["USEVARRESSWITCH"]; + // if (args.options.find("FAULTRESISTANCE") != args.options.end()) + // SwitchClosed = args.options["FAULTRESISTANCE"]; } - - DP_SynGenTrStab_3Bus_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); + + DP_SynGenTrStab_3Bus_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, useVarResSwitch, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); } diff --git a/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_SteadyState.cpp b/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_SteadyState.cpp index eb68eb086c..fee7e453e6 100644 --- a/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_SteadyState.cpp +++ b/dpsim/examples/cxx/Circuits/DP_SynGenTrStab_3Bus_SteadyState.cpp @@ -28,13 +28,13 @@ void DP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT auto n3PF = SimNode::make("n3", PhaseType::Single); //Synchronous generator 1 - auto gen1PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); gen1PF->setBaseVoltage(ThreeBus.Vnom); //Synchronous generator 2 - auto gen2PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); gen2PF->setBaseVoltage(ThreeBus.Vnom); @@ -97,7 +97,7 @@ void DP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT // Components //Synchronous generator 1 - auto gen1DP = DP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen1DP = DP::Ph1::SynchronGeneratorTrStab::make("SynGen1", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen1DP->setStandardParametersPU(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.nomFreq_G1, ThreeBus.Xpd_G1*std::pow(ThreeBus.t1_ratio,2), cmdInertia_G1*ThreeBus.H_G1, ThreeBus.Rs_G1, cmdDamping_G1*ThreeBus.D_G1); // Get actual active and reactive power of generator's Terminal from Powerflow solution @@ -105,7 +105,7 @@ void DP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT gen1DP->setInitialValues(initApparentPower_G1, ThreeBus.initMechPower_G1); //Synchronous generator 2 - auto gen2DP = DP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen2DP = DP::Ph1::SynchronGeneratorTrStab::make("SynGen2", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen2DP->setStandardParametersPU(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.nomFreq_G2, ThreeBus.Xpd_G2*std::pow(ThreeBus.t2_ratio,2), cmdInertia_G2*ThreeBus.H_G2, ThreeBus.Rs_G2, cmdDamping_G2*ThreeBus.D_G2); // Get actual active and reactive power of generator's Terminal from Powerflow solution diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_Fault.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_Fault.cpp new file mode 100644 index 0000000000..96a3f29169 --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_Fault.cpp @@ -0,0 +1,285 @@ +/* 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 "../Examples.h" + +using namespace DPsim; +using namespace CPS; +using namespace CIM::Examples::Grids::ThreeBus; + +ScenarioConfig ThreeBus; + +//Switch to trigger fault at generator terminal +Real SwitchOpen = 1e6; +Real SwitchClosed = 1; + +void EMT_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, Bool startFaultEvent, Bool endFaultEvent, Real startTimeFault, Real endTimeFault, Bool useVarResSwitch, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { + // ----- POWERFLOW FOR INITIALIZATION ----- + Real timeStepPF = finalTime; + Real finalTimePF = finalTime+timeStepPF; + String simNamePF = simName + "_PF"; + Logger::setLogDir("logs/" + simNamePF); + + // Components + auto n1PF = SimNode::make("n1", PhaseType::Single); + auto n2PF = SimNode::make("n2", PhaseType::Single); + auto n3PF = SimNode::make("n3", PhaseType::Single); + + //Synchronous generator 1 + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); + // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side + gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); + gen1PF->setBaseVoltage(ThreeBus.Vnom); + + //Synchronous generator 2 + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); + // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side + gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); + gen2PF->setBaseVoltage(ThreeBus.Vnom); + + //use Shunt as Load for powerflow + auto loadPF = SP::Ph1::Shunt::make("Load", Logger::Level::debug); + loadPF->setParameters(ThreeBus.activePower_L / std::pow(ThreeBus.Vnom, 2), - ThreeBus.reactivePower_L / std::pow(ThreeBus.Vnom, 2)); + loadPF->setBaseVoltage(ThreeBus.Vnom); + + //Line12 + auto line12PF = SP::Ph1::PiLine::make("PiLine12", Logger::Level::debug); + line12PF->setParameters(ThreeBus.lineResistance12, ThreeBus.lineInductance12, ThreeBus.lineCapacitance12, ThreeBus.lineConductance12); + line12PF->setBaseVoltage(ThreeBus.Vnom); + //Line13 + auto line13PF = SP::Ph1::PiLine::make("PiLine13", Logger::Level::debug); + line13PF->setParameters(ThreeBus.lineResistance13, ThreeBus.lineInductance13, ThreeBus.lineCapacitance13, ThreeBus.lineConductance13); + line13PF->setBaseVoltage(ThreeBus.Vnom); + //Line23 + auto line23PF = SP::Ph1::PiLine::make("PiLine23", Logger::Level::debug); + line23PF->setParameters(ThreeBus.lineResistance23, ThreeBus.lineInductance23, ThreeBus.lineCapacitance23, ThreeBus.lineConductance23); + line23PF->setBaseVoltage(ThreeBus.Vnom); + //Switch + auto faultPF = CPS::SP::Ph1::Switch::make("Br_fault", Logger::Level::debug); + faultPF->setParameters(SwitchOpen, SwitchClosed); + faultPF->open(); + + // Topology + gen1PF->connect({ n1PF }); + gen2PF->connect({ n2PF }); + loadPF->connect({ n3PF }); + line12PF->connect({ n1PF, n2PF }); + line13PF->connect({ n1PF, n3PF }); + line23PF->connect({ n2PF, n3PF }); + // faultPF->connect({SP::SimNode::GND , n1PF }); //terminal of generator 1 + faultPF->connect({SP::SimNode::GND , n2PF }); //terminal of generator 2 + // faultPF->connect({SP::SimNode::GND , n3PF }); //Load bus + auto systemPF = SystemTopology(60, + SystemNodeList{n1PF, n2PF, n3PF}, + SystemComponentList{gen1PF, gen2PF, loadPF, line12PF, line13PF, line23PF, faultPF}); + + // Logging + auto loggerPF = DataLogger::make(simNamePF); + loggerPF->logAttribute("v_bus1", n1PF->attribute("v")); + loggerPF->logAttribute("v_bus2", n2PF->attribute("v")); + loggerPF->logAttribute("v_bus3", n3PF->attribute("v")); + + // Simulation + Simulation simPF(simNamePF, Logger::Level::debug); + simPF.setSystem(systemPF); + simPF.setTimeStep(timeStepPF); + simPF.setFinalTime(finalTimePF); + simPF.setDomain(Domain::SP); + simPF.setSolverType(Solver::Type::NRP); + simPF.doInitFromNodesAndTerminals(false); + simPF.addLogger(loggerPF); + simPF.run(); + + // ----- Dynamic simulation ------ + String simNameEMT = simName + "_EMT"; + Logger::setLogDir("logs/"+simNameEMT); + + // Nodes + auto n1EMT = SimNode::make("n1", PhaseType::ABC); + auto n2EMT = SimNode::make("n2", PhaseType::ABC); + auto n3EMT = SimNode::make("n3", PhaseType::ABC); + + // Components + //Synchronous generator 1 + auto gen1EMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen1", Logger::Level::debug); + // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side + gen1EMT->setStandardParametersPU(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.nomFreq_G1, ThreeBus.Xpd_G1*std::pow(ThreeBus.t1_ratio,2), cmdInertia_G1*ThreeBus.H_G1, ThreeBus.Rs_G1, cmdDamping_G1*ThreeBus.D_G1); + // Get actual active and reactive power of generator's Terminal from Powerflow solution + Complex initApparentPower_G1= gen1PF->getApparentPower(); + gen1EMT->setInitialValues(initApparentPower_G1, ThreeBus.initMechPower_G1); + + //Synchronous generator 1 + auto gen2EMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen2", Logger::Level::debug); + // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side + gen2EMT->setStandardParametersPU(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.nomFreq_G2, ThreeBus.Xpd_G2*std::pow(ThreeBus.t2_ratio,2), cmdInertia_G2*ThreeBus.H_G2, ThreeBus.Rs_G2, cmdDamping_G2*ThreeBus.D_G2); + // Get actual active and reactive power of generator's Terminal from Powerflow solution + Complex initApparentPower_G2= gen2PF->getApparentPower(); + gen2EMT->setInitialValues(initApparentPower_G2, ThreeBus.initMechPower_G2); + + gen2EMT->setModelFlags(true); + gen2EMT->setReferenceOmega(gen1EMT->attributeTyped("w_r"), gen1EMT->attributeTyped("delta_r")); + + //Load + auto loadEMT = EMT::Ph3::RXLoad::make("Load", Logger::Level::debug); + loadEMT->setParameters(CPS::Math::singlePhasePowerToThreePhase(ThreeBus.activePower_L), CPS::Math::singlePhasePowerToThreePhase(ThreeBus.reactivePower_L), ThreeBus.Vnom); + + // Line12 + auto line12EMT = EMT::Ph3::PiLine::make("PiLine12", Logger::Level::debug); + line12EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance12)); + + // Line13 + auto line13EMT = EMT::Ph3::PiLine::make("PiLine13", Logger::Level::debug); + line13EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance13)); + + // Line23 + auto line23EMT = EMT::Ph3::PiLine::make("PiLine23", Logger::Level::debug); + line23EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance23)); + + + // //Switch + auto faultEMT = CPS::EMT::Ph3::Switch::make("Br_fault", Logger::Level::debug); + faultEMT->setParameters(Math::singlePhaseParameterToThreePhase(SwitchOpen), + Math::singlePhaseParameterToThreePhase(SwitchClosed)); + faultEMT->openSwitch(); + + // // Variable resistance Switch + // auto faultEMT = EMT::Ph3::varResSwitch::make("Br_fault", Logger::Level::debug); + // faultEMT->setParameters(Math::singlePhaseParameterToThreePhase(SwitchOpen), + // Math::singlePhaseParameterToThreePhase(SwitchClosed)); + // faultEMT->setInitParameters(timeStep); + // faultEMT->openSwitch(); + +// Topology + gen1EMT->connect({ n1EMT }); + gen2EMT->connect({ n2EMT }); + loadEMT->connect({ n3EMT }); + line12EMT->connect({ n1EMT, n2EMT }); + line13EMT->connect({ n1EMT, n3EMT }); + line23EMT->connect({ n2EMT, n3EMT }); + // faultEMT->connect({EMT::SimNode::GND , n1EMT }); //terminal of generator 1 + faultEMT->connect({EMT::SimNode::GND , n2EMT }); //terminal of generator 2 + // faultEMT->connect({EMT::SimNode::GND , n3EMT }); //Load bus + auto systemEMT = SystemTopology(60, + SystemNodeList{n1EMT, n2EMT, n3EMT}, + SystemComponentList{gen1EMT, gen2EMT, loadEMT, line12EMT, line13EMT, line23EMT, faultEMT}); + + // Initialization of dynamic topology + systemEMT.initWithPowerflow(systemPF); + + + // Logging + auto loggerEMT = DataLogger::make(simNameEMT); + loggerEMT->logAttribute("v1", n1EMT->attribute("v")); + loggerEMT->logAttribute("v2", n2EMT->attribute("v")); + loggerEMT->logAttribute("v3", n3EMT->attribute("v")); + loggerEMT->logAttribute("v_line12", line12EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line12", line12EMT->attribute("i_intf")); + loggerEMT->logAttribute("v_line13", line13EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line13", line13EMT->attribute("i_intf")); + loggerEMT->logAttribute("v_line23", line23EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line23", line23EMT->attribute("i_intf")); + loggerEMT->logAttribute("Ep_gen1", gen1EMT->attribute("Ep_mag")); + loggerEMT->logAttribute("v_gen1", gen1EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_gen1", gen1EMT->attribute("i_intf")); + loggerEMT->logAttribute("wr_gen1", gen1EMT->attribute("w_r")); + loggerEMT->logAttribute("delta_gen1", gen1EMT->attribute("delta_r")); + loggerEMT->logAttribute("Ep_gen2", gen2EMT->attribute("Ep_mag")); + loggerEMT->logAttribute("v_gen2", gen2EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_gen2", gen2EMT->attribute("i_intf")); + loggerEMT->logAttribute("wr_gen2", gen2EMT->attribute("w_r")); + loggerEMT->logAttribute("wref_gen2", gen2EMT->attribute("w_ref")); + loggerEMT->logAttribute("delta_gen2", gen2EMT->attribute("delta_r")); + loggerEMT->logAttribute("v_load", loadEMT->attribute("v_intf")); + loggerEMT->logAttribute("i_load", loadEMT->attribute("i_intf")); + loggerEMT->logAttribute("P_mech1", gen1EMT->attribute("P_mech")); + loggerEMT->logAttribute("P_mech2", gen2EMT->attribute("P_mech")); + loggerEMT->logAttribute("P_elec1", gen1EMT->attribute("P_elec")); + loggerEMT->logAttribute("P_elec2", gen2EMT->attribute("P_elec")); + + Simulation simEMT(simNameEMT, Logger::Level::debug); + simEMT.setSystem(systemEMT); + simEMT.setTimeStep(timeStep); + simEMT.setFinalTime(finalTime); + simEMT.setDomain(Domain::EMT); + simEMT.addLogger(loggerEMT); + simEMT.setDirectLinearSolverImplementation(DPsim::DirectLinearSolverImpl::SparseLU); + + if (useVarResSwitch == true) { + simEMT.doSystemMatrixRecomputation(true); + } + + // Events + if (startFaultEvent){ + auto sw1 = SwitchEvent3Ph::make(startTimeFault, faultEMT, true); + + simEMT.addEvent(sw1); + } + + if(endFaultEvent){ + + auto sw2 = SwitchEvent3Ph::make(endTimeFault, faultEMT, false); + simEMT.addEvent(sw2); + + } + + simEMT.run(); +} + +int main(int argc, char* argv[]) { + + + //Simultion parameters + String simName="EMT_SynGenTrStab_3Bus_Fault"; + Real finalTime = 30; + Real timeStep = 0.001; + Bool startFaultEvent=true; + Bool endFaultEvent=true; + Bool useVarResSwitch=false; + Real startTimeFault=10; + Real endTimeFault=10.1; + Real cmdInertia_G1= 1.0; + Real cmdInertia_G2= 1.0; + Real cmdDamping_G1=1.0; + Real cmdDamping_G2=1.0; + + + CommandLineArgs args(argc, argv); + if (argc > 1) { + timeStep = args.timeStep; + finalTime = args.duration; + if (args.name != "dpsim") + simName = args.name; + if (args.options.find("SCALEINERTIA_G1") != args.options.end()) + cmdInertia_G1 = args.getOptionReal("SCALEINERTIA_G1"); + if (args.options.find("SCALEINERTIA_G2") != args.options.end()) + cmdInertia_G2 = args.getOptionReal("SCALEINERTIA_G2"); + if (args.options.find("SCALEDAMPING_G1") != args.options.end()) + cmdDamping_G1 = args.getOptionReal("SCALEDAMPING_G1"); + if (args.options.find("SCALEDAMPING_G2") != args.options.end()) + cmdDamping_G2 = args.getOptionReal("SCALEDAMPING_G2"); + if (args.options.find("STARTTIMEFAULT") != args.options.end()) + startTimeFault = args.getOptionReal("STARTTIMEFAULT"); + if (args.options.find("ENDTIMEFAULT") != args.options.end()) + endTimeFault = args.getOptionReal("ENDTIMEFAULT"); + // if (args.options.find("USEVARRESSWITCH") != args.options.end()) + // useVarResSwitch = args.options["USEVARRESSWITCH"]; + // if (args.options.find("FAULTRESISTANCE") != args.options.end()) + // SwitchClosed = args.options["FAULTRESISTANCE"]; + } + + EMT_SynGenTrStab_3Bus_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, useVarResSwitch, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); +} \ No newline at end of file diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_SteadyState.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_SteadyState.cpp new file mode 100644 index 0000000000..cecda85ca8 --- /dev/null +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_3Bus_SteadyState.cpp @@ -0,0 +1,222 @@ +/* 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 "../Examples.h" + +using namespace DPsim; +using namespace CPS; +using namespace CIM::Examples::Grids::ThreeBus; + +ScenarioConfig ThreeBus; + +void EMT_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalTime, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { + // ----- POWERFLOW FOR INITIALIZATION ----- + Real timeStepPF = finalTime; + Real finalTimePF = finalTime+timeStepPF; + String simNamePF = simName + "_PF"; + Logger::setLogDir("logs/" + simNamePF); + + // Components + auto n1PF = SimNode::make("n1", PhaseType::Single); + auto n2PF = SimNode::make("n2", PhaseType::Single); + auto n3PF = SimNode::make("n3", PhaseType::Single); + + //Synchronous generator 1 + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); + // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side + gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); + gen1PF->setBaseVoltage(ThreeBus.Vnom); + + //Synchronous generator 2 + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); + // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side + gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); + gen2PF->setBaseVoltage(ThreeBus.Vnom); + + //use Shunt as Load for powerflow + auto loadPF = SP::Ph1::Shunt::make("Load", Logger::Level::debug); + loadPF->setParameters(ThreeBus.activePower_L / std::pow(ThreeBus.Vnom, 2), - ThreeBus.reactivePower_L / std::pow(ThreeBus.Vnom, 2)); + loadPF->setBaseVoltage(ThreeBus.Vnom); + + //Line12 + auto line12PF = SP::Ph1::PiLine::make("PiLine12", Logger::Level::debug); + line12PF->setParameters(ThreeBus.lineResistance12, ThreeBus.lineInductance12, ThreeBus.lineCapacitance12, ThreeBus.lineConductance12); + line12PF->setBaseVoltage(ThreeBus.Vnom); + //Line13 + auto line13PF = SP::Ph1::PiLine::make("PiLine13", Logger::Level::debug); + line13PF->setParameters(ThreeBus.lineResistance13, ThreeBus.lineInductance13, ThreeBus.lineCapacitance13, ThreeBus.lineConductance13); + line13PF->setBaseVoltage(ThreeBus.Vnom); + //Line23 + auto line23PF = SP::Ph1::PiLine::make("PiLine23", Logger::Level::debug); + line23PF->setParameters(ThreeBus.lineResistance23, ThreeBus.lineInductance23, ThreeBus.lineCapacitance23, ThreeBus.lineConductance23); + line23PF->setBaseVoltage(ThreeBus.Vnom); + + // Topology + gen1PF->connect({ n1PF }); + gen2PF->connect({ n2PF }); + loadPF->connect({ n3PF }); + line12PF->connect({ n1PF, n2PF }); + line13PF->connect({ n1PF, n3PF }); + line23PF->connect({ n2PF, n3PF }); + auto systemPF = SystemTopology(60, + SystemNodeList{n1PF, n2PF, n3PF}, + SystemComponentList{gen1PF, gen2PF, loadPF, line12PF, line13PF, line23PF}); + + // Logging + auto loggerPF = DataLogger::make(simNamePF); + loggerPF->logAttribute("v_bus1", n1PF->attribute("v")); + loggerPF->logAttribute("v_bus2", n2PF->attribute("v")); + loggerPF->logAttribute("v_bus3", n3PF->attribute("v")); + + // Simulation + Simulation simPF(simNamePF, Logger::Level::debug); + simPF.setSystem(systemPF); + simPF.setTimeStep(timeStepPF); + simPF.setFinalTime(finalTimePF); + simPF.setDomain(Domain::SP); + simPF.setSolverType(Solver::Type::NRP); + simPF.setSolverAndComponentBehaviour(Solver::Behaviour::Initialization); + simPF.doInitFromNodesAndTerminals(false); + simPF.addLogger(loggerPF); + simPF.run(); + + // ----- Dynamic simulation ------ + String simNameEMT = simName + "_EMT"; + Logger::setLogDir("logs/"+simNameEMT); + + // Nodes + auto n1EMT = SimNode::make("n1", PhaseType::ABC); + auto n2EMT = SimNode::make("n2", PhaseType::ABC); + auto n3EMT = SimNode::make("n3", PhaseType::ABC); + + // Components + //Synchronous generator 1 + auto gen1EMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen1", Logger::Level::debug); + // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side + gen1EMT->setStandardParametersPU(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.nomFreq_G1, ThreeBus.Xpd_G1*std::pow(ThreeBus.t1_ratio,2), cmdInertia_G1*ThreeBus.H_G1, ThreeBus.Rs_G1, cmdDamping_G1*ThreeBus.D_G1); + // Get actual active and reactive power of generator's Terminal from Powerflow solution + Complex initApparentPower_G1= gen1PF->getApparentPower(); + gen1EMT->setInitialValues(initApparentPower_G1, ThreeBus.initMechPower_G1); + + //Synchronous generator 2 + auto gen2EMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen2", Logger::Level::debug); + // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side + gen2EMT->setStandardParametersPU(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.nomFreq_G2, ThreeBus.Xpd_G2*std::pow(ThreeBus.t2_ratio,2), cmdInertia_G2*ThreeBus.H_G2, ThreeBus.Rs_G2, cmdDamping_G2*ThreeBus.D_G2); + // Get actual active and reactive power of generator's Terminal from Powerflow solution + Complex initApparentPower_G2= gen2PF->getApparentPower(); + gen2EMT->setInitialValues(initApparentPower_G2, ThreeBus.initMechPower_G2); + + //Load + auto loadEMT = EMT::Ph3::RXLoad::make("Load", Logger::Level::debug); + loadEMT->setParameters(CPS::Math::singlePhasePowerToThreePhase(ThreeBus.activePower_L), CPS::Math::singlePhasePowerToThreePhase(ThreeBus.reactivePower_L), ThreeBus.Vnom); + + // Line12 + auto line12EMT = EMT::Ph3::PiLine::make("PiLine12", Logger::Level::debug); + line12EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance12), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance12)); + + // Line13 + auto line13EMT = EMT::Ph3::PiLine::make("PiLine13", Logger::Level::debug); + line13EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance13), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance13)); + + // Line23 + auto line23EMT = EMT::Ph3::PiLine::make("PiLine23", Logger::Level::debug); + line23EMT->setParameters(Math::singlePhaseParameterToThreePhase(ThreeBus.lineResistance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineInductance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineCapacitance23), + Math::singlePhaseParameterToThreePhase(ThreeBus.lineConductance23)); + +// Topology + gen1EMT->connect({ n1EMT }); + gen2EMT->connect({ n2EMT }); + loadEMT->connect({ n3EMT }); + line12EMT->connect({ n1EMT, n2EMT }); + line13EMT->connect({ n1EMT, n3EMT }); + line23EMT->connect({ n2EMT, n3EMT }); + auto systemEMT = SystemTopology(60, + SystemNodeList{n1EMT, n2EMT, n3EMT}, + SystemComponentList{gen1EMT, gen2EMT, loadEMT, line12EMT, line13EMT, line23EMT}); + + // Initialization of dynamic topology + systemEMT.initWithPowerflow(systemPF); + + // Logging + auto loggerEMT = DataLogger::make(simNameEMT); + loggerEMT->logAttribute("v1", n1EMT->attribute("v")); + loggerEMT->logAttribute("v2", n2EMT->attribute("v")); + loggerEMT->logAttribute("v3", n3EMT->attribute("v")); + loggerEMT->logAttribute("v_line12", line12EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line12", line12EMT->attribute("i_intf")); + loggerEMT->logAttribute("v_line13", line13EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line13", line13EMT->attribute("i_intf")); + loggerEMT->logAttribute("v_line23", line23EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_line23", line23EMT->attribute("i_intf")); + loggerEMT->logAttribute("Ep_gen1", gen1EMT->attribute("Ep_mag")); + loggerEMT->logAttribute("v_gen1", gen1EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_gen1", gen1EMT->attribute("i_intf")); + loggerEMT->logAttribute("wr_gen1", gen1EMT->attribute("w_r")); + loggerEMT->logAttribute("delta_gen1", gen1EMT->attribute("delta_r")); + loggerEMT->logAttribute("Ep_gen2", gen2EMT->attribute("Ep_mag")); + loggerEMT->logAttribute("v_gen2", gen2EMT->attribute("v_intf")); + loggerEMT->logAttribute("i_gen2", gen2EMT->attribute("i_intf")); + loggerEMT->logAttribute("wr_gen2", gen2EMT->attribute("w_r")); + loggerEMT->logAttribute("wref_gen2", gen2EMT->attribute("w_ref")); + loggerEMT->logAttribute("delta_gen2", gen2EMT->attribute("delta_r")); + loggerEMT->logAttribute("v_load", loadEMT->attribute("v_intf")); + loggerEMT->logAttribute("i_load", loadEMT->attribute("i_intf")); + loggerEMT->logAttribute("P_mech1", gen1EMT->attribute("P_mech")); + loggerEMT->logAttribute("P_mech2", gen2EMT->attribute("P_mech")); + loggerEMT->logAttribute("P_elec1", gen1EMT->attribute("P_elec")); + loggerEMT->logAttribute("P_elec2", gen2EMT->attribute("P_elec")); + + Simulation simEMT(simNameEMT, Logger::Level::debug); + simEMT.setSystem(systemEMT); + simEMT.setTimeStep(timeStep); + simEMT.setFinalTime(finalTime); + simEMT.setDomain(Domain::EMT); + simEMT.addLogger(loggerEMT); + + simEMT.run(); +} + +int main(int argc, char* argv[]) { + + + //Simultion parameters + String simName="EMT_SynGenTrStab_3Bus_SteadyState"; + Real finalTime = 30; + Real timeStep = 0.001; + Real cmdInertia_G1= 1.0; + Real cmdInertia_G2= 1.0; + Real cmdDamping_G1=1.0; + Real cmdDamping_G2=1.0; + + CommandLineArgs args(argc, argv); + if (argc > 1) { + timeStep = args.timeStep; + finalTime = args.duration; + if (args.name != "dpsim") + simName = args.name; + if (args.options.find("SCALEINERTIA_G1") != args.options.end()) + cmdInertia_G1 = args.getOptionReal("SCALEINERTIA_G1"); + if (args.options.find("SCALEINERTIA_G2") != args.options.end()) + cmdInertia_G2 = args.getOptionReal("SCALEINERTIA_G2"); + if (args.options.find("SCALEDAMPING_G1") != args.options.end()) + cmdDamping_G1 = args.getOptionReal("SCALEDAMPING_G1"); + if (args.options.find("SCALEDAMPING_G2") != args.options.end()) + cmdDamping_G2 = args.getOptionReal("SCALEDAMPING_G2"); + } + + EMT_SynGenTrStab_3Bus_SteadyState(simName, timeStep, finalTime, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); +} \ No newline at end of file diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_Fault.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_Fault.cpp index 228bdee887..e2ed52bacf 100644 --- a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_Fault.cpp +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_Fault.cpp @@ -7,48 +7,20 @@ *********************************************************************************/ #include +#include "../Examples.h" using namespace DPsim; using namespace CPS; +using namespace CIM::Examples::Grids::SMIB; -//-----------Power system-----------// -//Voltage level as Base Voltage -Real Vnom = 230e3; - -//-----------Generator-----------// -Real nomPower = 500e6; -Real nomPhPhVoltRMS = 22e3; -Real nomFreq = 60; -Real nomOmega= nomFreq* 2*PI; -Real H = 5; -Real Xpd=0.31; -Real Rs = 0.003*0; -Real D = 1*50; -// Initialization parameters -Real initMechPower= 300e6; -Real initActivePower = 300e6; -Real setPointVoltage=nomPhPhVoltRMS + 0.05*nomPhPhVoltRMS; - -//-----------Transformer-----------// -Real t_ratio=Vnom/nomPhPhVoltRMS; - -//PiLine parameters calculated from CIGRE Benchmark system -Real lineResistance = 6.7; -Real lineInductance = 47./nomOmega; -// Real lineCapacitance = 3.42e-4/nomOmega; -Real lineCapacitance = 0; -Real lineConductance =0; - -// Parameters for powerflow initialization -// Slack voltage: 1pu -Real Vslack = Vnom; +ScenarioConfig smib; //Switch to trigger fault at generator terminal -Real SwitchOpen = 3e3; +Real SwitchOpen = 1e6; Real SwitchClosed = 0.1; -void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, bool startFaultEvent, bool endFaultEvent, Real startTimeFault, Real endTimeFault, Real cmdInertia) { - // // ----- POWERFLOW FOR INITIALIZATION ----- +void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, bool startFaultEvent, bool endFaultEvent, Real startTimeFault, Real endTimeFault, Real cmdInertia, Real cmdDamping) { + // ----- POWERFLOW FOR INITIALIZATION ----- Real timeStepPF = finalTime; Real finalTimePF = finalTime+timeStepPF; String simNamePF = simName + "_PF"; @@ -61,28 +33,34 @@ void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, b //Synchronous generator ideal model auto genPF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side - genPF->setParameters(nomPower, nomPhPhVoltRMS, initActivePower, setPointVoltage*t_ratio, PowerflowBusType::PV); - genPF->setBaseVoltage(Vnom); + genPF->setParameters(smib.nomPower, smib.nomPhPhVoltRMS, smib.initActivePower, smib.setPointVoltage*smib.t_ratio, PowerflowBusType::PV); + genPF->setBaseVoltage(smib.Vnom); genPF->modifyPowerFlowBusType(PowerflowBusType::PV); //Grid bus as Slack auto extnetPF = SP::Ph1::NetworkInjection::make("Slack", Logger::Level::debug); - extnetPF->setParameters(Vslack); - extnetPF->setBaseVoltage(Vnom); + extnetPF->setParameters(smib.Vnom); + extnetPF->setBaseVoltage(smib.Vnom); extnetPF->modifyPowerFlowBusType(PowerflowBusType::VD); //Line auto linePF = SP::Ph1::PiLine::make("PiLine", Logger::Level::debug); - linePF->setParameters(lineResistance, lineInductance, lineCapacitance, lineConductance); - linePF->setBaseVoltage(Vnom); + linePF->setParameters(smib.lineResistance, smib.lineInductance, smib.lineCapacitance, smib.lineConductance); + linePF->setBaseVoltage(smib.Vnom); + + //Switch + auto faultPF = CPS::SP::Ph1::Switch::make("Br_fault", Logger::Level::debug); + faultPF->setParameters(SwitchOpen, SwitchClosed); + faultPF->open(); // Topology genPF->connect({ n1PF }); + faultPF->connect({SP::SimNode::GND , n1PF }); linePF->connect({ n1PF, n2PF }); extnetPF->connect({ n2PF }); auto systemPF = SystemTopology(60, SystemNodeList{n1PF, n2PF}, - SystemComponentList{genPF, linePF, extnetPF}); + SystemComponentList{genPF, linePF, extnetPF, faultPF}); // Logging auto loggerPF = DataLogger::make(simNamePF); @@ -112,31 +90,38 @@ void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, b // Components auto genEMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side - genEMT->setStandardParametersPU(nomPower, nomPhPhVoltRMS, nomFreq, Xpd*std::pow(t_ratio,2), cmdInertia*H, Rs, D ); - + genEMT->setStandardParametersPU(smib.nomPower, smib.nomPhPhVoltRMS, smib.nomFreq, smib.Xpd*std::pow(smib.t_ratio,2), cmdInertia*smib.H, smib.Rs, cmdDamping*smib.D ); // Get actual active and reactive power of generator's Terminal from Powerflow solution Complex initApparentPower= genPF->getApparentPower(); - genEMT->setInitialValues(initApparentPower, initMechPower); + genEMT->setInitialValues(initApparentPower, smib.initMechPower); //Grid bus as Slack auto extnetEMT = EMT::Ph3::NetworkInjection::make("Slack", Logger::Level::debug); + extnetEMT->setParameters(CPS::Math::singlePhaseVariableToThreePhase(smib.Vnom), 60); //frequency must be set to 60 otherwise the subvoltage source will set it to 50 per default // Line auto lineEMT = EMT::Ph3::PiLine::make("PiLine", Logger::Level::debug); - lineEMT->setParameters(Math::singlePhaseParameterToThreePhase(lineResistance), - Math::singlePhaseParameterToThreePhase(lineInductance), - Math::singlePhaseParameterToThreePhase(lineCapacitance), - Math::singlePhaseParameterToThreePhase(lineConductance)); - - // Fault Switch - auto faultEMT = CPS::EMT::Ph3::Switch::make("Br_fault", Logger::Level::debug); + lineEMT->setParameters(Math::singlePhaseParameterToThreePhase(smib.lineResistance), + Math::singlePhaseParameterToThreePhase(smib.lineInductance), + Math::singlePhaseParameterToThreePhase(smib.lineCapacitance), + Math::singlePhaseParameterToThreePhase(smib.lineConductance)); + + // //Switch + // auto faultEMT = CPS::EMT::Ph3::Switch::make("Br_fault", Logger::Level::debug); + // faultEMT->setParameters(Math::singlePhaseParameterToThreePhase(SwitchOpen), + // Math::singlePhaseParameterToThreePhase(SwitchClosed)); + // faultEMT->openSwitch(); + + // Variable resistance Switch + auto faultEMT = EMT::Ph3::varResSwitch::make("Br_fault", Logger::Level::debug); faultEMT->setParameters(Math::singlePhaseParameterToThreePhase(SwitchOpen), - Math::singlePhaseParameterToThreePhase(SwitchClosed)); + Math::singlePhaseParameterToThreePhase(SwitchClosed)); + faultEMT->setInitParameters(timeStep); faultEMT->openSwitch(); // Topology genEMT->connect({ n1EMT }); - faultEMT->connect({EMT::SimNode::GND, n1EMT}); + faultEMT->connect({EMT::SimNode::GND , n1EMT }); lineEMT->connect({ n1EMT, n2EMT }); extnetEMT->connect({ n2EMT }); auto systemEMT = SystemTopology(60, @@ -159,6 +144,8 @@ void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, b loggerEMT->logAttribute("delta_r_gen", genEMT->attribute("delta_r")); loggerEMT->logAttribute("P_elec", genEMT->attribute("P_elec")); loggerEMT->logAttribute("P_mech", genEMT->attribute("P_mech")); + //Switch + loggerEMT->logAttribute("i_fault", faultEMT->attribute("i_intf")); //line loggerEMT->logAttribute("v_line", lineEMT->attribute("v_intf")); loggerEMT->logAttribute("i_line", lineEMT->attribute("i_intf")); @@ -174,8 +161,10 @@ void EMT_1ph_SynGenTrStab_Fault(String simName, Real timeStep, Real finalTime, b simEMT.setFinalTime(finalTime); simEMT.setDomain(Domain::EMT); simEMT.addLogger(loggerEMT); + simEMT.doSystemMatrixRecomputation(true); //for varres switch + simEMT.setDirectLinearSolverImplementation(DPsim::DirectLinearSolverImpl::SparseLU); - // Events + // Events if (startFaultEvent){ auto sw1 = SwitchEvent3Ph::make(startTimeFault, faultEMT, true); @@ -197,13 +186,30 @@ int main(int argc, char* argv[]) { //Simultion parameters String simName="EMT_SynGenTrStab_SMIB_Fault"; - Real finalTime = 1; - Real timeStep = 0.00001; - Bool startFaultEvent=false; - Bool endFaultEvent=false; - Real startTimeFault=0.5; - Real endTimeFault=0.6; - Real cmdInertia= 1; - - EMT_1ph_SynGenTrStab_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, cmdInertia); + Real finalTime = 30; + Real timeStep = 0.001; + Bool startFaultEvent=true; + Bool endFaultEvent=true; + Real startTimeFault=10; + Real endTimeFault=10.2; + Real cmdInertia= 1.0; + Real cmdDamping=1.0; + + CommandLineArgs args(argc, argv); + if (argc > 1) { + timeStep = args.timeStep; + finalTime = args.duration; + if (args.name != "dpsim") + simName = args.name; + if (args.options.find("SCALEINERTIA") != args.options.end()) + cmdInertia = args.getOptionReal("SCALEINERTIA"); + if (args.options.find("SCALEDAMPING") != args.options.end()) + cmdDamping = args.getOptionReal("SCALEDAMPING"); + if (args.options.find("STARTTIMEFAULT") != args.options.end()) + startTimeFault = args.getOptionReal("STARTTIMEFAULT"); + if (args.options.find("ENDTIMEFAULT") != args.options.end()) + endTimeFault = args.getOptionReal("ENDTIMEFAULT"); + } + + EMT_1ph_SynGenTrStab_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, cmdInertia, cmdDamping); } diff --git a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_SteadyState.cpp b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_SteadyState.cpp index 2f5714e8b7..0548a8a4e2 100644 --- a/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_SteadyState.cpp +++ b/dpsim/examples/cxx/Circuits/EMT_SynGenTrStab_SMIB_SteadyState.cpp @@ -5,46 +5,17 @@ * 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 "../Examples.h" using namespace DPsim; using namespace CPS; +using namespace CIM::Examples::Grids::SMIB; + +ScenarioConfig smib; -//-----------Power system-----------// -//Voltage level as Base Voltage -Real Vnom = 230e3; - -//-----------Generator-----------// -Real nomPower = 500e6; -Real nomPhPhVoltRMS = 22e3; -Real nomFreq = 60; -Real nomOmega= nomFreq* 2*PI; -Real H = 5; -Real Xpd=0.31; -Real Rs = 0.003*0; -Real D = 1*50; -// Initialization parameters -Real initMechPower= 300e6; -Real initActivePower = 300e6; -Real setPointVoltage=nomPhPhVoltRMS + 0.05*nomPhPhVoltRMS; - -//-----------Transformer-----------// -Real t_ratio=Vnom/nomPhPhVoltRMS; - -//PiLine parameters calculated from CIGRE Benchmark system -Real lineResistance = 6.7; -Real lineInductance = 47./nomOmega; -// Real lineCapacitance = 3.42e-4/nomOmega; -Real lineCapacitance = 0; -Real lineConductance =0; - -// Parameters for powerflow initialization -// Slack voltage: 1pu -Real Vslack = Vnom; - -void EMT_1ph_SynGenTrStab_SteadyState(String simName, Real timeStep, Real finalTime, bool startFaultEvent, bool endFaultEvent, Real startTimeFault, Real endTimeFault, Real cmdInertia) { - // // ----- POWERFLOW FOR INITIALIZATION ----- +void EMT_1ph_SynGenTrStab_SteadyState(String simName, Real timeStep, Real finalTime, Real cmdInertia, Real cmdDamping) { + // ----- POWERFLOW FOR INITIALIZATION ----- Real timeStepPF = finalTime; Real finalTimePF = finalTime+timeStepPF; String simNamePF = simName + "_PF"; @@ -57,20 +28,20 @@ void EMT_1ph_SynGenTrStab_SteadyState(String simName, Real timeStep, Real finalT //Synchronous generator ideal model auto genPF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side - genPF->setParameters(nomPower, nomPhPhVoltRMS, initActivePower, setPointVoltage*t_ratio, PowerflowBusType::PV); - genPF->setBaseVoltage(Vnom); + genPF->setParameters(smib.nomPower, smib.nomPhPhVoltRMS, smib.initActivePower, smib.setPointVoltage*smib.t_ratio, PowerflowBusType::PV); + genPF->setBaseVoltage(smib.Vnom); genPF->modifyPowerFlowBusType(PowerflowBusType::PV); //Grid bus as Slack auto extnetPF = SP::Ph1::NetworkInjection::make("Slack", Logger::Level::debug); - extnetPF->setParameters(Vslack); - extnetPF->setBaseVoltage(Vnom); + extnetPF->setParameters(smib.Vnom); + extnetPF->setBaseVoltage(smib.Vnom); extnetPF->modifyPowerFlowBusType(PowerflowBusType::VD); //Line auto linePF = SP::Ph1::PiLine::make("PiLine", Logger::Level::debug); - linePF->setParameters(lineResistance, lineInductance, lineCapacitance, lineConductance); - linePF->setBaseVoltage(Vnom); + linePF->setParameters(smib.lineResistance, smib.lineInductance, smib.lineCapacitance, smib.lineConductance); + linePF->setBaseVoltage(smib.Vnom); // Topology genPF->connect({ n1PF }); @@ -108,21 +79,20 @@ void EMT_1ph_SynGenTrStab_SteadyState(String simName, Real timeStep, Real finalT // Components auto genEMT = EMT::Ph3::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side - genEMT->setStandardParametersPU(nomPower, nomPhPhVoltRMS, nomFreq, Xpd*std::pow(t_ratio,2), cmdInertia*H, Rs, D ); - + genEMT->setStandardParametersPU(smib.nomPower, smib.nomPhPhVoltRMS, smib.nomFreq, smib.Xpd*std::pow(smib.t_ratio,2), cmdInertia*smib.H, smib.Rs, cmdDamping*smib.D ); // Get actual active and reactive power of generator's Terminal from Powerflow solution Complex initApparentPower= genPF->getApparentPower(); - genEMT->setInitialValues(initApparentPower, initMechPower); + genEMT->setInitialValues(initApparentPower, smib.initMechPower); //Grid bus as Slack auto extnetEMT = EMT::Ph3::NetworkInjection::make("Slack", Logger::Level::debug); - + extnetEMT->setParameters(CPS::Math::singlePhaseVariableToThreePhase(smib.Vnom), 60); //frequency must be set to 60 otherwise the subvoltage source will set it to 50 per default // Line auto lineEMT = EMT::Ph3::PiLine::make("PiLine", Logger::Level::debug); - lineEMT->setParameters(Math::singlePhaseParameterToThreePhase(lineResistance), - Math::singlePhaseParameterToThreePhase(lineInductance), - Math::singlePhaseParameterToThreePhase(lineCapacitance), - Math::singlePhaseParameterToThreePhase(lineConductance)); + lineEMT->setParameters(Math::singlePhaseParameterToThreePhase(smib.lineResistance), + Math::singlePhaseParameterToThreePhase(smib.lineInductance), + Math::singlePhaseParameterToThreePhase(smib.lineCapacitance), + Math::singlePhaseParameterToThreePhase(smib.lineConductance)); // Topology genEMT->connect({ n1EMT }); @@ -172,13 +142,22 @@ int main(int argc, char* argv[]) { //Simultion parameters String simName="EMT_SynGenTrStab_SMIB_SteadyState"; - Real finalTime = 0.1; - Real timeStep = 0.000001; - Bool startFaultEvent=false; - Bool endFaultEvent=false; - Real startTimeFault=10; - Real endTimeFault=10.1; + Real finalTime = 20; + Real timeStep = 0.001; Real cmdInertia= 1.0; - - EMT_1ph_SynGenTrStab_SteadyState(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, cmdInertia); + Real cmdDamping=1.0; + + CommandLineArgs args(argc, argv); + if (argc > 1) { + timeStep = args.timeStep; + finalTime = args.duration; + if (args.name != "dpsim") + simName = args.name; + if (args.options.find("SCALEINERTIA") != args.options.end()) + cmdInertia = args.getOptionReal("SCALEINERTIA"); + if (args.options.find("SCALEDAMPING") != args.options.end()) + cmdDamping = args.getOptionReal("SCALEDAMPING"); + } + + EMT_1ph_SynGenTrStab_SteadyState(simName, timeStep, finalTime, cmdInertia, cmdDamping); } diff --git a/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_Fault.cpp b/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_Fault.cpp index b16c3611aa..68d4051cd6 100644 --- a/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_Fault.cpp +++ b/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_Fault.cpp @@ -16,10 +16,10 @@ using namespace CIM::Examples::Grids::ThreeBus; ScenarioConfig ThreeBus; //Switch to trigger fault at generator terminal -Real SwitchOpen = 1e12; -Real SwitchClosed = 0.1; +Real SwitchOpen = 1e6; +Real SwitchClosed = 1; -void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, bool startFaultEvent, bool endFaultEvent, Real startTimeFault, Real endTimeFault, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { +void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, Bool startFaultEvent, Bool endFaultEvent, Real startTimeFault, Real endTimeFault, Bool useVarResSwitch, Real cmdInertia_G1, Real cmdInertia_G2, Real cmdDamping_G1, Real cmdDamping_G2) { // ----- POWERFLOW FOR INITIALIZATION ----- Real timeStepPF = finalTime; Real finalTimePF = finalTime+timeStepPF; @@ -32,13 +32,13 @@ void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b auto n3PF = SimNode::make("n3", PhaseType::Single); //Synchronous generator 1 - auto gen1PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); gen1PF->setBaseVoltage(ThreeBus.Vnom); //Synchronous generator 2 - auto gen2PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); gen2PF->setBaseVoltage(ThreeBus.Vnom); @@ -108,7 +108,7 @@ void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b // Components //Synchronous generator 1 - auto gen1SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen1SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen1", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen1SP->setStandardParametersPU(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.nomFreq_G1, ThreeBus.Xpd_G1*std::pow(ThreeBus.t1_ratio,2), cmdInertia_G1*ThreeBus.H_G1, ThreeBus.Rs_G1, cmdDamping_G1*ThreeBus.D_G1); // Get actual active and reactive power of generator's Terminal from Powerflow solution @@ -116,7 +116,7 @@ void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b gen1SP->setInitialValues(initApparentPower_G1, ThreeBus.initMechPower_G1); //Synchronous generator 2 - auto gen2SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen2SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen2", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen2SP->setStandardParametersPU(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.nomFreq_G2, ThreeBus.Xpd_G2*std::pow(ThreeBus.t2_ratio,2), cmdInertia_G2*ThreeBus.H_G2, ThreeBus.Rs_G2, cmdDamping_G2*ThreeBus.D_G2); // Get actual active and reactive power of generator's Terminal from Powerflow solution @@ -140,17 +140,17 @@ void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b auto line23SP = SP::Ph1::PiLine::make("PiLine23", Logger::Level::debug); line23SP->setParameters(ThreeBus.lineResistance23, ThreeBus.lineInductance23, ThreeBus.lineCapacitance23, ThreeBus.lineConductance23); - // //Switch - // auto faultSP = SP::Ph1::Switch::make("Br_fault", Logger::Level::debug); - // faultSP->setParameters(SwitchOpen, SwitchClosed); - // faultSP->open(); - - //Variable resistance Switch - auto faultSP = SP::Ph1::varResSwitch::make("Br_fault", Logger::Level::debug); + //Switch + auto faultSP = SP::Ph1::Switch::make("Br_fault", Logger::Level::debug); faultSP->setParameters(SwitchOpen, SwitchClosed); - faultSP->setInitParameters(timeStep); faultSP->open(); + // //Variable resistance Switch + // auto faultSP = SP::Ph1::varResSwitch::make("Br_fault", Logger::Level::debug); + // faultSP->setParameters(SwitchOpen, SwitchClosed); + // faultSP->setInitParameters(timeStep); + // faultSP->open(); + // Topology gen1SP->connect({ n1SP }); gen2SP->connect({ n2SP }); @@ -204,10 +204,12 @@ void SP_SynGenTrStab_3Bus_Fault(String simName, Real timeStep, Real finalTime, b simSP.setFinalTime(finalTime); simSP.setDomain(Domain::SP); simSP.addLogger(loggerSP); - simSP.doSystemMatrixRecomputation(true); simSP.setDirectLinearSolverImplementation(DPsim::DirectLinearSolverImpl::SparseLU); + if (useVarResSwitch == true) { + simSP.doSystemMatrixRecomputation(true); + } - // Events + // Events if (startFaultEvent){ auto sw1 = SwitchEvent::make(startTimeFault, faultSP, true); @@ -233,8 +235,9 @@ int main(int argc, char* argv[]) { Real timeStep = 0.001; Bool startFaultEvent=true; Bool endFaultEvent=true; + Bool useVarResSwitch=false; Real startTimeFault=10; - Real endTimeFault=10.2; + Real endTimeFault=10.1; Real cmdInertia_G1= 1.0; Real cmdInertia_G2= 1.0; Real cmdDamping_G1=1.0; @@ -258,7 +261,11 @@ int main(int argc, char* argv[]) { startTimeFault = args.getOptionReal("STARTTIMEFAULT"); if (args.options.find("ENDTIMEFAULT") != args.options.end()) endTimeFault = args.getOptionReal("ENDTIMEFAULT"); + // if (args.options.find("USEVARRESSWITCH") != args.options.end()) + // useVarResSwitch = args.options["USEVARRESSWITCH"]; + // if (args.options.find("FAULTRESISTANCE") != args.options.end()) + // SwitchClosed = args.options["FAULTRESISTANCE"]; } - SP_SynGenTrStab_3Bus_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); + SP_SynGenTrStab_3Bus_Fault(simName, timeStep, finalTime, startFaultEvent, endFaultEvent, startTimeFault, endTimeFault, useVarResSwitch, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); } diff --git a/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_SteadyState.cpp b/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_SteadyState.cpp index 941f7337ac..a62cf8a379 100644 --- a/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_SteadyState.cpp +++ b/dpsim/examples/cxx/Circuits/SP_SynGenTrStab_3Bus_SteadyState.cpp @@ -28,13 +28,13 @@ void SP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT auto n3PF = SimNode::make("n3", PhaseType::Single); //Synchronous generator 1 - auto gen1PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen1PF = SP::Ph1::SynchronGenerator::make("SynGen1", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen1PF->setParameters(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.initActivePower_G1, ThreeBus.setPointVoltage_G1*ThreeBus.t1_ratio, PowerflowBusType::VD); gen1PF->setBaseVoltage(ThreeBus.Vnom); //Synchronous generator 2 - auto gen2PF = SP::Ph1::SynchronGenerator::make("Generator", Logger::Level::debug); + auto gen2PF = SP::Ph1::SynchronGenerator::make("SynGen2", Logger::Level::debug); // setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side gen2PF->setParameters(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.initActivePower_G2, ThreeBus.setPointVoltage_G2*ThreeBus.t2_ratio, PowerflowBusType::PV); gen2PF->setBaseVoltage(ThreeBus.Vnom); @@ -97,7 +97,7 @@ void SP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT // Components //Synchronous generator 1 - auto gen1SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen1SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen1", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen1SP->setStandardParametersPU(ThreeBus.nomPower_G1, ThreeBus.nomPhPhVoltRMS_G1, ThreeBus.nomFreq_G1, ThreeBus.Xpd_G1*std::pow(ThreeBus.t1_ratio,2), cmdInertia_G1*ThreeBus.H_G1, ThreeBus.Rs_G1, cmdDamping_G1*ThreeBus.D_G1); // Get actual active and reactive power of generator's Terminal from Powerflow solution @@ -105,7 +105,7 @@ void SP_SynGenTrStab_3Bus_SteadyState(String simName, Real timeStep, Real finalT gen1SP->setInitialValues(initApparentPower_G1, ThreeBus.initMechPower_G1); //Synchronous generator 2 - auto gen2SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen", Logger::Level::debug); + auto gen2SP = SP::Ph1::SynchronGeneratorTrStab::make("SynGen2", Logger::Level::debug); // Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side gen2SP->setStandardParametersPU(ThreeBus.nomPower_G2, ThreeBus.nomPhPhVoltRMS_G2, ThreeBus.nomFreq_G2, ThreeBus.Xpd_G2*std::pow(ThreeBus.t2_ratio,2), cmdInertia_G2*ThreeBus.H_G2, ThreeBus.Rs_G2, cmdDamping_G2*ThreeBus.D_G2); // Get actual active and reactive power of generator's Terminal from Powerflow solution @@ -191,7 +191,7 @@ int main(int argc, char* argv[]) { Real cmdDamping_G1=1.0; Real cmdDamping_G2=1.0; - CommandLineArgs args(argc, argv); + CommandLineArgs args(argc, argv); if (argc > 1) { timeStep = args.timeStep; finalTime = args.duration; @@ -208,4 +208,4 @@ int main(int argc, char* argv[]) { } SP_SynGenTrStab_3Bus_SteadyState(simName, timeStep, finalTime, cmdInertia_G1, cmdInertia_G2, cmdDamping_G1, cmdDamping_G2); -} +} \ No newline at end of file diff --git a/dpsim/examples/cxx/Examples.h b/dpsim/examples/cxx/Examples.h index ed7af6ae35..7cc9a7f23f 100644 --- a/dpsim/examples/cxx/Examples.h +++ b/dpsim/examples/cxx/Examples.h @@ -413,7 +413,7 @@ namespace ThreeBus { struct ScenarioConfig { //-----------Network-----------// - Real Vnom = 230e3; + Real Vnom = 250e3; Real nomFreq = 60; Real nomOmega= nomFreq* 2*PI; @@ -422,23 +422,23 @@ namespace ThreeBus { Real nomPhPhVoltRMS_G1 = 25e3; Real nomFreq_G1 = 60; Real H_G1 = 6; - Real Xpd_G1=0.3; //in p.u - Real Rs_G1 = 0.003*0; //in p.u - Real D_G1 = 1.5; //in p.u - // Initialization parameters + Real Xpd_G1=0.3 +0.15; //in p.u of generator base + transformer 1 reactance in p.u of Tr bases !!!check value!!! + Real Rs_G1 = 0.003*0; //in p.u of generator base + Real D_G1 = 1.5; //in p.u + // Initialization parameters Real initActivePower_G1 = 270e6; Real initMechPower_G1 = 270e6; - Real setPointVoltage_G1=nomPhPhVoltRMS_G1+0.05*nomPhPhVoltRMS_G1; - + Real setPointVoltage_G1=nomPhPhVoltRMS_G1; + //-----------Generator 2 (bus2)-----------// Real nomPower_G2 = 50e6; Real nomPhPhVoltRMS_G2 = 13.8e3; Real nomFreq_G2 = 60; Real H_G2 = 2; - Real Xpd_G2=0.1; //in p.u - Real Rs_G2 = 0.003*0; //in p.u - Real D_G2 =1.5; //in p.u - // Initialization parameters + Real Xpd_G2=0.3 + 0.12; //in p.u of generator base + transformer 2 reactance in p.u of Tr bases + Real Rs_G2 = 0.003*0; //in p.u of generator base + Real D_G2 =1; //in p.u + // Initialization parameters Real initActivePower_G2 = 45e6; Real initMechPower_G2 = 45e6; Real setPointVoltage_G2=nomPhPhVoltRMS_G2-0.05*nomPhPhVoltRMS_G2; @@ -452,23 +452,21 @@ namespace ThreeBus { Real reactivePower_L= 150e6; // -----------Transmission Lines-----------// - // CIGREHVAmerican (230 kV) - Grids::CIGREHVAmerican::LineParameters lineCIGREHV; //line 1-2 (180km) - Real lineResistance12 = lineCIGREHV.lineResistancePerKm*180; - Real lineInductance12 = lineCIGREHV.lineReactancePerKm/nomOmega*180; - Real lineCapacitance12 = lineCIGREHV.lineSusceptancePerKm/nomOmega*180; - Real lineConductance12 = lineCIGREHV.lineConductancePerKm*180; + Real lineResistance12 = 7.2; + Real lineInductance12 = 72./nomOmega; + Real lineCapacitance12 = 774*1e-6/nomOmega; + Real lineConductance12 = 0; //line 1-3 (150km) - Real lineResistance13 = lineCIGREHV.lineResistancePerKm*150; - Real lineInductance13 = lineCIGREHV.lineReactancePerKm/nomOmega*150; - Real lineCapacitance13 = lineCIGREHV.lineSusceptancePerKm/nomOmega*150; - Real lineConductance13 = lineCIGREHV.lineConductancePerKm*150; + Real lineResistance13 = 4; + Real lineInductance13 = 40./nomOmega; + Real lineCapacitance13 = 645*1e-6/nomOmega;; + Real lineConductance13 = 0; //line 2-3 (80km) - Real lineResistance23 = lineCIGREHV.lineResistancePerKm*80; - Real lineInductance23 = lineCIGREHV.lineReactancePerKm/nomOmega*80; - Real lineCapacitance23 = lineCIGREHV.lineSusceptancePerKm/nomOmega*80; - Real lineConductance23 = lineCIGREHV.lineConductancePerKm*80; + Real lineResistance23 = 3.2; + Real lineInductance23 = 32./nomOmega; + Real lineCapacitance23 = 344*1e-6/nomOmega;; + Real lineConductance23 = 0; }; } diff --git a/dpsim/src/Event.cpp b/dpsim/src/Event.cpp index 2b05c2bbd2..e1c3717f27 100644 --- a/dpsim/src/Event.cpp +++ b/dpsim/src/Event.cpp @@ -20,8 +20,13 @@ void EventQueue::handleEvents(Real currentTime) { while (!mEvents.empty()) { e = mEvents.top(); - // if current time larger or equal to event time, execute event - if ( currentTime > e->mTime || (e->mTime - currentTime) < 100e-9) { + // // if current time larger or equal to event time, execute event + // if ( currentTime > e->mTime || (e->mTime - currentTime) < 1e-12) { + + int currentTime_us = (int)((currentTime + 1e-7)*1e6); //current time in microseconds as integer + int eTime_us = (int)((e->mTime+ 1e-7)*1e6); //event time in microseconds as integer + // if current time in us equal to event time in us + if ( (currentTime_us == eTime_us)) { e->execute(); std::cout << std::scientific << currentTime << ": Handle event time" << std::endl; //std::cout << std::scientific << e->mTime << ": Original event time" << std::endl; diff --git a/dpsim/src/pybind/EMTComponents.cpp b/dpsim/src/pybind/EMTComponents.cpp index 6c524f0f4e..24fce6831c 100644 --- a/dpsim/src/pybind/EMTComponents.cpp +++ b/dpsim/src/pybind/EMTComponents.cpp @@ -121,7 +121,29 @@ void addEMTPh3Components(py::module_ mEMTPh3) { .def("open", &CPS::EMT::Ph3::Switch::openSwitch) .def("close", &CPS::EMT::Ph3::Switch::closeSwitch) .def("connect", &CPS::EMT::Ph3::Switch::connect); - + + py::class_, CPS::SimPowerComp, CPS::Base::Ph3::Switch>(mEMTPh3, "varResSwitch", py::multiple_inheritance()) + .def(py::init(), "name"_a, "loglevel"_a = CPS::Logger::Level::off) + .def("set_parameters", &CPS::EMT::Ph3::varResSwitch::setParameters, "open_resistance"_a, "closed_resistance"_a, "closed"_a = false) // cppcheck-suppress assignBoolToPointer + .def("open", &CPS::EMT::Ph3::varResSwitch::openSwitch) + .def("close", &CPS::EMT::Ph3::varResSwitch::closeSwitch) + .def("set_init_parameters", &CPS::EMT::Ph3::varResSwitch::setInitParameters, "time_step"_a) + .def("connect", &CPS::EMT::Ph3::varResSwitch::connect); + + py::class_, CPS::SimPowerComp>(mEMTPh3, "SynchronGeneratorTrStab", py::multiple_inheritance()) + .def(py::init(), "name"_a, "loglevel"_a = CPS::Logger::Level::off) + .def("set_standard_parameters_PU", &CPS::EMT::Ph3::SynchronGeneratorTrStab::setStandardParametersPU, + "nom_power"_a, "nom_volt"_a, "nom_freq"_a, "Xpd"_a, "inertia"_a, "Rs"_a=0, "D"_a=0) + .def("set_fundamental_parameters_PU", &CPS::EMT::Ph3::SynchronGeneratorTrStab::setFundamentalParametersPU, + "nom_power"_a, "nom_volt"_a, "nom_freq"_a, "Ll"_a, "Lmd"_a, "Llfd"_a, "H"_a, "D"_a = 0) + .def("set_initial_values", &CPS::EMT::Ph3::SynchronGeneratorTrStab::setInitialValues, "elec_power"_a, "mech_power"_a) + .def("connect", &CPS::EMT::Ph3::SynchronGeneratorTrStab::connect) + .def("set_model_flags", &CPS::EMT::Ph3::SynchronGeneratorTrStab::setModelFlags, "convert_with_omega_mech"_a) + .def("set_reference_omega", [](CPS::EMT::Ph3::SynchronGeneratorTrStab &gen, std::string refOmegaName, CPS::IdentifiedObject::Ptr refOmegaComp, + std::string refDeltaName, CPS::IdentifiedObject::Ptr refDeltaComp) { + gen.setReferenceOmega(refOmegaComp->attributeTyped(refOmegaName), refDeltaComp->attributeTyped(refDeltaName)); + }, "ref_omega_name"_a="w_r", "ref_omage_comp"_a, "ref_delta_name"_a="delta_r", "ref_delta_comp"_a); + py::class_, CPS::SimPowerComp>(mEMTPh3, "SynchronGeneratorDQTrapez", py::multiple_inheritance()) .def(py::init(), "name"_a, "loglevel"_a = CPS::Logger::Level::off) .def("set_parameters_operational_per_unit", &CPS::EMT::Ph3::SynchronGeneratorDQTrapez::setParametersOperationalPerUnit, diff --git a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_Fault.ipynb b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_Fault.ipynb similarity index 63% rename from examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_Fault.ipynb rename to examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_Fault.ipynb index a4aae2ed07..1a66a6e9bb 100644 --- a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_Fault.ipynb +++ b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_Fault.ipynb @@ -17,7 +17,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "import villas.dataprocessing.readtools as rt\n", @@ -35,7 +37,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## DP Simulation" + "## EMT Simulation" ] }, { @@ -48,14 +50,352 @@ { "cell_type": "code", "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "#-----------Power system-----------\n", + "#Voltage level as Base Voltage\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", + "\n", + "#-----------Generator 1 (bus1)-----------\n", + "#Machine parameters\n", + "nom_power_G1 = 300e6\n", + "nom_ph_ph_volt_RMS_G1=25e3\n", + "nom_freq_G1 = 60\n", + "H_G1 = 6\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", + "# Initialization parameters \n", + "init_active_power_G1 = 270e6\n", + "init_mech_power_G1 = 270e6\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", + "\n", + "#-----------Generator 2 (bus2)-----------\n", + "#Machine parameters\n", + "nom_power_G2 = 50e6\n", + "nom_ph_ph_volt_RMS_G2=13.8e3\n", + "nom_freq_G2 = 60\n", + "H_G2 = 2\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", + "# Initialization parameters \n", + "init_active_power_G2 = 45e6\n", + "init_mech_power_G2 = 45e6\n", + "set_point_voltage_G2 = nom_ph_ph_volt_RMS_G2 * 0.95\n", + "\n", + "#-----------Transformers-----------\n", + "t1_ratio = V_nom / nom_ph_ph_volt_RMS_G1\n", + "t2_ratio = V_nom / nom_ph_ph_volt_RMS_G2\n", + "\n", + "#-----------Load (bus3)-----------\n", + "active_power_L = 310e6\n", + "reactive_power_L = 150e6\n", + "\n", + "#-----------Transmission Lines-----------\n", + "#PiLine parameters\n", + "#line 1-2 (180km)\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", + "line_conductance12 = 0\n", + "#line 1-3 (150km)\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", + "line_conductance13 = 0\n", + "#line 2-3 (80km)\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", + "line_conductance23 = 0\n", + "\n", + "#Switch to trigger fault at generator terminal\n", + "switch_open = 1e12\n", + "switch_closed = 0.1\n", + "\n", + "#Simulation parameters\n", + "sim_name=\"EMT_SynGenTrStab_3Bus_Fault\"\n", + "final_time = 30\n", + "time_step = 0.001\n", + "start_fault_event=True\n", + "end_fault_event=True\n", + "start_time_fault=10\n", + "end_time_fault=10.1\n", + "cmd_inertia_G1 = 1\n", + "cmd_inertia_G2 = 1\n", + "cmd_damping_G1 = 1\n", + "cmd_damping_G2 = 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Powerflow for Initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "time_step_pf = final_time\n", + "final_time_pf = final_time + time_step_pf\n", + "sim_name_pf = sim_name + \"_PF\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", + "\n", + "#Components\n", + "n1_pf = dpsimpy.sp.SimNode(\"n1\", dpsimpy.PhaseType.Single)\n", + "n2_pf = dpsimpy.sp.SimNode(\"n2\", dpsimpy.PhaseType.Single)\n", + "n3_pf = dpsimpy.sp.SimNode(\"n3\", dpsimpy.PhaseType.Single)\n", + "\n", + "#Synchronous generator 1\n", + "gen1_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "#setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side\n", + "gen1_pf.set_parameters(nom_power_G1, nom_ph_ph_volt_RMS_G1, init_active_power_G1,\n", + " set_point_voltage_G1 * t1_ratio, dpsimpy.PowerflowBusType.VD)\n", + "gen1_pf.set_base_voltage(V_nom)\n", + "\n", + "#Synchronous generator 2\n", + "gen2_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "#setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side\n", + "gen2_pf.set_parameters(nom_power_G2, nom_ph_ph_volt_RMS_G2, init_active_power_G2,\n", + " set_point_voltage_G2 * t2_ratio, dpsimpy.PowerflowBusType.PV)\n", + "gen2_pf.set_base_voltage(V_nom)\n", + "\n", + "#use Shunt as Load for powerflow\n", + "load_pf = dpsimpy.sp.ph1.Shunt(\"Load\", dpsimpy.LogLevel.debug)\n", + "load_pf.set_parameters(active_power_L / V_nom**2, - reactive_power_L / V_nom**2)\n", + "load_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line12\n", + "line12_pf = dpsimpy.sp.ph1.PiLine(\"PiLine12\", dpsimpy.LogLevel.debug)\n", + "line12_pf.set_parameters(line_resistance12, line_inductance12, line_capacitance12, line_conductance12)\n", + "line12_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line13\n", + "line13_pf = dpsimpy.sp.ph1.PiLine(\"PiLine13\", dpsimpy.LogLevel.debug)\n", + "line13_pf.set_parameters(line_resistance13, line_inductance13, line_capacitance13, line_conductance13)\n", + "line13_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line23\n", + "line23_pf = dpsimpy.sp.ph1.PiLine(\"PiLine23\", dpsimpy.LogLevel.debug)\n", + "line23_pf.set_parameters(line_resistance23, line_inductance23, line_capacitance23, line_conductance23)\n", + "line23_pf.set_base_voltage(V_nom)\n", + "\n", + "# Topology\n", + "gen1_pf.connect([n1_pf])\n", + "gen2_pf.connect([n2_pf])\n", + "load_pf.connect([n3_pf])\n", + "line12_pf.connect([n1_pf, n2_pf])\n", + "line13_pf.connect([n1_pf, n3_pf])\n", + "line23_pf.connect([n2_pf, n3_pf])\n", + "system_pf = dpsimpy.SystemTopology(60, [n1_pf, n2_pf, n3_pf], [gen1_pf, gen2_pf, load_pf, line12_pf, line13_pf, line23_pf])\n", + "\n", + "# Logging\n", + "logger_pf = dpsimpy.Logger(sim_name_pf)\n", + "logger_pf.log_attribute(\"v_bus1\", \"v\", n1_pf)\n", + "logger_pf.log_attribute(\"v_bus2\", \"v\", n2_pf)\n", + "logger_pf.log_attribute(\"v_bus3\", \"v\", n3_pf)\n", + "\n", + "# Simulation\n", + "sim_pf = dpsimpy.Simulation(sim_name_pf, dpsimpy.LogLevel.debug)\n", + "sim_pf.set_system(system_pf)\n", + "sim_pf.set_time_step(time_step_pf)\n", + "sim_pf.set_final_time(final_time_pf)\n", + "sim_pf.set_domain(dpsimpy.Domain.SP)\n", + "sim_pf.set_solver(dpsimpy.Solver.NRP)\n", + "sim_pf.set_solver_component_behaviour(dpsimpy.SolverBehaviour.Initialization)\n", + "sim_pf.do_init_from_nodes_and_terminals(False)\n", + "sim_pf.add_logger(logger_pf)\n", + "sim_pf.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Dynamic Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sim_name_emt = sim_name + \"_EMT\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_emt)\n", + "\n", + "#Components\n", + "n1_emt = dpsimpy.emt.SimNode(\"n1\", dpsimpy.PhaseType.ABC)\n", + "n2_emt = dpsimpy.emt.SimNode(\"n2\", dpsimpy.PhaseType.ABC)\n", + "n3_emt = dpsimpy.emt.SimNode(\"n3\", dpsimpy.PhaseType.ABC)\n", + "\n", + "#Synchronous generator 1\n", + "gen1_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen1_emt.set_standard_parameters_PU(nom_power_G1, nom_ph_ph_volt_RMS_G1, nom_freq_G1,\n", + " Xpd_G1 * t1_ratio**2, cmd_inertia_G1 * H_G1, Rs_G1,\n", + " cmd_damping_G1 *Kd_G1)\n", + "init_apparent_power_G1 = gen1_pf.get_apparent_power()\n", + "gen1_emt.set_initial_values(init_apparent_power_G1, init_mech_power_G1)\n", + "\n", + "#Synchronous generator 2\n", + "gen2_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen2_emt.set_standard_parameters_PU(nom_power_G2, nom_ph_ph_volt_RMS_G2, nom_freq_G2,\n", + " Xpd_G2 * t2_ratio**2, cmd_inertia_G2 * H_G2, Rs_G2,\n", + " cmd_damping_G2 *Kd_G2)\n", + "init_apparent_power_G2 = gen2_pf.get_apparent_power()\n", + "gen2_emt.set_initial_values(init_apparent_power_G2, init_mech_power_G2)\n", + "\n", + "gen2_emt.set_model_flags(convert_with_omega_mech=True)\n", + "gen2_emt.set_reference_omega(\"w_r\", gen1_emt, \"delta_r\", gen1_emt)\n", + "\n", + "# Load\n", + "load_emt = dpsimpy.emt.ph3.RXLoad(\"Load\", dpsimpy.LogLevel.debug)\n", + "load_emt.set_parameters(dpsimpy.Math.single_phase_power_to_three_phase(active_power_L),\n", + " dpsimpy.Math.single_phase_power_to_three_phase(reactive_power_L), V_nom)\n", + "\n", + "#Line12\n", + "line12_emt = dpsimpy.emt.ph3.PiLine(\"PiLine12\", dpsimpy.LogLevel.debug)\n", + "line12_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance12),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance12),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance12),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance12))\n", + "\n", + "#Line13\n", + "line13_emt = dpsimpy.emt.ph3.PiLine(\"PiLine13\", dpsimpy.LogLevel.debug)\n", + "line13_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance13),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance13),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance13),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance13))\n", + "\n", + "#Line23\n", + "line23_emt = dpsimpy.emt.ph3.PiLine(\"PiLine23\", dpsimpy.LogLevel.debug)\n", + "line23_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance23),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance23),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance23),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance23))\n", + "\n", + "# Switch\n", + "fault_emt = dpsimpy.emt.ph3.varResSwitch(\"Br_fault\", dpsimpy.LogLevel.debug)\n", + "fault_emt.set_parameters(dpsimpy.Math.single_phase_parameter_to_three_phase(switch_open), dpsimpy.Math.single_phase_parameter_to_three_phase(switch_closed))\n", + "fault_emt.set_init_parameters(time_step)\n", + "fault_emt.open()\n", + "\n", + "# Topology\n", + "gen1_emt.connect([n1_emt])\n", + "gen2_emt.connect([n2_emt])\n", + "load_emt.connect([n3_emt])\n", + "line12_emt.connect([n1_emt, n2_emt])\n", + "line13_emt.connect([n1_emt, n3_emt])\n", + "line23_emt.connect([n2_emt, n3_emt])\n", + "fault_emt.connect([dpsimpy.emt.SimNode.gnd, n2_emt]) # terminal of generator 2\n", + "system_emt = dpsimpy.SystemTopology(60, [n1_emt, n2_emt, n3_emt], [gen1_emt, gen2_emt, load_emt, line12_emt, line13_emt, line23_emt, fault_emt])\n", + "\n", + "# Initialization of dynamic topology\n", + "system_emt.init_with_powerflow(system_pf)\n", + "\n", + "# Logging\n", + "logger_emt = dpsimpy.Logger(sim_name_emt)\n", + "logger_emt.log_attribute(\"v1\", \"v\", n1_emt)\n", + "logger_emt.log_attribute(\"v2\", \"v\", n2_emt)\n", + "logger_emt.log_attribute(\"v3\", \"v\", n3_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_line12\", \"v_intf\", line12_emt)\n", + "logger_emt.log_attribute(\"v_line13\", \"v_intf\", line13_emt)\n", + "logger_emt.log_attribute(\"v_line23\", \"v_intf\", line23_emt)\n", + "logger_emt.log_attribute(\"i_line12\", \"i_intf\", line12_emt)\n", + "logger_emt.log_attribute(\"i_line13\", \"i_intf\", line13_emt)\n", + "logger_emt.log_attribute(\"i_line23\", \"i_intf\", line23_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep_gen1\", \"Ep_mag\", gen1_emt)\n", + "logger_emt.log_attribute(\"v_gen1\", \"v_intf\", gen1_emt)\n", + "logger_emt.log_attribute(\"i_gen1\", \"i_intf\", gen1_emt)\n", + "logger_emt.log_attribute(\"wr_gen1\", \"w_r\", gen1_emt)\n", + "logger_emt.log_attribute(\"delta_gen1\", \"delta_r\", gen1_emt)\n", + "logger_emt.log_attribute(\"P_mech1\", \"P_mech\", gen1_emt)\n", + "logger_emt.log_attribute(\"P_elec1\", \"P_elec\", gen1_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep_gen2\", \"Ep_mag\", gen2_emt)\n", + "logger_emt.log_attribute(\"v_gen2\", \"v_intf\", gen2_emt)\n", + "logger_emt.log_attribute(\"i_gen2\", \"i_intf\", gen2_emt)\n", + "logger_emt.log_attribute(\"wr_gen2\", \"w_r\", gen2_emt)\n", + "logger_emt.log_attribute(\"delta_gen2\", \"delta_r\", gen2_emt)\n", + "logger_emt.log_attribute(\"wref_gen2\", \"w_ref\", gen2_emt)\n", + "logger_emt.log_attribute(\"P_mech2\", \"P_mech\", gen2_emt)\n", + "logger_emt.log_attribute(\"P_elec2\", \"P_elec\", gen2_emt)\n", + "\n", + "logger_emt.log_attribute(\"i_fault\", \"i_intf\", fault_emt)\n", + "logger_emt.log_attribute(\"v_load\", \"v_intf\", load_emt)\n", + "logger_emt.log_attribute(\"i_load\", \"i_intf\", load_emt)\n", + "\n", + "# Simulation\n", + "sim_emt = dpsimpy.Simulation(sim_name_emt, dpsimpy.LogLevel.debug)\n", + "sim_emt.set_system(system_emt)\n", + "sim_emt.set_time_step(time_step)\n", + "sim_emt.set_final_time(final_time)\n", + "sim_emt.set_domain(dpsimpy.Domain.EMT)\n", + "sim_emt.add_logger(logger_emt)\n", + "sim_emt.do_system_matrix_recomputation(True)\n", + "sim_emt.set_direct_solver_implementation(dpsimpy.DirectLinearSolverImpl.SparseLU)\n", + "\n", + "\n", + "# Events\n", + "if start_fault_event:\n", + " sw1 = dpsimpy.event.SwitchEvent3Ph(start_time_fault, fault_emt, True)\n", + " sim_emt.add_event(sw1)\n", + " \n", + "if end_fault_event:\n", + " sw2 = dpsimpy.event.SwitchEvent3Ph(end_time_fault, fault_emt, False)\n", + " sim_emt.add_event(sw2)\n", + "\n", + "sim_emt.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## DP Simulation" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "### Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "#-----------Power system-----------\n", "#Voltage level as Base Voltage\n", - "V_nom = 230e3\n", - "nom_freq = 60\n", - "nom_omega = nom_freq * 2 * np.pi\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", "\n", "#-----------Generator 1 (bus1)-----------\n", "#Machine parameters\n", @@ -63,13 +403,13 @@ "nom_ph_ph_volt_RMS_G1=25e3\n", "nom_freq_G1 = 60\n", "H_G1 = 6\n", - "Xpd_G1 = 0.3 # in p.u\n", - "Rs_G1 = 0\n", - "Kd_G1 = 1.5\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", "# Initialization parameters \n", "init_active_power_G1 = 270e6\n", "init_mech_power_G1 = 270e6\n", - "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1 * 1.05\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", "\n", "#-----------Generator 2 (bus2)-----------\n", "#Machine parameters\n", @@ -77,9 +417,9 @@ "nom_ph_ph_volt_RMS_G2=13.8e3\n", "nom_freq_G2 = 60\n", "H_G2 = 2\n", - "Xpd_G2 = 0.1 # in p.u\n", - "Rs_G2 = 0\n", - "Kd_G2 = 1.5\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", "# Initialization parameters \n", "init_active_power_G2 = 45e6\n", "init_mech_power_G2 = 45e6\n", @@ -96,19 +436,19 @@ "#-----------Transmission Lines-----------\n", "#PiLine parameters\n", "#line 1-2 (180km)\n", - "line_resistance12 = 1.27e-4 * 529 *180\n", - "line_inductance12 = 9.05e-4 * 529 / nom_omega *180\n", - "line_capacitance12 = (1.81e-3 / 529) / nom_omega *180\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", "line_conductance12 = 0\n", "#line 1-3 (150km)\n", - "line_resistance13 = 1.27e-4 * 529 *150\n", - "line_inductance13 = 9.05e-4 * 529 / nom_omega *150\n", - "line_capacitance13 = (1.81e-3 / 529) / nom_omega *150\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", "line_conductance13 = 0\n", "#line 2-3 (80km)\n", - "line_resistance23 = 1.27e-4 * 529 *80\n", - "line_inductance23 = 9.05e-4 * 529 / nom_omega *80\n", - "line_capacitance23 = (1.81e-3 / 529) / nom_omega *80\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", "line_conductance23 = 0\n", "\n", "#Switch to trigger fault at generator terminal\n", @@ -122,7 +462,7 @@ "start_fault_event=True\n", "end_fault_event=True\n", "start_time_fault=10\n", - "end_time_fault=10.2\n", + "end_time_fault=10.1\n", "cmd_inertia_G1 = 1\n", "cmd_inertia_G2 = 1\n", "cmd_damping_G1 = 1\n", @@ -139,7 +479,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "time_step_pf = final_time\n", @@ -230,7 +572,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "sim_name_dp = sim_name + \"_DP\"\n", @@ -370,14 +714,16 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "#-----------Power system-----------\n", "#Voltage level as Base Voltage\n", - "V_nom = 230e3\n", - "nom_freq = 60\n", - "nom_omega = nom_freq * 2 * np.pi\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", "\n", "#-----------Generator 1 (bus1)-----------\n", "#Machine parameters\n", @@ -385,13 +731,13 @@ "nom_ph_ph_volt_RMS_G1=25e3\n", "nom_freq_G1 = 60\n", "H_G1 = 6\n", - "Xpd_G1 = 0.3 # in p.u\n", - "Rs_G1 = 0\n", - "Kd_G1 = 1.5\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", "# Initialization parameters \n", "init_active_power_G1 = 270e6\n", "init_mech_power_G1 = 270e6\n", - "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1 * 1.05\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", "\n", "#-----------Generator 2 (bus2)-----------\n", "#Machine parameters\n", @@ -399,9 +745,9 @@ "nom_ph_ph_volt_RMS_G2=13.8e3\n", "nom_freq_G2 = 60\n", "H_G2 = 2\n", - "Xpd_G2 = 0.1 # in p.u\n", - "Rs_G2 = 0\n", - "Kd_G2 = 1.5\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", "# Initialization parameters \n", "init_active_power_G2 = 45e6\n", "init_mech_power_G2 = 45e6\n", @@ -418,19 +764,19 @@ "#-----------Transmission Lines-----------\n", "#PiLine parameters\n", "#line 1-2 (180km)\n", - "line_resistance12 = 1.27e-4 * 529 *180\n", - "line_inductance12 = 9.05e-4 * 529 / nom_omega *180\n", - "line_capacitance12 = (1.81e-3 / 529) / nom_omega *180\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", "line_conductance12 = 0\n", "#line 1-3 (150km)\n", - "line_resistance13 = 1.27e-4 * 529 *150\n", - "line_inductance13 = 9.05e-4 * 529 / nom_omega *150\n", - "line_capacitance13 = (1.81e-3 / 529) / nom_omega *150\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", "line_conductance13 = 0\n", "#line 2-3 (80km)\n", - "line_resistance23 = 1.27e-4 * 529 *80\n", - "line_inductance23 = 9.05e-4 * 529 / nom_omega *80\n", - "line_capacitance23 = (1.81e-3 / 529) / nom_omega *80\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", "line_conductance23 = 0\n", "\n", "#Switch to trigger fault at generator terminal\n", @@ -444,7 +790,7 @@ "start_fault_event=True\n", "end_fault_event=True\n", "start_time_fault=10\n", - "end_time_fault=10.2\n", + "end_time_fault=10.1\n", "cmd_inertia_G1 = 1\n", "cmd_inertia_G2 = 1\n", "cmd_damping_G1 = 1\n", @@ -461,7 +807,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "time_step_pf = final_time\n", @@ -552,7 +900,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "sim_name_sp = sim_name + \"_SP\"\n", @@ -685,7 +1035,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "work_dir = 'logs/SP_SynGenTrStab_3Bus_Fault_SP/'\n", @@ -704,7 +1056,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "work_dir = 'logs/DP_SynGenTrStab_3Bus_Fault_DP/'\n", @@ -713,6 +1067,29 @@ "ts_dp1ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" ] }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "## EMT Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "work_dir = 'logs/EMT_SynGenTrStab_3Bus_Fault_EMT/'\n", + "log_name = 'EMT_SynGenTrStab_3Bus_Fault_EMT'\n", + "print(work_dir + log_name + '.csv')\n", + "ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -723,7 +1100,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "timestep=50e-6;\n", @@ -741,6 +1120,7 @@ "for name in ['v_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -749,9 +1129,8 @@ "for name in ['v_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend()" ] }, { @@ -764,7 +1143,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", @@ -775,6 +1156,7 @@ "for name in ['i_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -783,9 +1165,8 @@ "for name in ['i_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend()" ] }, { @@ -798,7 +1179,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", @@ -808,11 +1191,12 @@ "plt.ylabel('Generator 1 Rotor frequency (Hz)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(55,65)\n", + "plt.ylim(55,65)\n", "\n", "for name in ['wr_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='DP', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -820,14 +1204,13 @@ "plt.ylabel('Generator 2 Rotor frequency (Hz)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(55,65)\n", + "plt.ylim(55,65)\n", "\n", "for name in ['wr_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='DP', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", + "plt.legend()" ] }, { @@ -840,7 +1223,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", @@ -850,11 +1235,12 @@ "plt.ylabel('Generator 1 Rotor angular velocity (1/s)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(360,400)\n", + "plt.ylim(360,400)\n", "\n", "for name in ['wr_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -862,14 +1248,13 @@ "plt.ylabel('Generator 2 Rotor angular velocity (1/s)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(360,400)\n", + "plt.ylim(360,400)\n", "\n", "for name in ['wr_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", + "plt.legend()" ] }, { @@ -882,44 +1267,27 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", "\n", "f.add_subplot(1, 2, 1)\n", "plt.xlabel('time (s)', fontsize=20)\n", - "plt.ylabel('Generator 1 Rotor angle (degree)', fontsize=20)\n", + "plt.ylabel('$\\delta_2$ - $\\delta_1$ (degree)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(-20,20)\n", - "\n", - "for name in ['delta_gen1']:\n", - " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", - " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", "\n", - "f.add_subplot(1, 2, 2)\n", - "plt.xlabel('time (s)', fontsize=20)\n", - "plt.ylabel('Generator 2 Rotor angle (degree)', fontsize=20)\n", - "plt.xticks(fontsize=18)\n", - "plt.yticks(fontsize=18)\n", - "#plt.ylim(-20,20)\n", + "# plt.ylim(-50,50)\n", "\n", - "for name in ['delta_gen2']:\n", - " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", - " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" + "for name in ['delta_gen']:\n", + " plt.plot(ts_sp1ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_sp1ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", + " plt.plot(ts_dp1ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_dp1ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_emt3ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' EMT')\n", + "plt.legend()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -941,7 +1309,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.2" + "version": "3.9.13" }, "tests": { "skip": false diff --git a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb similarity index 61% rename from examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb rename to examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb index 7083c9a790..62ece97bce 100644 --- a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb +++ b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_3Bus_SteadyState.ipynb @@ -2,14 +2,16 @@ "cells": [ { "cell_type": "markdown", + "metadata": {}, "source": [ - "# 3 Bus (2 Generators, 1 Load) Steady state SP vs DP" - ], - "metadata": {} + "# 3 Bus (2 Generators, 1 Load) Steady state EMT vs DP vs SP" + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "import villas.dataprocessing.readtools as rt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", @@ -20,31 +22,331 @@ "import dpsimpy\n", "\n", "#%matplotlib widget" - ], + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## EMT Simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#-----------Power system-----------\n", + "#Voltage level as Base Voltage\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", + "\n", + "#-----------Generator 1 (bus1)-----------\n", + "#Machine parameters\n", + "nom_power_G1 = 300e6\n", + "nom_ph_ph_volt_RMS_G1=25e3\n", + "nom_freq_G1 = 60\n", + "H_G1 = 6\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", + "# Initialization parameters \n", + "init_active_power_G1 = 270e6\n", + "init_mech_power_G1 = 270e6\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", + "\n", + "#-----------Generator 2 (bus2)-----------\n", + "#Machine parameters\n", + "nom_power_G2 = 50e6\n", + "nom_ph_ph_volt_RMS_G2=13.8e3\n", + "nom_freq_G2 = 60\n", + "H_G2 = 2\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", + "# Initialization parameters \n", + "init_active_power_G2 = 45e6\n", + "init_mech_power_G2 = 45e6\n", + "set_point_voltage_G2 = nom_ph_ph_volt_RMS_G2 * 0.95\n", + "\n", + "#-----------Transformers-----------\n", + "t1_ratio = V_nom / nom_ph_ph_volt_RMS_G1\n", + "t2_ratio = V_nom / nom_ph_ph_volt_RMS_G2\n", + "\n", + "#-----------Load (bus3)-----------\n", + "active_power_L = 310e6\n", + "reactive_power_L = 150e6\n", + "\n", + "#-----------Transmission Lines-----------\n", + "#PiLine parameters\n", + "#line 1-2 (180km)\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", + "line_conductance12 = 0\n", + "#line 1-3 (150km)\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", + "line_conductance13 = 0\n", + "#line 2-3 (80km)\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", + "line_conductance23 = 0\n", + "\n", + "#Switch to trigger fault at generator terminal\n", + "switch_open = 1e12\n", + "switch_closed = 0.1\n", + "\n", + "#Simulation parameters\n", + "sim_name='EMT_SynGenTrStab_3Bus_SteadyState'\n", + "final_time = 10\n", + "time_step = 0.001\n", + "scale_inertia= 1.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Powerflow for Initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time_step_pf = final_time\n", + "final_time_pf = final_time + time_step_pf\n", + "sim_name_pf = sim_name + \"_PF\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", + "\n", + "#Components\n", + "n1_pf = dpsimpy.sp.SimNode(\"n1\", dpsimpy.PhaseType.Single)\n", + "n2_pf = dpsimpy.sp.SimNode(\"n2\", dpsimpy.PhaseType.Single)\n", + "n3_pf = dpsimpy.sp.SimNode(\"n3\", dpsimpy.PhaseType.Single)\n", + "\n", + "#Synchronous generator 1\n", + "gen1_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "#setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side\n", + "gen1_pf.set_parameters(nom_power_G1, nom_ph_ph_volt_RMS_G1, init_active_power_G1,\n", + " set_point_voltage_G1 * t1_ratio, dpsimpy.PowerflowBusType.VD)\n", + "gen1_pf.set_base_voltage(V_nom)\n", + "\n", + "#Synchronous generator 2\n", + "gen2_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "#setPointVoltage is defined as the voltage at the transfomer primary side and should be transformed to network side\n", + "gen2_pf.set_parameters(nom_power_G2, nom_ph_ph_volt_RMS_G2, init_active_power_G2,\n", + " set_point_voltage_G2 * t2_ratio, dpsimpy.PowerflowBusType.PV)\n", + "gen2_pf.set_base_voltage(V_nom)\n", + "\n", + "#use Shunt as Load for powerflow\n", + "load_pf = dpsimpy.sp.ph1.Shunt(\"Load\", dpsimpy.LogLevel.debug)\n", + "load_pf.set_parameters(active_power_L / V_nom**2, - reactive_power_L / V_nom**2)\n", + "load_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line12\n", + "line12_pf = dpsimpy.sp.ph1.PiLine(\"PiLine12\", dpsimpy.LogLevel.debug)\n", + "line12_pf.set_parameters(line_resistance12, line_inductance12, line_capacitance12, line_conductance12)\n", + "line12_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line13\n", + "line13_pf = dpsimpy.sp.ph1.PiLine(\"PiLine13\", dpsimpy.LogLevel.debug)\n", + "line13_pf.set_parameters(line_resistance13, line_inductance13, line_capacitance13, line_conductance13)\n", + "line13_pf.set_base_voltage(V_nom)\n", + "\n", + "#Line23\n", + "line23_pf = dpsimpy.sp.ph1.PiLine(\"PiLine23\", dpsimpy.LogLevel.debug)\n", + "line23_pf.set_parameters(line_resistance23, line_inductance23, line_capacitance23, line_conductance23)\n", + "line23_pf.set_base_voltage(V_nom)\n", + "\n", + "# Topology\n", + "gen1_pf.connect([n1_pf])\n", + "gen2_pf.connect([n2_pf])\n", + "load_pf.connect([n3_pf])\n", + "line12_pf.connect([n1_pf, n2_pf])\n", + "line13_pf.connect([n1_pf, n3_pf])\n", + "line23_pf.connect([n2_pf, n3_pf])\n", + "system_pf = dpsimpy.SystemTopology(60, [n1_pf, n2_pf, n3_pf], [gen1_pf, gen2_pf, load_pf, line12_pf, line13_pf, line23_pf])\n", + "\n", + "# Logging\n", + "logger_pf = dpsimpy.Logger(sim_name_pf)\n", + "logger_pf.log_attribute(\"v_bus1\", \"v\", n1_pf)\n", + "logger_pf.log_attribute(\"v_bus2\", \"v\", n2_pf)\n", + "logger_pf.log_attribute(\"v_bus3\", \"v\", n3_pf)\n", + "\n", + "# Simulation\n", + "sim_pf = dpsimpy.Simulation(sim_name_pf, dpsimpy.LogLevel.debug)\n", + "sim_pf.set_system(system_pf)\n", + "sim_pf.set_time_step(time_step_pf)\n", + "sim_pf.set_final_time(final_time_pf)\n", + "sim_pf.set_domain(dpsimpy.Domain.SP)\n", + "sim_pf.set_solver(dpsimpy.Solver.NRP)\n", + "sim_pf.set_solver_component_behaviour(dpsimpy.SolverBehaviour.Initialization)\n", + "sim_pf.do_init_from_nodes_and_terminals(False)\n", + "sim_pf.add_logger(logger_pf)\n", + "sim_pf.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Dynamic Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "sim_name_emt = sim_name + \"_EMT\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_emt)\n", + "\n", + "#Components\n", + "n1_emt = dpsimpy.emt.SimNode(\"n1\", dpsimpy.PhaseType.ABC)\n", + "n2_emt = dpsimpy.emt.SimNode(\"n2\", dpsimpy.PhaseType.ABC)\n", + "n3_emt = dpsimpy.emt.SimNode(\"n3\", dpsimpy.PhaseType.ABC)\n", + "\n", + "#Synchronous generator 1\n", + "gen1_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen1_emt.set_standard_parameters_PU(nom_power_G1, nom_ph_ph_volt_RMS_G1, nom_freq_G1,\n", + " Xpd_G1 * t1_ratio**2, scale_inertia * H_G1, Rs_G1,\n", + " Kd_G1)\n", + "init_apparent_power_G1 = gen1_pf.get_apparent_power()\n", + "gen1_emt.set_initial_values(init_apparent_power_G1, init_mech_power_G1)\n", + "\n", + "#Synchronous generator 2\n", + "gen2_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen2_emt.set_standard_parameters_PU(nom_power_G2, nom_ph_ph_volt_RMS_G2, nom_freq_G2,\n", + " Xpd_G2 * t2_ratio**2, scale_inertia * H_G2, Rs_G2,\n", + " Kd_G2)\n", + "init_apparent_power_G2 = gen2_pf.get_apparent_power()\n", + "gen2_emt.set_initial_values(init_apparent_power_G2, init_mech_power_G2)\n", + "\n", + "gen2_emt.set_model_flags(convert_with_omega_mech=True)\n", + "gen2_emt.set_reference_omega(\"w_r\", gen1_emt, \"delta_r\", gen1_emt)\n", + "\n", + "# Load\n", + "load_emt = dpsimpy.emt.ph3.RXLoad(\"Load\", dpsimpy.LogLevel.debug)\n", + "load_emt.set_parameters(dpsimpy.Math.single_phase_power_to_three_phase(active_power_L),\n", + " dpsimpy.Math.single_phase_power_to_three_phase(reactive_power_L), V_nom)\n", + "\n", + "#Line12\n", + "line12_emt = dpsimpy.emt.ph3.PiLine(\"PiLine12\", dpsimpy.LogLevel.debug)\n", + "line12_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance12),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance12),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance12),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance12))\n", + "\n", + "#Line13\n", + "line13_emt = dpsimpy.emt.ph3.PiLine(\"PiLine13\", dpsimpy.LogLevel.debug)\n", + "line13_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance13),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance13),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance13),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance13))\n", + "\n", + "#Line23\n", + "line23_emt = dpsimpy.emt.ph3.PiLine(\"PiLine23\", dpsimpy.LogLevel.debug)\n", + "line23_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance23),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance23),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance23),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance23))\n", + "\n", + "# Topology\n", + "gen1_emt.connect([n1_emt])\n", + "gen2_emt.connect([n2_emt])\n", + "load_emt.connect([n3_emt])\n", + "line12_emt.connect([n1_emt, n2_emt])\n", + "line13_emt.connect([n1_emt, n3_emt])\n", + "line23_emt.connect([n2_emt, n3_emt])\n", + "system_emt = dpsimpy.SystemTopology(60, [n1_emt, n2_emt, n3_emt], [gen1_emt, gen2_emt, load_emt, line12_emt, line13_emt, line23_emt])\n", + "\n", + "# Initialization of dynamic topology\n", + "system_emt.init_with_powerflow(system_pf)\n", + "\n", + "# Logging\n", + "logger_emt = dpsimpy.Logger(sim_name_emt)\n", + "logger_emt.log_attribute(\"v1\", \"v\", n1_emt)\n", + "logger_emt.log_attribute(\"v2\", \"v\", n2_emt)\n", + "logger_emt.log_attribute(\"v3\", \"v\", n3_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_line12\", \"v_intf\", line12_emt)\n", + "logger_emt.log_attribute(\"v_line13\", \"v_intf\", line13_emt)\n", + "logger_emt.log_attribute(\"v_line23\", \"v_intf\", line23_emt)\n", + "logger_emt.log_attribute(\"i_line12\", \"i_intf\", line12_emt)\n", + "logger_emt.log_attribute(\"i_line13\", \"i_intf\", line13_emt)\n", + "logger_emt.log_attribute(\"i_line23\", \"i_intf\", line23_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep_gen1\", \"Ep_mag\", gen1_emt)\n", + "logger_emt.log_attribute(\"v_gen1\", \"v_intf\", gen1_emt)\n", + "logger_emt.log_attribute(\"i_gen1\", \"i_intf\", gen1_emt)\n", + "logger_emt.log_attribute(\"wr_gen1\", \"w_r\", gen1_emt)\n", + "logger_emt.log_attribute(\"delta_gen1\", \"delta_r\", gen1_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep_gen2\", \"Ep_mag\", gen2_emt)\n", + "logger_emt.log_attribute(\"v_gen2\", \"v_intf\", gen2_emt)\n", + "logger_emt.log_attribute(\"i_gen2\", \"i_intf\", gen2_emt)\n", + "logger_emt.log_attribute(\"wr_gen2\", \"w_r\", gen2_emt)\n", + "logger_emt.log_attribute(\"delta_gen2\", \"delta_r\", gen2_emt)\n", + "\n", + "logger_emt.log_attribute(\"P_elec2\", \"P_elec\", gen2_emt)\n", + "\n", + "\n", + "# Simulation\n", + "sim_emt = dpsimpy.Simulation(sim_name_emt, dpsimpy.LogLevel.debug)\n", + "sim_emt.set_system(system_emt)\n", + "sim_emt.set_time_step(time_step)\n", + "sim_emt.set_final_time(final_time)\n", + "sim_emt.set_domain(dpsimpy.Domain.EMT)\n", + "sim_emt.add_logger(logger_emt)\n", + "\n", + "sim_emt.run()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## DP Simulation" - ], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Parameters" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "#-----------Power system-----------\n", "#Voltage level as Base Voltage\n", - "V_nom = 230e3\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", "\n", "#-----------Generator 1 (bus1)-----------\n", "#Machine parameters\n", @@ -52,13 +354,13 @@ "nom_ph_ph_volt_RMS_G1=25e3\n", "nom_freq_G1 = 60\n", "H_G1 = 6\n", - "Xpd_G1 = 0.3 # in p.u\n", - "Rs_G1 = 0\n", - "Kd_G1 = 1.5\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", "# Initialization parameters \n", "init_active_power_G1 = 270e6\n", "init_mech_power_G1 = 270e6\n", - "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1 * 1.05\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", "\n", "#-----------Generator 2 (bus2)-----------\n", "#Machine parameters\n", @@ -66,9 +368,9 @@ "nom_ph_ph_volt_RMS_G2=13.8e3\n", "nom_freq_G2 = 60\n", "H_G2 = 2\n", - "Xpd_G2 = 0.1 # in p.u\n", - "Rs_G2 = 0\n", - "Kd_G2 = 1.5\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", "# Initialization parameters \n", "init_active_power_G2 = 45e6\n", "init_mech_power_G2 = 45e6\n", @@ -85,19 +387,19 @@ "#-----------Transmission Lines-----------\n", "#PiLine parameters\n", "#line 1-2 (180km)\n", - "line_resistance12 = 0.04*180\n", - "line_inductance12 = (0.4/377)*180\n", - "line_capacitance12 = (4.3e-6/377)*180\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", "line_conductance12 = 0\n", "#line 1-3 (150km)\n", - "line_resistance13 = 0.0267*150\n", - "line_inductance13 = (0.267/377)*150\n", - "line_capacitance13 = (4.3e-6/377)*150\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", "line_conductance13 = 0\n", "#line 2-3 (80km)\n", - "line_resistance23 = 0.04*80\n", - "line_inductance23 = (0.267/377)*80\n", - "line_capacitance23 = (4.3e-6/377)*80\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", "line_conductance23 = 0\n", "\n", "#Switch to trigger fault at generator terminal\n", @@ -109,20 +411,20 @@ "final_time = 10\n", "time_step = 0.001\n", "scale_inertia= 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "time_step_pf = final_time\n", "final_time_pf = final_time + time_step_pf\n", @@ -194,20 +496,20 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic Simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_dp = sim_name + \"_DP\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_dp)\n", @@ -235,6 +537,9 @@ "init_apparent_power_G2 = gen2_pf.get_apparent_power()\n", "gen2_dp.set_initial_values(init_apparent_power_G2, init_mech_power_G2)\n", "\n", + "gen2_dp.set_model_flags(convert_with_omega_mech=True)\n", + "gen2_dp.set_reference_omega(\"w_r\", gen1_dp, \"delta_r\", gen1_dp)\n", + "\n", "# Load\n", "load_dp = dpsimpy.dp.ph1.RXLoad(\"Load\", dpsimpy.LogLevel.debug)\n", "load_dp.set_parameters(active_power_L, reactive_power_L, V_nom)\n", @@ -300,31 +605,33 @@ "sim_dp.add_logger(logger_dp)\n", "\n", "sim_dp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## SP Simulation" - ], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Parameters" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "#-----------Power system-----------\n", "#Voltage level as Base Voltage\n", - "V_nom = 230e3\n", + "V_nom = 250e3\n", + "nom_freq = 60;\n", + "nom_omega= nom_freq*2*np.pi;\n", "\n", "#-----------Generator 1 (bus1)-----------\n", "#Machine parameters\n", @@ -332,13 +639,13 @@ "nom_ph_ph_volt_RMS_G1=25e3\n", "nom_freq_G1 = 60\n", "H_G1 = 6\n", - "Xpd_G1 = 0.3 # in p.u\n", - "Rs_G1 = 0\n", - "Kd_G1 = 1.5\n", + "Xpd_G1 = 0.3 + 0.15 # in p.u\n", + "Rs_G1 = 0.003*0\n", + "Kd_G1 = 1.5 # in p.u\n", "# Initialization parameters \n", "init_active_power_G1 = 270e6\n", "init_mech_power_G1 = 270e6\n", - "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1 * 1.05\n", + "set_point_voltage_G1 = nom_ph_ph_volt_RMS_G1\n", "\n", "#-----------Generator 2 (bus2)-----------\n", "#Machine parameters\n", @@ -346,9 +653,9 @@ "nom_ph_ph_volt_RMS_G2=13.8e3\n", "nom_freq_G2 = 60\n", "H_G2 = 2\n", - "Xpd_G2 = 0.1 # in p.u\n", - "Rs_G2 = 0\n", - "Kd_G2 = 1.5\n", + "Xpd_G2 = 0.3 + 0.12 # in p.u\n", + "Rs_G2 = 0.003*0\n", + "Kd_G2 = 1\n", "# Initialization parameters \n", "init_active_power_G2 = 45e6\n", "init_mech_power_G2 = 45e6\n", @@ -365,19 +672,19 @@ "#-----------Transmission Lines-----------\n", "#PiLine parameters\n", "#line 1-2 (180km)\n", - "line_resistance12 = 0.04*180\n", - "line_inductance12 = (0.4/377)*180\n", - "line_capacitance12 = (4.3e-6/377)*180\n", + "line_resistance12 = 7.2\n", + "line_inductance12 = 72/nom_omega\n", + "line_capacitance12 = 774*1e-6/nom_omega\n", "line_conductance12 = 0\n", "#line 1-3 (150km)\n", - "line_resistance13 = 0.0267*150\n", - "line_inductance13 = (0.267/377)*150\n", - "line_capacitance13 = (4.3e-6/377)*150\n", + "line_resistance13 = 4\n", + "line_inductance13 = 40/nom_omega\n", + "line_capacitance13 = 645*1e-6/nom_omega\n", "line_conductance13 = 0\n", "#line 2-3 (80km)\n", - "line_resistance23 = 0.04*80\n", - "line_inductance23 = (0.267/377)*80\n", - "line_capacitance23 = (4.3e-6/377)*80\n", + "line_resistance23 = 3.2\n", + "line_inductance23 = 32/nom_omega\n", + "line_capacitance23 = 344*1e-6/nom_omega\n", "line_conductance23 = 0\n", "\n", "#Switch to trigger fault at generator terminal\n", @@ -389,20 +696,20 @@ "final_time = 10\n", "time_step = 0.001\n", "scale_inertia= 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "time_step_pf = final_time\n", "final_time_pf = final_time + time_step_pf\n", @@ -473,20 +780,20 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic Simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_sp = sim_name + \"_SP\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_sp)\n", @@ -513,6 +820,8 @@ " Kd_G2)\n", "init_apparent_power_G2 = gen2_pf.get_apparent_power()\n", "gen2_sp.set_initial_values(init_apparent_power_G2, init_mech_power_G2)\n", + "gen2_sp.set_model_flags(convert_with_omega_mech=True)\n", + "gen2_sp.set_reference_omega(\"w_r\", gen1_sp, \"delta_r\", gen1_sp)\n", "\n", "# Load\n", "load_sp = dpsimpy.sp.ph1.Load(\"Load\", dpsimpy.LogLevel.debug)\n", @@ -579,58 +888,79 @@ "sim_sp.add_logger(logger_sp)\n", "\n", "sim_sp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## SP Results" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/SP_SynGenTrStab_3Bus_SteadyState_SP/'\n", "log_name = 'SP_SynGenTrStab_3Bus_SteadyState_SP'\n", "print(work_dir + log_name + '.csv')\n", "ts_sp1ph_TrStab_dl= rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## DP Results" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/DP_SynGenTrStab_3Bus_SteadyState_DP/'\n", "log_name = 'DP_SynGenTrStab_3Bus_SteadyState_DP'\n", "print(work_dir + log_name + '.csv')\n", "ts_dp1ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "## EMT Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "work_dir = 'logs/EMT_SynGenTrStab_3Bus_SteadyState_EMT/'\n", + "log_name = 'EMT_SynGenTrStab_3Bus_SteadyState_EMT'\n", + "print(work_dir + log_name + '.csv')\n", + "ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator 1&2 terminal voltage" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "timestep=50e-6;\n", "t_begin=0\n", @@ -647,6 +977,7 @@ "for name in ['v_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -655,23 +986,22 @@ "for name in ['v_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator 1&2 terminal Current" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", "\n", @@ -681,6 +1011,7 @@ "for name in ['i_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -689,23 +1020,22 @@ "for name in ['i_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator 1&2 Rotor frequency" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", "\n", @@ -719,6 +1049,7 @@ "for name in ['wr_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='DP', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -731,23 +1062,22 @@ "for name in ['wr_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='DP', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", + "plt.legend()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator 1&2 Rotor angular velocity $\\omega _r$" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", "\n", @@ -761,6 +1091,7 @@ "for name in ['wr_gen1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", "plt.legend()\n", "\n", "f.add_subplot(1, 2, 2)\n", @@ -773,66 +1104,45 @@ "for name in ['wr_gen2']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", + "plt.legend()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator 1&2 Rotor angle $\\delta _r$" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "f= plt.figure(figsize=(16,8))\n", "\n", "f.add_subplot(1, 2, 1)\n", "plt.xlabel('time (s)', fontsize=20)\n", - "plt.ylabel('Generator 1 Rotor angle (degree)', fontsize=20)\n", + "plt.ylabel('$\\delta_2$ - $\\delta_1$ (degree)', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "plt.ylim(-20,20)\n", - "\n", - "for name in ['delta_gen1']:\n", - " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", - " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", "\n", - "f.add_subplot(1, 2, 2)\n", - "plt.xlabel('time (s)', fontsize=20)\n", - "plt.ylabel('Generator 2 Rotor angle (degree)', fontsize=20)\n", - "plt.xticks(fontsize=18)\n", - "plt.yticks(fontsize=18)\n", - "plt.ylim(-20,20)\n", + "plt.ylim(-50,50)\n", "\n", - "for name in ['delta_gen2']:\n", - " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", - " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", - "plt.legend()\n", - "\n", - "f.show()" - ], - "outputs": [], - "metadata": {} - }, - { - "cell_type": "code", - "execution_count": null, - "source": [], - "outputs": [], - "metadata": {} + "for name in ['delta_gen']:\n", + " plt.plot(ts_sp1ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_sp1ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' SP backshift')\n", + " plt.plot(ts_dp1ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_dp1ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'1'].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name+'2'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14 - ts_emt3ph_TrStab_dl[name+'1'].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' EMT')\n", + "plt.legend()" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", + "language": "python", "name": "python3" }, "language_info": { @@ -845,9 +1155,8 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.7" + "version": "3.9.13" }, - "orig_nbformat": 3, "tests": { "skip": false } diff --git a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_Fault.ipynb b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_Fault.ipynb similarity index 65% rename from examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_Fault.ipynb rename to examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_Fault.ipynb index 2b6dd3b650..4f6fa5ec91 100644 --- a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_Fault.ipynb +++ b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_Fault.ipynb @@ -2,14 +2,18 @@ "cells": [ { "cell_type": "markdown", + "metadata": {}, "source": [ - "# SMIB Fault SP vs DP" - ], - "metadata": {} + "# SMIB Fault EMT vs SP vs DP" + ] }, { "cell_type": "code", "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], "source": [ "import villas.dataprocessing.readtools as rt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", @@ -21,27 +25,278 @@ "import dpsimpy\n", "\n", "#%matplotlib widget" - ], + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## EMT Simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Parametrization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Power System\n", + "V_nom = 230e3\n", + "\n", + "# Generator\n", + "nom_power = 500e6\n", + "nom_ph_ph_volt_RMS = 22e3\n", + "nom_freq = 60\n", + "nom_omega = nom_freq * 2 * math.pi\n", + "H = 5\n", + "Xpd = 0.31\n", + "Rs = 0.003 * 0\n", + "D=1.5\n", + "\n", + "# Initialization parameters\n", + "init_mech_power = 300e6\n", + "init_active_power = 300e6\n", + "set_point_voltage = nom_ph_ph_volt_RMS * 1.05\n", + "\n", + "# Transformer\n", + "t_ratio = V_nom / nom_ph_ph_volt_RMS\n", + "\n", + "# PiLine parameters calculated from CIGRE Benchmark system\n", + "line_length = 100\n", + "line_resistance = 1.27e-4 * 529 * line_length\n", + "line_inductance = 9.05e-4 * 529 / nom_omega * line_length\n", + "line_capacitance = (1.81e-3 / 529) / nom_omega * line_length\n", + "line_conductance = 0\n", + "\n", + "# Parameters for powerflow initialization\n", + "# Slack voltage: 1pu\n", + "V_slack = V_nom\n", + "\n", + "# Switch to trigger fault at generator terminal\n", + "switch_open = 1e6\n", + "switch_closed = 0.1\n", + "\n", + "# Simulation parameters\n", + "sim_name = \"EMT_SynGenTrStab_SMIB_Fault\"\n", + "final_time = 30\n", + "time_step = 0.001\n", + "start_fault_event = True\n", + "end_fault_event = True\n", + "start_time_fault = 10\n", + "end_time_fault = 10.2\n", + "cmd_Inertia = 1.0\n", + "cmd_Damping = 1.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Powerflow for Initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sim_name_pf = sim_name + \"_PF\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", + "time_step_pf = final_time\n", + "final_time_pf = final_time + time_step_pf\n", + "\n", + "# Components\n", + "n1_pf = dpsimpy.sp.SimNode(\"n1\", dpsimpy.PhaseType.Single)\n", + "n2_pf = dpsimpy.sp.SimNode(\"n2\", dpsimpy.PhaseType.Single)\n", + "\n", + "# Synchronous generator ideal model\n", + "gen_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "gen_pf.set_parameters(rated_apparent_power=nom_power, rated_voltage=nom_ph_ph_volt_RMS,\n", + " set_point_active_power=init_active_power,\n", + " set_point_voltage=set_point_voltage*t_ratio,\n", + " powerflow_bus_type=dpsimpy.PowerflowBusType.PV)\n", + "gen_pf.set_base_voltage(V_nom)\n", + "gen_pf.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.PV)\n", + "\n", + "# Grid bus as Slack\n", + "extnet_pf = dpsimpy.sp.ph1.NetworkInjection(\"Slack\", dpsimpy.LogLevel.debug)\n", + "extnet_pf.set_parameters(V_nom)\n", + "extnet_pf.set_base_voltage(V_nom)\n", + "extnet_pf.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.VD)\n", + "\n", + "# Line\n", + "line_pf = dpsimpy.sp.ph1.PiLine(\"PiLine\", dpsimpy.LogLevel.debug)\n", + "line_pf.set_parameters(line_resistance, line_inductance, line_capacitance, line_conductance)\n", + "line_pf.set_base_voltage(V_nom)\n", + "\n", + "# Switch\n", + "fault_pf = dpsimpy.sp.ph1.Switch(\"Br_fault\", dpsimpy.LogLevel.debug)\n", + "fault_pf.set_parameters(switch_open, switch_closed)\n", + "fault_pf.open()\n", + "\n", + "# Topology\n", + "gen_pf.connect([n1_pf])\n", + "fault_pf.connect([dpsimpy.sp.SimNode.gnd, n1_pf])\n", + "line_pf.connect([n1_pf, n2_pf])\n", + "extnet_pf.connect([n2_pf])\n", + "system_pf = dpsimpy.SystemTopology(nom_freq, [n1_pf, n2_pf], [gen_pf, line_pf, extnet_pf, fault_pf])\n", + "\n", + "# Logging\n", + "logger_pf = dpsimpy.Logger(sim_name_pf)\n", + "logger_pf.log_attribute(\"v1\", \"v\", n1_pf)\n", + "logger_pf.log_attribute(\"v2\", \"v\", n2_pf)\n", + "\n", + "# Simulation\n", + "sim_pf = dpsimpy.Simulation(sim_name_pf, dpsimpy.LogLevel.debug)\n", + "sim_pf.set_system(system_pf)\n", + "sim_pf.set_time_step(time_step_pf)\n", + "sim_pf.set_final_time(final_time_pf)\n", + "sim_pf.set_domain(dpsimpy.Domain.SP)\n", + "sim_pf.set_solver(dpsimpy.Solver.NRP)\n", + "sim_pf.set_solver_component_behaviour(dpsimpy.SolverBehaviour.Initialization)\n", + "sim_pf.do_init_from_nodes_and_terminals(False)\n", + "sim_pf.add_logger(logger_pf)\n", + "sim_pf.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Dynamic simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, "outputs": [], - "metadata": {} + "source": [ + "sim_name_emt = sim_name + \"_EMT\" \n", + "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_emt)\n", + "\n", + "# Nodes\n", + "n1_emt = dpsimpy.emt.SimNode(\"n1\", dpsimpy.PhaseType.ABC)\n", + "n2_emt = dpsimpy.emt.SimNode(\"n2\", dpsimpy.PhaseType.ABC)\n", + "\n", + "# Components\n", + "gen_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen_emt.set_standard_parameters_PU(\n", + " nom_power=nom_power, nom_volt = nom_ph_ph_volt_RMS, nom_freq=nom_freq,\n", + " Xpd = Xpd*t_ratio**2, inertia=cmd_Inertia*H, Rs=Rs, D=cmd_Damping*D)\n", + "\n", + "init_apparent_power = gen_pf.get_apparent_power()\n", + "gen_emt.set_initial_values(init_apparent_power, init_mech_power)\n", + "\n", + "# Grid bus as Slack\n", + "extnet_emt = dpsimpy.emt.ph3.NetworkInjection(\"Slack\", dpsimpy.LogLevel.debug)\n", + "extnet_emt.set_parameters(dpsimpy.Math.single_phase_variable_to_three_phase(complex(V_slack, 0)), 60)\n", + "\n", + "# Line\n", + "line_emt = dpsimpy.emt.ph3.PiLine(\"PiLine\", dpsimpy.LogLevel.debug)\n", + "line_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance))\n", + "\n", + "# Switch\n", + "fault_emt = dpsimpy.emt.ph3.varResSwitch(\"Br_fault\", dpsimpy.LogLevel.debug)\n", + "fault_emt.set_parameters(dpsimpy.Math.single_phase_parameter_to_three_phase(switch_open), dpsimpy.Math.single_phase_parameter_to_three_phase(switch_closed))\n", + "fault_emt.set_init_parameters(time_step)\n", + "fault_emt.open()\n", + "\n", + "# Topology\n", + "gen_emt.connect([n1_emt])\n", + "line_emt.connect([n1_emt, n2_emt])\n", + "extnet_emt.connect([n2_emt])\n", + "fault_emt.connect([dpsimpy.emt.SimNode.gnd, n1_emt])\n", + "system_emt = dpsimpy.SystemTopology(nom_freq, [n1_emt, n2_emt], [gen_emt, line_emt, fault_emt, extnet_emt])\n", + "\n", + "\n", + "# Initialization of dynamic topology\n", + "system_emt.init_with_powerflow(system_pf)\n", + "\n", + "# Logging\n", + "logger_emt = dpsimpy.Logger(sim_name_emt)\n", + "logger_emt.log_attribute(\"v1\", \"v\", n1_emt)\n", + "logger_emt.log_attribute(\"v2\", \"v\", n2_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_line\", \"v_intf\", line_emt)\n", + "logger_emt.log_attribute(\"i_line\", \"i_intf\", line_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_slack\", \"v_intf\", extnet_emt)\n", + "logger_emt.log_attribute(\"i_slack\", \"i_intf\", extnet_emt)\n", + "\n", + "logger_emt.log_attribute(\"i_fault\", \"i_intf\", fault_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep\", \"Ep\", gen_emt)\n", + "logger_emt.log_attribute(\"v_gen\", \"v_intf\", gen_emt)\n", + "logger_emt.log_attribute(\"i_gen\", \"i_intf\", gen_emt)\n", + "logger_emt.log_attribute(\"wr_gen\", \"w_r\", gen_emt)\n", + "logger_emt.log_attribute(\"delta_r_gen\", \"delta_r\", gen_emt)\n", + "logger_emt.log_attribute(\"P_elec\", \"P_elec\", gen_emt)\n", + "logger_emt.log_attribute(\"P_mech\", \"P_mech\", gen_emt)\n", + "\n", + "# Simulation\n", + "sim_emt = dpsimpy.Simulation(sim_name, dpsimpy.LogLevel.debug)\n", + "sim_emt.set_system(system_emt)\n", + "sim_emt.set_time_step(time_step)\n", + "sim_emt.set_final_time(final_time)\n", + "sim_emt.set_domain(dpsimpy.Domain.EMT)\n", + "sim_emt.add_logger(logger_emt)\n", + "sim_emt.do_system_matrix_recomputation(True)\n", + "sim_emt.set_direct_solver_implementation(dpsimpy.DirectLinearSolverImpl.SparseLU)\n", + "\n", + "\n", + "# Events\n", + "if start_fault_event:\n", + " sw1 = dpsimpy.event.SwitchEvent3Ph(start_time_fault, fault_emt, True)\n", + " sim_emt.add_event(sw1)\n", + " \n", + "if end_fault_event:\n", + " sw2 = dpsimpy.event.SwitchEvent3Ph(end_time_fault, fault_emt, False)\n", + " sim_emt.add_event(sw2)\n", + " \n", + "sim_emt.run()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## DP Simulation" - ], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Parametrization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], "source": [ "# Power System\n", "V_nom = 230e3\n", @@ -85,20 +340,20 @@ "end_time_fault = 10.2\n", "cmd_Inertia = 1.0\n", "cmd_Damping = 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_pf = sim_name + \"_PF\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", @@ -157,20 +412,20 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_dp = sim_name + \"_DP\" \n", "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_dp)\n", @@ -256,27 +511,27 @@ " sim_dp.add_event(sw2)\n", " \n", "sim_dp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## SP Simulation" - ], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Parametrization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# Power System\n", "V_nom = 230e3\n", @@ -320,20 +575,20 @@ "end_time_fault = 10.2\n", "cmd_Inertia = 1.0\n", "cmd_Damping = 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_pf = sim_name + \"_PF\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", @@ -392,20 +647,20 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_sp = sim_name + \"_SP\" \n", "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_sp)\n", @@ -490,78 +745,77 @@ " sim_sp.add_event(sw2)\n", " \n", "sim_sp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Results 1ph SP" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/SP_SynGenTrStab_SMIB_Fault_SP/'\n", "log_name = 'SP_SynGenTrStab_SMIB_Fault_SP'\n", "print(work_dir + log_name + '.csv')\n", "ts_sp1ph_TrStab_dl= rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Results 1ph DP" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/DP_SynGenTrStab_SMIB_Fault_DP/' \n", "log_name = 'DP_SynGenTrStab_SMIB_Fault_DP'\n", "print(work_dir + log_name + '.csv')\n", "ts_dp1ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Results 3ph EMT Reference" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "# work_dir = 'logs/EMT_SynGenTrStab_SMIB_Fault_EMT/'\n", - "# log_name = 'EMT_SynGenTrStab_SMIB_Fault_EMT'\n", - "# print(work_dir + log_name + '.csv')\n", - "# #os.path.getsize(work_dir + log_name + '.csv')\n", - "# ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "work_dir = 'logs/EMT_SynGenTrStab_SMIB_Fault_EMT/'\n", + "log_name = 'EMT_SynGenTrStab_SMIB_Fault_EMT'\n", + "print(work_dir + log_name + '.csv')\n", + "ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator emf" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator emf (V)')\n", @@ -577,233 +831,208 @@ "for name in ['Ep']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' EMT')" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' EMT', linestyle='-.')\n", + "plt.legend(fontsize=18)\n", + "plt.show()" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator terminal voltage" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator terminal voltage (V)')\n", "\n", - "for name in ['v_gen']:\n", + "for name in ['v1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": { + "tags": [] + }, "source": [ "## Genrerator terminal Current" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator terminal current (A)')\n", "\n", "for name in ['i_gen']:\n", - " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], -ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", - " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], -ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "#plt.plot(ts_dp1ph_TrStab_dl['Ep'].interpolate(timestep).time[begin_idx:end_idx], 0.05*(ts_dp1ph_TrStab_dl['Ep'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]- ts_dp1ph_TrStab_dl['v_gen'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]))\n", - "#plt.plot(ts_sp1ph_TrStab_dl['Ep'].interpolate(timestep).time[begin_idx:end_idx], 0.05*(ts_sp1ph_TrStab_dl['Ep'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]- ts_sp1ph_TrStab_dl['v_gen'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]))\n", - "plt.legend()\n", - "plt.show()" - ], - "outputs": [], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "Fault Current" - ], - "metadata": {} - }, - { - "cell_type": "code", - "execution_count": null, - "source": [ - "plt.figure(figsize=(12,8))\n", - "plt.ylabel('Fault current (A)')\n", - "\n", - "for name in ['i_fault']:\n", - " #plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", + " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "#plt.plot(ts_dp1ph_TrStab_dl['Ep'].interpolate(timestep).time[begin_idx:end_idx], 0.05*(ts_dp1ph_TrStab_dl['Ep'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]- ts_dp1ph_TrStab_dl['v_gen'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]))\n", - "#plt.plot(ts_sp1ph_TrStab_dl['Ep'].interpolate(timestep).time[begin_idx:end_idx], 0.05*(ts_sp1ph_TrStab_dl['Ep'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]- ts_sp1ph_TrStab_dl['v_gen'].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx]))\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Voltage across line" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "\n", "for name in ['v_line']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Current through line" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "\n", "for name in ['i_line']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator electrical & mechanical energy" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", - "#plt.ylim(250e6,350e6)\n", + "# plt.ylim(250e6,350e6)\n", "\n", "for name in ['P_elec','P_mech']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor frequency" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "#plt.figure(figsize=(12,8))\n", - "#plt.xlabel('time (s)', fontsize=20)\n", - "#plt.ylabel('Rotor frequency (Hz)', fontsize=20)\n", - "#plt.xticks(fontsize=18)\n", - "#plt.yticks(fontsize=18)\n", - "##plt.ylim(55,65)\n", - "##plt.xlim(5.1,10)\n", - "\n", - "##if x_axis limits are changed above, change them again to consider the complete duration\n", - "\n", - "#for name in ['wr_gen']:\n", - "# plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$\\f_r$' + ' SP backshift')\n", - "# plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$\\f_r$' + ' DP backshift', linestyle='--')\n", - "# #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "#plt.legend(fontsize=18)\n", - "##plt.title('Post-Fault',fontsize=20)\n", - "#plt.show()\n" - ], + "metadata": { + "tags": [] + }, "outputs": [], - "metadata": {} + "source": [ + "plt.figure(figsize=(12,8))\n", + "plt.xlabel('time (s)', fontsize=20)\n", + "plt.ylabel('Rotor frequency (Hz)', fontsize=20)\n", + "plt.xticks(fontsize=18)\n", + "plt.yticks(fontsize=18)\n", + "plt.ylim(55,65)\n", + "\n", + "for name in ['wr_gen']:\n", + " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' SP backshift')\n", + " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", + "plt.legend(fontsize=18)\n", + "plt.show()\n" + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor angular velocity $\\omega _r$" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.xlabel('time (s)', fontsize=20)\n", "plt.ylabel('Rotor angular velocity', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(360,400)\n", - "#fix y axis\n", + "# plt.ylim(360,400)\n", "\n", "for name in ['wr_gen']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", "plt.legend(fontsize=18)\n", - "#plt.title('Post-Fault',fontsize=20)\n", "plt.show()\n" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor angle $\\delta _r$" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.xlabel('time (s)', fontsize=20)\n", - "#plt.ylim(-10,10)\n", + "# plt.ylim(0,50)\n", "plt.ylabel('Rotor angle in degrees', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", @@ -811,26 +1040,16 @@ "for name in ['delta_r_gen']: \n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label='$\\delta _r$' + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label='$\\delta _r$' + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' EMT')\n", " \n", "plt.legend(fontsize=18)\n", - "#plt.title('Post-Fault',fontsize=20)\n", "plt.show()\n" - ], - "outputs": [], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "Frequency of load angle in 1hz-2hz range" - ], - "metadata": {} + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -844,7 +1063,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.9.13" }, "tests": { "skip": false diff --git a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb similarity index 64% rename from examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb rename to examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb index 87a001197c..029fb69078 100644 --- a/examples/Notebooks/Circuits/DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb +++ b/examples/Notebooks/Circuits/EMT_DP_SP_SynGenTrStab_SMIB_SteadyState.ipynb @@ -1,15 +1,18 @@ { "cells": [ { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ - "# SMIB Steady state SP vs DP" - ], - "metadata": {} + "# SMIB Steady state EMT vs DP vs SP" + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "import villas.dataprocessing.readtools as rt\n", "from villas.dataprocessing.timeseries import TimeSeries as ts\n", @@ -21,27 +24,244 @@ "import dpsimpy\n", "\n", "#%matplotlib widget" - ], + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## EMT Simulation" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Parametrization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "# Power System\n", + "V_nom = 230e3\n", + "\n", + "# Generator\n", + "nom_power = 500e6\n", + "nom_ph_ph_volt_RMS = 22e3\n", + "nom_freq = 60\n", + "nom_omega = nom_freq * 2 * math.pi\n", + "H = 5\n", + "Xpd = 0.31\n", + "Rs = 0.003 * 0\n", + "D=1.5\n", + "\n", + "# Initialization parameters\n", + "init_mech_power = 300e6\n", + "init_active_power = 300e6\n", + "set_point_voltage = nom_ph_ph_volt_RMS * 1.05\n", + "\n", + "# Transformer\n", + "t_ratio = V_nom / nom_ph_ph_volt_RMS\n", + "\n", + "# PiLine parameters calculated from CIGRE Benchmark system\n", + "line_resistance = 1.27e-4 * 529 * 100\n", + "line_inductance = 9.05e-4 * 529 * 100/nom_omega\n", + "line_capacitance = (1.81e-3 / 529) * 100 /nom_omega\n", + "line_conductance = 0\n", + "\n", + "# Parameters for powerflow initialization\n", + "# Slack voltage: 1pu\n", + "V_slack = V_nom\n", + "\n", + "# Switch to trigger fault at generator terminal\n", + "switch_open = 1e6\n", + "switch_closed = 0.1\n", + "\n", + "# Simulation parameters\n", + "sim_name = \"EMT_SynGenTrStab_SMIB_SteadyState\"\n", + "final_time = 10\n", + "time_step = 0.001\n", + "cmd_Inertia = 1.0\n", + "cmd_Damping = 1.0" + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Powerflow for Initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim_name_pf = sim_name + \"_PF\"\n", + "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", + "time_step_pf = final_time\n", + "final_time_pf = final_time + time_step_pf\n", + "\n", + "# Components\n", + "n1_pf = dpsimpy.sp.SimNode(\"n1\", dpsimpy.PhaseType.Single)\n", + "n2_pf = dpsimpy.sp.SimNode(\"n2\", dpsimpy.PhaseType.Single)\n", + "\n", + "# Synchronous generator ideal model\n", + "gen_pf = dpsimpy.sp.ph1.SynchronGenerator(\"Generator\", dpsimpy.LogLevel.debug)\n", + "gen_pf.set_parameters(rated_apparent_power=nom_power, rated_voltage=nom_ph_ph_volt_RMS,\n", + " set_point_active_power=init_active_power,\n", + " set_point_voltage=set_point_voltage*t_ratio,\n", + " powerflow_bus_type=dpsimpy.PowerflowBusType.PV)\n", + "gen_pf.set_base_voltage(V_nom)\n", + "gen_pf.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.PV)\n", + "\n", + "# Grid bus as Slack\n", + "extnet_pf = dpsimpy.sp.ph1.NetworkInjection(\"Slack\", dpsimpy.LogLevel.debug)\n", + "extnet_pf.set_parameters(V_slack)\n", + "extnet_pf.set_base_voltage(V_nom)\n", + "extnet_pf.modify_power_flow_bus_type(dpsimpy.PowerflowBusType.VD)\n", + "\n", + "# Line\n", + "line_pf = dpsimpy.sp.ph1.PiLine(\"PiLine\", dpsimpy.LogLevel.debug)\n", + "line_pf.set_parameters(line_resistance, line_inductance, line_capacitance, line_conductance)\n", + "line_pf.set_base_voltage(V_nom)\n", + "\n", + "# Topology\n", + "gen_pf.connect([n1_pf])\n", + "line_pf.connect([n1_pf, n2_pf])\n", + "extnet_pf.connect([n2_pf])\n", + "system_pf = dpsimpy.SystemTopology(nom_freq, [n1_pf, n2_pf], [gen_pf, line_pf, extnet_pf])\n", + "\n", + "# Logging\n", + "logger_pf = dpsimpy.Logger(sim_name_pf)\n", + "logger_pf.log_attribute(\"v1\", \"v\", n1_pf)\n", + "logger_pf.log_attribute(\"v2\", \"v\", n2_pf)\n", + "\n", + "# Simulation\n", + "sim_pf = dpsimpy.Simulation(sim_name_pf, dpsimpy.LogLevel.debug)\n", + "sim_pf.set_system(system_pf)\n", + "sim_pf.set_time_step(time_step_pf)\n", + "sim_pf.set_final_time(final_time_pf)\n", + "sim_pf.set_domain(dpsimpy.Domain.SP)\n", + "sim_pf.set_solver(dpsimpy.Solver.NRP)\n", + "sim_pf.set_solver_component_behaviour(dpsimpy.SolverBehaviour.Initialization)\n", + "sim_pf.do_init_from_nodes_and_terminals(False)\n", + "sim_pf.add_logger(logger_pf)\n", + "sim_pf.run()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Dynamic simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim_name_emt = sim_name + \"_EMT\" \n", + "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_emt)\n", + "\n", + "# Nodes\n", + "n1_emt = dpsimpy.emt.SimNode(\"n1\", dpsimpy.PhaseType.ABC)\n", + "n2_emt = dpsimpy.emt.SimNode(\"n2\", dpsimpy.PhaseType.ABC)\n", + "\n", + "# Components\n", + "gen_emt = dpsimpy.emt.ph3.SynchronGeneratorTrStab(\"SynGen\", dpsimpy.LogLevel.debug)\n", + "# Xpd is given in p.u of generator base at transfomer primary side and should be transformed to network side\n", + "gen_emt.set_standard_parameters_PU(\n", + " nom_power=nom_power, nom_volt = nom_ph_ph_volt_RMS, nom_freq=nom_freq,\n", + " Xpd = Xpd*t_ratio**2, inertia=cmd_Inertia*H, Rs=Rs, D=cmd_Damping*D)\n", + "\n", + "init_apparent_power = gen_pf.get_apparent_power()\n", + "gen_emt.set_initial_values(init_apparent_power, init_mech_power)\n", + "\n", + "# Grid bus as Slack\n", + "extnet_emt = dpsimpy.emt.ph3.NetworkInjection(\"Slack\", dpsimpy.LogLevel.debug)\n", + "extnet_emt.set_parameters(dpsimpy.Math.single_phase_variable_to_three_phase(complex(V_slack, 0)), 60)\n", + "\n", + "# Line\n", + "line_emt = dpsimpy.emt.ph3.PiLine(\"PiLine\", dpsimpy.LogLevel.debug)\n", + "line_emt.set_parameters(\n", + " series_resistance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_resistance),\n", + " series_inductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_inductance),\n", + " parallel_capacitance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_capacitance),\n", + " parallel_conductance=dpsimpy.Math.single_phase_parameter_to_three_phase(line_conductance))\n", + "\n", + "# Topology\n", + "gen_emt.connect([n1_emt])\n", + "line_emt.connect([n1_emt, n2_emt])\n", + "extnet_emt.connect([n2_emt])\n", + "system_emt = dpsimpy.SystemTopology(nom_freq, [n1_emt, n2_emt], [gen_emt, line_emt, extnet_emt])\n", + "\n", + "\n", + "# Initialization of dynamic topology\n", + "system_emt.init_with_powerflow(system_pf)\n", + "\n", + "# Logging\n", + "logger_emt = dpsimpy.Logger(sim_name_emt)\n", + "logger_emt.log_attribute(\"v1\", \"v\", n1_emt)\n", + "logger_emt.log_attribute(\"v2\", \"v\", n2_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_line\", \"v_intf\", line_emt)\n", + "logger_emt.log_attribute(\"i_line\", \"i_intf\", line_emt)\n", + "\n", + "logger_emt.log_attribute(\"v_slack\", \"v_intf\", extnet_emt)\n", + "logger_emt.log_attribute(\"i_slack\", \"i_intf\", extnet_emt)\n", + "\n", + "logger_emt.log_attribute(\"Ep\", \"Ep\", gen_emt)\n", + "logger_emt.log_attribute(\"v_gen\", \"v_intf\", gen_emt)\n", + "logger_emt.log_attribute(\"i_gen\", \"i_intf\", gen_emt)\n", + "logger_emt.log_attribute(\"wr_gen\", \"w_r\", gen_emt)\n", + "logger_emt.log_attribute(\"delta_r_gen\", \"delta_r\", gen_emt)\n", + "logger_emt.log_attribute(\"P_elec\", \"P_elec\", gen_emt)\n", + "logger_emt.log_attribute(\"P_mech\", \"P_mech\", gen_emt)\n", + "\n", + "# Simulation\n", + "sim_emt = dpsimpy.Simulation(sim_name, dpsimpy.LogLevel.debug)\n", + "sim_emt.set_system(system_emt)\n", + "sim_emt.set_time_step(time_step)\n", + "sim_emt.set_final_time(final_time)\n", + "sim_emt.set_domain(dpsimpy.Domain.EMT)\n", + "sim_emt.add_logger(logger_emt)\n", + " \n", + "sim_emt.run()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, "source": [ "## DP Simulation" - ], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Parametrization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# Power System\n", "V_nom = 230e3\n", @@ -84,20 +304,21 @@ "time_step = 0.001\n", "cmd_Inertia = 1.0\n", "cmd_Damping = 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_pf = sim_name + \"_PF\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", @@ -150,20 +371,21 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_dp = sim_name + \"_DP\" \n", "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_dp)\n", @@ -228,27 +450,29 @@ "sim_dp.add_logger(logger_dp)\n", " \n", "sim_dp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## SP Simulation" - ], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Parametrization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# Power System\n", "V_nom = 230e3\n", @@ -293,20 +517,21 @@ "end_time_fault = 10.2\n", "cmd_Inertia = 1.0\n", "cmd_Damping = 1.0" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Powerflow for Initialization" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_pf = sim_name + \"_PF\"\n", "dpsimpy.Logger.set_log_dir(\"logs/\" + sim_name_pf)\n", @@ -358,20 +583,21 @@ "sim_pf.do_init_from_nodes_and_terminals(False)\n", "sim_pf.add_logger(logger_pf)\n", "sim_pf.run()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "### Dynamic simulation" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "sim_name_sp = sim_name + \"_SP\" \n", "dpsimpy.Logger.set_log_dir(\"logs/\"+sim_name_sp)\n", @@ -443,78 +669,84 @@ "sim_sp.add_logger(logger_sp)\n", " \n", "sim_sp.run()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Results 1ph SP" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/SP_SynGenTrStab_SMIB_SteadyState_SP/'\n", "log_name = 'SP_SynGenTrStab_SMIB_SteadyState_SP'\n", "print(work_dir + log_name + '.csv')\n", "ts_sp1ph_TrStab_dl= rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Results 1ph DP" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "work_dir = 'logs/DP_SynGenTrStab_SMIB_SteadyState_DP/' \n", "log_name = 'DP_SynGenTrStab_SMIB_SteadyState_DP'\n", "print(work_dir + log_name + '.csv')\n", "ts_dp1ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": { + "tags": [] + }, "source": [ - "## Results 3ph EMT Reference" - ], - "metadata": {} + "## Results 3ph EMT" + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "# work_dir = 'logs/EMT_SynGenTrStab_SMIB_SteadyState_EMT/'\n", - "# log_name = 'EMT_SynGenTrStab_SMIB_SteadyState_EMT'\n", - "# print(work_dir + log_name + '.csv')\n", - "# #os.path.getsize(work_dir + log_name + '.csv')\n", - "# ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" - ], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [ + "work_dir = 'logs/EMT_SynGenTrStab_SMIB_SteadyState_EMT/'\n", + "log_name = 'EMT_SynGenTrStab_SMIB_SteadyState_EMT'\n", + "print(work_dir + log_name + '.csv')\n", + "#os.path.getsize(work_dir + log_name + '.csv')\n", + "ts_emt3ph_TrStab_dl = rt.read_timeseries_dpsim(work_dir + log_name + '.csv')" + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator emf" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator emf (V)')\n", @@ -530,45 +762,51 @@ "for name in ['Ep']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' EMT')" - ], - "outputs": [], - "metadata": {} + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' EMT', linestyle='-.')\n", + "plt.legend(fontsize=18)\n", + "plt.show()" + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator terminal voltage" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator terminal voltage (V)')\n", "\n", - "for name in ['v_gen']:\n", + "for name in ['v1']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": { + "tags": [] + }, "source": [ "## Genrerator terminal Current" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylabel('Generator terminal current (A)')\n", @@ -576,70 +814,73 @@ "for name in ['i_gen']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", "\n", - "plt.legend()\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Voltage across line" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "\n", "for name in ['v_line']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Current through line" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "\n", "for name in ['i_line']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + 'SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).frequency_shift(60).values[begin_idx:end_idx], label=name + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name+'_0'].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Generator electrical & mechanical energy" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.ylim(250e6,350e6)\n", @@ -647,91 +888,88 @@ "for name in ['P_elec','P_mech']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' SP')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' DP', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], np.sqrt(3/2)*ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "plt.legend()\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + "plt.legend(fontsize=18)\n", "plt.show()" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor frequency" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [ - "# plt.figure(figsize=(12,8))\n", - "# plt.xlabel('time (s)', fontsize=20)\n", - "# plt.ylabel('Rotor frequency (Hz)', fontsize=20)\n", - "# plt.xticks(fontsize=18)\n", - "# plt.yticks(fontsize=18)\n", - "# #plt.ylim(55,65)\n", - "# #plt.xlim(5.1,10)\n", - "\n", - "# #if x_axis limits are changed above, change them again to consider the complete duration\n", - "\n", - "# for name in ['wr_gen']:\n", - "# plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$\\f_r$' + ' SP backshift')\n", - "# plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$\\f_r$' + ' DP backshift', linestyle='--')\n", - "# #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", - "# plt.legend(fontsize=18)\n", - "# #plt.title('Post-Fault',fontsize=20)\n", - "# plt.show()\n" - ], + "metadata": { + "tags": [] + }, "outputs": [], - "metadata": {} + "source": [ + "plt.figure(figsize=(12,8))\n", + "plt.xlabel('time (s)', fontsize=20)\n", + "plt.ylabel('Rotor frequency (Hz)', fontsize=20)\n", + "plt.xticks(fontsize=18)\n", + "plt.yticks(fontsize=18)\n", + "plt.ylim(55,65)\n", + "\n", + "for name in ['wr_gen']:\n", + " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' SP backshift')\n", + " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' DP backshift', linestyle='--')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*60/377, label='$f_r$' + ' EMT')\n", + "plt.legend(fontsize=18)\n", + "plt.show()\n" + ] }, { - "cell_type": "code", - "execution_count": null, + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor angular velocity $\\omega _r$" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.xlabel('time (s)', fontsize=20)\n", "plt.ylabel('Rotor angular velocity', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", - "#plt.ylim(360,400)\n", - "#fix y axis\n", + "plt.ylim(360,400)\n", "\n", "for name in ['wr_gen']:\n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx],ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label='$\\omega _r$' + ' EMT')\n", "plt.legend(fontsize=18)\n", - "#plt.title('Post-Fault',fontsize=20)\n", "plt.show()\n" - ], - "outputs": [], - "metadata": {} + ] }, { + "attachments": {}, "cell_type": "markdown", + "metadata": {}, "source": [ "## Rotor angle $\\delta _r$" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure(figsize=(12,8))\n", "plt.xlabel('time (s)', fontsize=20)\n", - "#plt.ylim(-10,10)\n", + "plt.ylim(0,50)\n", "plt.ylabel('Rotor angle in degrees', fontsize=20)\n", "plt.xticks(fontsize=18)\n", "plt.yticks(fontsize=18)\n", @@ -739,26 +977,16 @@ "for name in ['delta_r_gen']: \n", " plt.plot(ts_sp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_sp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label='$\\delta _r$' + ' SP backshift')\n", " plt.plot(ts_dp1ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_dp1ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label='$\\delta _r$' + ' DP backshift', linestyle='--')\n", - " #plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx], label=name + ' EMT')\n", + " plt.plot(ts_emt3ph_TrStab_dl[name].interpolate(timestep).time[begin_idx:end_idx], ts_emt3ph_TrStab_dl[name].interpolate(timestep).values[begin_idx:end_idx]*180/3.14, label=name + ' EMT')\n", " \n", "plt.legend(fontsize=18)\n", - "#plt.title('Post-Fault',fontsize=20)\n", "plt.show()\n" - ], - "outputs": [], - "metadata": {} - }, - { - "cell_type": "markdown", - "source": [ - "Frequency of load angle in 1hz-2hz range" - ], - "metadata": {} + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -772,11 +1000,15 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.7" + "version": "3.9.13" }, - "orig_nbformat": 3, "tests": { "skip": false + }, + "vscode": { + "interpreter": { + "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90" + } } }, "nbformat": 4,