Fixing the ACOPF rectangular formulation and updating datasets

src/gurobi_optimods/opf/api.py

@@ -78,7 +78,7 @@ def solve_opf(
     if opftype.lower() == "ac":
         opftype = "ac"
         useef = True
-        usejabr = True
+        usejabr = False
         default_solver_params = {"MIPGap": 1e-3, "OptimalityTol": 1e-3}
     # AC relaxation using the JABR inequality
     elif opftype.lower() == "acrelax":

src/gurobi_optimods/opf/converters.py

@@ -174,19 +174,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 +214,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

@@ -1,5 +1,6 @@
 import logging
 import math
+import sys
 import time
 from enum import Enum
 
@@ -98,7 +99,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.
@@ -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,7 +379,7 @@ def fill_result_fields(alldata, model, opftype, result):
 
                 alldata["LP"]["cvar"] = cvar
 
-            compute_voltage_angles(alldata, result)
+            # compute_voltage_angles(alldata, result)
 
         if alldata["dopolar"] and not alldata["doiv"]:
             for busindex in result["bus"]:

src/gurobi_optimods/opf/grbformulator_ac.py

@@ -100,9 +100,6 @@ def lpformulator_ac_create_vars(alldata, model):
 
             lower = gen.Qmin * gen.status
             upper = gen.Qmax * gen.status
-            if bus.nodetype == 3:
-                lower = -GRB.INFINITY
-                upper = GRB.INFINITY
             GenQvar[gen] = model.addVar(
                 lb=lower, ub=upper, name=f"GQ_{gen.count}_{gen.nodeID}"
             )
@@ -476,6 +473,8 @@ def lpformulator_ac_create_constraints(alldata, model):
     IDtoCountmap = alldata["IDtoCountmap"]
     cvar = alldata["LP"]["cvar"]
     svar = alldata["LP"]["svar"]
+    evar = alldata["LP"]["evar"]
+    fvar = alldata["LP"]["fvar"]
     Pvar_f = alldata["LP"]["Pvar_f"]
     Pvar_t = alldata["LP"]["Pvar_t"]
     Qvar_f = alldata["LP"]["Qvar_f"]
@@ -517,15 +516,16 @@ def lpformulator_ac_create_constraints(alldata, model):
             continue
 
         # Gff cff + Gft cft + Bft sft
-        expr = gp.LinExpr(
-            [branch.Gff, branch.Gft, branch.Bft],
-            [cvar[busf], cvar[branch], svar[branch]],
-        )
         expP = Pvar_f[branch]
         if alldata["branchswitching_mip"]:
             expP += twinPvar_f[branch]
-        branch.Pdeffconstr = model.addConstr(expr == expP, name=branch.Pfcname)
-
+        branch.Pdeffconstr = model.addConstr(
+            branch.Gff * (evar[busf] * evar[busf] + fvar[busf] * fvar[busf])
+            + branch.Gft * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust])
+            + branch.Bft * (-evar[busf] * fvar[bust] + fvar[busf] * evar[bust])
+            == expP,
+            name=branch.Pfcname,
+        )
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
                 coeff = branch.limit
@@ -554,20 +554,16 @@ def lpformulator_ac_create_constraints(alldata, model):
             count += 4
 
         # Gtt ctt + Gtf cft + Btf stf = Gtt ctt + Gtf cft - Btf sft
-        expr = gp.LinExpr(
-            [
-                branch.Gtt,
-                branch.Gtf,
-                -branch.Btf,
-            ],  # minus because svarft = -svartf
-            [cvar[bust], cvar[branch], svar[branch]],
-        )
-
         expP = Pvar_t[branch]
         if alldata["branchswitching_mip"]:
             expP += twinPvar_t[branch]
-
-        branch.Pdeftconstr = model.addConstr(expr == expP, name=branch.Ptcname)
+        branch.Pdeftconstr = model.addConstr(
+            branch.Gtt * (evar[bust] * evar[bust] + fvar[bust] * fvar[bust])
+            + branch.Gtf * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust])
+            - branch.Btf * (-evar[busf] * fvar[bust] + fvar[busf] * evar[bust])
+            == expP,
+            name=branch.Ptcname,
+        )
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -623,16 +619,16 @@ def lpformulator_ac_create_constraints(alldata, model):
             continue
 
         # -Bff cff - Bft cft + Gft sft
-        expr = gp.LinExpr(
-            [-branch.Bff, -branch.Bft, branch.Gft],
-            [cvar[busf], cvar[branch], svar[branch]],
-        )
-
         expQ = Qvar_f[branch]
         if alldata["branchswitching_mip"]:
             expQ += twinQvar_f[branch]
-        branch.Qdeffconstr = model.addConstr(expr == expQ, name=branch.Qfcname)
-
+        branch.Qdeffconstr = model.addConstr(
+            -branch.Bff * (evar[busf] * evar[busf] + fvar[busf] * fvar[busf])
+            - branch.Bft * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust])
+            + branch.Gft * (-evar[busf] * fvar[bust] + fvar[busf] * evar[bust])
+            == expQ,
+            name=branch.Qfcname,
+        )
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
                 coeff = branch.limit
