wip

Signed-off-by: p-arvy <pierre.arvy@artelys.com>
powsybl · Nov 27, 2024 · 90a84c8 · 90a84c8
1 parent 7bdbea1
commit 90a84c8
Show file tree

Hide file tree

Showing 9 changed files with 267 additions and 77 deletions.
diff --git a/open-reac/README.md b/open-reac/README.md
@@ -407,3 +407,8 @@ of the DCOPF fails (see [6](#6-direct-current-optimal-power-flow)),
 the problem is considered as inconsistent.
 Then, the script `reactiveopfexit.run` is executed and the file `reactiveopf_results_indic.txt` described in [8.1](#81-in-case-of-convergence)
  is exported, without the information on the calculated angles.
+
+TODO 
+Should ratio tap changers on branch with one side opened be optimized?
+Clarify if they should be taken into account on branches with side 1 opened
+Surement une erreur dans le code AMPL. QUand G1 est non nul alors il y a des problèmes
diff --git a/open-reac/src/main/resources/openreac/acopf.mod b/open-reac/src/main/resources/openreac/acopf.mod
@@ -132,48 +132,34 @@ var Red_Tran_Rea_Inv{(qq,m,n) in BRANCHCC } =
   + V[n] * (branch_admi[qq,m,n]*cos(branch_angper[qq,m,n])-branch_Bex_mod[qq,m,n])
   ;
 
-# Penalized active/reactive power on branches with one side opened
-var act_power_bus2_opened{(qq,m,n) in BRANCH_WITH_SIDE_2_OPENED} =
-(branch_Ror[qq,m,n])**2 * V[m] *
-(branch_Gor_mod[qq,m,n] + (branch_admi[qq,m,n])**2 * branch_Gex_mod[qq,m,n]
-/ ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-+ ((branch_Bex_mod[qq,m,n])**2 + (branch_Gex_mod[qq,m,n])**2) * branch_admi[qq,m,n] * sin(branch_angper[qq,m,n])
-/ ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-);
-
-var rea_power_bus2_opened{(qq,m,n) in BRANCH_WITH_SIDE_2_OPENED} =
-- (branch_Ror[qq,m,n])**2 * V[m] *
-(branch_Bor_mod[qq,m,n] + (branch_admi[qq,m,n])**2 * branch_Bex_mod[qq,m,n]
-/ ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-- ((branch_Bex_mod[qq,m,n])**2 + (branch_Gex_mod[qq,m,n])**2) * branch_admi[qq,m,n] * cos(branch_angper[qq,m,n])
-/ ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (-branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-);
-
-var act_power_bus1_opened{(qq,m,n) in BRANCH_WITH_SIDE_1_OPENED} =
-V[n] *
-(branch_Gex_mod[qq,m,n]
-+ (branch_admi[qq,m,n])**2 * branch_Gor_mod[qq,m,n]
-/ ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-+ ((branch_Bor_mod[qq,m,n])**2 + (branch_Gor_mod[qq,m,n])**2) * branch_admi[qq,m,n] * sin(branch_angper[qq,m,n])
-/ ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-);
-
-var rea_power_bus1_opened{(qq,m,n) in BRANCH_WITH_SIDE_1_OPENED} =
-- V[n] *
-(branch_Bex_mod[qq,m,n]
-+ (branch_admi[qq,m,n])**2 * branch_Bor_mod[qq,m,n]
-/ ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-- ((branch_Bor_mod[qq,m,n])**2 + (branch_Gor_mod[qq,m,n])**2) * branch_admi[qq,m,n] * cos(branch_angper[qq,m,n])
-/ ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
-+ (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) # Shunt
-);
+var Red_Tran_Act_Dir_Side_2_Opened{(qq,m,n) in BRANCHCC_WITH_SIDE_2_OPENED} =
+    (branch_Ror[qq,m,n])**2 * V[m] * (branch_Gor_mod[qq,m,n] + (branch_admi[qq,m,n])**2 * branch_Gex_mod[qq,m,n] / ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ) 
+    + ((branch_Bex_mod[qq,m,n])**2 + (branch_Gex_mod[qq,m,n])**2) * branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]) / ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ))
+  ;
+
+var Red_Tran_Rea_Dir_Side_2_Opened{(qq,m,n) in BRANCHCC_WITH_SIDE_2_OPENED} =
+    - (branch_Ror[qq,m,n])**2 * V[m] * (branch_Bor_mod[qq,m,n] + (branch_admi[qq,m,n])**2 * branch_Bex_mod[qq,m,n] / ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 )
+    - ((branch_Bex_mod[qq,m,n])**2 + (branch_Gex_mod[qq,m,n])**2) * branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]) / ( (branch_Gex_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (-branch_Bex_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ))
+  ;
+
+var Red_Tran_Act_Inv_Side_1_Opened{(qq,m,n) in BRANCHCC_WITH_SIDE_1_OPENED} =
+    V[n] * (branch_Gex_mod[qq,m,n] + (branch_admi[qq,m,n])**2 * branch_Gor_mod[qq,m,n] / ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 )
+    + ((branch_Bor_mod[qq,m,n])**2 + (branch_Gor_mod[qq,m,n])**2) * branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]) / ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ))
+  ;
+
+var Red_Tran_Rea_Inv_Side_1_Opened{(qq,m,n) in BRANCHCC_WITH_SIDE_1_OPENED} =
+    - V[n] * (branch_Bex_mod[qq,m,n]
+    + (branch_admi[qq,m,n])**2 * branch_Bor_mod[qq,m,n] / ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 )
+    - ((branch_Bor_mod[qq,m,n])**2 + (branch_Gor_mod[qq,m,n])**2) * branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]) / ( (branch_Gor_mod[qq,m,n] + branch_admi[qq,m,n] * sin(branch_angper[qq,m,n]))**2
+    + (- branch_Bor_mod[qq,m,n] + branch_admi[qq,m,n] * cos(branch_angper[qq,m,n]))**2 ))
+  ;
 
 
 #
