Why is IPOPT executed repeatedly for OPF on MultiNetworks from PowerModels.jl?

beccy_julia · January 19, 2025, 11:26pm

Dear community, I am using Multi Networks from PowerModels.jl on which I would like to run Optimal Power Flow (OPF). I use the MINLP solver Juniper with IPOPT as the nlp-solver. As a test grid I take the IEEE case14 and added two storage units at buses 6 and 12, and introduce negative loads (renewables = RES) as well as vary the load data over time. If I then replicate that grid and use “solve_mn_opf_strg” for OPF, and set the print output for IPOPT on 3, I see that IPOPT is executed repeatedly, each time finding an optimal solution. The number of variables varies a little (±10), the number of executions increases exponentially with the number of time steps (=replica), and also with how volatile the data is.

Storage units data:

%% storage data
%   storage_bus ps qs energy  energy_rating charge_rating  discharge_rating  charge_efficiency  discharge_efficiency  thermal_rating  qmin  qmax  r  x  p_loss  q_loss  status
mpc.storage = [
	 6	 0.0	 0.0	 20.0	 100.0	 50.0	 70.0	 0.8	 0.9	 100.0	 -50.0	 70.0	0.1	 0.0	 0.0	 0.0	 1;
	 12	 0.0	 0.0	 30.0	 100.0	 50.0	 70.0	 0.9	 0.8	 100.0	 -50.0	 70.0	0.1	 0.0	 0.0	 0.0	 1;
];

Further modifications (for completeness):

delete loads at buses 6,9,10,11
move load at bus 5 to bus 3
load at bus 2 (Pd: 45, Qd: 15), bus 3 (Pd: 30, Qd:10), bus 4 (Pd: 70, Qd:25)
neg. load at bus 12 (Pd: -30, Qd: -10), bus 14 (Pd: -20, Qd: -7)
Line ratings are 100, and 30 on lines 9,10,12,14, and 10 on branch 19.

These modifications should not be relevant to the described problem.

The optimizer is declared as follows:

using JuMP, Juniper, IPOPT
optimizer = Juniper.Optimizer
    nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3)
    minlp_solver = optimizer_with_attributes(optimizer, "nl_solver" => nl_solver, "log_levels" => [])

The execution of the OPF boils down to:

using PowerModels
T=10
mpc = PowerModels.parse_file("data/base_grids/IEEE/case14_strg.m") # load
mpc = replicate(mpc, T) # Multi Network
res = PowerModels.solve_mn_opf_strg(mpc, PowerModels.ACPPowerModel, minlp_solver) # OPF

Does anybody know why IPOPT runs so often?

odow · January 19, 2025, 11:34pm

Hi @beccy_julia, welcome to the forum

You’re using Juniper.jl as the solver.

Juniper is a solver for mixed-integer nonlinear programs. Juniper works via a branch-and-bound algorithm that searches over values of the discrete variables and, at each node, it solves a continuous nonlinear subproblem. The nl_solver attribute to Juniper controls the subsolver. In this example, you have chosen Ipopt. Juniper might also solve nonlinear subproblems as part of a heuristic, or an attempt to construct a feasible starting solution.

You can disable the printing from Ipopt with:

    nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

If you can provide a fully reproducible code example (with the data), then we might be able to help more.

beccy_julia · January 20, 2025, 1:08am

Dear @odow, thanks for the quick reply. I guessed I should provide the whole code, so here is a simple (silly) example that shows what I mean.

