globalex.m

optimization toolbox

clc
echo on
%*********************************************************
%
% Global bilinear programming
%
%*********************************************************
%
% YALMIP comes with a simple branch-and-bound code
% for solving problems with bilinear/quadratic constraints.
%
% The code is under development, and is currently only 
% applicable to small problems. 
%
% For the examples below to work, you must have an efficient
% LP solver (cplex,clp,nag,...) insrtalled and a solver for 
% general polynomial problems (currently only fmincon or PENBMI)


pause
yalmip('clear')
clc
%*********************************************************
% The first example is the problem from the moment 
% relaxation demo (concave quadratic constraint).
%*********************************************************
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);

objective = -2*x1+x2-x3;

F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0);
F = F + set(4-(x1+x2+x3)>0);
F = F + set(6-(3*x2+x3)>0);
F = F + set(x1>0);
F = F + set(2-x1>0);
F = F + set(x2>0);
F = F + set(x3>0);
F = F + set(3-x3>0);
pause

% Explicitly specify the global solver (this is currently recommended)
ops = sdpsettings('solver','bmibnb');                 % Global solver
pause

% Solve!
solvesdp(F,objective,ops)
pause
yalmip('clear')
clc
%***********************************************************
% The second example is slightly larger. It is a problem 
% with a concave quadratic objective function, and two 
% concave quadratic constraints.
%***********************************************************
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
x4 = sdpvar(1,1);
x5 = sdpvar(1,1);
x6 = sdpvar(1,1);

objective = -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) + 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);

pause

% The global solver in YALMIP can currently only deal with linear objectives.
% A simple epi-graph formulation solves this.
t = sdpvar(1,1);
F = F + set(objective<t);
pause

% Solve!
solvesdp(F,t,ops);
pause
clc
%***********************************************************
% Random bound-constrained indefinte quadratic program.
%***********************************************************

n = 5;
x = sdpvar(n,1);
t = sdpvar(1,1);

A = randn(n,n);A = A*A';
B = randn(n,n);B = B*B';

objective = x'*A*x-x'*B*x;

F = set(-1 < x < 1) + set(objective<t);
pause

% Solve!
solvesdp(F,t,ops);
pause
yalmip('clear')
clc
%***********************************************************
% Random boolean quadratic programming.
%
% By adding a nonlinear constraint x(x-1)==0, we can solve
% boolean programming using global nonlinear programming.
%***********************************************************

n = 10;
x = sdpvar(n,1);
t = sdpvar(1,1);

Q = randn(n,n);Q = Q*Q'/norm(Q)^2;
c = randn(n,1);

objective = x'*Q*x+c'*x;

F = set(x.*(x-1)==0) + set(objective < t);
pause

% Note, all complicating variables (variables involved in nonlinear terms)
% must be bounded (either explicitely or implicetely via linear and relaxed 
% nonlinear constraints).
%
% The easible set for the relaxed problem is not bounded
% so YALMIP will not be able to derive variable bounds.
% We add some valid bounds to fix this
F = F + set(0<x<1);
pause

% Solve!
solvesdp(F,t,ops);
double(x'*Q*x+c'*x)
pause

% Let us compare with an integer solver (use we use the internal MIQP solver)
solvesdp(set(binary(x)),x'*Q*x+c'*x,sdpsettings(ops,'solver','bnb'));
double(x'*Q*x+c'*x)
pause

yalmip('clear')
clc
%***********************************************************
% This example is a somewhat larger indefinite quadratic
% programming problem, taken from the GAMS problem library.
%***********************************************************
pause

% The problem has 11 variables
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);
x11 = sdpvar(1,1);
t = sdpvar(1,1);
x = [x1;x2;x3;x4;x5;x6;x7;x8;x9;x10;x11];

pause

