Repository: Gurobi/gurobi-optimods

Issue:
OPF problem
**Describe the bug**
The AC OPF is giving an infeasible solution for a case.

**To Reproduce**

- Run an OPF for attached case in MATPOWER and check infeasibility in power balance equations. 
- Save the case as .MAT file (It is attached as well).
- Load it in a python distribution with the command opf.read_case_matpower("file.mat") after importing optimod modules
- Run the AC OPF. 
- Check if solutions match.
- Check infeasibility the same way as before.

MATPOWER case file (.m) and Optimod case (.mat)
[IEEE69.zip](https://github.com/Gurobi/gurobi-optimods/files/13241332/IEEE69.zip)

**Expected behavior**
If the optimization problem was solved and the optimal solution was found, the solution must be feasible within tolerance. I expected the same solution as MATPOWER.

**Additional context**

With MATPOWER I run the OPF for a case (69 node Radial network, 1 reference node (node 1), 68 PQ nodes (2-69)) , whose cost function is polynomial (default in MATPOWER) with zero coefficients, except for the first order one (linear cost function). With this setup, we are minimizing indirectly the power in the slack generator, thus the losses.

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/3c68b78d-3d5d-4aec-9cd5-5796d06a6b75) 
![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/71d6a8c1-f428-4ead-8b95-a3398ab31775)


To verify feasibility, we use the obtained solution and check if the voltages are within limits and if the power balance equations are satisfied (in this case the line power limits are unconstrained). 

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/afc71c93-8767-4efb-b1bb-2f76d4923a70)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/d3acebe3-db67-4f4c-9d73-a558422ae8fb)

To check balance equations I use 2 tests:

1. First, create a zero array (size Nx1) for the active $p_{g}$ and reactive power $q_{g}$ at the slack bus and the demand ($p_{d}, q_{d}$) for each node, and assign the values from the solution obtained accordingly. Calculate the power losses $p_{loss},q_{loss}$. Then check power balance equations:

$$ \sum_{i=1}^{N} p_{g} - \sum_{i=1}^{N} p_{d} -p_{loss} =0$$

$$ \sum_{i=1}^{N} q_{g} - \sum_{i=1}^{N} q_{d} -q_{loss} =0$$

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/658eaeb9-8530-4445-b6f1-3d49001a5cfb)
![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/77d5041b-4183-468c-a489-8eb3de808463)


2.  Choose a power flow model, compute to each corresponding variable its values from the solution and calculate for each node the power balance equation. For instance, if the complex bus injection model is chosen, the complex admittance matrix $Y$, voltage $V_{i}$, generation $S_{i}^{g}$ and load $S_{i}^{d}$ vectors should be constructed.

$$ S_{i}^{g} -S_{i}^{d} = V_{i} \sum_{j=1}^{N} Y_{ij}^{\ast}V_{j}^{\ast} \quad \forall i,j \in N$$

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/d5d3198a-1fdb-469d-a5a8-2d8f44e10613)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/5d4bf136-95d5-47d4-a62c-bbb93655d788)

As observed, the MATPOWER solution is feasible, and is often found as reference value in the literature.


Applying the same procedure in python. First, the modules are imported, and the case is loaded. Then the AC OPF is excecuted:

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/83d51e43-c8c6-483c-b6ba-c022dc8f91c1)


The objective differs from MATPOWER, but the solution found is optimal.

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/65488656-3642-4f5c-81b6-d449c521d1d0)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/64bda85e-b082-4bc8-be00-b3514d88ddfd)
![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/1a87876a-f50d-4f2f-88fd-bbfd27d738f4)


This solution pass the first feasibility check:

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/cee7b4b1-9315-472e-b98f-241238f05371)
![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/63b33a7f-3bd5-4c03-8e06-5b98113ffbed)


However the second feasibility check is failed:

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/14d9d215-b153-481f-b7cc-0b21026843cc)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/6e814632-6aff-4ff7-b778-cadd2738d854)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/ef443684-4663-4ce0-a19f-16691310fb14)

Discussion:
@simonbowly: Hi @dimendozao, thanks for raising this. Could you please attach the code necessary to reproduce this result from your input files as a python script? We appreciate you sharing a working example, but transcribing code from screenshots is time consuming and error-prone.

