Nonlinear problem solved by different solvers

[kevin qiao]

Hi Prof. Johan.

    I'm solving a nonlinear problem with yalmip in Matlab, and different solvers, gurobi & ipopt, are used separately. However, the solution is quiet different. I use python+gurobi to checkout which solution is correct (as I think is correct), the solution from py+gurobi is the same as the one ipopt solves. I doubt why this happens, and I hope to get your help, thank you.

    The following is the code in Matlab.

clear
clc
tic
close all;
yalmip('clear')

%% ------------------------------- %%
%% Define a variable
%% ------------------------------- %%
x = sdpvar(5,1);
sdpvar Profit;

%% ------------------------------- %%
%% The constraints
%% ------------------------------- %%
Constraints = [
    x >= 0
    x(1)+x(2)^2-2 == 0;
    x(2)+2*x(3)^2-3 == 0;
    x(1)^2 - x(2) + x(3)^2 >= 0;
    x(1) + x(2)^2 + x(4)*x(5) <=20;
    x(4) == log(x(3));
    x(5) == x(3)*x(3);
];

%% ------------------------------- %%
%% My objective
%% ------------------------------- %%
Constraints = [Constraints,
    Profit == x(1)^2 + x(2)^2 + x(3)^2 + 8
];
Objective = Profit;

%% ------------------------------- %%
%% We do have many options
%% ------------------------------- %%
% Options = sdpsettings( 'verbose', 1, 'debug', 1 , 'solver', 'gurobi+', 'gurobi.NonConvex', '2', 'gurobi.MIPGap', '0.00000001');
Options = sdpsettings('solver', 'ipopt');

sol = optimize( Constraints, Objective , Options);

% Analyze error flags
if sol.problem == 0
    % Extract and display value
    display('========================================')
    display('The min objective is ')
    value(Objective)
    display('========================================')
    display('Thev variables are')
    value(x)
    display('========================================')
else
    display('Hmm, something went wrong!');
    sol.info
    yalmiperror(sol.problem)
end

[Johan Löfberg]

Gurobi could never solve this model, so I have no idea what you are talking about.

(adding a nonconvex equality to represent a convex quadratic objective is generally a bad idea (although it works ok here, indeed the solution found by ipopt can be confirmed to be the globally optimal one even))

[Johan Löfberg]

'gurobi.NonConvex' does not have to be touched by the user. YALMIP automatically invokes the nonconvex quadratic solver of gurobi when needed

lördag 16 oktober 2021 kl. 15:14:36 UTC+2 skrev wenji...@gmail.com:

[kevin qiao]

Gurobi indeed solve this model when I set sdpsettings('solver', 'gurobi+'), though the solution is not the optimal. While, as I change sdpsettings to sdpsettings('solver', 'gurobi'), the model can't be solved.

py+gurobi has the command like addGenConstrLog() to add log constraints, is there any command like this in YALMIP to add log constraints which can be solved by gurobi ? or I should linearize the nonlinear constraints like log func and x^n (n >= 2) ?

[Johan Löfberg]

Ah, with 'gurobi+' you re basically saying "I know what I am doing, I don't care if you think you cannot solve this, take the numerical data and just remove everything you think is wrong"

Gurobi does not support logarithms, and the high-level addGenConstrLog silently introduces a piecewise affine approximation to the model without your direct control. You can do that explicitly with with https://yalmip.github.io/command/interp1/. However, why not solve it using a global solver, bmibnb solves it trivially

[kevin qiao]

Thank you, I will try to solve it by a global solver.

[Johan Löfberg]

all you have to do is to change 'ipopt' to 'bmibnb'

[Johan Löfberg]

and if you want to use sos2 based interp1, i.e.  x(4) == interp1(xi,log(xi),x(3),'sos2');, you will have to install the develop branch as I found (and fixed) an issue in models combining nonconvex qudratics and sos2 constructs

söndag 17 oktober 2021 kl. 08:29:28 UTC+2 skrev wenji...@gmail.com:

