Model creation failed (Model not available in constraint #19376 at level 1)

[仝立杨]

Hi Dr.Johan,

    I encountered this problem when using yalmip to solve a bi-objective optimization problem. I used the episilon-constraint method. This constraint exists in Part A, but I can run Part A alone. and the results were obtained, although it took a whole day. BTW, the model construction time is far beyond my imagination. Therefore, I am not quite sure why such a problem occurs after integration. I know that the problem I want to solve is very tricky. I have posted complete codes of Part A alone and the integrated program on the forum. I hope you can give some suggestions for solving this problem. Thank you very much.

[Johan Löfberg]

https://yalmip.github.io/squareroots/

and I hope you realize you are constructing a very very (very!) hard optimization problem. Basically, you cannot construct a more complicated class of optimization problems

also, most likely a quadratic model will perform faster in the global solver, instead of a == sqrt(b^2+c^2) you write a^2 == b^2 + c^2, and I would not be surprised if cos(a)*b == c is better than a = acos(c/b) (avoids singularties in division and acos)

[仝立杨]

Thanks for your quick reply, Dr,Johan. I have taken you advices, it is creating model now, hope it can be successful created and solved, no matter the time consuming. 

Honestly, this program makes me mad, but i have to. I firstly tried several slovers with Python, like Gurobi, SCIP. but this problem cannot be solved. Maybe it's me, I don't know how to deal with the trig function in Python, It seems like you cannot use trig function in Python with a variable. and also the min and max. until i met Yalmip a couple days ago. Yalmip is so simple and convinient, amazing especiall for me because i don't know much complicated math knowledege about optimization. 

So you think this problem can be solved, right? even though it is very very complicated.

BTW, my supervisor suggested me to solve this problem with bi-level programming, but i don't think it can work out, because it involves too many binary variables so the reformulation seems impossible. am i right? or maybe it can be solved by genetic algorithm.

Again, thank you very much 

[Johan Löfberg]

Python is not a solver, so it makes no sense to say ", It seems like you cannot use trig function in Python with a variable"

Using acos on a decision variable, means you must use a general nonlinear solver (i.e. definitely not gurobi). Using that, and binary variables, means you must use a nonlinear nonconvex integer solver, i.e. the worst possible. YALMIPs global nonlinear solver, with integers, can be expected to work well on a few (like 1-10) variables. Sometimes you are lucky and it works on larger problems as the model looks complicated but it hides a simple structure. You model is theoretically intractable just based on the number of binary variables you have. It might work, but then you are just lucky

Your objective looks like a complicated way of saying Objective = min(sum(T)) (which will be implemented roughly as you've done it as it is nonconvex as you minimize a concave function, but why do that manually when you implement everything else using high-level operators)

Your best bet is to try to get rid of acos and quadratic stuff and get it down to a linear representation, possible using pwa approximations based on sos2-structures etc. Sure, you can try some random solver, but then you are outside the world of yalmip, and it's all down to luck and heuristics

[Johan Löfberg]

and if you use YALMIP, you must really vectorize your code because it awfully slow now with all the referencing

This

Constraints = [Constraints,d(i,j,m,n)== min(min(min(d1(i,j,m,n,1,1),d1(i,j,m,n,1,2)),min(d1(i,j,m,n,1,3),d1(i,j,m,n,1,4))),min(min(d1(i,j,m,n,2,1),d1(i,j,m,n,2,2)),min(d1(i,j,m,n,2,3),d1(i,j,m,n,2,4)))),min(min(min(d1(i,j,m,n,3,1),d1(i,j,m,n,3,2)),min(d1(i,j,m,n,3,3),d1(i,j,m,n,3,4))),min(min(d1(i,j,m,n,4,1),d1(i,j,m,n,4,2)),min(d1(i,j,m,n,4,3),d1(i,j,m,n,4,4))))];

takes forever, but it appears to be misplaced. It's inside a loop with indicies p  & q, but those are not used here

also, it is not clear if it really is what you intend it to be. The lhs is, of course, a scalar, but the rhs is a 1x2 vector

>> [min(min(min(d1(i,j,m,n,1,1),d1(i,j,m,n,1,2)),min(d1(i,j,m,n,1,3),d1(i,j,m,n,1,4))),min(min(d1(i,j,m,n,2,1),d1(i,j,m,n,2,2)),min(d1(i,j,m,n,2,3),d1(i,j,m,n,2,4)))),min(min(min(d1(i,j,m,n,3,1),d1(i,j,m,n,3,2)),min(d1(i,j,m,n,3,3),d1(i,j,m,n,3,4))),min(min(d1(i,j,m,n,4,1),d1(i,j,m,n,4,2)),min(d1(i,j,m,n,4,3),d1(i,j,m,n,4,4))))]

Nonlinear matrix variable 1x2 (full, real, 2 variables)

Coeffiecient range: 1 to 1

[仝立杨]

Thank you, Dr. Johan. You really help me a lot. I will take your advices seriously in the next day or two. really really appreciate your kind help !!

As for the constraint you posed, yes, it's a mistake. actually I want to get a minimum from the sixteen variables. and  i tried d(i,j,m,n)=min(min(d1(i,j,m,n,1:4,1:4))) before, but it doesn't work with a error message saying no defination of function min for the parameter with the type ndsdpvar. but when i try something simple like d(i) = d1(i,1:4), it can work. I am not sure if it is because the diamension of this variable is too high, so i adopted the stupid way to find a minimum, at least it can run.

[Johan Löfberg]

simply do  e = d1(i,j,m,n,1:4,1:4);d(i,j,m,n) = min(e(:))

[仝立杨]

Hi Dr. Johan,

I am a little confused to code  e = d1(i,j,m,n,1:4,1:4);d(i,j,m,n) = min(e(:)) you suggested. I mean how can i put this in my code correctly.

i added it this way:  Constraint = [Constraint, e == d1(i,j,m,n,1:4,1:4), d(i,j,m,n) = min(e(:)) ], but it seems weired. because e has no index and it can't run successfully. 

Actually, the logic is if a[m,n]=1 (a[m,n] is predetermined) and then if z[i,m]=1 and z[j,n]=1 (they are binvar) then we caculate the d1[i,j,m,n,1:4,1:4]. and we need to to find out the minimum of  d1[i,j,m,n,1:4,1:4], so how can i correct the below code? the BCX and BCY are predetermined, too.

for i = 1 : 11

    for j = i : 11

        if i~=j

            for m = 1 : 7

              for n = 1 : 7

                  if a(m,n)==1

                      for p = 1 : 4

                         for q  = 1 : 4

                            Constraints = [Constraints, 0 <= d1(i,j,m,n,p,q) <= (W+B)];

                            Constraints = [Constraints, implies(z(i,m)+z(j,n)==2, d1(i,j,m,n,p,q) == abs(cx(i)-BCX(m,p))+abs(cy(i)-BCY(m,p))+abs(cx(j)-BCX(n,q))+abs(cy(j)-BCY(n,q))+abs(BCX(m,p)-BCX(n,q))+abs(BCY(m,p)-BCY(n,q)))];

                      Constraint = [Constraint, implies(z(i,m)+z(j,n)==2, [e == d1(i,j,m,n,1:4,1:4), d(i,j,m,n) = min(e(:))] ]

                         end

                      end

                  end

              end

           end

      end

  end

end

[Johan Löfberg]

You seem to fail to understand the difference between assignment and equalities as you have flipped them around completely

The idea is that you simplify by first *assigning* the temporary e = d1(i,j,m,n,1:4,1:4) in the loops, and then use that expression in an equality in the constraints d(i,j,m,n) == min(e(:)). Using an assignment inside a constraint list does not work

In other words

for i = 

 e = x(i,:)

 Model = [Model, y(i) == min(e)]

end

is 100% equivalent to

for i = 

 Model = [Model, y(i) == min(x(i,:))]

end

but the version without a temporary variable e is not possible in your case since you must reshape it to become a standard 2d-array and not an ndsdpvar

[仝立杨]

Thank you, Dr. Johan. i kind of get it. so i write the codes below:

for i = 1 : 11

    for j = i : 11

        if i~=j

            for m = 1 : 7

              for n = 1 : 7

                  if a(m,n)==1

                      for p = 1 : 4

                         for q  = 1 : 4

                            Constraints = [Constraints, 0 <= d1(i,j,m,n,p,q) <= (W+B)];

                            Constraints = [Constraints, implies(z(i,m)+z(j,n)==2, d1(i,j,m,n,p,q) == abs(cx(i)-BCX(m,p))+abs(cy(i)-BCY(m,p))+abs(cx(j)-BCX(n,q))+abs(cy(j)-BCY(n,q))+abs(BCX(m,p)-BCX(n,q))+abs(BCY(m,p)-BCY(n,q)))];

                         end

                      end

                  end

              end

           end

      end

  end

end

e = sdpvar(4,4);

for i = 1 : 11

    for j = i : 11

        if i~=j

            for m = 1 : 7

              for n = 1 : 7

                  if a(m,n)==1

                     e=d1(i,j,m,n,1:4,1:4);

                     Constraints = [Constraints, d(i,j,m,n) = min(e(:))];

                  end

              end

            end

        end

    end

end

yalmip is creating model now, but i am a little worried that in my case the d1[i,j,m,n,p,q] is only cacuated when z[i,m]+z[j,n]=2. but it seems not considering this condition in the second part. will this casue some problems?

[仝立杨]

there is a typo in the above codes

Constraints = [Constraints, d(i,j,m,n) == min(e(:))]; 

[Johan Löfberg]

no reason t odefine e = sdpvar(4,4), since you redefine it later everytime in the loop by assigning it to be e=d1(i,j,m,n,1:4,1:4);

Not sure what you are talking about in the second issue. Calculated? As it is now, d1 is constrained to be equal to abs()... when z[i,m]+z[j,n]=2. It is also constrained to satisfy min(d1(i,j,m,n,1:4,1:4)) = d(i,j,m,n) when i~=j and a(m,n)=1

[仝立杨]

Ooo, i see, have removed it.

yes, it was 'caculated', sorry. 

i meant  the d1[i,j,m,n,p,q] is only calculated when z[i,m]+z[j,n]=2 as shown in the upper loop. 

and the d(i,j,m,n) should be the minimum of caculated d1(i,j,m,n,1:4,1:4). Actually, d1(i,j,m,n,p,q)  is not available for all m,n. for some m and n, there might be no value of d1(i,j,m,n,p,q), and it should return NaN according to my understanding.

but in the below loop to find the minimum, it seems not considering this condition that just assigns e when z[i,m]+z[j,n]=2. will this casue problems? or how to add this condition to assignment. Thank you 

[Johan Löfberg]

If a variables isn't involved in any constraint or in the objective, it has no value and will thus simply be nan if you apply value, just like you get nan if you define a variable and directly apply value

However, implies(d == 1, x == 5) does not mean x will be nan if d = 0. x is still part of the model, but is simply not constrained when d = 0, so the solver can return what ever it wants (unless it is constrained by other constraints of course)

[仝立杨]

Hi Dr. Johan,

Thanks for your previous kind help. I have simplified my codes, but I encountered a new problem. The result  was like:

Hmm, something went wrong!

ans =

Model creation failed (Model not available in objective at level 4)

ans =

Model creation failed

The code is attached. I have checked it. If I remove the Tower from my Objective, then it's ok to run. Could you help me with this? Thank you.

[Johan Löfberg]

Looks like you try to aim for a convex SOCP model in an obviously nonconvex model

https://yalmip.github.io/squareroots/

and don't do sqrt(stuff^2+bla^2)^2, you're just making the model overly complicated when you just as well could have written stuff^2+bla^2 and thus arrived at a lot less complex model representation

and the min/max/abs operator probably won't work either, as the arguments aren't convex and thus YALMIP will use MILP representation which is unnecessarily complicated. Model those manually by graphs instead, when possible

[仝立杨]

Thanks for your quick reply, Dr. Johan. 

Regarding the sqrt thing, do mean to write in this following way? 

dst^2   = (cx(1)-Tx)^2+(cy(1)-Ty)^2;

dsd^2  = (cx(1)-Demand(j,1))^2+(cy(1)-Demand(j,2))^2;

dtd^2   = (Tx-Demand(j,1))^2+(Ty-Demand(j,2))^2;

If so,  I got a problem that says the left hand expression is a invalid target for assignment....

[Johan Löfberg]

Of course, that's simply not valid matlab. It doesn't even make sense from any programming language to try to write *assignment* x^2 = 5 (what value should be assigned, 5 or -5?)

I'm saying you shouldn't write

dsd  = sqrt((cx(1)-Demand(j,1))^2+(cy(1)-Demand(j,2))^2);        
Tw   = acos((dsd^2...;

since that will introduce sqrt operators, and square operators on two levels, when you simply can write

whatever = ((cx(1)-Demand(j,1))^2+(cy(1)-Demand(j,2))^2);        

Tw   = acos((whatever...;

Having said that, no matter how you model this, you get an extremely bad nonlinear nonconvex mixed-integer model, and most likely we are just comparing very very very very hard with very very very hard

Also note that there are bugs in the acos implementation in the master branch of YALMIP, so you must install the develop version if you want to use that

[仝立杨]

Haha. How silly i am, ok I get it.  Thank you, Dr. Johan. And I have installed the develop version.  

Now it shows the problem is  infesible, I am checking it according to the Debugging Infeasible Models.