%Section 11.6.1 Reduction of a gas turbine aero-engine model
clear all
close all

% load the engine model
load aero0 				

minfo(G_eng);   % This is the 15 state 3 input 3 output model
				
% define frequency points	
%w1=logspace(-4,2,301);
%w2 = logspace(2, 3, 501);
%w3 = logspace(3, 4, 21);
%w = [w1 w2 w3];
w = logspace(-4, 4, 401);

G = G_eng;
Gjw = frsp(G, w);

Gd = strans(G);
Gdjw = frsp(Gd, w);
[Gb,hsig] = sysbal(G);
Gbjw = frsp(Gb, w);

%Figure 11.1
% a. Model reduction to 6 states using Rezidualization.
Gr = sresid(Gb, 6); 
Grjw = frsp(Gr, w); 
figure(1);subplot(231);
vplot('liv,lm', vsvd(Gjw), vsvd(Grjw), '--');
%set(gca, 'XTick', [1e-4 1e-2 1e0 1e2 1e4])
xlabel('Frequency'); 
axis([1e-2 1e4 1e-6 1e2]);drawnow


% b. Model reduction to 6 states using Truncation.
Gt = strunc(Gb, 6);
Gtjw = frsp(Gt, w);
subplot(232);
vplot('liv,lm', vsvd(Gjw), vsvd(Gtjw), '--');
xlabel('Frequency'); 
axis([1e-2 1e4 1e-6 1e2]);drawnow

% c. Model reduction using Optimal hankel norm approximation.
Gh = hankmr(Gb, hsig, 6, 'd');
Ghjw = frsp(Gh, w);
subplot(233);
vplot('liv,lm', vsvd(Gjw), vsvd(Ghjw), '--');
xlabel('Frequency'); 
axis([1e-2 1e4 1e-6 1e2]);drawnow

% Generate error systems.

Er = msub(G,Gr);
Et = msub(G,Gt);
Eh = msub(G,Gh);
Erjw = frsp(Er, w);
Etjw = frsp(Et, w);
Ehjw = frsp(Eh, w);

%Figure 11.2
figure(2);subplot(221);
vplot('liv,lm', vsvd(Erjw), ':', vsvd(Etjw), vsvd(Ehjw), '--')
xlabel('Frequency'); 
axis([1e-2 1e4 1e-10 1e1]);
drawnow

hinfnorm(Er, 1e-7)
hinfnorm(Et, 1e-7)
hinfnorm(Eh, 1e-7)
2*sum(hsig(7:15))

%Scaled truncation.
Gj0  = frsp(G, 0);
Gtj0 = frsp(Gt, 0);
h    = mmult(minv(Gtj0),Gj0);
St   = h(1:3, 1:3);
Gts = mmult(Gt, St);
Gtsjw = mmult(Gtjw, St);

%Scaled optimal Hankel-norm approximation.
Ghj0 = frsp(Gh, 0);
h    = mmult(minv(Ghj0),Gj0);
Sh   = h(1:3, 1:3);
Ghs = mmult(Gh, Sh);
Ghsjw = mmult(Ghjw, Sh);

%Error systems with scalings.

Ehsjw = msub(Gjw, Ghsjw);
Etsjw = msub(Gjw, Gtsjw);
subplot(222);
vplot('liv,lm', vsvd(Erjw), ':', vsvd(Etsjw), vsvd(Ehsjw), '--')
xlabel('Frequency'); 
axis([1e-2 1e4 1e-10 1e1]);
drawnow

Ehs = msub(G,Ghs);
Ets = msub(G,Gts);
hinfnorm(Ets, 1e-7)
hinfnorm(Ehs, 1e-7)

[A,B,C,D] = unpck(G);
t=0:0.0001:0.03;
[y,x,t] = impulse(A,B,C,D,2,t);
y = vpck(y,t);

% Impulse response for Balanced Residualization.
[A,B,C,D] = unpck(Gr);
[yr,x,t] = impulse(A,B,C,D,2,t);
yr = vpck(yr,t);
%Figure 11.3
figure(3);subplot(231)
vplot(y, yr, '--'); 
axis([0 0.03 -190 150])
xlabel('Time');
drawnow

% Impulse response for Scaled Balanced Truncation.
[A,B,C,D] = unpck(Gts);
[yt,x,t] = impulse(A,B,C,D,2,t);
yt = vpck(yt,t);

subplot(232);
vplot(y, yt, '--'); 
axis([0 0.03 -190 150])
xlabel('Time');
drawnow

% Impulse response for Scaled Hankel-norm approximation.
[A,B,C,D] = unpck(Ghs);
[yh,x,t] = impulse(A,B,C,D,2,t);
yh = vpck(yh,t);

subplot(233);
vplot(y, yh, '--'); 
axis([0 0.03 -190 150])
xlabel('Time');
drawnow

s = trsp(G, [0 1 0]', 1, 0.005);

% Step response for Balanced Residualization.
sr = trsp(Gr, [0 1 0]', 1, 0.005);
%Figure 11.4
figure(4);subplot(231);
vplot(s, sr, '--'); 
axis([0 1 -17.5 2.5]);
xlabel('Time');
drawnow

% Step response for Scaled Balanced Truncation.
st = trsp(Gts, [0 1 0]', 1, 0.005);
subplot(232);
vplot(s, st, '--'); 
axis([0 1 -17.5 2.5]);
xlabel('Time');
drawnow

% Step response for Scaled Optimal Hankel-norm Approximation.
sh = trsp(Ghs, [0 1 0]', 1, 0.005);
subplot(233);
vplot(s, sh, '--'); 
axis([0 1 -17.5 2.5]);
xlabel('Time');
drawnow












