[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 : Moment relaxations</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>Moment based relaxation of polynomials programs</h2> |
---|
| 22 | <hr noShade SIZE="1"> |
---|
| 23 | <p>YALMIP comes with a built-in module for polynomial programming using |
---|
| 24 | moment relaxations. This can be used for finding lower bounds on constrained |
---|
| 25 | polynomial programs (inequalities and equalities, element-wise and semidefinite), and to extract the corresponding optimizers. The implementation is entirely based on high-level |
---|
| 26 | YALMIP code, and can be somewhat inefficient for large problems (the |
---|
| 27 | inefficiency would then show in the setup of the problem, not actually |
---|
| 28 | solving the semidefinite resulting program). For very large problems, you might |
---|
| 29 | want to check out the dedicated software package |
---|
| 30 | <a target="_blank" href="http://www.laas.fr/~henrion/software/gloptipoly/">Gloptipoly</a> |
---|
| 31 | (the solution time will be the same, but the setup time might be reduced). |
---|
| 32 | For the underlying theory of moment relaxations, the reader is referred to |
---|
| 33 | <a href="readmore.htm#LASSERRE">[Lasserre]</a>.</p> |
---|
| 34 | <h3>Solving polynomial problems by relaxations</h3> |
---|
| 35 | <p>The following code calculates a lower bound on a concave quadratic |
---|
| 36 | optimization problem. As you can see, the only difference compared to |
---|
| 37 | solving the problem using a standard solver, such as |
---|
| 38 | <a href="solvers.htm#fmincon">fmincon</a> or <a href="solvers.htm#penbmi"> |
---|
| 39 | penbmi</a>, is that we call <a href="reference.htm#solvemoment">solvemoment</a> instead of |
---|
| 40 | <a href="reference.htm#solvesdp">solvesdp</a>.</p> |
---|
| 41 | <table cellPadding="10" width="100%"> |
---|
| 42 | <tr> |
---|
| 43 | <td class="xmpcode"> |
---|
| 44 | <pre>sdpvar x1 x2 x3 |
---|
| 45 | obj = -2*x1+x2-x3; |
---|
| 46 | F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0); |
---|
| 47 | F = F + set(4-(x1+x2+x3)>0); |
---|
| 48 | F = F + set(6-(3*x2+x3)>0); |
---|
| 49 | F = F + set(x1>0); |
---|
| 50 | F = F + set(2-x1>0); |
---|
| 51 | F = F + set(x2>0); |
---|
| 52 | F = F + set(x3>0); |
---|
| 53 | F = F + set(3-x3>0); |
---|
| 54 | solvemoment(F,obj); |
---|
| 55 | double(obj) |
---|
| 56 | <font color="#000000"> |
---|
| 57 | ans =</font></pre> |
---|
| 58 | <pre><font color="#000000"> -6.0000</font></pre> |
---|
| 59 | </td> |
---|
| 60 | </tr> |
---|
| 61 | </table> |
---|
| 62 | <p>Notice that YALMIP does not recover variables by default, a fact showing up in the |
---|
| 63 | difference between lifted variables and actual nonlinear variables (lifted |
---|
| 64 | variables are the variables used in the semidefinite relaxation to model |
---|
| 65 | nonlinear variables) The lifted variables can be obtained by using the |
---|
| 66 | command <code>relaxdouble</code> . The quadratic |
---|
| 67 | constraint above is satisfied in the lifted variables, but not in the true |
---|
| 68 | variables, as the following code illustrates.</p> |
---|
| 69 | <table cellPadding="10" width="100%"> |
---|
| 70 | <tr> |
---|
| 71 | <td class="xmpcode"> |
---|
| 72 | <pre>relaxdouble(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24) |
---|
| 73 | |
---|
| 74 | <font color="#000000">ans = |
---|
| 75 | |
---|
| 76 | 23.2648</font></pre> |
---|
| 77 | <pre>double(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24) |
---|
| 78 | |
---|
| 79 | <font color="#000000"> ans = |
---|
| 80 | |
---|
| 81 | -2.0000</font></pre> |
---|
| 82 | </td> |
---|
| 83 | </tr> |
---|
| 84 | </table> |
---|
| 85 | <p>An tighter relaxation can be obtained by using a higher order relaxation (the |
---|
| 86 | lowest possible is used if it is not specified).</p> |
---|
| 87 | <table cellPadding="10" width="100%"> |
---|
| 88 | <tr> |
---|
| 89 | <td class="xmpcode"> |
---|
| 90 | <pre>solvemoment(F,obj,[],2); |
---|
| 91 | double(obj) |
---|
| 92 | |
---|
| 93 | <font color="#000000"> ans =</font></pre> |
---|
| 94 | <pre><font color="#000000"> -5.6593</font></pre> |
---|
| 95 | </td> |
---|
| 96 | </tr> |
---|
| 97 | </table> |
---|
| 98 | <p>The obtained bound can be used iteratively to improve the bound by adding |
---|
| 99 | dynamically generated cuts.</p> |
---|
| 100 | <table cellPadding="10" width="100%"> |
---|
| 101 | <tr> |
---|
| 102 | <td class="xmpcode"> |
---|
| 103 | <pre>solvemoment(F+set(obj>double(obj)),obj,[],2); |
---|
| 104 | double(obj) |
---|
| 105 | |
---|
| 106 | <font color="#000000">ans =</font></pre> |
---|
| 107 | <pre><font color="#000000"> -5.3870</font></pre> |
---|
| 108 | <pre>solvemoment(F+set(obj>double(obj)),obj,[],2); |
---|
| 109 | double(obj) |
---|
| 110 | |
---|
| 111 | <font color="#000000">ans = |
---|
| 112 | |
---|
| 113 | -5.1270</font></pre> |
---|
| 114 | </td> |
---|
| 115 | </tr> |
---|
| 116 | </table> |
---|
| 117 | <p>The known true minima, -4, is found in the fourth order relaxation.</p> |
---|
| 118 | <table cellPadding="10" width="100%"> |
---|
| 119 | <tr> |
---|
| 120 | <td class="xmpcode"> |
---|
| 121 | <pre>solvemoment(F,obj,[],4); |
---|
| 122 | double(obj) |
---|
| 123 | |
---|
| 124 | <font color="#000000"> ans =</font></pre> |
---|
| 125 | <pre><font color="#000000"> -4.0000</font></pre> |
---|
| 126 | </td> |
---|
| 127 | </tr> |
---|
| 128 | </table> |
---|
| 129 | <p>The true global minima is however not recovered with the lifted variables, as we can see if we |
---|
| 130 | check the current solution (still violates the nonlinear constraint).</p> |
---|
| 131 | <table cellPadding="10" width="100%"> |
---|
| 132 | <tr> |
---|
| 133 | <td class="xmpcode"> |
---|
| 134 | <pre>checkset(F) |
---|
| 135 | <font color="#000000"> |
---|
| 136 | +++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 137 | | ID| Constraint | Type| Primal residual| |
---|
| 138 | +++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 139 | | #1| Numeric value| Element-wise| </font><font color="#FF0000">-0.88573</font><font color="#000000">| |
---|
| 140 | | #2| Numeric value| Element-wise| 1.834| |
---|
| 141 | | #3| Numeric value| Element-wise| 5.668| |
---|
| 142 | | #4| Numeric value| Element-wise| 1.834| |
---|
| 143 | | #5| Numeric value| Element-wise| 0.16599| |
---|
| 144 | | #6| Numeric value| Element-wise| 2.0873e-006| |
---|
| 145 | | #7| Numeric value| Element-wise| 0.33198| |
---|
| 146 | | #8| Numeric value| Element-wise| 2.668| |
---|
| 147 | +++++++++++++++++++++++++++++++++++++++++++++++++++</font></pre> |
---|
| 148 | </td> |
---|
| 149 | </tr> |
---|
| 150 | </table> |
---|
| 151 | <h3>Extracting solutions</h3> |
---|
| 152 | <p>To extract a (or several) globally optimal solution, we need two output |
---|
| 153 | arguments. The first output is a diagnostic structure (standard solution |
---|
| 154 | structure from the semidefinite solver), the second output is the |
---|
| 155 | (hopefully) extracted globally optimal solutions and the third output is |
---|
| 156 | a data structure containing all data that was needed to extract the |
---|
| 157 | solution.</p> |
---|
| 158 | <table cellPadding="10" width="100%"> |
---|
| 159 | <tr> |
---|
| 160 | <td class="xmpcode"> |
---|
| 161 | <pre>[sol,x] = solvemoment(F,obj,[],4); |
---|
| 162 | x{1} |
---|
| 163 | <font color="#000000"> |
---|
| 164 | ans = |
---|
| 165 | |
---|
| 166 | 0.5000 |
---|
| 167 | 0.0000 |
---|
| 168 | 3.0000</font> |
---|
| 169 | x{2} |
---|
| 170 | <font color="#000000"> |
---|
| 171 | ans = |
---|
| 172 | |
---|
| 173 | 2.0000 |
---|
| 174 | -0.0000 |
---|
| 175 | 0.0000</font></pre> |
---|
| 176 | </td> |
---|
| 177 | </tr> |
---|
| 178 | </table> |
---|
| 179 | <p>We can easily check that these are globally optimal solutions since |
---|
| 180 | they reach the lower bound -4 and satisfy the constraints (up to |
---|
| 181 | numerical precision).</p> |
---|
| 182 | <table cellPadding="10" width="100%"> |
---|
| 183 | <tr> |
---|
| 184 | <td class="xmpcode"> |
---|
| 185 | <pre>assign([x1;x2;x3],x{1}); |
---|
| 186 | double(obj) |
---|
| 187 | <font color="#000000">ans =</font></pre> |
---|
| 188 | <pre><font color="#000000"> -4.0000</font></pre> |
---|
| 189 | <pre>checkset(F) |
---|
| 190 | <font color="#000000">+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 191 | | ID| Constraint| Type| Primal residual| |
---|
| 192 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 193 | | #1| Numeric value| Element-wise (quadratic)| -1.1034e-006| |
---|
| 194 | | #2| Numeric value| Element-wise| 0.5| |
---|
| 195 | | #3| Numeric value| Element-wise| 3| |
---|
| 196 | | #4| Numeric value| Element-wise| 0.5| |
---|
| 197 | | #5| Numeric value| Element-wise| 1.5| |
---|
| 198 | | #6| Numeric value| Element-wise| 5.956e-007| |
---|
| 199 | | #7| Numeric value| Element-wise| 3| |
---|
| 200 | | #8| Numeric value| Element-wise| 4.6084e-007| |
---|
| 201 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++</font></pre> |
---|
| 202 | </td> |
---|
| 203 | </tr> |
---|
| 204 | </table> |
---|
| 205 | <h3>Extracting sum of squares decompositions.</h3> |
---|
| 206 | <p>Moment based approaches and <a href="sos.htm">sum of squares |
---|
| 207 | computations</a> are directly linked through duality. Given the solution |
---|
| 208 | to the moment problem, a sum of squares decomposition is easily derived |
---|
| 209 | from dual variables. Consider a simple unconstrained problem.</p> |
---|
| 210 | <table cellPadding="10" width="100%" id="table1"> |
---|
| 211 | <tr> |
---|
| 212 | <td class="xmpcode"> |
---|
| 213 | <pre>sdpvar x1 x2 |
---|
| 214 | p = x1^4+x2^4-x1*x2+1</pre> |
---|
| 215 | </td> |
---|
| 216 | </tr> |
---|
| 217 | </table> |
---|
| 218 | <p>A sum of squares decomposition of <b>p</b> can be obtained by using |
---|
| 219 | <a href="reference.htm#solvesos">solvesos</a></p> |
---|
| 220 | <table cellPadding="10" width="100%" id="table2"> |
---|
| 221 | <tr> |
---|
| 222 | <td class="xmpcode"> |
---|
| 223 | <pre>[sol,v,Q] = solvesos(set(sos(p))); |
---|
| 224 | sdisplay(clean(p-v{1}'*Q{1}*v{1},1e-6)) |
---|
| 225 | <font color="#000000">0</font></pre> |
---|
| 226 | </td> |
---|
| 227 | </tr> |
---|
| 228 | </table> |
---|
| 229 | <p>The decomposition can alternatively be computed from the moment |
---|
| 230 | results. If we minimize the polynomial using moments, |
---|
| 231 | <a href="reference.htm#solvemoment">solvemoment</a> will automatically |
---|
| 232 | extract <b>t</b>, <b>v</b> and <b>Q</b> in the decomposition <b>p(x)-t=v<sup>T</sup>Qv</b></p> |
---|
| 233 | <table cellPadding="10" width="100%" id="table3"> |
---|
| 234 | <tr> |
---|
| 235 | <td class="xmpcode"> |
---|
| 236 | <pre>[sol,x,moments,sosdec] = solvemoment([],p); |
---|
| 237 | t = sosdec.t; |
---|
| 238 | v = sosdec.v0; |
---|
| 239 | Q = sosdec.Q0; |
---|
| 240 | sdisplay(clean(p-t-v'*Q*v,1e-6)) |
---|
| 241 | <font color="#000000">0</font></pre> |
---|
| 242 | </td> |
---|
| 243 | </tr> |
---|
| 244 | </table> |
---|
| 245 | <p>In the more general constrained case, the polynomial multipliers for |
---|
| 246 | the constraints will also be returned.</p> |
---|
| 247 | <h3><a name="matrix"></a>Polynomial semidefinite constraints</h3> |
---|
| 248 | <p>Nonlinear semidefinite constraints can be |
---|
| 249 | added using exactly the same notation and syntax. The following example is taken from |
---|
| 250 | <a href="readmore.htm#HENRIONLASSERRE">[D. Henrion, J. B. Lasserre]</a>.</p> |
---|
| 251 | <table cellPadding="10" width="100%"> |
---|
| 252 | <tr> |
---|
| 253 | <td class="xmpcode"> |
---|
| 254 | <pre>sdpvar x1 x2 |
---|
| 255 | obj = -x1^2-x2^2; |
---|
| 256 | F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]); |
---|
| 257 | [sol,x] = solvemoment(F,obj,[],2); |
---|
| 258 | assign([x1;x2],x{1}); |
---|
| 259 | double(obj) |
---|
| 260 | <font color="#000000">ans =</font></pre> |
---|
| 261 | <pre><font color="#000000"> -4.00003</font></pre> |
---|
| 262 | <pre>checkset(F) |
---|
| 263 | <font color="#000000">++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 264 | | ID| Constraint| Type| Primal residual| |
---|
| 265 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 266 | | #1| Numeric value| LMI (quadratic)| -0.00034633| |
---|
| 267 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++</font></pre> |
---|
| 268 | </td> |
---|
| 269 | </tr> |
---|
| 270 | </table> |
---|
| 271 | <h3>Advanced features</h3> |
---|
| 272 | <p>A number of advanced features are available. We just briefly introduce |
---|
| 273 | these here by a quick example where we refine the extracted solution |
---|
| 274 | using a couple of Newton steps on an algebraic systems defining the global |
---|
| 275 | solutions given the optimal moment matrices, and change the numerical tolerance in the extraction |
---|
| 276 | algorithm. Finally, we calculate some different global solutions using |
---|
| 277 | the optimal moment matrices. Please check the moment relaxation example in |
---|
| 278 | <a href="reference.htm#yalmipdemo">yalmipdemo</a> for details.</p> |
---|
| 279 | <table cellPadding="10" width="100%"> |
---|
| 280 | <tr> |
---|
| 281 | <td class="xmpcode"> |
---|
| 282 | <pre>sdpvar x1 x2 |
---|
| 283 | obj = -x1^2-x2^2; |
---|
| 284 | F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]); |
---|
| 285 | ops = sdpsettings('moment.refine',5','moment.rceftol',1e-8); |
---|
| 286 | [sol,xe,data] = solvemoment(F,obj,ops); |
---|
| 287 | xopt1 = extractsolution(data,sdpsettings('moment.refine',0)); |
---|
| 288 | xopt2 = extractsolution(data,sdpsettings('moment.refine',1)); |
---|
| 289 | xopt3 = extractsolution(data,sdpsettings('moment.refine',10)); |
---|
| 290 | xopt4 = extractsolution(data,sdpsettings('moment.rceftol',1e-3,'moment.refine',5));</pre> |
---|
| 291 | </td> |
---|
| 292 | </tr> |
---|
| 293 | </table> |
---|
| 294 | <p>The moment relaxation solver can also be called using a more standard |
---|
| 295 | YALMIP notation, by simply defining the solver as <code>'moment'</code>.</p> |
---|
| 296 | <table cellPadding="10" width="100%"> |
---|
| 297 | <tr> |
---|
| 298 | <td class="xmpcode"> |
---|
| 299 | <pre>sdpvar x1 x2 |
---|
| 300 | obj = -x1^2-x2^2; |
---|
| 301 | F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]); |
---|
| 302 | sol = solvesdp(F,obj,sdpsettings('solver','moment','moment.order',2)); |
---|
| 303 | assign(sol.momentdata.x,sol.xoptimal{1}); |
---|
| 304 | double(obj) |
---|
| 305 | <font color="#000000">ans =</font></pre> |
---|
| 306 | <pre><font color="#000000"> -4.00003</font></pre> |
---|
| 307 | <pre>checkset(F) |
---|
| 308 | <font color="#000000">++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 309 | | ID| Constraint| Type| Primal residual| |
---|
| 310 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
| 311 | | #1| Numeric value| LMI (quadratic)| -0.00034633| |
---|
| 312 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++</font></pre> |
---|
| 313 | </td> |
---|
| 314 | </tr> |
---|
| 315 | </table> |
---|
| 316 | <p> |
---|
| 317 | <img border="0" src="demoicon.gif" width="16" height="16"> The |
---|
| 318 | semidefinite programs rapidly grows when the number of variables and the |
---|
| 319 | polynomial degree increase, so be careful when you model your problem.</p> |
---|
| 320 | <p> |
---|
| 321 | <img border="0" src="demoicon.gif" width="16" height="16"> The two |
---|
| 322 | most useful practical tips when working semidefinite relaxations seem to |
---|
| 323 | be to de-symmetrize your objective function, and compactify your |
---|
| 324 | feasible region. These two tricks typically increase the likelihood that |
---|
| 325 | you will be able to extract global solutions. By adding a perturbation |
---|
| 326 | to the polynomial, symmetry is lost, which generically means that there |
---|
| 327 | will not be infinitely many optima, and the extraction algorithm is more |
---|
| 328 | likely to work. Most theory in moment relaxations assumes that the |
---|
| 329 | feasible set is compact, and this is also affecting practical |
---|
| 330 | performance. By adding redundant compactifying constraints, you |
---|
| 331 | typically increase the likelihood of success. As an example, a simple |
---|
| 332 | redundant constraint which often work well in practice is an upper bound on |
---|
| 333 | the objective function based on a known feasible sub-optimal solution.</td> |
---|
| 334 | </tr> |
---|
| 335 | </table> |
---|
| 336 | |
---|
| 337 | <p> </div> |
---|
| 338 | |
---|
| 339 | </body> |
---|
| 340 | |
---|
| 341 | </html> |
---|