Constraint with a variable taking discrete values from a given set

[Monireh Mohebbi Moghaddam]

Hello everybody,

I am new to YALMIP. I need to solve an MILP optimization in attached file.

I post my question as a comment in YALMIP website, but I can't find it now.:-\

My problem has two decision variables, one of them is a binary ("alpha") and the other is a real variable ("u") in [0,1].

One of the Constrains is in a form like "e=q*u" ("q" is constant). In another Constraint, parameter "e" should take a discrete value from a given set. For example, "(e,b) ϵ

 {(e1,b1),(e2,b2),....,(en,bn)}". Optimization function is in a form like "b-e*k-alpha*p", where "k" and "p" are constant and "alpha" and "u" are desicion variables, as I mentioned.

Now, my question is how to define "e"? Which MATLAB function I use for "e" to take a discrete value from this set?

I will be very thankful if guide me.

[Johan Löfberg]

Haven't I already answered this?

To enforce [e;b] to take values from a given set, you use ismember

ismember([e;b],EB)

example

sdpvar e b

EB = randn(2,10);

optimize(ismember([e;b],EB),e^2+b^2)

Alternatively, you model it manually

sdpvar e b

EB = randn(2,10);
d = bnvar(10,1);

>> optimize([[e;b] == EB*d, sum(d)==1],e^2+b^2)

[Monireh Mohebbi Moghaddam]

Thank you so much.

Yes, you have answered but I couldn't find my comment and your answer! 

Sorry if my questions are basic.

Even when "e" is not an optimization decision variable, can I define it with sdvar? As the optimization variables are not "e" or "b". I thought sdvar just use for defining the optimization variables. Or, I should add "e" and "b" to the variables?

[Johan Löfberg]

How do you mean it is not a decision variable. If it is something the solver has to decide, it is a decision variable

unless it is a function of some other decision variable, then you don't define it using sdpvar, but simply create it from the other sdpvar variables

[Monireh Mohebbi Moghaddam]

Yes, "e" is not a decision variable. It is a function of a decision variable (e=q*u and "u" is a decision variable). But, in another constraint, "e" should take a discrete variable from a given set (as I explained in my first comment).

Hence, I have defined "u" by sdpvar and c1=[e==q*u] as first constraint. But, for the second constraint, is it true to define it in a way like this: c2=[ismemer([e;b],EB)], where EB={{1,2},{3,4}}?

Thank you.

[Johan Löfberg]

You either define e as a decision variable and add the constraint e == q*u, or you simply assign e = q*u. Does't really matter

EB in you case would be [[1;2] [3:4]], from the help

  F = ismember(x,P)

 

  If P is matrix of DOUBLE, a constraint constraining the vector x to equal

  one of the columns of P is created. This will introduce size(P,2) binary

  variables  

[Monireh Mohebbi Moghaddam]

Thank you so much again for your help :)

[Monireh Mohebbi Moghaddam]

I have another question about my problem. I would be very thankful if everyone help me.

I have a matrix "H" in my problem, for example with these values; H=[2,3,4;1,3,4;2,5,7], where H(i,j) represents number of hosts from a certain class "j" in a data center "i". Now, I need to define a real decision value for each host in each data center, for example 9 variables for row #1. I need to access to these variables in some constraints. 

I have defined optimization decision variables as follow:

H_i=sum(H,2);

H_class=sum(H,1);

H_s=sum(H_i);