@@ -185,8 +171,8 @@ subject to ctr_balance_P{PROBLEM_ACOPF,k in BUSCC}:
 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]
 # Flows on branches with one side opened
-+ sum{(qq,k,n) in BRANCH_WITH_SIDE_2_OPENED} base100MVA * V[k] * act_power_bus2_opened[qq,k,n]
-+ sum{(qq,m,k) in BRANCH_WITH_SIDE_1_OPENED} base100MVA * V[k] * act_power_bus1_opened[qq,m,k]
++ sum{(qq,k,n) in BRANCHCC_WITH_SIDE_2_OPENED} base100MVA * V[k] * Red_Tran_Act_Dir_Side_2_Opened[qq,k,n]
++ sum{(qq,m,k) in BRANCHCC_WITH_SIDE_1_OPENED} base100MVA * V[k] * Red_Tran_Act_Inv_Side_1_Opened[qq,m,k]
 # Generating units
 - sum{(g,k) in UNITON} P[g,k]
 # Batteries
@@ -217,8 +203,8 @@ subject to ctr_balance_Q{PROBLEM_ACOPF,k in BUSCC}:
 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]
 # Flows on branches with one side opened
-+ sum{(qq,k,n) in BRANCH_WITH_SIDE_2_OPENED} base100MVA * V[k] * rea_power_bus2_opened[qq,k,n]
-+ sum{(qq,m,k) in BRANCH_WITH_SIDE_1_OPENED} base100MVA * V[k] * rea_power_bus1_opened[qq,m,k]
++ sum{(qq,k,n) in BRANCHCC_WITH_SIDE_2_OPENED} base100MVA * V[k] * Red_Tran_Rea_Dir_Side_2_Opened[qq,k,n]
++ sum{(qq,m,k) in BRANCHCC_WITH_SIDE_1_OPENED} base100MVA * V[k] * Red_Tran_Rea_Inv_Side_1_Opened[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]

