[37] | 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 |
---|