In addition, can you please share the versions of installed packages required to reproduce the issue (python, gurobipy, gurobi-optimods) and the OS you are running?@dimendozao: Hi @simonbowly. Thanks for your reply. I used the following setup:

-  Windows 11 22H2
-  Python 3.11
-  Anaconda 3
-  Spyder IDE 5.4.3  
-  gurobipy version 10.0.2
-  gurobi-optimods version 1.1.0
-  numpy 1.24.3

Here is the .py file.

[IEEE69.zip](https://github.com/Gurobi/gurobi-optimods/files/13250484/IEEE69.zip)
@simonbowly: Hi @dimendozao, sorry for the delay! I could reproduce this with your code, thanks for sharing.

I believe this is linked to the rectangular coordinates representation currently used for the AC formulation. Gurobi doesn't report any large infeasibility in the solution itself, so the formulation would seem to be the issue here. We haven't added a polar form representation (which I believe MATPOWER uses?) to the optimod just yet, but this is planned for early next year.@100110110101: Hi @dimendozao & @simonbowly,

after talking with @JaromilNajman over the opf-optimod last week, I tried it myself and came across this issue here.

I was also able to reproduce the issue, although I found 3 slightly different solutions with Gurobi than dimendozao did:
Solution count 3: 80.886 80.9328 80.9432
Maybe this is because I am using Gurobi 11.0.1.
However, I think these different solutions are not the important point here. I think the important point here are the power differences at the buses, which show up for any arbitrary solution when the node powers are calculated from the generator and load powers on the one hand and from the grid equation system with the complex valued voltages on the other hand.
I did some more investigations because I couldn't believe that the use of rectangular coordinates could be the cause of such large deviations.

Long story short: I am pretty sure that the voltage angles are calculated incorrectly.
I didn't elaborate what the function `compute_voltage_angles(alldata, result)` in grbformulator.py is supposed to do in detail (I think I have a rough idea ...) but I think this function is way to complicated and maybe there is an error in it. However, I commented the according function-call in line 372 of grbformulator.py out and added the following simple line in the for-loop starting in line 345:
`resbus["Va"] = math.atan(alldata["LP"]["fvar"][databus].X/alldata["LP"]["evar"][databus].X)`
Furthermore, the angle is calculated with `ph[i]=np.radians(res_bus[i]['Va'])` in dimendozao's OptimodTest.py, but I think it simply needs to be `ph[i]=res_bus[i]['Va']`.
With both adjustments significantly better (and reasonable) values can be calculated and less differences result at each node. 

Can you reproduce the improvements with your code?

P.S.: It is still unsatisfactory that the remaining differences are so large. There is certainly still room for improvement.
From my experience with rectangular coordinates, I would generally recommend not to switch from rectangular coordinates to polar coordinates and back, as the (inverse) trigonometric functions entail numerical inaccuracies.@hhijazi: Hi folks, we believe that we found the issue and can get all violations down to zero. It has to do with using auxiliary variables in the power flow constraints, and this particular instance has numerically large coefficients (1e6) multiplying these variables.
The solution is to have the power flow constraints written using the original variables ((e,f) for complex V).
We are working on an updated model, stay tuned.@dimendozao: Hi everyone,

I am checking the response given by Thomas @100110110101 regarding the computation of phase angles.  However, I changed how the infeasibility was calculated. Instead of summing the absolute values I take the absolute value of the sum. The results change dramatically. 

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/240fed5c-a1ca-459e-b9c6-c0edbfed97f4)

If phase angles are  transformed to radians the infeasibility is:

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/0e616bd9-c7af-45b3-b676-6ea3b1b99fd8)

Notice that the absolute infeasibility (absinf) is considerably lower.

if phase angles are kept unchanged, as suggested by @100110110101, the infeasibility is  marginally greater:

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/6736c7e8-06c7-426c-8612-759a1ce58316)

This suggests that the phase angles are originally displayed in radians, and they are in general very close to zero.

Nevertheless, I agree with @100110110101 when he says that the infeasibility is still big (the apparent power error is around 26 KVA), specially when compared to the infeasibility shown using MATPOWER.

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/e1bf6f3e-8ae5-4eac-974f-4584cb953bc9).

Understanding that it is a non-convex problem (unless Jabr's relaxation is applied), the time it takes to solve the problem is also an issue, if again we compare it to MATPOWER.

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/128df361-0eae-474e-b4a0-192e6f4872ed)

