1 | clc |
---|
2 | echo on |
---|
3 | %********************************************************* |
---|
4 | % |
---|
5 | % Advanced YALMIP modeling, generalized eigenvalue problem |
---|
6 | % |
---|
7 | %********************************************************* |
---|
8 | % |
---|
9 | % The problem we solve here is estimatation of decay-rate |
---|
10 | % of a linear system x' = Ax. |
---|
11 | % |
---|
12 | % This can be formulated as |
---|
13 | % |
---|
14 | % max alpha |
---|
15 | % s.t A'P+PA < -2alphaP |
---|
16 | % P > I |
---|
17 | % |
---|
18 | pause |
---|
19 | |
---|
20 | % Define the variables |
---|
21 | A = [-1 2;-3 -4]; |
---|
22 | P = sdpvar(2,2); |
---|
23 | pause |
---|
24 | |
---|
25 | %To find a lower bound on alpha, we solve a standard Lyapunov stability problem. |
---|
26 | F = set(P>eye(2))+set(A'*P+P*A < -eye(2)); |
---|
27 | pause |
---|
28 | solvesdp(F,trace(P)); |
---|
29 | P0 = double(P) |
---|
30 | pause |
---|
31 | |
---|
32 | %For this particular solution, the decay-rate is |
---|
33 | alpha_lower = -max(eig(inv(P0)*(A'*P0+P0*A)))/2 |
---|
34 | pause |
---|
35 | |
---|
36 | % We now find an upper bound on the decay-rate by doubling alpha until |
---|
37 | % the problem is infeasible. To find out if the problem is infeasible, |
---|
38 | % we check the problem field in the solution structure. The |
---|
39 | % meaning of this variable is explained in the help text for the command |
---|
40 | % yalmiperror. Infeasibility is detected if the value is 1. To |
---|
41 | % reduce the amount of information written on the screen, we run the |
---|
42 | % solver in a completely silent mode. This can be accomplished by |
---|
43 | % using the verbose and warning options in sdpsettings. |
---|
44 | % |
---|
45 | % Note, on some solvers, just checking infeasibility flag is not enough. A |
---|
46 | % more careful bisection code should check also for numerical problems that |
---|
47 | % could indicate infeasibility |
---|
48 | pause |
---|
49 | |
---|
50 | options = sdpsettings('verbose',0,'warning',0); |
---|
51 | alpha_upper = alpha_lower*2; |
---|
52 | F = set(P>eye(2))+set(A'*P+P*A < -2*alpha_upper*P); |
---|
53 | sol = solvesdp(F,[],options); |
---|
54 | while ~(sol.problem==1) |
---|
55 | alpha_upper = alpha_upper*2; |
---|
56 | F = set(P>eye(2))+set(A'*P+P*A < -2*alpha_upper*P); |
---|
57 | sol = solvesdp(F,trace(P),options); |
---|
58 | end |
---|
59 | alpha_upper |
---|
60 | pause |
---|
61 | |
---|
62 | % Having both an upper bound and a lower bound allows us to perform a |
---|
63 | % bisection. (see code in decayex.m) |
---|
64 | pause |
---|
65 | echo off |
---|
66 | |
---|
67 | tol = 0.01; |
---|
68 | alpha_works = alpha_lower; |
---|
69 | clc |
---|
70 | disp('################################'); |
---|
71 | disp('# Lower Upper Current'); |
---|
72 | disp('################################'); |
---|
73 | while (alpha_upper-alpha_lower)>tol |
---|
74 | alpha_test = (alpha_upper+alpha_lower)/2; |
---|
75 | disp([alpha_lower alpha_upper alpha_test]) |
---|
76 | F = set(P>eye(2))+set(A'*P+P*A < -2*alpha_test*P); |
---|
77 | sol = solvesdp(F,trace(P),options); |
---|
78 | if sol.problem==1 |
---|
79 | alpha_upper = alpha_test; |
---|
80 | else |
---|
81 | alpha_lower = alpha_test; |
---|
82 | alpha_works = alpha_test; |
---|
83 | end |
---|
84 | end |
---|
85 | |
---|
86 | |
---|
87 | % A lower bound on the optimal value is thus |
---|
88 | alpha_works |
---|
89 | pause |
---|
90 | echo off |
---|