% ...and 22 linear constraints
e1 =  - x1 - 2*x2 - 3*x3 - 4*x4 - 5*x5 - 6*x6 - 7*x7 - 8*x8 - 9*x9 - 10*x10      - 11*x11 	- 0;
e2 =   - 2*x1 - 3*x2 - 4*x3 - 5*x4 - 6*x5 - 7*x6 - 8*x7 - 9*x8 - 10*x9 - 11*x10      - x11 	- 0;
e3 =   - 3*x1 - 4*x2 - 5*x3 - 6*x4 - 7*x5 - 8*x6 - 9*x7 - 10*x8 - 11*x9 - x10      - 2*x11 	- 0;
e4 =   - 4*x1 - 5*x2 - 6*x3 - 7*x4 - 8*x5 - 9*x6 - 10*x7 - 11*x8 - x9 - 2*x10      - 3*x11 	- 0;
e5 =   - 5*x1 - 6*x2 - 7*x3 - 8*x4 - 9*x5 - 10*x6 - 11*x7 - x8 - 2*x9 - 3*x10      - 4*x11 	- 0;
e6 =   - 6*x1 - 7*x2 - 8*x3 - 9*x4 - 10*x5 - 11*x6 - x7 - 2*x8 - 3*x9 - 4*x10      - 5*x11 	- 0;
e7 =   - 7*x1 - 8*x2 - 9*x3 - 10*x4 - 11*x5 - x6 - 2*x7 - 3*x8 - 4*x9 - 5*x10      - 6*x11 	- 0;
e8 =   - 8*x1 - 9*x2 - 10*x3 - 11*x4 - x5 - 2*x6 - 3*x7 - 4*x8 - 5*x9 - 6*x10      - 7*x11 	- 0;
e9 =   - 9*x1 - 10*x2 - 11*x3 - x4 - 2*x5 - 3*x6 - 4*x7 - 5*x8 - 6*x9 - 7*x10      - 8*x11 	- 0;
e10 =   - 10*x1 - 11*x2 - x3 - 2*x4 - 3*x5 - 4*x6 - 5*x7 - 6*x8 - 7*x9 - 8*x10       - 9*x11 	- 0;
e11 =   - 11*x1 - x2 - 2*x3 - 3*x4 - 4*x5 - 5*x6 - 6*x7 - 7*x8 - 8*x9 - 9*x10       - 10*x11 	- 0;
e12 =     x1 + 2*x2 + 3*x3 + 4*x4 + 5*x5 + 6*x6 + 7*x7 + 8*x8 + 9*x9 + 10*x10       + 11*x11 	- 66;
e13 =     2*x1 + 3*x2 + 4*x3 + 5*x4 + 6*x5 + 7*x6 + 8*x7 + 9*x8 + 10*x9 + 11*x10       + x11 	- 66;
e14 =     3*x1 + 4*x2 + 5*x3 + 6*x4 + 7*x5 + 8*x6 + 9*x7 + 10*x8 + 11*x9 + x10       + 2*x11 	- 66;
e15 =     4*x1 + 5*x2 + 6*x3 + 7*x4 + 8*x5 + 9*x6 + 10*x7 + 11*x8 + x9 + 2*x10       + 3*x11 	- 66;
e16 =     5*x1 + 6*x2 + 7*x3 + 8*x4 + 9*x5 + 10*x6 + 11*x7 + x8 + 2*x9 + 3*x10       + 4*x11 	- 66;
e17 =     6*x1 + 7*x2 + 8*x3 + 9*x4 + 10*x5 + 11*x6 + x7 + 2*x8 + 3*x9 + 4*x10       + 5*x11 	- 66;
e18 =     7*x1 + 8*x2 + 9*x3 + 10*x4 + 11*x5 + x6 + 2*x7 + 3*x8 + 4*x9 + 5*x10       + 6*x11 	- 66;
e19 =     8*x1 + 9*x2 + 10*x3 + 11*x4 + x5 + 2*x6 + 3*x7 + 4*x8 + 5*x9 + 6*x10       + 7*x11 	- 66;
e20 =     9*x1 + 10*x2 + 11*x3 + x4 + 2*x5 + 3*x6 + 4*x7 + 5*x8 + 6*x9 + 7*x10       + 8*x11 	- 66;
e21 =     10*x1 + 11*x2 + x3 + 2*x4 + 3*x5 + 4*x6 + 5*x7 + 6*x8 + 7*x9 + 8*x10       + 9*x11 	- 66;
e22 =     11*x1 + x2 + 2*x3 + 3*x4 + 4*x5 + 5*x6 + 6*x7 + 7*x8 + 8*x9 + 9*x10       + 10*x11 	- 66;

pause


% Define the objective function and the whole program
obj =  (0.5*x1*x2 - x1*x1 + 0.5*x2*x1 - x2*x2 + 0.5*x2*x3 + 0.5*x3*x2 - x3*x3       + 0.5*x3*x4 + 0.5*x4*x3 - x4*x4 + 0.5*x4*x5 + 0.5*x5*x4 - x5*x5 + 0.5*x5      *x6 + 0.5*x6*x5 - x6*x6 + 0.5*x6*x7 + 0.5*x7*x6 - x7*x7 + 0.5*x7*x8 + 0.5      *x8*x7 - x8*x8 + 0.5*x8*x9 + 0.5*x9*x8 - x9*x9 + 0.5*x9*x10 + 0.5*x10*x9       - x10*x10 + 0.5*x10*x11 + 0.5*x11*x10 - x11*x11);
e = -[e1;e2;e3;e4;e5;e6;e7;e8;e9;e10;e11;e12;e13;e14;e15;e16;e17;e18;e19;e20;e21;e22];
F = set(10>x>0);
F = F + set(e>0);
F = F + set(obj==t);

