Fix three winding transformer calculation

CHANGELOG.md

@@ -20,5 +20,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Location of `vn_simona` test grid (was partially in Berlin and Dortmund)
 - Let `ParticipantAgent` die after failed registration with secondary services (prevents stuck simulation)
 - Fix default resolution of weather source wrapper [#78](https://github.com/ie3-institute/simona/issues/78)
+- Support for three winding transformers  [#63](https://github.com/ie3-institute/simona/issues/63)
+  - Handle incoming slack voltage accordingly
+  - Allow multiple sub grid gates at one node (also allows multiple two winding transformers at one node)
+  - Perform power flow calculation in highest grid, if a three winding transformer is apparent
+  - Write out results
 
 [Unreleased]: https://github.com/ie3-institute/simona/compare/a14a093239f58fca9b2b974712686b33e5e5f939...HEAD

src/main/scala/edu/ie3/simona/agent/grid/GridResultsSupport.scala

@@ -422,9 +422,8 @@ private[grid] trait GridResultsSupport {
     }
   }
 
-  /** NOT IMPLEMENTED YET AS NO TEST DATA IS AVAILABLE! Creates an instance of
-    * [[Transformer3WResult]] based on the provided grid and power flow result
-    * data
+  /** Creates an instance of [[Transformer3WResult]] based on the provided grid
+    * and power flow result data
     *
     * @param trafo3w
     *   the instance of the 3 winding transformer that should be processed
@@ -447,12 +446,16 @@ private[grid] trait GridResultsSupport {
       iNominal: ComparableQuantity[ElectricCurrent],
       timestamp: ZonedDateTime
   ): PartialTransformer3wResult = {
-    val (iMag, iAng) = calcPortCurrent(
-      trafo3w,
-      nodeStateData.voltage,
+    val (_, iComplexPu) = iIJComplexPu(
       internalNodeStateData.voltage,
-      iNominal
+      nodeStateData.voltage,
+      yij(trafo3w),
+      Complex.zero,
+      None
     )
+
+    val (iMag, iAng) = iMagAndAngle(iComplexPu, iNominal)
+
     trafo3w.powerFlowCase match {
       case Transformer3wPowerFlowCase.PowerFlowCaseA =>
         PartialTransformer3wResult.PortA(
@@ -479,47 +482,6 @@ private[grid] trait GridResultsSupport {
     }
   }
 
-  /** Calculate the port current of the transformer
-    *
-    * @param transformer
-    *   The transformer model
-    * @param v1
-    *   Nodal voltage at the port
-    * @param v2
-    *   Nodal voltage at internal node
-    * @param iNominal
-    *   Nominal current
-    * @return
-    *   Magnitude and angle of the current
-    */
-  def calcPortCurrent(
-      transformer: Transformer3wModel,
-      v1: Complex,
-      v2: Complex,
-      iNominal: ComparableQuantity[ElectricCurrent]
-  ): (ComparableQuantity[ElectricCurrent], ComparableQuantity[Angle]) = {
-    val y = yij(transformer)
-    val (de, df) = (v1, v2) match {
-      case (Complex(e1, f1), Complex(e2, f2)) =>
-        (e1 - e2, f1 - f2)
-    }
-    val iReal = (de * y.real - df * y.imag) / (pow(y.real, 2) + pow(y.imag, 2))
-    val iImag = (de * y.imag + df * y.real) / (pow(y.real, 2) + pow(y.imag, 2))
-
-    val iMag = iNominal.multiply(sqrt(pow(iReal, 2) + pow(iImag, 2)))
-    val iAng = atan(iImag / iReal) match {
-      case angle if angle.isNaN =>
-        Quantities
-          .getQuantity(copySign(90d, iImag), PowerSystemUnits.DEGREE_GEOM)
-          .to(StandardUnits.ELECTRIC_CURRENT_ANGLE)
-      case angle =>
-        Quantities
-          .getQuantity(angle, Units.RADIAN)
-          .to(StandardUnits.ELECTRIC_CURRENT_ANGLE)
-    }
-    (iMag, iAng)
-  }
-
   /** Calculate the voltage magnitude and the voltage angle in physical units
     * based on a provided electric current in p.u. and the nominal referenced
     * electric current. The arctangent "only" calculates the angle between the

src/main/scala/edu/ie3/simona/model/grid/Transformer3wModel.scala

@@ -6,11 +6,15 @@
 
 package edu.ie3.simona.model.grid
 
+import breeze.linalg.max
+
 import java.time.ZonedDateTime
 import java.util.UUID
 import breeze.math.Complex
 import breeze.numerics.pow
+import com.typesafe.scalalogging.LazyLogging
 import edu.ie3.datamodel.exceptions.InvalidGridException
+import edu.ie3.datamodel.models.StandardUnits
 import edu.ie3.datamodel.models.input.connector.Transformer3WInput
 import edu.ie3.datamodel.models.input.connector.`type`.Transformer3WTypeInput
 import edu.ie3.simona.exceptions.{
@@ -29,7 +33,7 @@ import edu.ie3.util.scala.OperationInterval
 
 import javax.measure.Quantity
 import javax.measure.quantity.{Dimensionless, ElectricPotential}
-import tech.units.indriya.ComparableQuantity
+import tech.units.indriya.{AbstractUnit, ComparableQuantity}
 import tech.units.indriya.quantity.Quantities
 
 import scala.math.BigDecimal.RoundingMode
@@ -159,7 +163,7 @@ final case class Transformer3wModel(
   }
 }
 
-case object Transformer3wModel {
+case object Transformer3wModel extends LazyLogging {
 
   def apply(
       transformer3wInput: Transformer3WInput,
@@ -319,66 +323,56 @@ case object Transformer3wModel {
     val transformerRefSystem =
       RefSystem(transformerType.getsRatedA, transformerType.getvRatedA)
 
-    /* Extract the equivalent circuit diagram parameters from type, with the perspective of the
-     * transformer */
+    /* Get the physical equivalent circuit diagram parameters from type. They come with reference to the highest
+     * voltage side, therefore, in power flow case B and C, they need to be adapted. */
     val (rTrafo, xTrafo, gTrafo, bTrafo) = powerFlowCase match {
       case PowerFlowCaseA =>
         (
-          transformerRefSystem.rInPu(transformerType.getrScA),
-          transformerRefSystem.xInPu(transformerType.getxScA),
-          transformerRefSystem.gInPu(transformerType.getgM),
-          transformerRefSystem.gInPu(transformerType.getbM)
+          transformerType.getrScA,
+          transformerType.getxScA,
+          transformerType.getgM,
+          transformerType.getbM
         )
       case PowerFlowCaseB =>
+        val nominalRatio = transformerType
+          .getvRatedA()
+          .divide(transformerType.getvRatedB())
+          .asType(classOf[Dimensionless])
+          .to(AbstractUnit.ONE)
+          .getValue
+          .doubleValue()
         (
-          transformerRefSystem.rInPu(transformerType.getrScB),
-          transformerRefSystem.xInPu(transformerType.getxScB),
-          Quantities.getQuantity(0d, PU),
-          Quantities.getQuantity(0d, PU)
+          transformerType.getrScB.divide(pow(nominalRatio, 2)),
+          transformerType.getxScB.divide(pow(nominalRatio, 2)),
+          Quantities.getQuantity(0d, StandardUnits.ADMITTANCE),
+          Quantities.getQuantity(0d, StandardUnits.ADMITTANCE)
         )
       case PowerFlowCaseC =>
+        val nominalRatio = transformerType
+          .getvRatedA()
+          .divide(transformerType.getvRatedC())
+          .asType(classOf[Dimensionless])
+          .to(AbstractUnit.ONE)
+          .getValue
+          .doubleValue()
         (
-          transformerRefSystem.rInPu(transformerType.getrScC),
-          transformerRefSystem.xInPu(transformerType.getxScC),
-          Quantities.getQuantity(0d, PU),
-          Quantities.getQuantity(0d, PU)
+          transformerType.getrScC.divide(pow(nominalRatio, 2)),
+          transformerType.getxScC.divide(pow(nominalRatio, 2)),
+          Quantities.getQuantity(0d, StandardUnits.ADMITTANCE),
+          Quantities.getQuantity(0d, StandardUnits.ADMITTANCE)
         )
     }
 
-    /* Translate the single parameters to the grid's reference system */
+    /* Translate the single parameters to dimensionless units based on the grid's reference system */
     (
       /* r */
-      Quantities.getQuantity(
-        RefSystem
-          .transferImpedance(rTrafo, transformerRefSystem, refSystem)
-          .getValue
-          .doubleValue(),
-        PU
-      ),
+      refSystem.rInPu(rTrafo),
       /* x */
-      Quantities.getQuantity(
-        RefSystem
-          .transferImpedance(xTrafo, transformerRefSystem, refSystem)
-          .getValue
-          .doubleValue(),
-        PU
-      ),
+      refSystem.xInPu(xTrafo),
       /* g */
-      Quantities.getQuantity(
-        RefSystem
-          .transferAdmittance(gTrafo, transformerRefSystem, refSystem)
-          .getValue
-          .doubleValue(),
-        PU
-      ),
+      refSystem.gInPu(gTrafo),
       /* b */
-      Quantities.getQuantity(
-        RefSystem
-          .transferAdmittance(bTrafo, transformerRefSystem, refSystem)
-          .getValue
-          .doubleValue(),
-        PU
-      )
+      refSystem.bInPu(bTrafo)
     )
   }
 
