[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 : 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); |
---|
| 36 | b = rand(15,1); |
---|
| 37 | E = randn(15,2); |
---|
| 38 | |
---|
| 39 | z = sdpvar(3,1); |
---|
| 40 | x = [0.1;0.2]; |
---|
| 41 | |
---|
| 42 | F = set(A*z <= b+E*x); |
---|
| 43 | obj = (z-1)'*(z-1); |
---|
| 44 | |
---|
| 45 | sol = solvesdp(F,obj); |
---|
| 46 | double(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); |
---|
| 62 | F = set(A*z <= b+E*x) + set(-1 <= x <= 1) |
---|
| 63 | sol = 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 |
---|
| 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) |
---|
| 75 | sol{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]); |
---|
| 99 | double(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; |
---|
| 122 | A = [2 -1;1 0]; |
---|
| 123 | B = [1;0]; |
---|
| 124 | C = [0.5 0.5]; |
---|
| 125 | [H,S] = create_CHS(A,B,C,N);</pre> |
---|
| 126 | <pre>x = sdpvar(2,1); |
---|
| 127 | U = sdpvar(N,1);</pre> |
---|
| 128 | </td> |
---|
| 129 | </tr> |
---|
| 130 | </table> |
---|
| 131 | <p>The future output predictions are linear in the current state 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 > U > -1); |
---|
| 155 | F = F+set(1 > Y(N) > -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|≤1</font></b>.</p> |
---|
| 162 | <table cellPadding="10" width="100%"> |
---|
| 163 | <tr> |
---|
| 164 | <td class="xmpcode"> |
---|
| 165 | <pre>F = F + set(5 > x > -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"> |
---|
| 187 | ans = |
---|
| 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 |
---|
| 195 | plot(sol{1}.Pn) |
---|
| 196 | figure |
---|
| 197 | plot(Valuefunction) |
---|
| 198 | figure |
---|
| 199 | plot(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">∞</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]; |
---|
| 210 | B = [1;0]; |
---|
| 211 | C = [0.5 0.5]; |
---|
| 212 | N = 5; |
---|
| 213 | [H,S] = create_CHS(A,B,C,N); |
---|
| 214 | t = sdpvar(1,1); |
---|
| 215 | x = sdpvar(2,1); |
---|
| 216 | U = sdpvar(N,1); |
---|
| 217 | Y = H*x+S*U; |
---|
| 218 | F = set(-1 < U < 1); |
---|
| 219 | F = F+set(-1 < Y(N) < 1); |
---|
| 220 | F = F + set(-5 < x < 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]; |
---|
| 239 | B = [1;0]; |
---|
| 240 | C = [0.5 0.5]; |
---|
| 241 | N = 5; |
---|
| 242 | [H,S] = create_CHS(A,B,C,N); |
---|
| 243 | t = sdpvar(1,1); |
---|
| 244 | x = sdpvar(2,1); |
---|
| 245 | U = sdpvar(N,1); |
---|
| 246 | Y = H*x+S*U; |
---|
| 247 | F = set(-1 < U < 1); |
---|
| 248 | F = F+set(-1 < Y(N) < 1); |
---|
| 249 | F = F + set(-5 < x < 5); |
---|
| 250 | F = F + set(-t < [Y;U] < t); |
---|
| 251 | |
---|
| 252 | F = 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}) |
---|
| 264 | pass = |
---|
| 265 | <font color="#000000"> 0</font> |
---|
| 266 | tol = |
---|
| 267 | <font color="#000000"> 1.3484e-014 |
---|
| 268 | |
---|
| 269 | </font>figure; |
---|
| 270 | plot(Valuefunction); |
---|
| 271 | plot(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); |
---|
| 297 | U = sdpvar(N,1);</pre> |
---|
| 298 | <pre>x0 = x; % Save parametric state |
---|
| 299 | F = set(-5 < x0 < 5); % Exploration space |
---|
| 300 | objective= 0; |
---|
| 301 | for k = 1:N |
---|
| 302 | |
---|
| 303 | % Feasible region |
---|
| 304 | F = F + set(-1 < U(k) < 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); |
---|
| 310 | end |
---|
| 311 | F = F + set(-1 < C*x < 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); |
---|
| 323 | U = sdpvar(N,1);</pre> |
---|
| 324 | <pre>F = set(-5 < x0 < 5); % Exploration space |
---|
| 325 | objective = 0; |
---|
| 326 | for k = 1:N |
---|
| 327 | |
---|
| 328 | % Feasible region |
---|
| 329 | F = F + set(-1 < U(k) < 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)); |
---|
| 335 | end |
---|
| 336 | F = F + set(-1 < C*x(:,end) < 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') |
---|
| 361 | clear all</pre> |
---|
| 362 | <pre>% Model data |
---|
| 363 | A = [2 -1;1 0]; |
---|
| 364 | B1 = [1;0]; % Small gain for x(1) > 0 |
---|
| 365 | B2 = [pi;0]; % Larger gain for x(1) < 0 |
---|
| 366 | C = [0.5 0.5];</pre> |
---|
| 367 | <pre>nx = 2; % Number of states |
---|
| 368 | nu = 1; % Number of inputs |
---|
| 369 | |
---|
| 370 | % Prediction horizon |
---|
| 371 | N = 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) |
---|
| 381 | x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre> |
---|
| 382 | <pre>% Inputs u(k), ..., u(k+N) (last one not used) |
---|
| 383 | u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre> |
---|
| 384 | <pre>% Binary for PWA selection |
---|
| 385 | d = 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([]); |
---|
| 397 | obj = 0; |
---|
| 398 | |
---|
| 399 | for 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 < u{k} < 1); |
---|
| 408 | F = F + set(-1 < C*x{k} < 1); |
---|
| 409 | F = F + set(-5 < x{k} < 5); |
---|
| 410 | F = F + set(-1 < C*x{k+1} < 1); |
---|
| 411 | F = F + set(-5 < x{k+1} < 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) > 0)); |
---|
| 416 | F = F + set(implies(d{k}(2),x{k}(1) < 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 | |
---|
| 425 | end</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]); |
---|
| 446 | double(Optimizer) |
---|
| 447 | ans =</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) |
---|
| 457 | ans =</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); |
---|
| 469 | double(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. 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') |
---|
| 509 | clear all</pre> |
---|
| 510 | <pre>% Model data |
---|
| 511 | A = [2 -1;1 0]; |
---|
| 512 | B = [1;0]; |
---|
| 513 | C = [0.5 0.5];</pre> |
---|
| 514 | <pre>nx = 2; % Number of states |
---|
| 515 | nu = 1; % Number of inputs |
---|
| 516 | |
---|
| 517 | % Prediction horizon |
---|
| 518 | N = 4;</pre> |
---|
| 519 | <pre>% States x(k), ..., x(k+N) |
---|
| 520 | x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre> |
---|
| 521 | <pre>% Inputs u(k), ..., u(k+N) (last one not used) |
---|
| 522 | u = 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 | |
---|
| 535 | for k = N-1:-1:1 |
---|
| 536 | |
---|
| 537 | % Feasible region |
---|
| 538 | F = set(-1 < u{k} < 1); |
---|
| 539 | F = F + set(-1 < C*x{k} < 1); |
---|
| 540 | F = F + set(-5 < x{k} < 5); |
---|
| 541 | F = F + set(-1 < C*x{k+1} < 1); |
---|
| 542 | F = F + set(-5 < x{k+1} < 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}); |
---|
| 550 | end</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 |
---|
| 569 | end</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 |
---|
| 589 | end</pre> |
---|
| 590 | </td> |
---|
| 591 | </tr> |
---|
| 592 | </table> |
---|
| 593 | <p>Why not solve the problem with a worst-case L<sub><font face="Arial">∞</font></sub> cost! |
---|
| 594 | Notice 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 | |
---|
| 603 | for k = N-1:-1:1 |
---|
| 604 | |
---|
| 605 | % Feasible region |
---|
| 606 | F = set(-1 < u{k} < 1); |
---|
| 607 | F = F + set(-1 < C*x{k} < 1); |
---|
| 608 | F = F + set(-5 < x{k} < 5); |
---|
| 609 | F = F + set(-1 < C*x{k+1} < 1); |
---|
| 610 | F = F + set(-5 < x{k+1} < 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}); |
---|
| 618 | end</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') |
---|
| 638 | clear all</pre> |
---|
| 639 | <pre>% Model data |
---|
| 640 | A = [2 -1;1 0]; |
---|
| 641 | B1 = [1;0]; |
---|
| 642 | B2 = [pi;0]; % Larger gain for negative first state |
---|
| 643 | C = [0.5 0.5];</pre> |
---|
| 644 | <pre>nx = 2; % Number of states |
---|
| 645 | nu = 1; % Number of inputs |
---|
| 646 | |
---|
| 647 | % Prediction horizon |
---|
| 648 | N = 4;</pre> |
---|
| 649 | <pre>% States x(k), ..., x(k+N) |
---|
| 650 | x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre> |
---|
| 651 | <pre>% Inputs u(k), ..., u(k+N) (last one not used) |
---|
| 652 | u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre> |
---|
| 653 | <pre>% Binary for PWA selection |
---|
| 654 | d = 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 | |
---|
| 665 | for 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 < u{k} < 1); |
---|
| 674 | F = F + set(-1 < C*x{k} < 1); |
---|
| 675 | F = F + set(-5 < x{k} < 5); |
---|
| 676 | F = F + set(-1 < C*x{k+1} < 1); |
---|
| 677 | F = F + set(-5 < x{k+1} < 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) > 0)); |
---|
| 682 | F = F + set(implies(d{k}(2),x{k}(1) < 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}); |
---|
| 688 | end</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 |
---|
| 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)); |
---|
| 737 | sol.Bi = sol.Fi; |
---|
| 738 | sol.Ci = sol.Gi; |
---|
| 739 | Optimizer = 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); |
---|
| 754 | Valuefunction = 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> </div> |
---|
| 782 | |
---|
| 783 | </body> |
---|
| 784 | |
---|
| 785 | </html> |
---|