Solver not applicable (gurobi does not support polynomial constraints)

[吴思仁]

I have two models named M1 and M2. 

M1 can be solved successfully with solver Gurobi, while M2 returns "Solver not applicable (gurobi does not support polynomial constraints)".

I checked all the constraints of M1 and M2 and found that they all contain bilinear and polynomial constraints. Why do they give different result when solving with the same solver?

[Johan Löfberg]

As it says, your M2 models has polynomial constraints. Gurobi does not support polynomial constraints but only bilinear/quadratic, hence M1 does not include any polynomial constraints.

[吴思仁]

But M1 also includes polynomial constraints:

```

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|      ID|                                 Constraint|   Coefficient range|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|      #1|     Element-wise inequality (bilinear) 1x1|          100 to 150|
|      #2|    Element-wise inequality (bilinear) 28x1|              1 to 1|
|      #3|    Element-wise inequality (bilinear) 28x1|              1 to 1|
|      #4|   Element-wise inequality (polynomial) 7x1|              1 to 1|
|      #5|   Element-wise inequality (polynomial) 7x1|              1 to 1|
|      #6|     Element-wise inequality (bilinear) 1x1|           10 to 100|
|      #7|     Element-wise inequality (bilinear) 1x1|           10 to 100|

...

```

[Johan Löfberg]

You would have to supply a reproducble example as that should not be possible (unless you are using an optimizer setup with e.g. yx^2 with a parameter y and decisoon variable x, i.e. polynomial as written, but quadratic when instantiated for parameter values)

[吴思仁]

Thank you very much for your kind help, Johan. 

Here I would like to provide more information about the two models. M1 and M2 are both robust optimization models described with YALMIP. I carefully analyzed their polynomial constraints and found out their possible difference:

 - M1 includes constraints of form [z == k * (x + k * y)] where k is uncertain variable and x, y, z are all decision variables.

 - M2 includes constraints of form [z == k * (x + k * y), k == w * d] where d is uncertain variable and x, y, z, w, k are all decision variables.

May I ask you:

1. Are these two constraints mentioned above the reason why M1 and M2 contain polynomial constraints?

2. For robust optimization, the uncertain variable k in M1 is finally treated as some parameter values, so the model M1 can be solved normally; while M2 contains term w^2*d^2*y, which leads to a polynomial constraint with respect to w^2*y, although d^2 is treated as a parameter. Am I understanding this right?

And 3. Can you tell me how YALMIP deal with a constraint including uncertain variables and which kinds of constraints can robust optimization of YALMIP handle?

[Johan Löfberg]

To begin with, that is not a reproducible example, so original question cannot be investigated.

However, the model is completely flawed as you have uncertain variables in equality constraints. If you think that is sensible, please give me values of w and d, such that w*d == k for any 0<=k<=1 

But then again, the last paragraph makes everything confusing, as you then talk about parameters (as in optimizer objects) instead of uncertain variables. If you actually had an uncertain model, YALMIP should fail during derivation of the robust model, and should most likely never even make it to try to use gurobi ,while it definitely works in an optimizer construct with k as parameter as there is nothing polynomial here (there is no w^2d^2 term, unless you are confused about having an equality in the model, or if you actually mean that you have assigned k to be w*d, but then why even talk about parameter varaibles and/or uncertain variables)

Hence: unclear what you are doing in this non-reproducible example. Supply reproducible example

[吴思仁]