First, the IEEE case14 casefile including storage - to be saved as “case14.m” (MatPower/MATLAB):

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1.06	0	0	1	1.06	0.94;
	2	2	21.7	12.7	0	0	1	1.045	-4.98	0	1	1.06	0.94;
	3	2	94.2	19	0	0	1	1.01	-12.72	0	1	1.06	0.94;
	4	1	47.8	-3.9	0	0	1	1.019	-10.33	0	1	1.06	0.94;
	5	1	7.6	1.6	0	0	1	1.02	-8.78	0	1	1.06	0.94;
	6	2	11.2	7.5	0	0	1	1.07	-14.22	0	1	1.06	0.94;
	7	1	0	0	0	0	1	1.062	-13.37	0	1	1.06	0.94;
	8	2	0	0	0	0	1	1.09	-13.36	0	1	1.06	0.94;
	9	1	29.5	16.6	0	19	1	1.056	-14.94	0	1	1.06	0.94;
	10	1	9	5.8	0	0	1	1.051	-15.1	0	1	1.06	0.94;
	11	1	3.5	1.8	0	0	1	1.057	-14.79	0	1	1.06	0.94;
	12	1	6.1	1.6	0	0	1	1.055	-15.07	0	1	1.06	0.94;
	13	1	13.5	5.8	0	0	1	1.05	-15.16	0	1	1.06	0.94;
	14	1	14.9	5	0	0	1	1.036	-16.04	0	1	1.06	0.94;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	232.4	-16.9	10	0	1.06	100	1	332.4	0	0	0	0	0	0	0	0	0	0	0	0;
	2	40	42.4	50	-40	1.045	100	1	140	0	0	0	0	0	0	0	0	0	0	0	0;
	3	0	23.4	40	0	1.01	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
	6	0	12.2	24	-6	1.07	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
	8	0	17.4	24	-6	1.09	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	2	0.01938	0.05917	0.0528	0	0	0	0	0	1	-360	360;
	1	5	0.05403	0.22304	0.0492	0	0	0	0	0	1	-360	360;
	2	3	0.04699	0.19797	0.0438	0	0	0	0	0	1	-360	360;
	2	4	0.05811	0.17632	0.034	0	0	0	0	0	1	-360	360;
	2	5	0.05695	0.17388	0.0346	0	0	0	0	0	1	-360	360;
	3	4	0.06701	0.17103	0.0128	0	0	0	0	0	1	-360	360;
	4	5	0.01335	0.04211	0	0	0	0	0	0	1	-360	360;
	4	7	0	0.20912	0	0	0	0	0.978	0	1	-360	360;
	4	9	0	0.55618	0	0	0	0	0.969	0	1	-360	360;
	5	6	0	0.25202	0	0	0	0	0.932	0	1	-360	360;
	6	11	0.09498	0.1989	0	0	0	0	0	0	1	-360	360;
	6	12	0.12291	0.25581	0	0	0	0	0	0	1	-360	360;
	6	13	0.06615	0.13027	0	0	0	0	0	0	1	-360	360;
	7	8	0	0.17615	0	0	0	0	0	0	1	-360	360;
	7	9	0	0.11001	0	0	0	0	0	0	1	-360	360;
	9	10	0.03181	0.0845	0	0	0	0	0	0	1	-360	360;
	9	14	0.12711	0.27038	0	0	0	0	0	0	1	-360	360;
	10	11	0.08205	0.19207	0	0	0	0	0	0	1	-360	360;
	12	13	0.22092	0.19988	0	0	0	0	0	0	1	-360	360;
	13	14	0.17093	0.34802	0	0	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0.01	10	0;
	2	0	0	3	0.07	20	0;
	2	0	0	3	0.05 	50	0;
	2	0	0	3	0.04 	40	0;
	2	0	0	3	0.02 	30	0;
];

% hours
mpc.time_elapsed = 1.0

%% storage data
%   storage_bus ps qs energy  energy_rating charge_rating  discharge_rating  charge_efficiency  discharge_efficiency  thermal_rating  qmin  qmax  r  x  p_loss  q_loss  status
mpc.storage = [
	 6	 0.0	 0.0	 20.0	 100.0	 50.0	 70.0	 0.8	 0.9	 100.0	 -50.0	 70.0	0.1	 0.0	 0.0	 0.0	 1;
	 7	 0.0	 0.0	 30.0	 100.0	 50.0	 70.0	 0.9	 0.8	 100.0	 -50.0	 70.0	0.1	 0.0	 0.0	 0.0	 1;
];

Then, the code that runs optimal power flow (OPF):

using PowerModels, JuMP, Ipopt, Juniper

# optimizer
optimizer = Juniper.Optimizer
nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3)
minlp_solver = optimizer_with_attributes(optimizer, "nl_solver" => nl_solver, "log_levels" => [])

# time horizon (replica)
# T = 5
T = 15
# T = 30

# Multinetwork
net_single = PowerModels.parse_file("case14.m")
net_multi = replicate(net_single , T)

# OPF
results = PowerModels.solve_mn_opf_strg(net_multi, PowerModels.ACPPowerModel, minlp_solver)

The output is that it does converge, however there are several iterations of IPOPT (with larger T more iterations), like this:

[...]

Total number of variables............................:     1585
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      457
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1404
Total number of inequality constraints...............:      312
        inequality constraints with only lower bounds:       48
   inequality constraints with lower and upper bounds:      240
        inequality constraints with only upper bounds:       24


Number of Iterations....: 17

                                   (scaled)                 (unscaled)
Objective...............:   8.2689288055489249e+02    4.1344643945171869e+04
Dual infeasibility......:   3.4937205317300841e-13    1.7468602623762453e-11
Constraint violation....:   9.1315843775419125e-15    9.1315843775419125e-15
Variable bound violation:   1.0499368707783674e-08    1.0499368707783674e-08
Complementarity.........:   9.3111609659656525e-10    4.6555804736847871e-08
Overall NLP error.......:   9.3111609659656525e-10    4.6555804736847871e-08


