Using Implies in deciding the range of a loop

[Yanyan]

Hi Johan,

I am using implies to decide the range of loops. The main purpose of the following piece of code is like this. Given a matrix a and xx, I want to find VA with minimal sum(VA). In detail, (1). For the outer loop (i.e., for n=1:N), I need to find the best stopping point of n, which is decided by CA(n). If CA(n)==1, then n==ceil(bb)), VA(N-n+1:N)==Va(1:n) and VA(1:N-n)==0. (2). For the inner loop (i.e., for m=1:M), I need to find the best stopping point of m, which is decided by Ca(n,m)==1, then Va(n)==tem(n,2), where tem(n,2)==sum(a(1:m)*xx(1)*sum(a(1:m)*xx(2), when Ca(n,m)==1, as annotated in the code. But after running the code, it always shows the infeasible solutions, i.e., VA is all zero and binary variables are decimals, etc. However, this problem indeed has solution, which can be derived by hand. Could you please look into it and let me know whether there is something wrong in the following code? 

Thank you so much!

Yan

N=3;M=4;xx=[1,-2];a=[-1 5 -2 1];VA=sdpvar(1,N);Va=sdpvar(1,N);temb=binvar(1,N);tem=sdpvar(N,2);

Ca=binvar(N,M,'full');va=sdpvar(2,N,'full');CA=binvar(1,N);bb=sdpvar(1);

Model=[sum(CA)==1,binary(CA),binary(Ca),0<bb<=2,0<Va<50];

for n=1:N

    Model=[Model,sum(Ca(n,:))==1];

    for l=1:2

   if l==1     

    for m=1:M

     Model=[Model,implies(Ca(n,m),temb(n)==1)];

     Model=[Model,implies(Ca(n,m),tem(n,l)==sum(a(1:m))*xx(l))];

    end

   else

       for m=1:M

 Model=[Model,implies(Ca(n,m),tem(n,l)==tem(n,l-1)*sum(a(1:m))*xx(l))];

       end

   end

end

Model=[Model,implies(temb(n),Va(n)==tem(n,2))];%tem(n,2)=sum(a(1:m)*xx(1)*sum(a(1:m)*xx(2), when Ca(n,m)==1

Model=[Model,implies(CA(n),n==ceil(bb))]; Model=[Model,implies(CA(n),VA(N-n+1:N)==Va(1:n))];

 if n<N

   Model=[Model,implies(CA(n),VA(1:N-n)==0)];

 end

end

options=sdpsettings('solver','bnb','verbose',0,'debug',1);%Model=[Model,VA(N)>0];

opt=solvesdp(Model,sum(VA),options);

VA=double(VA)

Va=double(Va)

CA=double(CA)

Ca=double(Ca)

 

[Johan Löfberg]

To begin with, some oddities in your code

1. You define variables as binary, but then add binary constraint also. That is redundant

2. You use both strict and non-strict inequalities. Only use non-strict. You should see warnings about this

3. Bound non-binary variables involved in implies operators. You should see warnings about this

4. The solution you get when it is declared as infeasible has no meaning. It is infeasible, no solution exist

Even if I coment out some constraints, it is still infeasible. To try to understand what is inconsistent, try to relax constraints etc. Here I use some slacks to see how I can relax and make it feasible, and it appears the m=1 case is impossible

N=3;M=4;xx=[1,-2];a=[-1 5 -2 1];

VA=sdpvar(1,N);

Va=sdpvar(1,N);

temb=binvar(1,N);tem=sdpvar(N,2);

Ca=binvar(N,M,'full');

va=sdpvar(2,N,'full');CA=binvar(1,N);bb=sdpvar(1);

slack = sdpvar(N,2,M);

Model=[sum(CA)==1,binary(CA),binary(Ca),0<=bb<=2,0<=Va<=50, -1000 <= tem <= 1000,0<= slack <= 1000 ];

for n=1:N

    Model=[Model,sum(Ca(n,:))==1];

    for l=1:2

        if l==1

            for m=1:M

                Model=[Model,implies(Ca(n,m),temb(n)==1)];

                Model=[Model,implies(Ca(n,m),tem(n,l)==sum(a(1:m))*xx(l)+slack(n,l,m)) ];

            end

        else

            for m=1:M

                Model=[Model,implies(Ca(n,m),tem(n,l)==tem(n,l-1)*sum(a(1:m))*xx(l)+slack(n,l,m))];

            end

        end

    end

    Model=[Model,implies(temb(n),Va(n)==tem(n,2))];

   % Model=[Model,implies(CA(n),n==ceil(bb))]; Model=[Model,implies(CA(n),VA(N-n+1:N)==Va(1:n))];

   % if n<N

   %     Model=[Model,implies(CA(n),VA(1:N-n)==0)];

   % end

end

opt=solvesdp(Model,max(slack(:)))

value(slack)

[Johan Löfberg]

andf you should install a MILP solver so you don't have to resort to bnb.

[Mark L. Stone]

Have you tried to assign your hand "solution", and use checkset to check for feasibility against your constraints as entered?

[Yanyan]

The hand solution can be this:   CA=[1 0 0];Ca=[0 1 0 0;0 1 0 0;0 1 0 0];bb=0.6;Va=[-32 -32 -32];VA=[0 0 -32]; But could you specify how to use the checkset? I have tried to use checkset(Model), but it gives me the error: Conversion to double from cell is not possible. I also tried to use this:  value(CA)=[1 0 0];value(Ca)=[0 1 0 0;0 1 0 0;0 1 0 0];value(bb)=0.6;value(Va)=[-32 -32 -32];value(VA)=[0 0 -32]; or double(CA)=[1 0 0];, both of this assignment is not allowed. 

I would highly appreciate if you can give me some hint.

Thanks a lot,

[Yanyan]

Hi Johan,

I have installed mosek and used the non-strict inequalities, unfortunately, still the problem looks infeasible. After running your posted code, how do you conclude that when m==1, it is impossible? I could not figure it out. 

Also, if I have hand solution, such as: CA=[1 0 0];Ca=[0 1 0 0;0 1 0 0;0 1 0 0];bb=0.6;Va=[-32 -32 -32];VA=[0 0 -32]; Could you specify how to use the checkset?

Thanks,

Yan

[Johan Löfberg]

Yes it is infeasible, I'm just saying you should have a good solver. The proposed solution is trviially infeasible (Va!)

[Yanyan]

Thanks for your quick response.

How about my second question?  I.e., if I have hand solution, such as: CA=[1 0 0];Ca=[0 1 0 0;0 1 0 0;0 1 0 0];bb=0.6;Va=[-32 -32 -32];VA=[0 0 -32]; Could you specify how to use the checkset?

Thanks,

[Johan Löfberg]

You cannot use check (previously called checkset) when you have big-M represented models (such as implies). You simply have to check the constraints manually, which immediately reveals that Va is wrong.

[Yanyan]

Sorry to reply so late.

I find that whenever we have Model=[Model,implies(Ca(n,m),tem(n,l)==tem(n,l-1)*sum(a(1:m))*xx(l))]; the problem is infeasible. And once I remove this command, everything is fine. Also, I tested if we change it to  Model=[Model,implies(Ca(n,m),tem(n,l)==sum(a(1:m))*xx(l))]; the problem is also feasible. I am wondering, there must be something wrong with tem(n,l)==tem(n,l-1)*sum(a(1:m))*xx(l)) in implies. But why? I just want to use the multiplication of tem(n,1),tem(n,2),…,tem(n,l-1). If not put in this way, how could I realize this?

Thank you so much!

[Johan Löfberg]

what about the trivial error I told you about in your known solution? you claimed you knew a solution, which is it?

[Yanyan]

You didn't talk about the error in my known solution. Based on the current simplified code, a hand-solution is: Ca=[1 0 0 0;1 0 0 0;1 0 0 0]; tem=[-1 -2; -1 -2; -1 -2]; temb=[1 1 1]; Va=[-2 -2 -2]; The code we are testing now is:

N=3;M=4;xx=[1,-2];a=[-1 5 -2 1];

VA=sdpvar(1,N);

Va=sdpvar(1,N);

temb=binvar(1,N);tem=sdpvar(N,2);

Ca=binvar(N,M,'full');

va=sdpvar(2,N,'full');CA=binvar(1,N);bb=sdpvar(1);

% slack = sdpvar(N,2,M);

Model=[sum(CA)==1,0<=bb<=2,0<=Va<=50, -1000 <= tem <= 1000];%,0<= slack <= 1000 ];

for n=1:N

    Model=[Model,sum(Ca(n,:))==1];

    for l=1:2

        if l==1

            for m=1:M

                Model=[Model,implies(Ca(n,m),temb(n)==1)];

                Model=[Model,implies(Ca(n,m),tem(n,l)==sum(a(1:m))*xx(l))];

            end

        else

            for m=1:M

                Model=[Model,implies(Ca(n,m),tem(n,l)==sum(a(1:m))*xx(l)*tem(n,l-1))];

            end

        end

    end

    Model=[Model,implies(temb(n),Va(n)==tem(n,2))];

end

options=sdpsettings('solver','mosek','verbose',0,'debug',1);

opt=solvesdp(Model,-sum(Va),options)

value(Va)

[Johan Löfberg]

I've already told you that Va is wrong

[Yanyan]

Alright, but I don’t see why it is wrong. Could you please point out which constraint indicates it is wrong?

Thanks,

[Johan Löfberg]

It's negative and your model says it is positive

[Yanyan]

That is the problem! I am so embarrassed by the mistake I made. Now it is solvable. My apologies! 

BTW, if in implies, I want to use two conditions, i.e., binary variables to generate a constraint. For example, implies(a,b, c==1) (which is definitely wrong here); where a,b and c are all binary variables. So can I use: implies (a,b==1); implies(b,c==1); to achieve the function I want? 

Thank you!

[Johan Löfberg]

depends on what you actually mean with implies(a,b,c==1)

either

implies(a &&b,c)

or

implies(a| b,c)

which is equivalent to

implies(a,c)+implies(b,c)

[Yanyan]

I don't think so. What I want is: only if both a==1 and b==1, then c=2, where a,b are all binvar and c is spdvar. When I use implies(a&&b, c==1); it shows me the error: Conversion to logical from sdpvar is not possible. When I use implies(a,b); implies(b,c==2), it gives me: infeasible problem under both mosek and bnb. 

Thanks a lot!

[Johan Löfberg]

Sorry, typo,

>> implies(a & b,c==2)

+++++++++++++++++++++

|   ID|   Constraint|

+++++++++++++++++++++

|   #1|    (implies)|

+++++++++++++++++++++