[kevin qiao]

I use sos2 to linearize the log func, it seems well done. But it remains a problem how to determine the upper and lower bound. Is sos2 linearization method used by bmibnb in solving this kind of problems? If so, I just use bmibnb to solve considering simplify my model.

[kevin qiao]

plus, is there any method to directly call command in gurobi like addGenConstrLog() ? maybe the next version of YALMIP can add command to call special command in certain solver (it's just my thought)

[Johan Löfberg]

you have to derive bounds by simple bound propagation (either by using gurobi on a subset of equations, or by simply looking at the equations. it is easily seen from the first that x2<=sqrt(2)  from which it immediately follows from the second one that x3 >= sqrt((3-sqrt(2))/2) (and trivially x3<=sqrt(3/2), or something like that

bmibnb works in the continuous space, no approximations https://yalmip.github.io/tutorial/globaloptimization/

No, you have to model non-quadratic nonlinear operators using interp1, if you want to use gurobi for non-quadratic models

[kevin qiao]

hi Johan,

     I wonder that yalmip does not support the use of the symbol ‘...’ for line breaks ? my objective function is so long, I use the symbol, but the solution is quite different with not using the symbol.

[Johan Löfberg]

YALMIP is not involved in that parsing, so your matlab code in general is incorrect

[kevin qiao]

OK. Thank you!

[kevin qiao]

hello Johan

      I'm solving a nonlinear problem with an obj contains x.^5 ( x is a var matrix ) with linear constraints, and I want to solve it with gurobi, so I linearize the item by  sos2 based interp1, and my code is    lin_constrs = [lin_x == interp1(xi1, (xi1).^5, x, 'sos2')]; but the solution is something different with it solved by bmibnb (and with bmibnb solves it without linearization by hand), the following first solution is got by gurobi and the second is got by bmibnb. It can be seen that many solutions are similar, but the latter ones have a large gap, and I'm so confused what lead to this. btw, my obj contains y.^2 (y is a var matrix).

[Johan Löfberg]

You are approximating the problem, so why wouldn't you expect to see solutions which are slightly wrong?

[Johan Löfberg]

and I hope you mean x^.5 entering with a positive sign, as x^0.5 is concave and socp-representable, and thus no problems in gurobi when entering with negative sign

söndag 7 november 2021 kl. 09:00:17 UTC+1 skrev wenji...@gmail.com:

[kevin qiao]

well, it's x.^5, and I can accept there is a small gap between solutions got by the two solvers, however, the solutions in the back have a large gap

[kevin qiao]

I like solving problems with different solvers, because many problems I solve are nonlinear, to solve it, and to contrast solutions, I think in this way I can determine whether the solution found is available. sometimes it takes too long to solve a nonlinear problem, so I try to linearize it to get solutions with a faster time.

[Johan Löfberg]

All you can do then is to use more discretization points, there is nothing magic that will help you here (unless you simply have a bug of course)

If x is non-negative and x^5 enters positively, it is SOCP-representable via cpower

[kevin qiao]

okey, I'll try, thank you!

[kevin qiao]

hello Johan

       I'm solving a SOC problem contains binary constraints, and I want to extract the dual variable of my constraints, while mosek and gurobi seem not support this, so what can be done to extract the dual ?

[Johan Löfberg]

https://yalmip.github.io/faq/mipdual/

https://yalmip.github.io/faq/extractdual/

[kevin qiao]

hello Johan,

       I have a question about the dual variable, my model has a constrain whose variable is the dual of the other constrain, I can simply use dual() to extract the dual variable of the constrain and use the variable where I need? or what should I do? Thank you in advance!

[Johan Löfberg]

No, you will have to define the whole KKT system and work with dual variables that way. The command dual only returns the numerical value of a dual, after having solved a problem.

Effectively, you will have to solve a bilevel program, for which there is support on the LP/QP ide via kkt and solvebilevel.

[kevin qiao]

ok, I'll try, thank you

[kevin qiao]

hello Johan,

        I'm solving a two-stage RO problem, I wonder is there any lazy method that allows me to dualize the min problem of the second stage max-min subproblem, or extract the KKT condition of the min problem ? what I mean is that, I have the primal problem, and I want to transform it into the dual problem or the KKT condition  automatically.   Thank you in advance

[Johan Löfberg]

https://yalmip.github.io/tutorial/bilevelprogramming/

[kevin qiao]

aha, I forgot to solve it as bilevel optimization using this command, what if I want to solve max-min problem by using dual problem of min problem, thanks again!

[Johan Löfberg]

you have the linked kkt operator, that's it.

[kevin qiao]

OK， it seems I have to find the dual problem by hand, thank you, Dr. Johan

[Johan Löfberg]

why doesn't kkt work?

[kevin qiao]

it's a good question, the answer is kkt does work well, many papers, however, solve two-stage robust optimization using dual problem of the second stage min problem, with claim to be decoupled, I used to think kkt only can be used in LP (the min problem), I was wrong, of course. Now I think it's because the KKT condition is a necessary and insufficient condition for the optimal solution, while the dual problem is a sufficient and necessary condition for optimal solution. In fact, I'm confusd till now.

[kevin qiao]

hello Johan

      In my model, I use two classes of constraint name which are constrsA and constrsB to give the constraints, and two classes of obj name which are objA and objB to give my obj, and I set optimize([constrsA, constrsB], objA + objB) to optimize my model,  it runs and the solution is obtained, while, it's interesting that a variable named x, which is used in constrsA and constrsB both, satisfies constrsA but not satisfies constrsB (I checkout the solution of x and I find that constrsB seems to treat x as 0). and I wonder how to settle the problem. Thank you in advance.

[kevin qiao]

I have solved the problem, it's my fault.

[kevin qiao]

Hi Johan,

     I have a confusing question, my model is convex and nonlinear, it can be solved by ipopt with yalmip, and the solution of variable is marked as solVarIpopt and the obj of the problem as solObjIpopt, meantime, I use gurobi to solve the problem by linearize the model by SOS2 method, I take as many segments as possible as I think, it's lucky that the obj of the problem symboled as solObjGurobi is smaller than solObjIpopt (my problem is min-obj), is this mean that the quality of solution obtained by linearizing the model solved by gurobi better than the original nonlinear solved by ipopt (solVarGurobi is better than solVarIpopt)? cause I don't really know how ipopt solve nonlinear problem, I guess there are some linear method used by ipopt.

     Thank you in advance.

[Johan Löfberg]

you have not defined in what sense it is better and how sos2 is used. if an approximate objective is smaller, it could simply be that the approximation is a lower bound (tangent model as illustrated here https://yalmip.github.io/command/interp1/). if you mean the true objective is lower, it could be that the solution is not feasible due to outer approximation of feasible set, or that ipopt simply did not converge well

Why are you confused ipopt can solve the problem? ipopt is a solver for nonlinear programs

[kevin qiao]

Thank you for your reply, Johan, I know ipopt is a nonlinear slover, in my model, it is too time-consuming  for solving it, so I linearize it,  it becomes really faster to be solved, and I obtain a small obj, I think the quality of solution is better than the one I got from ipopt, because in my mind, ipopt use some kinds of method to linearize the nonlinear into linear constraints, I just replace the work it does. but I forget sos2 is just approximate to the original constraints, and I use 'check()' command to exam the numerical values of the primal and dual constraint violations, it makes me disappointed. as what can be seen below.

|   #2410|   Elementwise inequality|                 0|             NaN|
|   #2411|   Elementwise inequality|        16024.2028|             NaN|
|   #2412|   Elementwise inequality|         -1.93e-07|             NaN|
|   #2413|   Elementwise inequality|        17149.3904|             NaN|
|   #2414|   Elementwise inequality|                 0|             NaN|
|   #2415|   Elementwise inequality|        16264.7426|             NaN|
|   #2416|   Elementwise inequality|                 0|             NaN|
|   #2417|   Elementwise inequality|         16264.742|             NaN|
|   #2418|   Elementwise inequality|                 0|             NaN|
|   #2419|   Elementwise inequality|            1.1675|             NaN|
|   #2420|   Elementwise inequality|                 0|             NaN|
|   #2421|   Elementwise inequality|         14231.908|             NaN|
|   #2422|   Elementwise inequality|                 0|             NaN|
|   #2423|   Elementwise inequality|        19504.5174|             NaN|

however, my model contains y = x^4，z = ln(a) and some other kind of nonlinear constraints, I'll keep finding some methods to do with it. btw, ipopt don't solve nonlinear model with its approximation ?

[Johan Löfberg]

I don't understand what you re comparing or saying. 

If you re looking at the linearized constraints, then they are of course they are satisfied. You have to look at the objective and constraints of the original model, if you want to compare the solution you have obtained using your approximation, and the solution ipopt computes to the original model.

All nonlinear solvers work by iteratively approximating models, but the crucial part is that they do it iteratively and with theoretically justified steps and convergence properties meaning they (under various theoretical conditions) are guaranteed (sometimes...) to solve the original problem.

[kevin qiao]

got it, thank you again, Johan

[kevin qiao]

hello Johan,

       Does yalmip support LINGO(https://www.lindo.com/index.php/products/lingo-and-optimization-modeling) solver? btw, I find yalmip support LINDO(https://www.lindo.com/index.php/products/lindo-api-for-custom-optimization-application), I follow the step to install it, while it makes me sad cause I failed.

[Johan Löfberg]

No, obsolete remainings in current version, and then completely removed in develop branch. No users

[kevin qiao]

Hi Johan,

     I got a question and I have no idea what to do, so I long for your response. the obj is obj = x*ln(x) - x, I just use sos2 to linearize the total func. but the model is infeasible, and I try to linearize x*ln(x) as L1 and the obj transfers into obj = L1 - x, the model is infeasible, either. I wonder it's necessary to envelop x*ln(x) by McCormick envelop or not? because I think it can be linearized by sos2, while the truth is not, and I don't know why. Thanks  in advanced.

[Johan Löfberg]

How you model the objective can never influence feasibility, so you must be doing something wrong, such as only modelling it over a region where the problem is infeasible or something like that

Not sure why you talk about McCormick in this context. That is something which is used when solving the problem globally using general nonlinear global programming, and is done automatically 

[kevin qiao]

well, the solvers does reminder that the model is infeasible. I changed the obj func, it goes well, it means that my constraints are feasible, when I use sos2 method, the L1 = x*ln(x) is added as constraint, and the obj now is obj = L1 - x, so the problem must be this constraint when the model is feasible, can I take it this way ?

I think the obj = x*ln(x) - 1 is the bilinear problem ? if it does, the McCormick envelop is a method for me to solve it, by the way, I can take the locally feasible solution if there is some method to solve it.

[kevin qiao]

hi Johan

    I have addressed my problem. it seems that command sos2(variables) only acts on that variables belong to one certain original variable, if the original variable is a matrix, and I want to use the command, if will fail. for example, f(u_a) = ln(u_a) and I cut the original var u_a into m segments, and each additional variable corresponding to u_a^m is delta_a^m, if will act if I use for loop to code,  sos2(delta(:, a)), but fail when I code, sos2(delta(:, :)). Best wishes!

[Johan Löfberg]

I really don't see what you are talking about. x*ln(x) is not bilinear, it is nonlinear

McCormick is about over-approximating the set z = f(x), i.e. it does not give you *one* approximation of a function, but it allows the approximation to be in a set. It typically only works well when it is used in a b&b strategy where it is used over increasingly smaller domains guided by lower bounds obtained from these outer-approximations

Nothing prevents you from using a local solver directly with xlnx (or better yet, as you appear to minimize this, simply use -entropy(x) and perhaps you have the rest of the model expcone-representable too)