Export Yalmip model as a Cplex Standard form in a correct way?

[Rui Li]

When I tried to export the Yalmip model to the cplex standard form, I found that the constraints lb <= x <= ub, has been recognized as the constraints Ax<=b, as a result, the lb, and ub in the exported model are empty. Is it possible to handle this issue?

Here is the code.

x = sdpvar(4,1,'full');
Q = [2 0.4 0.1 0; 0.4 4 3 -1; 0.1 3 6 1; 0 -1 1 10];
c = [0 0 0 0]';
A = -[7 8 10 14];
b = -[8.5];
Al = [1 1 1 1];
bl = 1;
v = [0 0 0 0]';
u = [1 1 1 1]';
obj = 0.5*x'*Q*x + c'*x;
F = [A*x <= b, Al*x == bl, v <= x <= u];

[model,recoverymodel] = export(F,obj,sdpsettings('solver','cplex'));

model.lb

model.ub

[Johan Löfberg]

no simple way to move those. YALMIP does not care about detecting this, as the solver will do that much faster during presolve

[Rui Li]

Hi all, 

I found a mistake that the optimal decision variable provided by the primal problem and its KKT system is not equal when I solve the primal problem and its KKT system with the matric form exported from Yalmip. Here is one example of this issue. Could you help me figure it out?

Ricky

-----------------------------------------------------

clc; clear all;
% prim-question

x = sdpvar(4,1,'full');
Q = [2 0.4 0.1 0; 0.4 4 3 -1; 0.1 3 6 1; 0 -1 1 10];
c = [0 0 0 0]';
A = -[7 8 10 14];
b = -[8.5];
Al = [1 1 1 1];
bl = 1;
v = [0 0 0 0]';
u = [1 1 1 1]';
obj = 0.5*x'*Q*x + c'*x;
F = [A*x <= b, Al*x == bl, v <= x <= u];

optimize(F,obj,sdpsettings('solver','cplex'))
prim_x = value(x) % primal optimal x

% KKT-system in matric form
[model,recoverymodel] = export(F,obj,sdpsettings('solver','cplex')); % matric form
x_kkt = sdpvar(4,1,'full')
miu =  sdpvar(size(model.Aineq,1),1,'full');
h = binvar(size(model.Aineq,1),1,'full');
lambda = sdpvar(size(model.Aeq,1),1,'full');
big_M = 100;

KKT = [];
KKT = [KKT, model.Aineq*x_kkt <= model.bineq];
KKT = [KKT, model.Aeq*x_kkt == model.beq];
KKT = [KKT, model.H*x_kkt + model.f == -model.Aineq'*miu + model.Aeq'*lambda];
KKT = [KKT, -model.Aineq*x_kkt + model.bineq <= big_M*h]; % linear complementary slackness constraint
KKT = [KKT, miu <= big_M*(1-h)]; %

optimize(KKT,[],sdpsettings('solver','cplex'))
kkt_x = value(x_kkt) % the optimal x from kkt system
prim_x - kkt_x  % compare the results from prim and its kkt system

[Johan Löfberg]

you've set up the wrong kkt system of some sort. the data is correct

optimize([model.Aineq*x_kkt <= model.bineq,model.Aeq*x_kkt == model.beq],0.5*x_kkt'*model.H*x_kkt+model.f'*x_kkt)

value(x_kkt)

note thare is a cmmand kkt in yalmip which can do all this for you

[Rui Li]

Thank you professor.

Actually, I found that the kkt command was effective for small-scale problems. When I proceed the kkt command on a large-scale problem, containing 48 binary variables and 5200 other variables, it fails to derive the KKT system with kkt command in 2-hour computation time.

[Rui Li]

Dear professor,

Here is a linear program example which has 3696 continuous variables. 

It takes a long time to derive the KKT system with the kkt command in Yalmip. Is there any approaches can speed-up the generation of the kkt system?

Attached files are the code and related data file.

Best,

Ricky

[Johan Löfberg]

So which LP solver are you using? Solving the LP once takes  around 0.05 seconds with mosek on my machine, so performing bound propagation of the whole linear portion of the KKT model will take in the order of 8000*0.05 seconds, which tells me you have a slow LP solver.

I hope you realize you never use KKT to simply derive the data matrices. KKT is used to derive an strong MILP representation of the KKT conditions, which involves boundp propagation etc. So in priciple, an LP with n variables and m inequalities will lead to KKT solving 2*(n+m) LPs. I also hope you understand that you never derive the KKT model to simply solve the LP. It is used when you have a nonconvex problem of some kind where the KKT conditions gives you a handle on solving the problem (nonconvex QP, bilevel LP etc) (derivation of the dual bounds can be turned off though as described  in the print out)

[Johan Löfberg]

Appears to be another performance bug in kkt which shows up for large problems like this. 

Anyhow, you probably don't want a strong MILP representation (because you will never be able to solve a problem involving this large KKT model anyway), so just turn off propagation of dual bounds and you have the model in 0s, or simply add dual bounds manually as explained in the kkt/bilevel examples

[Rui Li]

Thank you very much. 

Actually, I am solving a bi-level programming which has a large-scale linear problem at the lower level. Thus, I need to derive the KKT system of the lower problem and add them to the upper-level problem.

 

By the way, how can I turn off the propagations of dual bounds when I use the kkt command?     I found [K, details] = kkt(C,obj, 'dualbounds', 'off') is not corrrect.

[Johan Löfberg]

you set it to 0

and then you add constraints manually

https://yalmip.github.io/command/kkt/

[Rui Li]

'Thanks a lot.  Indeed, the solution time has been decreased significantly when I set the dual.bounds  as 0 and add the details.dual <= 1e3 manually.  However, the results are not satisfactory since the fact that the solved value(dual) equals to 0.

Besides, I have another question on the kkt command. As indicated in the bi-level example https://yalmip.github.io/command/kkt/,  [K,details] = kkt(CI,OI,[x1 x2]), the term [x1, x2] are parameters for the lower level and they are variables for the upper level. When the lower level objective has the form of x1*y1 + x2*y2 where y1 and y2 are variables for the lower level, the kkt command fails to derive the KKT system by showing 'KKT system can only be derived for LPs or QPs'.

My question is that the term x1*y1 + x2*y2 is linear for [y1, y2] since [x1, x2] are parameters for the lower level, why the command fail to give the kkt system for this level?

[Johan Löfberg]

so is 1e3 a correct bound? and has the problem actually been solved to optimality?

for the other question , you must supply reproducible code

[Ricky Lee]

In fact, 1e3 is not the correct bound. Besides, the result is not optimal. I am wondering how can I decide the bounds for details.dual manually? Could you offer any suggestions?

[Johan Löfberg]

If you are using a too small bound and thus cut away the true dual, and you're not solving the kkt condition (if that is what you mean with not optimal), then of course you cannot solve the problem

If you don't have any a-priori bound on the dual, it is very hard. Solve, check if dual hits bound. If it does, increase bound and solve again.