Add Twtdata and unit tests

Signed-off-by: Sébastien LAIGRE <[email protected]>
powsybl · May 6, 2022 · 0b44f0a · 0b44f0a
1 parent 1a7835b
commit 0b44f0a
Show file tree

Hide file tree

Showing 8 changed files with 1,461 additions and 4 deletions.
diff --git a/include/powsybl/iidm/util/LinkData.hpp b/include/powsybl/iidm/util/LinkData.hpp
@@ -10,6 +10,8 @@
 
 #include <complex>
 
+#include <powsybl/iidm/Branch.hpp>
+
 namespace powsybl {
 
 namespace iidm {
@@ -25,6 +27,10 @@ struct Flow {
 
 struct BranchAdmittanceMatrix {
 public:
+    BranchAdmittanceMatrix() = default;
+
+    BranchAdmittanceMatrix(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22);
+
     std::complex<double> y11;
 
     std::complex<double> y12;
@@ -37,6 +43,23 @@ struct BranchAdmittanceMatrix {
 BranchAdmittanceMatrix calculateBranchAdmittance(double r, double x, double ratio1, double angle1,
                                                  double ratio2, double angle2, const std::complex<double>& ysh1, const std::complex<double>& ysh2);
 
+double getFixedX(double x, double epsilonX, bool applyReactanceCorrection);
+
+Flow flowBothEnds(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22,
+                  double u1, double theta1, double u2, double theta2);
+
+Flow flowBothEnds(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22,
+                  const std::complex<double>& v1, const std::complex<double>& v2);
+
+std::complex<double> flowYshunt(const std::complex<double>& 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<double> kronAntenna(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22, bool isOpenFrom);
+
 }  // namespace LinkData
 
 }  // namespace iidm

diff --git a/include/powsybl/iidm/util/TwtData.hpp 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 <complex>
+#include <string>
+
+#include <powsybl/iidm/ThreeWindingsTransformer.hpp>
+#include <powsybl/iidm/util/LinkData.hpp>
+
+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<double> calculateOneConnectedLegFlow(double u, double theta,
+                                                      const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg,
+                                                      const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg,
+                                                      const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg);
+
+    std::complex<double> calculateOneConnectedLegShunt(const LinkData::BranchAdmittanceMatrix& closeLeg,
+                                                       const LinkData::BranchAdmittanceMatrix& firstOpenLeg,
+                                                       const LinkData::BranchAdmittanceMatrix& secondOpenLeg);
+
+    std::complex<double> calculateOneConnectedLegStarBusVoltage(double u, double theta,
+                                                                const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg,
+                                                                const LinkData::BranchAdmittanceMatrix& admittanceMatrixFirstOpenLeg,
+                                                                const LinkData::BranchAdmittanceMatrix& admittanceMatrixSecondOpenLeg);
+
+    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);
+
+    LinkData::Flow calculateTwoConnectedLegsFlow(double u1, double theta1, double u2, double theta2,
+                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1,
+                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2,
+                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg);
+
+    std::complex<double> calculateTwoConnectedLegsStarBusVoltage(double u1, double theta1, double u2, double theta2,
+                                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg1,
+                                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixLeg2,
+                                                                 const LinkData::BranchAdmittanceMatrix& admittanceMatrixOpenLeg);
+
+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
@@ -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
@@ -7,12 +7,21 @@
 
 #include <powsybl/iidm/util/LinkData.hpp>
 
+#include <powsybl/math/ComplexUtils.hpp>
+
 namespace powsybl {
 
 namespace iidm {
 
 namespace LinkData {
 
+BranchAdmittanceMatrix::BranchAdmittanceMatrix(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& 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<double>& ysh1, const std::complex<double>& ysh2) {
     std::complex<double> 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<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22,
+                  double u1, double theta1, double u2, double theta2) {
+
+    const std::complex<double>& v1 = math::ComplexUtils::polar2Complex(u1, theta1);
+    const std::complex<double>& v2 = math::ComplexUtils::polar2Complex(u2, theta2);
+
+    return flowBothEnds(y11, y12, y21, y22, v1, v2);
+}
+
+Flow flowBothEnds(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22,
+                  const std::complex<double>& v1, const std::complex<double>& v2) {
+
+    Flow flow;
+    const std::complex<double>& ift = y12 * v2 + y11 * v1;
+    flow.fromTo = std::conj(ift) * v1;
+
+    const std::complex<double>& itf = y21 * v1 + y22 * v2;
+    flow.toFrom = std::conj(itf) * v2;
+
+    return flow;
+}
+
+std::complex<double> flowYshunt(const std::complex<double>& ysh, double u, double theta) {
+    std::complex<double> 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<double> yFirst11;
+    std::complex<double> yFirst1C;
+    std::complex<double> yFirstC1;
+    std::complex<double> 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<double> ySecond22;
+    std::complex<double> ySecond2C;
+    std::complex<double> ySecondC2;
+    std::complex<double> 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<double> kronAntenna(const std::complex<double>& y11, const std::complex<double>& y12, const std::complex<double>& y21, const std::complex<double>& y22, bool isOpenFrom) {
+    std::complex<double> ysh(0.0, 0.0);
+
+    if (isOpenFrom) {
+        if (y11 != std::complex<double>(0.0, 0.0)) {
+            ysh = y22 - (y21 * y12 / y11);
+        }
+    } else {
+        if (y22 != std::complex<double>(0.0, 0.0)) {
+            ysh = y11 - (y12 * y21 / y22);
+        }
+    }
+    return ysh;
+}
+
 }  // namespace LinkData
 
 }  // namespace iidm