nested if-cycle using implies or possible problem relaxation?

[Luca]

Dear Johan ,

First of all, thank you for the great job you are doing with Yalmip.

 

I have an optimization problem that I am not sure how to pose, and if it can be solved in a linear or quadratic way using YALMIP + Gurobi. I write the problem without objective because the main issue is whether it is possible to write some constraints or not. 

I have 3 products (P1, P2, and P3). Their sum can be maximum Pmax e.g. Pmax=30. Assume their initial value is given e.g. P1=10, P2=10, P3=10. 

At every moment in time (e.g. from 1 to 10 ), I have to check how many products I have and:

 - if I have enough P1 (P1>=P1_take), I have to take a certain quantity of P1_take.

-  if I have not enough P1(P1<=(P1_take-small_number)), and if I have enough P2 (P2>=P2_take) I have to take the quantity of P2_take.

-  if I have not enough P1(P1<=(P1_take-small_number)), and if I have not enough P2 (P2<=(P2_take-small_number)), and if I have enough P3(P3>=P3_take)  I have to take the quantity of P3_take.

I remark that, for example, if I have enough P1 I have to take only P1_take, even if P2>=P2_take and P3>=P3_take.

Also, I know that the use of "small_number" does not make the problem very robust (every alternative is welcome).

I do not know how (and if it is possible) to write the above mentioned constraints. 

I worked out a code draft without the if-constraints, to show how I am approaching the problem, but I am not sure how to continue.

T=10; %Time steps max
P_max=30;
P1_take=2;
P2_take=1;
P3_take=0.5;

P1=sdpvar(T,1,'full');
P1_bounds=[(0 <= P1):'LB_P1', (P1 <=P_max):'UB_P1'] ;
P2=sdpvar(T,1,'full');
P2_bounds=[(0 <= P2):'LB_P2', (P2 <=P_max):'UB_P2'] ;
P3=sdpvar(T,1,'full');
P3_bounds=[(0 <= P3):'LB_P3', (P3 <=P_max):'UB_P3'] ;

Take1_bin=binvar(T,1); %Binary variable deciding if P1 is taken
Take2_bin=binvar(T,1); %Binary variable deciding if P2 is taken
Take3_bin=binvar(T,1); %Binary variable deciding if P3 is taken
Take=sdpvar(T,1,'full'); %variable determining the quantity taken
Take_bounds=[(0 <= Take):'LB_Take', (Take <=P1_take):'UB_Take'] ; %P1_take is the biggest quantity that can be taken, so it is the upper bound

Bounds=[P1_bounds,P2_bounds,P3_bounds,Take_bounds];

Cons0=[P1(1,1) == 10 , P2(1,1) == 10 ,P3(1,1) == 10 ]:'P_initial';
Cons1=[(P1 + P2 + P3) <= P_max.*ones(T,1)]:'Sum_P';
Cons2=[Take == Take1_bin.*P1_take + Take2_bin.*P2_take + Take3_bin.*P3_take ]:'Take';
Cons3=[(Take1_bin + Take2_bin + Take3_bin) <= ones(T,1) ]:'Sum_Take_bin'; %Maximum one quantity should be taken

Constraints=[Bounds, Cons0,...
    Cons1,...
    Cons2,...
    Cons3...
    ];
Options =sdpsettings('verbose',2,'solver','gurobi','debug',1);
Solution=optimize(Constraints,[],Options);

Thank you very much for your time. 

[Johan Löfberg]

If you already managed to model it, but end up with a bilinear model,

Take == Take1_bin.*P1_take + Take2_bin.*P2_take + Take3_bin.*P3_take

then all you have to do now is to linearize those binary*continuous expressions

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

Alternatively, do it directly by using clean logic disjunction Take1_bin implies P1>=P1_take and P1_take =..., Take2_bin implies P1<=P1_take and P2>P2_take and P2_take = ...

https://yalmip.github.io/modellingif

[Luca]

Thanks to your suggestion I managed to formulate it properly using logic disjunctions with implies(). The amount of constraints that arose, however, were too many for the amount of time-steps of the original problem (T>>>>>10). Basically it takes really a lot of time to "assemble" the problem i.e. the time between the moment when I start optimize() and the time when I see the first message of the solve appearing on Matlab command window using verbose=2.

This makes the problem solvable in theory, but it requires quite some time in practice. But I guess there is nothing to do about it. 

