From 82d556b70b5e5c8da872f8fb1d3d1a6611e1781a Mon Sep 17 00:00:00 2001 From: "VIDAL Didier (Externe)" Date: Tue, 2 Apr 2024 16:44:15 +0200 Subject: [PATCH 1/3] Avoid angle divergence with HVDC AC emulation and PMax Signed-off-by: VIDAL Didier (Externe) --- ...stractHvdcAcEmulationFlowEquationTerm.java | 13 ++++++++++-- ...cEmulationSide1ActiveFlowEquationTerm.java | 20 ++++++++++++------- ...cEmulationSide2ActiveFlowEquationTerm.java | 19 +++++++++++------- 3 files changed, 36 insertions(+), 16 deletions(-) diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java index f9604b7804..64c287cbea 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java @@ -19,6 +19,8 @@ */ public abstract class AbstractHvdcAcEmulationFlowEquationTerm extends AbstractElementEquationTerm { + static final double EPSILON = 0.1; + protected final Variable ph1Var; protected final Variable ph2Var; @@ -37,6 +39,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 variableSet) { super(hvdc); ph1Var = variableSet.getVariable(bus1.getNum(), AcVariableType.BUS_PHI); @@ -48,13 +54,16 @@ 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); } - protected double boundedP(double rawP) { + protected double boundedP(double rawP, double ph1, double ph2) { // If there is a maximal active power // it is applied at the entry of the controller VSC station // on the AC side of the network. diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java index e9dc0be77e..c58291c0e0 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java @@ -23,8 +23,8 @@ public HvdcAcEmulationSide1ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus } private double p1(double ph1, double ph2) { - double rawP = rawP(p0, k, ph1, ph2); - double boundedP = boundedP(rawP); + double rawP = rawP(ph1, ph2); + double boundedP = boundedP(rawP, ph1, ph2); return (isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP; } @@ -37,12 +37,18 @@ 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; + double rawP = rawP(ph1, ph2); + double boundedP = boundedP(rawP, ph1, ph2); + 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 - p0) / (teta - tetaZero); + } else if (teta < tetaMin) { + factor = teta > tetaMin * 2 ? 0 : (p0 - boundedP) / (teta - tetaZero); } + return (isController(rawP) ? 1 : getVscLossMultiplier()) * factor; + } protected double dp1dph2(double ph1, double ph2) { diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java index e00b4837ea..17caac993a 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java @@ -23,8 +23,8 @@ public HvdcAcEmulationSide2ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus } private double p2(double ph1, double ph2) { - double rawP = rawP(p0, k, ph1, ph2); - double boundedP = boundedP(rawP); + double rawP = rawP(ph1, ph2); + double boundedP = boundedP(rawP, ph1, ph2); return -(isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP; } @@ -37,12 +37,17 @@ 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; + double rawP = rawP(ph1, ph2); + double boundedP = boundedP(rawP, ph1, ph2); + 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 - p0) / (teta - tetaZero); + } else if (teta < tetaMin) { + factor = teta > tetaMin * 2 ? 0 : (p0 - boundedP) / (teta - tetaZero); } + return -(isController(rawP) ? 1 : getVscLossMultiplier()) * factor; } private double dp2dph2(double ph1, double ph2) { From 4f7ea581d6db0dd237e0684458ee7ffef7f98292 Mon Sep 17 00:00:00 2001 From: "VIDAL Didier (Externe)" Date: Tue, 2 Apr 2024 17:09:41 +0200 Subject: [PATCH 2/3] factorize code Signed-off-by: VIDAL Didier (Externe) --- ...stractHvdcAcEmulationFlowEquationTerm.java | 26 ++++++++++++++++--- ...cEmulationSide1ActiveFlowEquationTerm.java | 12 ++------- ...cEmulationSide2ActiveFlowEquationTerm.java | 12 ++------- 3 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java index 64c287cbea..ec51462693 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java @@ -19,8 +19,6 @@ */ public abstract class AbstractHvdcAcEmulationFlowEquationTerm extends AbstractElementEquationTerm { - static final double EPSILON = 0.1; - protected final Variable ph1Var; protected final Variable ph2Var; @@ -63,7 +61,7 @@ protected double rawP(double ph1, double ph2) { return p0 + k * (ph1 - ph2); } - protected double boundedP(double rawP, double ph1, double ph2) { + protected double boundedP(double rawP) { // If there is a maximal active power // it is applied at the entry of the controller VSC station // on the AC side of the network. @@ -74,6 +72,28 @@ protected double boundedP(double rawP, double ph1, double ph2) { } } + /** + * 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 pseudoK(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()); } diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java index c58291c0e0..55103d7ff3 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java @@ -24,7 +24,7 @@ public HvdcAcEmulationSide1ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus private double p1(double ph1, double ph2) { double rawP = rawP(ph1, ph2); - double boundedP = boundedP(rawP, ph1, ph2); + double boundedP = boundedP(rawP); return (isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP; } @@ -38,15 +38,7 @@ private boolean isInOperatingRange(double rawP) { protected double dp1dph1(double ph1, double ph2) { double rawP = rawP(ph1, ph2); - double boundedP = boundedP(rawP, ph1, ph2); - 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 - p0) / (teta - tetaZero); - } else if (teta < tetaMin) { - factor = teta > tetaMin * 2 ? 0 : (p0 - boundedP) / (teta - tetaZero); - } + double factor = pseudoK(rawP, ph1, ph2); return (isController(rawP) ? 1 : getVscLossMultiplier()) * factor; } diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java index 17caac993a..2353c96d8e 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java @@ -24,7 +24,7 @@ public HvdcAcEmulationSide2ActiveFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus private double p2(double ph1, double ph2) { double rawP = rawP(ph1, ph2); - double boundedP = boundedP(rawP, ph1, ph2); + double boundedP = boundedP(rawP); return -(isController(rawP) ? 1 : getVscLossMultiplier()) * boundedP; } @@ -38,15 +38,7 @@ private boolean isInOperatingRange(double rawP) { private double dp2dph1(double ph1, double ph2) { double rawP = rawP(ph1, ph2); - double boundedP = boundedP(rawP, ph1, ph2); - 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 - p0) / (teta - tetaZero); - } else if (teta < tetaMin) { - factor = teta > tetaMin * 2 ? 0 : (p0 - boundedP) / (teta - tetaZero); - } + double factor = pseudoK(rawP, ph1, ph2); return -(isController(rawP) ? 1 : getVscLossMultiplier()) * factor; } From 33c06d3c48d11ee7a418aa06f80aeed9b1084b86 Mon Sep 17 00:00:00 2001 From: "VIDAL Didier (Externe)" Date: Wed, 3 Apr 2024 07:58:57 +0200 Subject: [PATCH 3/3] add test coverage and better comment Signed-off-by: VIDAL Didier (Externe) --- ...stractHvdcAcEmulationFlowEquationTerm.java | 4 +- ...cEmulationSide1ActiveFlowEquationTerm.java | 4 +- ...cEmulationSide2ActiveFlowEquationTerm.java | 4 +- .../openloadflow/ac/AcLoadFlowVscTest.java | 38 +++++++++++++++++++ 4 files changed, 46 insertions(+), 4 deletions(-) diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java index ec51462693..c6b95bcd5a 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/AbstractHvdcAcEmulationFlowEquationTerm.java @@ -53,7 +53,7 @@ protected AbstractHvdcAcEmulationFlowEquationTerm(LfHvdc hvdc, LfBus bus1, LfBus pMaxFromCS1toCS2 = hvdc.getPMaxFromCS1toCS2(); pMaxFromCS2toCS1 = hvdc.getPMaxFromCS2toCS1(); tetaMax = (pMaxFromCS1toCS2 - p0) / k; - tetaMin = -(pMaxFromCS2toCS1 - p0) / k; + tetaMin = (-pMaxFromCS2toCS1 - p0) / k; tetaZero = -p0 / k; } @@ -81,7 +81,7 @@ protected double boundedP(double rawP) { * @param ph2 * @return */ - protected double pseudoK(double rawP, double ph1, double ph2) { + protected double dummySlope(double rawP, double ph1, double ph2) { double boundedP = boundedP(rawP); double factor = k; double teta = ph1 - ph2; diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java index 55103d7ff3..4f20108df7 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide1ActiveFlowEquationTerm.java @@ -37,8 +37,10 @@ private boolean isInOperatingRange(double rawP) { } protected double dp1dph1(double ph1, double ph2) { + // 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 = pseudoK(rawP, ph1, ph2); + double factor = dummySlope(rawP, ph1, ph2); return (isController(rawP) ? 1 : getVscLossMultiplier()) * factor; } diff --git a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java index 2353c96d8e..92159593f9 100644 --- a/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java +++ b/src/main/java/com/powsybl/openloadflow/ac/equations/HvdcAcEmulationSide2ActiveFlowEquationTerm.java @@ -37,8 +37,10 @@ private boolean isInOperatingRange(double rawP) { } private double dp2dph1(double ph1, double ph2) { + // 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 = pseudoK(rawP, ph1, ph2); + double factor = dummySlope(rawP, ph1, ph2); return -(isController(rawP) ? 1 : getVscLossMultiplier()) * factor; } diff --git a/src/test/java/com/powsybl/openloadflow/ac/AcLoadFlowVscTest.java b/src/test/java/com/powsybl/openloadflow/ac/AcLoadFlowVscTest.java index c6776d5f6d..4ee930062e 100644 --- a/src/test/java/com/powsybl/openloadflow/ac/AcLoadFlowVscTest.java +++ b/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);