Add sensitivities on intensity

src/main/java/com/powsybl/openloadflow/OpenLoadFlowProvider.java

@@ -47,10 +47,7 @@
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 
-import java.util.ArrayList;
-import java.util.Collections;
-import java.util.List;
-import java.util.Objects;
+import java.util.*;
 import java.util.concurrent.CompletableFuture;
 import java.util.concurrent.TimeUnit;
 
@@ -120,11 +117,12 @@ private static SlackBusSelector getSlackBusSelector(Network network, LoadFlowPar
 
     public static AcLoadFlowParameters createAcParameters(Network network, MatrixFactory matrixFactory, LoadFlowParameters parameters,
                                                           OpenLoadFlowParameters parametersExt, boolean breakers) {
-        return createAcParameters(network, matrixFactory, parameters, parametersExt, breakers, false);
+        return createAcParameters(network, matrixFactory, parameters, parametersExt, breakers, false, null);
     }
 
     public static AcLoadFlowParameters createAcParameters(Network network, MatrixFactory matrixFactory, LoadFlowParameters parameters,
-                                                   OpenLoadFlowParameters parametersExt, boolean breakers, boolean forceA1Var) {
+                                                          OpenLoadFlowParameters parametersExt, boolean breakers, boolean forceA1Var,
+                                                          Set<String> branchesWithCurrent) {
 
         SlackBusSelector slackBusSelector = getSlackBusSelector(network, parameters, parametersExt);
 
@@ -173,7 +171,8 @@ public static AcLoadFlowParameters createAcParameters(Network network, MatrixFac
                                         breakers,
                                         parametersExt.getPlausibleActivePowerLimit(),
                                         forceA1Var,
-                                        parametersExt.isAddRatioToLinesWithDifferentNominalVoltageAtBothEnds());
+                                        parametersExt.isAddRatioToLinesWithDifferentNominalVoltageAtBothEnds(),
+                                        branchesWithCurrent);
     }
 
     private LoadFlowResult runAc(Network network, LoadFlowParameters parameters, OpenLoadFlowParameters parametersExt) {

src/main/java/com/powsybl/openloadflow/ac/equations/AbstractBranchAcFlowEquationTerm.java

@@ -19,6 +19,8 @@
  */
 abstract class AbstractBranchAcFlowEquationTerm extends AbstractNamedEquationTerm {
 
+    public static final double CURRENT_NORMALIZATION_FACTOR = 1000d / Math.sqrt(3d);
+
     protected final LfBranch branch;
 
     protected final double b1;

src/main/java/com/powsybl/openloadflow/ac/equations/AbstractOpenBranchAcFlowEquationTerm.java → src/main/java/com/powsybl/openloadflow/ac/equations/AbstractOpenSide1BranchAcFlowEquationTerm.java

@@ -16,15 +16,15 @@
 import java.util.List;
 
 /**
- * @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
+ * @author Gael Macherel <gael.macherel at artelys.com>
  */
-abstract class AbstractOpenBranchAcFlowEquationTerm extends AbstractBranchAcFlowEquationTerm {
+abstract class AbstractOpenSide1BranchAcFlowEquationTerm extends AbstractBranchAcFlowEquationTerm {
 
     protected final List<Variable> variables;
 
     protected double shunt;
 
-    protected AbstractOpenBranchAcFlowEquationTerm(LfBranch branch, VariableType variableType,
+    protected AbstractOpenSide1BranchAcFlowEquationTerm(LfBranch branch, VariableType variableType,
                                                    LfBus bus, VariableSet variableSet,
                                                    boolean deriveA1, boolean deriveR1) {
         super(branch);

src/main/java/com/powsybl/openloadflow/ac/equations/AbstractOpenSide2BranchAcFlowEquationTerm.java

@@ -0,0 +1,43 @@
+/**
+ * Copyright (c) 2019, 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/.
+ */
+package com.powsybl.openloadflow.ac.equations;
+
+import com.powsybl.openloadflow.equations.Variable;
+import com.powsybl.openloadflow.equations.VariableSet;
+import com.powsybl.openloadflow.equations.VariableType;
+import com.powsybl.openloadflow.network.LfBranch;
+import com.powsybl.openloadflow.network.LfBus;
+
+import java.util.Collections;
+import java.util.List;
+
+/**
+ * @author Gael Macherel <gael.macherel at artelys.com>
+ */
+abstract class AbstractOpenSide2BranchAcFlowEquationTerm extends AbstractBranchAcFlowEquationTerm {
+
+    protected final List<Variable> variables;
+
+    protected double shunt;
+
+    protected AbstractOpenSide2BranchAcFlowEquationTerm(LfBranch branch, VariableType variableType,
+                                                        LfBus bus, VariableSet variableSet,
+                                                        boolean deriveA1, boolean deriveR1) {
+        super(branch);
+        variables = Collections.singletonList(variableSet.getVariable(bus.getNum(), variableType));
+        shunt = (g2 + y * sinKsi) * (g2 + y * sinKsi) + (-b2 + y * cosKsi) * (-b2 + y * cosKsi);
+
+        if (deriveA1 || deriveR1) {
+            throw new IllegalArgumentException("Variable A1 or R1 on open branch not supported: " + branch.getId());
+        }
+    }
+
+    @Override
+    public List<Variable> getVariables() {
+        return variables;
+    }
+}

src/main/java/com/powsybl/openloadflow/ac/equations/AcEquationSystem.java

@@ -264,21 +264,34 @@ private static void createImpedantBranch(LfBranch branch, LfBus bus1, LfBus bus2
         EquationTerm q1 = null;
         EquationTerm p2 = null;
         EquationTerm q2 = null;
+        EquationTerm i1 = null;
+        EquationTerm i2 = null;
         boolean deriveA1 = creationParameters.isPhaseControl() && branch.isPhaseController()
                 && branch.getDiscretePhaseControl().getMode() == DiscretePhaseControl.Mode.CONTROLLER;
         deriveA1 = deriveA1 || (creationParameters.isForceA1Var() && branch.hasPhaseControlCapability());
+        boolean createCurrent = creationParameters.getBranchesWithCurrent() == null || creationParameters.getBranchesWithCurrent().contains(branch.getId());
         boolean deriveR1 = creationParameters.isTransformerVoltageControl() && branch.isVoltageController();
         if (bus1 != null && bus2 != null) {
             p1 = new ClosedBranchSide1ActiveFlowEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
             q1 = new ClosedBranchSide1ReactiveFlowEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
             p2 = new ClosedBranchSide2ActiveFlowEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
             q2 = new ClosedBranchSide2ReactiveFlowEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
+            if (createCurrent) {
+                i1 = new ClosedBranchSide1CurrentMagnitudeEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
+                i2 = new ClosedBranchSide2CurrentMagnitudeEquationTerm(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
+            }
         } else if (bus1 != null) {
             p1 = new OpenBranchSide2ActiveFlowEquationTerm(branch, bus1, variableSet, deriveA1, deriveR1);
             q1 = new OpenBranchSide2ReactiveFlowEquationTerm(branch, bus1, variableSet, deriveA1, deriveR1);
+            if (createCurrent) {
+                i1 = new OpenBranchSide2CurrentMagnitudeEquationTerm(branch, bus1, variableSet, deriveA1, deriveR1);
+            }
         } else if (bus2 != null) {
             p2 = new OpenBranchSide1ActiveFlowEquationTerm(branch, bus2, variableSet, deriveA1, deriveR1);
             q2 = new OpenBranchSide1ReactiveFlowEquationTerm(branch, bus2, variableSet, deriveA1, deriveR1);
+            if (createCurrent) {
+                i2 = new OpenBranchSide1CurrentMagnitudeEquationTerm(branch, bus2, variableSet, deriveA1, deriveR1);
+            }
         }
 
         if (p1 != null) {
@@ -308,6 +321,18 @@ private static void createImpedantBranch(LfBranch branch, LfBus bus1, LfBus bus2
             equationSystem.createEquation(branch.getNum(), EquationType.BRANCH_ALPHA1)
                 .addTerm(new BranchA1EquationTerm(branch, variableSet));
         }
+
+        if (i1 != null) {
+            Equation i =  equationSystem.createEquation(bus1.getNum(), EquationType.BUS_I).addTerm(i1);
+            i.setUpdateType(EquationSystem.EquationUpdateType.AFTER_NR); // only update those equations after the newton raphson
+            branch.setI1(i1);
+        }
+
+        if (i2 != null) {
+            Equation i =  equationSystem.createEquation(bus2.getNum(), EquationType.BUS_I).addTerm(i2);
+            i.setUpdateType(EquationSystem.EquationUpdateType.AFTER_NR); // only update those equations after the newton raphson
+            branch.setI2(i2);
+        }
     }
 
     private static void createBranchEquations(LfNetwork network, VariableSet variableSet, AcEquationSystemCreationParameters creationParameters,

src/main/java/com/powsybl/openloadflow/ac/equations/AcEquationSystemCreationParameters.java

@@ -6,6 +6,8 @@
  */
 package com.powsybl.openloadflow.ac.equations;
 
+import java.util.Set;
+
 /**
  * @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
  */
@@ -17,14 +19,21 @@ public class AcEquationSystemCreationParameters {
 
     private final boolean forceA1Var;
 
+    private final Set<String> branchesWithCurrent;
+
     public AcEquationSystemCreationParameters(boolean phaseControl, boolean transformerVoltageControl) {
         this(phaseControl, transformerVoltageControl, false);
     }
 
     public AcEquationSystemCreationParameters(boolean phaseControl, boolean transformerVoltageControl, boolean forceA1Var) {
+        this(phaseControl, transformerVoltageControl, forceA1Var, null);
+    }
+
+    public AcEquationSystemCreationParameters(boolean phaseControl, boolean transformerVoltageControl, boolean forceA1Var, Set<String> branchesWithCurrent) {
         this.phaseControl = phaseControl;
         this.transformerVoltageControl = transformerVoltageControl;
         this.forceA1Var = forceA1Var;
+        this.branchesWithCurrent = branchesWithCurrent;
     }
 
     public boolean isPhaseControl() {
@@ -38,4 +47,8 @@ public boolean isTransformerVoltageControl() {
     public boolean isForceA1Var() {
         return forceA1Var;
     }
+
+    public Set<String> getBranchesWithCurrent() {
+        return branchesWithCurrent;
+    }
 }

src/main/java/com/powsybl/openloadflow/ac/equations/ClosedBranchSide1CurrentMagnitudeEquationTerm.java

@@ -0,0 +1,117 @@
+/**
+ * Copyright (c) 2019, 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/.
+ */
+package com.powsybl.openloadflow.ac.equations;
+
+import com.powsybl.openloadflow.equations.Variable;
+import com.powsybl.openloadflow.equations.VariableSet;
+import com.powsybl.openloadflow.network.LfBranch;
+import com.powsybl.openloadflow.network.LfBus;
+import net.jafama.FastMath;
+
+import java.util.Objects;
+
+import static com.powsybl.openloadflow.network.PiModel.A2;
+import static com.powsybl.openloadflow.network.PiModel.R2;
+
+/**
+ * @author Gael Macherel <gael.macherel at artelys.com>
+ */
+public class ClosedBranchSide1CurrentMagnitudeEquationTerm extends AbstractClosedBranchAcFlowEquationTerm {
+
+    private double i1;
+
+    private double di1dv1;
+
+    private double di1dv2;
+
+    private double di1dph1;
+
+    private double di1dph2;
+
+    private double di1da1;
+
+    public ClosedBranchSide1CurrentMagnitudeEquationTerm(LfBranch branch, LfBus bus1, LfBus bus2, VariableSet variableSet,
+                                                         boolean deriveA1, boolean deriveR1) {
+        super(branch, bus1, bus2, variableSet, deriveA1, deriveR1);
+    }
+
+    @Override
+    protected double calculateDer(double dph1, double dph2, double dv1, double dv2, double a1, double r1) {
+        return di1dph1 * dph1 + di1dph2 * dph2 + di1dv1 * dv1 + di1dv2 * dv2;
+    }
+
+    @Override
+    public void update(double[] x) {
+        Objects.requireNonNull(x);
+        double v1 = x[v1Var.getRow()];
+        double v2 = x[v2Var.getRow()];
+        double ph1 = x[ph1Var.getRow()];
+        double ph2 = x[ph2Var.getRow()];
+        double r1 = r1Var != null && r1Var.isActive() ? x[r1Var.getRow()] : branch.getPiModel().getR1();
+        double w1 = r1 * v1;
+        double w2 = y * R2 * v2;
+        double cosPh1 = FastMath.cos(ph1);
+        double sinPh1 = FastMath.sin(ph1);
+        double cosPh1Ksi = FastMath.cos(ph1 + ksi);
+        double sinPh1Ksi = FastMath.sin(ph1 + ksi);
+        double theta = ksi - (a1Var != null && a1Var.isActive() ? x[a1Var.getRow()] : branch.getPiModel().getA1())
+            + A2 + ph2;
+        double sinTheta = FastMath.sin(theta);
+        double cosTheta = FastMath.cos(theta);
+
+        double reI1 = r1 * (w1 * (g1 * cosPh1 - b1 * sinPh1 + y * sinPh1Ksi) - w2 * sinTheta) * CURRENT_NORMALIZATION_FACTOR;
+        double imI1 = r1 * (w1 * (g1 * sinPh1 + b1 * cosPh1 - y * cosPh1Ksi) + w2 * cosTheta) * CURRENT_NORMALIZATION_FACTOR;
+        i1 = Math.hypot(reI1, imI1);
+
+        double dreI1dv1 = r1 * r1 * (g1 * cosPh1 - b1 * sinPh1 + y * sinPh1Ksi) * CURRENT_NORMALIZATION_FACTOR;
+        double dreI1dv2 = r1 * (-y * R2 * sinTheta) * CURRENT_NORMALIZATION_FACTOR;
+        double dreI1dph1 = r1 * w1 * (-g1 * sinPh1 - b1 * cosPh1 + y * cosPh1Ksi) * CURRENT_NORMALIZATION_FACTOR;
+        double dreI1dph2 = r1 * (-w2 * cosTheta) * CURRENT_NORMALIZATION_FACTOR;
+
+        double dimI1dv1 = r1 * r1 * (g1 * sinPh1 + b1 * cosPh1 - y * cosPh1Ksi) * CURRENT_NORMALIZATION_FACTOR;
+        double dimI1dv2 = r1 * (y * R2 * cosTheta) * CURRENT_NORMALIZATION_FACTOR;
+        double dimI1dph1 = r1 * w1 * (g1 * cosPh1 - b1 * sinPh1 + y * sinPh1Ksi) * CURRENT_NORMALIZATION_FACTOR;
+        double dimI1dph2 = r1 * (-w2 * sinTheta) * CURRENT_NORMALIZATION_FACTOR;
+
+        di1dv1 = (reI1 * dreI1dv1 + imI1 * dimI1dv1) / i1;
+        di1dv2 = (reI1 * dreI1dv2 + imI1 * dimI1dv2) / i1;
+        di1dph1 = (reI1 * dreI1dph1 + imI1 * dimI1dph1) / i1;
+        di1dph2 = (reI1 * dreI1dph2 + imI1 * dimI1dph2) / i1;
+
+        if (a1Var != null) {
+            di1da1 = -di1dph2;
+        }
+    }
+
+    @Override
+    public double eval() {
+        return i1;
+    }
+
+    @Override
+    public double der(Variable variable) {
+        Objects.requireNonNull(variable);
+        if (variable.equals(v1Var)) {
+            return di1dv1;
+        } else if (variable.equals(v2Var)) {
+            return di1dv2;
+        } else if (variable.equals(ph1Var)) {
+            return di1dph1;
+        } else if (variable.equals(ph2Var)) {
+            return di1dph2;
+        } else if (variable.equals(a1Var)) {
+            return di1da1;
+        } else {
+            throw new IllegalStateException("Unknown variable: " + variable);
+        }
+    }
+
+    @Override
+    protected String getName() {
+        return "ac_i_closed_1";
+    }
+}