OPF problem

[dimendozao]

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

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.

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).

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$$

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$$

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:

The objective differs from MATPOWER, but the solution found is optimal.

This solution pass the first feasibility check:

However the second feasibility check is failed:

[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

[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.

If phase angles are transformed to radians the infeasibility is:

Notice that the absolute infeasibility (absinf) is considerably lower.

if phase angles are kept unchanged, as suggested by @100110110101, the infeasibility is marginally greater:

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.

.

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.