source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/htmldata/moment.htm @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 14.1 KB
Line 
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
45obj = -2*x1+x2-x3;
46F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24&gt;0);
47F = F + set(4-(x1+x2+x3)&gt;0);
48F = F + set(6-(3*x2+x3)&gt;0);
49F = F + set(x1&gt;0);
50F = F + set(2-x1&gt;0);
51F = F + set(x2&gt;0);
52F = F + set(x3&gt;0);
53F = F + set(3-x3&gt;0);
54solvemoment(F,obj);
55double(obj)
56<font color="#000000">
57 ans =</font></pre>
58        <pre><font color="#000000">&nbsp;&nbsp;&nbsp;&nbsp; -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);
91double(obj)
92
93<font color="#000000"> ans =</font></pre>
94        <pre><font color="#000000">&nbsp;&nbsp;&nbsp;&nbsp; -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&gt;double(obj)),obj,[],2);
104double(obj)
105
106 <font color="#000000">ans =</font></pre>
107        <pre><font color="#000000">&nbsp;&nbsp;&nbsp;&nbsp; -5.3870</font></pre>
108        <pre>solvemoment(F+set(obj&gt;double(obj)),obj,[],2);
109double(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);
122double(obj)
123
124<font color="#000000"> ans =</font></pre>
125        <pre><font color="#000000">&nbsp;&nbsp;&nbsp;&nbsp; -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);
162x{1}
163<font color="#000000">
164ans =
165
166    0.5000
167    0.0000
168    3.0000</font>
169x{2}
170<font color="#000000">
171ans =
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});
186double(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
214p = 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)));
224sdisplay(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);
237t = sosdec.t;
238v = sosdec.v0;
239Q = sosdec.Q0;
240sdisplay(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
255obj = -x1^2-x2^2;
256F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
257[sol,x] = solvemoment(F,obj,[],2);
258assign([x1;x2],x{1});
259double(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
283obj = -x1^2-x2^2;
284F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
285ops = sdpsettings('moment.refine',5','moment.rceftol',1e-8);
286[sol,xe,data] = solvemoment(F,obj,ops);
287xopt1 = extractsolution(data,sdpsettings('moment.refine',0));
288xopt2 = extractsolution(data,sdpsettings('moment.refine',1));
289xopt3 = extractsolution(data,sdpsettings('moment.refine',10));
290xopt4 = 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
300obj = -x1^2-x2^2;
301F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
302sol = solvesdp(F,obj,sdpsettings('solver','moment','moment.order',2));
303assign(sol.momentdata.x,sol.xoptimal{1});
304double(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>&nbsp;</div>
338
339</body>
340
341</html>
Note: See TracBrowser for help on using the repository browser.