# Gurobi.jl numerical trouble in optimal power flow

**URL:** <https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization, gurobi, modelling\
**Created:** [January 27, 2025, 2:49pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267 "2025-01-27T14:49:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [January 27, 2025, 2:49pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/1 "2025-01-27T14:49:00Z")

</div>

Hi. I am trying to model one of the ‘classic’ conic relaxations of the optimal powerflow model (see [Branch Flow Model: Relaxations and Convexification—Part I](https://ieeexplore.ieee.org/abstract/document/6507355))

The problem I am running into is that depending on the loads / size of the grid solving becomes very unreliable. I am using Gurobi. I will frequently see:

> **Error**
>
> ```julia
> Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (win64 - Windows 11.0 (22631.2))
> 
> CPU model: AMD Ryzen 9 5950X 16-Core Processor, instruction set [SSE2|AVX|AVX2]
> Thread count: 16 physical cores, 32 logical processors, using up to 32 threads
> 
> Optimize a model with 3264 rows, 3456 columns and 8064 nonzeros
> Model fingerprint: 0x5ec5409d
> Model has 384 quadratic constraints
> Coefficient statistics:
> Matrix range [1e+00, 7e+00]
> QMatrix range [1e+00, 1e+00]
> Objective range [1e+00, 1e+00]
> Bounds range [8e-01, 1e+00]
> RHS range [1e-02, 1e+00]
> Presolve removed 1920 rows and 1824 columns
> Presolve time: 0.00s
> Presolved: 1344 rows, 1728 columns, 4224 nonzeros
> Presolved model has 384 second-order cone constraints
> Ordering time: 0.00s
> 
> Barrier statistics:
> AA' NZ : 6.528e+03
> Factor NZ : 1.008e+04 (roughly 1 MB of memory)
> Factor Ops : 9.744e+04 (less than 1 second per iteration)
> Threads : 1
> 
> Objective Residual
> Iter Primal Dual Primal Dual Compl Time
> 0 1.70141978e+02 -1.24205980e+02 1.70e+00 7.06e-01 4.13e-01 0s
> 1 1.70141978e+02 -1.24205980e+02 1.70e+00 7.06e-01 4.13e-01 0s
> 2 2.00214585e+05 -1.83305591e+05 2.09e+03 1.00e+03 6.37e+05 0s
> 3 2.00214585e+05 -1.83305591e+05 2.09e+03 1.00e+03 6.37e+05 0s
> 4 2.00214585e+05 -1.83305591e+05 2.09e+03 1.00e+03 6.37e+05 0s
> 5 2.00214585e+05 -1.83305591e+05 2.09e+03 1.00e+03 6.37e+05 0s
> 
> Barrier performed 5 iterations in 0.02 seconds (0.01 work units)
> Numerical trouble encountered
> 
> Model may be infeasible or unbounded. Consider using the
> homogeneous algorithm (through parameter 'BarHomogeneous')
> 
> ```

To improven numerical stability I have rescaled everything to p.u. dimensions. But since I am new to this there’s a chance I made a typo/mistake somewhere.

Here is the model code. I set up a small test grid of just 4 buses, 2 identical loads. The real problem has about 100 buses.

```julia
test network with 4 buses 
      0 (slack)
      |
      1 (junction)
     / \
    2 3 (loads)

```

> **Model**
>
> ```julia
> using JuMP
> using Gurobi
> using Clarabel
> using Random
> 
> Random.seed!(333)
> 
> # p.u. base values
> S_base = 1E5 # [VA]
> V_base = 230.0 # [V]
> 
> Z_base = V_base^2 / S_base
> I_base = S_base / V_base
> S_max = 1E6 * 0.23 * 0.16 * 3 / S_base
> 
> # voltage constraints V_base ± 10% -> (207V, 253V) at 230V base
> V_lb = 0.90 # [pu]
> V_ub = 1.10
> 
> T = 1:96
> 
> function add_grid_model(;
> model::Model,
> T=1:96
> )
> # test network with 4 buses 
> # 0 (slack)
> # |
> # 1 (junction)
> # / \
> # 2 3 (loads)
> 
> B = [0, 1, 2, 3]
> H = [2, 3]
> L = [(0, 1), (1, 2), (1, 3)]
> 
> bus_out = Dict(0 => [1], 1 => [2, 3], 2 => [], 3 => [])
> bus_in = Dict(0 => [], 1 => [0], 2 => [1], 3 => [1])
> 
> # base loads 
> # P_base = 1e4 * ones(length(T)) / S_base
> # P_base = round.(1e3 .* abs.(randn(length(T))) / S_base .+ 3e4 / S_base, digits=3)
> P_base = 1e3 .* abs.(randn(length(T))) / S_base .+ 2e4 / S_base
> Q_base = zeros(length(T))
> 
> # slack bus ID
> SB = 0
> 
> # real and reactive line impedances
> R = Dict((i, j) => 0.1 / Z_base for (i, j) in L)
> X = Dict((i, j) => 0.1 / Z_base for (i, j) in L)
> Inom = Dict((i, j) => 1.0 for (i, j) in L)
> 
> # variables
> @variables(model, begin
> # voltage squared eq. (7)
> V_lb^2 <= v[B, T] <= V_ub^2
> 
> # bus power injections
> P[B, T], (base_name = "PBusInjection")
> Q[B, T], (base_name = "QBusInjection")
> 
> # power flows
> P_line[L, T], (base_name = "PLineFlow")
> Q_line[L, T], (base_name = "QLineFlow")
> 
> # line current
> 0 <= I_line[L, T], (base_name = "CurrentSquare")
> end)
> 
> # slack bus constraints
> @constraints(model, begin
> TrafoLimit[t in T], [S_max, P[SB, t], Q[SB, t]] in SecondOrderCone()
> end)
> 
> # bus / line constraints
> @constraints(model, begin
> # real power balance eq. (1)
> RealPowerBalance[j in B, t in T],
> P[j, t] == sum(P_line[(j, k), t] for k in bus_out[j]) -
> sum(P_line[(i, j), t] - R[(i, j)] * I_line[(i, j), t] for i in bus_in[j])
> 
> # reactive power balance eq. (2)
> ReactivePowerBalance[j in B, t in T],
> Q[j, t] == sum(Q_line[(j, k), t] for k in bus_out[j]) -
> sum(Q_line[(i, j), t] - X[(i, j)] * I_line[(i, j), t] for i in bus_in[j])
> 
> # voltage drop eq. (3)
> VoltageDrop[(i, j) in L, t in T],
> v[j, t] == v[i, t] - 2 * (P_line[(i, j), t] * R[(i, j)] + Q_line[(i, j), t] * X[(i, j)]) +
> (R[(i, j)]^2 + X[(i, j)]^2) * I_line[(i, j), t]
> 
> # conic OPF eq. (4)
> ConicOPF[(i, j) in L, t in T],
> [
> I_line[(i, j), t] + v[i, t],
> 2 * P_line[(i, j), t],
> 2 * Q_line[(i, j), t],
> I_line[(i, j), t] - v[i, t],
> ] in SecondOrderCone()
> 
> # line current limit eq. (6)
> LineCurrentLimit[(i, j) in L, t in T], I_line[(i, j), t] <= Inom[(i, j)]
> end)
> 
> # load constraints
> @constraints(model, begin
> Transmission[i in B, t in T; i ∉ H && i ≠ SB], P[i, t] == 0
> RealBaseLoad[i in H, t in T], P[i, t] == -P_base[t]
> ReactiveBaseLoad[i in H, t in T], Q[i, t] == -Q_base[t]
> end)
> end
> 
> # create the model
> # model = Model(Clarabel.Optimizer)
> model = Model(Gurobi.Optimizer)
> set_attribute(model, "BarHomogeneous", 1)
> 
> # add grid model 
> add_grid_model(model=model, T=T)
> 
> P = model[:P]
> 
> # define objective functions
> @expressions(model, begin
> J_gen, sum(P[0, T])
> end)
> 
> @objective(model, Min, J_gen)
> 
> # optimize the model and return it
> optimize!(model)
> @assert is_solved_and_feasible(model)
> 
> ```

And here is some plotting code I used, in case it’s useful.

> **Plotting**
>
> ```julia
> using CairoMakie
> 
> function plot_line_currents(model::Model, lines::Vector{Tuple{Int,Int}}=[(34, 0), (2, 74)])
> fig = Figure(; size=(1000, 600))
> 
> I_line = value.(model[:I_line])
> for (idx, (i, j)) in enumerate(lines)
> ax = Axis(fig[idx, 1], xlabel="Time [hours]",
> ylabel="Current [A]",
> title="($i -> $j)",
> xticks=(1:4:97, string.(0:1:24))
> )
> scatterlines!(ax, Vector(I_line[(i, j), :]), color=:blue, label="$i->$j", linewidth=2)
> end
> 
> fig
> end
> 
> function plot_pq(model::Model, bus::Int)
> fig = Figure(; size=(1000, 300))
> ax = Axis(fig[1, 1], xlabel="Time [hours]",
> ylabel="Power [p.u.]",
> title="Power flow",
> xticks=(1:4:97, string.(0:1:24))
> )
> 
> # index of slack bus (transformer)
> # SB = argmin(value.(model[:P]).axes[2])
> P_trafo = Matrix(value.(model[:P]))[bus, :]
> Q_trafo = Matrix(value.(model[:Q]))[bus, :]
> 
> # lines!(ax, sol[:P], color = :blue, label = "Transformer", linestyle = :so
> scatterlines!(ax, P_trafo, color=:blue, label="P trafo", linewidth=2)
> scatterlines!(ax, Q_trafo, color=:red, label="Q trafo", linewidth=2)
> axislegend(ax, position=:lb)
> 
> # display
> fig
> end
> 
> ```

I have read pretty much everything in [Guidelines for Numerical Issues - Gurobi Optimizer Reference Manual](https://docs.gurobi.com/projects/optimizer/en/current/concepts/numericguide.html). Based on the advice there I rescaled everything to p.u. domain and the coefficient ranges are quite decent now, but the problems remain.

Any advice would be welcome. Also, if any PowerModels.jl users / experts are reading this please let me know if you see any mistakes in the formulations.

ps. I know I could use PowerModels.jl instead of doing this myself, but this is a better learning experience.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [January 27, 2025, 7:33pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/2 "2025-01-27T19:33:08Z")

</div>

For testing purposes I also ran the model using Clarabel.jl, the only change to the code being:

```julia
model = Model(Clarabel.Optimizer)

```

The Clarabel log, in case anyone spots something suspicious:

```julia
-------------------------------------------------------------
           Clarabel.jl v0.9.0 - Clever Acronym
                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables = 2016
  constraints = 4320
  nnz(P) = 0
  nnz(A) = 7680
  cones (total) = 387
    : Zero = 1, numel = 1536
    : Nonnegative = 2, numel = (672,672)
    : SecondOrder = 384, numel = (3,3,3,3,...,4)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf, max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
  0 4.7111e+01 -7.3618e+02 1.66e+01 5.89e-01 3.78e-01 1.00e+00 1.79e+00 ------
  1 5.3901e+01 -6.3766e+01 2.18e+00 1.44e-01 6.25e-02 2.55e-01 5.40e-01 8.50e-01
  2 4.6217e+01 3.4599e+01 3.36e-01 1.64e-02 8.08e-03 5.80e-02 8.34e-02 8.87e-01
  3 4.4820e+01 4.2589e+01 5.24e-02 3.18e-03 1.56e-03 9.80e-03 1.68e-02 8.34e-01
  4 4.4579e+01 4.4241e+01 7.63e-03 4.70e-04 2.34e-04 1.19e-03 2.54e-03 8.90e-01
  5 4.4573e+01 4.4548e+01 5.67e-04 3.52e-05 1.74e-05 8.60e-05 1.91e-04 9.28e-01
  6 4.4573e+01 4.4572e+01 2.05e-05 1.27e-06 6.28e-07 1.32e-06 6.87e-06 9.90e-01
  7 4.4573e+01 4.4573e+01 5.93e-07 3.68e-08 1.82e-08 3.82e-08 1.99e-07 9.71e-01  
  8 4.4573e+01 4.4573e+01 2.99e-08 1.85e-09 9.16e-10 1.92e-09 1.00e-08 9.50e-01
  9 4.4573e+01 4.4573e+01 1.70e-09 1.05e-10 5.20e-11 1.04e-10 5.69e-10 9.48e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 17.7ms

```

Oddly enough Clarabel has no issues with this model at all. I didn’t expect Clarabel to succeed where Gurobi fails.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [January 27, 2025, 9:08pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/3 "2025-01-27T21:08:39Z")

</div>

> Model may be infeasible or unbounded. Consider using the  
> homogeneous algorithm (through parameter ‘BarHomogeneous’)

Did you try this suggestion?

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [January 27, 2025, 9:16pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/4 "2025-01-27T21:16:53Z")

</div>

Yep. I edited the code to reflect that. I also fixed the random number generator to a seed that will always cause Gurobi to fail.

---

<div class="post-metadata">

**Author:** ![torressa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torressa/32/202204_2.png) [@torressa](https://discourse.julialang.org/u/torressa)\
**Post date:** [January 27, 2025, 9:20pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/5 "2025-01-27T21:20:35Z")

</div>

I tested this a bit earlier, it’s only reproducible on windows/mac (tested over 500 seeds on linux). For now, the best I could come up with is setting Presolve=0, please test this. But this won’t be great for larger models.

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [January 27, 2025, 9:29pm UTC](https://discourse.julialang.org/t/gurobi-jl-numerical-trouble-in-optimal-power-flow/125267/6 "2025-01-27T21:29:06Z")

</div>

Very odd. Thank you for looking into this.

`set_attribute(model, "Presolve", 0)` seems to work for both the larger and the smaller model.

I am still not sure where the numerical instability is coming from though.