Min max Ipopt opti : error in step computation

Specific Domains Optimization (Mathematical)
, ,
1

Hello everybody.

I want to solve a min max problem with Ipopt solver. I used the trick described in JuMP documentation (section tips and tricks).

I managed to write my problem but the optimization do not work. I hade the message EXIT: Restoration Failed!. I read that this could come from division by 0 issue. I then added large upper bounds on my variables and indicated a starting value. I obtained another error : EXIT: Error in step computation!

My problem has one analytical solution which is x=[cos (theta); sin(theta)] with theta=arctan(2.2);

The code is written in order to be general as the size of the vector may be larger than 2, once I will have this fixed.

Here is the code.

using  JuMP, Clarabel, Ipopt;

A=[  6.6  -3.0;-3.0   6.6];
nc=1;
tol=1.e-9
ub=1.e9;

model = Model(Ipopt.Optimizer)

@variable(model, x[1:2*nc], start=1);
@constraint(model,sum( x[i]*x[i] for i=1:2*nc) == 1 );
@constraint(model,[i=1:nc],tol<=sum( x[2*(i-1)+1]*x[2*(i-1)+1] + x[2*i]*x[2*i] ) <=ub);
@constraint(model,[i=1:nc],tol<=sum( x[2*(i-1)+1]*x[2*(i-1)+1] + x[2*i]*x[2*i] ) <=ub);

@variable(model, u[1:2*nc], start=1);
@constraint(model, [i=1:2*nc], sum( A[i,k]*x[k] for k=1:nc) == u[i] ); 
@constraint(model,[i=1:nc], u[2*(i-1)+1] == 0);
@constraint(model,[i=1:nc],tol<=u[2*i]*u[2*i]<=ub);

@variable(model, ratio[1:nc], start=1);
@constraint(model,[i=1:nc]  , ((x[2*i]  *u[2*i])/  ((sum( x[2*(i-1)+1]*x[2*(i-1)+1] + x[2*i]*x[2*i] ) ).^0.5 * (( u[2*i]*u[2*i] ).^0.5 )        )   ==ratio[i]));


@variable(model,t, start=1);


@constraint(model, [t; ratio] in MOI.NormInfinityCone(1 + length(ratio)))
@objective(model, Min,t);

optimize!(model)
objective_value(model)
JuMP.value.(x)

and here is what I obtained


Number of nonzeros in equality constraint Jacobian...:       11
Number of nonzeros in inequality constraint Jacobian.:        9
Number of nonzeros in Lagrangian Hessian.............:       13

Total number of variables............................:        6
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        5
Total number of inequality constraints...............:        5
        inequality constraints with only lower bounds:        2
   inequality constraints with lower and upper bounds:        2
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+00 5.60e+00 3.33e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.1658618e+00 2.83e+00 1.03e+00  -1.0 2.00e+00    -  8.01e-01 4.95e-01h  1
   2  1.1648725e+00 2.72e+00 5.60e+01  -1.0 5.05e-01    -  9.90e-01 3.74e-02h  1
   3  1.1648651e+00 2.72e+00 1.40e+05  -1.0 4.86e-01    -  9.90e-01 4.19e-04h  1
   4r 1.1648651e+00 2.72e+00 1.00e+03   0.4 0.00e+00    -  0.00e+00 2.62e-07R  5
   5r 1.2586349e+00 7.28e-01 9.83e+02   0.4 2.56e+03    -  1.82e-02 1.05e-03f  1
   6  1.2586322e+00 7.28e-01 8.02e+04  -1.0 4.62e-01    -  9.90e-01 1.79e-05h  1
   7r 1.2586322e+00 7.28e-01 1.00e+03  -0.1 0.00e+00    -  0.00e+00 3.10e-07R  3
   8r 1.1935496e+00 3.91e-01 9.91e+02  -0.1 5.91e+02    -  2.10e-02 1.22e-03f  1
   9r 1.1935496e+00 3.91e-01 9.99e+02  -0.4 0.00e+00    -  0.00e+00 2.62e-07R  6
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10r 1.1179795e+00 1.95e-01 9.93e+02  -0.4 1.11e+02    -  1.71e-02 1.77e-03f  1
  11  1.1176704e+00 1.91e-01 1.80e+02  -1.0 1.68e-01    -  9.90e-01 2.37e-02h  1
  12  1.1176675e+00 1.90e-01 7.23e+05  -1.0 1.64e-01    -  9.90e-01 2.55e-04h  1
  13r 1.1176675e+00 1.90e-01 1.00e+03  -0.7 0.00e+00    -  0.00e+00 3.20e-07R  4
  14r 1.0799620e+00 9.62e-02 1.15e+03  -0.7 4.42e+01    -  1.17e-01 2.13e-03f  1
  15  1.0799636e+00 9.62e-02 2.23e+05  -1.0 8.24e-02    -  9.90e-01 6.92e-05h  1
  16  1.0799642e+00 9.62e-02 9.16e+09  -1.0 4.24e-02    -  9.90e-01 2.39e-05h  1
  17  1.1049964e+00 8.85e-04 9.15e+07  -1.0 8.12e-02    -  9.90e-01 1.00e+00h  1
  18  1.1049964e+00 8.85e-04 2.45e+12  -1.0 1.25e-03    -  9.90e-01 3.72e-05h  1
  19  1.1049963e+00 8.82e-04 2.54e+14  -1.0 6.29e-04    -  1.00e+00 9.48e-03h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  1.1049963e+00 8.82e-04 1.33e+11  -1.0 3.08e-07  14.0 9.99e-01 1.00e+00h  1
  21  1.1049963e+00 8.81e-04 3.62e+07  -1.0 1.09e-06  13.5 1.00e+00 1.00e+00f  1
  22  1.1049963e+00 8.75e-04 6.57e+07  -3.8 5.91e-06  13.0 1.00e+00 1.00e+00h  1
