From ecfd366ff08f123cb2e71b5ee6b65250b8aa7948 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20LAIGRE?= Date: Fri, 29 Apr 2022 18:54:03 +0200 Subject: [PATCH] Add Twtdata and unit tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Sébastien LAIGRE --- include/powsybl/iidm/util/LinkData.hpp | 23 + include/powsybl/iidm/util/TwtData.hpp | 211 ++++++++ src/CMakeLists.txt | 1 + src/iidm/util/LinkData.cpp | 110 +++++ src/iidm/util/TwtData.cpp | 610 ++++++++++++++++++++++++ test/iidm/CMakeLists.txt | 1 + test/iidm/util/DanglingLineDataTest.cpp | 5 +- test/iidm/util/TwtDataTest.cpp | 496 +++++++++++++++++++ 8 files changed, 1453 insertions(+), 4 deletions(-) create mode 100644 include/powsybl/iidm/util/TwtData.hpp create mode 100644 src/iidm/util/TwtData.cpp create mode 100644 test/iidm/util/TwtDataTest.cpp diff --git a/include/powsybl/iidm/util/LinkData.hpp b/include/powsybl/iidm/util/LinkData.hpp index 400bf35b..59cb3897 100644 --- a/include/powsybl/iidm/util/LinkData.hpp +++ b/include/powsybl/iidm/util/LinkData.hpp @@ -10,6 +10,8 @@ #include +#include + namespace powsybl { namespace iidm { @@ -25,6 +27,10 @@ struct Flow { struct BranchAdmittanceMatrix { public: + BranchAdmittanceMatrix() = default; + + BranchAdmittanceMatrix(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22); + std::complex y11; std::complex y12; @@ -37,6 +43,23 @@ struct BranchAdmittanceMatrix { BranchAdmittanceMatrix calculateBranchAdmittance(double r, double x, double ratio1, double angle1, double ratio2, double angle2, const std::complex& ysh1, const std::complex& ysh2); +double getFixedX(double x, double epsilonX, bool applyReactanceCorrection); + +Flow flowBothEnds(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, + double u1, double theta1, double u2, double theta2); + +Flow flowBothEnds(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, + const std::complex& v1, const std::complex& v2); + +std::complex flowYshunt(const std::complex& ysh, double u, double theta); + +double getPhaseAngleClockDegrees(int phaseAngleClock); + +BranchAdmittanceMatrix kronChain(const BranchAdmittanceMatrix& firstAdm, const Branch::Side& firstChainNodeSide, + const BranchAdmittanceMatrix& secondAdm, const Branch::Side& secondChainNodeSide); + +std::complex kronAntenna(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, bool isOpenFrom); + } // namespace LinkData } // namespace iidm diff --git a/include/powsybl/iidm/util/TwtData.hpp b/include/powsybl/iidm/util/TwtData.hpp new file mode 100644 index 00000000..ca2dd193 --- /dev/null +++ b/include/powsybl/iidm/util/TwtData.hpp @@ -0,0 +1,211 @@ +/** + * Copyright (c) 2022, RTE (http://www.rte-france.com) + * 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 http://mozilla.org/MPL/2.0/. + */ + +#ifndef POWSYBL_IIDM_UTIL_TWTDATA_HPP +#define POWSYBL_IIDM_UTIL_TWTDATA_HPP + +#include +#include + +#include +#include + +namespace powsybl { + +namespace iidm { + +class ThreeWindingsTransformer; + +class TwtData { +public: + TwtData(const ThreeWindingsTransformer& twt, double epsilonX, bool applyReactanceCorrection); + + TwtData(const ThreeWindingsTransformer& twt, double epsilonX, bool applyReactanceCorrection, + bool twtSplitShuntAdmittance); + + TwtData(const ThreeWindingsTransformer& twt, int phaseAngleClock2, int phaseAngleClock3, double epsilonX, + bool applyReactanceCorrection, bool twtSplitShuntAdmittance); + + double getB1(const ThreeWindingsTransformer::Side& side) const; + + double getB2(const ThreeWindingsTransformer::Side& side) const; + + double getComputedP(const ThreeWindingsTransformer::Side& side) const; + + double getComputedQ(const ThreeWindingsTransformer::Side& side) const; + + double getG1(const ThreeWindingsTransformer::Side& side) const; + + double getG2(const ThreeWindingsTransformer::Side& side) const; + + const std::string& getId() const; + + double getP(const ThreeWindingsTransformer::Side& side) const; + + double getQ(const ThreeWindingsTransformer::Side& side) const; + + double getR(const ThreeWindingsTransformer::Side& side) const; + + double getRatedU(const ThreeWindingsTransformer::Side& side) const; + + double getStarTheta() const; + + double getStarU() const; + + double getTheta(const ThreeWindingsTransformer::Side& side) const; + + double getU(const ThreeWindingsTransformer::Side& side) const; + + double getX(const ThreeWindingsTransformer::Side& side) const; + + bool isConnected(const ThreeWindingsTransformer::Side& side) const; + + bool isMainComponent(const ThreeWindingsTransformer::Side& side) const; + + int getPhaseAngleClock2() const; + + int getPhaseAngleClock3() const; + + double getRatedU0() const; + +private: + static double alpha(const ThreeWindingsTransformer::Leg& leg); + + static double getB1(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance); + + static double getB2(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance); + + static double getG1(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance); + + static double getG2(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance); + + static double getR(const ThreeWindingsTransformer::Leg& leg); + + static double getTheta(const ThreeWindingsTransformer::Leg& leg); + + static double getV(const ThreeWindingsTransformer::Leg& leg); + + static double getValue(double initialValue, double rtcStepValue, double ptcStepValue); + + static double getX(const ThreeWindingsTransformer::Leg& leg); + + static bool isMainComponent(const ThreeWindingsTransformer::Leg& leg); + + static double rho(const ThreeWindingsTransformer::Leg& leg, double ratedU0); + + static bool valid(double voltage, double theta); + +private: + std::complex calculateOneConnectedLegFlow(double u, double theta, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg) const; + + std::complex calculateOneConnectedLegShunt(const LinkData::BranchAdmittanceMatrix& closeLeg, + const LinkData::BranchAdmittanceMatrix& firstOpenLeg, + const LinkData::BranchAdmittanceMatrix& secondOpenLeg) const; + + std::complex calculateOneConnectedLegStarBusVoltage(double u, double theta, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg) const; + + void calculateThreeConnectedLegsFlowAndStarBusVoltage(double u1, double theta1, double u2, double theta2, + double u3, double theta3, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg1, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg2, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg3); + + LinkData::BranchAdmittanceMatrix calculateTwoConnectedLegsAdmittance(const LinkData::BranchAdmittanceMatrix& firstCloseLeg, + const LinkData::BranchAdmittanceMatrix& secondCloseLeg, + const LinkData::BranchAdmittanceMatrix& openLeg) const; + + LinkData::Flow calculateTwoConnectedLegsFlow(double u1, double theta1, double u2, double theta2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg) const; + + std::complex calculateTwoConnectedLegsStarBusVoltage(double u1, double theta1, double u2, double theta2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg) const; + +private: + std::string m_id; + + double m_p1; + double m_q1; + double m_p2; + double m_q2; + double m_p3; + double m_q3; + + double m_u1; + double m_theta1; + double m_u2; + double m_theta2; + double m_u3; + double m_theta3; + + double m_r1; + double m_x1; + double m_r2; + double m_x2; + double m_r3; + double m_x3; + + double m_g11; + double m_b11; + double m_g12; + double m_b12; + double m_g21; + double m_b21; + double m_g22; + double m_b22; + double m_g31; + double m_b31; + double m_g32; + double m_b32; + + double m_rho1; + double m_alpha1; + double m_rho2; + double m_alpha2; + double m_rho3; + double m_alpha3; + + double m_ratedU1; + double m_ratedU2; + double m_ratedU3; + + bool m_connected1; + bool m_connected2; + bool m_connected3; + bool m_mainComponent1; + bool m_mainComponent2; + bool m_mainComponent3; + + double m_computedP1; + double m_computedQ1; + double m_computedP2; + double m_computedQ2; + double m_computedP3; + double m_computedQ3; + + double m_starU; + double m_starTheta; + + int m_phaseAngleClock2; + int m_phaseAngleClock3; + double m_ratedU0; +}; + +} // namespace iidm + +} // namespace powsybl + +#endif // POWSYBL_IIDM_UTIL_TWTDATA_HPP diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cb9da53b..6e93ced1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -209,6 +209,7 @@ set(IIDM_SOURCES iidm/util/Substations.cpp iidm/util/SV.cpp iidm/util/TerminalFinder.cpp + iidm/util/TwtData.cpp iidm/util/VoltageLevels.cpp logging/ConsoleLogger.cpp diff --git a/src/iidm/util/LinkData.cpp b/src/iidm/util/LinkData.cpp index 0a7478d3..caf5c205 100644 --- a/src/iidm/util/LinkData.cpp +++ b/src/iidm/util/LinkData.cpp @@ -7,12 +7,21 @@ #include +#include + namespace powsybl { namespace iidm { namespace LinkData { +BranchAdmittanceMatrix::BranchAdmittanceMatrix(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22) : + y11(y11), + y12(y12), + y21(y21), + y22(y22) { +} + BranchAdmittanceMatrix calculateBranchAdmittance(double r, double x, double ratio1, double angle1, double ratio2, double angle2, const std::complex& ysh1, const std::complex& ysh2) { std::complex a1 = std::polar(ratio1, angle1); @@ -33,6 +42,107 @@ BranchAdmittanceMatrix calculateBranchAdmittance(double r, double x, double rati return branchAdmittance; } +double getFixedX(double x, double epsilonX, bool applyReactanceCorrection) { + return std::abs(x) < epsilonX && applyReactanceCorrection ? epsilonX : x; +} + +Flow flowBothEnds(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, + double u1, double theta1, double u2, double theta2) { + + const std::complex& v1 = math::ComplexUtils::polar2Complex(u1, theta1); + const std::complex& v2 = math::ComplexUtils::polar2Complex(u2, theta2); + + return flowBothEnds(y11, y12, y21, y22, v1, v2); +} + +Flow flowBothEnds(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, + const std::complex& v1, const std::complex& v2) { + + Flow flow; + const std::complex& ift = y12 * v2 + y11 * v1; + flow.fromTo = std::conj(ift) * v1; + + const std::complex& itf = y21 * v1 + y22 * v2; + flow.toFrom = std::conj(itf) * v2; + + return flow; +} + +std::complex flowYshunt(const std::complex& ysh, double u, double theta) { + std::complex v = math::ComplexUtils::polar2Complex(u, theta); + + return std::conj(ysh) * std::conj(v) *v; +} + +double getPhaseAngleClockDegrees(int phaseAngleClock) { + double phaseAngleClockDegree = 0.0; + phaseAngleClockDegree += phaseAngleClock * 30.0; + phaseAngleClockDegree = std::remainder(phaseAngleClockDegree, 360.0); + if (phaseAngleClockDegree > 180.0) { + phaseAngleClockDegree -= 360.0; + } + return phaseAngleClockDegree; +} + +BranchAdmittanceMatrix kronChain(const BranchAdmittanceMatrix& firstAdm, const Branch::Side& firstChainNodeSide, + const BranchAdmittanceMatrix& secondAdm, const Branch::Side& secondChainNodeSide) { + BranchAdmittanceMatrix admittance; + + std::complex yFirst11; + std::complex yFirst1C; + std::complex yFirstC1; + std::complex yFirstCC; + if (firstChainNodeSide == Branch::Side::TWO) { + yFirst11 = firstAdm.y11; + yFirst1C = firstAdm.y12; + yFirstC1 = firstAdm.y21; + yFirstCC = firstAdm.y22; + } else { + yFirst11 = firstAdm.y22; + yFirst1C = firstAdm.y21; + yFirstC1 = firstAdm.y12; + yFirstCC = firstAdm.y11; + } + + std::complex ySecond22; + std::complex ySecond2C; + std::complex ySecondC2; + std::complex ySecondCC; + if (secondChainNodeSide == Branch::Side::TWO) { + ySecond22 = secondAdm.y11; + ySecond2C = secondAdm.y12; + ySecondC2 = secondAdm.y21; + ySecondCC = secondAdm.y22; + } else { + ySecond22 = secondAdm.y22; + ySecond2C = secondAdm.y21; + ySecondC2 = secondAdm.y12; + ySecondCC = secondAdm.y11; + } + + admittance.y11 = yFirst11 - (yFirst1C * yFirstC1 / (yFirstCC + ySecondCC)); + admittance.y12 = yFirst1C * ySecondC2 / -(yFirstCC + ySecondCC); + admittance.y21 = ySecond2C * yFirstC1 / -(yFirstCC + ySecondCC); + admittance.y22 = ySecond22 - (ySecond2C * ySecondC2 / (yFirstCC + ySecondCC)); + + return admittance; +} + +std::complex kronAntenna(const std::complex& y11, const std::complex& y12, const std::complex& y21, const std::complex& y22, bool isOpenFrom) { + std::complex ysh(0.0, 0.0); + + if (isOpenFrom) { + if (y11 != std::complex(0.0, 0.0)) { + ysh = y22 - (y21 * y12 / y11); + } + } else { + if (y22 != std::complex(0.0, 0.0)) { + ysh = y11 - (y12 * y21 / y22); + } + } + return ysh; +} + } // namespace LinkData } // namespace iidm diff --git a/src/iidm/util/TwtData.cpp b/src/iidm/util/TwtData.cpp new file mode 100644 index 00000000..23340eaf --- /dev/null +++ b/src/iidm/util/TwtData.cpp @@ -0,0 +1,610 @@ +/** + * Copyright (c) 2022, RTE (http://www.rte-france.com) + * 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 http://mozilla.org/MPL/2.0/. + */ + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace powsybl { + +namespace iidm { + +const std::string UNEXPECTED_SIDE = "Unexpected side"; + +TwtData::TwtData(const ThreeWindingsTransformer& twt, double epsilonX, bool applyReactanceCorrection) : + TwtData(twt, 0, 0, epsilonX, applyReactanceCorrection, false) { +} + +TwtData::TwtData(const ThreeWindingsTransformer& twt, double epsilonX, bool applyReactanceCorrection, bool twtSplitShuntAdmittance) : + TwtData(twt, 0, 0, epsilonX, applyReactanceCorrection, twtSplitShuntAdmittance) { +} + +TwtData::TwtData(const ThreeWindingsTransformer& twt, int phaseAngleClock2, int phaseAngleClock3, double epsilonX, + bool applyReactanceCorrection, bool twtSplitShuntAdmittance) : + m_id(twt.getId()), + m_p1(twt.getLeg1().getTerminal().getP()), + m_q1(twt.getLeg1().getTerminal().getQ()), + m_p2(twt.getLeg2().getTerminal().getP()), + m_q2(twt.getLeg2().getTerminal().getQ()), + m_p3(twt.getLeg3().getTerminal().getP()), + m_q3(twt.getLeg3().getTerminal().getQ()), + m_u1(getV(twt.getLeg1())), + m_theta1(getTheta(twt.getLeg1())), + m_u2(getV(twt.getLeg2())), + m_theta2(getTheta(twt.getLeg2())), + m_u3(getV(twt.getLeg3())), + m_theta3(getTheta(twt.getLeg3())), + m_r1(getR(twt.getLeg1())), + m_x1(LinkData::getFixedX(getX(twt.getLeg1()), epsilonX, applyReactanceCorrection)), + m_r2(getR(twt.getLeg2())), + m_x2(LinkData::getFixedX(getX(twt.getLeg2()), epsilonX, applyReactanceCorrection)), + m_r3(getR(twt.getLeg3())), + m_x3(LinkData::getFixedX(getX(twt.getLeg3()), epsilonX, applyReactanceCorrection)), + m_g11(getG1(twt.getLeg1(), twtSplitShuntAdmittance)), + m_b11(getB1(twt.getLeg1(), twtSplitShuntAdmittance)), + m_g12(getG2(twt.getLeg1(), twtSplitShuntAdmittance)), + m_b12(getB2(twt.getLeg1(), twtSplitShuntAdmittance)), + m_g21(getG1(twt.getLeg2(), twtSplitShuntAdmittance)), + m_b21(getB1(twt.getLeg2(), twtSplitShuntAdmittance)), + m_g22(getG2(twt.getLeg2(), twtSplitShuntAdmittance)), + m_b22(getB2(twt.getLeg2(), twtSplitShuntAdmittance)), + m_g31(getG1(twt.getLeg3(), twtSplitShuntAdmittance)), + m_b31(getB1(twt.getLeg3(), twtSplitShuntAdmittance)), + m_g32(getG2(twt.getLeg3(), twtSplitShuntAdmittance)), + m_b32(getB2(twt.getLeg3(), twtSplitShuntAdmittance)), + m_alpha1(alpha(twt.getLeg1())), + m_alpha2(alpha(twt.getLeg2())), + m_alpha3(alpha(twt.getLeg3())), + m_ratedU1(twt.getLeg1().getRatedU()), + m_ratedU2(twt.getLeg2().getRatedU()), + m_ratedU3(twt.getLeg3().getRatedU()), + m_connected1(twt.getLeg1().getTerminal().isConnected()), + m_connected2(twt.getLeg2().getTerminal().isConnected()), + m_connected3(twt.getLeg3().getTerminal().isConnected()), + m_mainComponent1(isMainComponent(twt.getLeg1())), + m_mainComponent2(isMainComponent(twt.getLeg2())), + m_mainComponent3(isMainComponent(twt.getLeg3())), + m_phaseAngleClock2(phaseAngleClock2), + m_phaseAngleClock3(phaseAngleClock3), + m_ratedU0(twt.getRatedU0()) { + m_rho1 = rho(twt.getLeg1(), m_ratedU0); + m_rho2 = rho(twt.getLeg2(), m_ratedU0); + m_rho3 = rho(twt.getLeg3(), m_ratedU0); + + double rhof = 1.0; + double alphaf = 0.0; + + double angle1 = -m_alpha1; + double angle2 = -m_alpha2 - stdcxx::toRadians * LinkData::getPhaseAngleClockDegrees(phaseAngleClock2); + double angle3 = -m_alpha3 - stdcxx::toRadians * LinkData::getPhaseAngleClockDegrees(phaseAngleClock3); + double anglef = -alphaf; + + LinkData::BranchAdmittanceMatrix branchAdmittanceLeg1 = LinkData::calculateBranchAdmittance(m_r1, m_x1, 1 / m_rho1, angle1, 1 / rhof, anglef, std::complex(m_g11, m_b11), std::complex(m_g12, m_b12)); + + LinkData::BranchAdmittanceMatrix branchAdmittanceLeg2 = LinkData::calculateBranchAdmittance(m_r2, m_x2, 1 / m_rho2, angle2, 1 / rhof, anglef, std::complex(m_g21, m_b21), std::complex(m_g22, m_b22)); + + LinkData::BranchAdmittanceMatrix branchAdmittanceLeg3 = LinkData::calculateBranchAdmittance(m_r3, m_x3, 1 / m_rho3, angle3, 1 / rhof, anglef, std::complex(m_g31, m_b31), std::complex(m_g32, m_b32)); + + // Assume the ratedU at the star bus is equal to ratedU of Leg1 + + if (m_connected1 && m_connected2 && m_connected3 && valid(m_u1, m_theta1) && valid(m_u2, m_theta2) && valid(m_u3, m_theta3)) { + + calculateThreeConnectedLegsFlowAndStarBusVoltage(m_u1, m_theta1, m_u2, m_theta2, m_u3, m_theta3, branchAdmittanceLeg1, + branchAdmittanceLeg2, branchAdmittanceLeg3); + } else if (m_connected1 && m_connected2 && valid(m_u1, m_theta1) && valid(m_u2, m_theta2)) { + + LinkData::Flow flow = calculateTwoConnectedLegsFlow(m_u1, m_theta1, m_u2, m_theta2, + branchAdmittanceLeg1, branchAdmittanceLeg2, branchAdmittanceLeg3); + m_computedP1 = std::real(flow.fromTo); + m_computedQ1 = std::imag(flow.fromTo); + m_computedP2 = std::real(flow.toFrom); + m_computedQ2 = std::imag(flow.toFrom); + m_computedP3 = 0.0; + m_computedQ3 = 0.0; + + std::complex v0 = calculateTwoConnectedLegsStarBusVoltage(m_u1, m_theta1, m_u2, m_theta2, + branchAdmittanceLeg1, branchAdmittanceLeg2, branchAdmittanceLeg3); + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else if (m_connected1 && m_connected3 && valid(m_u1, m_theta1) && valid(m_u3, m_theta3)) { + + LinkData::Flow flow = calculateTwoConnectedLegsFlow(m_u1, m_theta1, m_u3, m_theta3, + branchAdmittanceLeg1, branchAdmittanceLeg3, branchAdmittanceLeg2); + m_computedP1 = std::real(flow.fromTo); + m_computedQ1 = std::imag(flow.fromTo); + m_computedP2 = 0.0; + m_computedQ2 = 0.0; + m_computedP3 = std::real(flow.toFrom); + m_computedQ3 = std::imag(flow.toFrom); + + std::complex v0 = calculateTwoConnectedLegsStarBusVoltage(m_u1, m_theta1, m_u3, m_theta3, + branchAdmittanceLeg1, branchAdmittanceLeg3, branchAdmittanceLeg2); + + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else if (m_connected2 && m_connected3 && valid(m_u2, m_theta2) && valid(m_u3, m_theta3)) { + + LinkData::Flow flow = calculateTwoConnectedLegsFlow(m_u2, m_theta2, m_u3, m_theta3, + branchAdmittanceLeg2, branchAdmittanceLeg3, branchAdmittanceLeg1); + m_computedP1 = 0.0; + m_computedQ1 = 0.0; + m_computedP2 = std::real(flow.fromTo); + m_computedQ2 = std::imag(flow.fromTo); + m_computedP3 = std::real(flow.toFrom); + m_computedQ3 = std::imag(flow.toFrom); + + std::complex v0 = calculateTwoConnectedLegsStarBusVoltage(m_u2, m_theta2, m_u3, m_theta3, + branchAdmittanceLeg2, branchAdmittanceLeg3, branchAdmittanceLeg1); + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else if (m_connected1 && valid(m_u1, m_theta1)) { + + std::complex flow = calculateOneConnectedLegFlow(m_u1, m_theta1, branchAdmittanceLeg1, + branchAdmittanceLeg2, branchAdmittanceLeg3); + m_computedP1 = std::real(flow); + m_computedQ1 = std::imag(flow); + m_computedP2 = 0.0; + m_computedQ2 = 0.0; + m_computedP3 = 0.0; + m_computedQ3 = 0.0; + + std::complex v0 = calculateOneConnectedLegStarBusVoltage(m_u1, m_theta1, branchAdmittanceLeg1, + branchAdmittanceLeg2, branchAdmittanceLeg3); + + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else if (m_connected2 && valid(m_u2, m_theta2)) { + + std::complex flow = calculateOneConnectedLegFlow(m_u2, m_theta2, branchAdmittanceLeg2, + branchAdmittanceLeg1, branchAdmittanceLeg3); + + m_computedP1 = 0.0; + m_computedQ1 = 0.0; + m_computedP2 = std::real(flow); + m_computedQ2 = std::imag(flow); + m_computedP3 = 0.0; + m_computedQ3 = 0.0; + + std::complex v0 = calculateOneConnectedLegStarBusVoltage(m_u2, m_theta2, branchAdmittanceLeg2, + branchAdmittanceLeg1, branchAdmittanceLeg3); + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else if (m_connected3 && valid(m_u3, m_theta3)) { + + std::complex flow = calculateOneConnectedLegFlow(m_u3, m_theta3, branchAdmittanceLeg3, + branchAdmittanceLeg1, branchAdmittanceLeg2); + + m_computedP1 = 0.0; + m_computedQ1 = 0.0; + m_computedP2 = 0.0; + m_computedQ2 = 0.0; + m_computedP3 = std::real(flow); + m_computedQ3 = std::imag(flow); + + std::complex v0 = calculateOneConnectedLegStarBusVoltage(m_u3, m_theta3, branchAdmittanceLeg3, + branchAdmittanceLeg1, branchAdmittanceLeg2); + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); + } else { + + m_computedP1 = stdcxx::nan(); + m_computedQ1 = stdcxx::nan(); + m_computedP2 = stdcxx::nan(); + m_computedQ2 = stdcxx::nan(); + m_computedP3 = stdcxx::nan(); + m_computedQ3 = stdcxx::nan(); + + m_starU = stdcxx::nan(); + m_starTheta = stdcxx::nan(); + } +} + +double TwtData::alpha(const ThreeWindingsTransformer::Leg& leg) { + return leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getAlpha() : 0.0; +} + +std::complex TwtData::calculateOneConnectedLegFlow(double u, double theta, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg) const { + + std::complex ysh = calculateOneConnectedLegShunt(admittanceMatrixLeg, + admittanceMatrixFirstOpenLeg, admittanceMatrixSecondOpenLeg); + + return LinkData::flowYshunt(ysh, u, theta); +} + +std::complex TwtData::calculateOneConnectedLegShunt(const LinkData::BranchAdmittanceMatrix& closeLeg, + const LinkData::BranchAdmittanceMatrix& firstOpenLeg, + const LinkData::BranchAdmittanceMatrix& secondOpenLeg) const { + std::complex ysh1 = LinkData::kronAntenna(firstOpenLeg.y11, firstOpenLeg.y12, firstOpenLeg.y21, firstOpenLeg.y22, true); + std::complex ysh2 = LinkData::kronAntenna(secondOpenLeg.y11, secondOpenLeg.y12, secondOpenLeg.y21, secondOpenLeg.y22, true); + std::complex y22 = closeLeg.y22 + ysh1 + ysh2; + + return LinkData::kronAntenna(closeLeg.y11, closeLeg.y12, closeLeg.y21, y22, false); +} + +std::complex TwtData::calculateOneConnectedLegStarBusVoltage(double u, double theta, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg) const { + + std::complex v = math::ComplexUtils::polar2Complex(u, theta); + + std::complex ysh1O = LinkData::kronAntenna(admittanceMatrixFirstOpenLeg.y11, admittanceMatrixFirstOpenLeg.y12, + admittanceMatrixFirstOpenLeg.y21, admittanceMatrixFirstOpenLeg.y22, true); + std::complex ysh2O = LinkData::kronAntenna(admittanceMatrixSecondOpenLeg.y11, admittanceMatrixSecondOpenLeg.y12, + admittanceMatrixSecondOpenLeg.y21, admittanceMatrixSecondOpenLeg.y22, true); + + return (-admittanceMatrixLeg.y21 * v) / (admittanceMatrixLeg.y22 + ysh1O + ysh2O); +} + +void TwtData::calculateThreeConnectedLegsFlowAndStarBusVoltage(double u1, double theta1, double u2, double theta2, + double u3, double theta3, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg1, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg2, + const LinkData::BranchAdmittanceMatrix& branchAdmittanceLeg3) { + std::complex v1 = math::ComplexUtils::polar2Complex(u1, theta1); + std::complex v2 = math::ComplexUtils::polar2Complex(u2, theta2); + std::complex v3 = math::ComplexUtils::polar2Complex(u3, theta3); + + std::complex v0 = -(branchAdmittanceLeg1.y21 * v1 + branchAdmittanceLeg2.y21 * v2 + branchAdmittanceLeg3.y21 * v3) + / (branchAdmittanceLeg1.y22 + branchAdmittanceLeg2.y22 + branchAdmittanceLeg3.y22); + + LinkData::Flow flowLeg1 = LinkData::flowBothEnds(branchAdmittanceLeg1.y11, branchAdmittanceLeg1.y12, + branchAdmittanceLeg1.y21, branchAdmittanceLeg1.y22, v1, v0); + + LinkData::Flow flowLeg2 = LinkData::flowBothEnds(branchAdmittanceLeg2.y11, branchAdmittanceLeg2.y12, + branchAdmittanceLeg2.y21, branchAdmittanceLeg2.y22, v2, v0); + + LinkData::Flow flowLeg3 = LinkData::flowBothEnds(branchAdmittanceLeg3.y11, branchAdmittanceLeg3.y12, + branchAdmittanceLeg3.y21, branchAdmittanceLeg3.y22, v3, v0); + + m_computedP1 = std::real(flowLeg1.fromTo); + m_computedQ1 = std::imag(flowLeg1.fromTo); + m_computedP2 = std::real(flowLeg2.fromTo); + m_computedQ2 = std::imag(flowLeg2.fromTo); + m_computedP3 = std::real(flowLeg3.fromTo); + m_computedQ3 = std::imag(flowLeg3.fromTo); + + m_starU = std::abs(v0); + m_starTheta = std::arg(v0); +} + +LinkData::BranchAdmittanceMatrix TwtData::calculateTwoConnectedLegsAdmittance(const LinkData::BranchAdmittanceMatrix& firstCloseLeg, + const LinkData::BranchAdmittanceMatrix& secondCloseLeg, + const LinkData::BranchAdmittanceMatrix& openLeg) const { + + std::complex ysh = LinkData::kronAntenna(openLeg.y11, openLeg.y12, openLeg.y21, openLeg.y22, true); + LinkData::BranchAdmittanceMatrix secondCloseLegMod(secondCloseLeg.y11, secondCloseLeg.y12, secondCloseLeg.y21, secondCloseLeg.y22 + ysh); + return LinkData::kronChain(firstCloseLeg, Branch::Side::TWO, secondCloseLegMod, Branch::Side::TWO); +} + +LinkData::Flow TwtData::calculateTwoConnectedLegsFlow(double u1, double theta1, double u2, double theta2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg) const { + + std::complex v1 = math::ComplexUtils::polar2Complex(u1, theta1); + std::complex v2 = math::ComplexUtils::polar2Complex(u2, theta2); + + LinkData::BranchAdmittanceMatrix admittance = calculateTwoConnectedLegsAdmittance(admittanceMatrixLeg1, + admittanceMatrixLeg2, admittanceMatrixOpenLeg); + + return LinkData::flowBothEnds(admittance.y11, admittance.y12, admittance.y21, admittance.y22, v1, v2); +} + +std::complex TwtData::calculateTwoConnectedLegsStarBusVoltage(double u1, double theta1, double u2, double theta2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2, + const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg) const { + + std::complex v1 = math::ComplexUtils::polar2Complex(u1, theta1); + std::complex v2 = math::ComplexUtils::polar2Complex(u2, theta2); + + std::complex yshO = LinkData::kronAntenna(admittanceMatrixOpenLeg.y11, admittanceMatrixOpenLeg.y12, admittanceMatrixOpenLeg.y21, admittanceMatrixOpenLeg.y22, true); + return -(admittanceMatrixLeg1.y21 * v1 + (admittanceMatrixLeg2.y21 * v2)) + / (admittanceMatrixLeg1.y22 + admittanceMatrixLeg2.y22 + yshO); +} + +double TwtData::getB1(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_b11; + case ThreeWindingsTransformer::Side::TWO: + return m_b21; + case ThreeWindingsTransformer::Side::THREE: + return m_b31; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getB1(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance) { + return getValue(twtSplitShuntAdmittance ? leg.getB() / 2 : leg.getB(), + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getB() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getB() : 0.0); +} + +double TwtData::getB2(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_b12; + case ThreeWindingsTransformer::Side::TWO: + return m_b22; + case ThreeWindingsTransformer::Side::THREE: + return m_b32; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getB2(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance) { + return getValue(twtSplitShuntAdmittance ? leg.getB() / 2 : 0.0, + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getB() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getB() : 0.0); +} + +double TwtData::getComputedP(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_computedP1; + case ThreeWindingsTransformer::Side::TWO: + return m_computedP2; + case ThreeWindingsTransformer::Side::THREE: + return m_computedP3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getComputedQ(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_computedQ1; + case ThreeWindingsTransformer::Side::TWO: + return m_computedQ2; + case ThreeWindingsTransformer::Side::THREE: + return m_computedQ3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getG1(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_g11; + case ThreeWindingsTransformer::Side::TWO: + return m_g21; + case ThreeWindingsTransformer::Side::THREE: + return m_g31; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getG1(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance) { + return getValue(twtSplitShuntAdmittance ? leg.getG() / 2 : leg.getG(), + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getG() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getG() : 0.0); +} + +double TwtData::getG2(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_g12; + case ThreeWindingsTransformer::Side::TWO: + return m_g22; + case ThreeWindingsTransformer::Side::THREE: + return m_g32; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getG2(const ThreeWindingsTransformer::Leg& leg, bool twtSplitShuntAdmittance) { + return getValue(twtSplitShuntAdmittance ? leg.getG() / 2 : 0.0, + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getG() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getG() : 0.0); +} + +const std::string& TwtData::getId() const { + return m_id; +} + +double TwtData::getP(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_p1; + case ThreeWindingsTransformer::Side::TWO: + return m_p2; + case ThreeWindingsTransformer::Side::THREE: + return m_p3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +int TwtData::getPhaseAngleClock2() const { + return m_phaseAngleClock2; +} + +int TwtData::getPhaseAngleClock3() const { + return m_phaseAngleClock3; +} + +double TwtData::getQ(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_q1; + case ThreeWindingsTransformer::Side::TWO: + return m_q2; + case ThreeWindingsTransformer::Side::THREE: + return m_q3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getR(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_r1; + case ThreeWindingsTransformer::Side::TWO: + return m_r2; + case ThreeWindingsTransformer::Side::THREE: + return m_r3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getR(const ThreeWindingsTransformer::Leg& leg) { + return getValue(leg.getR(), + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getR() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getR() : 0.0); +} + +double TwtData::getRatedU(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_ratedU1; + case ThreeWindingsTransformer::Side::TWO: + return m_ratedU2; + case ThreeWindingsTransformer::Side::THREE: + return m_ratedU3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getRatedU0() const { + return m_ratedU0; +} + +double TwtData::getStarTheta() const { + return m_starTheta; +} + +double TwtData::getStarU() const { + return m_starU; +} + +double TwtData::getTheta(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_theta1; + case ThreeWindingsTransformer::Side::TWO: + return m_theta2; + case ThreeWindingsTransformer::Side::THREE: + return m_theta3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getTheta(const ThreeWindingsTransformer::Leg& leg) { + return leg.getTerminal().isConnected() ? stdcxx::toRadians * leg.getTerminal().getBusView().getBus().get().getAngle() : stdcxx::nan(); +} + +double TwtData::getU(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_u1; + case ThreeWindingsTransformer::Side::TWO: + return m_u2; + case ThreeWindingsTransformer::Side::THREE: + return m_u3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getV(const ThreeWindingsTransformer::Leg& leg) { + return leg.getTerminal().isConnected() ? leg.getTerminal().getBusView().getBus().get().getV() : stdcxx::nan(); +} + +double TwtData::getValue(double initialValue, double rtcStepValue, double ptcStepValue) { + return initialValue * (1 + rtcStepValue / 100) * (1 + ptcStepValue / 100); +} + +double TwtData::getX(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_x1; + case ThreeWindingsTransformer::Side::TWO: + return m_x2; + case ThreeWindingsTransformer::Side::THREE: + return m_x3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +double TwtData::getX(const ThreeWindingsTransformer::Leg& leg) { + return getValue(leg.getX(), + leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getX() : 0.0, + leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getX() : 0.0); +} + +bool TwtData::isConnected(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_connected1; + case ThreeWindingsTransformer::Side::TWO: + return m_connected2; + case ThreeWindingsTransformer::Side::THREE: + return m_connected3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +bool TwtData::isMainComponent(const ThreeWindingsTransformer::Side& side) const { + switch (side) { + case ThreeWindingsTransformer::Side::ONE: + return m_mainComponent1; + case ThreeWindingsTransformer::Side::TWO: + return m_mainComponent2; + case ThreeWindingsTransformer::Side::THREE: + return m_mainComponent3; + default: + throw AssertionError(UNEXPECTED_SIDE + ": " + Enum::toString(side)); + } +} + +bool TwtData::isMainComponent(const ThreeWindingsTransformer::Leg& leg) { + const auto& bus = leg.getTerminal().getBusView().getBus(); + const auto& connectableBus = leg.getTerminal().getBusView().getConnectableBus(); + bool connectableMainComponent = connectableBus && connectableBus.get().isInMainConnectedComponent(); + return bus ? bus.get().isInMainConnectedComponent() : connectableMainComponent; +} + +double TwtData::rho(const ThreeWindingsTransformer::Leg& leg, double ratedU0) { + double rho = ratedU0 / leg.getRatedU(); + rho *= leg.getOptionalRatioTapChanger() ? leg.getRatioTapChanger().getCurrentStep().getRho() : 1.0; + rho *= leg.getOptionalPhaseTapChanger() ? leg.getPhaseTapChanger().getCurrentStep().getRho() : 1.0; + return rho; +} + +bool TwtData::valid(double voltage, double theta) { + if (std::isnan(voltage) || voltage <= 0.0) { + return false; + } + return !std::isnan(theta); +} + +} // namespace iidm + +} // namespace powsybl diff --git a/test/iidm/CMakeLists.txt b/test/iidm/CMakeLists.txt index 150dc3d2..f858e342 100644 --- a/test/iidm/CMakeLists.txt +++ b/test/iidm/CMakeLists.txt @@ -61,6 +61,7 @@ set(UNIT_TEST_SOURCES util/DanglingLineDataTest.cpp util/SVTest.cpp util/TerminalFinderTest.cpp + util/TwtDataTest.cpp ) add_executable(unit-tests-iidm ${UNIT_TEST_SOURCES}) diff --git a/test/iidm/util/DanglingLineDataTest.cpp b/test/iidm/util/DanglingLineDataTest.cpp index e1d62296..a8b86791 100644 --- a/test/iidm/util/DanglingLineDataTest.cpp +++ b/test/iidm/util/DanglingLineDataTest.cpp @@ -71,10 +71,7 @@ Network createDanglingDataNetwork() { bool dlCompareBoundaryBusVoltage(const DanglingLineData& dlData, double boundaryBusU, double boundaryBusAngle) { double tol = 0.00001; - if (std::abs(dlData.getBoundaryBusU() - boundaryBusU) > tol || std::abs(stdcxx::toDegrees * dlData.getBoundaryBusTheta() - boundaryBusAngle) > tol) { - return false; - } - return true; + return !(std::abs(dlData.getBoundaryBusU() - boundaryBusU) > tol || std::abs(stdcxx::toDegrees * dlData.getBoundaryBusTheta() - boundaryBusAngle) > tol); } BOOST_AUTO_TEST_SUITE(DanglingLineDataTestSuite) diff --git a/test/iidm/util/TwtDataTest.cpp b/test/iidm/util/TwtDataTest.cpp new file mode 100644 index 00000000..2f8f9e05 --- /dev/null +++ b/test/iidm/util/TwtDataTest.cpp @@ -0,0 +1,496 @@ +/** + * Copyright (c) 2022, RTE (http://www.rte-france.com) + * 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 http://mozilla.org/MPL/2.0/. + */ + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace powsybl { + +namespace iidm { + +double P1 = 99.218431; +double Q1 = 2.7304328; +double P2 = -216.19819; +double Q2 = -85.368180; +double P3 = 118; +double Q3 = 92.612077; + +double U1 = 412.989001; +double ANGLE1 = -6.78071; +double U2 = 224.315268; +double ANGLE2 = -8.77012; +double U3 = 21.987; +double ANGLE3 = -6.6508; + +double STAR_U = 412.662007016922757; +double STAR_ANGLE = -7.353686938578365; + +double R1 = 0.898462; +double X1 = 17.204128; +double G11 = 0; +double B11 = 2.4375E-6; +double G12 = 0; +double B12 = 0; +double RATED_U1 = 400; +double R2 = 1.070770247933884; +double X2 = 19.6664; +double G21 = 0; +double B21 = 0; +double G22 = 0; +double B22 = 0; +double RATED_U2 = 220; +double R3 = 4.837006802721089; +double X3 = 21.76072562358277; +double G31 = 0; +double B31 = 0; +double G32 = 0; +double B32 = 0; +double RATED_U3 = 21; +int PHASE_ANGLE_CLOCK_2 = 0; +int PHASE_ANGLE_CLOCK_3 = 0; +double RATED_U0 = RATED_U1; + +Network createTwtDataTestNetwork() { + Network network("test", "test"); + Substation& substation = network.newSubstation() + .setId("S1") + .setName("S1_NAME") + .setCountry(Country::FR) + .setTso("TSO") + .add(); + + VoltageLevel& vl1 = substation.newVoltageLevel() + .setId("VL1") + .setName("VL1_NAME") + .setTopologyKind(TopologyKind::BUS_BREAKER) + .setNominalV(380.0) + .setLowVoltageLimit(340.0) + .setHighVoltageLimit(420.0) + .add(); + + Bus& vl1Bus1 = vl1.getBusBreakerView().newBus() + .setId("VL1_BUS1") + .add(); + + vl1.newLoad() + .setId("LOAD1") + .setName("LOAD1_NAME") + .setBus(vl1Bus1.getId()) + .setConnectableBus(vl1Bus1.getId()) + .setLoadType(LoadType::UNDEFINED) + .setP0(50.0) + .setQ0(40.0) + .add(); + + VoltageLevel& vl2 = substation.newVoltageLevel() + .setId("VL2") + .setName("VL2_NAME") + .setTopologyKind(TopologyKind::BUS_BREAKER) + .setNominalV(225.0) + .setLowVoltageLimit(200.0) + .setHighVoltageLimit(260.0) + .add(); + + Bus& vl2Bus1 = vl2.getBusBreakerView().newBus() + .setId("VL2_BUS1") + .add(); + + vl2.newLoad() + .setId("LOAD2") + .setName("LOAD2_NAME") + .setBus(vl2Bus1.getId()) + .setConnectableBus(vl2Bus1.getId()) + .setLoadType(LoadType::UNDEFINED) + .setP0(60.0) + .setQ0(70.0) + .add(); + + VoltageLevel& vl3 = substation.newVoltageLevel() + .setId("VL3") + .setName("VL3_NAME") + .setTopologyKind(TopologyKind::BUS_BREAKER) + .setNominalV(380.0) + .setLowVoltageLimit(340.0) + .setHighVoltageLimit(420.0) + .add(); + + Bus& vl3Bus1 = vl3.getBusBreakerView().newBus() + .setId("VL3_BUS1") + .add(); + + Substation& substation2 = network.newSubstation() + .setId("S2") + .setName("S2_NAME") + .setCountry(Country::FR) + .setTso("TSO") + .add(); + + VoltageLevel& vl4 = substation2.newVoltageLevel() + .setId("VL4") + .setName("VL4_NAME") + .setTopologyKind(TopologyKind::BUS_BREAKER) + .setNominalV(225.0) + .setLowVoltageLimit(200.0) + .setHighVoltageLimit(260.0) + .add(); + + vl4.getBusBreakerView().newBus() + .setId("VL4_BUS1") + .add(); + + substation.newThreeWindingsTransformer() + .setId("3WT_VL1_VL2_VL3") + .setName("3WT_VL1_VL2_VL3_NAME") + .setRatedU0(RATED_U0) + .newLeg1() + .setR(1.3) + .setX(1.4) + .setG(1.6) + .setB(1.7) + .setRatedU(1.1) + .setRatedS(2.2) + .setVoltageLevel(vl1.getId()) + .setBus(vl1Bus1.getId()) + .setConnectableBus(vl1Bus1.getId()) + .add() + .newLeg2() + .setR(2.3) + .setX(2.4) + .setG(0.0) + .setB(0.0) + .setRatedU(2.1) + .setVoltageLevel(vl2.getId()) + .setBus(vl2Bus1.getId()) + .setConnectableBus(vl2Bus1.getId()) + .add() + .newLeg3() + .setR(3.3) + .setX(3.4) + .setG(0.0) + .setB(0.0) + .setRatedU(3.1) + .setVoltageLevel(vl3.getId()) + .setBus(vl3Bus1.getId()) + .setConnectableBus(vl3Bus1.getId()) + .add() + .add(); + + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + + twt.getLeg1().getTerminal().getBusBreakerView().getBus().get().setV(U1); + twt.getLeg1().getTerminal().getBusBreakerView().getBus().get().setAngle(ANGLE1); + twt.getLeg1().getTerminal().setP(P1); + twt.getLeg1().getTerminal().setQ(Q1); + twt.getLeg1().setR(R1); + twt.getLeg1().setX(X1); + twt.getLeg1().setRatedU(RATED_U1); + twt.getLeg1().setB(B11 + B12); + twt.getLeg1().setG(G11 + G12); + + twt.getLeg2().getTerminal().getBusBreakerView().getBus().get().setV(U2); + twt.getLeg2().getTerminal().getBusBreakerView().getBus().get().setAngle(ANGLE2); + twt.getLeg2().getTerminal().setP(P2); + twt.getLeg2().getTerminal().setQ(Q2); + twt.getLeg2().setR(R2); + twt.getLeg2().setX(X2); + twt.getLeg2().setRatedU(RATED_U2); + twt.getLeg2().setB(B21 + B22); + twt.getLeg2().setG(G21 + G22); + + twt.getLeg3().getTerminal().getBusBreakerView().getBus().get().setV(U3); + twt.getLeg3().getTerminal().getBusBreakerView().getBus().get().setAngle(ANGLE3); + twt.getLeg3().getTerminal().setP(P3); + twt.getLeg3().getTerminal().setQ(Q3); + twt.getLeg3().setR(R3); + twt.getLeg3().setX(X3); + twt.getLeg3().setRatedU(RATED_U3); + twt.getLeg3().setB(B31 + B32); + twt.getLeg3().setG(G31 + G32); + + return network; +} + +struct T3xFlow { + double p1 = stdcxx::nan(); + double q1 = stdcxx::nan(); + double p2 = stdcxx::nan(); + double q2 = stdcxx::nan(); + double p3 = stdcxx::nan(); + double q3 = stdcxx::nan(); +}; + +bool sameFlow(const T3xFlow& expected, const T3xFlow& actual) { + double tol = 0.00001; + return !((!std::isnan(expected.p1) && std::isnan(actual.p1)) || + (!std::isnan(expected.q1) && std::isnan(actual.q1)) || + (!std::isnan(expected.p2) && std::isnan(actual.p2)) || + (!std::isnan(expected.q2) && std::isnan(actual.q2)) || + (!std::isnan(expected.p3) && std::isnan(actual.p3)) || + (!std::isnan(expected.q3) && std::isnan(actual.q3)) || + std::abs(expected.p1 - actual.p1) > tol || + std::abs(expected.q1 - actual.q1) > tol || + std::abs(expected.p2 - actual.p2) > tol || + std::abs(expected.q2 - actual.q2) > tol || + std::abs(expected.p3 - actual.p3) > tol || + std::abs(expected.q3 - actual.q3) > tol); +} + +bool t3xCompareFlow(const TwtData& twtData, double p1, double q1, double p2, double q2, double p3, double q3) { + T3xFlow actual; + actual.p1 = twtData.getComputedP(ThreeWindingsTransformer::Side::ONE); + actual.q1 = twtData.getComputedQ(ThreeWindingsTransformer::Side::ONE); + actual.p2 = twtData.getComputedP(ThreeWindingsTransformer::Side::TWO); + actual.q2 = twtData.getComputedQ(ThreeWindingsTransformer::Side::TWO); + actual.p3 = twtData.getComputedP(ThreeWindingsTransformer::Side::THREE); + actual.q3 = twtData.getComputedQ(ThreeWindingsTransformer::Side::THREE); + + T3xFlow expected; + expected.p1 = p1; + expected.q1 = q1; + expected.p2 = p2; + expected.q2 = q2; + expected.p3 = p3; + expected.q3 = q3; + + return sameFlow(expected, actual); +} + +bool t3xCompareStarBusVoltage(const TwtData& twtData, double starU, double starAngle) { + double tol = 0.00001; + return !(std::isnan(twtData.getStarU()) && !std::isnan(starU)) || + (std::isnan(twtData.getStarTheta()) && !std::isnan(starAngle)) || + std::abs(twtData.getStarU() - starU) > tol || + std::abs(stdcxx::toDegrees * twtData.getStarTheta() - starAngle) > tol; +} + +BOOST_AUTO_TEST_SUITE(TwtDataTestSuite) + +BOOST_AUTO_TEST_CASE(test) { + Network network = createTwtDataTestNetwork(); + const ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + TwtData twtData(twt, 0, false, false); + + BOOST_CHECK_CLOSE(P1, twtData.getComputedP(ThreeWindingsTransformer::Side::ONE), 0.3); + BOOST_CHECK_CLOSE(Q1, twtData.getComputedQ(ThreeWindingsTransformer::Side::ONE), 0.7); // 2.7304328 != 2.7471471852083678 + BOOST_CHECK_CLOSE(P2, twtData.getComputedP(ThreeWindingsTransformer::Side::TWO), 0.3); + BOOST_CHECK_CLOSE(Q2, twtData.getComputedQ(ThreeWindingsTransformer::Side::TWO), 0.3); + BOOST_CHECK_CLOSE(P3, twtData.getComputedP(ThreeWindingsTransformer::Side::THREE), 0.3); + BOOST_CHECK_CLOSE(Q3, twtData.getComputedQ(ThreeWindingsTransformer::Side::THREE), 0.3); + + BOOST_CHECK_CLOSE(STAR_U, twtData.getStarU(), 0.0001); + BOOST_CHECK_CLOSE(STAR_ANGLE, stdcxx::toDegrees * twtData.getStarTheta(), 0.0001); +} + +BOOST_AUTO_TEST_CASE(testSplitShuntAdmittance) { + Network network = createTwtDataTestNetwork(); + const ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + TwtData twtData(twt, 0, false, true); + + bool ok = t3xCompareFlow(twtData, 99.231950, 2.876479, -216.194348, -85.558437, 117.981856, 92.439531); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd1End2End3Connected) { + Network network = createTwtDataTestNetwork(); + const ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + TwtData twtData(twt, 0, false); + + bool ok = t3xCompareFlow(twtData, 99.227288294050, 2.747147185208, -216.195866533486, -85.490493190353, 117.988318295633, 92.500849015581); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 412.66200701692287, -7.353686938578365); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd2End3Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg1().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, 0.0, 0.0, -164.099476216398, -81.835885442800, 165.291731946141, 89.787051339157); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 412.29478568401856, -7.700275244269859); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd1End3Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg2().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, -18.723067158829, -59.239225729782, 0.0, 0.0, 18.851212571411, 59.694062940578); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 415.4806896701992, -6.690799426080698); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd1End2Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg3().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, 161.351352526949, 51.327798049323, -161.019856627996, -45.536840365345, 0.0, 0.0); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 410.53566804098494, -7.703116461849692); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd1Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg2().getTerminal().disconnect(); + twt.getLeg3().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, 0.0, -0.415739792683, 0.0, 0.0, 0.0, 0.0); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 412.9890009999999, -6.78071000000000); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd2Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg1().getTerminal().disconnect(); + twt.getLeg3().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, 0.0, 0.0, 0.000001946510, -0.405486077928, 0.0, 0.0); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 407.8654944214268, -8.77026956158324); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd3Connected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg1().getTerminal().disconnect(); + twt.getLeg2().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, 0.0, 0.0, 0.0, 0.0, 0.000005977974, -0.427562118410); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, 418.82221596280823, -6.65147559975559); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(testEnd1End2End3Disconnected) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + twt.getLeg1().getTerminal().disconnect(); + twt.getLeg2().getTerminal().disconnect(); + twt.getLeg3().getTerminal().disconnect(); + + TwtData twtData(twt, 0, false); + bool ok = t3xCompareFlow(twtData, stdcxx::nan(), stdcxx::nan(), stdcxx::nan(), stdcxx::nan(), stdcxx::nan(), stdcxx::nan()); + BOOST_CHECK(ok); + + ok = t3xCompareStarBusVoltage(twtData, stdcxx::nan(), stdcxx::nan()); + BOOST_CHECK(ok); +} + +BOOST_AUTO_TEST_CASE(checkInputs) { + Network network = createTwtDataTestNetwork(); + ThreeWindingsTransformer& twt = network.getThreeWindingsTransformer("3WT_VL1_VL2_VL3"); + + BOOST_CHECK_CLOSE(2.4375e-6, TwtData(twt, 0, false).getB1(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getB1(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getB1(ThreeWindingsTransformer::Side::THREE)); + + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getB2(ThreeWindingsTransformer::Side::ONE)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getB2(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getB2(ThreeWindingsTransformer::Side::THREE)); + + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG1(ThreeWindingsTransformer::Side::ONE)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG1(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG1(ThreeWindingsTransformer::Side::THREE)); + + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG2(ThreeWindingsTransformer::Side::ONE)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG2(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK_EQUAL(0.0, TwtData(twt, 0, false).getG2(ThreeWindingsTransformer::Side::THREE)); + + BOOST_CHECK_EQUAL("3WT_VL1_VL2_VL3", TwtData(twt, 0, false).getId()); + + BOOST_CHECK_EQUAL(0, TwtData(twt, 0, false).getPhaseAngleClock2()); + BOOST_CHECK_EQUAL(0, TwtData(twt, 0, false).getPhaseAngleClock3()); + + BOOST_CHECK_CLOSE(2.7304328, TwtData(twt, 0, false).getQ(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(-85.36818, TwtData(twt, 0, false).getQ(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(92.612077, TwtData(twt, 0, false).getQ(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(0.898462, TwtData(twt, 0, false).getR(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(1.070770247933884, TwtData(twt, 0, false).getR(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(4.837006802721089, TwtData(twt, 0, false).getR(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(400.0, TwtData(twt, 0, false).getRatedU(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(220.0, TwtData(twt, 0, false).getRatedU(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(21.0, TwtData(twt, 0, false).getRatedU(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(400, TwtData(twt, 0, false).getRatedU0(), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(-0.1183457151229047, TwtData(twt, 0, false).getTheta(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(-0.1530674697950051, TwtData(twt, 0, false).getTheta(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(-0.11607835789163888, TwtData(twt, 0, false).getTheta(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(412.989001, TwtData(twt, 0, false).getU(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(224.315268, TwtData(twt, 0, false).getU(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(21.987, TwtData(twt, 0, false).getU(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(17.204128, TwtData(twt, 0, false).getX(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(19.6664, TwtData(twt, 0, false).getX(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(21.76072562358277, TwtData(twt, 0, false).getX(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK_CLOSE(99.218431, TwtData(twt, 0, false).getP(ThreeWindingsTransformer::Side::ONE), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(-216.19819, TwtData(twt, 0, false).getP(ThreeWindingsTransformer::Side::TWO), std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(118.0, TwtData(twt, 0, false).getP(ThreeWindingsTransformer::Side::THREE), std::numeric_limits::epsilon()); + + BOOST_CHECK(TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::ONE)); + BOOST_CHECK(TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::ONE)); + twt.getLeg1().getTerminal().disconnect(); + BOOST_CHECK(!TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::ONE)); + BOOST_CHECK(!TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::ONE)); + + BOOST_CHECK(TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK(TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::TWO)); + twt.getLeg2().getTerminal().disconnect(); + BOOST_CHECK(!TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::TWO)); + BOOST_CHECK(!TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::TWO)); + + BOOST_CHECK(TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::THREE)); + BOOST_CHECK(TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::THREE)); + twt.getLeg3().getTerminal().disconnect(); + BOOST_CHECK(!TwtData(twt, 0, false).isConnected(ThreeWindingsTransformer::Side::THREE)); + BOOST_CHECK(!TwtData(twt, 0, false).isMainComponent(ThreeWindingsTransformer::Side::THREE)); +} + +BOOST_AUTO_TEST_SUITE_END() + +} // namespace iidm + +} // namespace powsybl