Tolerance in BMIBNB solver

[Frank Zhu]

Hi,

I am using Yalmip to solve a complicated MPC problem. This is a non-linear non-convex problem, but BMIBNB most of the time gets pretty good solutions, though it takes some time. However sometimes it raises infeasible error. But the problem is definitely feasible. Then I assigned one known feasible solution and checked it. It turns out the constraints voilated are the equality constraints, which are the dynamics of my model, and only by a magnitude of ~1e-15 (See below). I think this is a computational problem and could be ignored, but I don't know how I could set the solver to increase its tolerance on equality constraints. Could you give me some hint on this? Thanks a lot!

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|     ID|                           Constraint|   Primal residual|   Dual residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|     #1|                  Equality constraint|                 0|             NaN|
|     #2|               Elementwise inequality|         2.623e-08|             NaN|
|     #3|               Elementwise inequality|           0.20241|             NaN|
|     #4|   Elementwise inequality (quadratic)|           0.46535|             NaN|
|     #5|   Elementwise inequality (quadratic)|           0.17905|             NaN|
|     #6|   Elementwise inequality (quadratic)|           12.8228|             NaN|
|     #7|   Elementwise inequality (quadratic)|            2.0782|             NaN|
|     #8|   Elementwise inequality (quadratic)|            8.3035|             NaN|
|     #9|   Elementwise inequality (quadratic)|           11.4088|             NaN|
|    #10|   Elementwise inequality (quadratic)|            4.4986|             NaN|
|    #11|   Elementwise inequality (quadratic)|            4.4281|             NaN|
|    #12|                  Equality constraint|       -1.3878e-17|             NaN|
|    #13|               Elementwise inequality|        3.2323e-08|             NaN|
|    #14|               Elementwise inequality|           0.37527|             NaN|
|    #15|   Elementwise inequality (quadratic)|           0.21573|             NaN|
|    #16|   Elementwise inequality (quadratic)|           0.17905|             NaN|
|    #17|   Elementwise inequality (quadratic)|           10.8426|             NaN|
|    #18|   Elementwise inequality (quadratic)|            1.1116|             NaN|
|    #19|   Elementwise inequality (quadratic)|            7.4401|             NaN|
|    #20|   Elementwise inequality (quadratic)|            9.2982|             NaN|
|    #21|   Elementwise inequality (quadratic)|            3.3686|             NaN|
|    #22|   Elementwise inequality (quadratic)|             2.893|             NaN|
|    #23|                  Equality constraint|                 0|             NaN|
|    #24|               Elementwise inequality|        4.1026e-08|             NaN|
|    #25|               Elementwise inequality|           0.45876|             NaN|
|    #26|   Elementwise inequality (quadratic)|          0.045608|             NaN|
|    #27|   Elementwise inequality (quadratic)|           0.17905|             NaN|
|    #28|   Elementwise inequality (quadratic)|            8.7156|             NaN|
|    #29|   Elementwise inequality (quadratic)|           0.35038|             NaN|
|    #30|   Elementwise inequality (quadratic)|            6.8927|             NaN|
|    #31|   Elementwise inequality (quadratic)|            7.0154|             NaN|
|    #32|   Elementwise inequality (quadratic)|             2.275|             NaN|
|    #33|   Elementwise inequality (quadratic)|            1.7173|             NaN|
|    #34|                  Equality constraint|       -8.6736e-19|             NaN|
|    #35|               Elementwise inequality|        5.4237e-08|             NaN|
|    #36|               Elementwise inequality|           0.48737|             NaN|
|    #37|   Elementwise inequality (quadratic)|         0.0034645|             NaN|
|    #38|   Elementwise inequality (quadratic)|           0.17905|             NaN|
|    #39|   Elementwise inequality (quadratic)|            6.8074|             NaN|
|    #40|   Elementwise inequality (quadratic)|          0.069253|             NaN|
|    #41|   Elementwise inequality (quadratic)|              6.72|             NaN|
|    #42|   Elementwise inequality (quadratic)|            4.9472|             NaN|
|    #43|   Elementwise inequality (quadratic)|            1.4527|             NaN|
|    #44|   Elementwise inequality (quadratic)|           0.99034|             NaN|