@@ -661,15 +657,16 @@ def lpformulator_ac_create_constraints(alldata, model):
             count += 4
 
         # -Btt ctt - Btf cft + Gtf stf = -Btt ctt - Btf cft - Gtf sft
-        expr = gp.LinExpr(
-            [-branch.Btt, -branch.Btf, -branch.Gtf],  # again, same minus
-            [cvar[bust], cvar[branch], svar[branch]],
-        )
-
         expQ = Qvar_t[branch]
         if alldata["branchswitching_mip"]:
             expQ += twinQvar_t[branch]
-        branch.Qdeftconstr = model.addConstr(expr == expQ, name=branch.Qtcname)
+        branch.Qdeftconstr = model.addConstr(
+            -branch.Btt * (evar[bust] * evar[bust] + fvar[bust] * fvar[bust])
+            - branch.Btf * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust])
+            - branch.Gtf * (-evar[busf] * fvar[bust] + fvar[busf] * evar[bust])
+            == expQ,
+            name=branch.Qtcname,
+        )
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -708,11 +705,11 @@ def lpformulator_ac_create_constraints(alldata, model):
     )
     count = 0
     for j, bus in buses.items():
-        expr = gp.LinExpr()
+        expr = gp.QuadExpr()
         qexpr = gp.QuadExpr()
         doquad = False
         if bus.Gs != 0:
-            expr.add(bus.Gs * cvar[bus])
+            expr.add(bus.Gs * (evar[bus] * evar[bus] + fvar[bus] * fvar[bus]))
 
         if alldata["branchswitching_comp"] is False:
             for branchid in bus.frombranchids.values():
@@ -739,10 +736,10 @@ def lpformulator_ac_create_constraints(alldata, model):
         count += 1
 
     for j, bus in buses.items():
-        expr = gp.LinExpr()
+        expr = gp.QuadExpr()
 
         if bus.Bs != 0:
-            expr.add(-bus.Bs * cvar[bus])
+            expr.add(-bus.Bs * (evar[bus] * evar[bus] + fvar[bus] * fvar[bus]))
 
         if alldata["branchswitching_comp"] is False:
             for branchid in bus.frombranchids.values():
@@ -769,6 +766,52 @@ def lpformulator_ac_create_constraints(alldata, model):
 
     logger.info(f"    {count} balance constraints added.")
 
+    if alldata["skipjabr"] is True and alldata["use_ef"]:
+        # voltage magnitude bounds
+        logger.info("  Adding voltage magnitude constraints.")
+        count = 0
+        for j, bus in buses.items():
+            expr = gp.QuadExpr()
+            expr.add(evar[bus] * evar[bus] + fvar[bus] * fvar[bus])
+            model.addConstr(
+                expr <= bus.Vmax * bus.Vmax, name="VUb%d_%d" % (j, bus.nodeID)
+            )
+            model.addConstr(
+                expr >= bus.Vmin * bus.Vmin, name="VLb%d_%d" % (j, bus.nodeID)
+            )
+            count += 2
+        logger.info(f"    {count} voltage magnitude constraints added.")
+
+        # phase angle bounds
+        logger.info("  Adding phase angle bounds constraints.")
+        count = 0
+        for j, branch in branches.items():
+            if not branch.status or branch.unboundedlimit is True:
+                continue
+
+            f = branch.f
+            t = branch.t
+            count_of_f = IDtoCountmap[f]
+            count_of_t = IDtoCountmap[t]
+            busf = buses[count_of_f]
+            bust = buses[count_of_t]
+            expr = gp.QuadExpr()
+            expr.add(fvar[busf] * evar[bust] - evar[busf] * fvar[bust])
+            model.addConstr(
+                expr
+                <= math.tan(branch.maxangle_rad)
+                * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust]),
+                name="AngleUb_%d_%d_%d" % (j, f, t),
+            )
+            model.addConstr(
+                expr
+                >= math.tan(branch.minangle_rad)
+                * (evar[busf] * evar[bust] + fvar[busf] * fvar[bust]),
+                name="AngleLb_%d_%d_%d" % (j, f, t),
+            )
+            count += 2
+        logger.info(f"    {count} phase angle bounds constraints added.")
+
     # Injection defs
     logger.info("  Adding injection definition constraints.")
     count = 0
@@ -905,7 +948,7 @@ def lpformulator_ac_create_constraints(alldata, model):
         lpformulator_ac_add_polarconstraints(alldata, model)
 
     # nonconvex e, f representation
-    if alldata["use_ef"]:
+    if alldata["skipjabr"] is False and alldata["use_ef"]:
         lpformulator_ac_add_nonconvexconstraints(alldata, model)
 
     if alldata["branchswitching_mip"] or alldata["branchswitching_comp"]:

src/gurobi_optimods/opf/io.py

@@ -5,6 +5,8 @@
 - write_case_matpower: write case (or solution) to MATPOWER format .mat file
 """
 
+import os
+
 import numpy as np
 import pandas as pd
 import scipy
@@ -153,6 +155,7 @@ def fix_shape(arr, ncols=None):
         "gen": gen.to_dict("records"),
         "branch": branch.to_dict("records"),
         "gencost": gencost,
+        "casename": os.path.splitext(file_path)[0],
     }