Always busy and can't find the errors

[haliaetos]

Dear friends,

        My code can't run,the matlab is always busy.If I change the M to smaller ,the problem will get the solution ,I can't find the errors,whether my code has too many circulations,how to improve it?the following is my code.

tic

load data

[m,n]=size(d);

M=n;

y0=binvar(m,n,'full');

y1=binvar(m,n,'full');

mult=sdpvar(m,n,'full');

Obj1=0;Obj2=0;

for i=1:M

    for j=1:M

        Obj1 = Obj1 + (h(i)*d(i,j)*(1-q(j))+mult(i,j))*y0(i,j);%objective=h'*d*(1-q);

        for r=1:M

            if r~=j

                qq(r)=q(r)*y0(i,r);

            end

        end

        Obj2 = Obj2 + (h(i)*d(i,j)*sum(qq)+mult(i,j))*y1(i,j);

    end

end

ObjValue=Obj1+Obj2;

%%modeling the constraints

Constraints=[];

for i=1:M

    Constraints=[Constraints,sum(y0(i,:))==1,sum(y1(i,:))==1];

end

Constraints = [Constraints,mult>=0];

%set parameters

options = sdpsettings('verbose',0,'solver','cplex','cplex.mip',1);

sol = solvesdp(Constraints,ObjValue,options);

if sol.problem == 0

    y0= double(y0);

    y1= double(y1);

    mult= double(mult);

else

    display('Hmm, something went wrong!');

    sol.info

    yalmiperror(sol.problem)

end

toc

[haliaetos]

First,I can't get the result, so I think the problem  has too many  circulations,but the matlab  run the   "ObjValue=Obj1+Obj2;" Elapsed time is 204.072802 seconds., however if I put the d instead of the y0,y1,mult,the matlab will run fast.Obviously,the  sdpvar in yalmip affect the effcet, when I wait 4 minutes ,it run the solversdp,but it says  Error using eig ，Out of memory.  Error in lmi/categorizeproblem 

[Johan Löfberg]

To begin with,why do you expect cplex to solve this? You are setting up a bilinear program

>> sol = solvesdp(Constraints,ObjValue,options)

sol =

    solvertime: 0
          info: 'Solver not applicable (cplex)'
       problem: -4
    yalmiptime: 0.0960

[haliaetos]

A bilinear program？It has only one objective,but I rewrite it two , I am sure it can be solved by cplex,because the modeling is built in the paper "Reliable facility location design under disruption"(you can refer to the 3.2.1 in this paper).He or she  had solve it by cplex.My modeling is similar with it.SO,whether my code has some problems or not.Thank you very much.

And why the elapsed time is long  when the matlab run to the "ObjValue"

[Johan Löfberg]

It is bilinear since you multiply decision variables, e.g., mult with y0. However, it is easily transformed to a MILP by using the fact that continuous*binary can be modelled using logic constraints

a*b where b is binary is replaced with continuous variable c, and the constraints

[-M1*b <= c <= M1*b, -M2*(1-b) <= c-a <= M2*(1-b)]

where M1 and M2 are suitable selected small-but-sufficiently-large constants (called big-M). YALMIP has automatic support for such reformulations through the command binmodel.

However, there seem to be some performance issue in setting up the model at the moment. The best thing would be if you vectorized the definition of the objective. The for-loops makes it extremely slow.

[Johan Löfberg]

And there are logical bugs in your code

   for r=1:M
            if r~=j
                qq(r)=q(r)*y0(i,r);
            end
        end

When j = 1, you will get something like [0 value value value . . ]. Then, when j = 2, you want, I presume [value 0 value value . . ]. However, you will get [value value value ...] since the data in the second position will be kept from the last iteration

[haliaetos]

      On the basis of your analysis,if my problem is really the bilinear problem,it will not be solved by cplex,May I have to solve it through the command binmodel.,can I refer to some examples?Mealwhile I know the for-loops will make the my problem slow,but I  can"t code the part of r~=j.and the 4minutes is too long.

      And you called the logical bugs in my code,I entirely follow the modeling to code.if r=j,the y1(i,j)=1;the y0(i,j)=0,so the qq(r)=0;