![imagen](https://github.com/Gurobi/gurobi-optimods/assets/138025704/ef9bccac-d197-4174-b544-f8f4a67e404e)





Solution patch:
diff --git a/.github/workflows/code-quality.yml b/.github/workflows/code-quality.yml
index edc4fc5..b4194d6 100644
--- a/.github/workflows/code-quality.yml
+++ b/.github/workflows/code-quality.yml
@@ -27,7 +27,7 @@ jobs:
       - name: Lint with flake8
         run: |
           # check for Python syntax errors, undefined names, unused imports
-          python -m flake8 . --count --select=E9,F63,F7,F82,F811,F401 --show-source --statistics
+          python -m flake8 . --count --select=E9,F63,F7,F82,F811,F401,F811 --show-source --statistics
       - name: Verify pre-commit checks passed
         run: |
           pre-commit run --all-files
diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml
index aec83c1..15f21f9 100644
--- a/.github/workflows/python-tests.yml
+++ b/.github/workflows/python-tests.yml
@@ -1,4 +1,5 @@
-# Run the unit tests against python 3.8 + gurobipy 10.
+# Run the unit tests against oldest and newest supported python versions, and
+# major gurobipy versions.
 
 name: Unit tests
 
@@ -14,16 +15,11 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python: ["3.8", "3.12"]
-        gurobipy: ["10.0", "11.0"]
-        scipy: ["1.8.1", "1.13"]
+        python: ["3.9", "3.12"]
+        gurobipy: ["~=10.0", "~=11.0", "==12.0.0b3"]
         exclude:
           - python: "3.12"
-            gurobipy: "10.0"
-          - python: "3.8"
-            scipy: "1.13"
-          - python: "3.12"
-            scipy: "1.8.1"
+            gurobipy: "~=10.0"
 
     steps:
       - uses: actions/checkout@v4
@@ -37,8 +33,7 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install . \
-            gurobipy~=${{ matrix.gurobipy }} \
-            scipy~=${{ matrix.scipy }}
+            gurobipy${{ matrix.gurobipy }}
       - name: Run unittest
         run: python -m unittest discover -b
       - name: Install example dependencies
diff --git a/.github/workflows/wheel-tests.yml b/.github/workflows/wheel-tests.yml
index 97b9ea8..6446070 100644
--- a/.github/workflows/wheel-tests.yml
+++ b/.github/workflows/wheel-tests.yml
@@ -18,10 +18,11 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, windows-latest, macos-latest]
-        python: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python: ["3.9", "3.10", "3.11", "3.12"]
+        gurobipy: ["~=10.0", "~=11.0", "==12.0.0b3"]
         exclude:
-          - os: macos-latest
-            python: "3.8"
+          - python: "3.12"
+            gurobipy: "~=10.0"
 
     steps:
       - uses: actions/checkout@v4
@@ -37,11 +38,11 @@ jobs:
           python -m build
       - name: Fetch dependencies for offline install
         run: |
-          pip wheel --wheel-dir=dist gurobipy pandas numpy scipy gurobipy-pandas
+          pip wheel --wheel-dir=dist gurobipy${{ matrix.gurobipy }} pandas numpy scipy gurobipy-pandas
           pip wheel --wheel-dir=dist "numpy<2"
       - name: Install from built wheel
         run: |
-          python -m pip install --find-links dist --no-index --only-binary=:all: gurobi-optimods
+          python -m pip install --find-links dist --no-index --only-binary=:all: gurobi-optimods gurobipy${{ matrix.gurobipy }}
       - name: Run unit tests
         run: |
           python -m unittest discover -b
diff --git a/docs/source/mods/opf/opf.rst b/docs/source/mods/opf/opf.rst
index f3d9dc1..a74eec9 100644
--- a/docs/source/mods/opf/opf.rst
+++ b/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.
diff --git a/src/gurobi_optimods/opf/api.py b/src/gurobi_optimods/opf/api.py
index 9425de5..bab64e2 100644
--- a/src/gurobi_optimods/opf/api.py
+++ b/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
diff --git a/src/gurobi_optimods/opf/converters.py b/src/gurobi_optimods/opf/converters.py
index 936c841..d02dee1 100644
--- a/src/gurobi_optimods/opf/converters.py
+++ b/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,6 +178,16 @@ 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
@@ -181,12 +195,21 @@ def convert_case_to_internal_format(case_dict):
     # 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
