powsybl · p-arvy · Oct 30, 2023
diff --git a/open-reac/src/main/java/com/powsybl/openreac/OpenReacModel.java b/open-reac/src/main/java/com/powsybl/openreac/OpenReacModel.java
@@ -52,7 +52,8 @@ public Charset getFileEncoding() {
     public static OpenReacModel buildModel() {
         return new OpenReacModel(OUTPUT_FILE_PREFIX, "openreac",
                 List.of("reactiveopf.run"),
-                List.of("reactiveopf.mod", "reactiveopf.dat", "reactiveopfoutput.run", "reactiveopfexit.run"));
+                List.of("reactiveopf.mod", "ccomputation.mod", "dcopf.mod", "acopf.mod",
+                        "reactiveopf.dat", "reactiveopfoutput.run", "reactiveopfexit.run"));
     }
 
     private static final String NETWORK_DATA_PREFIX = "ampl";

diff --git a/open-reac/src/main/resources/openreac/acopf.mod b/open-reac/src/main/resources/openreac/acopf.mod
@@ -0,0 +1,224 @@
+###############################################################################
+#
+# Copyright (c) 2022 2023, 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/.
+#
+###############################################################################
+
+###############################################################################
+# Reactive OPF
+# Author:  Jean Maeght 2022 2023
+###############################################################################
+
+# Definition of optimization problem
+set PROBLEM_ACOPF default { };
+
+
+###############################################################################
+# Voltage variables and constraints 
+###############################################################################
+
+## Complex voltage = V*exp(i*teta). (with i**2=-1)
+
+# Phase of voltage
+var teta{BUSCC} 
+  <= teta_max, 
+  >= teta_min;
+
+# Modulus of voltage
+var V{n in BUSCC}
+  <=
+  if substation_Vnomi[1,bus_substation[1,n]] <= ignore_voltage_bounds then 2-epsilon_min_voltage else
+  voltage_upper_bound[1,bus_substation[1,n]],
+  >=
+  if substation_Vnomi[1,bus_substation[1,n]] <= ignore_voltage_bounds then epsilon_min_voltage else
+  voltage_lower_bound[1,bus_substation[1,n]];
+
+subject to ctr_null_phase_bus{PROBLEM_ACOPF}: teta[null_phase_bus] = 0;
+
+
+###############################################################################
+# Active power variables and constraints
+###############################################################################
+
+# General idea: Generation is an input data, but as voltage may vary, generation may vary a little.
+# Variations of generation is totally controlled by unique scalar variable alpha
+# Before and after optimization, there is no waranty that P is within
+# its "bounds" [corrected_unit_Pmin;corrected_unit_Pmax]
+
+# Active generation
+var alpha <=1, >=-1; # If alpha==1 then all units are at Pmax
+var P_bounded{(g,n) in UNITON} <= max(unit_Pc[1,g,n],corrected_unit_Pmax[g,n]), >= min(unit_Pc[1,g,n],corrected_unit_Pmin[g,n]);
+# If coeff_alpha == 1 then all P are defined by the single variable alpha
+# If coeff_alpha == 0 then all P are free within their respective bounds
+# todo faire des tests avec les valeurs de coeff_alpha
+var P{(g,n) in UNITON} =
+  if   ( unit_Pc[1,g,n] < (corrected_unit_Pmax[g,n] - Pnull) and unit_Pc[1,g,n] > Pnull )
+  then (     coeff_alpha   * ( unit_Pc[1,g,n] + alpha*(corrected_unit_Pmax[g,n]- unit_Pc[1,g,n]) )
+         + (1-coeff_alpha) *   P_bounded[g,n] )
+  else unit_Pc[1,g,n] ;
+
+# Ratios of transformers
+var branch_Ror_var{(qq,m,n) in BRANCHCC_REGL_VAR}
+  >= regl_ratio_min[1,branch_ptrRegl[1,qq,m,n]],
+  <= regl_ratio_max[1,branch_ptrRegl[1,qq,m,n]];
+
+# Active flows on branches
+var Red_Tran_Act_Dir{(qq,m,n) in BRANCHCC } =
+    V[n] * branch_admi[qq,m,n] * sin(teta[m]-teta[n]+branch_dephor[qq,m,n]-branch_angper[qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])
+  + V[m] * (branch_admi[qq,m,n]*sin(branch_angper[qq,m,n])+branch_Gor[1,qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])**2
+  ;
+
+var Red_Tran_Act_Inv{(qq,m,n) in BRANCHCC } =
+    V[m] * branch_admi[qq,m,n] * sin(teta[n]-teta[m]-branch_dephor[qq,m,n]-branch_angper[qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])
+  + V[n] * (branch_admi[qq,m,n]*sin(branch_angper[qq,m,n])+branch_Gex[1,qq,m,n])
+  ;
+
+# Active power balance
+subject to ctr_balance_P{PROBLEM_ACOPF,k in BUSCC}:
+  # Flows
+    sum{(qq,k,n) in BRANCHCC} base100MVA * V[k] * Red_Tran_Act_Dir[qq,k,n]
+  + sum{(qq,m,k) in BRANCHCC} base100MVA * V[k] * Red_Tran_Act_Inv[qq,m,k]
+  # Generating units
+  - sum{(g,k) in UNITON} P[g,k]
+  # Batteries
+  - sum{(b,k) in BATTERYCC} battery_p0[1,b,k]
+  # Loads
+  + sum{(c,k) in LOADCC} load_PFix[1,c,k]     # Fixed value
+  # VSC converters
+  + sum{(v,k) in VSCCONVON} vscconv_P0[1,v,k] # Fixed value
+  # LCC converters
+  + sum{(l,k) in LCCCONVON} lccconv_P0[1,l,k] # Fixed value
+  = 0; # No slack variables for active power. If data are really too bad, may not converge.
+
+
+###############################################################################
+# Reactive power variables and constraints
+###############################################################################
+
+# Reactive generation of units
+var Q{(g,n) in UNITON} 
+  <= corrected_unit_Qmax[g,n], 
+  >= corrected_unit_Qmin[g,n];
+# todo: add trapeze or hexagone constraints for reactive power
+
+# Variable shunts
+var shunt_var{(shunt,n) in SHUNT_VAR}
+  >= min{(1,shunt,k) in SHUNT} shunt_valmin[1,shunt,k],
+  <= max{(1,shunt,k) in SHUNT} shunt_valmax[1,shunt,k];
+
+# SVC reactive generation
+var svc_qvar{(svc,n) in SVCON} 
+  >= svc_bmin[1,svc,n], 
+  <= svc_bmax[1,svc,n];
+
+# VSCCONV reactive generation
+var vscconv_qvar{(v,n) in VSCCONVON}
+  >= min(vscconv_qP[1,v,n],vscconv_qp0[1,v,n],vscconv_qp[1,v,n]),
+  <= max(vscconv_QP[1,v,n],vscconv_Qp0[1,v,n],vscconv_Qp[1,v,n]);
+# todo: add trapeze or hexagone constraints for reactive power
+
+# Reactive flows on branches
+var Red_Tran_Rea_Dir{(qq,m,n) in BRANCHCC } =
+  - V[n] * branch_admi[qq,m,n] * cos(teta[m]-teta[n]+branch_dephor[qq,m,n]-branch_angper[qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])
+  + V[m] * (branch_admi[qq,m,n]*cos(branch_angper[qq,m,n])-branch_Bor[1,qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])^2
+  ;
+
+var Red_Tran_Rea_Inv{(qq,m,n) in BRANCHCC } =
+  - V[m] * branch_admi[qq,m,n] * cos(teta[n]-teta[m]-branch_dephor[qq,m,n]-branch_angper[qq,m,n])
+    * (if (qq,m,n) in BRANCHCC_REGL_VAR then branch_Ror_var[qq,m,n]*branch_cstratio[1,qq,m,n] else branch_Ror[qq,m,n])
+  + V[n] * (branch_admi[qq,m,n]*cos(branch_angper[qq,m,n])-branch_Bex[1,qq,m,n])
+  ;
+
+# Reactive balance slack variables only if there is a load or a shunt connected
+# If there is a unit, or SVC, or VSC, they already have reactive power generation, so no need to add slack variables
+set BUSCC_SLACK :=  {n in BUSCC:(card{(g,n) in UNITON: (g,n) not in UNIT_FIXQ}==0
+                                    and card{(svc,n) in SVCON}==0
+                                    and card{(vscconv,n) in VSCCONVON}==0)} ;
+var slack1_balance_Q{BUSCC_SLACK} >=0, <= 500; # 500 Mvar is already HUGE
+var slack2_balance_Q{BUSCC_SLACK} >=0, <= 500;
+#subject to ctr_compl_slack_Q{PROBLEM_ACOPF,k in BUSCC_SLACK}: slack1_balance_Q[k] >= 0 complements slack2_balance_Q[k] >= 0;
+
+# Reactive power balance
+subject to ctr_balance_Q{PROBLEM_ACOPF,k in BUSCC}:
+  # Flows
+    sum{(qq,k,n) in BRANCHCC} base100MVA * V[k] * Red_Tran_Rea_Dir[qq,k,n]
+  + sum{(qq,m,k) in BRANCHCC} base100MVA * V[k] * Red_Tran_Rea_Inv[qq,m,k]
+  # Generating units
+  - sum{(g,k) in UNITON: (g,k) not in UNIT_FIXQ } Q[g,k]
+  - sum{(g,k) in UNIT_FIXQ} unit_Qc[1,g,k]
+  # Batteries
+  - sum{(b,k) in BATTERYCC} battery_q0[1,b,k]
+  # Loads
+  + sum{(c,k) in LOADCC} load_QFix[1,c,k]
+  # Shunts
+  - sum{(shunt,k) in SHUNT_FIX} base100MVA * shunt_valnom[1,shunt,k] * V[k]^2
+  - sum{(shunt,k) in SHUNT_VAR} base100MVA * shunt_var[shunt,k] * V[k]^2
+  # SVC
+  - sum{(svc,k) in SVCON} base100MVA * svc_qvar[svc,k] * V[k]^2
+  # VSC converters
+  - sum{(v,k) in VSCCONVON} vscconv_qvar[v,k]
+  # LCC converters
+  + sum{(l,k) in LCCCONVON} lccconv_Q0[1,l,k] # Fixed value
+  # Slack variables
+  + if k in BUSCC_SLACK then
+  (- slack1_balance_Q[k]  # Homogeneous to a generation of reactive power (condensator)
+   + slack2_balance_Q[k]) # homogeneous to a reactive load (self)
+  = 0;
+
+
+###############################################################################
+# Objective function and penalties
+###############################################################################
+
+# Voltage target : ratio between Vmin and Vmax
+var target_voltage_ratio = sum{n in BUSCC: substation_Vnomi[1,bus_substation[1,n]] > ignore_voltage_bounds}
+  ( V[n] - (1-ratio_voltage_target)*voltage_lower_bound[1,bus_substation[1,n]] + ratio_voltage_target*voltage_upper_bound[1,bus_substation[1,n]] )**2;
+
+# Voltage target : value V0 in input data
+var target_voltage_data = sum{n in BUSVV} (V[n] - bus_V0[1,n])**2;
+
+param penalty_invest_rea_pos := 10;
+param penalty_invest_rea_neg := 10;
+param penalty_units_reactive := 0.1;
+param penalty_transfo_ratio  := 0.1;
+
+param penalty_active_power_high := 1;
+param penalty_active_power_low  := 0.01;
+
+param penalty_voltage_target_high := 1;
+param penalty_voltage_target_low  := 0.01;
+
+minimize problem_acopf_objective:
+  sum{n in BUSCC_SLACK} (
+      penalty_invest_rea_pos * slack1_balance_Q[n]
+    + penalty_invest_rea_neg * slack2_balance_Q[n]
+    )
+
+  # coeff_alpha == 1 : minimize sum of generation, all generating units vary with 1 unique variable alpha
+  # coeff_alpha == 0 : minimize sum of squared difference between target and value
+  + (if objective_choice==1 or objective_choice==2 then penalty_active_power_low else penalty_active_power_high)
+  * sum{(g,n) in UNITON} (coeff_alpha * P[g,n] + (1-coeff_alpha)*( (P[g,n]-unit_Pc[1,g,n])/max(1,abs(unit_Pc[1,g,n])) )**2 )
+
+  # Voltage for busses, ratio between Vmin and Vmax
+  + (if objective_choice==1 then penalty_voltage_target_high else penalty_voltage_target_low)
+  * target_voltage_ratio
+
+  # Voltage target : value V0 in input data
+  + (if objective_choice==2 then penalty_voltage_target_high else penalty_voltage_target_low)
+  * target_voltage_data
+
+  # Reactive power of units
+  + penalty_units_reactive * sum{(g,n) in UNITON} (Q[g,n]/max(1,abs(corrected_unit_Qmin[g,n]),abs(corrected_unit_Qmax[g,n])))**2
+
+  # Ratio of transformers
+  + penalty_transfo_ratio * sum{(qq,m,n) in BRANCHCC_REGL_VAR} (branch_Ror[qq,m,n]-branch_Ror_var[qq,m,n])**2
+  ;
+  #
diff --git a/open-reac/src/main/resources/openreac/ccomputation.mod b/open-reac/src/main/resources/openreac/ccomputation.mod
@@ -0,0 +1,40 @@
+###############################################################################
+#
+# Copyright (c) 2022 2023, 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/.
+#
+###############################################################################
+
+###############################################################################
+# Reactive OPF
+# Author:  Jean Maeght 2022 2023
+###############################################################################
+
+# Definition of optimization problem
+set PROBLEM_CCOMP default { };
+
+
+###############################################################################
+# Variables
+###############################################################################
+
+var teta_ccomputation{BUS2} >=0, <=1;
+
+
+###############################################################################
+# Constraints
+###############################################################################
+
+subject to ctr_null_phase_bus_cccomputation{PROBLEM_CCOMP}: teta_ccomputation[null_phase_bus] = 0;
+
+subject to ctr_flow_cccomputation{PROBLEM_CCOMP, (qq,m,n) in BRANCH2}: teta_ccomputation[m]-teta_ccomputation[n]=0;
+# All busses AC-connected to null_phase_bus will have '0' as optimal value, other will have '1'
+
+
+###############################################################################
+# Objective function 
+###############################################################################
+
+maximize cccomputation_objective: sum{n in BUS2} teta_ccomputation[n];
diff --git a/open-reac/src/main/resources/openreac/dcopf.mod b/open-reac/src/main/resources/openreac/dcopf.mod
@@ -0,0 +1,71 @@
+###############################################################################
+#
+# Copyright (c) 2022 2023, 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/.
+#
+###############################################################################
+
+###############################################################################
+# Reactive OPF
+# Author:  Jean Maeght 2022 2023
+###############################################################################
+
+# Definition of optimization problem
+set PROBLEM_DCOPF default { };
+
+
+###############################################################################
+# Voltage variables and constraints 
+###############################################################################
+
+# Phase of voltage
+var teta_dc{n in BUSCC} 
+  <= teta_max, 
+  >= teta_min;
+
+subject to ctr_null_phase_bus_dc{PROBLEM_DCOPF}: teta_dc[null_phase_bus] = 0;
+
+
+###############################################################################
+# Active power variables and constraints 
+###############################################################################
+
+# Variable flow is the flow from bus 1 to bus 2
+var activeflow{BRANCHCC};
+subject to ctr_activeflow{PROBLEM_DCOPF, (qq,m,n) in BRANCHCC}:
+  activeflow[qq,m,n] = base100MVA * (teta_dc[m]-teta_dc[n]) / branch_X_mod[qq,m,n];#* branch_X_mod[qq,m,n] / (branch_X_mod[qq,m,n]**2+branch_R[1,qq,m,n]**2);
+
+# Generation for DCOPF
+var P_dcopf{(g,n) in UNITON}; # >= unit_Pmin[1,g,n], <= unit_Pmax[1,g,n];
+
+# Slack variable for each bus
+# >=0 if too much generation in bus, <=0 if missing generation
+var balance_pos{BUSCC} >= 0;
+var balance_neg{BUSCC} >= 0;
+
+# Active power balance at each bus
+subject to ctr_balance{PROBLEM_DCOPF, n in BUSCC}:
+  - sum{(g,n) in UNITON} P_dcopf[g,n]
+  - sum{(b,n) in BATTERYCC} battery_p0[1,b,n]
+  + sum{(c,n) in LOADCC} load_PFix[1,c,n]
+  + sum{(qq,n,m) in BRANCHCC} activeflow[qq,n,m] # active power flow outgoing on branch qq at bus n
+  - sum{(qq,m,n) in BRANCHCC} activeflow[qq,m,n] # active power flow entering in bus n on branch qq
+  + sum{(vscconv,n) in VSCCONVON} vscconv_P0[1,vscconv,n]
+  + sum{(l,n) in LCCCONVON} lccconv_P0[1,l,n]
+  =
+  balance_pos[n] - balance_neg[n];
+
+
+###############################################################################
+# Objective function and penalties
+###############################################################################
+
+param penalty_gen     := 1;
+param penalty_balance := 1000;
+
+minimize problem_dcopf_objective:
+    penalty_gen     * sum{(g,n) in UNITON} ((P_dcopf[g,n]-unit_Pc[1,g,n])/max(0.01*abs(unit_Pc[1,g,n]),1))**2
+  + penalty_balance * sum{n in BUSCC} ( balance_pos[n] + balance_neg[n] )
+  ;