% Section 12.3 Aero-engine control, Part I
%  This file is used to produce the plots regarding 
%  Analysis of output selection problem for the aero 
%  engie.
% 
%  Kjetil Havre, 19/12-1995.
%

clear all
close all

% load the engine model
load aero1

w=logspace(-3,3,301);			

% scale the plant G_eng (Si and So are the input & output scaling matrices)
%G=mmult(So,G_eng,Si);
G = Gall;
Gjw = frsp(G,w);
Gj0 = frsp(G, 0);
[U0,S0,V0] = vsvd(Gj0);

 
Goh = round(1000*Gj0(1:6,1:3))/1000;
disp('Gall(0)=');
fprintf(1, '     %6.3f  %6.3f  %6.3f\n', vtp(Goh))
disp('RGA(Gall)=');
fprintf(1, '     %6.3f  %6.3f  %6.3f\n', vtp(rga(Goh)))
disp('RGA_sum=');
fprintf(1, '     %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n',sum(rga(Goh')))
U0h = round(1000*U0(1:6,1:6))/1000;
disp('U0=');
fprintf(1, '     %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n', vtp(U0h))

S0h = round(1000*S0(1:6,1:3))/1000;
disp('S0=');
fprintf(1, '     %6.3f  %6.3f  %6.3f\n', vtp(S0h))

V0h = round(1000*V0(1:3,1:3))/1000;
disp('V0=');
fprintf(1, '     %6.3f  %6.3f  %6.3f\n', vtp(V0h))

% Zeros.

Z   = szeros(G);
Z1  = szeros(sel(G,S1,':')); 
Z2  = szeros(sel(G,S2,':'));
Z3  = szeros(sel(G,S3,':'));
Z4  = szeros(sel(G,S4,':'));
Z5  = szeros(sel(G,S5,':'));
Z6  = szeros(sel(G,S6,':'));


% We have RHP zeros < 100 rad/s for the following sets:
%    In set 3, z = 30.853 ==> not use both NL and LPEMN measurements
%    In set 6, z = 27.664 ==> not use both OPR2 and LPEMN measurements
%    In set 7, z = 28.633 ==> not use both LPEMN and LPPR measurements
%
%    NL, LP-compressor spool speed is a function of LPPR the same is 
%    LPEMN ++> interactin effects which results in RHP zero
%
%    The OPR2 is a function of both LPPR and HPPR. LPEMN is a function of
%    LPPR ++> inteacting effects, which results in RHP zero.
%
%    Conclusin do not need to consider set s3 and s6 but we will consider set 7.
%clear Z3 Z6 s3 s6


disp(' ');
disp( '----------------------------------------');
disp( 'Frequency dependent relative gain array.');
disp( '----------------------------------------');

disp('Frequecy dependent RGA for full process.');
rga = vrga(Gjw);

RGA1j0 = vrga( sel(Gj0, [S1(2) S1(1) S1(3)], ':') ); 
RGA2j0 = vrga( sel(Gj0, [S2(2) S2(1) S2(3)], ':') );
RGA3j0 = vrga( sel(Gj0, [S3(2) S3(1) S3(3)], ':') );
RGA4j0 = vrga( sel(Gj0, S4, ':') );
RGA5j0 = vrga( sel(Gj0, S5, ':') );
RGA6j0 = vrga( sel(Gj0, S6, ':') );

RGA1 = vrga( sel(Gjw, [S1(2) S1(1) S1(3)], ':') );
RGA2 = vrga( sel(Gjw, [S2(2) S2(1) S2(3)], ':') );
RGA3 = vrga( sel(Gjw, [S3(2) S3(1) S3(3)], ':') );
RGA4 = vrga( sel(Gjw, S4, ':') );
RGA5 = vrga( sel(Gjw, S5, ':') );
RGA6 = vrga( sel(Gjw, S6, ':') );




rn1 = veval('sum', veval('sum', vabs(msub(RGA1,eye(3)))));
rn2 = veval('sum', veval('sum', vabs(msub(RGA2,eye(3)))));
rn3 = veval('sum', veval('sum', vabs(msub(RGA3,eye(3)))));
rn4 = veval('sum', veval('sum', vabs(msub(RGA4,eye(3)))));
rn5 = veval('sum', veval('sum', vabs(msub(RGA5,eye(3)))));
rn6 = veval('sum', veval('sum', vabs(msub(RGA6,eye(3)))));

%Figure 12.10
figure(1);
vplot('liv,d', rn1, rn2, rn3, rn4, '--', rn5,'-.', rn6)
axis([1e-3 1e2 0 7])
xlabel('Frequency')
text(6, 4.0,'S1');
text(10, 5.5,'S2');
text(4, 6.2,'S3');
text(20,  2.2,'S4');
text(20,  2.75,'S5');
text(31,  6.2,'S6');
drawnow

G0 = frsp( G, 0 );

% Smallest singular values:
Smsv1 = sel(vsvd(sel(G0,S1, ':')),3,1)
Smsv2 = sel(vsvd(sel(G0,S2, ':')),3,1)
Smsv3 = sel(vsvd(sel(G0,S3, ':')),3,1)
Smsv4 = sel(vsvd(sel(G0,S4, ':')),3,1)
Smsv5 = sel(vsvd(sel(G0,S5, ':')),3,1)
Smsv6 = sel(vsvd(sel(G0,S6, ':')),3,1)

clear all
load aero2
Gs = mmult(So, G_eng, Si );
P =    [0 0 0 1 0 0; 
        1 0 0 0 0 0;
        0 0 0 0 0 1;
        0 0 0 0 1 0;
        0 1 0 0 0 0;
        0 0 1 0 0 0];

S4 = [1 5 6];
S5 = [2 5 6];

Gall = mmult(P,Gs);
G4 = sel(Gall, S4, ':');
G5 = sel(Gall, S5, ':');

G1i = sel(G4, 1, ':'); % Output y1 = NL
G2i = sel(G5, 1, ':'); % Output y2 = OPR1
[G2ib, H4] = sysbal(G1i);
[G2ib, H5] = sysbal(G2i);

k = 1:1:15
%Figure 12.11
figure(2);
plot(k, H4, ':', k, H4, '+', k, H5, ':', k, H5, 'x' )
axis([0 16 0 0.4])
text(10, 0.3, 'Set 4: +');
text(10, 0.2, 'Set 5: x');
xlabel('n-state number');drawnow