diff --git a/open-reac/src/main/resources/openreac/commons.mod b/open-reac/src/main/resources/openreac/commons.mod
@@ -47,10 +47,13 @@ set BUS2:= setof {(1,n) in BUS:
 set BRANCH2:= setof {(1,qq,m,n) in BRANCH: m in BUS2 and n in BUS2} (qq,m,n);
 
 set BUSCC dimen 1 default {};
+# Branches with bus on side 1 and 2 in CC
 set BRANCHCC  := {(qq,m,n) in BRANCH2: m in BUSCC and n in BUSCC};
-set BRANCH_WITH_SIDE_2_OPENED := setof {(1,qq,m,n) in BRANCH: m in BUSCC and n == -1 and m != n} (qq,m,n);
-set BRANCH_WITH_SIDE_1_OPENED := setof {(1,qq,m,n) in BRANCH: m == -1 and n in BUSCC and m != n} (qq,m,n);
-set ALL_BRANCH_TO_CONSIDER:= BRANCHCC union BRANCH_WITH_SIDE_2_OPENED union BRANCH_WITH_SIDE_1_OPENED;
+# Branches with disconnected bus on side 2
+set BRANCHCC_WITH_SIDE_2_OPENED := setof {(1,qq,m,n) in BRANCH: m in BUSCC and n == -1 and m != n} (qq,m,n);
+# Branches with disconnected bus on side 1
+set BRANCHCC_WITH_SIDE_1_OPENED := setof {(1,qq,m,n) in BRANCH: m == -1 and n in BUSCC and m != n} (qq,m,n);
+set ALL_BRANCHCC := BRANCHCC union BRANCHCC_WITH_SIDE_2_OPENED union BRANCHCC_WITH_SIDE_1_OPENED;
 
 
 ###############################################################################
@@ -97,29 +100,29 @@ set LCCCONVON := setof{(t,l,n) in LCCCONV:
 
 # Branches with zero or near zero impedances
 # Notice: module of Z is equal to square root of (R^2+X^2)
-set BRANCHZNULL := {(qq,m,n) in ALL_BRANCH_TO_CONSIDER: branch_R[1,qq,m,n]^2+branch_X[1,qq,m,n]^2 <= Znull^2};
+set BRANCHZNULL := {(qq,m,n) in ALL_BRANCHCC: branch_R[1,qq,m,n]^2+branch_X[1,qq,m,n]^2 <= Znull^2};
 
 
 # If in BRANCHZNULL, then set X to ZNULL
-param branch_X_mod{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} :=
+param branch_X_mod{(qq,m,n) in ALL_BRANCHCC} :=
   if (qq,m,n) in BRANCHZNULL then Znull
   else branch_X[1,qq,m,n];
-check {(qq,m,n) in ALL_BRANCH_TO_CONSIDER}: abs(branch_X_mod[qq,m,n]) > 0;
+check {(qq,m,n) in ALL_BRANCHCC}: abs(branch_X_mod[qq,m,n]) > 0;
 
 # If in BRANCHZNULL, then set Gor/Gex/Bor/Bex to 0
-param branch_Gor_mod{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} :=
+param branch_Gor_mod{(qq,m,n) in ALL_BRANCHCC} :=
     if (qq,m,n) in BRANCHCC and (qq,m,n) in BRANCHZNULL then 0
     else branch_Gor[1,qq,m,n];
 
-param branch_Gex_mod{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} :=
+param branch_Gex_mod{(qq,m,n) in ALL_BRANCHCC} :=
     if (qq,m,n) in BRANCHCC and (qq,m,n) in BRANCHZNULL then 0
     else branch_Gex[1,qq,m,n];
 
-param branch_Bor_mod{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} :=
+param branch_Bor_mod{(qq,m,n) in ALL_BRANCHCC} :=
     if (qq,m,n) in BRANCHCC and (qq,m,n) in BRANCHZNULL then 0
     else branch_Bor[1,qq,m,n];
 
-param branch_Bex_mod{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} :=
+param branch_Bex_mod{(qq,m,n) in ALL_BRANCHCC} :=
     if (qq,m,n) in BRANCHCC and (qq,m,n) in BRANCHZNULL then 0
     else branch_Bex[1,qq,m,n];
 
@@ -187,18 +190,18 @@ param branch_Rdeph{(qq,m,n) in BRANCHCC_DEPH} =
   else branch_R[1,qq,m,n]
   ;
 
-param branch_angper{(qq,m,n) in ALL_BRANCH_TO_CONSIDER} =
+param branch_angper{(qq,m,n) in ALL_BRANCHCC} =
   if (qq,m,n) in BRANCHCC_DEPH
   then atan2(branch_Rdeph[qq,m,n], branch_Xdeph[qq,m,n])
   else atan2(branch_R[1,qq,m,n]  , branch_X_mod[qq,m,n]  );
 
-param branch_admi {(qq,m,n) in ALL_BRANCH_TO_CONSIDER} =
+param branch_admi {(qq,m,n) in ALL_BRANCHCC} =
   if (qq,m,n) in BRANCHCC_DEPH
   then 1./sqrt(branch_Rdeph[qq,m,n]^2 + branch_Xdeph[qq,m,n]^2 )
   else 1./sqrt(branch_R[1,qq,m,n]^2   + branch_X_mod[qq,m,n]^2   );
 
 # Later in this file, a variable branch_Ror_var will be created, to replace branch_Ror when it is not variable
-param branch_Ror {(qq,m,n) in ALL_BRANCH_TO_CONSIDER} =
+param branch_Ror {(qq,m,n) in ALL_BRANCHCC} =
     ( if ((qq,m,n) in BRANCHCC_REGL)
       then tap_ratio[1,regl_table[1,branch_ptrRegl[1,qq,m,n]],regl_tap0[1,branch_ptrRegl[1,qq,m,n]]]
       else 1.0
@@ -208,13 +211,13 @@ param branch_Ror {(qq,m,n) in ALL_BRANCH_TO_CONSIDER} =
       else 1.0
     )
   * (branch_cstratio[1,qq,m,n]);
-param branch_Rex {(q,m,n) in ALL_BRANCH_TO_CONSIDER} = 1; # In IIDM, everything is in bus1 so ratio at bus2 is always 1
+param branch_Rex {(q,m,n) in ALL_BRANCHCC} = 1; # In IIDM, everything is in bus1 so ratio at bus2 is always 1
 
-param branch_dephor {(qq,m,n) in ALL_BRANCH_TO_CONSIDER} =
+param branch_dephor {(qq,m,n) in ALL_BRANCHCC} =
   if ((qq,m,n) in BRANCHCC_DEPH)
   then tap_angle [1,deph_table[1,branch_ptrDeph[1,qq,m,n]],deph_tap0[1,branch_ptrDeph[1,qq,m,n]]]
   else 0;
-param branch_dephex {(qq,m,n) in ALL_BRANCH_TO_CONSIDER} = 0; # In IIDM, everything is in bus1 so dephase at bus2 is always 0
+param branch_dephex {(qq,m,n) in ALL_BRANCHCC} = 0; # In IIDM, everything is in bus1 so dephase at bus2 is always 0
 
 
 ###############################################################################

diff --git a/open-reac/src/main/resources/openreac/reactiveopfoutput.run b/open-reac/src/main/resources/openreac/reactiveopfoutput.run
@@ -286,7 +286,9 @@ printf "%s %i\n","nb_bus_with_voltage_value",card(BUSVV) > (fileOut);
 printf "%s %i\n","nb_bus_with_reactive_slacks",card(BUSCC_SLACK) > (fileOut);
 printf "%s %i\n","nb_bus_without_reactive_slacks",card(BUSCC diff BUSCC_SLACK) > (fileOut);
 printf "%s %i\n","nb_branch_in_data_file",card(BRANCH) > (fileOut);
-printf "%s %i\n","nb_branch_in_AC_CC",card(BRANCHCC) > (fileOut);
+printf "%s %i\n","nb_branch_in_AC_CC",card(ALL_BRANCHCC) > (fileOut);
+printf "%s %i\n","nb_branch_in_AC_CC_side_1_opened",card(BRANCHCC_WITH_SIDE_1_OPENED) > (fileOut);
+printf "%s %i\n","nb_branch_in_AC_CC_side_2_opened",card(BRANCHCC_WITH_SIDE_2_OPENED) > (fileOut);
 printf "%s %i\n","nb_branch_with_nonsmall_impedance",card(BRANCHCC diff BRANCHZNULL) > (fileOut);
 printf "%s %i\n","nb_branch_with_zero_or_small_impedance",card(BRANCHZNULL) > (fileOut);
 printf "%s %i\n","nb_unit_in_data_file",card(UNIT) > (fileOut);
@@ -376,3 +378,37 @@ for {(g,n) in UNITCC diff UNITON}
     '"' & bus_id[1,n] & '"'
     > (fileOut);
 close (fileOut);
+
+# write flows on the branches
+let fileOut := "reactiveopf_results_branches.csv";
+printf{LOG_DEBUG} "#qq;m;n;P1;Q1;P2;Q2;branch_id;\n" > (fileOut);
+for {(qq,m,n) in ALL_BRANCHCC}
+  printf{LOG_DEBUG} "%i;%i;%i;%.3f;%.3f;%.3f;%.3f;%s;\n",
+    qq,m,n,
+    # p1 flow
+    if (qq,m,n) in BRANCHCC 
+    then base100MVA * V[m] * Red_Tran_Act_Dir[qq,m,n]
+    else if (qq,m,n) in BRANCHCC_WITH_SIDE_2_OPENED
+    then base100MVA * V[m] * Red_Tran_Act_Dir_Side_2_Opened[qq,m,n]
+    else 0, # no flow as branch is opened on side 1
+    # q1 flow
+    if (qq,m,n) in BRANCHCC 
+    then base100MVA * V[m] * Red_Tran_Rea_Dir[qq,m,n]
+    else if (qq,m,n) in BRANCHCC_WITH_SIDE_2_OPENED
+    then base100MVA * V[m] * Red_Tran_Rea_Dir_Side_2_Opened[qq,m,n]
+    else 0, # no flow as branch is opened on side 1
+    # p2 flow
+    if (qq,m,n) in BRANCHCC 
+    then base100MVA * V[n] * Red_Tran_Act_Inv[qq,m,n]
+    else if (qq,m,n) in BRANCHCC_WITH_SIDE_1_OPENED
+    then base100MVA * V[n] * Red_Tran_Act_Inv_Side_1_Opened[qq,m,n]
+    else 0, # no flow as branch is opened on side 2
+    # q2 flow
+    if (qq,m,n) in BRANCHCC 
+    then base100MVA * V[n] * Red_Tran_Rea_Inv[qq,m,n]
+    else if (qq,m,n) in BRANCHCC_WITH_SIDE_1_OPENED
+    then base100MVA * V[n] * Red_Tran_Rea_Inv_Side_1_Opened[qq,m,n]
+    else 0, # no flow as branch is opened on side 2
+    '"' & branch_id[1,qq,m,n] & '"'
+    > (fileOut);
+close (fileOut);
diff --git a/open-reac/src/test/java/com/powsybl/openreac/OpenReacBranchOneSideOpenTest.java b/open-reac/src/test/java/com/powsybl/openreac/OpenReacBranchOneSideOpenTest.java
@@ -0,0 +1,6 @@
+package com.powsybl.openreac;
+
+import org.junit.jupiter.api.Test;
+
+public class OpenReacBranchOneSideOpenTest {
+}