network problem

[Tung]

Dear johan,

   I have a question about linear programming in terms of a network. I hope you could spend a little time to help me transforming the problem to yalmip model. 

Suppose that I have 4 nodes, and each node will connect to others through the links. Let see in the scenario below,

A ---- C

 |   X  | 

B ---- D

- Let say A transmits data packets (traffic demand) to D, where A is a source, and D is a destination. We have several ways for transmissions, including AD, ACD, ABD, ACBD, and ABCD. Note that packets transmission can be done for traffic allocation on one link or summation of traffic allocation in several links in order to reach destination D. In this case, traffic allocation (TA) is the output, and the inputs are traffic demand of source A and the number of links (L). 

- My goal is how to minimize the number of links when we have a given traffic demand from source A to transmit to destination D.

- Let x is a variable which is corresponding to each link. (x =1 means there are a link established btw two nodes, otherwise =0).

- Let a is traffic allocation variable on each link.

Now, my model is for a traffic demand from A-D:

-- model

Minimize c'*x   

   s.t. {Traffic demand (A-D) {sum of paths from (A-D) { a_{links on path}*{production of x}}}} <= Traffic capacity on a link;

   s.t. {sum of paths from (A-D) { a_{links on path}*{production of x}}} >= Traffic_demand;   % this constraint means that the left side has to meet the traffic demand from source A on the right side.

-- end model

-- Note:

The production of x means each link is corresponding to a x, when a path has several links like path ACD has two links AC and CD. x on each link, which belongs to the path, will be times together to see if this production is 1 which means the path is valid, otherwise invalid. The production of x is binary variable.

-- 

In matlab, I suppose that the inputs are the number of links and the number of paths from single source-destination A-D. The outputs are x and a. 

How do we form the model in yalmip and after running it, how can I get the specific value from each x, and traffic allocation value on each link ? 

One more thing is because the constraints are added through the return values from functions, so I don't know how to get the exactly values of x and a from functions. 

Please help me to figure out how to solve it. Thank you.

Best regards,

Tung 

 

[Johan Löfberg]

This sounds like a very standard network problem (which I know close to nothing about). Start by formulating it as a problem involving matrices and linear algebra (your network will define some sort of an incidence matrix). Turning on/off links will turn the problem into a binary problem as you indicate. At that point, we can start discussing suitable YALMIP models.

Just google network optimization and you get loads of hits. First one http://www.ee.ucla.edu/ee236a/lectures/networks.pdf

[Tung]

