[37] | 1 | <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN"> |
---|
| 2 | <html> |
---|
| 3 | |
---|
| 4 | <head> |
---|
| 5 | <meta http-equiv="Content-Language" content="en-us"> |
---|
| 6 | <title>YALMIP Example : Solving BMIs globally</title> |
---|
| 7 | <meta http-equiv="Content-Type" content="text/html; charset=windows-1251"> |
---|
| 8 | <meta content="Microsoft FrontPage 6.0" name="GENERATOR"> |
---|
| 9 | <meta name="ProgId" content="FrontPage.Editor.Document"> |
---|
| 10 | <link href="yalmip.css" type="text/css" rel="stylesheet"> |
---|
| 11 | <base target="_self"> |
---|
| 12 | </head> |
---|
| 13 | |
---|
| 14 | <body leftMargin="0" topMargin="0"> |
---|
| 15 | |
---|
| 16 | <div align="left"> |
---|
| 17 | |
---|
| 18 | <table border="0" cellpadding="4" cellspacing="3" style="border-collapse: collapse" bordercolor="#000000" width="100%" align="left" height="100%"> |
---|
| 19 | <tr> |
---|
| 20 | <td width="100%" align="left" height="100%" valign="top"> |
---|
| 21 | <h2>Global solution of polynomial programs</h2> |
---|
| 22 | <hr noShade SIZE="1"> |
---|
| 23 | <p> |
---|
| 24 | <img border="0" src="exclamationmark.jpg" align="left" width="16" height="16">This |
---|
| 25 | example requires the solver <a href="solvers.htm#penbmi">PENBMI</a> as a |
---|
| 26 | local nonlinear solver (<a href="solvers.htm#fmincon">fmincon</a> |
---|
| 27 | is sufficient for the first examples). |
---|
| 28 | Moreover, the examples below are all coded to use <a href="solvers.htm#pensdp">PENSDP</a> |
---|
| 29 | and <a href="solvers.htm#glpk">GLPK</a> for solving the relaxations.</p> |
---|
| 30 | <p>Global solutions! Well, almost... don't expect too much at this stage. |
---|
| 31 | The functionality described here is under development. The code |
---|
| 32 | is fairly robust on small bilinear and indefinite quadratic |
---|
| 33 | optimization problems, and a couple of small problems |
---|
| 34 | with bilinear matrix inequalities have been solved successfully.<br> |
---|
| 35 | <br> |
---|
| 36 | The algorithm is based on a simple spatial branch and bound strategy, |
---|
| 37 | using McCormick's convex envelopes for bounding bilinear terms. In |
---|
| 38 | addition to this, LP-based bound tightening can be applied iteratively to |
---|
| 39 | improve |
---|
| 40 | variable bounds. Relaxed problems are solved using either an LP solver or an SDP |
---|
| 41 | solver, depending on the problem, while upper bounds are found using the |
---|
| 42 | local solver <a href="solvers.htm#penbmi">PENBMI</a>. A more detailed |
---|
| 43 | description of the algorithm will be presented in the future.</p> |
---|
| 44 | <h3>Options</h3> |
---|
| 45 | <p>In the examples below, we are mainly using the default settings when solving |
---|
| 46 | the problem, but there are several options that can be changed to |
---|
| 47 | alter the computations. The most important are</p> |
---|
| 48 | <table border="1" cellspacing="1" style="border-collapse: collapse" width="100%" bordercolor="#000000" bgcolor="#FFFFE0" id="table1"> |
---|
| 49 | <tr> |
---|
| 50 | <td width="367"><code>sdpsettings('bmibnb.roottight',[0|1])</code></td> |
---|
| 51 | <td>Improve variable bounds at root-node by performing bound |
---|
| 52 | strengthening based on the full relaxed model (Can be very expensive, |
---|
| 53 | but lead to improved branching)</td> |
---|
| 54 | </tr> |
---|
| 55 | <tr> |
---|
| 56 | <td width="367"><code>sdpsettings('bmibnb.lpreduce',[0|1])</code></td> |
---|
| 57 | <td>Improve variable bounds in all nodes by performing bound |
---|
| 58 | strengthening using only the scalar constraints (including scalar cut |
---|
| 59 | constraints) in the model (Can be very expensive, but lead to |
---|
| 60 | improved branching, in particular for semidefinite problems)</td> |
---|
| 61 | </tr> |
---|
| 62 | <tr> |
---|
| 63 | <td width="367"><code>sdpsettings('bmibnb.lowersolver', solvertag)</code></td> |
---|
| 64 | <td>Solver for relaxed problems.</td> |
---|
| 65 | </tr> |
---|
| 66 | <tr> |
---|
| 67 | <td width="367"><code>sdpsettings('bmibnb.uppersolver', solvertag)</code></td> |
---|
| 68 | <td>Local solver for computing upper bounds.</td> |
---|
| 69 | </tr> |
---|
| 70 | <tr> |
---|
| 71 | <td width="367"><code>sdpsettings('bmibnb.lpsolver', solvertag)</code></td> |
---|
| 72 | <td>LP solver for bound strengthening |
---|
| 73 | (only used if <code>bmibnb.lpreduce</code> is set to 1)</td> |
---|
| 74 | </tr> |
---|
| 75 | <tr> |
---|
| 76 | <td width="367"><code>sdpsettings('bmibnb.numglobal', [int])</code></td> |
---|
| 77 | <td>A major computational burden along the branching process is to |
---|
| 78 | solver the upper bound problems using a local nonlinear solver. By |
---|
| 79 | setting this value to a finite value, the local solver will no |
---|
| 80 | longer be used when the upper bound has been improved <code>numglobal</code> times. |
---|
| 81 | This option is useful if you believe your local solver quickly gives |
---|
| 82 | globally optimal solutions. It can also be useful if you have a |
---|
| 83 | feasible solution and just want to compute a gap. In this case, use |
---|
| 84 | the <code>usex0</code> option and set <code>numglobal</code> to 0.</td> |
---|
| 85 | </tr> |
---|
| 86 | </table> |
---|
| 87 | <h3>Nonconvex quadratic programming</h3> |
---|
| 88 | <p>The first example is a problem with a concave quadratic constraint (this |
---|
| 89 | is the example addressed in the <a href="moment.htm">moment |
---|
| 90 | relaxation section</a>). Three different optimization problems are solved |
---|
| 91 | during the branching: Upper bounds using a local nonlinear solver (<code>bmibnb.uppersolver</code>), |
---|
| 92 | lower bounds (<code>bmibnb.lowersolver</code>) and bound tightening using |
---|
| 93 | linear programming (<code>bmibnb.lpsolver</code>).</p> |
---|
| 94 | <table cellPadding="10" width="100%"> |
---|
| 95 | <tr> |
---|
| 96 | <td class="xmpcode"> |
---|
| 97 | <pre>clear all |
---|
| 98 | x1 = sdpvar(1,1); |
---|
| 99 | x2 = sdpvar(1,1); |
---|
| 100 | x3 = sdpvar(1,1); |
---|
| 101 | p = -2*x1+x2-x3; |
---|
| 102 | F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0); |
---|
| 103 | F = F + set(4-(x1+x2+x3)>0); |
---|
| 104 | F = F + set(6-(3*x2+x3)>0); |
---|
| 105 | F = F + set(x1>0); |
---|
| 106 | F = F + set(2-x1>0); |
---|
| 107 | F = F + set(x2>0); |
---|
| 108 | F = F + set(x3>0); |
---|
| 109 | F = F + set(3-x3>0); |
---|
| 110 | options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk'); |
---|
| 111 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 112 | options = sdpsettings(options,'verbose',2); |
---|
| 113 | solvesdp(F,p,options); |
---|
| 114 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 115 | * Lower solver : glpk |
---|
| 116 | * Upper solver : penbmi |
---|
| 117 | Node Upper Gap(%) Lower Open |
---|
| 118 | 1 : Inf NaN -6.000E+000 2 |
---|
| 119 | 2 : Inf NaN -6.000E+000 3 |
---|
| 120 | 3 : -4.000E+000 40.00 -6.000E+000 4 Improved solution |
---|
| 121 | 4 : -4.000E+000 40.00 -6.000E+000 3 Infeasible |
---|
| 122 | 5 : -4.000E+000 37.98 -5.899E+000 4 Poor bound |
---|
| 123 | 6 : -4.000E+000 37.98 -5.899E+000 5 |
---|
| 124 | 7 : -4.000E+000 25.80 -5.290E+000 6 |
---|
| 125 | 8 : -4.000E+000 25.80 -5.290E+000 5 Infeasible |
---|
| 126 | 9 : -4.000E+000 7.08 -4.354E+000 4 Infeasible |
---|
| 127 | 10 : -4.000E+000 7.08 -4.354E+000 3 Infeasible |
---|
| 128 | 11 : -4.000E+000 0.06 -4.003E+000 4 |
---|
| 129 | + 11 Finishing. Cost: -4 Gap: 0.06486%</font></pre> |
---|
| 130 | </td> |
---|
| 131 | </tr> |
---|
| 132 | </table> |
---|
| 133 | <p>The second example is a slightly larger problem indefinite quadratic |
---|
| 134 | programming problem. The problem is easily solved to a gap of |
---|
| 135 | less than 1%. </p> |
---|
| 136 | <table cellPadding="10" width="100%"> |
---|
| 137 | <tr> |
---|
| 138 | <td class="xmpcode"> |
---|
| 139 | <pre>clear all |
---|
| 140 | x1 = sdpvar(1,1); |
---|
| 141 | x2 = sdpvar(1,1); |
---|
| 142 | x3 = sdpvar(1,1); |
---|
| 143 | x4 = sdpvar(1,1); |
---|
| 144 | x5 = sdpvar(1,1); |
---|
| 145 | x6 = sdpvar(1,1); |
---|
| 146 | p = -25*(x1-2)^2-(x2-2)^2-(x3-1)^2-(x4-4)^2-(x5-1)^2-(x6-4)^2; |
---|
| 147 | F = set((x3-3)^2+x4>4)+set((x5-3)^2+x6>4); |
---|
| 148 | F = F + set(x1-3*x2<2)+set(-x1+x2<2); |
---|
| 149 | F = F + set(x1-3*x2 <2)+set(x1+x2>2); |
---|
| 150 | F = F + set(6>x1+x2>2); |
---|
| 151 | F = F + set(1<x3<5) + set(0<x4<6)+set(1<x5<5)+set(0<x6<10)+set(x1>0)+set(x2>0); |
---|
| 152 | options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk'); |
---|
| 153 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 154 | options = sdpsettings(options,'verbose',2,'solver','bmibnb'); |
---|
| 155 | solvesdp(F,p,options); |
---|
| 156 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 157 | * Lower solver : glpk |
---|
| 158 | * Upper solver : penbmi |
---|
| 159 | Node Upper Gap(%) Lower Open |
---|
| 160 | 1 : Inf NaN -3.130E+002 2 |
---|
| 161 | 2 : -3.100E+002 0.96 -3.130E+002 3 Improved solution |
---|
| 162 | + 2 Finishing. Cost: -310 Gap: 0.96463%</font></pre> |
---|
| 163 | </td> |
---|
| 164 | </tr> |
---|
| 165 | </table> |
---|
| 166 | <p>Quadratic equality constraints can also be dealt with. This can be used |
---|
| 167 | for, e.g., boolean programming. </p> |
---|
| 168 | <table cellPadding="10" width="100%"> |
---|
| 169 | <tr> |
---|
| 170 | <td class="xmpcode"> |
---|
| 171 | <pre>n = 10; |
---|
| 172 | x = sdpvar(n,1); |
---|
| 173 | |
---|
| 174 | Q = randn(n,n);Q = Q*Q'/norm(Q)^2; |
---|
| 175 | c = randn(n,1); |
---|
| 176 | |
---|
| 177 | objective = x'*Q*x+c'*x; |
---|
| 178 | |
---|
| 179 | F = set(x.*(x-1)==0) + set(0<x<1); |
---|
| 180 | options = sdpsettings('bmibnb.uppersolver','penbmi','bmibnb.lowersolver','glpk'); |
---|
| 181 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 182 | options = sdpsettings(options,'verbose',2,'solver','bmibnb'); |
---|
| 183 | solvesdp(F,objective,options); |
---|
| 184 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 185 | * Lower solver : glpk |
---|
| 186 | * Upper solver : penbmi |
---|
| 187 | Node Upper Gap(%) Lower Open |
---|
| 188 | 1 : -2.967E+000 0.00 -2.967E+000 2 Improved solution |
---|
| 189 | + 1 Finishing. Cost: -2.9671 Gap: 2.5207e-009%</font></pre> |
---|
| 190 | </td> |
---|
| 191 | </tr> |
---|
| 192 | </table> |
---|
| 193 | <h3>Nonconvex polynomial programming</h3> |
---|
| 194 | <p>The global solver in YALMIP can only solve bilinear programs, but |
---|
| 195 | a pre-processor is capable of transforming higher order problems to |
---|
| 196 | bilinear programs. As an example, the variable <b>x<sup>3</sup>y<sup>2</sup></b> |
---|
| 197 | will be replaced with the the variable <b>w</b> and the constraints |
---|
| 198 | <b>w == uv</b>, <b>u == zx</b>, <b>z == x<sup>2</sup></b>, <b>v == y<sup>2</sup></b>. |
---|
| 199 | The is done automatically if the global solver, or <a href="solvers.htm#penbmi">PENBMI</a>, |
---|
| 200 | is called with a higher order polynomial problem. Note that this |
---|
| 201 | conversion is rather inefficient, and only very small problems can |
---|
| 202 | be addressed using this simple approach. Additionally, this module |
---|
| 203 | has not been tested in any serious sense.</p> |
---|
| 204 | <table cellPadding="10" width="100%"> |
---|
| 205 | <tr> |
---|
| 206 | <td class="xmpcode"> |
---|
| 207 | <pre>sdpvar x y |
---|
| 208 | |
---|
| 209 | F = set(x^3+y^5 < 5) + set(y > 0); |
---|
| 210 | F = F + set(-100 < [x y] < 100) % Always bounded domain |
---|
| 211 | solvesdp(F,-x,options) |
---|
| 212 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 213 | * Lower solver : glpk |
---|
| 214 | * Upper solver : penbmi |
---|
| 215 | Node Upper Gap(%) Lower Open |
---|
| 216 | 1 : Inf NaN -2.990E+001 2 |
---|
| 217 | 2 : Inf NaN -2.990E+001 1 Infeasible |
---|
| 218 | 3 : Inf NaN -2.763E+001 2 |
---|
| 219 | 4 : Inf NaN -2.763E+001 3 |
---|
| 220 | 5 : Inf NaN -1.452E+001 4 |
---|
| 221 | 6 : Inf NaN -1.452E+001 5 |
---|
| 222 | 7 : -1.710E+000 171.54 -6.359E+000 4 Improved solution |
---|
| 223 | 8 : -1.710E+000 171.54 -6.359E+000 3 Infeasible |
---|
| 224 | 9 : -1.710E+000 127.11 -5.155E+000 4 Poor bound |
---|
| 225 | 10 : -1.710E+000 127.11 -5.155E+000 3 Infeasible |
---|
| 226 | 11 : -1.710E+000 0.00 -1.710E+000 2 Infeasible |
---|
| 227 | + 11 Finishing. Cost: -1.71 Gap: 1.247e-005%</font></pre> |
---|
| 228 | </td> |
---|
| 229 | </tr> |
---|
| 230 | </table> |
---|
| 231 | <h3>Nonconvex semidefinite programming</h3> |
---|
| 232 | <p>The following problem is a classical BMI problem</p> |
---|
| 233 | <table cellPadding="10" width="100%"> |
---|
| 234 | <tr> |
---|
| 235 | <td class="xmpcode"> |
---|
| 236 | <pre>yalmip('clear') |
---|
| 237 | x = sdpvar(1,1); |
---|
| 238 | y = sdpvar(1,1); |
---|
| 239 | t = sdpvar(1,1); |
---|
| 240 | A0 = [-10 -0.5 -2;-0.5 4.5 0;-2 0 0]; |
---|
| 241 | A1 = [9 0.5 0;0.5 0 -3 ; 0 -3 -1]; |
---|
| 242 | A2 = [-1.8 -0.1 -0.4;-0.1 1.2 -1;-0.4 -1 0]; |
---|
| 243 | K12 = [0 0 2;0 -5.5 3;2 3 0]; |
---|
| 244 | F = lmi(x>-0.5)+lmi(x<2) + lmi(y>-3)+lmi(y<7); |
---|
| 245 | F = F + lmi(A0+x*A1+y*A2+x*y*K12-t*eye(3)<0); |
---|
| 246 | options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi'); |
---|
| 247 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 248 | options = sdpsettings(options,'verbose',2,'solver','bmibnb'); |
---|
| 249 | solvesdp(F,t,options); |
---|
| 250 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 251 | * Lower solver : pensdp |
---|
| 252 | * Upper solver : penbmi |
---|
| 253 | Node Upper Gap(%) Lower Open |
---|
| 254 | 1 : -9.565E-001 2.22 -1.000E+000 2 Improved solution |
---|
| 255 | 2 : -9.565E-001 2.22 -1.000E+000 1 Infeasible |
---|
| 256 | 3 : -9.565E-001 2.22 -1.000E+000 2 |
---|
| 257 | 4 : -9.565E-001 2.22 -1.000E+000 3 |
---|
| 258 | 5 : -9.565E-001 0.24 -9.613E-001 2 Infeasible |
---|
| 259 | + 5 Finishing. Cost: -0.95653 Gap: 0.24126%</font></pre> |
---|
| 260 | </td> |
---|
| 261 | </tr> |
---|
| 262 | </table> |
---|
| 263 | <p>The constrained LQ problem in the section <a href="bmi.htm">Local BMIs</a> |
---|
| 264 | can actually easily be solved to global optimality in just a couple of |
---|
| 265 | iterations. For the global code to work, global lower and upper and bound on |
---|
| 266 | all complicating variables (involved in nonlinear terms) must be supplied, either explicitly or implicitly in the linear |
---|
| 267 | constraints. This was the case for all problems above. In this example, the variable <b> |
---|
| 268 | K</b> is already bounded in |
---|
| 269 | the original problem, but the elements in <b>P</b> have to be bounded |
---|
| 270 | artificially. </p> |
---|
| 271 | <table cellPadding="10" width="100%"> |
---|
| 272 | <tr> |
---|
| 273 | <td class="xmpcode"> |
---|
| 274 | <pre>yalmip('clear') |
---|
| 275 | A = [-1 2;-3 -4];B = [1;1]; |
---|
| 276 | P = sdpvar(2,2); |
---|
| 277 | K = sdpvar(1,2); |
---|
| 278 | F = set(P>0)+set((A+B*K)'*P+P*(A+B*K) < -eye(2)-K'*K); |
---|
| 279 | F = F + set(K<0.1)+set(K>-0.1); |
---|
| 280 | F = F + set(1000> P(:)>-1000) |
---|
| 281 | options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi'); |
---|
| 282 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 283 | options = sdpsettings(options,'verbose',2,'solver','bmibnb'); |
---|
| 284 | solvesdp(F,trace(P),options); |
---|
| 285 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 286 | * Lower solver : pensdp |
---|
| 287 | * Upper solver : penbmi |
---|
| 288 | Node Upper Gap(%) Lower Open |
---|
| 289 | 1 : 4.834E-001 17.70 2.209E-001 2 Improved solution |
---|
| 290 | 2 : 4.834E-001 17.70 2.209E-001 1 Infeasible |
---|
| 291 | 3 : 4.834E-001 1.93 4.548E-001 2 |
---|
| 292 | 4 : 4.834E-001 1.93 4.548E-001 1 Infeasible |
---|
| 293 | 5 : 4.834E-001 0.25 4.797E-001 2 |
---|
| 294 | + 5 Finishing. Cost: 0.48341 Gap: 0.25032%</font></pre> |
---|
| 295 | </td> |
---|
| 296 | </tr> |
---|
| 297 | </table> |
---|
| 298 | <p>Beware, the BMI solver is absolutely not that efficient in general, |
---|
| 299 | this was just a lucky case. Here is the <a href="gevp.htm">decay-rate |
---|
| 300 | example</a> instead (with some additional constraints on the elements |
---|
| 301 | in <b>P</b> to bound the feasible set). </p> |
---|
| 302 | <table cellPadding="10" width="100%"> |
---|
| 303 | <tr> |
---|
| 304 | <td class="xmpcode"> |
---|
| 305 | <pre>yalmip('clear'); |
---|
| 306 | A = [-1 2;-3 -4]; |
---|
| 307 | t = sdpvar(1,1); |
---|
| 308 | P = sdpvar(2,2); |
---|
| 309 | F = set(P>eye(2))+set(A'*P+P*A < -2*t*P); |
---|
| 310 | F = F + set(-1e4 < P(:) < 1e4) + set(100 > t > 0) ; |
---|
| 311 | options = sdpsettings('bmibnb.lowersolver','pensdp','bmibnb.uppersolver','penbmi'); |
---|
| 312 | options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk'); |
---|
| 313 | options = sdpsettings(options,'verbose',2,'solver','bmibnb'); |
---|
| 314 | solvesdp(F,-t,options); |
---|
| 315 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 316 | * Lower solver : pensdp |
---|
| 317 | * Upper solver : penbmi |
---|
| 318 | Node Upper Gap(%) Lower Open |
---|
| 319 | 1 : Inf NaN -9.965E+001 2 |
---|
| 320 | 2 : -2.500E+000 2775.83 -9.965E+001 3 Improved solution |
---|
| 321 | 3 : -2.500E+000 2721.00 -9.773E+001 4 |
---|
| 322 | 4 : -2.500E+000 2721.00 -9.773E+001 3 Infeasible |
---|
| 323 | 5 : -2.500E+000 2634.59 -9.471E+001 4 |
---|
| 324 | 6 : -2.500E+000 2634.59 -9.471E+001 3 Infeasible |
---|
| 325 | 7 : -2.500E+000 2577.61 -9.272E+001 4 |
---|
| 326 | 8 : -2.500E+000 2577.61 -9.272E+001 3 Infeasible |
---|
| 327 | 9 : -2.500E+000 2501.03 -9.004E+001 4 |
---|
| 328 | 10 : -2.500E+000 2501.03 -9.004E+001 3 Infeasible |
---|
| 329 | ... |
---|
| 330 | 94 : -2.500E+000 4.30 -2.650E+000 49 |
---|
| 331 | 95 : -2.500E+000 0.47 -2.516E+000 50 |
---|
| 332 | + 95 Finishing. Cost: -2.5 Gap: 0.46843%</font></pre> |
---|
| 333 | </td> |
---|
| 334 | </tr> |
---|
| 335 | </table> |
---|
| 336 | <p>The linear relaxations give very poor lower bounds on this problem, |
---|
| 337 | leading to poor convergence of the branching process. However, for this |
---|
| 338 | particular problem, the reason is easy to find. The original BMI is |
---|
| 339 | homogeneous w.r.t <b>P</b>, and to guarantee a somewhat reasonable |
---|
| 340 | solution, we artificially added the constraint <b>P<font face="Tahoma">≥I</font></b>. |
---|
| 341 | A better model is obtained if we instead fix the trace of <b>P</b>. This will |
---|
| 342 | make the feasible set w.r.t P much smaller, but the problems are |
---|
| 343 | equivalent. Note also that we no longer need any artificial constraints |
---|
| 344 | on the elements in <b>P</b>.</p> |
---|
| 345 | <table cellPadding="10" width="100%"> |
---|
| 346 | <tr> |
---|
| 347 | <td class="xmpcode"> |
---|
| 348 | <pre>F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ; |
---|
| 349 | F = F + set(trace(P)==1); |
---|
| 350 | solvesdp(F,-t,options); |
---|
| 351 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 352 | * Lower solver : pensdp |
---|
| 353 | * Upper solver : penbmi |
---|
| 354 | Node Upper Gap(%) Lower Open |
---|
| 355 | 1 : -2.500E+000 1396.51 -5.138E+001 2 Improved solution |
---|
| 356 | 2 : -2.500E+000 1396.51 -5.138E+001 1 Infeasible |
---|
| 357 | 3 : -2.500E+000 1.41 -2.549E+000 2 |
---|
| 358 | 4 : -2.500E+000 1.41 -2.549E+000 1 Infeasible |
---|
| 359 | 5 : -2.500E+000 0.17 -2.506E+000 2 |
---|
| 360 | + 5 Finishing. Cost: -2.5 Gap: 0.1746%</font></pre> |
---|
| 361 | </td> |
---|
| 362 | </tr> |
---|
| 363 | </table> |
---|
| 364 | <p>For this problem, we can easily find a very efficient additional cutting |
---|
| 365 | plane. The |
---|
| 366 | decay-rate BMI together with the constant trace on <b>P</b> implies <b>trace(A<sup>T</sup>P+PA)<font face="Times New Roman">≤</font>-2t</b>. |
---|
| 367 | Adding this cut reduce the number of iterations needed.</p> |
---|
| 368 | <table cellPadding="10" width="100%"> |
---|
| 369 | <tr> |
---|
| 370 | <td class="xmpcode"> |
---|
| 371 | <pre>F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0); |
---|
| 372 | F = F + set(trace(P)==1); |
---|
| 373 | F = F + set(trace(A'*P+P*A)<-2*t); |
---|
| 374 | solvesdp(F,-t,options); |
---|
| 375 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 376 | * Lower solver : pensdp |
---|
| 377 | * Upper solver : penbmi |
---|
| 378 | Node Upper Gap(%) Lower Open |
---|
| 379 | 1 : -2.500E+000 26.60 -3.431E+000 2 Improved solution |
---|
| 380 | 2 : -2.500E+000 26.60 -3.431E+000 3 |
---|
| 381 | 3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible |
---|
| 382 | + 3 Finishing. Cost: -2.5 Gap: 0.5461%</font></pre> |
---|
| 383 | </td> |
---|
| 384 | </tr> |
---|
| 385 | </table> |
---|
| 386 | <p> |
---|
| 387 | A Schur complement on the decay-rate BMI gives us yet another cut |
---|
| 388 | which improves the node relaxation even more.</p><table cellPadding="10" width="100%"> |
---|
| 389 | <tr> |
---|
| 390 | <td class="xmpcode"> |
---|
| 391 | <pre>F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ; |
---|
| 392 | F = F + set(trace(P)==1); |
---|
| 393 | F = F + set(trace(A'*P+P*A)<-2*t) |
---|
| 394 | F = F + set([-A'*P-P*A P*t;P*t P*t/2]); |
---|
| 395 | solvesdp(F,-t,options); |
---|
| 396 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 397 | * Lower solver : pensdp |
---|
| 398 | * Upper solver : penbmi |
---|
| 399 | Node Upper Gap(%) Lower Open |
---|
| 400 | 1 : -2.500E+000 17.84 -3.124E+000 2 Improved solution |
---|
| 401 | 2 : -2.500E+000 17.84 -3.124E+000 3 |
---|
| 402 | 3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible |
---|
| 403 | + 3 Finishing. Cost: -2.5 Gap: 0.55275%</font></pre> |
---|
| 404 | </td> |
---|
| 405 | </tr> |
---|
| 406 | </table> |
---|
| 407 | <p> |
---|
| 408 | By adding valid cuts, the relaxations are possibly tighter, leading to |
---|
| 409 | better lower bounds. A problem however is that we add additional |
---|
| 410 | burden to the local solver used for the upper bounds. The additional |
---|
| 411 | cuts are redundant for the local solver, and most likely detoriate the |
---|
| 412 | performance. To avoid this, cuts can be exlicitely specified by using |
---|
| 413 | the command <a target="topic" href="reference.htm#cut">cut</a>. |
---|
| 414 | Constraints defined using this command (instead of |
---|
| 415 | <a target="topic" href="reference.htm#set">set</a>) will only be used |
---|
| 416 | in the solution of relaxations, and omitted when the local solver is |
---|
| 417 | called. </p><table cellPadding="10" width="100%"> |
---|
| 418 | <tr> |
---|
| 419 | <td class="xmpcode"> |
---|
| 420 | <pre>F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ; |
---|
| 421 | F = F + set(trace(P)==1); |
---|
| 422 | F = F + cut(trace(A'*P+P*A)<-2*t); |
---|
| 423 | F = F + cut([-A'*P-P*A P*t;P*t P*t/2]); |
---|
| 424 | solvesdp(F,-t,options); |
---|
| 425 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 426 | * Lower solver : pensdp |
---|
| 427 | * Upper solver : penbmi |
---|
| 428 | Node Upper Gap(%) Lower Open |
---|
| 429 | 1 : -2.500E+000 17.84 -3.124E+000 2 Improved solution |
---|
| 430 | 2 : -2.500E+000 17.84 -3.124E+000 3 |
---|
| 431 | 3 : -2.500E+000 0.55 -2.519E+000 2 Infeasible |
---|
| 432 | + 3 Finishing. Cost: -2.5 Gap: 0.55275%</font></pre> |
---|
| 433 | </td> |
---|
| 434 | </tr> |
---|
| 435 | </table> |
---|
| 436 | <p> |
---|
| 437 | Upper bounds were obtained above by solving the BMI locally using |
---|
| 438 | <a href="solvers.htm#penbmi">PENBMI</a>. If no local BMI solver is available, an alternative is to |
---|
| 439 | check if the relaxed solution is a feasible solution. If so, |
---|
| 440 | the upper bound can be updated. This scheme |
---|
| 441 | can be obtained by specifying <code>'none'</code> as the upper bound solver.</p> |
---|
| 442 | <table cellPadding="10" width="100%"> |
---|
| 443 | <tr> |
---|
| 444 | <td class="xmpcode"> |
---|
| 445 | <pre>options = sdpsettings(options,'bmibnb.uppersolver','none'); |
---|
| 446 | F = set(P>0)+set(A'*P+P*A < -2*t*P) + set(100 > t > 0) ; |
---|
| 447 | F = F + set(trace(P)==1); |
---|
| 448 | F = F + cut(trace(A'*P+P*A)<-2*t); |
---|
| 449 | F = F + cut([-A'*P-P*A P*t;P*t P*t/2]); |
---|
| 450 | solvesdp(F,-t,options); |
---|
| 451 | <font color="#000000">* Starting YALMIP bilinear branch & bound. |
---|
| 452 | * Lower solver : pensdp |
---|
| 453 | * Upper solver : none |
---|
| 454 | Node Upper Gap(%) Lower Open |
---|
| 455 | 1 : Inf NaN -3.124E+000 2 |
---|
| 456 | 2 : Inf NaN -3.124E+000 3 |
---|
| 457 | 3 : -9.045E-001 104.47 -2.894E+000 4 Improved solution |
---|
| 458 | 4 : -9.045E-001 104.47 -2.894E+000 5 |
---|
| 459 | 5 : -1.467E+000 52.43 -2.761E+000 4 Improved solution |
---|
| 460 | 6 : -1.467E+000 52.43 -2.761E+000 5 |
---|
| 461 | 7 : -1.831E+000 29.87 -2.677E+000 4 Improved solution |
---|
| 462 | 8 : -1.831E+000 29.87 -2.677E+000 5 |
---|
| 463 | 9 : -2.071E+000 17.77 -2.616E+000 4 Improved solution |
---|
| 464 | 10 : -2.071E+000 17.77 -2.616E+000 5 |
---|
| 465 | 11 : -2.229E+000 10.46 -2.567E+000 4 Improved solution |
---|
| 466 | 12 : -2.229E+000 10.46 -2.567E+000 5 |
---|
| 467 | 13 : -2.332E+000 5.70 -2.522E+000 4 Improved solution |
---|
| 468 | 14 : -2.332E+000 5.70 -2.522E+000 5 |
---|
| 469 | 15 : -2.397E+000 3.25 -2.508E+000 4 Improved solution |
---|
| 470 | 16 : -2.397E+000 3.25 -2.508E+000 5 |
---|
| 471 | 17 : -2.435E+000 2.01 -2.504E+000 4 Improved solution |
---|
| 472 | 18 : -2.435E+000 2.01 -2.504E+000 5 |
---|
| 473 | 19 : -2.459E+000 1.25 -2.503E+000 4 Improved solution |
---|
| 474 | 20 : -2.459E+000 1.25 -2.503E+000 5 |
---|
| 475 | 21 : -2.478E+000 0.70 -2.502E+000 4 Improved solution |
---|
| 476 | + 21 Finishing. Cost: -2.4776 Gap: 0.69631%</font></pre> |
---|
| 477 | </td> |
---|
| 478 | </tr> |
---|
| 479 | </table> |
---|
| 480 | <p> |
---|
| 481 | More examples can be found in <a href="reference.htm#yalmipdemo"> |
---|
| 482 | yalmipdemo</a>.<br> |
---|
| 483 | <br> |
---|
| 484 | <img border="0" src="demoicon.gif" width="16" height="16"> The global |
---|
| 485 | branch and bound algorithm will hopefully be improved upon |
---|
| 486 | significantly in future releases. This includes both algorithmic |
---|
| 487 | improvements and code optimization. If you have any ideas, you are welcome to contribute.</td> |
---|
| 488 | </tr> |
---|
| 489 | </table> |
---|
| 490 | |
---|
| 491 | </div> |
---|
| 492 | |
---|
| 493 | </body> |
---|
| 494 | |
---|
| 495 | </html> |
---|