function BurlCE52

% Computer exercise 5.2 in Burl's book

close all;

% Define system matrices
s=tf('s');
A = [-10 0 0 0; 0 -10 0 0; 6 -3 -5 0; 0 0.5 0 -5];
B = [10 0; 0 10; 0 0; 0 0];
C = [0 0 1 0; 0 0 0 1];
D = [0 0; 0 0];
K = (2 + 10/s) * eye(2);
% Check for nominal stability
% If G and K are stable and detIpGK does not have any zeros in the RHP, 
% then the system is internally stable (sufficient condition)
[numG1, denG1] = ss2tf(A, B, C, D, 1);
[numG2, denG2] = ss2tf(A, B, C, D, 2);
G = [numG1/denG1 numG2/denG2];
IpGK = [1 0; 0 1] + G*K;
detIpGK = IpGK(1,1)*IpGK(2,2) - IpGK(1,2)*IpGK(2,1);

% If Gall is stable, then then system is internally stable
% (necessary and sufficient condition)
Gall = [inv(eye(2)+K*G) -K*inv(eye(2)+G*K); G*inv(eye(2)+K*G) inv(eye(2)+G*K)];

% Check for robust stability. If the max SSV is less than 1, 
% then the system is robustly stable.
Acl = [-10 0 -20 0 10 0; 0 -10 0 -20 0 10; 6 -3 -5 0 0 0; 0 0.5 0 -5 0 0; 0 0 -10 0 0 0; 0 0 0 -10 0 0];
Bdelta = [0 0 0; 0 0 0; 1 1 0; 0 0 1; 0 0 0; 0 0 0];
Cdelta = [2 0 0 0 0 0; 0 1 0 0 0 0; 0 0.4 0 0 0 0];
Ddelta = [0 0 0; 0 0 0; 0 0 0];
Nydud = pck(Acl, Bdelta, Cdelta, Ddelta);
w = logspace(-2,2);
fr = frsp(Nydud, w);
blk = [1 1; 1 1; 1 1];
ssv = mu(fr, blk);
figure; loglog(ssv(:,3), ssv(:,1), ssv(:,3), ssv(:,2));
xlabel('frequency (rad/s)'); ylabel('SSV(Nydud)');

% Check for nominal performance. If the max principal gain is less than 0.2
% for frequencies < 0.1 rad/s, then the system meets the performance specs.
Bcl = [20 0; 0 20; 0 0; 0 0; 10 0; 0 10];
C = [0 0 1 0 0 0; 0 0 0 1 0 0];
D = [-1 0; 0 -1];
F = pck(Acl, Bcl, C, D);
fr = frsp(F, w);
sigma = vsvd(fr);
figure; vplot('liv,lm', sigma);
xlabel('frequency (rad/s)'); ylabel('principal gains (nominal system)');

% Check for robust performance. If the max SSV is less than 1,
% then the system performs robustly.
Ap = [Acl zeros(6,2); C [-0.1 0; 0 -0.1]];
Bp = [Bdelta Bcl; zeros(2,3) D];
Cp = [Cdelta zeros(3,2); zeros(2,6) [0.5 0; 0 0.5]];
Dp = zeros(5,5);
N = pck(Ap, Bp, Cp, Dp);
fr = frsp(N, w);
blk = [1 1; 1 1; 1 1; 2 2];
ssv = mu(fr, blk);
figure; vplot('liv,lm', ssv);
xlabel('frequency (rad/s)'); ylabel('SSV(N)');

% Alternate approach for robust performance test
Integrator = tf(1, [1 0]) * eye(6);
sigmaxinv = tf(0.5, [1 0.1]) * eye(2);
systemnames = 'Integrator Acl Bcl Bdelta C Cdelta D sigmaxinv';
inputvar = '[udel(3); w(2)]';
outputvar = '[Cdelta; sigmaxinv]';
input_to_Integrator = '[Acl+Bcl+Bdelta]';
input_to_Acl = '[Integrator]';
input_to_Bcl = '[w]';
input_to_Bdelta = '[udel]';
input_to_C = '[Integrator]';
input_to_Cdelta = '[Integrator]';
input_to_D = '[w]';
input_to_sigmaxinv = '[C+D]';
sysoutname = 'N';
sysic; % N = transfer function matrix
Nss = minreal(ss(N)); % Nss = state space model
Nf = frd(N, w); % frequency response data model of N
[SSV, info] = mussv(Nf, blk); % SSV bounds of N
muN = norm(SSV(1,1), inf); % upper bound of SSV of N