Dear Johan,

   Thanks for your reply. I got the result for my problem formula. But I need your help to interpret my results. Let say x is a matrix vector (e.g, x = binvar(8,1); ) which are decision variables. Note that my problem needs to be minimized. (Min cost'*x; ) 

When I put the input value with small amount, the output x=[1 0 0 0 0 0 0 0]; and cost =1; 

When increasing the value with larger amount, x = [1 0 0 1 1 0 0 0]; and cost =3;

When increasing the value with more larger amount, cost =Inf; and objective = 4; 

Note that the maximum cost=8; (define cost=ones(8,1); )

As I guess, cost = inf and objective = only 4 but no more higher is because the problem is minimizing. Is it correct?

Thanks.

Best regards,

Tung  

[Johan Löfberg]

I don't understand what you are asking. The second argument to solvesdp is minimized if that is what you wonder. Using inf in the coefficients of the objective is a very bad idea (since all x then yield the same objective, inf, -inf, or undefined)

[Tung]

Dear Johan,

   Thanks for reply. My question is that why the cost (output) can not be higher than 4, when I do increasing the input value? In my guess, the cost could be 5 or 6 when I increase the input value. Why x variables can not be more than 4 decision values? Why the objective stops at 4, and the cost is infinite? I don't understand this point. Please explain to me about that.

Thanks. 

Best regards,

Tung

[Johan Löfberg]

I think you would have to post some code to show what you mean

[Tung]

Dear Johan,

   Here is my code. Note that the traffic as below is from T(1,5) and T(2,5);. Suppose that you can add the traffic T(1,6) to 100, and T(3,5) to 50. The cost is increased higher. But when I try to increase higher, the cost goes to the infinite, and the objective obj stops at 4. My intention is that the cost should be increased over 4, to 5, or 6 if I increase the traffic correspondingly. Is my intention correct? or it is alright?

W=100; % capacity of type A

V=200; % capacity of type B

% % traffic

 T = [ 0 0 0 0 100 0;

       0 0 0 0 0 100;

       0 0 0 0 0 0;

       0 0 0 0 0 0;

       0 0 0 0 0 0;

       0 0 0 0 0 0;];

% % matrix links

N = [0 1 1 0 1 1;

     1 0 0 1 1 1;

     1 0 0 1 1 1;

     0 1 1 0 1 1;

     1 1 1 1 0 1;

     1 1 1 1 1 0;]; 

cost = ones(8,1);

% def variables

x = binvar(8,1); % x = 0 | 1;

a = sdpvar(8,1); % 0 <= a <= 1;

% type A constraints

l1 = l1 = T(3,5)*a(1)*x(1) + T(3,6)*a(1)*N(5,6)*x(1) + T(1,5)*a(1)*N(1,3)*x(1) + .. + T(4,6)*a(1)*N(4,2)*N(2,1)*N(1,3)*N(5,6)*x(1);  

l2 = T(3,6)*a(2)*x(2) + T(3,5)*a(2)*N(6,5)*x(2) + T(1,6)*a(2)*N(1,3)*x(2) + .. + T(4,5)*a(2)*N(4,2)*N(2,1)*N(1,3)*N(6,5)*x(2);  

...

l8 = ..;

% type B constraints

s1 = T(1,5)*a(1)*N(1,2)*N(2,4)*N(4,3)*x(1) + T(1,6)*a(1)*N(1,2)*N(2,4)*N(4,3)*N(5,6)*x(1) +..+ T(4,5)*a(6)*N(4,2)*N(2,1)*N(6,5)*x(6);

...

s5 = ...;

% type C constraints

tr1 = T(3,5)*a(1)*x(1) + T(3,5)*a(2)*N(6,5)*x(2) + .. + T(3,5)*a(8)*N(3,4)*N(4,2)*N(6,5)*x(8);

...

tr8 = ...;

% def constraints

F = [0<=x<=1, 0<=a<=1];

F = [F, l1<=W, l2<=W, l3<=W, l4<=W, l5<=W, l6<=W, l7<=W, l8<=W];

F = [F, s1<=V, s2<=V, s3<=V, s4<=V, s5<=V];

F = [F, tr1 >= T(1,5), tr2 >= T(1,6), tr3 >= T(2,5), tr4 >= T(2,6), tr5 >= T(3,5), tr6 >= T(3,6), tr7 >= T(4,5), tr8 >= T(4,6)];

% objective function

obj = cost'*x;

% solver

options = sdpsettings('solver','bnb');

% solve

solution = solvesdp(F,obj,options);

% get output

obj = double(obj);

x = double(x);

a = double(a);

[Tung]

I would like to show 3 cases that I am talking about the problem.

Case 1: This is the output with a small amount of the input T(1,5)

-----

Starting YALMIP integer branch & bound.

* Lower solver   : fmincon-standard

* Upper solver   : rounder

* Max iterations : 300

 

Warning : The continuous relaxation may be nonconvex. This means 

that the branching process is not guaranteed to find a

globally optimal solution, since the lower bound can be

invalid. Hence, do not trust the bound or the gap...

 Node       Upper       Gap(%)      Lower    Open

    1 :          Inf      Inf      1.000E+00   2  Maximum iterations or time limit exceeded 

    2 :          Inf      Inf      1.000E+00   3  Successfully solved 

    3 :          Inf      Inf      1.000E+00   4  Successfully solved 

    4 :          Inf      Inf      1.000E+00   5  Successfully solved 

    5 :          Inf      Inf      1.000E+00   6  Maximum iterations or time limit exceeded 

    6 :    1.000E+00     0.00      1.000E+00   5  Successfully solved 

+   6 Finishing.  Cost: 1

x =

     0

     1

     0

     0

     0

     0

     0

     0

------

Case 2: When I increase the input with T(1,5); T(1,6); T(2,5); T(3,5); T(3,6); T(4,5); the output is following as below

------

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.

RCOND =  7.955966e-17. 

> In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\backsolveSys.p>backsolveSys at 18

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\solveAugSystem.p>solveAugSystem at 23

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\projConjGrad.p>computeProjResidual at 159

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\projConjGrad.p>projConjGrad at 95

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\tangentialStep.p>tangentialStep at 36

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\computeTrialStep.p>computeTrialStep at 104

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\barrier.p>barrier at 348

  In fmincon at 841

  In callfmincon at 64

  In bnb_solvelower at 136

  In bnb>branch_and_bound at 603

  In bnb at 279

  In solvesdp at 358

  In test8 at 127 

    8 :          Inf      Inf      4.000E+00   1  Infeasible problem 

    9 :    4.000E+00     0.00      4.000E+00   0  Successfully solved 

+   9 Finishing.  Cost: 4

>> x

x =

     0

     1

     1

     0

     1

     0

     1

     0

------

Case 3: Increasing with more traffic, it turns out as follows.

-----

Warning: Matrix is singular to working precision. 

> In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\backsolveSys.p>backsolveSys at 18

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\solveAugSystem.p>solveAugSystem at 23

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\projConjGrad.p>computeProjResidual at 159

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\projConjGrad.p>projConjGrad at 95

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\tangentialStep.p>tangentialStep at 36

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\private\computeTrialStep.p>computeTrialStep at 104

  In C:\Program Files\MATLAB\R2012a\toolbox\optim\optim\barrier.p>barrier at 348

  In fmincon at 841

  In callfmincon at 64

  In bnb_solvelower at 136

  In bnb>branch_and_bound at 603

  In bnb at 279

  In solvesdp at 358

  In test8 at 127 

    1 :          Inf      Inf            NaN   0  Infeasible problem 

+   1 Finishing.  Cost: Inf

>> x

x =

    0.5000

    0.5000

    0.5000

    0.5000

    0.5000

    0.5000

    0.5000

    0.5000

-----

My question is that could the cost in case 3 be over 4? Why it cannot be over 4 and turns out to the infinite?

Thanks.

Best regards,

Tung

[Johan Löfberg]

For future reference, please supply short succint code which actually runs.

The problem here is that you setup a very hard mixed-integer nonconvex program since you multiply the decision variables a and x. Hence, the integer relaxations are not simple LPs or QPs, but nonlinear programs. This is a very bad way to model the problem.

Instead, you should exploit the fact that a*x where 0<=a<=1 and x binary can be linearized by introducing a new variable z (replacing a*x) and

[0 <= z <= x,-(1-x) <= z-a <= 1-x]

This you have to do for every bilinear product. YALMIP can help you through the command binmodel.

[Tung]

Dear Johan,

   Thanks for reply. I tried convert non-LP to LP but it seems there are somethings more to be concerned. My code is attached in this post, please take a look at it and please explain to me if it does a wrong way. In addition, can you explain to me how we have this constraint [-(1-x) <= z-a <= 1-x]. Thanks.

Best regards,

Tung  

[Johan Löfberg]

You seem to have missed the whole idea. You should not add the constraint z==a*x. This bilinear nasty constraint is guaranteed to hold if you add the linear constraints [0<= z <= x, -(1-x) <= z-a <= 1-x] (right?, check!). That is the whole idea in introducing the new variable z.

The problem is infeasible though, already this is infeasible

F = [0<=x<=1, 0<=a<=1];

F = [F, 0 <= z(1) <= x(1) ,-(1-x(1)) <= z(1)-a(1) <= 1-x(1)];
F = [F, 0 <= z(2) <= x(2) ,-(1-x(2)) <= z(2)-a(2) <= 1-x(2)];
F = [F, 0 <= z(3) <= x(3) ,-(1-x(3)) <= z(3)-a(3) <= 1-x(3)];
F = [F, 0 <= z(4) <= x(4) ,-(1-x(4)) <= z(4)-a(4) <= 1-x(4)];
F = [F, 0 <= z(5) <= x(5) ,-(1-x(5)) <= z(5)-a(5) <= 1-x(5)];
F = [F, 0 <= z(6) <= x(6) ,-(1-x(6)) <= z(6)-a(6) <= 1-x(6)];
F = [F, 0 <= z(7) <= x(7) ,-(1-x(7)) <= z(7)-a(7) <= 1-x(7)];
F = [F, 0 <= z(8) <= x(8) ,-(1-x(8)) <= z(8)-a(8) <= 1-x(8)];

F = [F, l1<=Wn, l2<=Wn, l3<=Wn, l4<=Wn, l5<=Wn, l6<=Wn, l7<=Wn, l8<=Wn];
F = [F, s1<=We, s2<=We, s3<=We, s4<=We, s5<=We];
F = [F, tr1 >= T(1,5)]

and don't use bnb. You must install a real MILP solver, such as cplex, gurobi or mosek (all with academic licenses), or a non-commercial alternative, of which I prefer scip.

[Tung]

Dear Johan,

   Thanks for reply. I got your point. In addition, in order to satisfy this constraint [0<= z <= x, -(1-x) <= z-a <= 1-x], a must be -1<=a<=1, not 0<=a<=1 as in the code. I used CPLEX (academic version) to solve the problem and it works perfectly. Again, thanks very much for your help.

Have a nice day,

Best regards,

Tung  

[Johan Löfberg]

a is already constrained in your other code segments to 0<=a<=1, but that is not interesting anyway. The model I give is not supposed to constrain a (on the contrary, it should not constrain a at all, constraints on a should be declared outside and are only used to derive the big-M model)

The constraint

[0<= z <= x, -(1-x) <= z-a <= 1-x]

Says that either
x = 0 implying z = 0 (thus correctly multiplying a*x when x is 0)
x = 1 implying z=a   (thus correctly multiplying a*x when x is 1)

a can never become negative in the model

[Tung]

Dear Johan,

   Thanks for your reply. You are right on statement of a. 

   In addition, I have a question here is that if I want to insert an algorithm to solve the problem instead of a solver (in the case I have a better algorithm to solve the problem), is it possible to do in yalmip? Then, how/where to implement my algorithm in order to solve the abovementioned model ? Please advice to me this point. I really appreciate that. Thanks.

Best regards,

Tung

  

[Johan Löfberg]

Your only possibility would be to write a function which takes the same inputs as a solver in YALMIP, and then export the model to that specific format and send the data to your algorithm

[Tung]

Dear Johan,

   Thanks for your reply. I mean that in order to solve the problem through the algorithm, I have to relax the MILP to LP, but how can I get the LP relaxation in YALMIP? thereby, I can add these relaxed variables into my algorithm to find an optimal solution. In addition, how can I get the upper bound and lower bound values of the solver used. Please explain to me. Thanks.

[Johan Löfberg]

What do you mean with "get the LP relaxation in YALMIP". Do you mean how to solve the relaxation in YALMIP (define them as sdpvars instead of binaries, or use the option 'relax' in sdpsettings) or do you mean how to get the numerical data which defines the LP relaxation (use export, and export the model to, e.g., cplex numerical format, and simply skip the binary constraints)

[Tung]

Dear Johan,

   Thanks for your reply. My intention is to get the data after LP relaxation form the original problem, so that I can use those relaxed variables to run the algorithm (to check each relaxed variable is whether feasible or not, based on given constraints). It seems likely that the second method (export) is close to my question, but when I run the export, I am confused about which elements are relaxed variables that I need. For instance, I ran the opcode-2.m, which is attached here, and I got the model information as follows.(I used cplex solver here)

---------------------------------------------------------

model =         H: []

                     f: [24x1 double]

                    Aineq: [83x24 double]

                    bineq: [83x1 double]

                    Aeq: []

                    beq: []

                    lb: [24x1 double]

                    ub: [24x1 double]

                    x0: []

                    options: [1x1 struct]

                    verbose: 1

                    integer_variables: [1x0 double]

                    binary_variables: [1 2 3 4 5 6 7 8]

                    semicont_variables: [1x0 double]

                    K: [1x1 struct]

                    ctype: [24x1 char]

                    Qi: []

                    Li: []

                    ri: []

                    NegativeSemiVar: []   

--------------------------------------------------------------------------------

But when I use the options relax=1 instead of the cplex solver (without the solver), the result is following:

--

model =     A: [83x24 double]

                 obj: [24x1 double]

                 sense: [83x1 char]

                 rhs: [83x1 double]

                 lb: [24x1 double]

                 ub: [24x1 double]

                 objcon: 0

                 vtype: [24x1 char]

                 Q: [24x24 double]

                 params: [1x1 struct]

                 NegativeSemiVar: []

--   

So, in both cases, which element I can get the relaxed variables (after LP relaxation) from the original problem before the solver intervenes those variables? (of course, solvesdp is disable to check the relaxed variables in opcode-2.m). Please explain to me. Thanks.

Best regards,

Tung

[Johan Löfberg]

You are not clear. Do you want the model data, or the variables obtained when solving the relaxation?

In the first model above, you have clearly obtained data for

min obj'*x
s.t Aeq*x == beq, Aineq*x <= Bineq, lb <= x <= lb, x(1:8) binary.

Trivially, the LP relaxation is the model where you simply skip the binary condition on x(1:8)

The second model you have exported has been exported for gurobi instead it looks like. The model looks exactly the same apart from a different format, and, of course, that there is no binary condition (vtype should thus be a vector of 'C' which is Gurobis way of saying the variable should be binary)

[Tung]

Dear Dr. Johan,

   Thanks for your reply. I actually want the variable obtained when solving the relaxation. Based on those relaxed variables, I can apply them to my algorithm to check whether feasible or not. You make me confused about how to get variables obtained after relaxation. Is it just only set relax = 1 in sdpsettings and put this option in solvesdp? I think it would be better to get variables obtained after relaxation and before those variables going to solvesdp. Because after solvesdp, variables have already been solved by the solver, therefore, after this process, there is no job for my algorithm. Is it correct?  Please let me know how to get variables obtained after relaxation and before they go to solvesdp. Thanks.

[Tung]

Dear Dr. Johan,

   I came up with a simple algorithm as follows.

--

Relax MILP to LP => (*) solve LP (by a solver?!): if infeasible variable => stop

                                                                    if feasible variable => Store this variable   

                                                                    Check until all variables could be feasible variable? => if Yes:   optimal solution

                                                                                                                                                   if No: check other if possible

                                                                                                                                                   Establish new LP with stored variables (constraints) to solve ( go to (*) )

--                                                     

Please advise to me that this algorithm is possible to solve the LP or not. Thanks.

Best regards,

Tung 

[Johan Löfberg]

First, I don't see what in your algorithm would stop the solver from returning exactly the same solution as you just found and stored, when you sovle the problem the second time. You are not cutting or branching anything

Second: Why on earth would you want to try to develop this new (not functioning) version of a MILP solver, when you have loads of solvers available to solve this problem? You are reinventing a square wheel.

[Tung]

Dear Dr. Johan,

   Thanks for your reply. I think that there are some specific problems that could not be solved by such solvers, for example, a large-scale network. It could take a very long time to solve such problems because the solvers cannot know exactly specific purposes of the problems. Instead, designing an algorithm in polynomial-time complexity that could probably solve the problem efficiently for the large scale network, in terms of reducing in time computation. Anyway, thank you for your help.

[Johan Löfberg]

Yes, but when you have to develop your own algorithm, you will have to take care of the modelling manually, i.e., no support from YALMIP

[Tung]

Dear Dr. Johan,

    First of all, thank you for your support on this problem. I got stuck at this error as below and don't know how to solve it yet. Could you help me out? I also attached the running file into this post.  I appreciate it. Thank you 

-- error message--

Error using lmi/horzcat (line 11)

One of the constraints evaluates to a FALSE LOGICAL variable

Error in subs2cores16 (line 554)

F = [F, tr(1) >= (T(1,9) + T(9,1)), tr(2) >= (T(2,9) + T(9,2)), tr(3) >= (T(3,9) + T(9,3)), tr(4) >= (T(4,9) + T(9,4)), tr(5) >= (T(5,9) + T(9,5)), tr(6) >= (T(6,9) +

T(9,6)), tr(7) >= (T(7,9) + T(9, 

--

[Johan Löfberg]

You would have to supply a smaller example. I cannot run it on my machine (I cannot even open the file, it is too large for matlab)

You are effectively doing 0>=1 or similiar somewhere in that line, similiar to

T(1,9)=1;
T(9,1)=tr(1);

tr(1) >= (T(1,9) + T(9,1))...

or something trivial like that

[Tung]

Here is a simple example between two subnets (6 nodes). When you change the value in matrix T, binvar x will change (decision variables). However, when I run with a larger network (16 nodes) in two subnets, it works well with small workloads, but with the high workload, it is crashed at the error as above. I still not know why, whether my code is wrong or the solver's constraints limitation. I am thinking that when the number of connections increases, the constraints apparently are increased, so that it is bounded, otherwise, my constraints configuration has some problems. If so, why it can run in the small workload. 

[Johan Löfberg]

Once again: You are adding a garbage constraint

This line

One of the constraints evaluates to a FALSE LOGICAL variable

Error in subs2cores16 (line 554)

F = [F, tr(1) >= (T(1,9) + T(9,1)), tr(2) >= (T(2,9) + T(9,2)), tr(3) >= (T(3,9) + T(9,3)),...

tells you that one of those constraints is not a symbolic expression, but false logical (0). Hence, some of these expressions

tr(1)-(T(1,9) + T(9,1))
tr(2)-(T(2,9) + T(9,2))


evaluates to a negative number, and can thus never be non-negative (as it doesn't depend on any decision variable). You have to put a break on that line and investigate those expressions, and then understand why you are defining a trivially infeasible model

[Tung]

Dear Dr Johan,

    I am sorry for late replying. Thanks for your hints. I've already fixed that. But now I got another problem. When I increase the traffic demands in the traffic matrix, the number of decision variables x should be increased. For example, in the attached file here I have a small network with 2 subnets as follows.(subnet1 contains 1,2,3,4, and subnet2 contains 5,6)

1---3    5

|      |     |

2---4    6 

we have decision variables x (4x2=8) corresponding to the link established from one subnet to another subnet. When the traffic demands increases, it is apparent that decision variables x should be increased. With a small traffic, it should be x[1,0,0,0,0,0,0,0] since it only uses a link to connect communication between two subnets and meet all given constraints. When a larger traffic is injected (through the traffic matrix in the file), x should be increased correspondingly like x[1,0,1,0,1,1,1,0] and meet all constraints. 

However, when I try to increase the traffic demands, x does not increase appropriately. It always keeps only one "1" or nothing (means all zeros). I still not figure out how to fix it. Thanks for supporting me on this problem, hope you have a merry christmas and happy new year.

Best wishes for your holidays.

Best regards,

Tung    

[Tung]

I mean that could you advice me how to figure out this problem above? Many thanks.

Best regards,

Tung

[Johan Löfberg]

Impossible to answer the question with such a complex and verbose model.

[Tung]

Thanks Dr. Johan,

   I wrote a concise program that describes what I want to do as in the file attached. Now, the problem is that the result of decision variables x always keeps TWO "1" variables whatever I change the traffic demands. This means that when the traffic demands change, the decision variables x might change the position, like x[0 1 1 0 0 0 0 0] to x[1 0 0 0 1 0 0 0], but it does not change the number of variables. What I would like to expect here is that when I change the traffic demands T, for instance, T is increased, x should be more "1", and when T is decreased to a small traffic amount, x should be ONE "1". Can you figure out any bugs here?  Many thanks. 

[Johan Löfberg]

Absolutely impossible for us to see in such a elementwise unrolled model  Could be a typo in one single digit or anything.

The minimal number of non-zeros in the current model is 2 (seen by using obj = sum(x))

Perhaps you should start by creating a high-level model and let YALMIP take care of the linearization of the model (binmodel), to see if your basic model is wrong, or if your linearization is wrong.

[Tung]

Do you mean that instead of using z(i) in constraints, I should use x(i)*x(j)*x(k), and using a high-level model for inputs.

As I understand, it could be something like this:

l(1)=T(1,2)*z(15)+N(3,4)*N(2,4)*T(1,2)*z(16)+N(3,4)*T(1,2)*z(17)+...  will be => l(1)=T(1,2)*x(1)*x(3)+N(3,4)*N(2,4)*T(1,2)*x(1)*x(5)+N(3,4)*T(1,2)*x(1)*x(5)*x(8)*x(4)+...

so on so forth with s(i) and tr(i). And then yalmip can linearize the x variables' multiplication and extract to the constraints as I did. Note that all x variables are binary variables.

Then I create a high-level model like:

p = cost'*x;

[plin, F]=binmodel(p, [0<=x<=1, 0<=z<=1]);

F = [F, l(1)<=Wn, l(2)<=Wn, l(3)<=Wn, l(4)<=Wn, l(5)<=Wn, l(6)<=Wn, l(7)<=Wn, l(8)<=Wn];

F = [F, s(1)<=We, s(2)<=We, s(3)<=We, s(4)<=We, s(5)<=We];

F = [F, tr(1)>=T(1,5), tr(2)>=T(1,6), tr(3)>=T(2,5), tr(4)>=T(2,6), tr(5)>=T(3,5), tr(6)>=T(3,6), tr(7)>=T(4,5), tr(8)>=T(4,6), tr(9)>=T(5,1), tr(10)>=T(5,2), tr(11)>=T(5,3), tr(12)>=T(5,4), tr(13)>=T(6,1), tr(14)>=T(6,2), tr(15)>=T(6,3), tr(16)>=T(6,4)];

options = sdpsettings('solver','cplex');

solution = solvesdp([F,plin],obj,options);

Am I right? can you explain a bit more? Many thanks.

[Tung]

Or do you mean that I should extract each part in l(i) s(i) tr(i) like this

From this: l(1)=T(1,2)*x(1)*x(3)+N(3,4)*N(2,4)*T(1,2)*x(1)*x(5)+N(3,4)*T(1,2)*x(1)*x(5)*x(8)*x(4)+...

p1= T(1,2)*x(1)*x(3);

p2= N(3,4)*N(2,4)*T(1,2)*x(1)*x(5);

p3= N(3,4)*T(1,2)*x(1)*x(5)*x(8)*x(4);

..

pn=...

[plinear1, plinear2, plinear3,...,plinearn, F] = binmodel(p1,p2,p3,...,pn, [0<=x<=1]);

options = sdpsettings('solver','cplex'); 

obj = cost'*x;

solution=solvesdp([plinear1, plinear2, plinear3,...,plinearn, F], obj, options);

Many thanks.

[Johan Löfberg]

Yes, as a test that is what I would do

Note, the bounds on x and z are redundant as they are defined as binary

[Tung]

Thanks Dr. Johan,

   I wrote again as in the file, and have an error as follows.

--

reached traffic demands

...

reached pN constraints

reached a binmodel

Error using binmodel (line 57)

Arguments should be SDPVAR or SET objects

Error in test151 (line 112)

[plinear1, plinear2, plinear3, plinear4, plinear5, plinear6, plinear7, plinear8, plinear9,

plinear10, plinear11, plinear12, plinear13, plinear14, plinear15, plinear16, plinear17, plinear18,

plinear19 

--

Now the program only has x variables' multiplications, how can z variables be defined in yalmip?

I think there's something wrong here, right?

I attached the update file in this post. 

Many thanks.

[Tung]

I have checked the previous one with my linearization and it works with one variable instead of TWO as before. But now, when I increase the traffic demands, it shows either only one decision variable x, like x [0,0,0,0,1,0,0,0] or all-zeros like x [0,0,0,0,0,0,0,0] if it doesn't meet the constraints. What I expect that when the traffic demands are increased, one decision variable can't hold all traffics, it could increase decision variables x, with more than "1", e.g. x [0,1,0,1,0,0,0,0] or [0,1,0,1,0,0,1,0] to handle all given traffics. Do you see any bugs I could have ? The file is attached in this post. 

--

Note that the result of this file as follows.

--

z  =   54

Linear scalar (real, binary, 47 variables, current value : 0)

Linear scalar (real, binary, 47 variables, current value : 0)

Linear scalar (real, binary, 46 variables, current value : 0)

Linear scalar (real, binary, 46 variables, current value : 0)

Linear scalar (real, binary, 48 variables, current value : 24) <==== x(5)

Linear scalar (real, binary, 48 variables, current value : 0)

Linear scalar (real, binary, 49 variables, current value : 0)

Linear scalar (real, binary, 49 variables, current value : 0)

--

z(54) is corresponding to x(5)=1;

--

Why when the traffic demands are increased, the decision variables can not increase to TWO or more "1"?

Many thanks.

--  

[Johan Löfberg]

It crashes because the first argument into binmodel is 0. binmodel gets confused (why would you want to linear a constant 0, of course this should e handled, but at the moment the code does not check for trivial case).

This can be coded around though by first defining a vector of all things you want to linearize (why don't you work with vectors and/or vectorize your code? It looks horrible, and I hope you are using some kind of code generator to generate this, otherwise it is is almost sure that you are introducing errors)

p = [p1, p2,...]
[plinear,cuts] = binmodel(p);

Then when you solve the problem, it is not sure what it is you want with the vector p. Now, you simply use the model

[F,p]

This is short for

[F, p>=0]

Clearly that is not what you want. In addition to that, you don't add the linearizing constraints to the model. You should have something like

optimize([F,cuts, p>=(<=,==)?],obj)

 

[Tung]

Thanks Dr. Johan,

   Yes, I wrote another code generator for l(i), s(i), tr(i), and p. Now, I changed as you said, and added the cplex solver (optimize([F, cuts, plinear, p], obj, options); ), but it does not work, even with a small traffic demands. The result is as below.

--

...

reached an objective

reached a solver

Warning: Solver not applicable (cplex)

Polynomial scalar (real, binary, 8 variables)

Polynomial scalar (real, binary, 8 variables)

Polynomial scalar (real, binary, 8 variables)

Polynomial scalar (real, binary, 8 variables)

Polynomial scalar (real, binary, 7 variables)

Polynomial scalar (real, binary, 7 variables)

Polynomial scalar (real, binary, 7 variables)

Polynomial scalar (real, binary, 7 variables)

--

When I run without the solver ( optimize([F, cuts, plinear, p], obj); ), the result is as follows 

--

...

reached an objective

reached a solver

* Starting YALMIP integer branch & bound.

* Lower solver   : fmincon-standard

* Upper solver   : rounder

* Max iterations : 300

 

Warning : The continuous relaxation may be nonconvex. This means 

that the branching process is not guaranteed to find a

globally optimal solution, since the lower bound can be

invalid. Hence, do not trust the bound or the gap...

 Node       Upper       Gap(%)      Lower    Open

    1 :          Inf      Inf      5.000E+00   2  Successfully solved 

    2 :          Inf      Inf      5.000E+00   3  Successfully solved 

    3 :          Inf      Inf      5.000E+00   2  Infeasible problem 

    4 :          Inf      Inf      5.000E+00   1  Infeasible problem 

    5 :          Inf      Inf      1.800E+01   2  Successfully solved 

    6 :    1.000E+01     0.00      1.000E+01   0  Successfully solved 

+   6 Finishing.  Cost: 10

Polynomial scalar (real, binary, 8 variables, current value : 0)

Polynomial scalar (real, binary, 8 variables, current value : 0)

Polynomial scalar (real, binary, 8 variables, current value : 0)

Polynomial scalar (real, binary, 8 variables, current value : 0)

Polynomial scalar (real, binary, 7 variables, current value : 0)

Polynomial scalar (real, binary, 7 variables, current value : 20)

Polynomial scalar (real, binary, 7 variables, current value : 0)

Polynomial scalar (real, binary, 7 variables, current value : 0)

--  

Why I cannot run the program using an installed solver (cplex)? 

When running the program without the solver (without the options call), the result shows that the problem has one feasible solution, but in polynomial scalar, not in linear scalar. I don't understand here because the result should be linear scalars after linearizing process in yalmip, am I right? In addition, if it works, why when I increase the traffic demands a little bit, the result shows an infeasible solution, not increasing the number of feasible solutions, so it gets confusing.

I would expect when I increase the traffic demands, it should have more than one feasible solution, but for now, how to make it work with cplex solver in the program? 

PS: The updated file is attached in this post.

Many thanks.

   

[Johan Löfberg]

Your code with constraints

[F, cuts, plinear, p]

is a bad way of writing

[F, cuts, plinear >= 0, p>= 0]

Is that what you intend?

The fact that cplex isn't applicable is since you still constrain the nonlinear vector p. The whole idea with linearization is to skip that, and replace p>=0 with [cuts, plinear >=0]

[Tung]

Thanks Dr. Johan,

   Because after binmodeling, we have plinear which are the result of linearization of x variables. In addition, I intend to add the results after linearizing, so vector p will not be used, since we already have plinear. so it should be like this 

optimize([F, cuts, plinear>=0], obj);

So, at this point, can I use cplex to solve the problem?

Moreover, I wonder that the binmodel can handle x1*x2*x3*x4*x5*x6  to have z<=x1, ..., z<=x6, z>= ((x1+x2+...+x6) - 5) or something like that. I mean how about the maximum number of decision variables' multiplication in the binmodel?

--

I also run the program with this update, but it does not work.

--

...

reached a solver

* Starting YALMIP integer branch & bound.

* Lower solver   : fmincon-standard

* Upper solver   : rounder

* Max iterations : 300

 

Warning : The continuous relaxation may be nonconvex. This means 

that the branching process is not guaranteed to find a

globally optimal solution, since the lower bound can be

invalid. Hence, do not trust the bound or the gap...

 Node       Upper       Gap(%)      Lower    Open

    1 :          Inf      Inf            NaN   0  Infeasible problem 

+   1 Finishing.  Cost: Inf

Polynomial scalar (real, binary, 8 variables, current value : 59.375)

Polynomial scalar (real, binary, 8 variables, current value : 64.375)

Polynomial scalar (real, binary, 8 variables, current value : 60.9375)

Polynomial scalar (real, binary, 8 variables, current value : 55.9375)

Polynomial scalar (real, binary, 8 variables, current value : 66.5625)

Polynomial scalar (real, binary, 8 variables, current value : 61.5625)

Polynomial scalar (real, binary, 8 variables, current value : 45.625)

Polynomial scalar (real, binary, 7 variables, current value : 50.625)

--

Many thanks. 

[Johan Löfberg]

Yes, binmodel derives a linear model

You can have as many variables you want. The model might be bad though, and a manually derived model will most likely be cleaner (but the reason I am asking you to use binmodel temporarily is to make sure your problems don't come from a bug in your linearization)

[Johan Löfberg]

F contains nonlinearties too which you have to linearize

[Johan Löfberg]

My guess is that your whole introduction of p is useless, i.e, it makes no sense since it isn't used anywhere. However, it contains nonlinear expressions which also are used in F. Hence, my guess is what you really want to do is

Flinearized = binmodel(F);
solution = optimize(Flinearized , obj, options);

[Tung]

Thanks Dr. Johan,

   Yes, you are right. p will not be considered because l(i) s(i) tr(i) already contain all x variables multiplication for linearization. I have updated as you mentioned and i am still dealing with the issue that when the traffic demand is low, it shows that only One variable x is triggered to "1", all others are zeros. When the traffic demands increase higher, if One x variable can not handle the traffics, it should choose more x variables to handle the traffics and meet all predefined requirements. But in the file attached in this post, when traffic demands increase higher, it shows the infeasible problem. For instance, I have tried with various scenarios of traffic demands as in file attached, from T1 to T7. x variable is changed corresponding to different traffic demands, T1 to T6 show different positions of x variable,  but only One decision variable is made. In T7, when I increase the traffic demands, it turns to the infeasible solution, but I would expect that the result should be one more or two more decision variables which are made. I guess there are some possibilities as follows:

   1. Traffic demands T7 does not meet all constraints.  

   2. l(i) s(i) tr(i) do not offer enough possible solutions for cplex solver to solve the problem. 

In case 1, I have checked the solution, it shows like this:

--

solution = 

    yalmiptime: 0.4062

    solvertime: 0.0038

          info: 'Infeasible problem (CPLEX-IBM)'

       problem: 1

--

Therefore, may be the solver could not find out the solution for such a traffic demand.

In case 2, I have also checked my code generator. For each long link l(i), all traffic from sources to destinations can go through this link. For each short link s(i), all traffic from sources to destinations can go through this link. For all traffic paths from a source i to a destination j in the network are considered. 

Do you figure out where the error could be ? How can I change to make more decision variables when I increase the traffic demands like T7? Many thanks.   

[Johan Löfberg]

Impossible to sort this out

You have to create a minimal example where you trivially can illustrate that the solution is sub-optimal (i.e., by manually stating the optimal solution), and then see which constraint is violated by that solution

[Johan Löfberg]

some simple removal of constraints until feasibnility is retained reveals that these two constraints are inconsistent


S = [s(5)<=We, tr(5)>=T(1,6)]
optimize(binmodel(S))

[Johan Löfberg]

...and it is trivial to see why. I leave that as an exercise for you

[Tung]

Thanks Dr. Johan,

    I wrote a minimal example for debugging where the problem is. Here is the file that is attached in this post. The result shows as follows.

% traffic demands 

T=[  

0 5 5 5;

0 0 10 15;

0 0 0 5;

0 15 0 0;];

-- output

reached traffic demands

...

reached a solver

===== result: long links

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 50)

===== result: short links

Bilinear scalar (real, binary, 4 variables, current value : 14)

Bilinear scalar (real, binary, 4 variables, current value : 20)

===== result: traffic links

Bilinear scalar (real, binary, 4 variables, current value : 4)

Polynomial scalar (real, binary, 4 variables, current value : 5)

Polynomial scalar (real, binary, 4 variables, current value : 5)

ans =

     0

Polynomial scalar (real, binary, 4 variables, current value : 10)

Polynomial scalar (real, binary, 4 variables, current value : 15)

ans =

     0

ans =

     0

Bilinear scalar (real, binary, 4 variables, current value : 5)

ans =

     0

Polynomial scalar (real, binary, 4 variables, current value : 15)

ans =

     0

--

It shows that One decision variable x(4) is selected. Long link l(4)=50 meets its constraint requirement. Short links s(1)=14; s(2)=20 also meet their constraint requirements.

I would expect that when I increase the traffic demands with more traffics on sources to destinations, it should increase the number of long links to meet all constraints, for instance, the result shows two decision variables x(4) and x(2) to handle all increased traffic and guarantee to meet all constraint requirements. When the traffic demands increase higher, if x(2) and x(4) can't handle such a higher traffics, it creates more decision variable x, say x(1), so on so forth. 

The infeasible solution should occur when all resources long links capacity can not handle the traffic demands through the network, so it is apparent that there is not enough resources on long links to handle the traffic demands.   

   

However, the problem is that when I increase the traffic demand, it shows the infeasible solution. 

% traffic demands increased T(1,3)=6

T=[  

0 5 6 5;

0 0 10 15;

0 0 0 5;

0 15 0 0;];

-- output

reached a solver

===== result: long links

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

===== result: short links

Bilinear scalar (real, binary, 4 variables, current value : 5)

Bilinear scalar (real, binary, 4 variables, current value : 5)

===== result: traffic links

Bilinear scalar (real, binary, 4 variables, current value : 5)

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

ans =

     0

Polynomial scalar (real, binary, 4 variables, current value : 0)

Polynomial scalar (real, binary, 4 variables, current value : 0)

ans =

     0

ans =

     0

Bilinear scalar (real, binary, 4 variables, current value : 5)

ans =

     0

Polynomial scalar (real, binary, 4 variables, current value : 0)

ans =

     0

--

I have checked the long links, short links, traffic constraints from my code generator, it is correct as I expect. If so, do you guess where the error could be ? Many thanks.

[Johan Löfberg]

As I said, the model is trivially infeasible due to, e.g.,

s(5)<=We, tr(5)>=T(1,6)

Look at those constraints. They are inconsistent

[Tung]

Thanks Dr. Johan,

   I still not figure out where the bugs could be. Here is a smallest network with two subnets (file attached), each subnet consists of two nodes. It still has the same problem, where the result only has either one decision variable (x [1 0 0 0] ) or all zeros. Because the constraints s(i)<=We and tr(i)>=T(x,y) are for the general cases, since we have to cope with different types of traffic demands, so that, when s(i)=0<We, it still meets this constraint's requirement, and tr(i)>=T(x,y) where tr(i) always contains T(x,y), so if T(x,y)=0, it results in tr(i) at least is equal to zero, otherwise, tr(i) is a positive integer. Thus, I think the bugs may not be resided on these constraints. 

   Then I run this program and the result is as follows:

-----

x =

     0

     1

     0

     0

>> s(1)

Bilinear scalar (real, binary, 4 variables, current value : 15)

>> s(2)

Bilinear scalar (real, binary, 4 variables, current value : 20)

>> l(2)

Polynomial scalar (real, binary, 4 variables, current value : 30)

>> 

-----

From this result, we can see that the decision variable is established in x(2), corresponding to the long link l(2) ~ link 14 in the network. This long link meets its requirement where l(2)<=Wn where Wn=50. In addition, short links capacities meet their requirements (s(i)<=We where We=20). So now, long link and short links meet all requirements.

Regarding the traffics, the parameters are as follows:

-----

tr(1) ~ T(1,2): Bilinear scalar (real, binary, 4 variables, current value : 5)

tr(3) ~ T(1,4): Polynomial scalar (real, binary, 4 variables, current value : 20)

tr(6) ~ T(2,4): Polynomial scalar (real, binary, 4 variables, current value : 10)

tr(9) ~ T(3,4): Bilinear scalar (real, binary, 4 variables, current value : 20)

other traffics =0;

-----

This result shows that it is matched to the traffic demands at the input.

If so, the problem could be the variables linearization part. Then, I looked at the binmodel.m, you already have the cases of multiplication of binary decision variables in binmodel file. It's great since I don't have to take care of this linearization. 

However, I am thinking one possibility that it could be the case of this problem. When binmodel tries to linearize the multiplication of binary variables, does it take care of the case when the same multiplication of binary variables should have the same z_linearize? I mean because l(i) s(i) and tr(i) contain many of the same multiplication of binary variables. For example, l(i) has T(2,3)*x(2)*x(4)*x(1) + T(3,4)*x(4)*x(2)*x(1); they should have the same z_linearize(id) like T(2,3)*z_linearize(1) + T(3,4)*z_linearize(1); the same as in s(i) has T(1,2)*x(1)*x(2)*x(4) +...; and this multiplication of binary variables should be replaced by z_linearize(1) as the same in l(i). If binmodel generates another z_linearize rather than z_linearize(1), so, it could be a problem. 

Do you have any ideas or suggestions to solve this problem? Please take a look at this code to see, since everything has been simplified in this code. Many thanks.

       

[Johan Löfberg]

As I said, 1641 is trivially infeasible since those two constraints say

30*x(1)+36*x(2)+33*x(3)+26*x(4)+34*x(5)+31*x(6)+36*x(7)+33*x(8) + more positive monomials <= 20

and

8*x(1)+8*x(2)+16*x(3)+16*x(4)+16*x(5)+16*x(6)+16*x(7)+16*x(8) + more positive monomials >= 8

What does the first constraint tell you? What will happen in the second constraint?

[Tung]

Thanks Dr. Johan,

    I have fixed this problem. Now it works well with various traffic inputs. I have changed a bit in the system model since the traffic demands can not be fixed, it should be tuning for "bargaining" among solutions. I also don't use the binmodel to linearize x variables multiplication, instead, I have used my own code generator for linearization. Thus, my system model looks like as follows.

    Let x be the decision variables (type: binary). 

           z be the linearized variable (type: binary) from multiplication of x binary variables.

           a be the bandwidth allocation in each path from a source to a destination.

           y be the linearized variable (type: continuous) of a*z.

// the lines are for pseudo-algorithm description
0<=x<=1; %binvar
0<=a<=1; %sdpvar
0<=z<=1; %binvar
0<=y<=z;
-(1-z)<=y-a<=(1-z);
where
z=(x)*; % linearization of multiplication of binary variables.
y=a*z; % linearization of a continuous variable and a linearized binary variable. 



This means we have to linearize two times, the first one is for multiplication of binary variables, the second one is for both continuous and binary variables multiplication.

 

Again, thanks a lot for your supports. 

[Tung]

Hi Dr. Johan,

   Do you know how to extract all details from the results of linear scalar values? For instance, after running the solver, I got this result.

      Linear scalar (real, 4 variables, current value : 100)

      Linear scalar (real, 5 variables, current value : 10)

I want to see these in details, something like this:

30*x(1)+36*x(2)+33*x(3)+26*x(4)+34*x(5)+31*x(6)+36*x(7)+33*x(8) + ... <= 20 (this one is copied from your post above)

Because I want to see how the values are assigned into the linear scalar 100 or 10, for example. 

So, how can I get it? Please let me know asap. Thanks a lot.

Best,

Tung

[Johan Löfberg]

expression = f(x)
sdisplay(expression)
see(expression)
getvariables(expression)
getbase(expression)
value(expression)
value(x)

[Tung]

Thanks Dr. Johan,

    I ran this:

          exprs = l(i);

          sdisplay=exprs;

    results:

    -----

Linear scalar (real, 5 variables, current value : 31.25)

    6*internal(30273)+6*internal(30274)+110*internal(30277)+50*internal(30282)+30*internal(30284)

Linear scalar (real, 5 variables, current value : 100)

    6*internal(30275)+6*internal(30276)+110*internal(30280)+30*internal(30281)+30*internal(30284)

    -----

How to see inside the internal(*)? something like 30*x(1)+36*x(2)+33*x(3)+26*x(4)+34*x(5)+31*x(6)+36*x(7);

[Johan Löfberg]

YALMIP sdoes not have a notion of variable name, only indicies. Mst often the symbolic display will not figure out what you name is on a variable, and it will then only tell the index to the variable (internal(n))

yalmip('clear');
sdpvar x y

p = x+3*y
sdisplay(p)

clear x y
sdisplay(p)

>>x+3*y
>>internal(1)+3*internal(2)

[Tung]

Thanks Dr. Johan. Have a great day!

Best, 

[Tung]

Hi Dr. Johan,

   Do the solvers like CPLEX or Gurobi always give the optimal solutions in any problems?

As in my concern, I am trying to prove that these solvers will give the optimal solutions for my problem in networking, however, I have failed to convincing about it due to my heuristic algorithm always gives the better solution on a specific problem in networking. For instance, for a given set of traffic demands, if this set is a light traffic, both solver and heuristic algorithm give the same result, let's say the same number of decision variables x. Assume that we have 1 decision variable corresponding to the given light traffic. However, if a set of high traffic demands is given, solver (cplex) always gives more decision variables than heuristic algorithm, meanwhile, heuristic algorithm only gives lower number of decision variables, and I have checked the constraints in both cases, they both are met the constraints' requirements. I am doubting about it. How do you think about my situation? Please give me some suggestion or advice. Thanks a lot.

Best,

[Johan Löfberg]

Yes, Cplex and gurobi are global solvers (up to the tolerance set by the user, and of course numerical problems which might cause problems for any solver). You can see the claimed optimality in the display.

Cplex does not generate or give any decision variables. All variables are declared by YALMIP, and then cplex simply finds the optimal value. The fact that you have more variables when YALMIP sets up a MILP etc is not strange. YALMIP has to introduce a lot of variables to implement various high-level operators to make them MILP-representable etc.

If your heuristics gives better solutions than cplex/gurobi, something is wrong. You are not implementing the same model, the difference is simply within tolerances, there is a bug somewhere,...Show the model (as small and compact as possible, reproducible etc), together with you claimed better solution, if that is what you are saying is happening.

[Tung]

Thanks Dr. Johan,

   Yes, it might not probably be compared in terms of the modeling perspective. In heuristics, the cost of the flows is minimized, whereas, in the solver the cost of links between subnets is minimized. However, I don't know (1) how to let the solver considering about the cost of the links and (2)how to set a constraint of flows conservation (i.e. guarantee that a traffic from a source to a destination should be conservative, without packets dropped during transmission.). For instance, I would like to add the cost, let's say 100, on links within subnets, and another cost 10000 on links between subnets, each subnet has several nodes to communicate to other nodes either within its subnet or interconnecting to other subnets, it depends on the traffic demands. I would like to let the solver consider that it should choose the low cost for assigning traffics, high cost links should be considered for assigning iff there is no way to go through the low cost links to reach the destination node. I would also like to let the solver consider that when it uses a high cost link, it has to use it for later assigning traffics until reaching link capacity constraint. This link should be set to low cost link in order for assigning traffics as the same as other low cost links.

Please give me some suggestions if possible. Thanks a lot.

Best, 

[Omar-Sayan Karabayev]

but would that cause problems like for example if I take a Jacobian of a polynomial: F

x=sdpvar(2,1)

A = [x(1)^2 2 ; 1 3];

F = A*x;

so F is supposed to be     '2*x(2)+x(1)^3'

                                         'x(1)+3*x(2)'

However,

sdisplay(F)

ans= 

    '2*x(2)+internal(1)^2*x(1)'

    'x(1)+3*x(2)'

sdisplay(jacobian(F,x))

ans = 

    'internal(1)^2'    '2'

    '1'                '3'

where, internal(1) in F is supposed to be x(1) , so i should've gotten 3*x(1)^2 (or at least 3*internal(1)^2) in the jacobian at position (1,1)

I am running this in an .m file. However if I do this in command window it solves it right. 

Help would be very much appreciated

[Johan Löfberg]

Using an old version of YALMIP? Works here

>> x=sdpvar(2,1)

A = [x(1)^2 2 ; 1 3];

F = A*x;

Linear matrix variable 2x1 (full, real, 2 variables)

>> sdisplay(F)

ans = 

    '2*x(2)+x(1)^3'

    'x(1)+3*x(2)'

[Omar-Sayan Karabayev]

Thank you for the quick reply!

I'm using the current version. It works fine when I run these commands in the command window like you, but I need to have all that code in a function inside an m-file which I call from the command window. 

[Johan Löfberg]

Works here

>> testsdisplay

ans = 

    '2*x(2)+x(1)^3'

    'x(1)+3*x(2)'

function testdisplay

x=sdpvar(2,1);

A = [x(1)^2 2 ; 1 3];

F = A*x;

sdisplay(F)

[Johan Löfberg]

but note though that, as discussed in this thread, sdisplay is not guaranteed to work as it is a hack. yalmip does not have any notion of variable names, and sometimes it fails in figuring out what you have called the variables

[Omar-Sayan Karabayev]

Lots of thanks for your help!

[Omar-Sayan Karabayev]

Actually, one more thing, Dr. Lofberg

I ran the function you wrote and it worked fine, BUT, Could you try this:

in the testsdisplay function delete the declaration of the A matrix, and declare A in the command window instead, then call the function giving matrix A as input, so:

function testsdisplay(A)
x = sdpvar(2,1);


F = A*x;
sdisplay(F)
J = jacobian(F,x);
sdisplay(J)

>>x =sdpvar(2,1);

>>A = [x(1)^2, 2;4 5];

>> testsdisplay(A)

ans = 

    '2*x(2)+internal(1)^2*x(1)'

    '3*x(1)+3*x(2)'

ans = 

    'internal(1)^2'    '2'

    '3'                '3'

If the error is only in the displaying the variables, then it's fine. I'm more worried if this will affect when I'm actually solving SOS problems with constraints by solvers later.

Thank you again!

[Johan Löfberg]

That's entirely expected, as x defined outside the function is another variable (you've defined different vectors x)

[Omar-Sayan Karabayev]

OK, now it makes sense, thanks!