for i=1:(#of rows in H)

        U{i}=sdpvar(H_i(i),1,'full');

end

Is it true? If it is not true, how should I define the variables?  For example, definition such as follow can be useful? U=sdpvar((#of rows in H),H_s)

And, how can I access to one of these variable? For example, U(1,2), i.e. host number 2 from data center 1 in  constraint (2).

[Johan Löfberg]

U{i}(j)

'full' is redundant as it is a vector

[Monireh Mohebbi Moghaddam]

Thank you so much for fast reply :)

[Monireh Mohebbi Moghaddam]

I came back with another question about constraint defining :|

I have used "ismember" function in my problem and it works good. Now, I want to use it for some variables, not for one, but I don't know how use it.

I have defined the following variables:

e=sdpvar(2,1);

b=sdpvar(2,1);

Now, I need to define one "ismember" constraint for each [e;b]. Something like this: C=[ismember([e(i);b(i)],EB(i))]. But, this form does not work! :(

How should I define them?

Thank you so much for your reply in advance.

[Johan Löfberg]

works here

>> e=sdpvar(2,1);b = sdpvar(2,1);EB=randn(2);ismember([e(1);b(1)],EB(1))

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

|   ID|                    Constraint|

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

|   #1|                    (ismember)|

|   #2|   Element-wise inequality 4x1|

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

[Monireh Mohebbi Moghaddam]

Should I define one constraint for each variable pair [e(i),b(i)]? I mean something like:

C1=ismember([e(1);b(1)],EB(1))

C2=ismember([e(2);b(2)],EB(2))

As when I use a loop as follow, it generates some error.

for i=1:2

 C=ismember([e(i);b(i)],EB(i))

end

[Johan Löfberg]

Runs here without "some error"

Although I guess you mean

for i=1:2

 C=[C, ismember([e(i);b(i)],EB(i))];

end

[Monireh Mohebbi Moghaddam]

This is my code and I receive the following error!

N_CP=3;

price=[20,40,60];

EB=[[0.01;20] [0.015;200] [0.02;20]];

e=sdpvar(1,N_CP);

b=sdpvar(1,N_CP);

%%

C=ismember([e(1);b(1)],EB);

for i=2:N_CP

    C=[C,ismember([e(i);b(i)],EB)];

end

%%

ep=zeros(N_CP,1);

for i=1:N_CP

        ep(i)=e(i)*price(i);

end

E=sum(ep,2);

B=sum(b,2);

obj=B-E;

sol=optimize(C,-obj);

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

Error using sdpvar/minus (line 54)

Adding NaN to an SDPVAR makes no sense.

Error in OPt_Test3 (line 22)

obj=B-E;

I did not understand what wrong is it? :\

[Johan Löfberg]

https://yalmip.github.io/naninmodel/

but the ep code doesn't even run as it is not valid matlab code

[Johan Löfberg]

did you edit the code?

anyway, you still have the nan in model due to reasons explained in the link

However, the code is unnecessarily complex, ep = e.*price...

[Monireh Mohebbi Moghaddam]

Thank you for your help.

[Monireh Mohebbi Moghaddam]

Hi,

Unfortunately, I faced to the another problem which is related to NaN! I wrote the following code:

    mig_cost=[3,1;1,0];

    H=[1,2,3;2,3,2];

    W=[2,2,2;1,1,2];

    o_h=[1,1,1,1,1,1,2,2,2,2,2,2,2,2];

    z=[1,1,1,1,1,1,2,2,2,2];

    H_s=13;

    W_s=10;

    %%

    alpha=binvar(H_s,W_s);

    alpha_sum=sum(alpha,1);

    C1=[alpha_sum==1];

    %%

   for i=1:H_s

        for j=1:W_s

            O=o_h(i);    

            Z=z(j);

            mig1(i,j)=alpha(i,j)*mig_cost(Z,O);            

        end

    end

    mig2=sum(mig1,2);

    Mig_total=sum(mig2);

    obj=-Mig_total; 

    sol=optimize(C1,-obj); 

    

     %%

    if sol.problem == 0 

        disp('Solver thinks it is feasible.');

        solution = value(sol);

        Opt_Value=value(obj);

    elseif sol.problem == 1

        disp('Solver thinks it is infeasible');

    else            

        display('Something else happened');

        sol.info;

    end

This code works well without any problem. But, when I change the value of "mig_cost(1,1)" to zero, "Opt_Value" changes to NaN, because of the value of "mig1". I check the solutions have been suggested in https://yalmip.github.io/naninmodel/, but the problem did not resolve :|

Could you help me, please?

[Johan Löfberg]

Works here

You have to show us the print-out from the solver and the diagnostic code etc

[Monireh Mohebbi Moghaddam]

Here the print screen of code and the solver.

Thank you

[Johan Löfberg]

but this is precisely the issue that is discussed in the link. You're creating mig1 by placing sdpvars into it something that isn't an sdpvar from start. In your second case, the first element you put in is a double, hence it is created as a double, and then when you try to put in other elements being sdpvars, you will get nans instead. In your first case, you are lucky as it is empty and you start by inserting an sdpvar, hence it will be an sdpvar

Create mig1 using concatenation instead (or better, vectorize the code)

[Monireh Mohebbi Moghaddam]

Thank you again.

I'm sorry for asking so many questions.

I have used "double2sdpvar", but it did not work. How can I use concatenation? As mig1 consists 1 value for each (i,j) and is not a vector! As for each (i,j) I should find the value of some function (O and Z).

Thank you so much

[Johan Löfberg]

mig1 = [];
 for i=1:H_s
     temp= [];

       for j=1:W_s

           O=o_h(i);    

           Z=z(j);

           temp = [temp alpha(i,j)*mig_cost(Z,O)];

       end
      mig1 = [mig1;temp]

   end

however, since you only use mig1 to construct mig2, you can just as well construct mig2 directly instead of creating mig1 first and then summing it. You can also vectorize the code above

[Monireh Mohebbi Moghaddam]

It works.

Thank you so much :)

[Monireh Mohebbi Moghaddam]

Hi again, 

I'm sorry for asking so many questions :\

I have asked a question about "ismemeer" function and Mr. Lofberg answered me. Now, I have another question which is related to this function.

In my problem, I have a matrix as EB(2,:) with some variables. For example, EB=[1,2,3;0,1,3]. I have two decision variables "e" and "b", which take their values from EB(1,:) and EB(2,:), respectively.

I have defined a constraint by "ismember" as follows:

C1=ismember([e;b],EB(:,:));

My optimization problem has another variable "en", which is not a decision variable. It is a function of other decision variable "q". (assume en=f(q)). Now, I want to define another constraint for "en"to link it to "e", in such a way "en<=e". But, I need the value of "e" is the smallest value which is greater than "en". 

As the objective function is something like "max (b-e)",constraint C2=[en<=e] does not generate correct value in some cases. For example, if "en=1.5", the correct value for "e" is 2, subsequently 1 for "b". But, in some cases, as by selecting "b=3", the optimization generates bigger value for "obj", it selects "e=3" not "e=2" form EB.

How can I change the constraint to solve this issue?

P.S. I try to use "find" function to find the index of "e" in EB to say that C2=[en<=e(i) and en>=e(i-1)] , but "find" is not defined for input arguments of type 'constraint'.

Thank you so much for your help.

[Johan Löfberg]

I don't see any clever way directly.

To model x = min(e(e>=en)), introduce a binary vector d to indicate which element you should pick out, add the constraint that sum(d)==1, implies(d(i), e(i)>=en), implies(d(i),x = e(i)), e(i) <= e(setdiff(1:n,i))

[Monireh Mohebbi Moghaddam]

Thank you for reply. But, I did not understand what I do clearly.

Do you mean that I keep the constraint C1=ismember([e;b],EB(:,:)), and define a binary vector "d" and add other constraints C2=[sum(d)==1] as well as C3=[e(i) <= e(setdiff(1:n,i))]?

How do I link "en" to "e"? And, how does "d" link to "e"?

[Johan Löfberg]

You text was confusing since you used the word en to denote both a threshold (if less than en) and to denote a variable.

Hence, if you want to model that x is equal to the smallest value of a vector e which is larger than a constant threshold z, i.e. x = min(e(e>=z)), you use the proposed model. Call the variables what ever you want. z, i.e. threshold, is a fixed number, e is a vector, and x i a new variable to represent the result of the function min(e(e>=z))

ismember has nothing to do with this

[Monireh Mohebbi Moghaddam]

Excuse me for unclear description.

"en" is a variable which is a function of some decision variable. But, as its value is independent from "e", we can simply assume it is a fixed number (threshold).

I have a vector "EB" with some describe values. "e" is a decision variable, which should take its value from "EB", i.e., e=sdpvar. And, C1=[ismember(e,EB)]

 Now, I want to model x = min(e(e>=en)), under C1 condition. 

[Johan Löfberg]

x = min(e(e>=en)) is modelled precisely as I said

What ever other stuff you have in your model is completely irrelevant (besides what is important in big-M modelling, i.e, e and en are explicitly bounded)

[Monireh Mohebbi Moghaddam]

Thank you. I truly sorry for many questions :(

I really would be thankful if you said me how I exactly define the constraints. 

I have EB=[1 2 3; 1 3 4] as a vector in my problem, and values of "e" and "b" should be chosen from this vector, i.e. ismember([e;b],EB(:,:)).

I should find something like e==min(EB(1,EB(1,:)>=en) ).

I model my problem as follows,but  I get an error.

b=sdpvar;

e=sdpvar;

en=1.5;

C1=ismember([e;b],EB(:,:));

C2=[e==min(EB(1,EB(1,:)>=en))];

C=[C1,C2]

obj=b-e;

sol=optimize(C,-obj);

Error using subsindex

Function 'subsindex' is not defined for values of class 'constraint'.

[Monireh Mohebbi Moghaddam]

This is my code and I have attached an print screen of the error.

"a,M,wl,p_idle,p_peak,gama" are vectors with fixed values and "EB" is a 3-D matrix.

w=intvar(3,3,'full'); 

e=sdpvar(1,3); 

b=sdpvar(3,1);

%********************* C1 *************************************************

w_sum1=sum(w,1);

C1=[w_sum1==wl];

%********************* C2 *************************************************    

w_sum2=sum(w,2);

C2=[w_sum2<=(a.*M)'];

%********************* C3 *************************************************   

C3=[w>=0];

%**************************************************************************

energy=[];

for j=1:3

    temp=[p_idle(j)*M(j)+(p_peak(j)-p_idle(j))*(w_sum2(j)/a(j))]*gama(j);

    energy=[energy;temp];

end 

%********************* C4 *************************************************

C4=ismember([e(1);b(1)],EB(:,:,1));

for i=2:3

    C4=[C4,ismember([e(i);b(i)],EB(:,:,i))];

end

%********************* C5 *************************************************

C5=[e(1)==min(EB(1,EB(1,:,1)>=energy(1),1))];

for i=2:3

   C5=[C5,e(i)==min(EB(1,EB(1,:,i)>=energy(i),i))];

end

%*********************************************************************

C=[C1,C2,C3,C4,C5];

obj=sum(b,1)-sum(e,2);

sol=optimize(C,-obj);

[Johan Löfberg]

You cannot index variables using constraints. That's precisely why we derive a big-M model

Here is the completely generic example to implement e = min(values(values>=en))

values = [1 2 3 4 5 6 7]';

n = length(values);

e = sdpvar(1);

en = sdpvar(1);

larger = binvar(n,1);

position = binvar(n,1);

Model = [e >= en, implies(values >= en,larger), sum(position) == 1];

for i = 1:n

    % If we pick out this position, then e should have this value

    Model = [Model, implies(position(i), e == values(i))];

    % If we pick out this solution, then e(i) must be the smallest

    % among those that are larger than en

    theother = setdiff(1:n,i);

    Model = [Model, implies(position(i) + larger(theother) >= 2, values(theother) >= e)];

end

Model = [Model, -10 <= [e;en] <= 10]   

% This should lead to position 4

optimize([Model,en == 3.5], e) 

value(position)

value(e)

% This should lead to position 7

optimize([Model,en == 6.5], e) 

value(position)

value(e)

Additional constraint on e or en, or how those are used elsewhere is not interesting. This is just a minimal example to show that the modeling works

[Monireh Mohebbi Moghaddam]

Thank you very much for your attention and help.

I am sorry for having disturbed you and asking many questions. I try to rewrite my code according to your example, but, I have faced some other problems:(

You have explained how I should implement min(values(values>=en)) by an example. Now, if “values” is a matrix, not a vector, and “b” should have the values of "values(2,i)" if we pick up position(i), how do I revise the code? Following code leads to incorrect value for “b”, but, I did not why.

Model=[Model, implies(position(i), b==values(2,i))];  

 

Additionally, the objective function in my problem is not min(values(values>=en)). This is one of the constraints. In other words, I should solve max(b-e*p) subject to: e=min(values(values>=en)) and some other constraints. How can I define it as a constraint, not as an objective function in optimization?

I change some parts of the code to better illustrate the issue. But, I did not know how to revise the code to generate the right result.

Thank you so much.

%======================================================

values = [1,3,9,9.5,12,18;2,5,7,9,3,1];

n = size(values,2);

e = sdpvar(1);

b = sdpvar(1);

en = sdpvar(1);

p=5;

 

larger = binvar(n,1);

position = binvar(n,1);

 

Model=[e>=en, implies(values(1,:)>=en,larger), sum(position)==1];

for i=1:n

    Model=[Model, implies(position(i), e==values(1,i];

    Model=[Model, implies(position(i), b==values(2,i))];  

    theother =setdiff(1:n,i);

    Model=[Model, implies(position(i)+larger(theother)>=2, values(1,theother)>=e)];

end

  

optimize([Model,en==11.9],-(b-e*p))

value(position)

value(e)

value(b)

 

 

[Johan Löfberg]

Your objective has nothing to do with the way we model min(e(e>=en)). If you want to use objective max(b-e*p), do so.

The generic model I have given is for one column. If you have many columns etc, you simply apply the same basic idea on every column, etc.

[Monireh Mohebbi Moghaddam]

Thank you.

I changed the example code to vectorize it. But, it leads to some error :-\ What is the problem?

values(:,:,1) = [1,3,9,9.5,12,18;2,5,7,9,3,1];

values(:,:,2) = [2,4,5,7.5,10,11;1,2.3,4,8,1,0.5];

n = 6;

N_S=4;

e=sdpvar(1,N_S); 

b=sdpvar(N_S,1); 

larger = binvar(N_S,n);

position = binvar(N_S,n);

en=11.9;

%======================================

    C4=[e(1)>=en,implies(values(1,:,1)>=en,larger(1,:)),sum(position(1,:))==1];

    for i=1:n

        C4=[C4,implies(position(1,i),e==values(1,i,1))];

        theother=setdiff(1:n,i);

        C4=[C4,implies(position(1,i)+larger(theother)>=2,values(1,theother,1)>=e(1))];

    end

%===============================================================

         for i=2:N_S

            C4=[e(i)>=en,implies(values(1,:,i)>=en,larger(i,:)),sum(position(i,:))==1];

            for j=1:n

                C4=[C4,implies(position(i,j),e==values(1,j,i))];

                theother=setdiff(1:n,j);

                C4=[C4,implies(position(i,j)+larger(theother)>=2,values(1,theother,i)>=e(i))];

            end

        end

%=================================================================

optimize(C4,e);

[Johan Löfberg]

You appear to don't know how nd-arrays work in MATLAB, as values(1,:,i) doesn't make sense for i=3 as you have now

you have vectors of different length, you should use cells instead

[Monireh Mohebbi Moghaddam]

Excuse me.

I made a mistake in the initializing parameter N_S. Actually, N_S=2 not 4. Matrixes "values" have the same number of rows and columns for all "i"s. Additionally, I should optimize (C4, min(e,2)).

It generates the following error:

Error using implies (line 80)

Size mismatch in input arguments

Error in OPt_Test (line 16)

    C4=[e(1)>=en,implies(EB(1,:,1)>=en,larger(1,:)),sum(position(1,:))==1]; 

[Monireh Mohebbi Moghaddam]

Excuse me again. I am so sorry for disturbing you.

I have attached a file, in which I have described a simple version of my optimization problem.

I need to model e=min(e(e>=en)) and use big-M for constraints C3 and C4. As we can see, the objective function does not include "e". Even when I include "e" in the objective function, and change it to max(b-e*p-e) or max(b-en-e), reported result is not correct, yet.

As an example, I wrote a code, in the attached file, by using of your example code . But, the values reported for "w", "en" and "e" have  not the correct values. I did not know how can I revise it.

I would be very thankful if help me.

[Johan Löfberg]

Your code does not run

>> N_S

N_S =

     4

>> a

a =

     4     3     2

>> w

Linear matrix variable 4x1 (full, real, integer, 4 variables)

>> 

[Monireh Mohebbi Moghaddam]

Excuse me,

I have forgotten to initialize some variables.

I have revised and attached the correct version. Now, it works. 

The value of "en" in this setting is 30. But, instead of generating value (32,10) for (e,b), it gives (35,161), because of the value of parameter "b" is bigger in this case.

[Johan Löfberg]

You still haven't fixed this

Warning: You have unbounded variables in IFF leading to a lousy big-M relaxation. 

I get the numbers you talk about though. Do you really have a real MILP solver installed?

K>> [values;value(position)']

ans =

    21    25    28    32    35    41

     5     8     9    10    61     1

     0     0     0     0     1     0

K>> value(e)

ans =

   35.0000

K>> value(b)

ans =

    61

[Monireh M]

I have installed MOSEK on my system. When I type "yalmiptest", YALMIP can find it.

Indeed, I have installed CPLEX. But, despite adding its path to MATLAB, YALMIP cannot find it.

Unfortunately, I don't know how can I fixed this problem :\

[Johan Löfberg]

Which versions? Yalmip r20170921 detects up to cplex 12.7.1

however, you shouldn't use cplex (unless you are on an old matlab installation) as they have bugs causing it to crash randomly in matlab 016/17

[Monireh M]

YALMIP version is 20160930. And, I have MATLAB R2013a on my system.

[Johan Löfberg]

hence you should update yalmip

[Monireh M]

Thank you so much.

If I update YALMIP and use CPLEX, will my issue caused to that warning be solved?

[Johan Löfberg]

the warning about unbounded variables in big-m has nothing to do with the solver. you have to add explicit bounds on all variables involved in big-m modelled expression

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

[Monireh M]

If I understood correctly, I should define some bounds on "en", "e" and "b" in my problem, as they are involved in "implies" constraints. Is it true?

[Johan Löfberg]

yes, and those bounds are trivially derived from your vectors values1 and values2

[Monireh M]

Hi,

Excuse me for many questions.

Unfortunately, even though I have bounded the variables, I still get the previous warning and the result is not correct:\

I have changed that constraint as follows and keep the remaining unchanged. 

a=[4,3,1];

values=[21,25,28,32,35,41;5,8,9,10,61,1];

values1=values(1,:);

values2=values(2,:);

n=size(values,2);

p=5;

N_S=3;

WL=10;

%=======================================================

e=sdpvar(1); 

b=sdpvar(1);

w=intvar(N_S,1); 

en=sdpvar(1);

en1=sdpvar(1,N_S);

larger=binvar(n,1);

position=binvar(n,1);

 

%% **********************************************************

C1=[w>=0];

C2=[sum(w,1)==WL];

 

for i=1:N_S

    en1(i)=a(i)*w(i);

end

en=sum(en1,2);

 

C3=[e>=en,21<=e<=41,5<=b<=61,20<=en<=80,implies(values1>=en,larger),sum(position)==1];

for i=1:n

    C3=[C3,implies(position(i),e==values1(i))];

    C3=[C3,implies(position(i),b==values2(i))];  

    theother =setdiff(1:n,i);

    C3=[C3,implies(position(i)+larger(theother)>=2,values1(theother)>=e)];

end

%% ************************************************************  

ep=e*pp;

obj=b-ep-e;

sol=optimize([C1,C2,C3],-obj);

[values;value(position)']

value(en)

value(w)

 

if sol.problem == 0 

    disp('*****Solver thinks it is feasible******');

    solution = value(sol);

    disp('***** Optimum value******');

    value(solution)

 

elseif sol.problem == 1

    disp('Solver thinks it is infeasible');

    Opt_Value=-9999999999999999999999;

else            

    display('Something else happened');

    sol.info;

    Opt_Value=-1000000000000000000000;

end

[Johan Löfberg]

Undefined function or variable 'pp'.

 

[Johan Löfberg]

You don't have any bound on w (en is not a variable, it is an expression derived from other variables. YALMIP cannot derive big-M bounds from general inequalities as that would be way too expensive)

[Monireh M]

Thank you. 

after I bounded "w", it works and the error disappeared.

But, the result is not correct when I put some variable except "e" in the objective function (for e.g. when I define obj=b-e). It seems there is no way to solve this issue :(

a=[4,3,1];

values=[9,25,28,32,35,41;5,8,9,10,61,1];

values1=values(1,:);

values2=values(2,:);

n=size(values,2);

N_S=3;

WL=10;

pp=5;

%% **********************************************************

e=sdpvar(1); 

b=sdpvar(1);

w=intvar(N_S,1); 

larger=binvar(n,1);

position=binvar(n,1);

 

%% **********************************************************

C1=[w>=0];

C2=[sum(w,1)==WL];

en=27;

 

C3=[e>=en,9<=e<=41,1<=b<=61,0<=w<=10,implies(values1>=en,larger),sum(position)==1];

for i=1:n

    C3=[C3,implies(position(i),e==values1(i))];

    C3=[C3,implies(position(i),b==values2(i))];  

    theother =setdiff(1:n,i);

    C3=[C3,implies(position(i)+larger(theother)>=2,values1(theother)>=e)];

end

%% ************************************************************  

obj=b-e;

sol=optimize([C1,C2,C3],-obj);

[values;value(position)']

value(obj)

value(w)

[Johan Löfberg]

You will have to define what "incorrect" means and what a correct result would be

[Monireh M]

For example, when en=27, "e" should be 28. But, after running the code, it gives value "35".

Only, when I define obj=e and sol==optimize([C1,C2,C3],obj), this result generates.

But, in my problem, objective function should be obj=b-e*pp, not obj=e.

[Johan Löfberg]

Show the code. Your old code does not run as pp isn't defined

[Monireh M]

a=[4,3,1];

values=[9,25,28,32,35,41;5,8,9,10,61,1];

values1=values(1,:);

values2=values(2,:);

n=size(values,2);

N_S=3;

WL=10;

pp=5;

%% **********************************************************

e=sdpvar(1); 

b=sdpvar(1);

w=intvar(N_S,1); 

%en=sdpvar(1);

%en1=sdpvar(1,N_S);

larger=binvar(n,1);

position=binvar(n,1);

 

%% **********************************************************

C1=[w>=0];

C2=[sum(w,1)==WL];

for i=1:N_S

    en1(i)=a(i)*w(i); %Equal to energy formulation

end

en=sum(en1,2);

%en=27;

 

C3=[e>=en,9<=e<=41,1<=b<=61,0<=w<=10,implies(values1>=en,larger),sum(position)==1];

for i=1:n

    C3=[C3,implies(position(i),e==values1(i))];

    C3=[C3,implies(position(i),b==values2(i))];  

    theother =setdiff(1:n,i);

    C3=[C3,implies(position(i)+larger(theother)>=2,values1(theother)>=e)];

end

%% ************************************************************  

ep=e*pp;

obj=b-ep-e;

[Johan Löfberg]

Well, I get this. If it is wrong, you have to explain why and how it is wrong

>> value(b)

ans =

     8

>> value(e)

ans =

    25

>> [value(position)';values]

ans =

     0     1     0     0     0     0

     9    25    28    32    35    41

     5     8     9    10    61     1

>> value(en)

ans =

    10

>> value(en1)

ans =

     0     0    10

[Johan Löfberg]

i.e., it appears to correctly model e = min(values1(values1 >= en)) as you initially wanted (and b = min(values2(values1 >= en)))

[Monireh M]

Yes! It seems that it works correctly as I want :))

Thank you so much.

Excuse me for disturbing you and asking many questions.

[Monireh M]

Please run the code with the following values for the parameters. It generates incorrect result :(

a=[4,3,2];

values=[40,45,55,60,80,91;35,70,9,90,10,70];

values1=values(1,:);

values2=values(2,:);

n=size(values,2);

N_S=3;

WL=20;

pp=0.005;

ans =

    40    45    55    60    80    91

    35    70     9    90    10    70

     0     0     0     1     0     0

ans =

   89.7000

ans =

    46

ans =

     0

     6

    14

[Johan Löfberg]

It is infeasible. You probably have other bounds etc

[Monireh M]

I have attached a print screen of the result.

a=[4,3,2];

values=[40,45,55,60,80,91;35,70,9,90,10,90];

values1=values(1,:);

values2=values(2,:);

n=size(values,2);

N_S=3;

WL=20;

pp=0.005;

%% **********************************************************

e=sdpvar(1); 

b=sdpvar(1);

w=intvar(N_S,1); 

%en=sdpvar(1);

%en1=sdpvar(1,N_S);

larger=binvar(n,1);

position=binvar(n,1);

 

%% **********************************************************

C1=[w>=0];

C2=[sum(w,1)==WL];

for i=1:N_S

    en1(i)=a(i)*w(i); %Equal to energy formulation

end

en=sum(en1,2);

 

C3=[e>=en,40<=e<=91,9<=b<=90,0<=w<=20,implies(values1>=en,larger),sum(position)==1];

for i=1:n

    C3=[C3,implies(position(i),e==values1(i))];

    C3=[C3,implies(position(i),b==values2(i))];  

    theother =setdiff(1:n,i);

    C3=[C3,implies(position(i)+larger(theother)>=2,values1(theother)>=e)];

end

%% ************************************************************  

ep=e*pp;

obj=b-ep-e;

sol=optimize([C1,C2,C3],-obj);

[values;value(position)']

value(obj)

value(en)

value(w)

 

if sol.problem == 0 % Extract and display value

[Monireh M]

Why does the following piece of code not work? It leads to an infeasible solution.

C3=[e>=en,40<=e<=91,9<=b<=90,0<=w<=20,sum(position)==1];

for i=1:n

    C3=[C3,implies(values1(i)>=en,position(i))];

    C3=[C3,implies(position(i),e==values1(i))];

    C3=[C3,implies(position(i),b==values2(i))];

end

[Johan Löfberg]

It might be that it should be 

    for j = theother

        C3=[C3,implies(position(i)+larger(j)>=2, e <= values1(j))];

    end

instead of

C3=[C3,implies(position(i)+larger(theother)>=2,values1(theother)>=e)];

You will have to go through the logic though

[Johan Löfberg]

Change WL to an sdpvar and minimize it, and you'll see tat the minimal value is 21, hence your value 20 is not possible

[Monireh M]

Excuse me, I have a primitive question about definition of C3, especially about implies function. I try to trace a simple case for testing C3.

Assume that "en=53", values1=[40,45,55,60,80,91] and values2=[35,70,9,10,50,90].

When optimizer runs and comes to "implies(values1>=en,larger)", "larger" takes [0,0,1,1,1,1]. Is it true?

After that, when the for-loop runs, just for i>=3, "implies(position(i)+larger(theother)>=2" will be true. For example, for i=3, this expresion is true for theother=4,5,6. And values1(theother)=[10,50,90]. Is it true? 

Now, my question is what does "values1(theother)>=e" express? It give "e"s that their values is less than values1(theother)? As value of values1(theother) is specific. 

Why we does not use == instead of >=? 

And for your post about WL, I should tell that WL is a constant.

As it is obvious from the optimization problem, sum(w) should be equal to WL, where "w"s are decision variables. Solver allocates some values to "w" parameters. By doing this, the value of "en" will be determined. "en" determination leads to a value for "e", which should be e=min(e(e>=en)), as well as for "b". The solver by changing the values of "w"s changes  "en", subsequently "e" and "b". This lead to different value for obj=b-e*pp.

Now, why should I change a constant to a variable and minimize it?

Thank you so much for your time.

[Johan Löfberg]

implies(values1>=en,larger) says that if values1(i)>=en then larger(i) = 1, hence larger = [0,0,1,1,1,1]

values1(theother)>=e is not used anymore as I think it was wrong.

Yes, WL is a constant, but the model was infeasible and I debugged by simply showing that if WL is a constant, it only works if that constant is 21 or larger (or the infeasibility comes from some other constraint) as that is required for feasibility. Once feasible, you can start talking about what you actually want to minimize. 

[Monireh M]

I am confused :\

I have just one constraint related to WL. Solver should choose some values for "w"s in such a way sum(w)=WL.

Did you just change the definition of WL to an intvar instead of a constant and change obj=WL, and got an infeasibility?

[Johan Löfberg]

NO, I just showed you a way to debug an infeasible model

If I give you the model

sdpvar x y
Objective = x^2+y^2
Model=[x==1,x>=100]

It will trivially be infeasible. To find out what the cause could be (trivial here of course) I could change it and try to solve

sdpvar x y WL
Objective = -WL
Model=[x==1,x>=WL]

this turns out to be feasible,  but the largest possible WL wil be 1, hence we now that the constant 100 simply isn't possible to have (if we know that the constant 1 in the first equality constraint is correct). It gives us a way to show that the current model simply is far from being feasible.

In your case you had picked a constant as 20, but using similiar strategy, I showed that the model simply isn't possible to solve then, as that constant has to be at least 21 for it to work. THe solver cannot pick any value of w such that it sums to 20 with the other constraints being satisfied

[Monireh M]

Hi,

I am really sorry for disturbing you.

I have vectorized my code. Whan, I run it, I face some errors. But, I could not understand what is the problem :\

I have attached my code. I would be really thankful if guide me.

[Johan Löfberg]

Vectorized implications have to be consistent in size. optimize(implies(EB(1,:,1)>=energy(1),larger(1,:))) for instance is 5x1 vs 1x5

Next time, send .m files, not docx...

[Monireh M]

Thank you for the reply.

I am really sorry for bothering you a lot.

I have revised the code and that problem solved. But, the vectorized version of the code (Test2.m) does not work correctly. some "Position" variables take NaN values!! :( 

I have attached two .m files, named Test2 and Test3.

Test2.m is vectorized version of Test3.m file.  

[Johan Löfberg]

To begin with, you create C4 on line 71, then you recreate it on line 88...

[Monireh M]

Oh, yes! Thank you.

 I forgot to add C4 in the constraint in line 88.

The previous problem has been solved. But, the problem is infeasible, now! :|

[Monireh M]

It seems that the problem becomes infeasible when C4 is vectorized! 

As when I hold C4 for "e(1)" and define e(2)=energy(2) and give some constant to "b(2)", the problem becomes feasible (something like the code in the attached file).

I did not know what does cause this issue? As C4 just told to the solver that select a pair (e,b) for other "e_i"s than e(1), as it does for e_1.

[Johan Löfberg]

https://yalmip.github.io/debugginginfeasible/

[Monireh M]

Hi,

I am so sorry  to be disturbing you by asking many questions.

I have changed the way that I defined the constraint which was related to discrete values (C4).

But, when I ran the attached code now, I have faced the following error:

Error using sdpvar/minus (line 44)

Adding NaN to an SDPVAR makes no sense.

It seems that there is a problem in defining "ep" and "b". I studied the reasons that cause NaN in the models in the following link (https://yalmip.github.io/naninmodel/). I thought that maybe it is because of the inserting the sdpvar in a double variable. But, even though I have used concatenation, I see this error.

Could you please guide me?

Thanks a lot.

[Johan Löfberg]

First, this line

LargeNo=1000000000;

should be

% I accept that I am completely killing any kind of reasonable numerics in my model

LargeNo=1000000000;

Then you claim you concatenate. No you don't. You still do b_prime(i,k) = ..., and to make mattes worse, your code looks obviously wrong as you try to put a vector in that element. Luckily, it is never a vector since you have even more bugs causing the right-hand side to never be a vector...

[Monireh M]

Hi,

I am so sorry for asking many questions.

I have revised my code. The final version has been attached.

It works now. But, the problem is about the execution time. When N_S=3, it runs in a time less than 1 minute, but for N=4, it takes several hours to run (sometimes more than 12 hours). I know that IP problems take a lot of time to run. But, is it usual to take this amount of time? I  have installed MATLAB 2017a and the last version of yalmip and use MOSEK as the solver.

Is it because of the problem model which is an IP or because of the used solver or because of the modeling or something else?

Do you have any suggestion?

Thank you

[Johan Löfberg]

It can take inf hours if you are unlucky. You have to try to tweak the model and come up with a stronger formulation and/or tweak various solver options