checking inner variable in yalmip

[G]

Hello,

I am trying to solve a problem which involves a big-M relaxation. The model is quite complicated and at the end I get an error like

Warning: You have unbounded variables in IMPLIES leading to a lousy big-M relaxation.
...
...
Solver not applicable: sedumi

I guess (hope) I forgot some constraints regarding the variable inside the implies command, even if I really think having write everything.... There is an easy way to check the "inner status" of yalmip to understand which variable(s) are unbounded?

Thanks

[Johan Löfberg]

No simple way. If you supply some reproducible code, I can probably find it easily.

[G]

thanks I am more or less forced to give you the entire code

function uSol = MPCBigM(k0, betas, roads_o, roads_v, Tc, L, PhiIN, phi_max, v_max, w, x_max )

rho_c = 40;

x0 = zeros(12,1);

x0([1 3 5 7 9 11]) = [roads_o(1).rho(k0) ...
    roads_o(2).rho(k0) ...
    roads_o(3).rho(k0) ...
    roads_o(4).rho(k0) ...
    roads_o(5).rho(k0) ...
    roads_o(6).rho(k0)];

x0([2 4 6 8 10 12]) = [roads_v(1).rho(k0) ...
    roads_v(2).rho(k0) ...
    roads_v(3).rho(k0) ...
    roads_v(4).rho(k0) ...
    roads_v(5).rho(k0) ...
    roads_v(6).rho(k0)];

%%% receives x(0), betas, Tc

% T = 150 sec
% Tc = 15 sec
% #k = 10;

num = 150/15 -2;

xvar = sdpvar(12, num);
xvar(:,1) = x0;

constraints=[];

u = binvar(8, num); % for a network with 4 intersection 2x2, 12 roads

dem = sdpvar(12, num);

%constraints = [constraints dem>= 0];

supp = sdpvar(12, num);

%constraints = [constraints sup>= 0];

phi_o = sdpvar(12, num);

%constraints = [constraints phi_o<=phi_max phi_o>=0];

z = sdpvar(8, num);

%constraints = [constraints z<=phi_max z>=0];

cr = sdpvar(8,num);

alpha_sw = sdpvar(8,num);

obj = sdpvar(1);

for k = 2:num-1

constraints = [constraints xvar(:,k)<=x_max xvar(:,k)>=0 ...
    xvar(:,k-1)>=0 xvar(:,k-1)<=x_max dem(:,k-1)<=v_max*xvar(:,k-1) dem(:,k-1)<=phi_max dem(:,k-1)>=0];

constraints = [constraints supp(:,k-1)<=w*(x_max - xvar(:,k-1)) supp(:,k-1)<=phi_max supp(:,k-1) >= 0];

constraints = [constraints phi_o(:,k-1) <= dem(:,k-1) phi_o(:,k-1)>=0];

constraints = [constraints ...
    phi_o(1,k-1) <= supp(3,k-1)/betas(1,1) ...
    phi_o(1,k-1) <= supp(6,k-1)/betas(1,2) ...
    phi_o(2,k-1) <= supp(3,k-1)/betas(1,1) ...
    phi_o(2,k-1) <= supp(6,k-1)/betas(1,2) ...
    phi_o(3,k-1) <= supp(5,k-1)/betas(2,1) ...
    phi_o(3,k-1) <= supp(4,k-1)/betas(2,2) ...
    phi_o(6,k-1) <= supp(7,k-1)/betas(2,1) ...
    phi_o(6,k-1) <= supp(10,k-1)/betas(2,2) ...
    phi_o(8,k-1) <= supp(4,k-1)/betas(3,1) ...
    phi_o(8,k-1) <= supp(5,k-1)/betas(3,2) ...
    phi_o(9,k-1) <= supp(7,k-1)/betas(3,1) ...
    phi_o(9,k-1) <= supp(10,k-1)/betas(3,2) ...
    phi_o(11,k-1) <= supp(8,k-1)/betas(4,1) ...
    phi_o(11,k-1) <= supp(9,k-1)/betas(4,2) ...
    phi_o(12,k-1) <= supp(8,k-1)/betas(4,1) ...
    phi_o(12,k-1) <= supp(9,k-1)/betas(4,2)];

%phi_in = sdvpvar(12,1);

