[WIP] Move DCOPF in AMPL to DC initialization in Java

open-reac/pom.xml

@@ -101,6 +101,11 @@
             <groupId>com.powsybl</groupId>
             <artifactId>powsybl-loadflow-api</artifactId>
         </dependency>
+        <dependency>
+            <groupId>com.powsybl</groupId>
+            <artifactId>powsybl-open-loadflow</artifactId>
+            <version>${powsybl-open-loadflow.version}</version>
+        </dependency>
 
         <dependency>
             <groupId>org.junit.jupiter</groupId>
@@ -148,11 +153,5 @@
             <artifactId>powsybl-ieee-cdf-converter</artifactId>
             <scope>test</scope>
         </dependency>
-        <dependency>
-            <groupId>com.powsybl</groupId>
-            <artifactId>powsybl-open-loadflow</artifactId>
-            <version>${powsybl-open-loadflow.version}</version>
-            <scope>test</scope>
-        </dependency>
     </dependencies>
 </project>

open-reac/src/main/java/com/powsybl/openreac/OpenReacRunner.java

@@ -14,6 +14,18 @@
 import com.powsybl.computation.ComputationManager;
 import com.powsybl.computation.local.LocalComputationManager;
 import com.powsybl.iidm.network.Network;
+import com.powsybl.iidm.network.VariantManager;
+import com.powsybl.iidm.network.extensions.ActivePowerControlAdder;
+import com.powsybl.loadflow.LoadFlowParameters;
+import com.powsybl.math.matrix.MatrixFactory;
+import com.powsybl.math.matrix.SparseMatrixFactory;
+import com.powsybl.openloadflow.FullVoltageInitializer;
+import com.powsybl.openloadflow.ac.VoltageMagnitudeInitializer;
+import com.powsybl.openloadflow.dc.DcValueVoltageInitializer;
+import com.powsybl.openloadflow.dc.equations.DcApproximationType;
+import com.powsybl.openloadflow.network.*;
+import com.powsybl.openloadflow.network.impl.LfNetworkLoaderImpl;
+import com.powsybl.openloadflow.network.util.VoltageInitializer;
 import com.powsybl.openreac.parameters.OpenReacAmplIOFiles;
 import com.powsybl.openreac.parameters.input.OpenReacParameters;
 import com.powsybl.openreac.parameters.output.OpenReacResult;