@@ -461,15 +455,27 @@ case object Transformer3wModel {
       )
 
     // check if nominal power of winding A is able to supply winding B and C
-    val nomSwindingA = trafo3wType.getsRatedA.getValue.doubleValue
-    val nomSwindingB = trafo3wType.getsRatedB.getValue.doubleValue
-    val nomSwindingC = trafo3wType.getsRatedC.getValue.doubleValue
-
-    if (nomSwindingA < (nomSwindingB + nomSwindingC))
-      throw new InvalidParameterException(
-        s"The winding A of transformer type has a lower rating as both windings B and C together! " +
-          s"($nomSwindingA < $nomSwindingB + $nomSwindingC)"
-      )
+    val sNomA = trafo3wType.getsRatedA.getValue.doubleValue
+    val sNomB = trafo3wType.getsRatedB.getValue.doubleValue
+    val sNomC = trafo3wType.getsRatedC.getValue.doubleValue
+
+    if (sNomA < (sNomB + sNomC)) {
+      val maxSNomUnderlying = max(sNomB, sNomC)
+      if (sNomA < maxSNomUnderlying)
+        throw new InvalidParameterException(
+          s"The winding A of transformer type has a lower rating ($sNomA) as winding B ($sNomB) or C ($sNomC)!"
+        )
+      else
+        logger.warn(
+          "The port A of three winding transformer type {} ({}) has lower power rating ({}) than both underlying ports together ({} + {} = {})!",
+          trafo3wType.getUuid,
+          trafo3wType.getId,
+          trafo3wType.getsRatedA(),
+          trafo3wType.getsRatedB(),
+          trafo3wType.getsRatedC(),
+          trafo3wType.getsRatedB().add(trafo3wType.getsRatedC())
+        )
+    }
   }
 
   /** Calculates the current, tap dependent voltage ratio between the high
@@ -521,7 +527,7 @@ case object Transformer3wModel {
         val bij = transformer3wModel.bij().getValue.doubleValue()
         val gii = transformer3wModel.g0().getValue.doubleValue()
         val bii = transformer3wModel.b0().getValue.doubleValue()
-        amount * ((1 - transformer3wModel.tapRatio) * Complex(
+        amount * ((1 - 1 / transformer3wModel.tapRatio) * Complex(
           gij,
           bij
         ) + Complex(
@@ -550,7 +556,7 @@ case object Transformer3wModel {
     val bij = transformer3wModel.bij().getValue.doubleValue()
     transformer3wModel.powerFlowCase match {
       case PowerFlowCaseA =>
-        amount * pow(transformer3wModel.tapRatio, 2) * Complex(gij, bij)
+        amount * Complex(gij, bij) / transformer3wModel.tapRatio
       case _ => amount * Complex(gij, bij)
     }
   }

src/test/scala/edu/ie3/simona/agent/grid/GridResultsSupportSpec.scala

@@ -467,44 +467,6 @@ class GridResultsSupportSpec
         StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
       )
 
-      "calculating the port currents correctly" in {
-        forAll(
-          Table(
-            ("transformer", "expectedMagnitude", "expectedAngle"),
-            (transformerA, 0.61188, -45.00000),
-            (transformerB, 0.70710, -45.00000),
-            (transformerC, 0.70710, -45.00000)
-          )
-        ) {
-          (
-              transformer: Transformer3wModel,
-              expectedMagnitude: Double,
-              expectedAngle: Double
-          ) =>
-            {
-              val v1 = Complex(1.0, 0.0)
-              val v2 = Complex(0.97, -0.01)
-              calcPortCurrent(transformer, v1, v2, iNominal) match {
-                case (actualMagnitude, actualAngle) =>
-                  actualMagnitude should equalWithTolerance(
-                    Quantities.getQuantity(
-                      expectedMagnitude,
-                      StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
-                    ),
-                    1e-4
-                  )
-                  actualAngle should equalWithTolerance(
-                    Quantities.getQuantity(
-                      expectedAngle,
-                      StandardUnits.ELECTRIC_CURRENT_ANGLE
-                    ),
-                    1e-4
-                  )
-              }
-            }
-        }
-      }
-
       val timeStamp =
         TimeUtil.withDefaults.toZonedDateTime("2021-06-10 14:45:00")
       "assemble correct result for transformer at node A" in {
@@ -531,7 +493,7 @@ class GridResultsSupportSpec
             tapPos shouldBe transformerA.currentTapPos
             currentMagnitude should equalWithTolerance(
               Quantities.getQuantity(
-                0.6118825,
+                13.15547500,
                 StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
               ),
               1e-6
@@ -567,14 +529,14 @@ class GridResultsSupportSpec
             input shouldBe transformerB.uuid
             currentMagnitude should equalWithTolerance(
               Quantities.getQuantity(
-                0.7071067,
+                14.14213562,
                 StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
               ),
               1e-6
             )
             currentAngle should equalWithTolerance(
               Quantities
-                .getQuantity(-45.000000, StandardUnits.ELECTRIC_CURRENT_ANGLE),
+                .getQuantity(135.000000, StandardUnits.ELECTRIC_CURRENT_ANGLE),
               1e-6
             )
           case wrong => fail(s"Got wrong result: '$wrong'")
@@ -603,14 +565,14 @@ class GridResultsSupportSpec
             input shouldBe transformerC.uuid
             currentMagnitude should equalWithTolerance(
               Quantities.getQuantity(
-                0.7071067,
+                14.14213562,
                 StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
               ),
               1e-6
             )
             currentAngle should equalWithTolerance(
               Quantities
-                .getQuantity(-45.0000000, StandardUnits.ELECTRIC_CURRENT_ANGLE),
+                .getQuantity(135.0000000, StandardUnits.ELECTRIC_CURRENT_ANGLE),
               1e-6
             )
           case wrong => fail(s"Got wrong result: '$wrong'")
@@ -641,14 +603,14 @@ class GridResultsSupportSpec
             input shouldBe transformer3wInput.getUuid
             currentMagnitude should equalWithTolerance(
               Quantities.getQuantity(
-                0.00092397,
+                11.4542161,
                 StandardUnits.ELECTRIC_CURRENT_MAGNITUDE
               ),
               1e-6
             )
             currentAngle should equalWithTolerance(
               Quantities
-                .getQuantity(-72.648555, StandardUnits.ELECTRIC_CURRENT_ANGLE),
+                .getQuantity(-89.4475391, StandardUnits.ELECTRIC_CURRENT_ANGLE),
               1e-6
             )
           case wrong => fail(s"Got wrong result: '$wrong'")