Thank you again. 

Best regards

Luca 

[Johan Löfberg]

Note that implies works with vectorized arguments.

The setup should really be very small compared to the time it takes to solve the problem, as this is a theoretically intractable problem that should take ages to solve

[Luca]

Dear Johan,

I tried to use implies with vectors but I get a different results. I updated my code with an objective that is basically trying NOT to respect the constraint. i.e. It is more convenient to do not respect the priority (P1->P2->P3).

If I state Cons5 as in the first way (for loop), the code works (objective=15). If I state it in the second way, it does not reach the expected solution (the objective in this case is 5.5).

Thank you. 

Best regards

T=10; %Time steps max
P_max=30;
P1_take=2;
P2_take=1;
P3_take=0.5;

P1=sdpvar(T,1,'full');
P1_bounds=[(0 <= P1):'LB_P1', (P1 <=P_max):'UB_P1'] ;
P2=sdpvar(T,1,'full');
P2_bounds=[(0 <= P2):'LB_P2', (P2 <=P_max):'UB_P2'] ;
P3=sdpvar(T,1,'full');
P3_bounds=[(0 <= P3):'LB_P3', (P3 <=P_max):'UB_P3'] ;

Take1_bin=binvar(T,1); %Binary variable deciding if P1 is taken
Take2_bin=binvar(T,1); %Binary variable deciding if P2 is taken
Take3_bin=binvar(T,1); %Binary variable deciding if P3 is taken

Bounds=[P1_bounds,P2_bounds,P3_bounds ];

Cons0=[P1(1,1) == 10 , P2(1,1) == 10 ,P3(1,1) == 10 ]:'P_initial';
Cons1=[(P1 + P2 + P3) <= P_max.*ones(T,1)]:'Sum_P';

Cons3=[(Take1_bin + Take2_bin + Take3_bin) == ones(T,1) ]:'Sum_Take_bin'; %Maximum one quantity should be taken

Cons4=[Cons4, P1(2:T) == P1(1:T-1) - Take1_bin(2:T)*P1_take,...
      P2(2:T) == P2(1:T-1) - Take2_bin(2:T)*P2_take,...
      P3(2:T) == P3(1:T-1) - Take3_bin(2:T)*P3_take ]:'P_amount';


Cons5=[];
for i=1:T
    Cons5=[Cons5,...
        implies(Take1_bin(i)==1,P1(i)>=P1_take),...
        implies(Take2_bin(i)==1,P2(i)>=P2_take),...
        implies(Take2_bin(i)==1,P1(i)<=P1_take),...
        implies(Take3_bin(i)==1,P3(i)>=P3_take),...
        implies(Take3_bin(i)==1,P2(i)<=P2_take),...
        implies(Take3_bin(i)==1,P1(i)<=P1_take) ];
end

% Cons5=[implies(Take1_bin == ones(T,1), P1 >= P1_take*ones(T,1)),...

%     implies(Take2_bin == ones(T,1),P2 >= P2_take*ones(T,1)), ...

%     implies(Take2_bin== ones(T,1),P1 <= P1_take*ones(T,1)), ...

%     implies(Take3_bin == ones(T,1),P3 >= P3_take*ones(T,1)), ...

%     implies(Take3_bin == ones(T,1) ,P2 <= P2_take*ones(T,1)), ...

%     implies(Take3_bin == ones(T,1),P1 <= P1_take*ones(T,1))];

Constraints=[Bounds, Cons0,Cons1, Cons3,Cons4,Cons5];
Options = sdpsettings('verbose',2,'solver','gurobi','debug',1);

Objective=0;
for i=1:T
    Objective=Objective + Take1_bin(i).*P1_take + Take2_bin(i).*P2_take + Take3_bin(i).*P3_take;
end

Solution=optimize(Constraints,Objective,Options);

Sol_P1=value(P1);
Sol_P2=value(P2);
Sol_P3=value(P3);
Sol_Take1_bin=value(Take1_bin);
Sol_Take2_bin=value(Take2_bin);
Sol_Take3_bin=value(Take3_bin);
Sol_Objective=value(Objective)

[Johan Löfberg]

The vectorization occurs for the case implies(Take1_bin, P1 >= P1_take*ones(T,1))

[Luca]

Great thanks, now it works and it is ~x3 faster