infeasibility of LPACCPowerModel for PEGASE 89 bus case

[grecht]

Neither IPOPT nor Gurobi converge to a feasible solution on this test case. I tried generating an irreducible inconsistent subsystem using Gurobi, but apparently JuMP needs a GenericModel to do this.

using PowerModels
using Ipopt

nw = parse_file("pglib_opf_case89_pegase.m")
res = solve_opf(nw, LPACCPowerModel, Ipopt.Optimizer)

This is the output after IPOPT stops:

Number of Iterations....: 478

                                   (scaled)                 (unscaled)
Objective...............:   2.5403995113601827e+03    1.0744328603513892e+05
Dual infeasibility......:   9.9999999740854648e+01    4.2293853890397450e+03
Constraint violation....:   1.4724721180794376e-03    1.4724721180794376e-03
Variable bound violation:   8.5406254152209371e-09    8.5406254152209371e-09
Complementarity.........:   7.5127157862845708e-09    3.1774170460861485e-07
Overall NLP error.......:   9.9999999740854648e+01    4.2293853890397450e+03


Number of objective function evaluations             = 2329
Number of objective gradient evaluations             = 436
Number of equality constraint evaluations            = 2329
Number of inequality constraint evaluations          = 2329
Number of equality constraint Jacobian evaluations   = 491
Number of inequality constraint Jacobian evaluations = 491
Number of Lagrangian Hessian evaluations             = 479
Total seconds in IPOPT                               = 2.204

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

Interestingly, the objective value is quite close to the AC-OPF solution (1.07285e5). Maybe this just requires some different solver settings, but briefly playing around with lower IPOPT tolerances didn't do anything.

Using my own implementation of the model I was able to compute an IIS using CPLEX, in case it helps. Ignore the first argument, which is a snapshot index (thus always "0" here). It narrows the infeasibility down to the buses 7637, 8581 and the branches 95=(5416, 7637), 205=(7637, 8581). Looks to me like the demands cannot be met at these nodes.

Number of equations in conflict: 15
fixed: flow_pf(0,95,5416,7637) = 11.7706
fixed: flow_pf(0,205,7637,8581) = 0.856809
fixed: flow_qf(0,95,5416,7637) = 111.884
fixed: flow_qf(0,205,7637,8581) = 64.5153
fixed: flow_pt(0,95,5416,7637) = 11.7706
fixed: flow_pt(0,205,7637,8581) = -0.107518
fixed: flow_qt(0,95,5416,7637) = 111.884
fixed: flow_qt(0,205,7637,8581) = 64.5209
fixed: nodal_balance_p(0,7637) = 0.0029
fixed: nodal_balance_p(0,8581) = -12.9913
fixed: nodal_balance_q(0,7637) = -0.0618
fixed: nodal_balance_q(0,8581) = -1.4085
upper: cosine_relaxation(0,95,5416,7637) < 1
upper: cosine_relaxation(0,205,7637,8581) < 0.999973
upper: thermal_limit_f(0,95) < 169.781

Number of variables in conflict: 0

Removing e.g. the active power demand at node 8581 (load 33) makes it feasible in PowerModels.

[grecht]

Never mind: It is feasible using Gurobi with a rather low convergence tolerance of 1e-4.

res = solve_opf(nw, LPACCPowerModel, optimizer_with_attributes(Gurobi.Optimizer, "BarQCPConvTol"=>1e-4))

I am closing this issue, however it would still be interesting to me how this could be related to e.g. the scaling of the problem and if there are any best modeling practices when using the QPAC model in particular.

[ccoffrin]

Yes, I agree this is not a bug in PowerModels.

For posterity, I did want to follow up and say that I am able to reproduce the result you show here, the LPAC-Cold model is not feasible for this network. This result is not necessarily a bug because an approximation for the AC equations (which the LPAC is) can cut off feasible solutions, which is the case in this example. This is a good case for highlighting the benefits of going for a full AC OPF model.