[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 : Generalized eigenvalue problems</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>Decay-rate estimation</h2> |
---|
| 22 | <hr noShade SIZE="1"> |
---|
| 23 | <p>The problem we will solve is to estimate the decay-rate of a linear |
---|
| 24 | system <strong>x' = Ax</strong>. This can be formulated as a generalized |
---|
| 25 | eigenvalue problem.</p> |
---|
| 26 | <p><img border="0" src="gevp.h4.gif" hspace="45"></p> |
---|
| 27 | <p>Due to the product between <b><font face="Tahoma">t</font></b> and |
---|
| 28 | <b><font face="Tahoma">P</font></b>, the problem cannot be solved |
---|
| 29 | directly. However, it is easily solved by bisection in <b>t</b>.</p> |
---|
| 30 | <p>Define the variables.</p> |
---|
| 31 | <table cellPadding="10" width="100%"> |
---|
| 32 | <tr> |
---|
| 33 | <td class="xmpcode"> |
---|
| 34 | <pre>A = [-1 2;-3 -4]; |
---|
| 35 | P = sdpvar(2,2);</pre> |
---|
| 36 | </td> |
---|
| 37 | </tr> |
---|
| 38 | </table> |
---|
| 39 | <p>To find a lower bound on <b>t</b>, we solve a standard Lyapunov stability |
---|
| 40 | problem.</p> |
---|
| 41 | <table cellPadding="10" width="100%"> |
---|
| 42 | <tr> |
---|
| 43 | <td class="xmpcode"> |
---|
| 44 | <pre>F = set(P>eye(2))+set(A'*P+P*A < -eye(2)); |
---|
| 45 | solvesdp(F,trace(P)); |
---|
| 46 | P0 = double(P);</pre> |
---|
| 47 | </td> |
---|
| 48 | </tr> |
---|
| 49 | </table> |
---|
| 50 | <p>In the code above, we minimized the trace just to get a numerically sound |
---|
| 51 | solution. This solution gives us a lower bound on decay-rate</p> |
---|
| 52 | <table cellPadding="10" width="100%"> |
---|
| 53 | <tr> |
---|
| 54 | <td class="xmpcode"> |
---|
| 55 | <pre>t_lower = -max(eig(inv(P0)*(A'*P0+P0*A)))/2;</pre> |
---|
| 56 | </td> |
---|
| 57 | </tr> |
---|
| 58 | </table> |
---|
| 59 | <p>We now find an upper bound on the decay-rate by doubling t until the |
---|
| 60 | problem is infeasible. To find out if the problem is infeasible, we check |
---|
| 61 | the field <code>problem</code> in the solution structure. The meaning of |
---|
| 62 | this variable is explained in the help text for the command <a href="reference.htm#yalmiperror">yalmiperror</a>. Infeasibility has been detected by the solver if the value is 1. To reduce |
---|
| 63 | the amount of information written on the screen, we run the solver in a |
---|
| 64 | completely silent mode. This can be accomplished by using the <code>verbose</code> |
---|
| 65 | and <code>warning</code> options in |
---|
| 66 | <a href="reference.htm#sdpsettings"> |
---|
| 67 | sdpsettings</a>.</p> |
---|
| 68 | <table cellPadding="10" width="100%"> |
---|
| 69 | <tr> |
---|
| 70 | <td class="xmpcode"> |
---|
| 71 | <pre>t_upper = t_lower*2; |
---|
| 72 | F = set(P>eye(2))+set(A'*P+P*A < -2*t_upper*P); |
---|
| 73 | ops = sdpsettings('verbose',0,'warning',0); |
---|
| 74 | sol = solvesdp(F,[],ops);</pre> |
---|
| 75 | <pre>while ~(sol.problem==1) |
---|
| 76 | t_upper = t_upper*2; |
---|
| 77 | F = set(P>eye(2))+set(A'*P+P*A < -2*t_upper*P); |
---|
| 78 | sol = solvesdp(F,[],ops); |
---|
| 79 | end</pre> |
---|
| 80 | </td> |
---|
| 81 | </tr> |
---|
| 82 | </table> |
---|
| 83 | <p>Having both an upper bound and a lower bound allows us to perform a |
---|
| 84 | bisection.</p> |
---|
| 85 | <table cellPadding="10" width="100%"> |
---|
| 86 | <tr> |
---|
| 87 | <td class="xmpcode"> |
---|
| 88 | <pre>tol = 0.01; |
---|
| 89 | t_works = t_lower |
---|
| 90 | while (t_upper-t_lower)>tol |
---|
| 91 | t_test = (t_upper+t_lower)/2; |
---|
| 92 | disp([t_lower t_upper t_test]) |
---|
| 93 | F = set(P>eye(2))+set(A'*P+P*A < -2*t_test*P); |
---|
| 94 | sol = solvesdp(F,[],ops); |
---|
| 95 | if sol.problem==1 |
---|
| 96 | t_upper = t_test; |
---|
| 97 | else |
---|
| 98 | t_lower = t_test; |
---|
| 99 | t_works = t_test; |
---|
| 100 | end |
---|
| 101 | end</pre> |
---|
| 102 | </td> |
---|
| 103 | </tr> |
---|
| 104 | </table> |
---|
| 105 | <p>This example will be revisited later when we study <a href="bmi.htm"> |
---|
| 106 | BMIs</a></td> |
---|
| 107 | </tr> |
---|
| 108 | </table> |
---|
| 109 | |
---|
| 110 | </div> |
---|
| 111 | |
---|
| 112 | </body> |
---|
| 113 | |
---|
| 114 | </html> |
---|