globalex.m

optimization toolbox

% Nothing prevents us from adding an additional 
% constraint derived from trace(P)=1
F = F + set(trace(A'*P+P*A) < -2*t);
pause

solvesdp(F,-t,ops);
pause

% When we add redundant constraints, we improve the relaxations
% and obtain better lower bounds. For the upper bounds however,
% these cuts only add additional computational burden.
% To overcome this, the command CUT can be used.
% Constraints defined with CUT are not used when the upper 
% bound problems are solved, but are only used in relaxations
F = set(P>0)+set(A'*P+P*A < -2*t*P);
F = F + set(trace(P)==1) + set(100 > t > 0) ;
F = F + cut(trace(A'*P+P*A) < -2*t);
pause

solvesdp(F,-t,ops);
pause

% Finally, it should be mentioned that the branch & bound
% algorithm can run without any installed local BMI solver. 
% This version of the algorithms is obtained by specifying 
% the upper bound solver as 'none'
ops = sdpsettings(ops,'bmibnb.uppersolver','none');
F = set(P>0)+set(A'*P+P*A < -2*t*P);
F = F + set(trace(P)==1) + set(100 > t > 0) ;
F = F + cut(trace(A'*P+P*A) < -2*t);
pause

solvesdp(F,-t,ops);
pause

echo off



break

% ********************************
% Sahinidis (-6.666 in 1)
% ********************************
yalmip('clear');
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
F = set(0 < x1 < 6)+set(0 < x2 < 4) + set(x1*x2 < 4) 
solvesdp(F,-x1-x2,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bmibnb.maxiter',100,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','cdd','cdd.method','dual-simplex'))

% ****************************
% Decay-rate (-2.5)
% ****************************
yalmip('clear');
A = [-1 2;-3 -4];
t = sdpvar(1,1);
P = sdpvar(2,2);
F = set(P>eye(2))+set(A'*P+P*A < -2*t*P);
F = F + set(-1e4 < P(:) < 1e4) + set(100 > t > -100) ;
ops = sdpsettings('bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bmibnb.maxiter',100,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','cdd','cdd.method','dual-simplex');
solvesdp(F,-t,ops);

% ********************************
% Classical example from Safonov etc
% ********************************
yalmip('clear')
x = sdpvar(1,1);
y = sdpvar(1,1);
t = sdpvar(1,1);
A0 = [-10 -0.5 -2;-0.5 4.5 0;-2 0 0];
A1 = [9 0.5 0;0.5 0 -3 ; 0 -3 -1];
A2 = [-1.8 -0.1 -0.4;-0.1 1.2 -1;-0.4 -1 0];
K12 = [0 0 2;0 -5.5 3;2 3 0];
F = lmi(x>-0.5)+lmi(x<2) + lmi(y>-3)+lmi(y<7) + lmi(t>-1000)+lmi(t<1000);%set(x+y<2.47)+set(x+y>2.47);
F = F + lmi(A0+x*A1+y*A2+x*y*K12-t*eye(3)<0);
solvesdp(F,t,sdpsettings('bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% Lassere 1
% optimal : -4, takes 13
% ********************************
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x = [x1;x2;x3];
B = [0 0 1;0 -1 0;-2 1 -1];
b = [3;0;-4];
v = [0;-1;-6];
r = [1.5;-0.5;-5];
p = -2*x1+x2-x3;
F = set(x1+x2+x3 < 4) + set(x1<2)+set(x3<3) + set(3*x2+x3 < 6);
F = F + set(x>0) + set(x'*B'*B*x-2*r'*B*x+r'*r-0.25*(b-v)'*(b-v) >0);
solvesdp(F,p,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% Lassere 2 
% optimal : -310 in 3
% ********************************
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
x6 = sdpvar(1,1);
t = sdpvar(1,1);
p = -25*(x1-2)^2-(x2-2)^2-(x3-1)^2-(x4-4)^2-(x5-1)^2-(x6-4)^2;
F = set((x3-3)^2+x4>4)+set((x5-3)^2+x6>4);
F = F + set(x1-3*x2<2)+set(-x1+x2<2);
F = F + set(x1-3*x2 <2)+set(x1+x2>2);
F = F + set(6>x1+x2>2);
F = F + set(1<x3<5) + set(0<x4<6)+set(1<x5<5)+set(0<x6<10)+set(x1>0)+set(x2>0);

solvesdp(F,p,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))
F = F + set(p<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))
% ********************************
% Lassere 2.6
% optimal : -39 in 6
% ********************************
clear all
Q = 100*eye(10);
c = [48, 42, 48, 45, 44, 41, 47, 42, 45, 46]';
b = [-4, 22,-6,-23,-12]';
A =[-2 -6 -1 0 -3 -3 -2 -6 -2 -2;
6 -5 8 -3 0 1 3 8 9 -3;
-5 6 5 3 8 -8 9 2 0 -9;
9 5 0 -9 1 -8 3 -9 -9 -3;
-8 7 -4 -5 -9 1 -7 -1 3 -2];
x = sdpvar(10,1);
t = sdpvar(1,1);
p = c'*x-0.5*x'*Q*x;
F = set(0<x<1)+set(A*x<b);
solvesdp(F,p,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
F = F + set(p<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
% ********************************
% Lassere
% -0.37? >100
% ********************************
clear all
x = sdpvar(10,1);
t = sdpvar(1,1);
x(1) = 1-sum(x(2:end));
obj = 0;
for i = 1:9
    obj = obj+x(i)*x(i+1);
end
for i = 1:8
    obj = obj+x(i)*x(i+2);
end
obj = obj + x(1)*x(7)+ x(1)*x(9)+ x(2)*x(10)+ x(4)*x(7);
F = set(1>x>0);
F = F + set(-obj<t);
solvesdp(F,-obj,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
F = F + set(-obj<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% Is J co-positive ?
% Optimal : 0
% ********************************
clear all
t = sdpvar(1,1);
x = sdpvar(5,1);
J = [1 -1  1  1 -1;
    -1  1 -1  1  1;
    1 -1  1 -1  1;
    1  1 -1  1 -1;
    -1  1  1 -1  1];
F = set(x>0) + set(sum(x)<1);
F = F+set(-10<t<10);
solvesdp(F,x'*J*x,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
F = F + set(x'*J*x<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% BARON example problems
% -17 in 9
% ********************************   
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
t = sdpvar(1,1);
p = 42*x1-50*x1^2+44*x2-50*x2^2+45*x3-50*x3^2+47*x4-50*x4^2+47.5*x5-50*x5^2;
F = set([20 12 11 7 4]*[x1;x2;x3;x4;x5] < 40) + set(0<[x1;x2;x3;x4;x5]<1);
solvesdp(F,p,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','nag','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
F = F+set(p<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','nag','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% BARON example
% -13 in 1
% ********************************   
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
t = sdpvar(1,1);
p = x1 - x2 - x3 - x1*x3 + x1*x4 + x2*x3 - x2*x4;
F = set(x1+4*x2 <= 8) + set(4*x1+x2 <= 12) + set(3*x1+4*x2 <= 12) + set(2*x3+x4 <= 8) + set(x3+2*x4 <= 8) + set(x3+x4 <= 5)+set([x1;x2;x3;x4]>0);
F = F+set(p<t);
solvesdp(F,p,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))
F = F+set(p<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.09e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% Control design
% 5.46 in 12
% ********************************   
A = [1 2;-3 0];B = [1;1];
[K0,P0] = lqr(A,B,eye(2),1);
P = sdpvar(2,2);setsdpvar(P,2*P0);K0(K0>1)=1;K0(K0<-1)=-1;
K = sdpvar(1,2);setsdpvar(K,-K0);
F = set(K<1)+set(K>-1)+set(P>0)+set((A+B*K)'*P+P*(A+B*K) < -eye(2)-K'*K);
%F = F + set([-(A+B*K)'*P-P*(A+B*K)-eye(2) K';K eye(1)] > 0)
F = F+lmi(diag(P)>0)+lmi(P(:)>-151) + lmi(P(:)<150) + lmi(P>P0)+lmi(K>-100) + lmi(K<100);
solvesdp(F,trace(P),sdpsettings('usex0',1,'bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk','bmibnb.lpreduce',1))

% ********************************
% GAMS st_qpc_m1
% optimal -473.77 in 3
% ********************************
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
t = sdpvar(1,1);

F = set([x1;x2;x3;x4;x5]>0);
F = F + set(x1 + x2 + 2*x3 + x4 + x5 > 10);
F = F + set(2*x1 + 3*x2 + x5 > 8);
F = F + set(x2 + 4*x3 - x4 + 2*x5 > 12);
F = F + set(8*x1 - x2 - x3 + 6*x4 > 20);
F = F + set(- 2*x1 - x2 - 3*x3 - x4 - x5 > -30);
obj = -(- (10*x1 - 0.34*x1*x1 - 0.28*x1*x2 + 10*x2 - 0.22*x1*x3 + 10*x3 - 0.24*x1*x4 + 10*x4 - 0.51*x1*x5 + 10*x5 - 0.28*x2*x1 - 0.34*x2*x2 - 0.23*x2*x3 -0.24*x2*x4 - 0.45*x2*x5 - 0.22*x3*x1 - 0.23*x3*x2 - 0.35*x3*x3 - 0.22*x3*x4 - 0.34*x3*x5 - 0.24*x4*x1 - 0.24*x4*x2 - 0.22*x4*x3 - 0.2*x4*x4 - 0.38*x4*x5 - 0.51*x5*x1 - 0.45*x5*x2 - 0.34*x5*x3 - 0.38*x5*x4 - 0.99*x5*x5));

solvesdp(F,obj,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))
F = F + set(obj<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))

% ********************************
% GAMS st_qpk1
% optimal -3 in 7
% ********************************
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
t = sdpvar(1,1);

F = set([x1;x2]>0);
F = F + set(- x1 + x2 < 1);
F = F + set(x1 - x2 < 1);
F = F + set(- x1 + 2*x2 < 3);
F = F + set(2*x1 - x2 < 3);
obj = (2*x1 - 2*x1*x1 + 2*x1*x2 + 3*x2 - 2*x2*x2);
F = F + set(obj<t);
solvesdp(F,t,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))


% ********************************
% GAMS st_robot
% optimal -3.52 in 4
% ********************************
yalmip('clear')
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
x6 = sdpvar(1,1);
x7 = sdpvar(1,1);
x8 = sdpvar(1,1);
x = [x1;x2;x3;x4;x5;x6;x7;x8];
F = set(-1<x<1);
F = F + set( 0.004731*x1*x3 - 0.1238*x1 - 0.3578*x2*x3 - 0.001637*x2 - 0.9338*x4 + x7 == 0.3571);
F = F + set(0.2238*x1*x3 + 0.2638*x1 + 0.7623*x2*x3 - 0.07745*x2 - 0.6734*x4 - x7 == 0.6022);
F = F + set( x6*x8 + 0.3578*x1 + 0.004731*x2 == 0);
F = F + set( - 0.7623*x1 + 0.2238*x2 == -0.3461);
F = F + set(x1^2 + x2^2 == 1);
F = F + set(x3^2 + x4^2 == 1);
F = F + set(x5^2 + x6^2 == 1);
F = F + set(x7^2 + x8^2 == 1);

solvesdp(F,sum(x),sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','pensdp','verbose',2,'solver','bmibnb','bnb.maxiter',30,'penbmi.P0',0.9e-1,'bmibnb.lpsolver','glpk'))


% ********************************
% GAMS st_robot
% optimal -382 in 5
% ********************************
yalmip('clear')
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
x6 = sdpvar(1,1);
x7 = sdpvar(1,1);
x8 = sdpvar(1,1);
x9 = sdpvar(1,1);
x10 = sdpvar(1,1);
t = sdpvar(1,1);
x = [x1;x2;x3;x4;x5;x6;x7;x8;x9;x10];

e1=    20*x1 + 20*x2 + 60*x3 + 60*x4 + 60*x5 + 60*x6 + 5*x7 + 45*x8 + 55*x9   + 65*x10 - 600.1;

e2=    5*x1 + 7*x2 + 3*x3 + 8*x4 + 13*x5 + 13*x6 + 2*x7 + 14*x8 + 14*x9      + 14*x10 - 310.5;

e3=    100*x1 + 130*x2 + 50*x3 + 70*x4 + 70*x5 + 70*x6 + 20*x7 + 80*x8 + 80*x9      + 80*x10 - 1800;

e4=    200*x1 + 280*x2 + 100*x3 + 200*x4 + 250*x5 + 280*x6 + 100*x7 + 180*x8      + 200*x9 + 220*x10 - 3850;

e5=    2*x1 + 2*x2 + 4*x3 + 4*x4 + 4*x5 + 4*x6 + 2*x7 + 6*x8 + 6*x9 + 6*x10      - 18.6;

e6=    4*x1 + 8*x2 + 2*x3 + 6*x4 + 10*x5 + 10*x6 + 5*x7 + 10*x8 + 10*x9      + 10*x10 - 198.7;

e7=    60*x1 + 110*x2 + 20*x3 + 40*x4 + 60*x5 + 70*x6 + 10*x7 + 40*x8 + 50*x9      + 50*x10 - 882;

e8=    150*x1 + 210*x2 + 40*x3 + 70*x4 + 90*x5 + 105*x6 + 60*x7 + 100*x8      + 140*x9 + 180*x10 - 4200;

e9=    80*x1 + 100*x2 + 6*x3 + 16*x4 + 20*x5 + 22*x6 + 20*x8 + 30*x9 + 30*x10      - 40.25;

e10=    40*x1 + 40*x2 + 12*x3 + 20*x4 + 24*x5 + 28*x6 + 40*x9 + 50*x10       - 327;

obj =   (10*x1 - 6.8*x1*x1 - 4.6*x1*x2 + 10*x2 - 7.9*x1*x3 + 10*x3 - 5.1*x1*x4       + 10*x4 - 6.9*x1*x5 + 10*x5 - 6.8*x1*x6 + 10*x6 - 4.6*x1*x7 + 10*x7 -      7.9*x1*x8 + 10*x8 - 5.1*x1*x9 + 10*x9 - 6.9*x1*x10 + 10*x10 - 4.6*x2*x1       - 5.5*x2*x2 - 5.8*x2*x3 - 4.5*x2*x4 - 6*x2*x5 - 4.6*x2*x6 - 5.5*x2*x7 -      5.8*x2*x8 - 4.5*x2*x9 - 6*x2*x10 - 7.9*x3*x1 - 5.8*x3*x2 - 13.3*x3*x3 -      6.7*x3*x4 - 8.9*x3*x5 - 7.9*x3*x6 - 5.8*x3*x7 - 13.3*x3*x8 - 6.7*x3*x9 -      8.9*x3*x10 - 5.1*x4*x1 - 4.5*x4*x2 - 6.7*x4*x3 - 6.9*x4*x4 - 5.8*x4*x5 -      5.1*x4*x6 - 4.5*x4*x7 - 6.7*x4*x8 - 6.9*x4*x9 - 5.8*x4*x10 - 6.9*x5*x1 -      6*x5*x2 - 8.9*x5*x3 - 5.8*x5*x4 - 11.9*x5*x5 - 6.9*x5*x6 - 6*x5*x7 - 8.9*      x5*x8 - 5.8*x5*x9 - 11.9*x5*x10 - 6.8*x6*x1 - 4.6*x6*x2 - 7.9*x6*x3 - 5.1      *x6*x4 - 6.9*x6*x5 - 6.8*x6*x6 - 4.6*x6*x7 - 7.9*x6*x8 - 5.1*x6*x9 - 6.9*      x6*x10 - 4.6*x7*x1 - 5.5*x7*x2 - 5.8*x7*x3 - 4.5*x7*x4 - 6*x7*x5 - 4.6*x7      *x6 - 5.5*x7*x7 - 5.8*x7*x8 - 4.5*x7*x9 - 6*x7*x10 - 7.9*x8*x1 - 5.8*x8*      x2 - 13.3*x8*x3 - 6.7*x8*x4 - 8.9*x8*x5 - 7.9*x8*x6 - 5.8*x8*x7 - 13.3*x8      *x8 - 6.7*x8*x9 - 8.9*x8*x10 - 5.1*x9*x1 - 4.5*x9*x2 - 6.7*x9*x3 - 6.9*x9      *x4 - 5.8*x9*x5 - 5.1*x9*x6 - 4.5*x9*x7 - 6.7*x9*x8 - 6.9*x9*x9 - 5.8*x9*      x10 - 6.9*x10*x1 - 6*x10*x2 - 8.9*x10*x3 - 5.8*x10*x4 - 11.9*x10*x5 - 6.9      *x10*x6 - 6*x10*x7 - 8.9*x10*x8 - 5.8*x10*x9 - 11.9*x10*x10);

F = set(0<x<10000);
F = F + set([e1;e2;e3;e4;e5;e6;e7;e8;e9;e10]<0);
%F = F + set(obj==t);
solvesdp(F,obj,sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk','solver','bmibnb','bmibnb.maxiter',15,'penbmi.P0',0.01,'bmibnb.lpsolver','glpk'))