diff --git a/src/gurobi_optimods/opf/grbformulator.py b/src/gurobi_optimods/opf/grbformulator.py
index cdf016c..556881b 100644
--- a/src/gurobi_optimods/opf/grbformulator.py
+++ b/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]
diff --git a/src/gurobi_optimods/opf/grbformulator_ac.py b/src/gurobi_optimods/opf/grbformulator_ac.py
index 436e0e6..d51a78a 100644
--- a/src/gurobi_optimods/opf/grbformulator_ac.py
+++ b/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,9 @@ def lpformulator_ac_create_constraints(alldata, model):
     IDtoCountmap = alldata["IDtoCountmap"]
     cvar = alldata["LP"]["cvar"]
     svar = alldata["LP"]["svar"]
+    if alldata["use_ef"]:
+        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"]
@@ -515,16 +515,27 @@ def lpformulator_ac_create_constraints(alldata, model):
                 Pvar_t[branch] == 0, name=branch.Ptcname
             )
             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)
+        if alldata["use_ef"]:
+            # Gff cff + Gft cft + Bft sft
+            expP = Pvar_f[branch]
+            if alldata["branchswitching_mip"]:
+                expP += twinPvar_f[branch]
+            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,
+            )
+        else:
+            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)
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -554,20 +565,32 @@ 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]],
-        )
+        if alldata["use_ef"]:
+            expP = Pvar_t[branch]
+            if alldata["branchswitching_mip"]:
+                expP += twinPvar_t[branch]
+            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,
+            )
+        else:
+            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]
+            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(expr == expP, name=branch.Ptcname)
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -623,15 +646,27 @@ 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]],
-        )
+        if alldata["use_ef"]:
+            expQ = Qvar_f[branch]
+            if alldata["branchswitching_mip"]:
+                expQ += twinQvar_f[branch]
+            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,
+            )
+        else:
+            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)
+            expQ = Qvar_f[branch]
+            if alldata["branchswitching_mip"]:
+                expQ += twinQvar_f[branch]
+            branch.Qdeffconstr = model.addConstr(expr == expQ, name=branch.Qfcname)
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -661,15 +696,27 @@ 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]],
-        )
+        if alldata["use_ef"]:
+            expQ = Qvar_t[branch]
+            if alldata["branchswitching_mip"]:
+                expQ += twinQvar_t[branch]
+            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,
+            )
+        else:
+            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)
+            expQ = Qvar_t[branch]
+            if alldata["branchswitching_mip"]:
+                expQ += twinQvar_t[branch]
+            branch.Qdeftconstr = model.addConstr(expr == expQ, name=branch.Qtcname)
 
         if alldata["branchswitching_mip"]:
             if branch.constrainedflow:
@@ -708,13 +755,16 @@ 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])
+            if alldata["use_ef"]:
+                expr.add(bus.Gs * (evar[bus] * evar[bus] + fvar[bus] * fvar[bus]))
+            else:
+                expr.add(bus.Gs * cvar[bus])
 
-        if alldata["branchswitching_comp"] is False:
+        if not alldata["branchswitching_comp"]:
             for branchid in bus.frombranchids.values():
                 expr.add(Pvar_f[branches[branchid]])
 
@@ -739,12 +789,14 @@ 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])
+            if alldata["use_ef"]:
+                expr.add(-bus.Bs * (evar[bus] * evar[bus] + fvar[bus] * fvar[bus]))
+            else:
+                expr.add(-bus.Bs * cvar[bus])
 
-        if alldata["branchswitching_comp"] is False:
+        if not alldata["branchswitching_comp"]:
             for branchid in bus.frombranchids.values():
                 expr.add(Qvar_f[branches[branchid]])
 
@@ -769,6 +821,51 @@ def lpformulator_ac_create_constraints(alldata, model):
 
     logger.info(f"    {count} balance constraints added.")
 
+    if alldata["skipjabr"] and alldata["use_ef"]:
+        # voltage magnitude bounds
+        logger.info("  Adding voltage magnitude constraints.")
+        count = 0
+        for j, bus in buses.items():
+            expr = 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:
+                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
@@ -798,7 +895,7 @@ def lpformulator_ac_create_constraints(alldata, model):
     logger.info("  Adding branch limits.")
     count = 0
     for j, branch in branches.items():
