Fix ACOPF rectangular formulation and update datasets (#153)

* Fix rectangular formulation
* Add 'ACLocal' option for OPF
* Update pglib datasets

docs/source/mods/opf/opf.rst

@@ -224,7 +224,7 @@ argument specifies otherwise.
     :options: +NORMALIZE_WHITESPACE +ELLIPSIS
 
     ...
-    Optimize a model with 73 rows, 107 columns and 208 nonzeros
+    Optimize a model with 19 rows, 107 columns and 28 nonzeros
     ...
     Optimal solution found...
     ...
@@ -276,7 +276,7 @@ solution information, as specified below.
             :options: +NORMALIZE_WHITESPACE
 
             >>> result['bus'][0]
-            {... 'Vm': 1.09..., 'Va': 0, ...}
+            {... 'Vm': 1.09..., 'Va': 0.0, ...}
 
     .. tab:: Branches
 
@@ -410,7 +410,7 @@ whether branch switching allows a better solution.
     :options: +NORMALIZE_WHITESPACE +ELLIPSIS
 
     ...
-    Optimize a model with 278 rows, 185 columns and 694 nonzeros
+     Optimize a model with 212 rows, 185 columns and 424 nonzeros
     ...
 
 Plotting the resulting solution shows that one branch has been turned off in the
@@ -483,9 +483,9 @@ data:
 .. doctest:: opf
 
     >>> print(violations['branch'][6]['limitviol'])
-    66.33435016796234
+    66.33435...
     >>> print(violations['bus'][3]['Pviol'])
-    -318.8997836192236
+    -318.8997...
 
 In this case, the limit at branch 6 and the real power injection at bus 3 are
 violated by the given input voltages.

src/gurobi_optimods/opf/api.py

@@ -74,8 +74,18 @@ def solve_opf(
         fields
     """
 
-    # Exact cartesian AC (force jabr for performance reasons)
-    if opftype.lower() == "ac":
+    # use aclocal to run Gurobi as a local solver (stops at the first feasible solution)
+    if opftype.lower() == "aclocal":
+        opftype = "ac"
+        useef = True
+        usejabr = False
+        default_solver_params = {
+            "Presolve": 0,
+            "SolutionLimit": 1,
+            "NodeLimit": 0,
+            "GURO_PAR_NLBARSLOPPYLIMIT": 2000,
+        }
+    elif opftype.lower() == "ac":
         opftype = "ac"
         useef = True
         usejabr = True

src/gurobi_optimods/opf/converters.py

@@ -78,6 +78,10 @@ def build_internal_settings(
     else:
         raise ValueError(f"Unknown OPF type {opftype}")
 
+    if opftype == "dc":
+        settings["use_ef"] = False
+        settings["skipjabr"] = True
+
     # Sub-type for IV models
     ivtype = ivtype.lower()
     if ivtype in ["plain", "aggressive"]:
@@ -174,19 +178,38 @@ def convert_case_to_internal_format(case_dict):
     if any(gen["bus"] not in bus_ids for gen in case_dict["gen"]):
         raise ValueError("Unknown bus ID referenced in generator bus")
 
+    # For each field we create a key value for use internally.
+    # Note: index != nodeID (bus_i) for buses.
+    # Note: some parts of the code still rely on these indexes being 1..n and
+    # the keys being in ascending order in dictionary iterators. Be very careful
+    # when trying to change this.
+    # Remove isolated buses and their connected branches
+    remove_buses = {bus["bus_i"] for bus in case_dict["bus"] if bus["type"] == 4}
+    if remove_buses:
+        logger.info(f"Removing buses {remove_buses} (bustype=4)")
+
     # For each field we create a key value for use internally.
     # Note: index != nodeID (bus_i) for buses.
     # Note: some parts of the code still rely on these indexes being 1..n and
     # the keys being in ascending order in dictionary iterators. Be very careful
     # when trying to change this.
     case_dict = {
         "baseMVA": case_dict["baseMVA"],
-        "bus": {i + 1: dict(bus) for i, bus in enumerate(case_dict["bus"])},
-        "branch": {i + 1: dict(branch) for i, branch in enumerate(case_dict["branch"])},
+        "bus": {
+            i + 1: dict(bus)
+            for i, bus in enumerate(case_dict["bus"])
+            if bus["bus_i"] not in remove_buses
+        },
+        "branch": {
+            i + 1: dict(branch)
+            for i, branch in enumerate(case_dict["branch"])
+            if branch["fbus"] not in remove_buses and branch["tbus"] not in remove_buses
+        },
         "gen": {i + 1: dict(gen) for i, gen in enumerate(case_dict["gen"])},
         "gencost": {
             i + 1: dict(gencost) for i, gencost in enumerate(case_dict["gencost"])
         },
+        "casename": case_dict["casename"],
     }
 
     # Manual correction to ratios
@@ -195,7 +218,7 @@ def convert_case_to_internal_format(case_dict):
             branch["ratio"] = 1.0
 
     alldata = {"LP": {}, "MIP": {}}
-
+    alldata["casename"] = case_dict["casename"]
     logger.info("Building case data structures from dictionary.")
     baseMVA = alldata["baseMVA"] = case_dict["baseMVA"]
     # Buses

src/gurobi_optimods/opf/grbformulator.py

@@ -98,7 +98,17 @@ def lpformulator_optimize(alldata, model, opftype):
         for j in range(1, 1 + numbranches):
             branch = branches[j]
             zvar[branch].Start = 1.0
-
+    # Initialize evar to 1 for barrier to converge
+    if alldata["use_ef"]:
+        evar = alldata["LP"]["evar"]
+        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]
+            evar[bus].PStart = 1.0
+        model.update()
     model.optimize()
 
     # Check model status and re-optimize if numerical trouble or inconclusive results.
@@ -161,6 +171,7 @@ def turn_solution_into_result_dict(alldata, model, opftype, type):
     # Reconstruct case data from our data
     result = {}
     result["baseMVA"] = alldata["baseMVA"]
+    result["casename"] = alldata["casename"]
     baseMVA = result["baseMVA"]
     result["bus"] = {}
     result["gen"] = {}
@@ -347,14 +358,13 @@ def fill_result_fields(alldata, model, opftype, result):
                 databus = alldata["buses"][busindex]
                 # Override old values
                 # Voltage magnitude is root of cvar because cvar = square of voltage magnitude given as e^2 + f^2
-                if alldata["doiv"]:  # doiv makes sure that e, f variables are present
-                    resbus["Vm"] = math.sqrt(
-                        alldata["LP"]["evar"][databus].X ** 2
-                        + alldata["LP"]["fvar"][databus].X ** 2
-                    )
-                else:
-                    resbus["Vm"] = math.sqrt(alldata["LP"]["cvar"][databus].X)
-
+                resbus["Vm"] = math.sqrt(
+                    alldata["LP"]["evar"][databus].X ** 2
+                    + alldata["LP"]["fvar"][databus].X ** 2
+                )
+                resbus["Va"] = math.atan(
+                    alldata["LP"]["fvar"][databus].X / alldata["LP"]["evar"][databus].X
+                )
             if alldata["doiv"]:
                 # Need to fill cvar[branch] dictionary to compute angles for IV
                 # Note that there is no cvar dictionary for IV!!!
@@ -369,8 +379,6 @@ def fill_result_fields(alldata, model, opftype, result):
 
                 alldata["LP"]["cvar"] = cvar
 
-            compute_voltage_angles(alldata, result)
-
         if alldata["dopolar"] and not alldata["doiv"]:
             for busindex in result["bus"]:
                 resbus = result["bus"][busindex]