Not convex with MOSEK solver

GAMS
1

Dear all,

I am a beginner in the GAMS.

My optimization problem is a Mixed Integer Cone Programming problem and trying to solve with the MOSEK solver. It is showing like below after running code

MOSEK 32.2.0 rc62c018 Released Aug 26, 2020 WEI x86 64bit/MS Window

M O S E K version 9.2.18 (Build date: 2020-8-11 14:59:10)
Copyright (C) MOSEK ApS, Fruebjergvej 3, Box 16
DK-2100 Copenhagen, Denmark
https://www.mosek.com

GAMS/MOSEK Link license detected.

Setting up initial integer solution.
Problem
Name :
Objective sense : min
Type : QCQO (quadratically constrained optimization problem)
Constraints : 8836
Cones : 0
Scalar variables : 8904
Matrix variables : 0
Integer variables : 240

Optimizer started.
Quadratic to conic reformulation started.
The constraint ‘eq3(1,2,t1)’(3072) is not convex. Q should be negative semidefinite for a constraint with finite lower bound. (1294)
Quadratic to conic reformulation terminated. Time: 0.00
Optimizer terminated. Time: 0.00


Integer solution solution summary
Problem status : UNKNOWN
Solution status : UNKNOWN
Primal. obj: 0.0000000000e+00 nrm: 1e+03 Viol. con: 1e+03 var: 0e+00 itg: 0e+00

Return code - 1294 [MSK_RES_ERR_CON_Q_NOT_NSD]: The quadratic constraint matrix is not NSD.

The same optimization problem is running with the CPLEX and GUROBI solvers successfully. I would like to know, why MOSEK solver is showing like that? Please help me in this regard.

2

I am not entirely sure, but that’s because CPLEX and Gurobi can solve nonconvex quadratic programs, whereas MOSEK cannot.

Atharv

3

Mosek is pretty pure in terms of how your cones need to look like. Cplex and Gurobi try to transform your constraints and see if they can reformulate this into a convex cone. Without seeing the constraints/model it is hard to say, but the message from MOSEK about PSD indicates that MOSEK interpreted your constraint as a regular quadratic constraint and not as a SOCP constraint.

-Michael

4

Dear sir,

The following is my code of SOCP based optimization problem,

$Title Loadflow of IEEE 33 bus distribution system

Set
i ‘System Buses’ /1*33/
slack(i) ‘Substation Bus’ /1/;

set non_slack(i);
non_slack(i) = NOT slack(i);
Alias (i,j);

Table LN(i,j,*) ‘Line data’
r x limit
1 .2 0.0922 0.0470 4600
2 .3 0.4930 0.2511 4100
3 .4 0.3660 0.1864 2900
4 .5 0.3811 0.1941 2900
5 .6 0.8190 0.7070 2900
6 .7 0.1872 0.6188 1500
7 .8 0.7114 0.2351 1050
8 .9 1.0300 0.7400 1050
9 .10 1.0440 0.7400 1050
10.11 0.1966 0.0650 1050
11.12 0.3744 0.1298 1050
12.13 1.4680 1.1550 500
13.14 0.5416 0.7129 450
14.15 0.5910 0.5260 300
15.16 0.7463 0.5450 250
16.17 1.2890 1.7210 250
17.18 0.7320 0.5740 100
2 .19 0.1640 0.1565 500
19.20 1.5042 1.3554 500
20.21 0.4095 0.4784 210
21.22 0.7089 0.9373 110
3 .23 0.4512 0.3083 1050
23.24 0.8980 0.7091 1050
24.25 0.8960 0.7011 500
6 .26 0.2030 0.1034 1500
26.27 0.2842 0.1447 1500
27.28 1.0590 0.9337 1500
28.29 0.8042 0.7006 1500
29.30 0.5075 0.2585 1500
30.31 0.9744 0.9630 500
31.32 0.3105 0.3619 500
32.33 0.3410 0.5302 100;

Scalar
Sbase / 1000 /
Vbase / 12.66 /;