...

[Johan Löfberg]

Too much going in there to make any suggestion. You would have to supply a reproducible example with specified solvers

[Frank Zhu]

Hi,

You could try running the example here. If it's infeasible, it would check one known feasible solution and print out the Primal residuals. The one negative are equality constraints, but by a very small magnitude. The agents have been at the target places, so doing nothing will be a feasible solution. However, the solver I used, 'bmibnb', still gets infeasible solution. I just want to set the solver to increase its tolerance on theses equality constraints. 

Thanks.

[Johan Löfberg]

Which solvers are you using? It is not set in the bmibnb options, and you have not told us or shown a log where it is displayed

[Johan Löfberg]

and this is very far from a minimal example as it appears I have to run some simulation and then hope it will do exactly the same path and hitting the problematic state. Save the state where it fails and then crate a file which loads that state and calls the solver with everything setup

[Johan Löfberg]

or maybe that is precisely what you have done i see now

Anyway you have to tell us which solvers you are using. If you don't know, simply turn on verbose when you use optimize and it is listed by bmibnb

[Johan Löfberg]

well I tried to look at it and it just way to messy to follow

However, looking at the model you set up, it looks like something that easily could become compelete chaos in an optimizer setup as you will get super high dimension on the quadratics, which then has to be reduced once parameters are set. If here is any non-strictness in feasibility, numerical things will likely accumulate and cause issues

First step is to introducing extra parameters to avoid creating complex parametric expressions. Instead of having (parametera + parameter2)^2, simply call that expression parameter3 and send parameter3 as a parameter value (computed from the parameters a and b)

Instead of doing (x-parameter)^2, which thus will be complied as a quadratic in the variables (x,parameter), just introduce a new decision variable e, and use e^2 with the constraints e == x - parameter. Super simple sparse symbolic model, and trivially reduced once the parameter is fixed and no reduction of QP needed

[Frank Zhu]

Hi Johan,

Thanks a lot for your advice!

I turned on verbose and this is what I get:

* Starting YALMIP global branch & bound.
* Upper solver     : fmincon
* Lower solver     : QUADPROG
* LP solver        : LINPROG
* -Extracting bounds from model
* -Perfoming root-node bound propagation
* -Calling upper solver (no solution found)
* -Branch-variables : 144
* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation
* Timing: 7% spent in upper solver (1 problems solved)
*         0% spent in lower solver (0 problems solved)
*         0% spent in LP-based domain reduction (0 problems solved)
*         0% spent in upper heuristics (0 candidates tried)

Unfortunately, I tried your way of introducing extra parameters to avoid creating complex parametric expressions, but still got the same infeasible warning. May I know if there are other ways of solving this numerical problem? Could I solve it by somehow explicitly increasing the tolerance of the solver?

Thanksa great deal.

Best

[Johan Löfberg]

Ouch, solving nonconvex quadratic programs with 144 variables is basically destined to fail. If it ever works you are insanely lucky. To add to the misery, you don't have any efficient LP/QP solver

I dug out a minimal example which looks like

 load dummy
 check_val(ag,z, z_hat, mu, mu2, nu, u_star);

and probably illustrates your woes (dummy is the data required to make the call at the point of problems)

BTW, you have a bug in checkval as ver isn't declared, meaning it will be the output from the mtlab command ver, which will cause the call o optimize to crash

The problem is not related to optimizer at all, bmibnb simply thinks it is infeasible, which I will have to investigate

[Johan Löfberg]

and as far as I can see, a minimal example in check_val is

objective = 0;
for i = 1:N-1      
    objective = objective + param.alpha * (X{i}-param.Zdesired).'*(X{i}-param.Zdesired) + param.rho *(U{i}).'*U{i} ;
end
objective = objective + param.alpha * (X{end}-param.Zdesired).'*(X{end}-param.Zdesired);
constraints = [X{1} == z];

hence, you never tried to reduce to a minimal example

it appears something goes wrong in the presolve of quadratic objective

fredag 1 juli 2022 kl. 00:26:33 UTC+2 skrev Frank Zhu:

[Johan Löfberg]

which leads to the minimal example

x = sdpvar(10,1)

optimize([x(1:3) == 1e-4*randn(3,1)],x'*x,sdpsettings('solver','bmibnb'))

[Johan Löfberg]

If you uncomment line 62 in preprocess_bilinear_bounds at least the particular case runs. I think the problem is that you get a lot of stuff very close to zero, which is too aggressively presolved (and this is also why it is even possible to solve it, the problems are trivial as soo much is zero I guess)

>> load dummy
 
 check_val(ag,z, z_hat, mu, mu2, nu, u_star)
 

* Starting YALMIP global branch & bound.
* Upper solver     : fmincon

* Lower solver     : GUROBI
* LP solver        : GUROBI

* -Extracting bounds from model
* -Perfoming root-node bound propagation

* -Calling upper solver (found a solution, objective 1.2783e-08)
* -Branch-variables : 58

* -More root-node bound-propagation
* -Performing LP-based bound-propagation
* -And some more root-node bound-propagation

* Starting the b&b process
 Node       Upper       Gap(%)       Lower     Open   Time
    1 :   1.27831E-08     0.00    1.27831E-08    0     3s  Infeasible node in lower solver  
* Finished.  Cost: 1.2783e-08 (lower bound: 1.2783e-08, relative gap 0%)
* Termination with all nodes pruned
* Timing: 89% spent in upper solver (1 problems solved)

*         0% spent in lower solver (0 problems solved)

*         12% spent in LP-based domain reduction (108 problems solved)
*         1% spent in upper heuristics (0 candidates tried)

ans =

  struct with fields:

    yalmipversion: '20210331'
    matlabversion: '9.9.0.1524771 (R2020b) Update 2'
       yalmiptime: 0.179309343624414
       solvertime: 3.755690656375588
             info: 'Successfully solved (BMIBNB)'
          problem: 0

fredag 1 juli 2022 kl. 00:26:33 UTC+2 skrev Frank Zhu:

[Frank Zhu]

Hi,

But is there a way to solve this somehow? I just found out this kind of problem also appears when they have not yet been in the desired position, but is definitely feasible somehow... I understand this is a huge non convex problem which could be very difficult to solve. May I know if you have any suggestions for solving this kind of problem?

Many thanks!

[Johan Löfberg]

The way to solve it is to lure out the issues in bmibnb, which only can be done if you supply reproducible examples (in the form I did when I cleaned up your code and made a minimal example by stopping at the point at which the failure appeared, saved all relevant data, and then made an example from that, instead of forcing the full code to run to end up there)

[Frank Zhu]

Hi,

1. I tried to provide a minimal example including the needed constraints. I hope this could help solve the problem.

2. Besides, I have another question about how to specify using QUADPROG as the lower solver? Somtimes when I called optimize then it uses this solver and problem is feasible, but when I called optimizer it uses SeDuMi by default and the problem is infeasible... but by setting 'bmibnb.lowersolver' to 'QUADPROG, 'I got this error: "Failed exporting the model: No suitable solver". I don't know where it went wrong.

Many thanks

[Johan Löfberg]

I think this model in combination with optimizer+bmibnb simply won't work

% Works, yalmip understands that lower relaxations will be done in elementwise stuff and does not introduce socps for the convex stuff as linprog does not support it

x = sdpvar(2,1);
y = sdpvar(2,1);
ops = sdpsettings('solver','bmibnb','bmibnb.lowersolver','linprog')
optimize([-5>=[x,y] >= -5,x'*x >= 1,(x-y)'*(x-y) <= 1],sum(x),ops)

% Does not work, yalmip misses that solver will be relax the convex qps to elementwise stuff and models socp cones for the convex stuff, which linprog does not support

optimizer([-5>=[x,y] >= -5,x'*x >= 1,(x-y)'*(x-y) <= 1],sum(x),ops,y,x)

[Johan Löfberg]

Add

do_not_convert = do_not_convert | strcmpi(options.solver,'+bmibnb');

after line 220 in compileinterfacedata.m and solver selection in optimizer + bmibnb will be solved