Yes,  I actually mean that I have assigned k to be w*d, and  assigned z to be k * (x + k * y), and there are other constraints on z such as 0 <= z <= 1 (Actually the uncertain variables didn't only appear in equality constraints). 

Please forgive me for not giving a complete reproducible example, because the original model is too complicated for me. Actually, in a nutshell, I'm trying to understand why both robust optimization models M1 and M2 contain polynomial constraints, but M1 can be solved successfully but M2 cannot. As I said earlier, the polynomial constraints contained in M1 and M2 differ only in (here I try to express it more rigorously again) :

 - M1: [0 <= k * (x + k * y) <= 1] where k is uncertain variable that has a range such as 0<=k<=0.1 and x, y are all decision variables.

 - M2: [0 <= w * d * (x +  w * d * y) <= 1] where d is uncertain variable that has a range such as 0<=d<=0.1 and x, y, w are all decision variables.

Now the only question I want to ask is, can the above information I provided explain why both M1 and M2 contain polynomial constraints but M1 can be successfully solved but M2 cannot?

[Johan Löfberg]

You are not being consistent in your description. You are first saying that k has been *assigned* to be w*d and  z to be k * (x + k * y), i.e. the models

sdpvar x y w
k = w*d;
z = k*(x+k*y);
M1 = [0 <= k*(x+k*y) <= 1]
M2 = [0 <= w*d*(x+w*d*y) <= 1]

i.e. what you have defined are the same models

M1 = [0 <= x*w*d+y*w^2*d^2 <= 1]
M2 = [0 <= x*w*d+y*w^2*d^2 <= 1]

i.e. the same thing, and that model should not work with Gurobi with x, y, w, d as decision variables. The only time gurobi would work is if you used them in an optimizer with w and d as parameters as it then would instantiate as quadratic models

However, then you start talking about robust stuff and uncertainty. you say k and d are uncertain variables, i.e. you would have 

sdpvar x y k d

M1 = [0 <= k*(x+k*y) <= 1, uncertain(k), 0 <= k <= 0.1]

M2 = [0 <= w * d * (x +  w * d * y) <= 1,uncertain(d), 0 <= d <= 0.1] 

but that will run into troubles as quadratic stuff in the uncertainty isn't supported by the robustness engine, so effectively this will lead to the engine trying to eliminate k and d from the models, which means adding [y==0] to M1 and [w^2*y==0] to M2. The final robustified M1 will easily be solved by Gurobi, M2 will not as it will involve a cubic constraint

[吴思仁]

Yes, what I'm trying to say is exactly what you said in the second half. To be precise, the two models should look like this: 

sdpvar x y w k d

M1 = [0 <= k*(x+k*y) <= 1, uncertain(k), 0 <= k <= 0.1]

M2 = [0 <= w * d * (x +  w * d * y) <= 1,uncertain(d), 0 <= d <= 0.1] 

Now I want to make it clear:

1. Will both models M1 and M2 be reported by Yalmip to contain polynomial constraints?

2. Yalmip's robustness engine can accept and robustify models M1 and M2, but only the robustified M1 can be solved by Gurobi successfully while M2 will fail for the polynomial constraint?

3. The solution log of real M1 is as follows, can you help me judge whether it can be considered as a successful solution:

***** Starting YALMIP robustification module. *********************
 - Detected 48 uncertain variables
 - Detected 48 independent group(s) of uncertain variables
 - Complicating terms in w encountered. Trying to eliminate by forcing some decision variables to 0
 - Complicating terms in w encountered. Trying to eliminate by forcing some decision variables to 0
 - Eliminating uncertainty using explicit maximization of inf-norm
***** Derivation of robust counterpart done ***********************
Academic license - for non-commercial use only - expires 20xx-xx-xx
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 47653 rows, 21106 columns and 343533 nonzeros
Model fingerprint: xxxxxxx
Variable types: 21050 continuous, 56 integer (56 binary)
Coefficient statistics:
  Matrix range     [4e-05, 1e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [3e-14, 4e+03]
Presolve removed 35203 rows and 15560 columns
Presolve time: 3.57s
Presolved: 12450 rows, 5546 columns, 62908 nonzeros
Variable types: 5524 continuous, 22 integer (22 binary)

Root simplex log...

Iteration    Objective       Primal Inf.    Dual Inf.      Time
    7170    1.0557720e+04   7.604274e+05   0.000000e+00      5s
   17840    1.0703594e+04   5.011939e+06   0.000000e+00     10s
   21690    1.0741404e+04   5.357545e+05   0.000000e+00     15s
   25330    1.0779904e+04   1.569191e+06   0.000000e+00     20s
   28690    1.0803012e+04   7.330740e+06   0.000000e+00     25s
   32190    1.0835664e+04   5.400681e+05   0.000000e+00     30s
   36180    1.0868502e+04   1.365175e+05   0.000000e+00     35s
   40693    1.4235013e+04   5.781053e+07   0.000000e+00     40s
Warning: Markowitz tolerance tightened to 0.03125
   44201    9.0072367e+04   1.554219e+07   0.000000e+00     45s
   47701    9.5070696e+04   1.028131e+07   0.000000e+00     50s
   50926    9.9985530e+04   7.728223e+07   0.000000e+00     55s
   54483    1.0145411e+05   1.686782e+06   0.000000e+00     60s
Warning: Markowitz tolerance tightened to 0.125
   58858    1.0218533e+05   7.344124e+05   0.000000e+00     65s
   65421    1.0295841e+05   5.031727e+05   0.000000e+00     70s
   67601    1.0159633e+05   0.000000e+00   0.000000e+00     71s

Root relaxation: objective 1.015963e+05, 67601 iterations, 67.25 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 101596.329    0   17          - 101596.329      -     -   71s
H    0     0                    116883.95528 101596.329  13.1%     -   75s
     0     0 101596.329    0   17 116883.955 101596.329  13.1%     -   78s
     0     0 101596.877    0   17 116883.955 101596.877  13.1%     -   84s
     0     0 101596.877    0   17 116883.955 101596.877  13.1%     -   84s
     0     0 101596.877    0   17 116883.955 101596.877  13.1%     -   94s
     0     2 101596.877    0   17 116883.955 101596.877  13.1%     -   94s
     1     4 105364.479    1   19 116883.955 101596.877  13.1% 53752  206s
     3     4 106279.815    2   17 116883.955 101596.877  13.1% 61197  305s
     7     4 108790.046    3   18 116883.955 101596.877  13.1% 54390  360s
    11     6 109945.422    4   16 116883.955 107516.341  8.01% 42710  367s
    15     6 110657.957    5   14 116883.955 107997.764  7.60% 32513  406s
    19     6 113969.874    6    9 116883.955 107997.764  7.60% 29342  410s
*   20     6               6    111210.68895 107997.764  2.89% 27877  410s
*   25     0               6    110973.10238 110973.102  0.00% 22656  412s

Cutting planes:
  Implied bound: 3
  RLT: 8

Explored 26 nodes (634199 simplex iterations) in 412.90 seconds
Thread count was 8 (of 8 available processors)

Solution count 3: 110973 111211 116884

Optimal solution found (tolerance 1.00e-04)
Best objective 1.109731023838e+05, best bound 1.109731023838e+05, gap 0.0000%

[Johan Löfberg]

M1 is bilinear in the decision variables. It has a nonlinearity w.r.t k uncertainty, which YALMIP does not have a tactic to robustify except forcing decision variables to take values such that k is not seen, i.e. the extremely conservative robust counterpart of  [0 <= k*x+k^2*y <= 1] is to force y==0 leaving it to robustify  [y == 0, 0 <= k*x <= 1] leading to  [y == 0, 0 <= 0.1*x <= 1] 

M2 is cubic in the decision variables. It has a nonlinearity w.r.t d uncertainty, which YALMIP does not have a tactic to robustify except forcing decision variables to take values such that d is not seen, i.e. the extremely conservative robust counterpart of  0 <= w*d*x + w^2*d^2*y <= 1 is to force w^2*y==0 leaving it to robustify  [w^2*y == 0, 0 <= w'*d*x <= 1] leading to  [w^2*y == 0, 0 <= 0.1*w*x <= 1] . This is cubic thus gurobi not applicable

[吴思仁]

“ the extremely conservative robust counterpart of  [0 <= k*x+k^2*y <= 1] is to force y==0 ” 

Why can we robustify [0 <= k*x+k^2*y <= 1] through forcing decision variables y==0?  What's the reason for doing this? Is this a last resort approximation, or a strictly mathematical trick? (Please forgive my lack of a solid foundation for robust optimization)

[Johan Löfberg]

You want   [0 <= k*x+k^2*y <= 1] to hold for all 0 <= k<=0.1

with y==0 the term k^2 is eliminated, and robustification of the remaining linear part is simple

it is last resort that is done when there are no filters available for the model, in this case for nonlinear dependence of uncertainty. Brutally conservative, so if important you would have to either make an outer approximation using sampling, or attack it using a sum-of-squares constraint, i.e. a semidefinite programming model (or some other method to find the exact robust counter-part)

[吴思仁]

Oh, I see...Thanks. And are there any literature references for this brutal approach?

[吴思仁]

Also, regarding the last two methods you mentioned (1. make an outer approximation using sampling, 2. attack it using a sum-of-squares constraint), could you please provide some details or give some references? Thanks!!!

在2022年1月24日星期一 UTC+8 17:09:38<Johan Löfberg> 写道：

[Johan Löfberg]

there would not be any references on that, as it is the trivial "I give up and have no idea left" approach. You should not use it, but come up with the exact robust counter-part or at least a good approximation, which shouldn't be too hard

[Johan Löfberg]

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

sdpvar x y
K = [];
for i = 1:1000
    k = 0.1*rand(1);
    K = [K, [0 <= k*x+k^2*y <= 1]];
end
subplot(1,2,1);plot(K)

sdpvar x y k
p1 = 1-(k*x+k^2*y);
p2 = k*x+k^2*y;
g = (0.0025-(k-0.05)^2)
[s1,c1] = polynomial(k,2);
[s2,c2] = polynomial(k,2);

M = [sos(p1 - s1*g),
     sos(p2 - s2*g),
     sos(s1),sos(s2)]
F = compilesos(M,[],sdpsettings('sos.model','image'),[x;y;c1;c2]);
subplot(1,2,2);plot(F,[x;y])

[吴思仁]

Wonderful! Thanks!!!

[吴思仁]

For the SOS method, why does the following scheme not work?

sdpvar x y k
p1 = 1-(k*x+k^2*y);
p2 = k*x+k^2*y;

g = [k;0.1-k]
[s11,c11] = polynomial(k,2);
[s12,c12] = polynomial(k,2);
[s21,c21] = polynomial(k,2);
[s22,c22] = polynomial(k,2);

M = [sos(p1 - [s11 s12]*g),
     sos(p2 - [s21 s22]*g),
     sos(s11),sos(s12),sos(s21),sos(s22)]
F = compilesos(M,[],sdpsettings('sos.model','image'),[x;y;c11;c12;c21;c22]);
plot(F,[x;y])

在2022年1月25日星期二 UTC+8 17:57:36<Johan Löfberg> 写道：

[Johan Löfberg]

Don't know

[吴思仁]

And why does the following code draw a weird graph? Just change g to g = (0.0025-abs(k-0.05)^2). (Isn't it equivalent to g = (0.0025-(k-0.05)^2)? )

p1 = 1-(k*x+k^2*y);
p2 = k*x+k^2*y;

g = (0.0025-abs(k-0.05)^2)

[s1,c1] = polynomial(k,2);
[s2,c2] = polynomial(k,2);

M = [sos(p1 - s1*g),
     sos(p2 - s2*g),
     sos(s1),sos(s2)]
F = compilesos(M,[],sdpsettings('sos.model','image'),[x;y;c1;c2]);

plot(F,[x;y])

[Johan Löfberg]

It is not a polynomial

[吴思仁]

oh...yes