constraints = [constraints u(:,k-1)<=1 u(:,k-1)>=0 z(:,k-1)>=0 z(:,k-1)<= phi_max];
constr=[];
for i=1:8
    constr = [constr implies(u(i,k-1)==0, z(i,k-1)==0)];
end
    constraints = [constraints constr ...
        implies(u(1,k-1)==1, z(1,k-1)==phi_o(1,k-1)) ...
        implies(u(2,k-1)==1, z(2,k-1)==phi_o(2,k-1)) ...
        implies(u(3,k-1)==1, z(3,k-1)==phi_o(3,k-1)) ...
        implies(u(4,k-1)==1, z(4,k-1)==phi_o(6,k-1)) ...
        implies(u(5,k-1)==1, z(5,k-1)==phi_o(8,k-1)) ...
        implies(u(6,k-1)==1, z(6,k-1)==phi_o(9,k-1)) ...
        implies(u(7,k-1)==1, z(7,k-1)==phi_o(11,k-1)) ...
        implies(u(8,k-1)==1, z(8,k-1)==phi_o(12,k-1))];


constraints = [constraints alpha_sw(:,k-1) >= u(:,k) - u(:,k-1)];
constraints = [constraints alpha_sw(:,k-1) >= -(u(:,k) - u(:,k-1))];
cr(:,k) = cr(:,k-1) + alpha_sw(:,k-1);

constraints = [constraints alpha_sw(:,k-1)<=1 cr(:,k)<=2];

xvar(1,k) = xvar(1,k-1) + Tc/L*(PhiIN(1) - z(1,k-1));
xvar(2,k) = xvar(2,k-1) + Tc/L*(PhiIN(2) - z(2,k-1));
xvar(3,k) = xvar(3,k-1) + Tc/L*(beta(1,1)* (z(1,k-1) + z(2,k-1)) - z(3,k-1));
xvar(4,k) = xvar(4,k-1) + Tc/L*(beta(2,2)*(z(3,k-1) + z(5,k-1)) - phi_o(4,k-1));
xvar(5,k) = xvar(5,k-1) + Tc/L*(beta(2,1)*(z(3,k-1) + z(5,k-1)) - phi_o(5,k-1));
xvar(6,k) = xvar(6,k-1) + Tc/L*(beta(1,2)*(z(1,k-1) + z(2,k-1)) - z(4,k-1));
xvar(7,k) = xvar(7,k-1) + Tc/L*(beta(3,1)*(z(4,k-1) + z(6,k-1)) - phi_o(7,k-1));
xvar(8,k) = xvar(8,k-1) + Tc/L*(beta(4,2)*(z(7,k-1) + z(8,k-1)) - z(5,k-1));
xvar(9,k) = xvar(9,k-1) + Tc/L*(beta(4,1)*(z(7,k-1) + z(8,k-1)) - z(6,k-1));
xvar(10,k) = xvar(10,k-1) + Tc/L*(beta(3,2)*(z(6,k-1) + z(4,k-1)) - phi_o(10,k-1));
xvar(11,k) = xvar(11,k-1) + Tc/L*(PhiIN(3) - z(7,k-1));
xvar(12,k) = xvar(12,k-1) + Tc/L*(PhiIN(4) - z(8,k-1));

obj = obj + (xvar(:,k) - rho_c)'*(xvar(:,k) - rho_c);

end

ops = sdpsettings('solver','sedumi');
sol = optimize(constraints,obj,ops);
    if sol.problem == 0
        % Extract and display value
        uSol = value(u);
    else
        %display('Hmm, something went wrong!');
        sol.info
        yalmiperror(sol.problem)
    end

[G]

I don't know how to give you an input set for the function I wrote.....

[Johan Löfberg]

Just send me a zip file with all the data, and a ready-made command to run the model

load gdata
rungproblem

[Johan Löfberg]

Constraints of this type

 implies(u(i,k-1)==0, z(i,k-1)==0)

is extremly badly posed. Solvers work with finite precision of, say 1e-8. Hence, when the variables are around or smaller that tolerance, the implication is entirely undefined and you can get anything.

[Johan Löfberg]

Ah, u is binary. Then it is ok. However, you should write it as

constr = [constr implies(1-u(i,k-1), z(i,k-1)==0)];
implies(u(1,k-1), z(1,k-1)==phi_o(1,k-1))

