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

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

Added original make3d
File size: 29.5 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 : Multiparametric programming</title>
7<meta http-equiv="Content-Type" content="text/html; charset=windows-1251">
8<link href="yalmip.css" type="text/css" rel="stylesheet">
9<base target="_self">
10</head>
11
12<body leftMargin="0" topMargin="0">
13
14<div align="left">
15  <table border="0" cellpadding="4" cellspacing="3" style="border-collapse: collapse" bordercolor="#000000" width="100%" align="left" height="100%">
16    <tr>
17      <td width="100%" align="left" height="100%" valign="top">
18      <h2>Multiparametric programming</h2>
19      <hr noShade SIZE="1">
20    <p>
21    <img border="0" src="exclamationmark.jpg" align="left" width="16" height="16">This
22    example requires the <a href="solvers.htm#mpt">Multi-Parametric
23      Toolbox (MPT)</a>. </p>
24      <p>YALMIP can be used to calculate explicit solutions of linear and quadratic
25      programs by interfacing the <a href="solvers.htm#mpt">Multi-Parametric
26      Toolbox (MPT)</a>. This chapter assumes that the reader is familiar with
27                parametric programming and the basics of <a href="solvers.htm#mpt">MPT</a>.</p>
28                <h3>Standard example.</h3>
29                <p>Consider the following simple quadratic program in the<i> decision
30                variable</i> <b>z</b>, solved
31                for a particular value on a <i>parameter</i> <b>x</b>.</p>
32      <table cellPadding="10" width="100%" id="table1">
33        <tr>
34          <td class="xmpcode">
35          <pre>A = randn(15,3);
36b = rand(15,1);
37E = randn(15,2);
38
39z = sdpvar(3,1);
40x = [0.1;0.2];
41
42F = set(A*z &lt;= b+E*x);
43obj = (z-1)'*(z-1);
44
45sol = solvesdp(F,obj);
46double(z)
47<font color="#000000">ans =</font></pre>
48                        <pre><font color="#000000">   -0.1454
49   -0.1789
50   -0.0388</font></pre>
51          </td>
52        </tr>
53      </table>
54      <p>To obtain the parametric solution with respect to <b>x</b>, we call the
55                function <code>solvemp.m</code>, and tell the solver that <b>x</b> is a
56                parametric variable. Moreover, we must add constraints to declare a
57                bound on the region where we want to compute the parametric solution, the so called exploration set.</p>
58      <table cellPadding="10" width="100%">
59        <tr>
60          <td class="xmpcode">
61          <pre>x = sdpvar(2,1);
62F = set(A*z &lt;= b+E*x) + set(-1 &lt;= x &lt;= 1)
63sol = solvemp(F,obj,[],x);</pre>
64          </td>
65        </tr>
66      </table>
67      <p>The first output is an <a href="solvers.htm#mpt">MPT</a> structure.
68                In accordance with <a href="solvers.htm#mpt">MPT</a> syntax, the optimizer for the parametric value&nbsp; 
69                is given by the following code.</p>
70      <table cellPadding="10" width="100%" id="table2">
71        <tr>
72          <td class="xmpcode">
73          <pre>xx = [0.1;0.2];
74[i,j] = isinside(sol{1}.Pn,xx)
75sol{1}.Bi{j}*xx + sol{1}.Ci{j}
76<font color="#000000">ans =</font></pre>
77                        <pre><font color="#000000">   -0.1454
78   -0.1789
79   -0.0388</font></pre>
80          </td>
81        </tr>
82      </table>
83      <p>By using more outputs from <code>solvemp</code>, it is possible to simplify
84                things considerably.</p>
85      <table cellPadding="10" width="100%" id="table3">
86        <tr>
87          <td class="xmpcode">
88          <pre>[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,obj,[],x);</pre>
89          </td>
90        </tr>
91      </table>
92      <p>The function now returns solutions using YALMIPs
93                <a href="extoperators.htm">nonlinear operator</a> framework. To retrieve the solution for a particular
94                parameter, simply use <code>assign</code> and <code>double</code> in standard fashion.</p>
95      <table cellPadding="10" width="100%" id="table4">
96        <tr>
97          <td class="xmpcode">
98          <pre>assign(x,[0.1;0.2]);
99double(Optimizer)</pre>
100          </td>
101        </tr>
102      </table>
103      <p>Some of the plotting capabilities of <a href="solvers.htm#mpt">MPT</a> are overloaded for
104                the piecewise functions.</p>
105      <table cellPadding="10" width="100%" id="table25">
106        <tr>
107          <td class="xmpcode">
108          <pre>plot(Valuefunction)</pre>
109          </td>
110        </tr>
111      </table>
112      <p>More about the output later.</p>
113      <h3>MPC example</h3>
114                <p>Define numerical data for a linear system, prediction matrices, and
115                variables for current state <b>x</b> and the
116       
117        future control sequence<b> U</b>, for an MPC problem with horizon 5.</p>
118      <table cellPadding="10" width="100%">
119        <tr>
120          <td class="xmpcode">
121          <pre>N = 5;
122A = [2 -1;1 0];
123B = [1;0];
124C = [0.5 0.5];
125[H,S] = create_CHS(A,B,C,N);</pre>
126          <pre>x = sdpvar(2,1);
127U = sdpvar(N,1);</pre>
128          </td>
129        </tr>
130      </table>
131      <p>The future output predictions are linear in the current state&nbsp;and the
132      control sequence.</p>
133      <table cellPadding="10" width="100%">
134        <tr>
135          <td class="xmpcode">
136          <pre>Y = H*x+S*U;</pre>
137          </td>
138        </tr>
139      </table>
140      <p>We wish to minimize a  quadratic cost, compromising between small input
141                and outputs.</p>
142      <table cellPadding="10" width="100%">
143        <tr>
144          <td class="xmpcode">
145          <pre>objective = Y'*Y+U'*U;</pre>
146          </td>
147        </tr>
148      </table>
149      <p>The input variable has a hard constraint, and so does the output at the
150                terminal state.</p>
151      <table cellPadding="10" width="100%">
152        <tr>
153          <td class="xmpcode">
154          <pre>F = set(1 &gt; U &gt; -1);
155F = F+set(1 &gt; Y(N) &gt; -1); </pre>
156          </td>
157        </tr>
158      </table>
159      <p>We seek the explicit solution <b><font face="Tahoma,Arial,sans-serif">U(x)</font></b> over the
160                exploration set <b>
161      <font face="Tahoma,Arial,sans-serif">|x|&#8804;1</font></b>.</p>
162      <table cellPadding="10" width="100%">
163        <tr>
164          <td class="xmpcode">
165          <pre>F = F + set(5 &gt; x &gt; -5);</pre>
166          </td>
167        </tr>
168      </table>
169      <p>The explicit solution <b>U(x)</b> is obtained by calling
170      <code>solvemp</code> with the parametric variable <b>x</b> as the fourth argument.
171                Additionally, since we only are interested in the first element of the
172                solution <b>U</b>, we use a fifth input to communicate this.</p>
173      <table cellPadding="10" width="100%">
174        <tr>
175          <td class="xmpcode">
176          <pre>[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,objective,[],x,U(1));</pre>
177          </td>
178        </tr>
179      </table>
180      <p>The explicit solution can, e.g, be plotted (see the <a href="solvers.htm#mpt">MPT</a> 
181      manual informationon the solution structure)</p>
182      <table cellPadding="10" width="100%">
183        <tr>
184          <td class="xmpcode">
185          <pre>sol{1}
186<font color="#000000">
187ans =
188
189 Pn: [1x13 polytope]
190 Fi: {1x13 cell}
191 Gi: {1x13 cell}
192 activeConstraints: {1x13 cell}
193 Phard: [1x1 polytope]</font></pre>
194          <pre>figure
195plot(sol{1}.Pn)
196figure
197plot(Valuefunction)
198figure
199plot(Optimizer)</pre>
200          </td>
201        </tr>
202      </table>
203      <p>The following piece of code calculates the explicit MPC controller with
204      an worst-case L<sub><font face="Arial">&#8734;</font></sub> cost instead
205                (i.e. worst case control/output peak along the predictions).</p>
206      <table cellPadding="10" width="100%">
207        <tr>
208          <td class="xmpcode">
209          <pre>A = [2 -1;1 0];
210B = [1;0];
211C = [0.5 0.5];
212N = 5;
213[H,S] = create_CHS(A,B,C,N);
214t = sdpvar(1,1);
215x = sdpvar(2,1);
216U = sdpvar(N,1);
217Y = H*x+S*U;
218F = set(-1 &lt; U &lt; 1);
219F = F+set(-1 &lt; Y(N) &lt; 1);
220F = F + set(-5 &lt; x &lt; 5);
221[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,norm([Y;U],inf),[],x);</pre>
222          </td>
223        </tr>
224      </table>
225      <p>This example enables us to introduce the overloading of the projection
226      functionality in <a href="solvers.htm#mpt">MPT</a>. In MPC,
227      only the first input, <b>U(1)</b>, is of interest. What we can do is to
228      project the whole problem to the parametric variable <b>x</b>, the
229      objective function <b>t</b>, and the variable of interest, <b>U(1)</b>.
230                The following piece of code projects the problem to the reduced set of
231                variables, and solves the multi-parametric LP in this reduced space (not
232                necessarily an efficient approach though, projection is very expensive,
233                and we can easily end up with a problem with many constraints in the
234                reduced variables.)</p>
235      <table cellPadding="10" width="100%">
236        <tr>
237          <td class="xmpcode">
238          <pre>A = [2 -1;1 0];
239B = [1;0];
240C = [0.5 0.5];
241N = 5;
242[H,S] = create_CHS(A,B,C,N);
243t = sdpvar(1,1);
244x = sdpvar(2,1);
245U = sdpvar(N,1);
246Y = H*x+S*U;
247F = set(-1 &lt; U &lt; 1);
248F = F+set(-1 &lt; Y(N) &lt; 1);
249F = F + set(-5 &lt; x &lt; 5);
250F = F + set(-t &lt; [Y;U] &lt; t);
251
252F = projection(F,[x;t;U(1)]);
253
254[sol_p,diagnostics_p,Zsol_p,Valuefunction_p,Optimizer_p] = solvemp(F,t,[],x,[U(1);t]);</pre>
255          </td>
256        </tr>
257      </table>
258      <p>To check that the two solutions actually coincide, we compare the value
259                functions and conclude that they are essentially the same.</p>
260      <table cellPadding="10" width="100%" id="table30">
261        <tr>
262          <td class="xmpcode">
263          <pre>[pass,tol] = mpt_ispwabigger(sol{1},sol_p{1})
264pass =
265    <font color="#000000"> 0</font>
266tol =
267 <font color="#000000"> 1.3484e-014
268
269</font>figure;
270plot(Valuefunction);
271plot(Valuefunction_p):</pre>
272          </td>
273        </tr>
274      </table>
275      <p>Since the objective function is <b>t</b>, the optimal t should coincide
276                with the value function. The optimal <b>t</b> in the projected problem
277                was requested as the second optimizer</p>
278      <table cellPadding="10" width="100%" id="table31">
279        <tr>
280          <td class="xmpcode">
281          <pre>plot(Optimizer(2));</pre>
282          </td>
283        </tr>
284      </table>
285      <p>In the code above, we used the pre-constructed function <code>create_CHS.m</code>,
286                which generates numerical data for defining predictions from
287                current state and future input signals. The problem can however be
288                solved without using this function. We will show two alternative
289                methods. </p>
290                <p>The first approach explicitly generates future states using the given
291                current state and future inputs (note that the solution not is exactly
292                the same as above due to a small difference in indexing logic.)</p>
293      <table cellPadding="10" width="100%" id="table32">
294        <tr>
295          <td class="xmpcode">
296          <pre>x  = sdpvar(2,1);
297U  = sdpvar(N,1);</pre>
298                        <pre>x0 = x; % Save parametric state
299F  =  set(-5 &lt; x0 &lt; 5); % Exploration space
300objective= 0;
301for k = 1:N   
302
303    % Feasible region
304    F = F + set(-1 &lt; U(k) &lt; 1);</pre>
305                        <pre>    % Add stage cost to total cost
306    objective = objective + (C*x)'*(C*x) + U(k)'*U(k);
307
308    % Explicit state update
309    x = A*x + B*U(k); 
310end
311F = F + set(-1 &lt; C*x &lt; 1); % Terminal constraint
312
313[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,objective,[],x0,U(1));</pre>
314          </td>
315        </tr>
316      </table>
317      <p>The second method introduced new decision variables for both future
318                states and connect these using equality constraints.</p>
319      <table cellPadding="10" width="100%" id="table33">
320        <tr>
321          <td class="xmpcode">
322          <pre>x  = sdpvar(2,N+1);
323U  = sdpvar(N,1);</pre>
324                        <pre>F  =  set(-5 &lt; x0 &lt; 5); % Exploration space
325objective = 0;
326for k = 1:N
327   
328    % Feasible region
329    F = F + set(-1 &lt; U(k) &lt; 1);
330 
331    % Add stage cost to total cost
332    objective = objective + (C*x(:,k))'*(C*x(:,k)) + U(k)'*U(k);</pre>
333                        <pre>    % Implicit state update
334    F = F + set(x(:,k+1) == A*x(:,k) + B*U(k)); 
335end
336F = F + set(-1 &lt; C*x(:,end) &lt; 1); % Terminal constraint
337
338[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,objective,[],x(:,1),U(1));</pre>
339          </td>
340        </tr>
341      </table>
342      <p>This second method is less efficient, since it has a lot more decision
343                variables. However, in some cases, it is the most convenient way to
344                model a system (see the PWA examples below for similar approaches), and
345                the additional decision variables are removed anyway during a pre-solve
346                process YALMIP applies before calling the multi-parametric solver.</p>
347      <h3>Parametric programs with binary variables and equality constraints</h3>
348                <p>YALMIP extends the parametric algorithms in <a href="solvers.htm#mpt">MPT</a> 
349                by adding a layer to
350                enable binary variables and equality constraints. We can use this to
351                find explicit solutions to, e.g., predictive control of PWA (piecewise
352                affine) systems.</p>
353                <p>Let us find the explicit solution to a predictive control problem,
354                where the gain of the system depends on the sign of the first state.
355                This will be a pretty advanced example, so let us start slowly by defining the
356                some data.</p>
357      <table cellPadding="10" width="100%" id="table6">
358        <tr>
359          <td class="xmpcode">
360          <pre>yalmip('clear')
361clear all</pre>
362                        <pre>% Model data
363A = [2 -1;1 0];
364B1 = [1;0];  % Small gain for x(1) &gt; 0
365B2 = [pi;0]; % Larger gain for x(1) &lt; 0
366C = [0.5 0.5];</pre>
367                        <pre>nx = 2; % Number of states
368nu = 1; % Number of inputs
369
370% Prediction horizon
371N = 4;</pre>
372          </td>
373        </tr>
374      </table>
375      <p>To simplify the code and the notation, we create a bunch of state and
376                control vectors in cell arrays</p>
377      <table cellPadding="10" width="100%" id="table5">
378        <tr>
379          <td class="xmpcode">
380          <pre>% States x(k), ..., x(k+N)
381x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre>
382                        <pre>% Inputs u(k), ..., u(k+N) (last one not used)
383u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre>
384                        <pre>% Binary for PWA selection
385d = binvar(repmat(2,1,N),repmat(1,1,N));</pre>
386          </td>
387        </tr>
388      </table>
389      <p>We now run a loop to add constraints on all states and inputs. Note
390                the use of binary variables define the PWA dynamics, and the recommended
391                use of <a href="logic.htm#bounds">bounds</a> to improve big-M
392                relaxations)</p>
393      <table cellPadding="10" width="100%" id="table7">
394        <tr>
395          <td class="xmpcode">
396          <pre>F   = set([]);
397obj = 0;
398
399for k = N-1:-1:1   
400
401    % Strenghten big-M (improves numerics)
402    bounds(x{k},-5,5);
403    bounds(u{k},-1,1);
404    bounds(x{k+1},-5,5);
405   
406    % Feasible region
407    F = F + set(-1 &lt; u{k}     &lt; 1);
408    F = F + set(-1 &lt; C*x{k}   &lt; 1);
409    F = F + set(-5 &lt; x{k}     &lt; 5);
410    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
411    F = F + set(-5 &lt; x{k+1}   &lt; 5);</pre>
412                        <pre>    % PWA Dynamics
413    F = F + set(implies(d{k}(1),x{k+1} == (A*x{k}+B1*u{k})));
414    F = F + set(implies(d{k}(2),x{k+1} == (A*x{k}+B2*u{k})));
415    F = F + set(implies(d{k}(1),x{k}(1) &gt; 0));
416    F = F + set(implies(d{k}(2),x{k}(1) &lt; 0));
417
418    % It is EXTREMELY important to add as many
419    % constraints as possible to the binary variables
420    F = F + set(sum(d{k}) == 1);
421 
422    % Add stage cost to total cost
423    obj = obj + norm(x{k},1) + norm(u{k},1);
424
425end</pre>
426          </td>
427        </tr>
428      </table>
429      <p>The parametric variable here is the current state <b>x{1}</b>. In this
430                optimization problem, there are a lot of variables that we have no
431                interest in. To tell YALMIP that we only want the optimizer for the
432                current state <b>u{1}</b>, we use a fifth input argument.</p>
433      <table cellPadding="10" width="100%" id="table8">
434        <tr>
435          <td class="xmpcode">
436          <pre>[sol,diagnostics,Z,Valuefunction,Optimizer] = solvemp(F,obj,[],x{1},u{1});</pre>
437          </td>
438        </tr>
439      </table>
440      <p>To obtain the optimal control input for a specific state, we use <code>double</code> 
441                and <code>assign</code> as usual.</p>
442      <table cellPadding="10" width="100%" id="table9">
443        <tr>
444          <td class="xmpcode">
445          <pre>assign(x{1},[-1;1]);
446double(Optimizer)
447ans =</pre>
448                        <pre>    0.9549</pre>
449          </td>
450        </tr>
451      </table>
452      <p>The optimal cost at this state is available in the value function</p>
453      <table cellPadding="10" width="100%" id="table10">
454        <tr>
455          <td class="xmpcode">
456          <pre>double(Valuefunction)
457ans =</pre>
458                        <pre>    4.2732</pre>
459          </td>
460        </tr>
461      </table>
462      <p>To convince our self that we have a correct parametric solution, let us
463                compare it to the solution obtained by solving the problem for this
464                specific state. </p>
465      <table cellPadding="10" width="100%" id="table11">
466        <tr>
467          <td class="xmpcode">
468          <pre>sol = solvesdp(F+set(x{1}==[-1;1]),obj);
469double(u{1})
470<font color="#000000">ans =</font></pre>
471                        <pre><font color="#000000">    0.9549</font></pre>
472                        <pre>double(obj)
473<font color="#000000">ans =</font></pre>
474                        <pre><font color="#000000">    4.2732</font></pre>
475          </td>
476        </tr>
477      </table>
478      <h3><a name="dp"></a>Dynamic programming with LTI systems</h3>
479                <p>The capabilities in YALMIP to work with piecewise functions and
480                parametric programs enables easy coding of dynamic programming
481                algorithms. The value function with respect to the parametric variable
482                for a parametric linear program is a convex PWA function, and this is
483                the function returned in the fourth output. YALMIP creates this function
484                internally, and also saves information about convexity etc, and uses it
485                as any other <a href="extoperators.htm">nonlinear operator</a> (see more
486                details below). For
487                binary parametric linear programs, the value function is no longer
488                convex, but a so called overlapping PWA function. This means that, at
489                each point, it is defined as the minimum of a set of convex PWA
490                function. This information is also handled transparently in YALMIP, it
491                is simply another type of <a href="extoperators.htm">nonlinear operator</a>.
492                The main difference between the two function classes is that the second
493                class requires introduction of binary variables when used.</p>
494                <p>Note, the algorithms described in the following sections are mainly
495                intended for with (picewise) linear objectives. Dynamic programming with
496                quadratic
497                objective functions give rise to problems that are much harder to solve,
498                although it is supported.</p>
499                <p>To illustrate how easy it is to work with these PWA functions, we can
500                solve predictive control using dynamic programming, instead of setting
501                up the whole problem in one shot as we did above. As a first example, we solve a
502                standard linear predictive control problem.&nbsp; To fully understand
503                this example, it is required that you are familiar with predictive
504                control, dynamic programming and parametric optimization.</p>
505      <table cellPadding="10" width="100%" id="table12">
506        <tr>
507          <td class="xmpcode">
508          <pre>yalmip('clear')
509clear all</pre>
510                        <pre>% Model data
511A = [2 -1;1 0];
512B = [1;0];
513C = [0.5 0.5];</pre>
514                        <pre>nx = 2; % Number of states
515nu = 1; % Number of inputs
516
517% Prediction horizon
518N = 4;</pre>
519          <pre>% States x(k), ..., x(k+N)
520x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre>
521                        <pre>% Inputs u(k), ..., u(k+N) (last one not used)
522u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre>
523          </td>
524        </tr>
525      </table>
526      <p>Now, instead of setting up the problem for the whole prediction
527                horizon, we only set it up for one step, solve the problem
528                parametrically, take one step back, and perform a standard dynamic
529                programming value iteration.</p>
530      <table cellPadding="10" width="100%" id="table13">
531        <tr>
532          <td class="xmpcode">
533          <pre>J{N} = 0;
534
535for k = N-1:-1:1   
536   
537    % Feasible region
538    F =     set(-1 &lt; u{k}     &lt; 1);
539    F = F + set(-1 &lt; C*x{k}   &lt; 1);
540    F = F + set(-5 &lt; x{k}     &lt; 5);
541    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
542    F = F + set(-5 &lt; x{k+1}   &lt; 5);   </pre>
543                        <pre>    % Dynamics
544    F = F + set(x{k+1} == A*x{k}+B*u{k});</pre>
545                        <pre>    % Cost in value iteration
546    obj = norm(x{k},1) + norm(u{k},1) + J{k+1}
547
548    % Solve one-step problem   
549    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
550end</pre>
551          </td>
552        </tr>
553      </table>
554      <p>Notice the minor changes needed compared to the one shot solution.
555                Important to understand is that the value function at step <b>k</b> will
556                be a function in <b>x{k}</b>, hence when it is used at <b>k-1</b>, it
557                will be a function penalizing the predicted state. Note that YALMIP
558                automatically keep track of convexity of the value function. Hence, for
559                this example, no binary variables are introduced along the solution
560                process (this is why we safely can
561                omit the use of <a href="reference.htm#bounds">bounds</a>).</p>
562                <p>To study the development of the value function, we can easily plot
563                them. </p>
564      <table cellPadding="10" width="100%" id="table14">
565        <tr>
566          <td class="xmpcode">
567          <pre>for k = N-1:-1:1
568 plot(J{k});hold on
569end</pre>
570          </td>
571        </tr>
572      </table>
573      <p>The optimal control input can also be plotted.</p>
574      <table cellPadding="10" width="100%" id="table15">
575        <tr>
576          <td class="xmpcode">
577          <pre>clf;plot(Optimizer{1})</pre>
578          </td>
579        </tr>
580      </table>
581      <p>Of course, you can do the same thing by working directly with the
582                numerical MPT
583                data.</p>
584      <table cellPadding="10" width="100%" id="table23">
585        <tr>
586          <td class="xmpcode">
587          <pre>for k = N-1:-1:1
588 mpt_plotPWA(sol{k}{1}.Pn,sol{k}{1}.Bi,sol{k}{1}.Ci);hold on
589end</pre>
590          </td>
591        </tr>
592      </table>
593      <p>Why not solve the problem with a worst-case L<sub><font face="Arial">&#8734;</font></sub> cost!
594                Notice&nbsp;the non-additive value iteration! Although this might look
595                complicated using a mix of maximums, norms and piecewise value-functions
596                , the formulation fits perfectly into the <a href="extoperators.htm">nonlinear operator</a> 
597                framework in YALMIP.</p>
598      <table cellPadding="10" width="100%" id="table24">
599        <tr>
600          <td class="xmpcode">
601          <pre>J{N} = 0;
602
603for k = N-1:-1:1   
604   
605    % Feasible region
606    F =     set(-1 &lt; u{k}     &lt; 1);
607    F = F + set(-1 &lt; C*x{k}   &lt; 1);
608    F = F + set(-5 &lt; x{k}     &lt; 5);
609    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
610    F = F + set(-5 &lt; x{k+1}   &lt; 5);   </pre>
611                        <pre>    % Dynamics
612    F = F + set(x{k+1} == A*x{k}+B*u{k});</pre>
613                        <pre>    % Cost in value iteration
614    obj = max(norm([x{k};u{k}],inf),J{k+1})
615
616    % Solve one-step problem   
617    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
618end</pre>
619          </td>
620        </tr>
621      </table>
622      <h3>Dynamic programming with PWA systems</h3>
623                <p>As mentioned above, the difference between the value function from a
624                parametric LP and a parametric LP with binary variables is that
625                convexity of the value function no longer is guaranteed. In practice,
626                this means that (additional) binary variables are required when the the
627                value function is used in an optimization problem. This is done
628                automatically in YALMIP through the <a href="extoperators.htm">nonlinear
629                operator</a> framework, but it is extremely important to realize that
630                the value function from a binary problem is much more complex than the
631                one resulting from a standard parametric problem. Nevertheless, working
632                with the objects is easy, as the following example illustrates. In this
633                case, we solve the dynamic programming problem for our PWA system</p>
634      <table cellPadding="10" width="100%" id="table20">
635        <tr>
636          <td class="xmpcode">
637          <pre>yalmip('clear')
638clear all</pre>
639                        <pre>% Model data
640A = [2 -1;1 0];
641B1 = [1;0];
642B2 = [pi;0]; % Larger gain for negative first state
643C = [0.5 0.5];</pre>
644                        <pre>nx = 2; % Number of states
645nu = 1; % Number of inputs
646
647% Prediction horizon
648N = 4;</pre>
649          <pre>% States x(k), ..., x(k+N)
650x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre>
651                        <pre>% Inputs u(k), ..., u(k+N) (last one not used)
652u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre>
653                        <pre>% Binary for PWA selection
654d = binvar(repmat(2,1,N),repmat(1,1,N));</pre>
655          </td>
656        </tr>
657      </table>
658      <p>Just as in the LTI case, we set up the problem for the one step case,
659                and use the value function from the previous iteration.</p>
660      <table cellPadding="10" width="100%" id="table21">
661        <tr>
662          <td class="xmpcode">
663          <pre>J{N} = 0;
664
665for k = N-1:-1:1   
666
667    % Strenghten big-M (improves numerics)
668    bounds(x{k},-5,5);
669    bounds(u{k},-1,1);
670    bounds(x{k+1},-5,5);
671   
672    % Feasible region
673    F =     set(-1 &lt; u{k}     &lt; 1);
674    F = F + set(-1 &lt; C*x{k}   &lt; 1);
675    F = F + set(-5 &lt; x{k}     &lt; 5);
676    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
677    F = F + set(-5 &lt; x{k+1}   &lt; 5);  </pre>
678                        <pre>    % PWA Dynamics
679    F = F + set(implies(d{k}(1),x{k+1} == (A*x{k}+B1*u{k})));
680    F = F + set(implies(d{k}(2),x{k+1} == (A*x{k}+B2*u{k})));
681    F = F + set(implies(d{k}(1),x{k}(1) &gt; 0));
682    F = F + set(implies(d{k}(2),x{k}(1) &lt; 0));
683    F = F + set(sum(d{k}) == 1);</pre>
684                        <pre>    % Cost in value iteration
685    obj = norm(x{k},1) + norm(u{k},1) + J{k+1}</pre>
686                        <pre>    % Solve one-step problem
687    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
688end</pre>
689          </td>
690        </tr>
691      </table>
692      <p>Plotting the final value function clearly reveals the nonconvexity.</p>
693      <table cellPadding="10" width="100%" id="table22">
694        <tr>
695          <td class="xmpcode">
696          <pre>clf;plot(J{1})</pre>
697          </td>
698        </tr>
699      </table>
700      <h3>Behind the scenes and advanced use</h3>
701                <p>The first thing that might be a bit unusual to the advanced user is
702                the piecewise functions that YALMIP returns in the fourth and fifth
703                output from <code>solvemp</code> . In principle, they are specialized
704                <a href="reference.htm#pwf">pwf</a> objects. To create a PWA value
705                function after solving a multi-parametric LP, the following
706                command is used.</p>
707      <table cellPadding="10" width="100%" id="table26">
708        <tr>
709          <td class="xmpcode">
710          <pre>Valuefunction = pwf(sol,x,'convex')</pre>
711          </td>
712        </tr>
713      </table>
714      <p>The <a href="reference.htm#pwf">pwf</a> command will recognize the MPT
715                solution structure and create a PWA function based on the
716                fields <b>Pn</b>, <b>Bi</b> and <b>Ci</b>. The dedicated
717                <a href="extoperators.htm">nonlinear operator</a> is implemented in the file <code>pwa_yalmip.m</code> 
718                . The <a href="extoperators.htm">nonlinear operator</a> will exploit the
719                fact that the PWA function is convex and implements an efficient epi-graph
720                representation. In case the PWA function is used in a non-convex fashion
721                (i.e. the automatic <a href="extoperators.htm#propagation">convexity
722                propagation</a> fails), a
723                <a href="extoperators.htm#milp">MILP implementation</a> is also
724                available.<p>If the field <b>Ai</b> is non-empty (solution obtained from
725                a multi-parametric QP), a corresponding PWQ function is created (<code>pwq_yalmip.m</code> 
726                ). <p>To create a PWA function
727                representing the optimizer, two things have to be changed. To begin
728                with, YALMIP searches for the <b>Bi</b> and <b>Ci</b> fields, but since
729                we want to create a PWA function based on <b>Fi</b> and <b>Gi</b> 
730                fields, the field names have to be changed. Secondly, the piecewise&nbsp; 
731                optimizer is typically not convex, so a general PWA function is created
732                instead (requiring 1 binary variable per region if the variable later
733                actually is used in an optimization problem.)<table cellPadding="10" width="100%" id="table27">
734        <tr>
735          <td class="xmpcode">
736          <pre>sol.Ai = cell(1,length(sol.Ai));
737sol.Bi = sol.Fi;
738sol.Ci = sol.Gi;
739Optimizer = pwf(sol,x,'general')</pre>
740          </td>
741        </tr>
742      </table>
743      <p>A third important case is when the solution structure returned from <code>solvemp</code> is
744                a cell with several <a href="solvers.htm#mpt">MPT</a> structures. This means that a
745                multi-parametric problem with binary variables was
746                solved, and the different cells represent overlapping solutions. One way
747                to get rid of the overlapping solutions is to use the MPT command
748                <code>mpt_removeOverlaps.m</code> and create a PWA function based on the result. Since the
749                resulting PWA function typically is non-convex,
750                we must create a general function.<table cellPadding="10" width="100%" id="table28">
751        <tr>
752          <td class="xmpcode">
753          <pre>sol_single = mpt_removeOverlaps(sol);
754Valuefunction = pwf(sol_single,x,'general')</pre>
755          </td>
756        </tr>
757      </table>
758                <p>This is only recommended if you just intend to plot or investigate
759                the value function. Typically, if you want to continue computing with the
760                result, i.e. use it in another optimization problem, as in the dynamic
761                programming examples above, it is recommended to keep the overlapping
762                solutions. To model a general PWA function, 1 binary variable per region
763                is needed. However, by using a dedicated <a href="extoperators.htm">nonlinear operator</a> 
764                for overlapping convex PWA functions, only one binary per PWA function
765                is needed.</p>
766                <table cellPadding="10" width="100%" id="table29">
767        <tr>
768          <td class="xmpcode">
769          <pre>Valuefunction = pwf(sol,x,'overlapping')</pre>
770          </td>
771        </tr>
772      </table>
773                <p>As we have seen in the examples above, the PWA and PWQ functions can be
774                plotted. This is currently nothing but a simple hack. Normally, when you
775                apply the plot command, the corresponding double values are plotted.
776                However, if the input is a simple scalar PWA or PWQ variable, the
777                underlying MPT structures will be extracted and the plot commands in MPT
778                will be called.</td>
779    </tr>
780  </table>
781        <p>&nbsp;</div>
782
783</body>
784
785</html>
Note: See TracBrowser for help on using the repository browser.