Number of objective function evaluations             = 21
Number of objective gradient evaluations             = 18
Number of equality constraint evaluations            = 21
Number of inequality constraint evaluations          = 21
Number of equality constraint Jacobian evaluations   = 18
Number of inequality constraint Jacobian evaluations = 18
Number of Lagrangian Hessian evaluations             = 17
Total seconds in IPOPT                               = 0.827

EXIT: Optimal Solution Found.
Total number of variables............................:     1584
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      456
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1404
Total number of inequality constraints...............:      312
        inequality constraints with only lower bounds:       48
   inequality constraints with lower and upper bounds:      240
        inequality constraints with only upper bounds:       24


Number of Iterations....: 19

                                   (scaled)                 (unscaled)
Objective...............:   8.2689420825235823e+02    4.1344710330045025e+04
Dual infeasibility......:   3.3940203096770292e-13    1.6970101514492774e-11
Constraint violation....:   7.6327832942979512e-15    7.6327832942979512e-15
Variable bound violation:   1.0500180280814675e-08    1.0500180280814675e-08
Complementarity.........:   1.1374174929826561e-09    5.6870874535551351e-08
Overall NLP error.......:   1.1374174929826561e-09    5.6870874535551351e-08


Number of objective function evaluations             = 25
Number of objective gradient evaluations             = 20
Number of equality constraint evaluations            = 25
Number of inequality constraint evaluations          = 25
Number of equality constraint Jacobian evaluations   = 20
Number of inequality constraint Jacobian evaluations = 20
Number of Lagrangian Hessian evaluations             = 19
Total seconds in IPOPT                               = 0.538

EXIT: Optimal Solution Found.
Total number of variables............................:     1584
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      456
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1404
Total number of inequality constraints...............:      312
        inequality constraints with only lower bounds:       48
   inequality constraints with lower and upper bounds:      240
        inequality constraints with only upper bounds:       24


Number of Iterations....: 17

                                   (scaled)                 (unscaled)
Objective...............:   8.2689288055486543e+02    4.1344643945170516e+04
Dual infeasibility......:   4.6103210505637627e-13    2.3051605206780574e-11
Constraint violation....:   8.5487172896137054e-15    8.5487172896137054e-15
Variable bound violation:   1.0499368707783674e-08    1.0499368707783674e-08
Complementarity.........:   9.3111301238100338e-10    4.6555650526070085e-08
Overall NLP error.......:   9.3111301238100338e-10    4.6555650526070085e-08


Number of objective function evaluations             = 21
Number of objective gradient evaluations             = 18
Number of equality constraint evaluations            = 21
Number of inequality constraint evaluations          = 21
Number of equality constraint Jacobian evaluations   = 18
Number of inequality constraint Jacobian evaluations = 18
Number of Lagrangian Hessian evaluations             = 17
Total seconds in IPOPT                               = 0.492

EXIT: Optimal Solution Found.

