diff --git a/Examples/Cxx/CIM/Slack_Trafo3W_Gen_Load.cpp b/Examples/Cxx/CIM/Slack_Trafo3W_Gen_Load.cpp new file mode 100644 index 0000000000..682fdf3a17 --- /dev/null +++ b/Examples/Cxx/CIM/Slack_Trafo3W_Gen_Load.cpp @@ -0,0 +1,67 @@ +/* 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 +#include + +using namespace std; +using namespace DPsim; +using namespace CPS; +using namespace CPS::CIM; + +/* + * This example runs the powerflow for the CIGRE MV benchmark system (neglecting the tap changers of the transformers) + */ +int main(int argc, char** argv) { + + // Find CIM files + std::list filenames; + if (argc <= 1) { + filenames = DPsim::Utils::findFiles({ + "Slack_Trafo3W_PowerFlow_EQ.xml", + "Slack_Trafo3W_PowerFlow_TP.xml", + "Slack_Trafo3W_PowerFlow_SV.xml", + "Slack_Trafo3W_PowerFlow_SSH.xml" + }, "build/_deps/cim-data-src/BasicGrids/PowerFactory/Slack_Trafo3W_PowerFlow", "CIMPATH"); + } + else { + filenames = std::list(argv + 1, argv + argc); + } + + String simName = "Trafo3W_PFSolver"; + CPS::Real system_freq = 60; + + CIM::Reader reader(simName, Logger::Level::debug, Logger::Level::off); + SystemTopology system = reader.loadCIM( + system_freq, + filenames, + CPS::Domain::SP, + CPS::PhaseType::Single, + CPS::GeneratorType::PVNode + ); + system.renderToFile("build/_deps/cim-data-src/BasicGrids/PowerFactory/Slack_Trafo3W_PowerFlow/Trafo3W_PowerFlowTest.svg"); + + auto logger = DPsim::DataLogger::make(simName); + for (auto node : system.mNodes) + { + logger->addAttribute(node->name() + ".V", node->attribute("v")); + } + + Simulation sim(simName, Logger::Level::debug); + sim.setSystem(system); + sim.setTimeStep(1); + sim.setFinalTime(1); + sim.setDomain(Domain::SP); + sim.setSolverType(Solver::Type::NRP); + sim.doInitFromNodesAndTerminals(true); + sim.addLogger(logger); + + sim.run(); + + return 0; +} diff --git a/Examples/Cxx/CMakeLists.txt b/Examples/Cxx/CMakeLists.txt index bbcc668c67..9adb8ab79e 100644 --- a/Examples/Cxx/CMakeLists.txt +++ b/Examples/Cxx/CMakeLists.txt @@ -126,6 +126,7 @@ if(WITH_CIM) # PF(Power Flow) example CIM/Slack_Trafo_Load.cpp + CIM/Slack_Trafo3W_Gen_Load.cpp CIM/Slack_TrafoTapChanger_Load.cpp CIM/CIGRE_MV_PowerFlowTest.cpp CIM/CIGRE_MV_PowerFlowTest_LoadProfiles.cpp diff --git a/Include/dpsim/PFSolver.h b/Include/dpsim/PFSolver.h index fc813ac167..0dd46d19ac 100644 --- a/Include/dpsim/PFSolver.h +++ b/Include/dpsim/PFSolver.h @@ -57,6 +57,8 @@ namespace DPsim { CPS::SystemTopology mSystem; /// Vector of transformer components std::vector> mTransformers; + /// Vector of transformer components + std::vector> mTransformers3W; /// Vector of solid state transformer components std::vector> mSolidStateTransformers; /// Vector of synchronous generator components diff --git a/Source/PFSolver.cpp b/Source/PFSolver.cpp index 330792d66b..84bdbe436d 100644 --- a/Source/PFSolver.cpp +++ b/Source/PFSolver.cpp @@ -28,6 +28,8 @@ void PFSolver::initialize(){ mLoads.push_back(load); else if (std::shared_ptr trafo = std::dynamic_pointer_cast(comp)) mTransformers.push_back(trafo); + else if (std::shared_ptr trafo3W = std::dynamic_pointer_cast(comp)) + mTransformers3W.push_back(trafo3W); else if (std::shared_ptr line = std::dynamic_pointer_cast(comp)) mLines.push_back(line); else if (std::shared_ptr extnet = std::dynamic_pointer_cast(comp)) @@ -86,6 +88,9 @@ void PFSolver::initializeComponents(){ for(auto trans : mTransformers) { trans->calculatePerUnitParameters(mBaseApparentPower, mSystem.mSystemOmega); } + for(auto trans3w : mTransformers3W) { + trans3w->calculatePerUnitParameters(mBaseApparentPower, mSystem.mSystemOmega); + } for(auto shunt : mShunts) { shunt->calculatePerUnitParameters(mBaseApparentPower, mSystem.mSystemOmega); } @@ -113,6 +118,16 @@ void PFSolver::setBaseApparentPower() { if (trafo->attribute("S")->get() > maxPower) maxPower = trafo->attribute("S")->get(); } + else if (!mTransformers3W.empty()) { + for (auto trafo3w : mTransformers3W){ + if (trafo3w->attribute("S1")->get() > maxPower) + maxPower = trafo3w->attribute("S1")->get(); + if (trafo3w->attribute("S2")->get() > maxPower) + maxPower = trafo3w->attribute("S2")->get(); + if (trafo3w->attribute("S3")->get() > maxPower) + maxPower = trafo3w->attribute("S3")->get(); + } + } if (maxPower != 0.) mBaseApparentPower = pow(10, 1 + floor(log10(maxPower))); else @@ -196,6 +211,11 @@ void PFSolver::determinePFBusType() { mSLog->debug("{}: VD, PV and PQ type component connect -> set as VD bus", node->name()); mVDBusIndices.push_back(node->matrixNodeIndex()); mVDBuses.push_back(node); + } // VD and PQ type component connect -> set as VD bus + else if (!connectedPV && connectedPQ && connectedVD) { + mSLog->debug("{}: VD and PQ type component connect -> set as VD bus", node->name()); + mVDBusIndices.push_back(node->matrixNodeIndex()); + mVDBuses.push_back(node); } else { std::stringstream ss; @@ -263,13 +283,25 @@ void PFSolver::composeAdmittanceMatrix() { } trans->pfApplyAdmittanceMatrixStamp(mY); } + for(auto trans3W : mTransformers3W) { + //to check if this transformer could be ignored + if (trans3W->attribute("R1") == 0 && trans3W->attribute("L1") == 0 && + trans3W->attribute("R2") == 0 && trans3W->attribute("L2") == 0 && + trans3W->attribute("R3") == 0 && trans3W->attribute("L3") == 0) { + mSLog->info("{} {} ignored for R = 0 and L = 0 in all windings", trans3W->type(), trans3W->name()); + continue; + } + trans3W->pfApplyAdmittanceMatrixStamp(mY); + } for(auto shunt : mShunts) { shunt->pfApplyAdmittanceMatrixStamp(mY); } } - if(mLines.empty() && mTransformers.empty()) { + if(mLines.empty() && mTransformers.empty() && mTransformers3W.empty()) { throw std::invalid_argument("There are no bus"); } + + mSLog->info("#### System matrix:\n{}", mY); } CPS::Real PFSolver::G(int i, int j) { diff --git a/Source/PFSolverPowerPolar.cpp b/Source/PFSolverPowerPolar.cpp index ce59a683a1..4888c1e032 100644 --- a/Source/PFSolverPowerPolar.cpp +++ b/Source/PFSolverPowerPolar.cpp @@ -305,11 +305,38 @@ void PFSolverPowerPolar::setSolution() { baseVoltage_ = trans->attribute("nominal_voltage_end2")->get(); break; } - else if (std::shared_ptr gen = std::dynamic_pointer_cast(comp)) { - baseVoltage_ =gen->attribute("base_Voltage")->get(); + } + else if (std::shared_ptr trans3W = std::dynamic_pointer_cast(comp)) { + if (trans3W->terminal(0)->node()->name() == node->name()){ + baseVoltage_ = trans3W->attribute("nominal_voltage_end1")->get(); + break; + } + else if (trans3W->terminal(1)->node()->name() == node->name()){ + baseVoltage_ = trans3W->attribute("nominal_voltage_end2")->get(); + break; + } + else if (trans3W->terminal(2)->node()->name() == node->name()){ + baseVoltage_ = trans3W->attribute("nominal_voltage_end3")->get(); break; } - else + } + else if (std::shared_ptr gen = std::dynamic_pointer_cast(comp)) { + baseVoltage_ = gen->attribute("base_Voltage")->get(); + break; + } + else if (std::shared_ptr extnet = std::dynamic_pointer_cast(comp)) { + baseVoltage_ = extnet->attribute("base_Voltage")->get(); + break; + } + else if (std::shared_ptr shunt = std::dynamic_pointer_cast(comp)) { + baseVoltage_ = shunt->attribute("base_Voltage")->get(); + break; + } + else if (std::shared_ptr load = std::dynamic_pointer_cast(comp)) { + baseVoltage_ = load->attribute("V_nom")->get(); + break; + } + else { mSLog->warn("Unable to get base voltage at {}", node->name()); } } @@ -340,6 +367,17 @@ void PFSolverPowerPolar::calculateBranchFlow() { VectorComp flow_on_branch = v.array()*current.conjugate().array(); trafo->updateBranchFlow(current, flow_on_branch); } + for (auto trafo3W : mTransformers3W) { + VectorComp v(3); + v(0) = sol_V_complex.coeff(trafo3W->node(0)->matrixNodeIndex()); + v(1) = sol_V_complex.coeff(trafo3W->node(1)->matrixNodeIndex()); + v(2) = sol_V_complex.coeff(trafo3W->node(2)->matrixNodeIndex()); + /// I = Y * V + VectorComp current = trafo3W->Y_element() * v; + /// pf on branch [S_01; S_10] = [V_0 * conj(I_0); V_1 * conj(I_1)] + VectorComp flow_on_branch = v.array()*current.conjugate().array(); + trafo3W->updateBranchFlow(current, flow_on_branch); + } } void PFSolverPowerPolar::calculateNodalInjection() { @@ -358,6 +396,10 @@ void PFSolverPowerPolar::calculateNodalInjection() { trafo->storeNodalInjection(sol_S_complex.coeff(node->matrixNodeIndex())); break; } + else if (std::shared_ptr trafo3W = std::dynamic_pointer_cast(comp)) { + trafo3W->storeNodalInjection(sol_S_complex.coeff(node->matrixNodeIndex())); + break; + } } } } diff --git a/models/Include/cps/Base/Base_Ph1_Transformer.h b/models/Include/cps/Base/Base_Ph1_Transformer.h index a3dcace814..a7aec7f7c3 100644 --- a/models/Include/cps/Base/Base_Ph1_Transformer.h +++ b/models/Include/cps/Base/Base_Ph1_Transformer.h @@ -36,6 +36,59 @@ namespace Ph1 { mInductance = inductance; } }; + + class Transformer3W { + protected: + /// Nominal voltage of primary side + Real mNominalVoltageEnd1; + /// Nominal voltage of secondary side + Real mNominalVoltageEnd2; + /// Nominal voltage of secondary side + Real mNominalVoltageEnd3; + /// Transformer ratio primary side + Complex mRatio1; + /// Transformer ratio secondary side + Complex mRatio2; + /// Transformer ratio tetrary side + Complex mRatio3; + /// Resistance [Ohm] primary side + Real mResistance1; + /// Resistance [Ohm] secondary side + Real mResistance2; + /// Resistance [Ohm] tetrary side + Real mResistance3; + /// Inductance [H] primary side + Real mInductance1; + /// Inductance [H] secondary side + Real mInductance2; + /// Inductance [H] tetrary side + Real mInductance3; + + public: + /// + void setParameters( + Real nomVoltageEnd1, Real nomVoltageEnd2, Real nomVoltageEnd3, + Real ratioAbs1, Real ratioAbs2, Real ratioAbs3, + Real ratioPhase1, Real ratioPhase2, Real ratioPhase3, + Real resistance1, Real resistance2, Real resistance3, + Real inductance1, Real inductance2, Real inductance3 + ) { + + mNominalVoltageEnd1 = nomVoltageEnd1; + mNominalVoltageEnd2 = nomVoltageEnd2; + mNominalVoltageEnd3 = nomVoltageEnd3; + mRatio1 = std::polar(ratioAbs1, ratioPhase1); + mRatio2 = std::polar(ratioAbs2, ratioPhase2); + mRatio3 = std::polar(ratioAbs3, ratioPhase3); + mResistance1 = resistance1; + mResistance2 = resistance2; + mResistance3 = resistance3; + mInductance1 = inductance1; + mInductance2 = inductance2; + mInductance3 = inductance3; + } + }; + } } } diff --git a/models/Include/cps/CIM/Reader.h b/models/Include/cps/CIM/Reader.h index b101ede81a..fb56723355 100644 --- a/models/Include/cps/CIM/Reader.h +++ b/models/Include/cps/CIM/Reader.h @@ -47,7 +47,9 @@ namespace CIMPP { class ExternalNetworkInjection; class EnergyConsumer; class PowerTransformer; + class PowerTransformer3W; class EquivalentShunt; + class LinearShuntCompensator; class TopologicalNode; class ConductingEquipment; }; @@ -147,6 +149,7 @@ namespace CIM { TopologicalPowerComp::Ptr mapExternalNetworkInjection(CIMPP::ExternalNetworkInjection* extnet); /// Returns a shunt TopologicalPowerComp::Ptr mapEquivalentShunt(CIMPP::EquivalentShunt *shunt); + TopologicalPowerComp::Ptr mapEquivalentShunt(CIMPP::LinearShuntCompensator *shunt); // #### Helper Functions #### /// Determine base voltage associated with object diff --git a/models/Include/cps/Components.h b/models/Include/cps/Components.h index 051b70f03c..97914e3338 100644 --- a/models/Include/cps/Components.h +++ b/models/Include/cps/Components.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include diff --git a/models/Include/cps/SP/SP_Ph1_Transformer3W.h b/models/Include/cps/SP/SP_Ph1_Transformer3W.h new file mode 100644 index 0000000000..3e27a7192d --- /dev/null +++ b/models/Include/cps/SP/SP_Ph1_Transformer3W.h @@ -0,0 +1,192 @@ +/* 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 SP { +namespace Ph1 { + /// Transformer that includes an inductance and resistance + class Transformer3W : + public SimPowerComp, + public SharedFactory, + public PFSolverInterfaceBranch, + public MNAInterface, + public Base::Ph1::Transformer3W { + + private: + /// Internal inductors to model losses + std::shared_ptr mSubInductor1; + std::shared_ptr mSubInductor2; + std::shared_ptr mSubInductor3; + /// Internal resistors to model losses + std::shared_ptr mSubResistor1; + std::shared_ptr mSubResistor2; + std::shared_ptr mSubResistor3; + /// Internal parallel resistance as snubber + std::shared_ptr mSubSnubResistor; + + /// Snubber resistance added on the low voltage side + Real mSnubberResistance; + + /// Rated Apparent Powers [VA] + Real mRatedPower1 = 0; + Real mRatedPower2 = 0; + Real mRatedPower3 = 0; + /// Transformer ratios magnitude + Real mRatioAbs1 = 1; + Real mRatioAbs2 = 1; + Real mRatioAbs3 = 1; + /// Transformer ratios phase [deg] + Real mRatioPhase1 = 0; + Real mRatioPhase2 = 0; + Real mRatioPhase3 = 0; + + /// Nominal omega ?? + Real mNominalOmega; + + /// Voltage [V] ?? + Real mSvVoltage; + + /// Conductances [S] + Real mConductance1; + Real mConductance2; + Real mConductance3; + + /// Reactances [Ohm] + Real mReactance1; + Real mReactance2; + Real mReactance3; + + + /// Magnetizing reactance [Ohm] + Real mMagnetizingReactance=1e9; + /// Leakage + Complex mLeakage1; + Complex mLeakage2; + Complex mLeakage3; + /// Magnetizing impedance + Complex mMagnetizing; + + /// base apparent power[VA] + Real mBaseApparentPower; + /// base impedance [ohm] + Real mBaseImpedance; + /// base inductance [H] + Real mBaseInductance; + /// base admittance [S] + Real mBaseAdmittance; + ///base omega [1/s] + Real mBaseOmega; + /// base voltage [V] + Real mBaseVoltage; + ///base current [A] + Real mBaseCurrent; + + /// resistance + Real mResistancePerUnit1; + /// resistance + Real mResistancePerUnit2; + /// resistance + Real mResistancePerUnit3; + /// reactance + Real mReactancePerUnit1; + /// reactance + Real mReactancePerUnit2; + /// reactance + Real mReactancePerUnit3; + + /// inductance + Real mInductancePerUnit1; + /// inductance + Real mInductancePerUnit2; + /// inductance + Real mInductancePerUnit3; + /// leakage impedance + Complex mLeakagePerUnit1; + /// leakage impedance + Complex mLeakagePerUnit2; + /// leakage impedance + Complex mLeakagePerUnit3; + /// magnetizing impedance + Complex mMagnetizingPerUnit; + + // #### Admittance matrix stamp #### + MatrixComp mY_element; + + // #### Power flow results #### + /// branch Current flow [A] + Eigen::Matrix mCurrent; + /// branch active powerflow [W], coef(0) has data from node 0, coef(1) from node 1. + Eigen::Matrix mActivePowerBranch; + /// branch reactive powerflow [Var] + Eigen::Matrix mReactivePowerBranch; + /// whether the total power injection of its from node is stored in this line + Bool mStoreNodalPowerInjection = false; + /// nodal active power injection + Real mActivePowerInjection; + /// nodal reactive power injection + Real mReactivePowerInjection; + + /// Boolean for considering resistive losses with sub resistor + Bool mWithResistiveLosses; + public: + /// Defines UID, name and logging level + Transformer3W(String uid, String name, + Logger::Level logLevel = Logger::Level::off, Bool withResistiveLosses = false); + /// Defines name and logging level + Transformer3W(String name, Logger::Level logLevel = Logger::Level::off) + : Transformer3W(name, name, logLevel) { } + + SimPowerComp::Ptr clone(String name) override; + + // #### General #### + /// Set transformer specific parameters (without rated power) + void setParameters(Real nomVoltageEnd1, Real nomVoltageEnd2, Real nomVoltageEnd3, + Real ratioAbs1, Real ratioAbs2, Real ratioAbs3, + Real ratioPhase1, Real ratioPhase2, Real ratioPhase3, + Real resistance1, Real resistance2, Real resistance3, + Real inductance1, Real inductance2, Real inductance3); + /// Set transformer specific parameters + void setParameters(Real nomVoltageEnd1, Real nomVoltageEnd2, Real nomVoltageEnd3, + Real ratedPower1, Real ratedPower2, Real ratedPower3, + Real ratioAbs1, Real ratioAbs2, Real ratioAbs3, + Real ratioPhase1, Real ratioPhase2, Real ratioPhase3, + Real resistance1, Real resistance2, Real resistance3, + Real inductance1, Real inductance2, Real inductance3); + /// Initializes component from power flow data + void initializeFromNodesAndTerminals(Real frequency) override; + + // #### Powerflow section #### + /// Set base voltage + void setBaseVoltage(Real baseVoltage); + /// Initializes component from power flow data + void calculatePerUnitParameters(Real baseApparentPower, Real baseOmega); + /// Stamps admittance matrix + void pfApplyAdmittanceMatrixStamp(SparseMatrixCompRow & Y) override; + /// updates branch current and power flow, input pu value, update with real value + void updateBranchFlow(VectorComp& current, VectorComp& powerflow); + /// stores nodal injection power in this line object + void storeNodalInjection(Complex powerInjection); + + // #### Getter #### + /// get admittance matrix + MatrixComp Y_element(); + + }; +} +} +} diff --git a/models/Source/CIM/Reader.cpp b/models/Source/CIM/Reader.cpp index 496aeaaa70..05532647fb 100644 --- a/models/Source/CIM/Reader.cpp +++ b/models/Source/CIM/Reader.cpp @@ -81,6 +81,8 @@ TopologicalPowerComp::Ptr Reader::mapComponent(BaseClass* obj) { return mapExternalNetworkInjection(extnet); if (CIMPP::EquivalentShunt *shunt = dynamic_cast(obj)) return mapEquivalentShunt(shunt); + if (CIMPP::LinearShuntCompensator *shunt = dynamic_cast(obj)) + return mapEquivalentShunt(shunt); return nullptr; } @@ -370,124 +372,254 @@ TopologicalPowerComp::Ptr Reader::mapACLineSegment(CIMPP::ACLineSegment* line) { } TopologicalPowerComp::Ptr Reader::mapPowerTransformer(CIMPP::PowerTransformer* trans) { - if (trans->PowerTransformerEnd.size() != 2) { - mSLog->warn("PowerTransformer {} does not have exactly two windings, ignoring", trans->name); + if (trans->PowerTransformerEnd.size() != 2 && trans->PowerTransformerEnd.size() != 3) { + mSLog->warn("PowerTransformer {} does not have exactly two or three windings, ignoring", trans->name); return nullptr; } mSLog->info("Found PowerTransformer {}", trans->name); // assign transformer ends - CIMPP::PowerTransformerEnd* end1 = nullptr, *end2 = nullptr; + CIMPP::PowerTransformerEnd* end1 = nullptr, *end2 = nullptr, *end3 = nullptr; for (auto end : trans->PowerTransformerEnd) { if (end->Terminal->sequenceNumber == 1) end1 = end; else if (end->Terminal->sequenceNumber == 2) end2 = end; + else if (end->Terminal->sequenceNumber == 3) end3 = end; else return nullptr; } // setting default values for non-set resistances and reactances mSLog->info(" PowerTransformerEnd_1 {}", end1->name); mSLog->info(" Srated={} Vrated={}", (float) end1->ratedS.value, (float) end1->ratedU.value); - try{ + try { mSLog->info(" R={}", (float) end1->r.value); - }catch(ReadingUninitializedField* e1){ + } catch(ReadingUninitializedField* e1){ end1->r.value = 1e-12; mSLog->warn(" Uninitialized value for PowerTrafoEnd1 setting default value of R={}", (float) end1->r.value); } - try{ + try { mSLog->info(" X={}", (float) end1->x.value); - }catch(ReadingUninitializedField* e1){ + } catch(ReadingUninitializedField* e1){ end1->x.value = 1e-12; mSLog->warn(" Uninitialized value for PowerTrafoEnd1 setting default value of X={}", (float) end1->x.value); } mSLog->info(" PowerTransformerEnd_2 {}", end2->name); mSLog->info(" Srated={} Vrated={}", (float) end2->ratedS.value, (float) end2->ratedU.value); - try{ + try { mSLog->info(" R={}", (float) end2->r.value); - }catch(ReadingUninitializedField* e1){ + } catch(ReadingUninitializedField* e1){ end2->r.value = 1e-12; mSLog->warn(" Uninitialized value for PowerTrafoEnd2 setting default value of R={}", (float) end2->r.value); } - try{ + try { mSLog->info(" X={}", (float) end2->x.value); - }catch(ReadingUninitializedField* e1){ + } catch(ReadingUninitializedField* e1){ end2->x.value = 1e-12; mSLog->warn(" Uninitialized value for PowerTrafoEnd2 setting default value of X={}", (float) end2->x.value); } - if (end1->ratedS.value != end2->ratedS.value) { + if (end3) { + mSLog->info(" PowerTransformerEnd_3 {}", end3->name); + mSLog->info(" Srated={} Vrated={}", (float) end3->ratedS.value, (float) end3->ratedU.value); + try { + mSLog->info(" R={}", (float) end3->r.value); + } catch(ReadingUninitializedField* e1){ + end3->r.value = 1e-12; + mSLog->warn(" Uninitialized value for PowerTrafoEnd3 setting default value of R={}", (float) end3->r.value); + } + try { + mSLog->info(" X={}", (float) end3->x.value); + } catch(ReadingUninitializedField* e1){ + end3->x.value = 1e-12; + mSLog->warn(" Uninitialized value for PowerTrafoEnd3 setting default value of X={}", (float) end3->x.value); + } + } + + if (!end3 && end1->ratedS.value != end2->ratedS.value) { mSLog->warn(" PowerTransformerEnds of {} come with distinct rated power values. Using rated power of PowerTransformerEnd_1.", trans->name); } - Real ratedPower = unitValue(end1->ratedS.value, UnitMultiplier::M); + Real voltageNode1 = unitValue(end1->ratedU.value, UnitMultiplier::k); Real voltageNode2 = unitValue(end2->ratedU.value, UnitMultiplier::k); + + Real ratedPower1 = unitValue(end1->ratedS.value, UnitMultiplier::M); + Real ratedPower2 = unitValue(end2->ratedS.value, UnitMultiplier::M); - Real ratioAbsNominal = voltageNode1 / voltageNode2; - Real ratioAbs = ratioAbsNominal; + if (!end3) { + Real ratioAbsNominal = voltageNode1 / voltageNode2; + Real ratioAbs = ratioAbsNominal; - // use normalStep from RatioTapChanger - if (end1->RatioTapChanger) { - ratioAbs = voltageNode1 / voltageNode2 * (1 + (end1->RatioTapChanger->normalStep - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); - } + // use normalStep from RatioTapChanger + if (end1->RatioTapChanger && !end3) { + ratioAbs = voltageNode1 / voltageNode2 * (1 + (end1->RatioTapChanger->normalStep - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); + } - // if corresponding SvTapStep available, use instead tap position from there - if (end1->RatioTapChanger) { - for (auto obj : mModel->Objects) { - auto tapStep = dynamic_cast(obj); - if (tapStep && tapStep->TapChanger == end1->RatioTapChanger) { - ratioAbs = voltageNode1 / voltageNode2 * (1 + (tapStep->position - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); + // if corresponding SvTapStep available, use instead tap position from there + if (end1->RatioTapChanger && !end3) { + for (auto obj : mModel->Objects) { + auto tapStep = dynamic_cast(obj); + if (tapStep && tapStep->TapChanger == end1->RatioTapChanger) { + ratioAbsNominal = voltageNode1 / voltageNode2 * (1 + (tapStep->position - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); + } } } - } - - // TODO: To be extracted from cim class - Real ratioPhase = 0; - // Calculate resistance and inductance referred to high voltage side - Real resistance = 0; - Real inductance = 0; - if (voltageNode1 >= voltageNode2 && abs(end1->x.value) > 1e-12) { - inductance = end1->x.value / mOmega; - resistance = end1->r.value; - } else if (voltageNode1 >= voltageNode2 && abs(end2->x.value) > 1e-12) { - inductance = end2->x.value / mOmega * std::pow(ratioAbsNominal, 2); - resistance = end2->r.value * std::pow(ratioAbsNominal, 2); - } - else if (voltageNode2 > voltageNode1 && abs(end2->x.value) > 1e-12) { - inductance = end2->x.value / mOmega; - resistance = end2->r.value; - } - else if (voltageNode2 > voltageNode1 && abs(end1->x.value) > 1e-12) { - inductance = end1->x.value / mOmega / std::pow(ratioAbsNominal, 2); - resistance = end1->r.value / std::pow(ratioAbsNominal, 2); - } + // TODO: To be extracted from cim class + Real ratioPhase = 0; + + // Calculate resistance and inductance referred to high voltage side + Real resistance = 0; + Real inductance = 0; + + if (voltageNode1 >= voltageNode2 && abs(end1->x.value) > 1e-12) { + inductance = end1->x.value / mOmega; + resistance = end1->r.value; + } else if (voltageNode1 >= voltageNode2 && abs(end2->x.value) > 1e-12) { + inductance = end2->x.value / mOmega * std::pow(ratioAbsNominal, 2); + resistance = end2->r.value * std::pow(ratioAbsNominal, 2); + } else if (voltageNode2 > voltageNode1 && abs(end2->x.value) > 1e-12) { + inductance = end2->x.value / mOmega; + resistance = end2->r.value; + } else if (voltageNode2 > voltageNode1 && abs(end1->x.value) > 1e-12) { + inductance = end1->x.value / mOmega / std::pow(ratioAbsNominal, 2); + resistance = end1->r.value / std::pow(ratioAbsNominal, 2); + } - if (mDomain == Domain::EMT) { - if (mPhase == PhaseType::ABC) { - Matrix resistance_3ph = CPS::Math::singlePhaseParameterToThreePhase(resistance); - Matrix inductance_3ph = CPS::Math::singlePhaseParameterToThreePhase(inductance); + if (mDomain == Domain::EMT) { + if (mPhase == PhaseType::ABC) { + Matrix resistance_3ph = CPS::Math::singlePhaseParameterToThreePhase(resistance); + Matrix inductance_3ph = CPS::Math::singlePhaseParameterToThreePhase(inductance); + Bool withResistiveLosses = resistance > 0; + auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel, withResistiveLosses); + transformer->setParameters(voltageNode1, voltageNode2, ratioAbs, ratioPhase, resistance_3ph, inductance_3ph); + return transformer; + } else { + mSLog->info(" Transformer for EMT not implemented yet"); + return nullptr; + } + } else if (mDomain == Domain::SP) { + auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel); + transformer->setParameters(voltageNode1, voltageNode2, ratedPower1, ratioAbs, ratioPhase, resistance, inductance); + Real baseVolt = voltageNode1 >= voltageNode2 ? voltageNode1 : voltageNode2; + transformer->setBaseVoltage(baseVolt); + return transformer; + } else { Bool withResistiveLosses = resistance > 0; - auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel, withResistiveLosses); - transformer->setParameters(voltageNode1, voltageNode2, ratioAbs, ratioPhase, resistance_3ph, inductance_3ph); + auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel, withResistiveLosses); + transformer->setParameters(voltageNode1, voltageNode2, ratioAbs, ratioPhase, resistance, inductance); return transformer; } - else - { - mSLog->info(" Transformer for EMT not implemented yet"); - return nullptr; - } - } - else if (mDomain == Domain::SP) { - auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel); - transformer->setParameters(voltageNode1, voltageNode2, ratedPower, ratioAbs, ratioPhase, resistance, inductance); - Real baseVolt = voltageNode1 >= voltageNode2 ? voltageNode1 : voltageNode2; - transformer->setBaseVoltage(baseVolt); - return transformer; } + else { - Bool withResistiveLosses = resistance > 0; - auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel, withResistiveLosses); - transformer->setParameters(voltageNode1, voltageNode2, ratioAbs, ratioPhase, resistance, inductance); - return transformer; + Real ratedPower3 = unitValue(end3->ratedS.value, UnitMultiplier::M); + Real voltageNode3 = unitValue(end3->ratedU.value, UnitMultiplier::k); + + Real ratioAbsNominal12 = voltageNode1 / voltageNode2; + Real ratioAbsNominal23 = voltageNode2 / voltageNode3; + Real ratioAbsNominal13 = voltageNode1 / voltageNode3; + + Real ratioAbs1 = 1.0; + // use normalStep from RatioTapChanger + if (end1->RatioTapChanger) { + ratioAbs1 = (1 + (end1->RatioTapChanger->normalStep - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); + } + // if corresponding SvTapStep available, use instead tap position from there + if (end1->RatioTapChanger) { + for (auto obj : mModel->Objects) { + auto tapStep = dynamic_cast(obj); + if (tapStep && tapStep->TapChanger == end1->RatioTapChanger) { + ratioAbs1 = (1 + (tapStep->position - end1->RatioTapChanger->neutralStep) * end1->RatioTapChanger->stepVoltageIncrement.value / 100); + } + } + } + + Real ratioAbs2 = 1.0; + // use normalStep from RatioTapChanger + if (end2->RatioTapChanger) { + ratioAbs2 = (1 + (end2->RatioTapChanger->normalStep - end2->RatioTapChanger->neutralStep) * end2->RatioTapChanger->stepVoltageIncrement.value / 100); + } + // if corresponding SvTapStep available, use instead tap position from there + if (end2->RatioTapChanger) { + for (auto obj : mModel->Objects) { + auto tapStep = dynamic_cast(obj); + if (tapStep && tapStep->TapChanger == end2->RatioTapChanger) { + ratioAbs2 = (1 + (tapStep->position - end2->RatioTapChanger->neutralStep) * end2->RatioTapChanger->stepVoltageIncrement.value / 100); + } + } + } + + Real ratioAbs3 = 1.0; + // use normalStep from RatioTapChanger + if (end3->RatioTapChanger) { + ratioAbs3 = (1 + (end3->RatioTapChanger->normalStep - end3->RatioTapChanger->neutralStep) * end3->RatioTapChanger->stepVoltageIncrement.value / 100); + } + // if corresponding SvTapStep available, use instead tap position from there + if (end3->RatioTapChanger) { + for (auto obj : mModel->Objects) { + auto tapStep = dynamic_cast(obj); + if (tapStep && tapStep->TapChanger == end3->RatioTapChanger) { + ratioAbs3 = (1 + (tapStep->position - end3->RatioTapChanger->neutralStep) * end3->RatioTapChanger->stepVoltageIncrement.value / 100); + } + } + } + + // TODO: To be extracted from cim class + Real ratioPhase1 = 0; + Real ratioPhase2 = 0; + Real ratioPhase3 = 0; + + // Calculate resistance and inductance referred to high voltage side + Real resistance1 = 0; + Real resistance2 = 0; + Real resistance3 = 0; + Real inductance1 = 0; + Real inductance2 = 0; + Real inductance3 = 0; + + if (voltageNode1 >= voltageNode2 && voltageNode1 >= voltageNode3) { + inductance1 = end1->x.value / mOmega; + resistance1 = end1->r.value; + inductance2 = end2->x.value / mOmega * std::pow(ratioAbsNominal12, 2); + resistance2 = end2->r.value * std::pow(ratioAbsNominal12, 2); + inductance3 = end3->x.value / mOmega * std::pow(ratioAbsNominal13, 2); + resistance3 = end3->r.value * std::pow(ratioAbsNominal13, 2); + } else if (voltageNode2 >= voltageNode1 && voltageNode2 >= voltageNode3) { + inductance2 = end2->x.value / mOmega; + resistance2 = end2->r.value; + inductance1 = end1->x.value / mOmega / std::pow(ratioAbsNominal12, 2); + resistance1 = end1->r.value / std::pow(ratioAbsNominal12, 2); + inductance3 = end3->x.value / mOmega * std::pow(ratioAbsNominal23, 2); + resistance3 = end3->r.value * std::pow(ratioAbsNominal23, 2); + } else if (voltageNode3 >= voltageNode1 && voltageNode3 >= voltageNode2) { + inductance3 = end3->x.value / mOmega; + resistance3 = end3->r.value; + inductance1 = end1->x.value / mOmega / std::pow(ratioAbsNominal13, 2); + resistance1 = end1->r.value / std::pow(ratioAbsNominal13, 2); + inductance2 = end2->x.value / mOmega / std::pow(ratioAbsNominal23, 2); + resistance2 = end2->r.value / std::pow(ratioAbsNominal23, 2); + } + + if (mDomain == Domain::EMT) { + mSLog->info(" 3-Winding transformer for EMT not implemented yet"); + return nullptr; + } else if (mDomain == Domain::SP) { + auto transformer = std::make_shared(trans->mRID, trans->name, mComponentLogLevel); + transformer->setParameters( + voltageNode1, voltageNode2, voltageNode3, + ratedPower1, ratedPower2, ratedPower3, + ratioAbs1, ratioAbs2, ratioAbs3, + ratioPhase1, ratioPhase2, ratioPhase3, + resistance1, resistance2, resistance3, + inductance1, inductance2, inductance3 + ); + Real baseVolt = voltageNode1 >= voltageNode2 ? voltageNode1 : voltageNode2 >= voltageNode3 ? voltageNode2 : voltageNode3; + transformer->setBaseVoltage(baseVolt); + return transformer; + } else { + mSLog->info(" 3-Winding transformer for DP not implemented yet"); + return nullptr; + } + } } @@ -562,7 +694,7 @@ TopologicalPowerComp::Ptr Reader::mapSynchronousMachine(CIMPP::SynchronousMachin if (syncGen->mRID == machine->mRID){ // Check whether relevant input data are set, otherwise set default values Real setPointActivePower = 0; - Real setPointVoltage = 0; + Real setPointVoltage = unitValue(machine->ratedU.value, UnitMultiplier::k); Real maximumReactivePower = 1e12; try{ setPointActivePower = unitValue(genUnit->initialP.value, UnitMultiplier::M); @@ -630,13 +762,15 @@ TopologicalPowerComp::Ptr Reader::mapExternalNetworkInjection(CIMPP::ExternalNet try { if(extnet->RegulatingControl){ mSLog->info(" Voltage set-point={}", (float) extnet->RegulatingControl->targetValue); - cpsextnet->setParameters(extnet->RegulatingControl->targetValue*baseVoltage); // assumes that value is specified in CIM data in per unit + cpsextnet->setParameters(unitValue(extnet->RegulatingControl->targetValue, UnitMultiplier::k)); // assumes that value is specified in CIM data in per unit } else { mSLog->info(" No voltage set-point defined. Using 1 per unit."); cpsextnet->setParameters(1.*baseVoltage); } } catch (ReadingUninitializedField* e ) { - std::cerr << "Ignore incomplete RegulatingControl" << std::endl; + std::cerr << "Ignore incomplete RegulatingControl." << std::endl; + mSLog->info(" Incomplete RegulatingControl. Using 1 per unit."); + cpsextnet->setParameters(1.*baseVoltage); } return cpsextnet; @@ -666,6 +800,17 @@ TopologicalPowerComp::Ptr Reader::mapEquivalentShunt(CIMPP::EquivalentShunt* shu return cpsShunt; } +TopologicalPowerComp::Ptr Reader::mapEquivalentShunt(CIMPP::LinearShuntCompensator* shunt){ + mSLog->info("Found linear shunt compensator {}", shunt->name); + + Real baseVoltage = determineBaseVoltageAssociatedWithEquipment(shunt); + + auto cpsShunt = std::make_shared(shunt->mRID, shunt->name, mComponentLogLevel); + cpsShunt->setParameters(shunt->gPerSection.value, shunt->bPerSection.value); + cpsShunt->setBaseVoltage(baseVoltage); + return cpsShunt; +} + void Reader::initDynamicSystemTopologyWithPowerflow(SystemTopology& systemPF, SystemTopology& system) { for (auto nodePF : systemPF.mNodes) { if (auto node = system.node(nodePF->name())) { diff --git a/models/Source/CMakeLists.txt b/models/Source/CMakeLists.txt index fbcc08cdbd..e720df9b03 100644 --- a/models/Source/CMakeLists.txt +++ b/models/Source/CMakeLists.txt @@ -90,6 +90,7 @@ list(APPEND CPS_SOURCES SP/SP_Ph1_Load.cpp SP/SP_Ph1_Switch.cpp SP/SP_Ph1_Transformer.cpp + SP/SP_Ph1_Transformer3W.cpp SP/SP_Ph1_SolidStateTransformer.cpp SP/SP_Ph1_Shunt.cpp SP/SP_Ph1_SynchronGenerator.cpp diff --git a/models/Source/SP/SP_Ph1_NetworkInjection.cpp b/models/Source/SP/SP_Ph1_NetworkInjection.cpp index 4eebdd4ad2..46b8dff823 100644 --- a/models/Source/SP/SP_Ph1_NetworkInjection.cpp +++ b/models/Source/SP/SP_Ph1_NetworkInjection.cpp @@ -28,6 +28,7 @@ SP::Ph1::NetworkInjection::NetworkInjection(String uid, String name, for (auto subcomp: mSubComponents) mSLog->info("- {}", subcomp->name()); + addAttribute("base_Voltage", &mBaseVoltage, Flags::read | Flags::write); addAttribute("V_set", &mVoltageSetPoint, Flags::read | Flags::write); addAttribute("V_set_pu", &mVoltageSetPointPerUnit, Flags::read | Flags::write); addAttribute("p_inj", &mActivePowerInjection, Flags::read | Flags::write); diff --git a/models/Source/SP/SP_Ph1_Shunt.cpp b/models/Source/SP/SP_Ph1_Shunt.cpp index 9d4910ba0e..5480a5aa8f 100644 --- a/models/Source/SP/SP_Ph1_Shunt.cpp +++ b/models/Source/SP/SP_Ph1_Shunt.cpp @@ -16,6 +16,7 @@ SP::Ph1::Shunt::Shunt(String uid, String name, Logger::Level logLevel) mSLog->info("Create {} of type {}", this->type(), name); setTerminalNumber(1); + addAttribute("base_Voltage", &mBaseVoltage, Flags::read | Flags::write); addAttribute("G", &mConductance, Flags::write | Flags::read); addAttribute("B", &mSusceptance, Flags::write | Flags::read); } diff --git a/models/Source/SP/SP_Ph1_Transformer3W.cpp b/models/Source/SP/SP_Ph1_Transformer3W.cpp new file mode 100644 index 0000000000..8fb561202e --- /dev/null +++ b/models/Source/SP/SP_Ph1_Transformer3W.cpp @@ -0,0 +1,408 @@ +/* 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; + +// #### General #### +SP::Ph1::Transformer3W::Transformer3W(String uid, String name, Logger::Level logLevel, Bool withResistiveLosses) + : SimPowerComp(uid, name, logLevel) { + if (withResistiveLosses) + setVirtualNodeNumber(7); + else + setVirtualNodeNumber(4); + + mSLog->info("Create {} {}", this->type(), name); + mIntfVoltage = MatrixComp::Zero(3, 1); + mIntfCurrent = MatrixComp::Zero(3, 1); + setTerminalNumber(3); + + addAttribute("nominal_voltage_end1", &mNominalVoltageEnd1, Flags::read | Flags::write); + addAttribute("nominal_voltage_end2", &mNominalVoltageEnd2, Flags::read | Flags::write); + addAttribute("nominal_voltage_end3", &mNominalVoltageEnd3, Flags::read | Flags::write); + addAttribute("base_voltage", &mBaseVoltage, Flags::read | Flags::write); + addAttribute("S1", &mRatedPower1, Flags::write | Flags::read); + addAttribute("S2", &mRatedPower2, Flags::write | Flags::read); + addAttribute("S3", &mRatedPower3, Flags::write | Flags::read); + addAttribute("ratio1", &mRatio1, Flags::write | Flags::read); + addAttribute("ratio2", &mRatio2, Flags::write | Flags::read); + addAttribute("ratio3", &mRatio3, Flags::write | Flags::read); + addAttribute("R1", &mResistance1, Flags::write | Flags::read); + addAttribute("R2", &mResistance2, Flags::write | Flags::read); + addAttribute("R3", &mResistance3, Flags::write | Flags::read); + addAttribute("L1", &mInductance1, Flags::write | Flags::read); + addAttribute("L2", &mInductance2, Flags::write | Flags::read); + addAttribute("L3", &mInductance3, Flags::write | Flags::read); + addAttribute("nodal_injection_stored", &mStoreNodalPowerInjection, Flags::read | Flags::write); + addAttribute("p_inj", &mActivePowerInjection, Flags::read | Flags::write); + addAttribute("q_inj", &mReactivePowerInjection, Flags::read | Flags::write); + addAttribute("current_1", &mCurrent(0), Flags::read | Flags::write); + addAttribute("current_2", &mCurrent(1), Flags::read | Flags::write); + addAttribute("current_3", &mCurrent(2), Flags::read | Flags::write); + addAttribute("p_branch_1", &mActivePowerBranch(0), Flags::read | Flags::write); + addAttribute("q_branch_1", &mReactivePowerBranch(0), Flags::read | Flags::write); + addAttribute("p_branch_2", &mActivePowerBranch(1), Flags::read | Flags::write); + addAttribute("q_branch_2", &mReactivePowerBranch(1), Flags::read | Flags::write); + addAttribute("p_branch_3", &mActivePowerBranch(2), Flags::read | Flags::write); + addAttribute("q_branch_3", &mReactivePowerBranch(2), Flags::read | Flags::write); +} + + +void SP::Ph1::Transformer3W::setParameters( + Real nomVoltageEnd1, Real nomVoltageEnd2, Real nomVoltageEnd3, + Real ratioAbs1, Real ratioAbs2, Real ratioAbs3, + Real ratioPhase1, Real ratioPhase2, Real ratioPhase3, + Real resistance1, Real resistance2, Real resistance3, + Real inductance1, Real inductance2, Real inductance3 + ) { + + // Note: to be consistent impedance values must be referred to high voltage side (and base voltage set to higher voltage) + Base::Ph1::Transformer3W::setParameters( + nomVoltageEnd1, nomVoltageEnd2, nomVoltageEnd3, + ratioAbs1, ratioAbs2, ratioAbs3, + ratioPhase1, ratioPhase2, ratioPhase3, + resistance1, resistance2, resistance3, + inductance1, inductance2, inductance3 + ); + + mSLog->info("Nominal Voltage End 1={} [V] Nominal Voltage End 2={} [V] Nominal Voltage End 2={} [V]", + mNominalVoltageEnd1, mNominalVoltageEnd2, mNominalVoltageEnd3); + mSLog->info("Resistance={} [Ohm] Inductance={} [H] (primary side)", mResistance1, mInductance1); + mSLog->info("Resistance={} [Ohm] Inductance={} [H] (secondary side)", mResistance2, mInductance2); + mSLog->info("Resistance={} [Ohm] Inductance={} [H] (tetrary side)", mResistance3, mInductance3); + mSLog->info("Tap Ratio 1={} [/] Phase Shift 1={} [deg]", std::abs(mRatio1), std::arg(mRatio1)); + mSLog->info("Tap Ratio 2={} [/] Phase Shift 2={} [deg]", std::abs(mRatio2), std::arg(mRatio2)); + mSLog->info("Tap Ratio 3={} [/] Phase Shift 3={} [deg]", std::abs(mRatio3), std::arg(mRatio3)); + + mRatioAbs1 = std::abs(mRatio1); + mRatioAbs2 = std::abs(mRatio2); + mRatioAbs3 = std::abs(mRatio3); + mRatioPhase1 = std::arg(mRatio1); + mRatioPhase2 = std::arg(mRatio2); + mRatioPhase3 = std::arg(mRatio3); + mConductance1 = 1 / mResistance1; + mConductance2 = 1 / mResistance2; + mConductance3 = 1 / mResistance3; + + mParametersSet = true; +} + +void SP::Ph1::Transformer3W::setParameters( + Real nomVoltageEnd1, Real nomVoltageEnd2, Real nomVoltageEnd3, + Real ratedPower1, Real ratedPower2, Real ratedPower3, + Real ratioAbs1, Real ratioAbs2, Real ratioAbs3, + Real ratioPhase1, Real ratioPhase2, Real ratioPhase3, + Real resistance1, Real resistance2, Real resistance3, + Real inductance1, Real inductance2, Real inductance3 + ) { + + mRatedPower1 = ratedPower1; + mRatedPower2 = ratedPower2; + mRatedPower3 = ratedPower3; + mSLog->info("Rated Power ={} [W]", mRatedPower1); + mSLog->info("Rated Power ={} [W]", mRatedPower2); + mSLog->info("Rated Power ={} [W]", mRatedPower3); + + SP::Ph1::Transformer3W::setParameters( + nomVoltageEnd1, nomVoltageEnd2, nomVoltageEnd3, + ratioAbs1, ratioAbs2, ratioAbs3, + ratioPhase1, ratioPhase2, ratioPhase3, + resistance1, resistance2, resistance3, + inductance1, inductance2, inductance3); +} + + +SimPowerComp::Ptr SP::Ph1::Transformer3W::clone(String name) { + auto copy = Transformer3W::make(name, mLogLevel); + copy->setParameters( + mNominalVoltageEnd1, mNominalVoltageEnd2, mNominalVoltageEnd3, + mRatedPower1, mRatedPower2, mRatedPower3, + std::abs(mRatio1), std::abs(mRatio2), std::abs(mRatio3), + std::arg(mRatio1), std::arg(mRatio2), std::arg(mRatio3), + mResistance1, mResistance2, mResistance3, + mInductance1, mInductance2, mInductance3 + ); + return copy; +} + +void SP::Ph1::Transformer3W::initializeFromNodesAndTerminals(Real frequency) { + mNominalOmega = 2. * PI * frequency; + mReactance1 = mNominalOmega * mInductance1; + mReactance2 = mNominalOmega * mInductance2; + mReactance3 = mNominalOmega * mInductance3; + mSLog->info("Reactance={} [Ohm] (primary side)", mReactance1); + mSLog->info("Reactance={} [Ohm] (secondary side)", mReactance2); + mSLog->info("Reactance={} [Ohm] (tetrary side)", mReactance3); + + + // Component parameters are referred to high voltage side. + // Switch terminals if transformer is connected the other way around. + // Without modifying the mNominalVoltageEnd value. + if (mNominalVoltageEnd1 < mNominalVoltageEnd2 && + mNominalVoltageEnd1 < mNominalVoltageEnd3) { + std::shared_ptr> tmp = mTerminals[0]; + if (mNominalVoltageEnd2 > mNominalVoltageEnd3) { + mTerminals[0] = mTerminals[1]; + mTerminals[1] = mTerminals[2]; + } + else { + mTerminals[0] = mTerminals[2]; + } + mTerminals[2] = tmp; + } + else if (mNominalVoltageEnd1 < mNominalVoltageEnd2) { + std::shared_ptr> tmp = mTerminals[0]; + mTerminals[0] = mTerminals[1]; + mTerminals[1] = tmp; + } + else if (mNominalVoltageEnd1 < mNominalVoltageEnd3) { + std::shared_ptr> tmp = mTerminals[0]; + mTerminals[0] = mTerminals[2]; + mTerminals[2] = mTerminals[1]; + mTerminals[1] = tmp; + } + else if (mNominalVoltageEnd2 < mNominalVoltageEnd3) { + std::shared_ptr> tmp = mTerminals[1]; + mTerminals[1] = mTerminals[2]; + mTerminals[2] = tmp; + } + + // TODO: Temporary for some reason its resetting to 1x1 the matrix ignoring L22 and L23 + mIntfVoltage = MatrixComp::Zero(3, 1); + mIntfCurrent = MatrixComp::Zero(3, 1); + + // Static calculations from load flow data + Complex impedance1 = { mResistance1, mReactance1 }; + Complex impedance2 = { mResistance2, mReactance2 }; + Complex impedance3 = { mResistance3, mReactance3 }; + + // Set initial voltage of virtual node in between + mVirtualNodes[1]->setInitialVoltage(initialSingleVoltage(0)); + mVirtualNodes[2]->setInitialVoltage(initialSingleVoltage(1)); + mVirtualNodes[3]->setInitialVoltage(initialSingleVoltage(2)); + + mIntfVoltage(0, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(0) * mRatio1; + mIntfVoltage(1, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(1) * mRatio2; + mIntfVoltage(2, 0) = mVirtualNodes[0]->initialSingleVoltage() - initialSingleVoltage(2) * mRatio3; + + mIntfCurrent(0, 0) = mIntfVoltage(0, 0) / impedance1; + mIntfCurrent(1, 0) = mIntfVoltage(1, 0) / impedance2; + mIntfCurrent(2, 0) = mIntfVoltage(2, 0) / impedance3; + + // Create series sub components + mSubInductor1 = std::make_shared(mUID + "_ind1", mName + "_ind1", Logger::Level::off); + mSubInductor1->setParameters(mInductance1); + mSubInductor2 = std::make_shared(mUID + "_ind2", mName + "_ind2", Logger::Level::off); + mSubInductor2->setParameters(mInductance2); + mSubInductor3 = std::make_shared(mUID + "_ind3", mName + "_ind3", Logger::Level::off); + mSubInductor3->setParameters(mInductance3); + + + if (mNumVirtualNodes == 7) { + mVirtualNodes[4]->setInitialVoltage(initialSingleVoltage(0)); + mSubResistor1 = std::make_shared(mUID + "_res1", mName + "_res1", Logger::Level::off); + mSubResistor1->setParameters(mResistance1); + mSubResistor1->connect({ node(0), mVirtualNodes[4] }); + mSubResistor1->initialize(mFrequencies); + mSubResistor1->initializeFromNodesAndTerminals(frequency); + mSubInductor1->connect({ mVirtualNodes[4], mVirtualNodes[1] }); + + mVirtualNodes[5]->setInitialVoltage(initialSingleVoltage(1)); + mSubResistor2 = std::make_shared(mUID + "_res2", mName + "_res2", Logger::Level::off); + mSubResistor2->setParameters(mResistance1); + mSubResistor2->connect({ node(1), mVirtualNodes[5] }); + mSubResistor2->initialize(mFrequencies); + mSubResistor2->initializeFromNodesAndTerminals(frequency); + mSubInductor2->connect({ mVirtualNodes[5], mVirtualNodes[2] }); + + mVirtualNodes[6]->setInitialVoltage(initialSingleVoltage(2)); + mSubResistor2 = std::make_shared(mUID + "_res3", mName + "_res3", Logger::Level::off); + mSubResistor2->setParameters(mResistance1); + mSubResistor2->connect({ node(2), mVirtualNodes[6] }); + mSubResistor2->initialize(mFrequencies); + mSubResistor2->initializeFromNodesAndTerminals(frequency); + mSubInductor2->connect({ mVirtualNodes[6], mVirtualNodes[3] }); + + } else { + mSubInductor1->connect({ node(0), mVirtualNodes[1] }); + mSubInductor2->connect({ node(1), mVirtualNodes[2] }); + mSubInductor3->connect({ node(2), mVirtualNodes[3] }); + } + mSubInductor1->initialize(mFrequencies); + mSubInductor1->initializeFromNodesAndTerminals(frequency); + mSubInductor2->initialize(mFrequencies); + mSubInductor2->initializeFromNodesAndTerminals(frequency); + mSubInductor3->initialize(mFrequencies); + mSubInductor3->initializeFromNodesAndTerminals(frequency); + + // Create parallel sub components + // A snubber conductance is added on the low voltage side (resistance scaled with lower voltage side) + mSnubberResistance = mNominalVoltageEnd1 > mNominalVoltageEnd2 ? std::abs(mNominalVoltageEnd2)*1e6 : std::abs(mNominalVoltageEnd1)*1e6; + mSubSnubResistor = std::make_shared(mUID + "_snub_res", mName + "_snub_res", Logger::Level::off); + mSubSnubResistor->setParameters(mSnubberResistance); + mSubSnubResistor->connect({ node(2), SP::SimNode::GND }); + mSubSnubResistor->initialize(mFrequencies); + mSubSnubResistor->initializeFromNodesAndTerminals(frequency); + mSLog->info("Snubber Resistance={} [Ohm] (connected to LV side)", mSnubberResistance); + + mSLog->info( + "\n--- Initialization from powerflow ---" + "\nVoltage across 1: {:s}" + "\nCurrent 1: {:s}" + "\nVoltage across 2: {:s}" + "\nCurrent 2: {:s}" + "\nVoltage across 3: {:s}" + "\nCurrent 3: {:s}" + "\nTerminal 1 voltage: {:s}" + "\nTerminal 2 voltage: {:s}" + "\nTerminal 3 voltage: {:s}" + "\n--- Initialization from powerflow finished ---", + Logger::phasorToString(mIntfVoltage(0, 0)), + Logger::phasorToString(mIntfCurrent(0, 0)), + Logger::phasorToString(mIntfVoltage(1, 0)), + Logger::phasorToString(mIntfCurrent(1, 0)), + Logger::phasorToString(mIntfVoltage(2, 0)), + Logger::phasorToString(mIntfCurrent(2, 0)), + Logger::phasorToString(initialSingleVoltage(0)), + Logger::phasorToString(initialSingleVoltage(1)), + Logger::phasorToString(initialSingleVoltage(2))); +} + + +// #### Powerflow section #### + +void SP::Ph1::Transformer3W::setBaseVoltage(Real baseVoltage) { + // Note: to be consistent set base voltage to higher voltage (and impedance values must be referred to high voltage side) + // TODO: use attribute setter for setting base voltage + mBaseVoltage = baseVoltage; +} + +void SP::Ph1::Transformer3W::calculatePerUnitParameters(Real baseApparentPower, Real baseOmega) { + mSLog->info("#### Calculate Per Unit Parameters for {}", mName); + mBaseApparentPower = baseApparentPower; + mBaseOmega = baseOmega; + mSLog->info("Base Power={} [VA] Base Omega={} [1/s]", baseApparentPower, baseOmega); + + mBaseImpedance = mBaseVoltage * mBaseVoltage / mBaseApparentPower; + mBaseAdmittance = 1.0 / mBaseImpedance; + mBaseCurrent = baseApparentPower / (mBaseVoltage*sqrt(3)); // I_base=(S_threephase/3)/(V_line_to_line/sqrt(3)) + mSLog->info("Base Voltage={} [V] Base Impedance={} [Ohm]", mBaseVoltage, mBaseImpedance); + + mResistancePerUnit1 = mResistance1 / mBaseImpedance; + mReactancePerUnit1 = mReactance1 / mBaseImpedance; + mSLog->info("Resistance={} [pu] Reactance={} [pu]", mResistancePerUnit1, mReactancePerUnit1); + + mResistancePerUnit2 = mResistance2 / mBaseImpedance; + mReactancePerUnit2 = mReactance2 / mBaseImpedance; + mSLog->info("Resistance={} [pu] Reactance={} [pu]", mResistancePerUnit2, mReactancePerUnit2); + + mResistancePerUnit3 = mResistance3 / mBaseImpedance; + mReactancePerUnit3 = mReactance3 / mBaseImpedance; + mSLog->info("Resistance={} [pu] Reactance={} [pu]", mResistancePerUnit3, mReactancePerUnit3); + + mBaseInductance = mBaseImpedance / mBaseOmega; + mInductancePerUnit1 = mInductance1 / mBaseInductance; + mInductancePerUnit2 = mInductance2 / mBaseInductance; + mInductancePerUnit3 = mInductance3 / mBaseInductance; + mMagnetizing = Complex(0., mMagnetizingReactance); + // omega per unit=1, hence 1.0*mInductancePerUnit. + mLeakagePerUnit1 = Complex(mResistancePerUnit1, 1.*mInductancePerUnit1); + mLeakagePerUnit2 = Complex(mResistancePerUnit2, 1.*mInductancePerUnit2); + mLeakagePerUnit3 = Complex(mResistancePerUnit3, 1.*mInductancePerUnit3); + mMagnetizingPerUnit = mMagnetizing / mBaseImpedance; + + // set snubber resistance + if (mBehaviour == Behaviour::Initialization) { + Real snubberResistance = 1e3; + mSubSnubResistor = std::make_shared(mUID + "_snub_res", mName + "_snub_res", mLogLevel); + mSubSnubResistor->setParameters(snubberResistance); + mSubSnubResistor->connect({ node(2), SP::SimNode::GND }); + mSubSnubResistor->initializeFromNodesAndTerminals(mBaseOmega); + } + if (mSubSnubResistor) { + mSubSnubResistor->setBaseVoltage(mBaseVoltage); + mSubSnubResistor->calculatePerUnitParameters(baseApparentPower); + } + + mSLog->info("Leakage Impedance 12 Per Unit={} [pu] ", mLeakagePerUnit1 + mLeakagePerUnit2); + mSLog->info("Leakage Impedance 23 Per Unit={} [pu] ", mLeakagePerUnit2 + mLeakagePerUnit3); + mSLog->info("Leakage Impedance 31 Per Unit={} [pu] ", mLeakagePerUnit3 + mLeakagePerUnit1); +} + +void SP::Ph1::Transformer3W::pfApplyAdmittanceMatrixStamp(SparseMatrixCompRow & Y) { + // calculate matrix stamp + mY_element = MatrixComp(3, 3); + Complex a, b, c, d; + + a = mLeakagePerUnit2 * mLeakagePerUnit3 * mRatioAbs1; + b = mLeakagePerUnit3 * mLeakagePerUnit1 * mRatioAbs2; + c = mLeakagePerUnit1 * mLeakagePerUnit2 * mRatioAbs3; + + if ((std::abs(mMagnetizingPerUnit) > 1e9)) { + d = a * mRatioAbs1 + b * mRatioAbs2 + c * mRatioAbs3; + } + else { + d = a * mRatioAbs1 + b * mRatioAbs2 + c * mRatioAbs3 + + (mLeakagePerUnit1 * mLeakagePerUnit2 * mLeakagePerUnit3 / mMagnetizingPerUnit); + } + + mY_element(0, 0) = (d - mRatioAbs1 * a) / (mLeakagePerUnit1 * d); + mY_element(0, 1) = - mRatioAbs1 * b / (mLeakagePerUnit1 * d); + mY_element(0, 2) = - mRatioAbs1 * c / (mLeakagePerUnit1 * d); + + mY_element(1, 0) = - mRatioAbs2 * a / (mLeakagePerUnit2 * d); + mY_element(1, 1) = (d - mRatioAbs2 * b) / (mLeakagePerUnit2 * d); + mY_element(1, 2) = - mRatioAbs2 * c / (mLeakagePerUnit2 * d); + + mY_element(2, 0) = - mRatioAbs3 * a / (mLeakagePerUnit3 * d); + mY_element(2, 1) = - mRatioAbs3 * b / (mLeakagePerUnit3 * d); + mY_element(2, 2) = (d - mRatioAbs3 * c) / (mLeakagePerUnit3 * d); + + //check for inf or nan + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + if (std::isinf(mY_element.coeff(i, j).real()) || std::isinf(mY_element.coeff(i, j).imag())) { + std::cout << mY_element << std::endl; + std::cout << "Zm:" << mMagnetizing << std::endl; + std::stringstream ss; + ss << "Transformer>>" << this->name() << ": infinite or nan values in the element Y at: " << i << "," << j; + throw std::invalid_argument(ss.str()); + } + + //set the circuit matrix values + Y.coeffRef(this->matrixNodeIndex(0), this->matrixNodeIndex(0)) += mY_element.coeff(0, 0); + Y.coeffRef(this->matrixNodeIndex(0), this->matrixNodeIndex(1)) += mY_element.coeff(0, 1); + Y.coeffRef(this->matrixNodeIndex(0), this->matrixNodeIndex(2)) += mY_element.coeff(0, 2); + + Y.coeffRef(this->matrixNodeIndex(1), this->matrixNodeIndex(0)) += mY_element.coeff(1, 0); + Y.coeffRef(this->matrixNodeIndex(1), this->matrixNodeIndex(1)) += mY_element.coeff(1, 1); + Y.coeffRef(this->matrixNodeIndex(1), this->matrixNodeIndex(2)) += mY_element.coeff(1, 2); + + Y.coeffRef(this->matrixNodeIndex(2), this->matrixNodeIndex(0)) += mY_element.coeff(2, 0); + Y.coeffRef(this->matrixNodeIndex(2), this->matrixNodeIndex(1)) += mY_element.coeff(2, 1); + Y.coeffRef(this->matrixNodeIndex(2), this->matrixNodeIndex(2)) += mY_element.coeff(2, 2); + + mSLog->info("#### Y matrix stamping: {}", mY_element); +} + +void SP::Ph1::Transformer3W::updateBranchFlow(VectorComp& current, VectorComp& powerflow) { + mCurrent = current * mBaseCurrent; + mActivePowerBranch = powerflow.real()*mBaseApparentPower; + mReactivePowerBranch = powerflow.imag()*mBaseApparentPower; +} + +void SP::Ph1::Transformer3W::storeNodalInjection(Complex powerInjection) { + mActivePowerInjection = std::real(powerInjection)*mBaseApparentPower; + mReactivePowerInjection = std::imag(powerInjection)*mBaseApparentPower; + mStoreNodalPowerInjection = true; +} + +MatrixComp SP::Ph1::Transformer3W::Y_element() { + return mY_element; +}