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

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

Added original make3d
File size: 4.0 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 : 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&nbsp;<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];
35P = 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&gt;eye(2))+set(A'*P+P*A &lt; -eye(2));
45solvesdp(F,trace(P));
46P0 = 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&nbsp;<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;
72F = set(P&gt;eye(2))+set(A'*P+P*A &lt; -2*t_upper*P);
73ops = sdpsettings('verbose',0,'warning',0);
74sol = solvesdp(F,[],ops);</pre>
75        <pre>while ~(sol.problem==1)
76&nbsp;&nbsp;&nbsp; t_upper = t_upper*2;
77&nbsp;&nbsp;&nbsp; F = set(P&gt;eye(2))+set(A'*P+P*A &lt; -2*t_upper*P);
78&nbsp;&nbsp;&nbsp; sol = solvesdp(F,[],ops);
79end</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;
89t_works = t_lower
90while (t_upper-t_lower)&gt;tol
91&nbsp; t_test = (t_upper+t_lower)/2;
92&nbsp; disp([t_lower t_upper t_test])
93&nbsp; F = set(P&gt;eye(2))+set(A'*P+P*A &lt; -2*t_test*P);
94&nbsp; sol = solvesdp(F,[],ops);
95&nbsp; if sol.problem==1
96&nbsp;&nbsp;&nbsp; t_upper = t_test;
97&nbsp; else
98&nbsp;&nbsp;&nbsp; t_lower = t_test;
99&nbsp;&nbsp;&nbsp; t_works = t_test;
100 end
101end</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>
Note: See TracBrowser for help on using the repository browser.