[Johan Löfberg]

binmodel is explained on the Wiki. However, you should not use binmodel, since it requires the bilinear model as a start, and the bilinear model is absolutely huge. I think you have failed to appreciate how many terms there is in your model

% The following (almost vectorized) version still takes a long time to run, and in the end, Obj2 contains over 115000 bilinear terms. Absolutely ginormous in the context of mixed-integer programming

QQ = sdpvar(length(q));

for i=1:M    
    for j=1:

M  
        qq = q'.*y0(i,:);qq(j) = 0;
        QQ(i,j) = sum(qq);  
    end    
end

H = repmat(h,1,length(h));
Q = repmat(q',length(q),1);
Obj1fast = sum(sum((H.*d.*(1-Q) + mult).*y0));
Obj2fast = sum(sum((H.*d.*QQ + mult).*y1));

You have to do the linearization manually directly on the expressions (H.*d.*QQ).* y1, mult.*y1 etc. Very straightforward to do, use precisely the linearization strategy I explained above, but in a matrix format.


In your code qq(r) will be q(j-1)*y0(i,j-1) when r = j and j>1 (since you never clear qq, so qq is occupied with data from last loop)

[haliaetos]

  First, I run your code , my matlab can't run the results,I always doubt that using the sdpvar or binvar to define the variable will make the program run slower than number value ,so I replace the y0,y1,mult with  d in my data  in same for-loops,the matlab can run the results fast.

clc;clear all;

tic

load data

[m,n]=size(d);

M=n;

y0=binvar(m,n,'full');

y1=binvar(m,n,'full');

mult=sdpvar(m,n,'full');

QQ = sdpvar(length(q));

for i=1:M     

    for j=1:M   

        qq = q'.*y0(i,:);

        qq(j) = 0;

        QQ(i,j) = sum(qq);  

    end    

end
toc

    Second, since the  mult and y1 are variables so mult.*y1 is the bilinear program that you called.I need to build a  continuous variable c to replace the y1 like  a*b where b is binary is replaced with continuous variable c, and the constraints

[-M1*b <= c <= M1*b, -M2*(1-b) <= c-a <= M2*(1-b)]

.

[Johan Löfberg]

What do you mean you can't run the results? The code runs. It is horribly slow though, so I left it to you to vectorize the definition of QQ

I have no idea what you are saying in the next sentences.

You have to create logical models for

H.*d.*(1-Q) TIMES y0
mult TIMES y0
H.*d.*QQ TIMES y1
mult TIMES y1

I.e, 4 * (M^2) logical models


[haliaetos]

I  know your purpose,your code should have helped me to simple the code,so  first to calculate the QQ,but only calculating the QQ will last  over 30 minutes ,let alone to create logical models,I have not start to create the logical models.

Now I want to know how to get the QQ faster,why is it slow?Then I decide to start to create the logical models.

[haliaetos]

And if my mult is constant not the sdpvar,whether my problem will be easy or not?

[Johan Löfberg]

Yes, it was curiously slow. I'll try to understand why. Here are some fast alternatives

tic
QQ = [];
for i=1:M
    Qrow = [];

    for j=1:M  
        qq = q

'.*y0(i,:);
        qq(j) = 0;      
        Qrow = [Qrow sum(qq)];
    end    
    QQ = [QQ;Qrow];
end
toc

tic
QQ2 = [];
for i=1:M    
    Qrow = sum(q'.*y0(i,:))-q'.*y0(i,:);        
    QQ2 = [QQ2;Qrow];
end
toc

tic
QQ3 = repmat(sum(repmat(q(:)',M,1).*y0,2),1,M)-repmat(q(:)',M,1).*y0
toc

[Johan Löfberg]

You still have bilinear terms since you multiply QQ, which containts y0, with y1.

[haliaetos]

I approximately know some thing,only if I deal with some about over one binvar,eg, QQ.*y1 ,or QQ(i,j) = sum(qq);the matlab  will obviously slow.

So,if I want to work out my problem,I must overcome question about  using the binvar,If your yalmip had solve the problem, my question will turn to creat the logic models.

[Johan Löfberg]

No idea what you are trying to say

[haliaetos]

My idea is intend to transform my problem to a MILP according  to your help. I am trying to create logical models for "H.*d.*QQ TIMES y1",but it still can'trun.I don't know whether it is right or not.

clc;clear all;

load data

[m,n]=size(d);

M=n;

d_ave = sum(sum(d))/(M*M);

mult = repmat(h*d_ave/10,1,M)

y0=binvar(m,n,'full');

y1=binvar(m,n,'full');

QQ = [];

for i=1:M    

    Qrow = sum(q'.*y0(i,:))-q'.*y0(i,:);         

    QQ = [QQ;Qrow];

end

H = repmat(h,1,length(h));

Q = repmat(q',length(q),1);

Obj1fast = sum(sum((H.*d.*(1-Q) + mult).*y0));

% a*b where b is binary is replaced with continuous variable c, and the constraints

%  [-M1*b <= c <= M1*b, -M2*(1-b) <= c-a <= M2*(1-b)]

M1=1;M2=1;c=semivar(m,n);%c=sdpvar(m,n);

Obj2fast = sum(sum((H.*d.*QQ + mult).*c));

ObjValue = obj1fast + obj2fast ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Constraints = [];

for i = 1:M

    Constraints = [Constraints, sum(y0(i,:)) == 1,sum(y1(i,:)) == 1];

end

%　Constraints = [Constraints,mult >= 0];

Constraints = [Constraints, -M1*y1 <= c <= M1*y1, -M2*(1-y1) <= c-H.*d.*QQ <= M2*(1-y1)];

[Johan Löfberg]

I don't understand what you are trying to do, as you still multiply decision variables (QQ and c)

QQ*y1 should be modelled by a new variable so your new objective should be

sum(sum((H.*d.*QQ_TIMES_Y1 + mult.*y1));

[haliaetos]

I am sorry to trouble you  all the time ,but do you remember  your answer:

"It is bilinear since you multiply decision variables, e.g., mult with y0. However, it is easily transformed to a MILP by using the fact that continuous*binary can be modelled using logic constraints  a*b where b is binary is replaced with continuous variable c, and the constraints

[-M1*b <= c <= M1*b, -M2*(1-b) <= c-a <= M2*(1-b)]


where M1 and M2 are suitable selected small-but-sufficiently-large constants (called big-M)." I had not deal with this problem,so I imitate your methods.

Now  I follow you to deal with "H.*d.*QQ TIMES y1",I replace y1 to b like you ,but the c is the continuous variable, so I redefine the c ,the result is that the QQ and c is still the decision variables,I need to learn to how to deal with it,eg  the examples about your big -M.

[Johan Löfberg]

You never replaced the product a*b with the new variable c, you simply replaced b with c.

something+a*b <= 0

is written as

[something+c <= 0, -M1*b <= c <= M1*b, -M2*(1-b) <= c-a <= M2*(1-b)]


Think through the combinatorial cases in b and realize that it always holds that c=a*b if b is binary

[haliaetos]

With the your help,I have a little understand the theory about the big- M ,so I set the M1 and M2 according to the max（max(H.*d）），I get the the value of the ObjValue, but the results was different from the results with the genetic algorithm ,the former is 4573.2 ,the latter is 4581.6.So I had changed the M1,M2,the matlab can still run the  different resluts.Can I believe the results?what is the true value?

[Johan Löfberg]

Is the other solution feasible?

Otherwise, the gap between these two solutions are 0.2%, so you might have to decrease the solver tolerance to force the solver to continue its search for longer.

[Johan Löfberg]

And I presume you know what you do when you derive M1 and M2 (a bit odd the M2 doesn't depend on any known bound on mult...)