as that better signals the logic variable in the first argument in the implies expression (I think YALMIP checks for this simple equality case though and exploits the structure, but I am not sure, hence there might be redundant variables introduced to model the truth of the equality, if you use the equality version, i.e,  implies(a==0, b==0) turns into implies(a==0,binary)+implies(binary,b==0) as yalmip flattens everything to expressions using implies(variable,condition) and implies(condition,variable))

[Johan Löfberg]

phi_o which is used in an implies is not explicitly bounded from above

[Johan Löfberg]

BTW, why are you using SeDuMi? I cannot see any conic in this model. Looks like a QP to me

[Johan Löfberg]

and of course, trying to use SeDuMi makes no sense as you have a mixed-integer program. Hence the solver not applicable message

[G]

yeah sorry for that silly line, it was there from before, i was just checking the boundness right now (btw which is the best solver for this case in your opinion?).

I think phi_o is bounded becouse I wrote

phi_o(:,k-1) <= dem(:,k-1) phi_o(:,k-1)>=0

and

[Johan Löfberg]

That is an implied bound that YALMIP effectively would have to solve an LP for to find. That would be way to expensive to do for every single variable (LP-based bound propagation). The variable has to be explicitly bounded.

[Johan Löfberg]

Gurobi, Mosek and cplex are all good mixed-integer solvers

[G]

Wow....now the warning it's not more shown.

Thank you very much for having looked at my code!!

P.S. Right now it's a good moment, if you may, to answer my question about the solver

[Johan Löfberg]

Note that if num happens to be 8 or 12, you might get symmetric variables. Use

xvar = sdpvar(12, num,'full');

etc

[G]

thanks for that trick as well!

What do you think about SCIP solver? I was quickly looking on the web and I found that as suggestion for free solver. Please note that I am an academic researcher but I cannot buy a solver!

[Johan Löfberg]

scip is a good solver, but note that gurobi, mosek and cplex all are free for academia

[G]

Done, thank you once again.

Another question: I get this error in the second time (in the same simulation) I try to solve the problem

Error using constraint/horzcat (line 11)
One of the constraints evaluates to a FALSE LOGICAL variable. Your model is infeasible

Is there a way to check which constraint is it?

[Johan Löfberg]

Just step through with debug until you find it. You are adding a constraint like 0*x >= 1

[Johan Löfberg]

and the error message should show the offendling line where it happened

function dgd
sdpvar x

f = [x >= 1];
f = [f, 0*x >= 1]

>> dgd

Error using constraint/horzcat (line 11)
One of the constraints evaluates to a FALSE LOGICAL variable. Your
model is infeasible

Error in dgd (line 5)
f = [f, 0*x >= 1]

[G]

but at each debug step how can i check the constraints ?

[Johan Löfberg]

You simply look at it and realize why it is trivially infeasible, or evaluate it

K>> 0*x >= 1

ans =

     0

[G]

Hi,

i restart this post because I have a related question.

Once again I have a quite complex model, representing a network. The constraints are created automatically (there is a piece of code that analyzes the adjacency matrix and adds constraints).

An instance of my problem result unfeasible, and this is not necessarily a mistake, but I need to check which constraints give the unfeasibility, to understand what happens at that intersection.

There exists a way to check the constraints in yalmip, to perform this kind of analysis?


Thanks

[Johan Löfberg]

There is no simple standard way of doing that, partly because the problem is hard to begin with

The following model does not have any single constraint which yields infeasibility

sdpvar x
F = [x >= 1, x <= -1]

but obviously having them both at the same time yields infeasibility

To search for bad constraints ad-hoc, simply solve the problem while removing some constraint

optimize(F) % Infeasible
optimize(F(2:end)) % Feasible. Aha, F(1) is problematic in the context of the other

or add slacks and minimize those

sdpvar s1 s2
F = [x >= 1-s1, x-s2 <= -1,s1>=0,s2>=0]
optimize(F,s1+s2)
>> s1
Linear scalar (real, 1 variable, current value : 2)
>> s2
Linear scalar (real, 1 variable, current value : 0)

aha, at least it works when relaxing the first constraint

[G]

ok but there is a way to print F in readable form while doing debug?

[Johan Löfberg]

If you mean a nice symbolic display of the constraints, typically no

however, maybe you are helped by tags
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.Tagging