% We will nowe try to solve this using the same strategy as before
solvesdp(F,t,ops)

pause
clc

% As an alternative, we might use additional cuts using the CUT functionality.
% These additional constraints are only used in domain reduction and
% solution of lower bound programs, but not in the local solver.

% We square all linear constraints and add these as cuts
%
% Note that the triu operator leaves returns a matrix with
% zeros below the diagonal. These are however neglected by YALMIP
% in the call to SET.
%
f = sdpvar(F(find(is(F,'linear'))));
F = F + cut(triu(f*f') > 0);
pause

% ...and solve (To speed up things, we turn off LP based domain reduction)
%
% NOTE : The first ~15 iterations are rather slow, but things starts to 
%        speed up when YALMIP starts pruning redundnant linear constraints.
%
solvesdp(F,t,sdpsettings(ops,'bmibnb.lpreduce',0))
pause
clc

% The global solver in YALMIP can only solve bilinear programs, 
% but a pre-processor is capable of transforming higher order problems 
% to bilinear programs. As an example, the variable x^3y^2 will be replaced 
% with the the variable w and the constraints w == uv, u == zx, z == x^2, v == y^2. 
% The is done automatically if the global solver, or PENBMI, is called with 
% a higher order polynomial problem. Note that this conversion is rather 
% inefficient, and only very small problems can be addressed using this simple approach.
% Additionally, this module has not been tested in any serious sense.
pause

% Let us solve a small trivial polynomial program
sdpvar x y
F = set(x^3+y^5 < 5) + set(y > 0);
F = F + set(-100 < [x y] < 100) % Always bounded domain
pause

solvesdp(F,-x,ops)
pause

%***********************************************************
% The main motivation for implementing a bilinear solver in
% YALMIP is matrix-valued problems, i.e. BMI problems.
%
% Of-course, BMI problems are even harder than the nonconvex 
% quadratic problems solved above. However, in some cases, 
% it is actually possible to solve BMI problems globally.
%
% Don't expect too much though...
%***********************************************************

pause
clc

%***********************************************************
% The following is a classical problem
%***********************************************************
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];

x = sdpvar(1,1);
y = sdpvar(1,1);
t = sdpvar(1,1);

F = set(A0+x*A1+y*A2+x*y*K12-t*eye(3)<0);
F = F + set(-0.5 < x < 2) + set(-3 < y < 7);
pause

% Specify solvers (this requires PENBMI!)
ops = sdpsettings('solver','bmibnb');                 % Global solver
ops = sdpsettings(ops,'bmibnb.uppersolver','penbmi'); % Local solver

pause

% Solve!
solvesdp(F,t,ops);
pause
clc
%***********************************************************
% The next example shows how to solve a standard LQ control
% problem with constraints on the feedback gain
%***********************************************************

A = [-1 2;-3 -4];B = [1;1];

P = sdpvar(2,2);
K = sdpvar(1,2);

F = set(P>0)+set((A+B*K)'*P+P*(A+B*K) < -eye(2)-K'*K);
F = F + set(K<0.1)+set(K>-0.1);
F = F + set(1000> P(:)>-1000);
solvesdp(F,trace(P),ops);
pause
clc
%***********************************************************
% Decay-rate estimation as a global BMI example ...
%
% Note, this is a quasi-convex problem that PENBMI actually 
% solves to global optimality, so the use of a global solver
% is pretty stupid. It is only included here to show some 
% properties of the global solver.
%
% may run for up to 100 iterations depending on the solvers
%***********************************************************

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 > 0) ;
pause
solvesdp(F,-t,ops);
pause

% Now, that sucked, right! 100 iterations and still 10% gap...
%
% The way we model a problem can make a huge difference.
% The decay-rate problem can be substantially improved by 
% using a smarter model. The original problem is homogenous 
% w.r.t P, and a better way to make the feasible set bounded
% is to constrain the trace of P
A = [-1 2;-3 -4];
t = sdpvar(1,1);
P = sdpvar(2,2);
F = set(P>0)+set(A'*P+P*A < -2*t*P);
F = F + set(trace(P)==1) + set(100 > t > 0) ;
pause
solvesdp(F,-t,ops);
pause