-        if not branch.status or branch.unboundedlimit is True:
+        if not branch.status or branch.unboundedlimit:
             continue
 
         f = branch.f
@@ -822,7 +919,7 @@ def lpformulator_ac_create_constraints(alldata, model):
     logger.info(f"    {count} branch limits added.")
 
     # JABR.
-    if alldata["skipjabr"] is False:
+    if not alldata["skipjabr"]:
         logger.info("  Adding Jabr constraints.")
         count = 0
         for j, branch in branches.items():
@@ -845,7 +942,7 @@ def lpformulator_ac_create_constraints(alldata, model):
         logger.info("  Skipping Jabr inequalities.")
 
     # Active loss constraints.
-    if alldata["useactivelossineqs"] is True:
+    if alldata["useactivelossineqs"]:
         logger.info("  Adding active loss constraints.")
         count = 0
         for j, branch in branches.items():
@@ -905,7 +1002,7 @@ def lpformulator_ac_create_constraints(alldata, model):
         lpformulator_ac_add_polarconstraints(alldata, model)
 
     # nonconvex e, f representation
-    if alldata["use_ef"]:
+    if not alldata["skipjabr"] and alldata["use_ef"]:
         lpformulator_ac_add_nonconvexconstraints(alldata, model)
 
     if alldata["branchswitching_mip"] or alldata["branchswitching_comp"]:
diff --git a/src/gurobi_optimods/opf/io.py b/src/gurobi_optimods/opf/io.py
index 3a8b979..46be45c 100644
--- a/src/gurobi_optimods/opf/io.py
+++ b/src/gurobi_optimods/opf/io.py
@@ -153,6 +153,7 @@ def read_case_matpower(file_path):
         "gen": gen.to_dict("records"),
         "branch": branch.to_dict("records"),
         "gencost": gencost,
+        "casename": mpc["casename"] if "casename" in mpc.dtype.names else "case",
     }
 
 
@@ -185,6 +186,7 @@ def write_case_matpower(case, file_path):
                 "gen": gen[gen_field_names].values,
                 "branch": branch[branch_field_names].values,
                 "gencost": gencost_values,
+                "casename": case.get("casename", "Undefined"),
             }
         },
     )
diff --git a/src/gurobi_optimods/opf/structs.py b/src/gurobi_optimods/opf/structs.py
index 04731fb..0e3f29b 100644
--- a/src/gurobi_optimods/opf/structs.py
+++ b/src/gurobi_optimods/opf/structs.py
@@ -212,11 +212,12 @@ class Branch:
         self.minangle = minangle
         self.minangle_rad = math.pi * minangle / 180.0
 
+        # We need to have -pi/2 < maxangle_rad < pi/2 since we use tan(angle) in the rectangular formulation
         if maxangle == 360 or maxangle == 0:
-            self.maxangle_rad = 2 * math.pi
+            self.maxangle_rad = math.pi / 2 - 1e-6
 
         if minangle == -360 or minangle == 0:
-            self.minangle_rad = -2 * math.pi
+            self.minangle_rad = -math.pi / 2 + 1e-6
 
         self.invratio2 = invratio2 = 1 / ratio**2
         self.multtf = multtf = 1 / (ratio * cmath.exp(1j * self.angle_rad))
diff --git a/src/gurobi_optimods/opf/violations.py b/src/gurobi_optimods/opf/violations.py
index 5f07249..e768e58 100644
--- a/src/gurobi_optimods/opf/violations.py
+++ b/src/gurobi_optimods/opf/violations.py
@@ -130,6 +130,7 @@ def grbderive_xtra_sol_values_from_voltages(alldata, model):
     Qvar_f = alldata["LP"]["Qvar_f"]
     Qvar_t = alldata["LP"]["Qvar_t"]
 
+    # Recomputing the active and reactive power flows using the voltage variables
     for j, branch in branches.items():
         for item in [
             (branch.Pdeffconstr, Pvar_f[branch]),
@@ -139,17 +140,34 @@ def grbderive_xtra_sol_values_from_voltages(alldata, model):
         ]:
             constr = item[0]
             var = item[1]
-
-            row = model.getRow(constr)
-            sum = -constr.RHS
             leadcoeff = 0