LN(i,j,‘r’)=LN(i,j,‘r’)(1/Vbase**2);
LN(i,j,‘x’)=LN(i,j,‘x’)(1/Vbase**2);


Table BD(i,*) ‘demands of each bus in MW’
Pd Qd
1 0 0
2 100 60
3 90 40
4 120 80
5 60 30
6 60 20
7 200 100
8 200 100
9 60 20
10 60 20
11 45 30
12 60 35
13 60 35
14 120 80
15 60 10
16 60 20
17 60 20
18 90 40
19 90 40
20 90 40
21 90 40
22 90 40
23 90 50
24 420 200
25 420 200
26 60 25
27 60 25
28 60 20
29 120 70
30 200 600
31 150 70
32 210 100
33 60 40;


LN(i,j,‘x’)(LN(i,j,'x')=0) = LN(j,i,'x'); LN(i,j,'r')(LN(i,j,‘r’)=0) = LN(j,i,‘r’);
LN(i,j,‘Limit’)$(LN(i,j,‘Limit’)=0) = LN(j,i,‘Limit’);
LN(i,j,‘z’)LN(i,j,'Limit') = sqrt(sqr(LN(i,j,'x')) + sqr(LN(i,j,'r'))); LN(j,i,'z')(LN(i,j,‘z’)=0) = LN(i,j,‘z’);
LN(i,j,‘th’)(LN(i,j,'Limit') and LN(i,j,'x') and LN(i,j,'r')) = arctan(LN(i,j,'x')/(LN(i,j,'r'))); LN(i,j,'th')(LN(i,j,‘Limit’) and LN(i,j,‘x’) and LN(i,j,‘r’)=0) = pi/2;
LN(i,j,‘th’)$(LN(i,j,‘Limit’) and LN(i,j,‘r’) and LN(i,j,‘x’)=0) = 0;
LN(j,i,‘th’)$LN(i,j,‘Limit’) = LN(i,j,‘th’);
LN(i,j,‘Gij’)$LN(i,j,‘Limit’) = (LN(i,j,‘r’))/(sqr(LN(i,j,‘x’)) + sqr(LN(i,j,‘r’)));
LN(i,j,‘Bij’)LN(i,j,'Limit') = -((LN(i,j,'x'))/(sqr(LN(i,j,'x')) + sqr(LN(i,j,'r')))); LN(j,i,'Gij')(LN(i,j,‘Gij’)=0) = LN(i,j,‘Gij’);
LN(j,i,‘Bij’)$(LN(i,j,‘Bij’)=0) = LN(i,j,‘Bij’);

Parameter cx(i,j);
cx(i,j)(LN(i,j,'limit') and LN(j,i,'limit')) = 1; cx(i,j)(cx(j,i)) = 1;

Variable OF, Pij(i,j), Qij(i,j), V2(i), Rij(i,j), Uij(i,j), Pss(i), Qss(i);
Equation eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8;


eq1(i,j)$cx(i,j)… Pij(i,j) =e= V2(i)*LN(j,i,‘Gij’) - Rij(i,j)*LN(j,i,‘Gij’) - Uij(i,j)*LN(j,i,‘Bij’);

eq2(i,j)$cx(i,j)… Qij(i,j) =e= -V2(i)*LN(j,i,‘Bij’) + Rij(i,j)*LN(j,i,‘Bij’) - Uij(i,j)*LN(j,i,‘Gij’);

eq3(i,j)$cx(i,j)… (V2(i)V2(j)) =g= sqr(Rij(i,j)) + sqr(Uij(i,j));
eq3(i,j)$cx(i,j)… sqr(2Rij(i,j)) + sqr(2Uij(i,j))+ sqr(V2(i) - V2(j)) =l= sqr(V2(i) + V2(j));
eq4(i,j)$cx(i,j)… Uij(j,i) =e= -Uij(i,j);
eq5(j,i)$cx(j,i)… Rij(j,i) =e= Rij(i,j);