@@ -27,6 +39,8 @@
  */
 public final class OpenReacRunner {
 
+    private static final String WITH_USER_INITIALIZATION_ID_VARIANT = "WithDcLoadFlow";
+
     private OpenReacRunner() {
     }
 
@@ -74,7 +88,30 @@ public static OpenReacResult run(Network network, String variantId, OpenReacPara
         checkParameters(network, variantId, parameters, config, manager, reportNode);
         AmplModel reactiveOpf = OpenReacModel.buildModel();
         OpenReacAmplIOFiles amplIoInterface = new OpenReacAmplIOFiles(parameters, amplExportConfig, network, config.isDebug(), Reports.createOpenReacReporter(reportNode, network.getId(), parameters.getObjective()));
-        AmplResults run = AmplModelRunner.run(network, variantId, reactiveOpf, manager, amplIoInterface);
+
+        AmplResults run;
+        // current voltage values are used to initialize ACOPF optimization
+        if (parameters.getVoltageInitialization() == OpenReacParameters.OpenReacVoltageInitialization.PREVIOUS_VALUES) {
+            run = AmplModelRunner.run(network, variantId, reactiveOpf, manager, amplIoInterface);
+
+        // initialize optimization with user option
+        } else {
+            // create new variant to avoid modifications of the network
+            VariantManager variantManager = network.getVariantManager();
+            variantManager.cloneVariant(variantId, WITH_USER_INITIALIZATION_ID_VARIANT);
+            variantManager.setWorkingVariant(WITH_USER_INITIALIZATION_ID_VARIANT);
+
+            // initialize the optimization with given process
+            initializeVoltageBeforeOptimization(network, parameters);
+
+            // execute ampl code
+            run = AmplModelRunner.run(network, WITH_USER_INITIALIZATION_ID_VARIANT, reactiveOpf, manager, amplIoInterface);
+
+            // remove the variant created to store dc load flow results
+            variantManager.setWorkingVariant(variantId);
+            variantManager.removeVariant(WITH_USER_INITIALIZATION_ID_VARIANT);
+        }
+
         OpenReacResult result = new OpenReacResult(run.isSuccess() && amplIoInterface.checkErrors() ? OpenReacStatus.OK : OpenReacStatus.NOT_OK, amplIoInterface, run.getIndicators());
         Reports.createShuntModificationsReporter(reportNode, network.getId(), amplIoInterface.getNetworkModifications().getShuntsWithDeltaDiscreteOptimalOverThreshold());
         return result;
@@ -107,8 +144,29 @@ public static CompletableFuture<OpenReacResult> runAsync(Network network, String
         checkParameters(network, variantId, parameters, config, manager, reportNode);
         AmplModel reactiveOpf = OpenReacModel.buildModel();
         OpenReacAmplIOFiles amplIoInterface = new OpenReacAmplIOFiles(parameters, amplExportConfig, network, config.isDebug(), Reports.createOpenReacReporter(reportNode, network.getId(), parameters.getObjective()));
-        CompletableFuture<AmplResults> runAsync = AmplModelRunner.runAsync(network, variantId, reactiveOpf, manager, amplIoInterface);
+
+        CompletableFuture<AmplResults> runAsync;
+        if (parameters.getVoltageInitialization() == OpenReacParameters.OpenReacVoltageInitialization.PREVIOUS_VALUES) {
+            runAsync = AmplModelRunner.runAsync(network, variantId, reactiveOpf, manager, amplIoInterface);
+        } else {
+            // create new variant to avoid modifications of the network
+            VariantManager variantManager = network.getVariantManager();
+            variantManager.cloneVariant(variantId, WITH_USER_INITIALIZATION_ID_VARIANT);
+            variantManager.setWorkingVariant(WITH_USER_INITIALIZATION_ID_VARIANT);
+
+            // initialize the optimization with given process
+            initializeVoltageBeforeOptimization(network, parameters);
+
+            // execute ampl code
+            runAsync = AmplModelRunner.runAsync(network, WITH_USER_INITIALIZATION_ID_VARIANT, reactiveOpf, manager, amplIoInterface);
+        }
+
         return runAsync.thenApply(run -> {
+            if (parameters.getVoltageInitialization() != OpenReacParameters.OpenReacVoltageInitialization.PREVIOUS_VALUES) {
+                // remove the variant created to store dc load flow results
+                network.getVariantManager().setWorkingVariant(variantId);
+                network.getVariantManager().removeVariant(WITH_USER_INITIALIZATION_ID_VARIANT);
+            }
             OpenReacResult result = new OpenReacResult(run.isSuccess() && amplIoInterface.checkErrors() ? OpenReacStatus.OK : OpenReacStatus.NOT_OK, amplIoInterface, run.getIndicators());
             Reports.createShuntModificationsReporter(reportNode, network.getId(), amplIoInterface.getNetworkModifications().getShuntsWithDeltaDiscreteOptimalOverThreshold());
             return result;
@@ -124,4 +182,62 @@ private static void checkParameters(Network network, String variantId, OpenReacP
         Objects.requireNonNull(reportNode);
         parameters.checkIntegrity(network, Reports.createParameterIntegrityReporter(reportNode, network.getId()));
     }
+
+    public static void initializeVoltageBeforeOptimization(Network network, OpenReacParameters openReacParameters) {
+        // gsk to only distribute on generators
+        network.getGeneratorStream()
+                .forEach(generator -> generator.newExtension(ActivePowerControlAdder.class)
+                        .withParticipate(true)
+                        .withParticipationFactor(generator.getTargetP()));
+
+        // slack bus selection
+        SlackBusSelector slackBusSelector = new MostMeshedSlackBusSelector();
+
+        // get lfNetwork to apply voltage initialization
+        LfNetwork lfNetwork = LfNetwork.load(network, new LfNetworkLoaderImpl(), slackBusSelector).get(0);
+
+        // get initializer depending on user option
+        MatrixFactory matrixFactory = new SparseMatrixFactory();
+        switch (openReacParameters.getVoltageInitialization()) {
+            // uniform voltage initializer, for flat start voltage angles
+            case UNIFORM_VALUES -> {
+                for (LfBus lfBus : lfNetwork.getBuses()) {
+                    lfBus.setV(1);
+                    lfBus.setAngle(0);
+                }
+            }
+            // direct current initializer, to initialize voltage angles
+            case DC_VALUES -> {
+                VoltageInitializer initializer = new DcValueVoltageInitializer(new LfNetworkParameters().setSlackBusSelector(slackBusSelector),
+                        true,
+                        LoadFlowParameters.BalanceType.PROPORTIONAL_TO_GENERATION_PARTICIPATION_FACTOR,
+                        true,
+                        DcApproximationType.IGNORE_R,
+                        matrixFactory,
+                        1);
+                initializer.prepare(lfNetwork);
+            }
+            // full voltage initializer, to initialize voltage magnitudes and angles
+            case FULL_VOLTAGE -> {
+                VoltageInitializer initializer = new FullVoltageInitializer(
+                        new VoltageMagnitudeInitializer(false, matrixFactory, openReacParameters.getLowImpedanceThreshold()),
+                        new DcValueVoltageInitializer(new LfNetworkParameters().setSlackBusSelector(slackBusSelector),
+                                true,
+                                LoadFlowParameters.BalanceType.PROPORTIONAL_TO_GENERATION_PARTICIPATION_FACTOR,
+                                true,
+                                DcApproximationType.IGNORE_R,
+                                matrixFactory,
+                                1));
+                initializer.prepare(lfNetwork);
+            }
+            default -> throw new IllegalStateException("Unexpected value: " + openReacParameters.getVoltageInitialization());
+        }
+
+        // initialize voltage values
+
+        // update the network with initialization
+        LfNetworkStateUpdateParameters updateParameters = new LfNetworkStateUpdateParameters(false, false, true, false, false,
+                false, true, false, ReactivePowerDispatchMode.Q_EQUAL_PROPORTION, false, ReferenceBusSelectionMode.FIRST_SLACK, false);
+        lfNetwork.updateState(updateParameters);
+    }
 }

open-reac/src/main/java/com/powsybl/openreac/parameters/input/OpenReacParameters.java

@@ -127,6 +127,15 @@ public class OpenReacParameters {
 
     private double shuntCompensatorActivationAlertThreshold;
 
+    private OpenReacVoltageInitialization voltageInitialization = OpenReacVoltageInitialization.FULL_VOLTAGE;
+
+    public enum OpenReacVoltageInitialization {
+        PREVIOUS_VALUES,
+        UNIFORM_VALUES,
+        DC_VALUES,
+        FULL_VOLTAGE;
+    }
+
     /**
      * Override some voltage level limits in the network. This will NOT modify the network object.
      * <p>
@@ -548,6 +557,15 @@ public OpenReacParameters setShuntCompensatorActivationAlertThreshold(double shu
         return this;
     }
 
+    public OpenReacVoltageInitialization getVoltageInitialization() {
+        return voltageInitialization;
+    }
+
+    public OpenReacParameters setVoltageInitialization(OpenReacVoltageInitialization voltageInitialization) {
+        this.voltageInitialization = Objects.requireNonNull(voltageInitialization);
+        return this;
+    }
+
     public List<OpenReacAlgoParam> getAllAlgorithmParams() {
         ArrayList<OpenReacAlgoParam> allAlgoParams = new ArrayList<>();
         allAlgoParams.add(objective.toParam());

open-reac/src/main/java/com/powsybl/openreac/parameters/input/json/OpenReacParametersDeserializer.java

@@ -151,6 +151,10 @@ public OpenReacParameters deserialize(JsonParser parser, DeserializationContext
                     parser.nextToken();
                     parameters.setShuntVariableScalingFactor(parser.readValueAs(Double.class));
                 }
+                case "voltageInitialization" -> {
+                    parser.nextToken();
+                    parameters.setVoltageInitialization(OpenReacParameters.OpenReacVoltageInitialization.valueOf(parser.getText()));
+                }
                 default -> throw new IllegalStateException("Unexpected field: " + parser.getCurrentName());
             }
         }

open-reac/src/main/java/com/powsybl/openreac/parameters/input/json/OpenReacParametersSerializer.java

@@ -62,6 +62,7 @@ public void serialize(OpenReacParameters openReacParameters, JsonGenerator jsonG
         serializerProvider.defaultSerializeField("reactiveSlackVariableScalingFactor", openReacParameters.getReactiveSlackVariableScalingFactor(), jsonGenerator);
         serializerProvider.defaultSerializeField("twoWindingTransformerRatioVariableScalingFactor", openReacParameters.getTwoWindingTransformerRatioVariableScalingFactor(), jsonGenerator);
         serializerProvider.defaultSerializeField("shuntVariableScalingFactor", openReacParameters.getShuntVariableScalingFactor(), jsonGenerator);
+        serializerProvider.defaultSerializeField("voltageInitialization", openReacParameters.getVoltageInitialization(), jsonGenerator);
         jsonGenerator.writeEndObject();
     }
 }

open-reac/src/main/resources/openreac/acopf.mod

@@ -30,6 +30,8 @@ set PROBLEM_ACOPF default { };
 # Complex voltage = V*exp(i*teta). (with i**2=-1)
 
 # Phase of voltage
+param teta_min default -3; # roughly 3*pi
+param teta_max default  3; # roughly 3*pi
 var teta{BUSCC} <= teta_max, >= teta_min;
 subject to ctr_null_phase_bus{PROBLEM_ACOPF}: teta[null_phase_bus] = 0;
 

open-reac/src/main/resources/openreac/acopf.run

@@ -20,7 +20,12 @@
 ###############################################################################
 
 # DC values used as initialization for phases
-let {n in BUSCC} teta[n] := teta_dc[n];
+
+let teta_max := max({n in BUSCC : -10 <= bus_angl0[1,n] and bus_angl0[1,n] <= 10} bus_angl0[1,n]) + 3;
+let teta_min := min({n in BUSCC : -10 <= bus_angl0[1,n] and bus_angl0[1,n] <= 10} bus_angl0[1,n]) - 3;
+let {n in BUSCC} teta[n] := if teta_min <= bus_angl0[1,n] and bus_angl0[1,n] <= teta_max
+                            then bus_angl0[1,n]
+                            else 0;
 
 # Voltages
 let {n in BUSVV}            V[n] := bus_V0[1,n];
@@ -30,7 +35,7 @@ let {n in BUSCC diff BUSVV} V[n] := voltage_lower_bound[1,bus_substation[1,n]]
 let {n in BUSCC_SLACK} slack1_shunt_B[n] := 0;
 let {n in BUSCC_SLACK} slack2_shunt_B[n] := 0;
 let alpha := 0.01;
-let {(g,n) in UNITON} P_bounded[g,n] := max(P_dcopf[g,n],unit_Pc[1,g,n]);
+let {(g,n) in UNITON} P_bounded[g,n] := max(unit_P0[1,g,n],unit_Pc[1,g,n]);;
 let {(g,n) in UNITON} Q[g,n] := 0.5*(corrected_unit_Qmax[g,n] + corrected_unit_Qmin[g,n]);
 let {(shunt,n) in SHUNT_VAR} shunt_var[shunt,n] := min{(1,shunt,k) in SHUNT} shunt_valnom[1,shunt,k];
 let {(svc,n) in SVCON} svc_qvar[svc,n] := 0.1;
@@ -118,10 +123,10 @@ printf{LOG_INFO} "** ACopf results\n";
 if 1 in LOG_INFO then display
   nb_iter_last,nb_iter_total,
   max({(qq,m,n) in BRANCHCC} branch_R[1,qq,m,n]),max({(qq,m,n) in BRANCHCC} branch_X[1,qq,m,n]),
-  teta_max, max({n in BUSCC} teta[n]), max({n in BUSCC} teta_dc[n]),
-  teta_min, min({n in BUSCC} teta[n]), min({n in BUSCC} teta_dc[n]),
-  max({(qq,m,n) in BRANCHCC} (teta[m]-teta[n])), max({(qq,m,n) in BRANCHCC} (teta_dc[m]-teta_dc[n])),
-  min({(qq,m,n) in BRANCHCC} (teta[m]-teta[n])), min({(qq,m,n) in BRANCHCC} (teta_dc[m]-teta_dc[n])),
+  teta_max, max({n in BUSCC} teta[n]),# max({n in BUSCC} teta_dc[n]),
+  teta_min, min({n in BUSCC} teta[n]),# min({n in BUSCC} teta_dc[n]),
+  max({(qq,m,n) in BRANCHCC} (teta[m]-teta[n])),# max({(qq,m,n) in BRANCHCC} (teta_dc[m]-teta_dc[n])),
+  min({(qq,m,n) in BRANCHCC} (teta[m]-teta[n])),# min({(qq,m,n) in BRANCHCC} (teta_dc[m]-teta_dc[n])),
   min({n in BUSCC}V[n]),max({n in BUSCC}V[n])
   ;