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> |
---|