Dict{String, Any} with 8 entries:
  "solve_time"         => 90.632
  "optimizer"          => "Juniper"
  "termination_status" => LOCALLY_SOLVED
  "dual_status"        => NO_SOLUTION
  "primal_status"      => FEASIBLE_POINT
  "objective"          => 41344.7
  "solution"           => Dict{String, Any}("multiinfrastructure"=>false, "multinetwork"=>true, "nw"=>Dict{String, Any}("4"=>Dict{String, Any}("baseMVA"=>100, "branch"=>Dict{String, Any}("4"=>Dict{String, Any}("qf"=>-0.0149719, "qt"=>0.0252754, "pt"=>-0.511128, "pf"=>0.526193), "1"=>Dict{String, Any}("qf"=>-0.…  "objective_lb"       => 41344.7

My problem is that, if I have a more complex problem (volatile and negative loads, curtailment generators, more storage units), it takes a lot longer to compute even at T=12 (several hours to not finishing…).

Hence, my question is, whether this is as it should be or whether I overlooked something that makes my code faster (e.g. Multithreading doesn’t work here)?

Thanks a lot!

odow · January 20, 2025, 1:24am

This behavior is expected. If you turn off logging of Ipopt and turn on logging from Juniper, you will see:

julia> using PowerModels, JuMP, Ipopt, Juniper

julia> solver = optimizer_with_attributes(
           Juniper.Optimizer,
           "nl_solver" => optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0),
       );

julia> net_single = PowerModels.parse_file("case14.m");
[... lines of warnings removed ...]

julia> net_multi = replicate(net_single , 15);

julia> results = PowerModels.solve_mn_opf_strg(
           net_multi,
           PowerModels.ACPPowerModel,
           solver,
       )
nl_solver         : MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 0])
feasibility_pump  : false
log_levels        : [:Options, :Table, :Info]

#Variables: 2010
#IntBinVar: 60
Obj Sense: Min

Start values are not feasible.
Status of relaxation: LOCALLY_SOLVED
Time for relaxation: 0.8341710567474365
Relaxation Obj: 51860.75540370797

 ONodes   CLevel          Incumbent                   BestBound            Gap    Time   Restarts  GainGap  
============================================================================================================
    2       2                 -                        51860.76             -     19.8      0         -     
    3       3                 -                        51860.76             -     20.2      -       83.2%   
    4       4                 -                        51860.76             -     20.5      -       83.2%   
    5       5                 -                        51860.76             -     20.8      -       83.2%   
    6       6                 -                        51860.76             -     21.2      -       83.2%   
    7       7                 -                        51860.76             -     21.5      -       83.2%   
    8       8                 -                        51860.76             -     21.8      -       83.2%   
    9       9                 -                        51860.76             -     22.1      -       83.2%   
   10       10                -                        51860.76             -     22.5      -       83.2%   
   11       11                -                        51860.76             -     22.8      -       83.2%   
   12       12                -                        51860.76             -     23.1      -       83.2%   
   13       13                -                        51860.76             -     23.5      -       83.2%   
   14       14                -                        51860.76             -     23.8      -       83.2%   
   15       15                -                        51860.76             -     24.2      -       83.2%   
   16       16                -                        51860.76             -     24.5      -       83.2%   
   17       17                -                        51860.76             -     24.8      -       83.2%   
   18       18                -                        51860.76             -     25.2      -       83.2%   
   19       19                -                        51860.76             -     25.6      -       83.2%   
   20       20                -                        51860.76             -     25.9      -       83.2%   
   21       21                -                        51860.76             -     26.3      -       83.2%   
   22       22                -                        51860.76             -     26.7      -       83.2%   
   23       23                -                        51860.76             -     27.0      -       83.2%   
   24       24                -                        51860.76             -     27.4      -       83.2%   
   25       25                -                        51860.76             -     27.7      -       83.2%   
   26       26                -                        51860.76             -     28.1      -       83.2%   
   27       27                -                        51860.76             -     28.5      -       83.2%   
   28       28                -                        51860.76             -     28.8      -       83.2%   
   29       29                -                        51860.76             -     29.2      -       83.2%   
   30       30                -                        51860.76             -     29.5      -       83.2%   
    0       31             51860.8                     51860.76           0.0%    29.9      -       83.2%   

#branches: 30
Obj: 51860.79722002109
Dict{String, Any} with 8 entries:
  "solve_time"         => 30.6949
  "optimizer"          => "Juniper"
  "termination_status" => LOCALLY_SOLVED
  "dual_status"        => NO_SOLUTION
  "primal_status"      => FEASIBLE_POINT
  "objective"          => 51860.8
  "solution"           => Dict{String, Any}("multiinfrastructure"=>false, "multinetwork"=>true, "nw"=>Dict{String, Any}("4"=>Dict{String, Any}("baseMVA"=>100, "branch"=>Dict{String, Any}("4"=>Dict{String…
  "objective_lb"       => 51860.8

In general, Juniper may require hundreds or thousands (if not more!) of Ipopt subproblem solves to find a solution. For larger MINLPs, it may be impossible for Juniper to find a solution. MINLPs are very difficult to solve.

beccy_julia · January 20, 2025, 2:00am

I see, thanks for the explanation! And the tipp with Juniper output. That calms my nerves. ^^

odow · January 20, 2025, 2:09am

If you have a license, you might also consider trying Gurobi.jl for this problem. When I tried it quickly found a good incumbent, but it struggled to prove optimality, so consider running it with a time limit.

julia> solver = optimizer_with_attributes(Gurobi.Optimizer, MOI.TimeLimitSec() => 60.0)
MathOptInterface.OptimizerWithAttributes(Gurobi.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.TimeLimitSec() => 60.0])

julia> results = PowerModels.solve_mn_opf_strg(
                 net_multi,
                 PowerModels.ACPPowerModel,
                 solver,
             )

Topic		Replies	Views
Why is Ipopt executed repeatedly for MultiNetworks from PowerModels.jl? Optimization (Mathematical) question , jump	0	16	January 19, 2025
Newton solution from PowerModels.jl Optimization (Mathematical)	11	1174	July 23, 2021
PowerModels.jl multinetwork for Security Constrained AC-OPF Optimization (Mathematical) question	4	289	July 21, 2024
Sudden stop of parallel optimization with JuMP and Juniper Optimization (Mathematical) jump , parallel	8	789	July 23, 2021
ACOPF Problem Optimization (Mathematical) query	2	546	December 1, 2021

Why is IPOPT executed repeatedly for OPF on MultiNetworks from PowerModels.jl?

Related topics