Decay-rate estimation
The problem we will solve is to estimate the decay-rate of a linear
system x' = Ax. This can be formulated as a generalized
eigenvalue problem.
Due to the product between t and
P, the problem cannot be solved
directly. However, it is easily solved by bisection in t.
Define the variables.
A = [-1 2;-3 -4];
P = sdpvar(2,2);
|
To find a lower bound on t, we solve a standard Lyapunov stability
problem.
F = set(P>eye(2))+set(A'*P+P*A < -eye(2));
solvesdp(F,trace(P));
P0 = double(P);
|
In the code above, we minimized the trace just to get a numerically sound
solution. This solution gives us a lower bound on decay-rate
t_lower = -max(eig(inv(P0)*(A'*P0+P0*A)))/2;
|
We now find an upper bound on the decay-rate by doubling t until the
problem is infeasible. To find out if the problem is infeasible, we check
the field problem in the solution structure. The meaning of
this variable is explained in the help text for the command yalmiperror. Infeasibility has been detected by the solver if the value is 1. To reduce
the amount of information written on the screen, we run the solver in a
completely silent mode. This can be accomplished by using the verbose
and warning options in
sdpsettings.
t_upper = t_lower*2;
F = set(P>eye(2))+set(A'*P+P*A < -2*t_upper*P);
ops = sdpsettings('verbose',0,'warning',0);
sol = solvesdp(F,[],ops);
while ~(sol.problem==1)
t_upper = t_upper*2;
F = set(P>eye(2))+set(A'*P+P*A < -2*t_upper*P);
sol = solvesdp(F,[],ops);
end
|
Having both an upper bound and a lower bound allows us to perform a
bisection.
tol = 0.01;
t_works = t_lower
while (t_upper-t_lower)>tol
t_test = (t_upper+t_lower)/2;
disp([t_lower t_upper t_test])
F = set(P>eye(2))+set(A'*P+P*A < -2*t_test*P);
sol = solvesdp(F,[],ops);
if sol.problem==1
t_upper = t_test;
else
t_lower = t_test;
t_works = t_test;
end
end
|
This example will be revisited later when we study
BMIs |