eq6(i)… Pss(i)- BD(i,‘Pd’)/Sbase =e= sum(j$cx(j,i), Pij(i,j));

eq7(i)… Qss(i)- BD(i,‘Qd’)/Sbase =e= sum(j$cx(j,i), Qij(i,j));

eq8… OF =e= sum(i,Pss(i));

V2.lo(i) = 0.81;
V2.up(i) = 1.21;
V2.fx(slack) = 1;

Rij.lo(i,j)$cx(i,j) = 0;
Rij.up(i,j)$cx(i,j) = 1.21;
Uij.lo(i,j)$cx(i,j) = -1.21;
Uij.up(i,j)$cx(i,j) = 1.21;

Pij.up(i,j)((cx(i,j))) = 5; Pij.lo(i,j)((cx(i,j))) =-5;
Qij.up(i,j)((cx(i,j))) = 5; Qij.lo(i,j)((cx(i,j))) =-5;

option qcp = mosek;
Pss.lo(slack(i)) = 0;
Qss.lo(slack(i)) = 0;
Pss.up(slack(i)) = 5;
Qss.up(slack(i)) = 5;
Pss.fx(non_slack(i)) = 0;
Qss.fx(non_slack(i)) = 0;

Model loadflow / all /;

solve loadflow using qcp minimizing OF;

execute_unload ‘results.gdx’,
OF.l, Pss.l, Qss.l, V2.l;

execute ‘gdx2xls results.gdx results3.xlsx’;

5

Dear sir,

I am attaching my code. Please, help me.
loadflow_33bus.gms (4.51 KB)
.

6

As I said before, MOSEK is strict for what it expects as SOCP constraints. The definition of the various SOCPs can be found e.g. here: https://www.gams.com/latest/docs/UG_LanguageFeatures.html#UG_LanguageFeatures_ConicProgramming_Intro. So you have a rotated quadratic cone for MOSEK you can’t have constant (e.g. 2Uij(i,j)) and the you need the explicit 2x*y with x and y positive variables, not terms. Hence you need to make some new variables (MOSEK’s presolve will take care of substituting these out again):

variable vvm(i,j), rij2(i,j), uij2(i,j);
positive variables vvp(i,j),vvp2(i,j);
equation defvvm(i,j), defvvp(i,j), defvvp2(i,j), defrij2(i,j), defuij2(i,j);
defvvm(i,j)$cx(i,j).. vvm(i,j) =e= V2(i) - V2(j); [attachment=0]loadflow_33bus.gms[/attachment]
defvvp(i,j)$cx(i,j).. vvp(i,j) =e= V2(i) + V2(j);
defvvp2(i,j)$cx(i,j).. vvp2(i,j) =e= vvp(i,j)/2;
defuij2(i,j)$cx(i,j).. uij2(i,j) =e= 2*Uij(i,j);
defrij2(i,j)$cx(i,j).. rij2(i,j) =e= 2*Rij(i,j);
eq3(i,j)$cx(i,j).. sqr(rij2(i,j)) + sqr(uij2(i,j))+ sqr(vvm(i,j)) =l= 2*vvp(i,j)*vvp2(i,j);
*eq3(i,j)$cx(i,j).. sqr(2*Rij(i,j)) + sqr(2*Uij(i,j)) + sqr(V2(i) - V2(j)) =l= sqr(V2(i) + V2(j));

With that MOSEK solves the model just fine. The entire model is attached.

Good luck.
-Michael
loadflow_33bus.gms (5.01 KB)

7

Thank you bussieck. It’s working now.

8

Dear Michael,

I have reformulated my optimization problem. I try to solve with the mosek solver. It is showing like below

Variable ‘u(1)’ (128) is a member of cone ‘’ (0). (1307)
*** mosekgms.c:3263: MSK_GAMS_readtask: Error code MSK_RES_ERR_CONE_OVERLAP_APPEND (1307): The cone to be appended has one variable which is already member of another cone.

I attached my code. please help me.
loadflow_mosek.gms (4.54 KB)
Thanks in advance.

Sincerly,
Kasi