-            for i in range(row.size()):
-                v = row.getVar(i)
-                coeff = row.getCoeff(i)
-                if v.Varname != var.Varname:
-                    sum += coeff * xbuffer[v]
-                else:
-                    leadcoeff = coeff
+            # If the constraint is linear, e.g., Pvar_f = f(cosvar,sinvar), f linear
+            if type(constr) is gp.Constr:
+                row = model.getRow(constr)
+                sum = -constr.RHS
+                for i in range(row.size()):
+                    v = row.getVar(i)
+                    coeff = row.getCoeff(i)
+                    if v.Varname != var.Varname:
+                        sum += coeff * xbuffer[v]
+                    else:
+                        leadcoeff = coeff
+            else:  # Rectangular formulation with Pvar_f = f(evar,fvar), f quadratic
+                row = model.getQCRow(constr)
+                sum = -constr.QCRHS
+                for i in range(row.size()):
+                    v1 = row.getVar1(i)
+                    v2 = row.getVar2(i)
+                    coeff = row.getCoeff(i)
+                    sum += coeff * xbuffer[v1] * xbuffer[v2]
+                lterms = row.getLinExpr()
+                for i in range(lterms.size()):
+                    v = lterms.getVar(i)
+                    coeff = lterms.getCoeff(i)
+                    if v.Varname != var.Varname:
+                        sum += coeff * xbuffer[v]
+                    else:
+                        leadcoeff = coeff
 
             xbuffer[var] = -sum / leadcoeff
             # leadcoeff should be +1 or -1


