Enable ACOPF polar formulation (#176)

* Add ACPlocal (default), ACPglobal, ACRlocal, ACRglobal options for opftype
* Remove AC and AClocal options
* Update formulation to use nonlinear expressions

Author: hhijazi

docs/source/mods/opf/opf.rst

@@ -222,13 +222,13 @@ argument specifies otherwise.
     from gurobi_optimods import datasets
 
     case = datasets.load_opf_example("case9")
-    result = opf.solve_opf(case, opftype="AC")
+    result = opf.solve_opf(case, opftype="ACRGLOBAL")
 
 .. testoutput:: opf
     :options: +NORMALIZE_WHITESPACE +ELLIPSIS
 
     ...
-    Optimize a model with 19 rows, 107 columns and 28 nonzeros
+    Optimize a model with 18 rows, 106 columns and 24 nonzeros...
     ...
     Optimal solution found...
     ...
@@ -242,7 +242,7 @@ optimality within the time limit, the best known solution will be returned.
 
 .. testcode:: opf
 
-    result = opf.solve_opf(case, opftype="AC", time_limit=60)
+    result = opf.solve_opf(case, opftype="ACRGLOBAL", time_limit=60)
 
 .. testoutput:: opf
     :hide:
@@ -406,7 +406,7 @@ whether branch switching allows a better solution.
 
     case = datasets.load_opf_example("case9-switching")
     result = opf.solve_opf(
-        case, opftype="AC",
+        case, opftype="ACRGLOBAL",
         branch_switching=True,
         min_active_branches=0.1,
         time_limit=60,
@@ -417,7 +417,7 @@ whether branch switching allows a better solution.
     :options: +NORMALIZE_WHITESPACE +ELLIPSIS
 
     ...
-     Optimize a model with 212 rows, 185 columns and 424 nonzeros
+     Optimize a model with 211 rows, 184 columns and 420 nonzeros...
     ...
 
 Plotting the resulting solution shows that one branch has been turned off in the

src/gurobi_optimods/opf/api.py

@@ -7,6 +7,8 @@
 
 import logging
 
+import gurobipy as gp
+
 from gurobi_optimods.opf import converters, grbformulator, violations
 from gurobi_optimods.utils import optimod
 
@@ -16,7 +18,7 @@
 @optimod()
 def solve_opf(
     case,
-    opftype="AC",
+    opftype="ACPLOCAL",
     branch_switching=False,
     min_active_branches=0.9,
     use_mip_start=False,
@@ -54,7 +56,8 @@ def solve_opf(
     case : dict
         Dictionary holding case data
     opftype : str
-        Desired OPF model type. One of ``AC``, ``AClocal``, ``ACrelax``, or ``DC``
+        Desired OPF model type. One of ``ACRlocal``, ``ACPlocal``, ``ACRglobal``,
+        ``ACPglobal``, ``ACrelax``, or ``DC``. Defaults to ``ACPlocal``.
     branch_switching : bool, optional
         If set to True, enables branch switching.
     min_active_branches : float, optional
@@ -74,33 +77,68 @@ def solve_opf(
         fields
     """
 
-    # use aclocal to run Gurobi as a local solver (stops at the first feasible solution)
-    if opftype.lower() == "aclocal":
+    # use acrlocal to run Gurobi as a local solver with the rectangular formulation (QCQP)
+    if opftype.lower() == "acrlocal":
         opftype = "ac"
         useef = True
         usejabr = False
-        default_solver_params = {
-            "Presolve": 0,
-            "SolutionLimit": 1,
-            "NodeLimit": 0,
-            "GURO_PAR_NLBARSLOPPYLIMIT": 2000,
-        }
-    elif opftype.lower() == "ac":
+        polar = False
+        version = gp.gurobi.version()
+        if version >= (13, 0, 0):
+            default_solver_params = {
+                "OptimalityTarget": 1,
+            }
+        else:
+            default_solver_params = {
+                "Presolve": 0,
+                "SolutionLimit": 1,
+                "NodeLimit": 0,
+                "GURO_PAR_NLBARSLOPPYLIMIT": 2000,
+            }
+    # use acplocal to run Gurobi as a local solver with the polar formulation (trigonometric functions)
+    elif opftype.lower() == "acplocal":
+        opftype = "ac"
+        useef = False
+        usejabr = False
+        polar = True
+        version = gp.gurobi.version()
+        if version >= (13, 0, 0):
+            default_solver_params = {
+                "OptimalityTarget": 1,
+            }
+        else:
+            default_solver_params = {
+                "Presolve": 0,
+                "SolutionLimit": 1,
+                "NodeLimit": 0,
+                "GURO_PAR_NLBARSLOPPYLIMIT": 2000,
+            }
+    elif opftype.lower() == "acrglobal":
         opftype = "ac"
         useef = True
         usejabr = True
+        polar = False
+        default_solver_params = {"MIPGap": 1e-3, "OptimalityTol": 1e-3}
+    # Exact polar AC
+    elif opftype.lower() == "acpglobal":
+        opftype = "ac"
+        useef = False
+        usejabr = False
+        polar = True
         default_solver_params = {"MIPGap": 1e-3, "OptimalityTol": 1e-3}
     # AC relaxation using the JABR inequality
     elif opftype.lower() == "acrelax":
         opftype = "ac"
         useef = False
         usejabr = True
+        polar = False
         default_solver_params = {"MIPGap": 1e-3, "OptimalityTol": 1e-3}
     # DC linear approximation (ef & jabr are irrelevant)
     elif opftype.lower() == "dc":
         opftype = "dc"
         useef = False
         usejabr = False
+        polar = False
         default_solver_params = {"MIPGap": 1e-4, "OptimalityTol": 1e-4}
     else:
         raise ValueError(f"Unknown opftype '{opftype}'")
@@ -115,8 +153,8 @@ def solve_opf(
             branchswitching=branch_switching,
             usemipstart=use_mip_start,
             minactivebranches=min_active_branches,
-            polar=False,
-            ivtype="aggressive",
+            polar=polar,
+            ivtype="aggressive",  # ignored, no IV formulation used
             useactivelossineqs=False,
         )
 
@@ -223,5 +261,6 @@ def compute_violations(case, voltages, polar=False, *, create_env):
     converters.grbmap_volts_from_dict(alldata, voltages)
 
     # Compute model violations based on user input voltages
-    with create_env() as env:
-        return violations.compute_violations_from_voltages(env, alldata)
+    if not polar:
+        with create_env() as env:
+            return violations.compute_violations_from_voltages(env, alldata)

src/gurobi_optimods/opf/grbformulator.py

@@ -112,6 +112,16 @@ def lpformulator_optimize(alldata, model, opftype):
             for zvar in alldata["MIP"]["zvar"].values():
                 zvar.PStart = 1.0
         model.update()
+    elif alldata["dopolar"]:
+        vvar = alldata["LP"]["vvar"]
+        buses = alldata["buses"]
+        numbuses = alldata["numbuses"]
+        for var in model.getVars():
+            var.PStart = 0.0
+        for j in range(1, 1 + numbuses):
+            bus = buses[j]
+            vvar[bus].PStart = 1.0
+        model.update()
     model.optimize()
 
     # Check model status and re-optimize if numerical trouble or inconclusive results.