[WIP] Avoid angle divergence with HVDC AC emulation and PMax

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

@@ -37,6 +37,10 @@ public abstract class AbstractHvdcAcEmulationFlowEquationTerm extends AbstractEl
 
     protected final double pMaxFromCS2toCS1;
 
+    protected final double tetaMax;
+    protected final double tetaMin;
+    protected final double tetaZero;
+
     protected AbstractHvdcAcEmulationFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus bus2, VariableSet<AcVariableType> variableSet) {
         super(hvdc);
         ph1Var = variableSet.getVariable(bus1.getNum(), AcVariableType.BUS_PHI);
@@ -48,9 +52,12 @@ protected AbstractHvdcAcEmulationFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus
         lossFactor2 = hvdc.getConverterStation2().getLossFactor() / 100;
         pMaxFromCS1toCS2 = hvdc.getPMaxFromCS1toCS2();
         pMaxFromCS2toCS1 = hvdc.getPMaxFromCS2toCS1();
+        tetaMax = (pMaxFromCS1toCS2 - p0) / k;
+        tetaMin = (-pMaxFromCS2toCS1 - p0) / k;
+        tetaZero = -p0 / k;
     }
 
-    protected double rawP(double p0, double k, double ph1, double ph2) {
+    protected double rawP(double ph1, double ph2) {
         return p0 + k * (ph1 - ph2);
     }
 
@@ -65,6 +72,28 @@ protected double boundedP(double rawP) {
         }
     }
 
+    /**
+     * Returns a "corrected value for k that is
+     *     the droop factor between tetaMin and tetaMax
+     *     0 in a resonable range above tetaMin and tetaMax
+     *     the slot between pMax and P0 for larger angle differences to help convergence
+     * @param ph1
+     * @param ph2
+     * @return
+     */
+    protected double dummySlope(double rawP, double ph1, double ph2) {
+        double boundedP = boundedP(rawP);
+        double factor = k;
+        double teta = ph1 - ph2;
+        // for large values of teta return a value that helps convergence
+        if (teta > tetaMax) {
+            factor = teta < tetaMax * 2 ? 0 : boundedP / (teta - tetaZero);
+        } else if (teta < tetaMin) {
+            factor = teta > tetaMin * 2 ? 0 : boundedP / (teta - tetaZero);
+        }
+        return factor;
+    }
+
     protected double ph1() {
         return sv.get(ph1Var.getRow());
     }

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

@@ -23,7 +23,7 @@ public HvdcAcEmulationSide1ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus
     }
 
     private double p1(double ph1, double ph2) {
-        double rawP = rawP(p0, k, ph1, ph2);
+        double rawP = rawP(ph1, ph2);
         double boundedP = boundedP(rawP);
         return (isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP;
     }
@@ -37,12 +37,12 @@ private boolean isInOperatingRange(double rawP) {
     }
 
     protected double dp1dph1(double ph1, double ph2) {
-        double rawP = rawP(p0, k, ph1, ph2);
-        if (isInOperatingRange(rawP)) {
-            return (isController(rawP) ? 1 : getVscLossMultiplier()) * k;
-        } else {
-            return 0;
-        }
+        // return the exact derivative unless ph1-oh2 us ti large (twice value of P limit)
+        // in this case returns a dummy slope to heo th NR convergence towards desired P
+        double rawP = rawP(ph1, ph2);
+        double factor = dummySlope(rawP, ph1, ph2);
+        return (isController(rawP) ? 1 : getVscLossMultiplier()) * factor;
+
     }
 
     protected double dp1dph2(double ph1, double ph2) {

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

@@ -23,7 +23,7 @@ public HvdcAcEmulationSide2ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus
     }
 
     private double p2(double ph1, double ph2) {
-        double rawP = rawP(p0, k, ph1, ph2);
+        double rawP = rawP(ph1, ph2);
         double boundedP = boundedP(rawP);
         return -(isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP;
     }
@@ -37,12 +37,11 @@ private boolean isInOperatingRange(double rawP) {
     }
 
     private double dp2dph1(double ph1, double ph2) {
-        double rawP = rawP(p0, k, ph1, ph2);
-        if (isInOperatingRange(rawP)) {
-            return -(isController(rawP) ? 1 : getVscLossMultiplier()) * k;
-        } else {
-            return 0;
-        }
+        // return the exact derivative unless ph1-oh2 us ti large (twice value of P limit)
+        // in this case returns a dummy slope to heo th NR convergence towards desired P
+        double rawP = rawP(ph1, ph2);
+        double factor = dummySlope(rawP, ph1, ph2);
+        return -(isController(rawP) ? 1 : getVscLossMultiplier()) * factor;
     }
 
     private double dp2dph2(double ph1, double ph2) {

src/test/java/com/powsybl/openloadflow/ac/AcLoadFlowVscTest.java

@@ -502,6 +502,44 @@ void testAcEmuWithOperationalLimits() {
         assertEquals(180, network.getHvdcConverterStation("cs3").getTerminal().getP(), DELTA_POWER);
     }
 
+    @Test
+    void testAcEmuWithOperationalLimitsWithDCInit() {
+        Network network = HvdcNetworkFactory.createHvdcLinkedByTwoLinesAndSwitch(HvdcConverterStation.HvdcType.VSC);
+        // without limit p=195
+        network.getHvdcLine("hvdc23")
+                .newExtension(HvdcOperatorActivePowerRangeAdder.class)
+                .withOprFromCS2toCS1(180)
+                .withOprFromCS1toCS2(170)
+                .add();
+
+        LoadFlow.Runner loadFlowRunner = new LoadFlow.Runner(new OpenLoadFlowProvider(new DenseMatrixFactory()));
+        LoadFlowParameters p = new LoadFlowParameters();
+        p.setHvdcAcEmulation(true);
+        p.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.DC_VALUES);
+
+        // fool the DC initialization with large angle
+        network.getHvdcLine("hvdc23").setActivePowerSetpoint(0);
+        network.getLine("l14").setX(20);
+
+        LoadFlowResult result = loadFlowRunner.run(network, p);
+
+        assertTrue(result.isFullyConverged());
+
+        // Active flow capped at limit. Output has losses (due to VSC stations)
+        assertEquals(170, network.getHvdcConverterStation("cs2").getTerminal().getP(), DELTA_POWER);
+        assertEquals(-166.280, network.getHvdcConverterStation("cs3").getTerminal().getP(), DELTA_POWER);
+
+        // now invert the angle to test ill conditioned converence in the other direction
+        network.getLine("l14").setX(-20);
+        result = loadFlowRunner.run(network, p);
+        assertTrue(result.isFullyConverged());
+
+        // Active flow capped at other direction's limit. Output has losses (due to VSC stations)
+        assertEquals(-176.062, network.getHvdcConverterStation("cs2").getTerminal().getP(), DELTA_POWER);
+        assertEquals(180, network.getHvdcConverterStation("cs3").getTerminal().getP(), DELTA_POWER);
+
+    }
+
     @Test
     void testAcEmuAndPMax() {
         Network network = HvdcNetworkFactory.createHvdcLinkedByTwoLinesAndSwitch(HvdcConverterStation.HvdcType.VSC);