Test patch:
diff --git a/tests/opf/data/bus1-branch1-gen1.json b/tests/opf/data/bus1-branch1-gen1.json
index a61842e..373bd51 100644
--- a/tests/opf/data/bus1-branch1-gen1.json
+++ b/tests/opf/data/bus1-branch1-gen1.json
@@ -1,4 +1,5 @@
 {
+    "casename": "us1-branch1-gen1",
     "baseMVA": 100.0,
     "branch": [
         {
diff --git a/tests/opf/data/bus2-branch2-gen2.json b/tests/opf/data/bus2-branch2-gen2.json
index cb75121..293bf77 100644
--- a/tests/opf/data/bus2-branch2-gen2.json
+++ b/tests/opf/data/bus2-branch2-gen2.json
@@ -1,4 +1,5 @@
 {
+    "casename": "bus2-branch2-gen2",
     "baseMVA": 100.0,
     "branch": [
         {
diff --git a/tests/opf/data/case9_switching_solution.json b/tests/opf/data/case9_switching_solution.json
index 4b4f0ab..3aa547f 100644
--- a/tests/opf/data/case9_switching_solution.json
+++ b/tests/opf/data/case9_switching_solution.json
@@ -1,4 +1,5 @@
 {
+    "casename": "case9",
     "baseMVA": 100.0,
     "branch": [
         {
diff --git a/tests/opf/test_io.py b/tests/opf/test_io.py
index f6a38ea..e4ecf9d 100644
--- a/tests/opf/test_io.py
+++ b/tests/opf/test_io.py
@@ -58,10 +58,20 @@ class TestBadData(unittest.TestCase):
 class TestDatasets(unittest.TestCase):
     def test_load_opf_example(self):
         # Check that example cases can be loaded by name
-        for case in ["case9", "case14", "case57", "case118", "case300", "caseNY"]:
+        for case in [
+            "pglib_opf_case3_lmbd",
+            "pglib_opf_case5_pjm",
+            "case9",
+            "case14",
+            "case57",
+            "case118",
+            "case300",
+            "caseNY",
+        ]:
             case_data = load_opf_example(case)
             self.assertEqual(
-                set(case_data.keys()), {"baseMVA", "bus", "branch", "gen", "gencost"}
+                set(case_data.keys()),
+                {"baseMVA", "bus", "branch", "gen", "gencost", "casename"},
             )
 
     def test_load_opf_extras_coordinates(self):
@@ -157,7 +167,8 @@ class TestIO(unittest.TestCase):
         # Check that all example cases are read without errors
 
         case_mat_files = [
-            self.dataset_dir.joinpath(f"pglib_case{case}.mat") for case in ["14"]
+            self.dataset_dir.joinpath(f"pglib_opf_case{case}.mat")
+            for case in ["3_lmbd", "5_pjm", "14_ieee"]
         ]
 
         for file_path in case_mat_files:
diff --git a/tests/opf/test_solver.py b/tests/opf/test_solver.py
index e2868c4..a172589 100644
--- a/tests/opf/test_solver.py
+++ b/tests/opf/test_solver.py
@@ -1,9 +1,12 @@
 # Tests that solve models using the public API
 
 import collections
+import math
 import random
 import unittest
 
+from gurobipy import GRB
+
 from gurobi_optimods.datasets import load_opf_example
 from gurobi_optimods.opf import solve_opf
 
@@ -146,12 +149,10 @@ class TestAPICase9(unittest.TestCase):
         self.assert_solution_valid(solution)
 
         self.assertIsNotNone(solution["f"])
-        self.assert_approx_equal(solution["f"], 5296.6862, tol=1.0)
+        self.assert_approx_equal(solution["f"], 5296.76, tol=1.0)
         self.assert_approx_equal(solution["bus"][0]["Va"], 0.0, tol=1e-1)
         self.assert_approx_equal(solution["gen"][1]["Qg"], 0.0318441, tol=1e-1)
-        self.assert_approx_equal(
-            solution["branch"][2]["Pt"], 55.96906046691643, tol=1e-1
-        )
+        self.assert_approx_equal(solution["branch"][2]["Pt"], 56.23, tol=1e-1)
 
     def test_ac_relax(self):
         solution = solve_opf(self.case, opftype="ACRelax")
@@ -162,6 +163,101 @@ class TestAPICase9(unittest.TestCase):
         self.assert_approx_equal(solution["gen"][0]["Pg"], 89.803524, tol=1e-1)
         self.assert_approx_equal(solution["branch"][1]["Pt"], -34.1774, tol=1e-1)
 
+    @unittest.skipIf(GRB.VERSION_MAJOR < 12, "Needs Gurobi 12")
+    def test_aclocal(self):
+        solution = solve_opf(self.case, opftype="aclocal")
+        self.assertEqual(solution["success"], 1)
+        self.assert_solution_valid(solution)
+
+        # A local solution should get somewhere close to optimal
+        self.assert_approx_equal(solution["f"], 5296.66532, tol=1e1)
+
+    def test_ac_branchswitching_infeasible(self):
+        # Modify the case to make it infeasible
+        self.case["bus"][1]["Vmax"] = 0.8
+
+        # Solve model, expect failure
+        with self.assertRaisesRegex(ValueError, "Infeasible model"):
+            solve_opf(self.case, opftype="AC", branch_switching=True)
+
+
+class TestAPICase5_PJM(unittest.TestCase):
+    # Test configurations of case5_pjm with known optimal values. We should
+    # get reasonable reproducibility on this small case.
+
+    def setUp(self):
+        self.case = load_opf_example("pglib_opf_case5_pjm")
+
+    def assert_approx_equal(self, value, expected, tol):
+        self.assertLess(abs(value - expected), tol)
+
+    def assert_solution_valid(self, solution):
+        # The solution carries over the same structure as the case input data,
+        # and some values (e.g. bus refs) must match exactly in the ordering.
+
+        self.assertEqual(
+            set(solution.keys()), set(self.case.keys()) | {"et", "success", "f"}
+        )
+        self.assertEqual(solution["baseMVA"], self.case["baseMVA"])
+
+        for sol_bus, case_bus in zip(solution["bus"], self.case["bus"]):
+            self.assertTrue(set(sol_bus.keys()).issuperset((set(case_bus.keys()))))
+            self.assertEqual(sol_bus["bus_i"], case_bus["bus_i"])
+
+        for sol_branch, case_branch in zip(solution["branch"], self.case["branch"]):
+            self.assertTrue(
+                set(sol_branch.keys()).issuperset((set(case_branch.keys())))
+            )
+            self.assertEqual(sol_branch["fbus"], case_branch["fbus"])
+            self.assertEqual(sol_branch["tbus"], case_branch["tbus"])
+
+        for sol_gen, case_gen in zip(solution["gen"], self.case["gen"]):
+            self.assertTrue(set(sol_gen.keys()).issuperset((set(case_gen.keys()))))
+            self.assertEqual(sol_gen["bus"], case_gen["bus"])
+
+        for sol_gencost, case_gencost in zip(solution["gencost"], self.case["gencost"]):
+            self.assertTrue(
+                set(sol_gencost.keys()).issuperset((set(case_gencost.keys())))
+            )
+            self.assertEqual(sol_gencost["costtype"], case_gencost["costtype"])
+            self.assertEqual(sol_gencost["n"], case_gencost["n"])
+
+    def test_dc(self):
+        solution = solve_opf(self.case, opftype="DC")
+        self.assertEqual(solution["success"], 1)
+        self.assert_solution_valid(solution)
+
+        self.assert_approx_equal(solution["f"], 17479.89, tol=1e-1)
+
+    def test_ac(self):
+        solution = solve_opf(self.case, opftype="AC")
+        self.assertEqual(solution["success"], 1)
+        self.assert_solution_valid(solution)
+
+        self.assert_approx_equal(solution["f"], 17551.89, tol=1e-1)
+
+    @unittest.skip("shaky")
+    def test_ac_branchswitching(self):
+        solution = solve_opf(
+            self.case,
+            opftype="AC",
+            branch_switching=True,
+            use_mip_start=False,
+            solver_params={"MIPGap": 1e-4},
+        )
+        self.assertEqual(solution["success"], 1)
+        self.assert_solution_valid(solution)
+
+        self.assertIsNotNone(solution["f"])
+        self.assert_approx_equal(solution["f"], 15174, tol=1)
+
+    def test_ac_relax(self):
+        solution = solve_opf(self.case, opftype="ACRelax")
+        self.assertEqual(solution["success"], 1)
+        self.assert_solution_valid(solution)
+
+        self.assert_approx_equal(solution["f"], 14999.71, tol=1e-1)
+
     def test_ac_branchswitching_infeasible(self):
         # Modify the case to make it infeasible
         self.case["bus"][1]["Vmax"] = 0.8
@@ -191,12 +287,12 @@ class TestComputeVoltageAnglesBug(unittest.TestCase):
         self.assertEqual(num_zeros, 1)
 
 
-class TestAPICase9Reordered(unittest.TestCase):
+class TestAPICase5_PJMReordered(unittest.TestCase):
     # It should be possible to re-order the input data (i.e. the lists of buses,
     # branches, generators) and get the same result
 
     def setUp(self):
-        self.case = load_opf_example("case9")
+        self.case = load_opf_example("pglib_opf_case5_pjm")
 
         self.bus_reorder = [i for i, _ in enumerate(self.case["bus"])]
         self.branch_reorder = [i for i, _ in enumerate(self.case["branch"])]
@@ -208,6 +304,7 @@ class TestAPICase9Reordered(unittest.TestCase):
 
         self.case_reordered = {
             "baseMVA": self.case["baseMVA"],
+            "casename": self.case["casename"],
             "bus": [dict(self.case["bus"][ind]) for ind in self.bus_reorder],
             "branch": [dict(self.case["branch"][ind]) for ind in self.branch_reorder],
             "gen": [dict(self.case["gen"][ind]) for ind in self.gen_reorder],
@@ -245,7 +342,9 @@ class TestAPICase9Reordered(unittest.TestCase):
             gen_original = solution_original["gen"][orig_ind]
             gen_reordered = solution_reordered["gen"][ind]
             for key, value in gen_original.items():
-                self.assert_approx_equal(gen_reordered[key], value, tol=1e-3)
+                # Some Gen entries are not initialized, e.g., Qc2min, check for nan
+                if not math.isnan(value):
+                    self.assert_approx_equal(gen_reordered[key], value, tol=1e-3)
 
     def test_ac(self):
         # Results are too numerically unstable, just check the call goes through
@@ -275,6 +374,17 @@ class TestAPICase9Reordered(unittest.TestCase):
             solution_original["f"], solution_reordered["f"], tol=1e1
         )
 
+    @unittest.skipIf(GRB.VERSION_MAJOR < 12, "Needs Gurobi 12")
+    def test_aclocal(self):
+        # Test AC local option, should arrive at a similar result
+        kwargs = dict(opftype="aclocal")
+        solution_original = solve_opf(self.case, **kwargs)
+        solution_reordered = solve_opf(self.case_reordered, **kwargs)
+
+        self.assert_approx_equal(
+            solution_original["f"], solution_reordered["f"], tol=1e1
+        )
+
 
 @unittest.skipIf(size_limited_license(), "size-limited-license")
 class TestAPILargeModels(unittest.TestCase):