Question about five available positions choose three

[bin sun]

Hi, Johan,

I have a problem like what is shown in picture below. I have three cables, and five available positions. The current(I1,I2,I3) of each cable depends on the positions. My objective is to get the largest I1+I2+I3; The constraints are (c*(I.^2)')./d' <= 1(c and d are calculated by position). So I use below code to let each time two positions have I==0 and I get the largest I1+I2+I3 for different configurations.

c=[1,0.385740902327350,0.252086400923626,0.180051756213799,0.134338625773115;0.385740902327350,1,0.385740902327350,0.252086400923626,0.180051756213799;0.252086400923626,0.385740902327350,1,0.385740902327350,0.252086400923626;0.180051756213799,0.252086400923626,0.385740902327350,1,0.385740902327350;0.134338625773115,0.180051756213799,0.252086400923626,0.385740902327350,1];

d=[401933.568419812,401933.568419812,401933.568419812,401933.568419812,401933.568419812];

m=1

for i=1:n

    for j=1:n

        if j~=i

      I=intvar(1,5);

      intvar Ibase;

      desired = [ zeros(1,2) Ibase  Ibase Ibase];

      f=-sum(I);                      

      F=[I(i)==0]+[I(j)==0];   % each time two position have I=0, no cable

      F=F+[(c*(I.^2)')./d' <= 1]; 

      F=F+[sort(I)==desired];

      F=F+[0<=I<=600];   

      F=F+[1<=Ibase<=600];

      sol=optimize(F,f);

      I=value(I);

      Current(m,:)=I

     

        end

      m=m+1;

    end

end

I find that for different arrangements, I have different largest values. So I try to use below code to find the largest value and the smallest value shown in Figure.(516 and 476). So I use another decision variable pos.

C=[same as before]

D=[same as before]

I=intvar(1,5);

pos=binvar(1,5);

intvar Ibase;   

desired = [ zeros(1,2) Ibase  Ibase Ibase]; 

f=-sum(I);

 

F=[implies(pos,I== 0)];

F=F+[sum(pos)==2];   %so two I==0, no cables

 

 

F=F+[(c*(I.^2)')./d' <= 1];

F=F+[sort(I)==desired];    

F=F+[0<=I<=600];   

F=F+[1<=Ibase<=600];

 

sol=optimize(F,f);

I=value(I)  

pos=value(pos)

 

I find I only can get the value of 516. Could you tell me how can I change my code to get the value 476?

 

Thank you very much,

Best Regards,

Bin

[Johan Löfberg]

You will have to ask a much more specific YALMIP or modelling question, just not "here is my model I thought I would get some other solution"

[bin sun]

I use brute force method at first and I find that for different arrangements I have different largest current. 

So I want to automatically find these value.

So I introduce the decision variable pos, when pos==1, I==0, means no cable here. And my objective is -sum(I);

So I think I should get the value 516 and I exactly get it, which means the method is correct.

But when I try to get the value 476, I have no idea how to solve it. If I still use  -sum(I) it will still give me the largest value for best arrangement. But I prefer the largest value for worst arrangement.

I know you are very familiar with optimization. Do you have any idea about how to get the 476 value? Just an idea is enough. For example, should I introduce another decision variable?

Thank you very much,

Best,

Bin

[Johan Löfberg]

??

If you maximize, you will of course get the maximum solution. Why/how can you expect anything else? If you want something else, you will have to change the objective accordingly

[bin sun]

The problem is when I try to find 516, the final objective(find largest value among all configurations) and each step objective(find largest value for each configuration) is the same. Both the find the largest value.

But when I try to find the 476, the final objective(find smallest value among all configurations) and each step objective(find largest value for each configuration) is opposite.

[bin sun]

Because when I finish the maximize, the result is known to all. But the 476 value is nobody know and very useful for practice.

[Johan Löfberg]

How do you mean you try to find? At the moment, you maximize sum(f) (minimize -sum(f)) and thus you find the maximum solution.  What do you optmize to find the "476" solution, and what do you get instead?

[bin sun]

The 516 result is known to all. But the 476 value is nobody know and very useful for practice. So I want to use code to get this value. But I don't know how to find the 476 solution. 

The final objective(find smallest value among all configurations) and each step objective(find largest value for each configuration) is opposite. So I don't know how to write the code. 

I think three methods might be useful: 1. change objective to f=sum(I); but it give me -516. So it's wrong.

                                                          2. divide the code to two parts: 1. find find largest value for each configuration; 2. find the smallest value among step one. But have not idea how to do this.

                                                           3. add new decision variables

I think the method 2 should be best. But I can only think the method is brute force like the code at the beginning. But this is not optimization and it waste lots of time when I have more available positions.

[Johan Löfberg]

If you minimize sum(I), you cannot get the objective -516. You have a constraint I>=0, hence a trivial lower bound is 0.

[bin sun]

Yeah, that make sense. 

Do you have any idea about how to get the value 476? Like maybe use another decision variable will be useful? I really think about it the whole week but still have no idea until now.

Or

I can't have an optimization problem that the final objective(find smallest value among all configurations) and each step objective(find largest value for each configuration) is opposite. 

Thank you very much for your help,

Best,

Bin

[Johan Löfberg]

It is not clear what the problem is.

You wanted maximum so you maximized the objective and got what you wanted. Now you want to minimal solution instead, so what stops you from minimizing objective instead?

[bin sun]

Because in order to get 476, this optimization includes two steps.

The first step is to find the largest value for each configuration. which is finished by: f=-sum(I); and F=F+[(c*(I.^2)')./d' <= 1];. Using this objective function and this constraint, I can get the largest value for each configuration and even the largest value for all possible configurations.

But the second step is to find the smallest value from all these largest values I get from step one, which should be finished by f=sum(I); But this is not enough. 

This is where I feel confused.

[Johan Löfberg]

??

smallest = inf;

for i = 1:N

 solve optimization problem i

 if value < smallest

    smallest = value

 end

end

[bin sun]

So this part looks like a brute force method, right? 

for i = 1:N

 solve optimization problem i

 if value < smallest

    smallest = value

 end

 Just like the first code I list above: F=[I(i)==0]+[I(j)==0];

This is ok for only five available positions, but if I have 20 available positions, it will take days or weeks to finish N iterations. 

I do this for a 15 available positions and I want to choose 6 cables. The MATLAB calculates about 4500+ results but haven't finished.

Is there any method to calculate it faster? I will post my code below.

[Johan Löfberg]

your initial code you have looks strange though. Why use a for-loop to manually go through the all the combinatorial cases, instead of simply writing nnz(l)<= 3, or manually writing the big-M model  z = binvar(1,5)...l <= M*z, sum(z)<=3

[bin sun]

Thank you very much. This is very helpful. I will try to figure it out.

[Johan Löfberg]

aren't you then simply looking for the minimal value, i.e. minimize sum(f) subject to nnz(l)<=6.

[bin sun]

If I use nnz or big-M, the MATLAB will all give me only one result, but using the for-loop I can get all results and then I can use your code to find the smallest one.

for i = 1:N

 solve optimization problem i

 if value < smallest

    smallest = value

 end

end

But if I only get one final result, I can't use this loop code.

If I use sum(f) subject to nnz(I)<=6, Since I have a constraint:F=F+[1<=Ibase<=600];I will get 1,1,1; But I hope to find a method to get 476, which is  the smallest value from all these largest values I get from step one.

Do we have a command that can show all results like the for-loop? So I can get the smallest value like you mentioned before.

Thanks.

[bin sun]

Attached is my brute force method using for-loop for a 15 available positions. I will get all largest value for each configuration and I can use all these results find the smallest one among them.

c=[1,0.255860716455266,0.167208109795452,0.119427758544508,0.0891063841831220,0.215796124147757,0.196140596273289,0.159750264079153,0.127121232599849,0.101579754444553,0.147303771089326,0.142193424838124,0.129319129336975,0.113304068497061,0.0974536785386143;0.255860716455266,1,0.255860716455266,0.167208109795452,0.119427758544508,0.196140596273289,0.215796124147757,0.196140596273289,0.159750264079153,0.127121232599849,0.142193424838124,0.147303771089326,0.142193424838124,0.129319129336975,0.113304068497061;0.167208109795452,0.255860716455266,1,0.255860716455266,0.167208109795452,0.159750264079153,0.196140596273289,0.215796124147757,0.196140596273289,0.159750264079153,0.129319129336975,0.142193424838124,0.147303771089326,0.142193424838124,0.129319129336975;0.119427758544508,0.167208109795452,0.255860716455266,1,0.255860716455266,0.127121232599849,0.159750264079153,0.196140596273289,0.215796124147757,0.196140596273289,0.113304068497061,0.129319129336975,0.142193424838124,0.147303771089326,0.142193424838124;0.0891063841831220,0.119427758544508,0.167208109795452,0.255860716455266,1,0.101579754444553,0.127121232599849,0.159750264079153,0.196140596273289,0.215796124147757,0.0974536785386143,0.113304068497061,0.129319129336975,0.142193424838124,0.147303771089326;0.215796124147757,0.196140596273289,0.159750264079153,0.127121232599849,0.101579754444553,1,0.309401534633576,0.218425513520097,0.167208109795452,0.132807993278190,0.260910883769874,0.240787597022779,0.203052375128684,0.168358658633581,0.140245211353001;0.196140596273289,0.215796124147757,0.196140596273289,0.159750264079153,0.127121232599849,0.309401534633576,1,0.309401534633576,0.218425513520097,0.167208109795452,0.240787597022779,0.260910883769874,0.240787597022779,0.203052375128684,0.168358658633581;0.159750264079153,0.196140596273289,0.215796124147757,0.196140596273289,0.159750264079153,0.218425513520097,0.309401534633576,1,0.309401534633576,0.218425513520097,0.203052375128684,0.240787597022779,0.260910883769874,0.240787597022779,0.203052375128684;0.127121232599849,0.159750264079153,0.196140596273289,0.215796124147757,0.196140596273289,0.167208109795452,0.218425513520097,0.309401534633576,1,0.309401534633576,0.168358658633581,0.203052375128684,0.240787597022779,0.260910883769874,0.240787597022779;0.101579754444553,0.127121232599849,0.159750264079153,0.196140596273289,0.215796124147757,0.132807993278190,0.167208109795452,0.218425513520097,0.309401534633576,1,0.140245211353001,0.168358658633581,0.203052375128684,0.240787597022779,0.260910883769874;0.147303771089326,0.142193424838124,0.129319129336975,0.113304068497061,0.0974536785386143,0.260910883769874,0.240787597022779,0.203052375128684,0.168358658633581,0.140245211353001,1,0.347683397091418,0.255860716455266,0.203314404699171,0.167208109795452;0.142193424838124,0.147303771089326,0.142193424838124,0.129319129336975,0.113304068497061,0.240787597022779,0.260910883769874,0.240787597022779,0.203052375128684,0.168358658633581,0.347683397091418,1,0.347683397091418,0.255860716455266,0.203314404699171;0.129319129336975,0.142193424838124,0.147303771089326,0.142193424838124,0.129319129336975,0.203052375128684,0.240787597022779,0.260910883769874,0.240787597022779,0.203052375128684,0.255860716455266,0.347683397091418,1,0.347683397091418,0.255860716455266;0.113304068497061,0.129319129336975,0.142193424838124,0.147303771089326,0.142193424838124,0.168358658633581,0.203052375128684,0.240787597022779,0.260910883769874,0.240787597022779,0.203314404699171,0.255860716455266,0.347683397091418,1,0.347683397091418;0.0974536785386143,0.113304068497061,0.129319129336975,0.142193424838124,0.147303771089326,0.140245211353001,0.168358658633581,0.203052375128684,0.240787597022779,0.260910883769874,0.167208109795452,0.203314404699171,0.255860716455266,0.347683397091418,1]

d=[772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826,772899.067110826]

m=1

for pos1=1:n

    for pos2=1:n

        for pos3=1:n

            for pos4=1:n

                for pos5=1:n

                    for pos6=1:n

                       

        if pos1~=pos2~=pos3~=pos4~=pos5~=pos6

                          

      I=intvar(1,15);

      intvar Ibase; 

      desired = [ zeros(1,6) Ibase Ibase  Ibase Ibase Ibase Ibase Ibase Ibase  Ibase];

      f=-sum(I);                       % max (maximize a scalar function)

      

      

      F=[I(pos1)==0]+[I(pos2)==0]+[I(pos3)==0]+[I(pos4)==0]+[I(pos5)==0]+[I(pos6)==0];   % change 2 and 4 to specify the arrangement, which positions not choose.

      F=F+[(c*(I.^2)')./d' <= 1];  %key constraint of each arrangement

      F=F+[sort(I)==desired];

      F=F+[0<=I<=600];    

      F=F+[1<=Ibase<=600]; 

      sol=optimize(F,f);

      I=value(I);

      Current(m,:)=I

      

%         else

%             Curent(m,:)=[0,0,0,0,0]

                    end

      m=m+1;

                               

                                

                            end

                            

                        end

                        

                    end

                     

                end

                

            end

           

  

       end

This will get about 5000+ results and I try to find the smallest one among them. But this waste too many time So I hope maybe we have a command to save all results and then I can find the smallest  using :

smallest = inf;

for i = 1:N

 solve optimization problem i

 if value < smallest

    smallest = value

 end

end

[Johan Löfberg]

It feels incredibly confusing and unnecessarily complicated what you are doing.

intvar Ibase;

desired = [ zeros(1,2) Ibase  Ibase Ibase];

f=-sum(I);                      

F=[I(i)==0]+[I(j)==0];   % each time two position have I=0, no cable

F=F+[(c*(I.^2)')./d' <= 1]; 

F=F+[sort(I)==desired];

F=F+[0<=I<=600];   

F=F+[1<=Ibase<=600];

sol=optimize(F,f);

So effectively for fixed i and j , you only have one decision variable Ibase, which means that (c*(I.^2)')./d' will simplify to Ibase^2*(c*(q.^2)')./d'<=1 for some 0/1 vector q. Since sum(f) simply will be 3*Ibase, you are just maximizing Ibase. Hence, the optimal Ibase will be sqrt(1/(c*(q.^2)')./d'), and you're done.

[bin sun]

Yes. So I can get largest value of Ibase for every fixed i and j. 

But my problem is to find the smallest value from all possible largest values for each fixed configuration.

So I need to have two steps:

1. get each largest value of Ibase for every pair of i and j. 

2. find the smallest value from step 1.

But now the problem is how to get the smallest value from step one (Or how to get all values of step one so I can use a loop to find the smallest one among them)

[Johan Löfberg]

Basically the inner problem boils down to something like Ibase^2*q <= 1 where q is a sum of all elements but 2 in a given data vector, best possible Ibase^2 is 1/q, and the worst Ibase^2 you can get is is if q is large, meaning you select the n-2 largest elements. What I mean is that if you think through the model, I would not be surprised if you can solve it explicitly something like that

Generally speaking, you have what is called a bilvel problem,

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

which YALMIP has partial support for. However, YALMIP only solves bilevel problem where the inner prolem is a QP or LP, while the outer problem can be a MILP etc. Hence, if you can formulate the inner problem as an LP/QP, and the outer selection problem as a MILP, you're home. If you cannot, you have a much harder problem. But a I said, I think you can solve this explicitly