In iteration 22, 1 Slack too small, adjusting variable bound
  23  1.1049963e+00 8.75e-04 2.43e+22  -3.8 7.30e+01    -  1.00e+00 5.27e-11h  1
  24  1.1049963e+00 8.72e-04 3.65e+18  -3.8 3.34e-06  13.5 1.00e+00 1.00e+00h  1
  25  1.1049963e+00 3.98e-04 5.44e+14  -3.8 5.83e-04  13.0 1.00e+00 1.00e+00h  1
  26r 1.1049963e+00 3.98e-04 1.00e+03  -3.4 0.00e+00    -  0.00e+00 2.19e-07R  2
  27r 1.1051626e+00 1.48e-04 1.97e+03  -3.4 3.98e-01    -  4.00e-03 1.25e-03f  1
  28  1.1017412e+00 1.43e-04 2.91e+05  -3.8 1.00e-01    -  1.00e+00 3.42e-02h  1
  29r 1.1017412e+00 1.43e-04 1.00e+03  -3.8 0.00e+00    -  0.00e+00 3.74e-07R  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30r 1.1017407e+00 3.44e-05 1.95e+03  -3.8 1.31e-01    -  8.78e-03 1.08e-03f  1
  31  1.0048481e+00 2.00e+00 1.39e+02  -3.8 9.69e-02    -  1.00e+00 1.00e+00f  1
  32  1.4599240e-03 9.99e-01 4.54e+05  -3.8 4.00e+00    -  1.00e+00 5.01e-01f  1
  33  1.6267715e-03 9.98e-01 9.50e+05  -3.8 1.05e+00    -  2.54e-04 2.49e-03H  1
  34  1.0001504e+00 4.71e-10 3.87e-04  -3.8 2.00e+00    -  1.00e+00 1.00e+00H  1
  35  1.0001504e+00 1.11e-16 1.93e-06  -3.8 5.84e-09    -  1.00e+00 1.00e+00h  1
  36  9.9999999e-01 1.11e-16 2.81e+14  -8.6 1.50e-04    -  1.00e+00 1.00e+00f  1
WARNING: Problem in step computation; switching to emergency mode.
Cannot call restoration phase at point that is almost feasible (violation 1.110223e-16).
Abort in line search due to no other fall back.

Number of Iterations....: 36

                                   (scaled)                 (unscaled)
Objective...............:   9.9999999000037698e-01    9.9999999000037698e-01
Dual infeasibility......:   2.8147497034264009e+14    2.8147497034264009e+14
Constraint violation....:   1.1102230246251565e-16    1.1102230246251565e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.3819234855155923e-08    1.3819234855155923e-08
Overall NLP error.......:   2.8147497034264009e+14    2.8147497034264009e+14


Number of objective function evaluations             = 65
Number of objective gradient evaluations             = 37
Number of equality constraint evaluations            = 65
Number of inequality constraint evaluations          = 65
Number of equality constraint Jacobian evaluations   = 43
Number of inequality constraint Jacobian evaluations = 43
Number of Lagrangian Hessian evaluations             = 37
Total seconds in IPOPT                               = 0.025

EXIT: Error in step computation!
2

It looks like failure of constraint qualification, so the Lagrange multipliers blow up. This is consistent with what I’m seing in your model: constraints 2 and 3 are identical.

3

Hi @baroqueuse, welcome to the forum :smile:

Here’s how I would write your model. You’ll see that there are number of things that are different. As Charlie says, you have some repeated constraints, and there are also other issues with your algebraic expressions, like (... * u[2*i]) / (... * u[2*i]). In the JuMP side, sum(A + B) is just A + B, you can use linear algebra like A * x .== u, and it’s always better to provide explicit variable bounds.

This code finds your analytical solution:

using JuMP, Ipopt
A = [6.6 -3.0; -3.0 6.6]
nc = div(size(A, 2), 2)
tol, ub = 1e-6, 1e9
model = Model(Ipopt.Optimizer)
@variables(model, begin
    -1 <= x[i in 1:2nc] <= 1, (start = 1)
    tol <= u[i in 1:2nc] <= sqrt(ub)
    ratio[1:nc]
    t
end)
fix.(u[1:2:2nc], 0.0; force = true)  # [i in 1:nc], u[2i-1] == 0
@constraints(model, begin
    sum(x.^2) == 1
    [i in 1:nc], tol <= x[2i-1]^2 + x[2i]^2 <= 1
    A * x .== u
    [i in 1:nc], ratio[i] == x[2i] / sqrt(x[2i-1]^2 + x[2i]^2)
    [t; ratio] in MOI.NormInfinityCone(nc + 1)
end)
@objective(model, Min, t)
optimize!(model)
objective_value(model)
value.(x)
4

Even simpler

using JuMP, Ipopt
A = [6.6 -3.0; -3.0 6.6]
nc = div(size(A, 2), 2)
ub = maximum(sum(abs, A; dims = 2))
model = Model(Ipopt.Optimizer)
@variables(model, begin
    -1 <= x[i in 1:2nc] <= 1, (start = 1 / sqrt(nc))
    0 <= u[i in 1:2nc] <= (isodd(i) ? 0 : ub)
end)
@constraints(model, begin
    sum(x.^2) == 1
    A * x .== u
end)
@expression(model, ratio[i in 1:nc], x[2i] / sqrt(x[2i-1]^2 + x[2i]^2))
@objective(model, Min, maximum(abs, ratio))
optimize!(model)
objective_value(model)
value.(x)