Global solution of polynomial programs
 This
example requires the solver PENBMI as a
local nonlinear solver (fmincon
is sufficient for the first examples).
Moreover, the examples below are all coded to use PENSDP
and GLPK for solving the relaxations.
Global solutions! Well, almost... don't expect too much at this stage.
The functionality described here is under development. The code
is fairly robust on small bilinear and indefinite quadratic
optimization problems, and a couple of small problems
with bilinear matrix inequalities have been solved successfully.
The algorithm is based on a simple spatial branch and bound strategy,
using McCormick's convex envelopes for bounding bilinear terms. In
addition to this, LP-based bound tightening can be applied iteratively to
improve
variable bounds. Relaxed problems are solved using either an LP solver or an SDP
solver, depending on the problem, while upper bounds are found using the
local solver PENBMI. A more detailed
description of the algorithm will be presented in the future.
Options
In the examples below, we are mainly using the default settings when solving
the problem, but there are several options that can be changed to
alter the computations. The most important are
sdpsettings('bmibnb.roottight',[0|1]) |
Improve variable bounds at root-node by performing bound
strengthening based on the full relaxed model (Can be very expensive,
but lead to improved branching) |
sdpsettings('bmibnb.lpreduce',[0|1]) |
Improve variable bounds in all nodes by performing bound
strengthening using only the scalar constraints (including scalar cut
constraints) in the model (Can be very expensive, but lead to
improved branching, in particular for semidefinite problems) |
sdpsettings('bmibnb.lowersolver', solvertag) |
Solver for relaxed problems. |
sdpsettings('bmibnb.uppersolver', solvertag) |
Local solver for computing upper bounds. |
sdpsettings('bmibnb.lpsolver', solvertag) |
LP solver for bound strengthening
(only used if bmibnb.lpreduce is set to 1) |
sdpsettings('bmibnb.numglobal', [int]) |
A major computational burden along the branching process is to
solver the upper bound problems using a local nonlinear solver. By
setting this value to a finite value, the local solver will no
longer be used when the upper bound has been improved numglobal times.
This option is useful if you believe your local solver quickly gives
globally optimal solutions. It can also be useful if you have a
feasible solution and just want to compute a gap. In this case, use
the usex0 option and set numglobal to 0. |
Nonconvex quadratic programming
The first example is a problem with a concave quadratic constraint (this
is the example addressed in the moment
relaxation section). Three different optimization problems are solved
during the branching: Upper bounds using a local nonlinear solver (bmibnb.uppersolver),
lower bounds (bmibnb.lowersolver) and bound tightening using
linear programming (bmibnb.lpsolver).
clear all
x1 = sdpvar(1,1);
x2 = sdpvar(1,1);
x3 = sdpvar(1,1);
p = -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);
options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2);
solvesdp(F,p,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : glpk
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : Inf NaN -6.000E+000 2
2 : Inf NaN -6.000E+000 3
3 : -4.000E+000 40.00 -6.000E+000 4 Improved solution
4 : -4.000E+000 40.00 -6.000E+000 3 Infeasible
5 : -4.000E+000 37.98 -5.899E+000 4 Poor bound
6 : -4.000E+000 37.98 -5.899E+000 5
7 : -4.000E+000 25.80 -5.290E+000 6
8 : -4.000E+000 25.80 -5.290E+000 5 Infeasible
9 : -4.000E+000 7.08 -4.354E+000 4 Infeasible
10 : -4.000E+000 7.08 -4.354E+000 3 Infeasible
11 : -4.000E+000 0.06 -4.003E+000 4
+ 11 Finishing. Cost: -4 Gap: 0.06486%
|
The second example is a slightly larger problem indefinite quadratic
programming problem. The problem is easily solved to a gap of
less than 1%.
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);
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);
options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2,'solver','bmibnb');
solvesdp(F,p,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : glpk
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : Inf NaN -3.130E+002 2
2 : -3.100E+002 0.96 -3.130E+002 3 Improved solution
+ 2 Finishing. Cost: -310 Gap: 0.96463%
|
Quadratic equality constraints can also be dealt with. This can be used
for, e.g., boolean programming.
n = 10;
x = sdpvar(n,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(0<x<1);
options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2,'solver','bmibnb');
solvesdp(F,objective,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : glpk
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -2.967E+000 0.00 -2.967E+000 2 Improved solution
+ 1 Finishing. Cost: -2.9671 Gap: 2.5207e-009%
|
Nonconvex polynomial programming
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 x3y2
will be replaced with the the variable w and the constraints
w == uv, u == zx, z == x2, v == y2.
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.
sdpvar x y
F = set(x^3+y^5 < 5) + set(y > 0);
F = F + set(-100 < [x y] < 100) % Always bounded domain
solvesdp(F,-x,options)
* Starting YALMIP bilinear branch & bound.
* Lower solver : glpk
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : Inf NaN -2.990E+001 2
2 : Inf NaN -2.990E+001 1 Infeasible
3 : Inf NaN -2.763E+001 2
4 : Inf NaN -2.763E+001 3
5 : Inf NaN -1.452E+001 4
6 : Inf NaN -1.452E+001 5
7 : -1.710E+000 171.54 -6.359E+000 4 Improved solution
8 : -1.710E+000 171.54 -6.359E+000 3 Infeasible
9 : -1.710E+000 127.11 -5.155E+000 4 Poor bound
10 : -1.710E+000 127.11 -5.155E+000 3 Infeasible
11 : -1.710E+000 0.00 -1.710E+000 2 Infeasible
+ 11 Finishing. Cost: -1.71 Gap: 1.247e-005%
|
Nonconvex semidefinite programming
The following problem is a classical BMI problem
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);
F = F + lmi(A0+x*A1+y*A2+x*y*K12-t*eye(3)<0);
options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2,'solver','bmibnb');
solvesdp(F,t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -9.565E-001 2.22 -1.000E+000 2 Improved solution
2 : -9.565E-001 2.22 -1.000E+000 1 Infeasible
3 : -9.565E-001 2.22 -1.000E+000 2
4 : -9.565E-001 2.22 -1.000E+000 3
5 : -9.565E-001 0.24 -9.613E-001 2 Infeasible
+ 5 Finishing. Cost: -0.95653 Gap: 0.24126%
|
The constrained LQ problem in the section Local BMIs
can actually easily be solved to global optimality in just a couple of
iterations. For the global code to work, global lower and upper and bound on
all complicating variables (involved in nonlinear terms) must be supplied, either explicitly or implicitly in the linear
constraints. This was the case for all problems above. In this example, the variable
K is already bounded in
the original problem, but the elements in P have to be bounded
artificially.
yalmip('clear')
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)
options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2,'solver','bmibnb');
solvesdp(F,trace(P),options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : 4.834E-001 17.70 2.209E-001 2 Improved solution
2 : 4.834E-001 17.70 2.209E-001 1 Infeasible
3 : 4.834E-001 1.93 4.548E-001 2
4 : 4.834E-001 1.93 4.548E-001 1 Infeasible
5 : 4.834E-001 0.25 4.797E-001 2
+ 5 Finishing. Cost: 0.48341 Gap: 0.25032%
|
Beware, the BMI solver is absolutely not that efficient in general,
this was just a lucky case. Here is the decay-rate
example instead (with some additional constraints on the elements
in P to bound the feasible set).
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 > 0) ;
options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi');
options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');
options = sdpsettings(options,'verbose',2,'solver','bmibnb');
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : Inf NaN -9.965E+001 2
2 : -2.500E+000 2775.83 -9.965E+001 3 Improved solution
3 : -2.500E+000 2721.00 -9.773E+001 4
4 : -2.500E+000 2721.00 -9.773E+001 3 Infeasible
5 : -2.500E+000 2634.59 -9.471E+001 4
6 : -2.500E+000 2634.59 -9.471E+001 3 Infeasible
7 : -2.500E+000 2577.61 -9.272E+001 4
8 : -2.500E+000 2577.61 -9.272E+001 3 Infeasible
9 : -2.500E+000 2501.03 -9.004E+001 4
10 : -2.500E+000 2501.03 -9.004E+001 3 Infeasible
...
94 : -2.500E+000 4.30 -2.650E+000 49
95 : -2.500E+000 0.47 -2.516E+000 50
+ 95 Finishing. Cost: -2.5 Gap: 0.46843%
|
The linear relaxations give very poor lower bounds on this problem,
leading to poor convergence of the branching process. However, for this
particular problem, the reason is easy to find. The original BMI is
homogeneous w.r.t P, and to guarantee a somewhat reasonable
solution, we artificially added the constraint P≥I.
A better model is obtained if we instead fix the trace of P. This will
make the feasible set w.r.t P much smaller, but the problems are
equivalent. Note also that we no longer need any artificial constraints
on the elements in P.
F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ;
F = F + set(trace(P)==1);
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -2.500E+000 1396.51 -5.138E+001 2 Improved solution
2 : -2.500E+000 1396.51 -5.138E+001 1 Infeasible
3 : -2.500E+000 1.41 -2.549E+000 2
4 : -2.500E+000 1.41 -2.549E+000 1 Infeasible
5 : -2.500E+000 0.17 -2.506E+000 2
+ 5 Finishing. Cost: -2.5 Gap: 0.1746%
|
For this problem, we can easily find a very efficient additional cutting
plane. The
decay-rate BMI together with the constant trace on P implies trace(ATP+PA)≤-2t.
Adding this cut reduce the number of iterations needed.
F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0);
F = F + set(trace(P)==1);
F = F + set(trace(A'*P+P*A)<-2*t);
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -2.500E+000 26.60 -3.431E+000 2 Improved solution
2 : -2.500E+000 26.60 -3.431E+000 3
3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible
+ 3 Finishing. Cost: -2.5 Gap: 0.5461%
|
A Schur complement on the decay-rate BMI gives us yet another cut
which improves the node relaxation even more.
F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ;
F = F + set(trace(P)==1);
F = F + set(trace(A'*P+P*A)<-2*t)
F = F + set([-A'*P-P*A P*t;P*t P*t/2]);
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -2.500E+000 17.84 -3.124E+000 2 Improved solution
2 : -2.500E+000 17.84 -3.124E+000 3
3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible
+ 3 Finishing. Cost: -2.5 Gap: 0.55275%
|
By adding valid cuts, the relaxations are possibly tighter, leading to
better lower bounds. A problem however is that we add additional
burden to the local solver used for the upper bounds. The additional
cuts are redundant for the local solver, and most likely detoriate the
performance. To avoid this, cuts can be exlicitely specified by using
the command cut.
Constraints defined using this command (instead of
set) will only be used
in the solution of relaxations, and omitted when the local solver is
called.
F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ;
F = F + set(trace(P)==1);
F = F + cut(trace(A'*P+P*A)<-2*t);
F = F + cut([-A'*P-P*A P*t;P*t P*t/2]);
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : penbmi
Node Upper Gap(%) Lower Open
1 : -2.500E+000 17.84 -3.124E+000 2 Improved solution
2 : -2.500E+000 17.84 -3.124E+000 3
3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible
+ 3 Finishing. Cost: -2.5 Gap: 0.55275%
|
Upper bounds were obtained above by solving the BMI locally using
PENBMI. If no local BMI solver is available, an alternative is to
check if the relaxed solution is a feasible solution. If so,
the upper bound can be updated. This scheme
can be obtained by specifying 'none' as the upper bound solver.
options = sdpsettings(options,'bmibnb.uppersolver','none');
F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ;
F = F + set(trace(P)==1);
F = F + cut(trace(A'*P+P*A)<-2*t);
F = F + cut([-A'*P-P*A P*t;P*t P*t/2]);
solvesdp(F,-t,options);
* Starting YALMIP bilinear branch & bound.
* Lower solver : pensdp
* Upper solver : none
Node Upper Gap(%) Lower Open
1 : Inf NaN -3.124E+000 2
2 : Inf NaN -3.124E+000 3
3 : -9.045E-001 104.47 -2.894E+000 4 Improved solution
4 : -9.045E-001 104.47 -2.894E+000 5
5 : -1.467E+000 52.43 -2.761E+000 4 Improved solution
6 : -1.467E+000 52.43 -2.761E+000 5
7 : -1.831E+000 29.87 -2.677E+000 4 Improved solution
8 : -1.831E+000 29.87 -2.677E+000 5
9 : -2.071E+000 17.77 -2.616E+000 4 Improved solution
10 : -2.071E+000 17.77 -2.616E+000 5
11 : -2.229E+000 10.46 -2.567E+000 4 Improved solution
12 : -2.229E+000 10.46 -2.567E+000 5
13 : -2.332E+000 5.70 -2.522E+000 4 Improved solution
14 : -2.332E+000 5.70 -2.522E+000 5
15 : -2.397E+000 3.25 -2.508E+000 4 Improved solution
16 : -2.397E+000 3.25 -2.508E+000 5
17 : -2.435E+000 2.01 -2.504E+000 4 Improved solution
18 : -2.435E+000 2.01 -2.504E+000 5
19 : -2.459E+000 1.25 -2.503E+000 4 Improved solution
20 : -2.459E+000 1.25 -2.503E+000 5
21 : -2.478E+000 0.70 -2.502E+000 4 Improved solution
+ 21 Finishing. Cost: -2.4776 Gap: 0.69631%
|
More examples can be found in
yalmipdemo.
 The global
branch and bound algorithm will hopefully be improved upon
significantly in future releases. This includes both algorithmic
improvements and code optimization. If you have any ideas, you are welcome to contribute. |