%% Application of model order reduction techniques 
% Based on the example from the monograph by Obinata and Anderson (2001), page 54.

%% The original full-order model of the order 8
w(1) = 0.56806689746895; zeta(1) = 0.00096819582773; k(1) = 0.01651378989774; 
w(2) = 3.94093897440699; zeta(2) = 0.00100229920475; k(2) = 0.00257034576009;
w(3) = 10.58229653714164;zeta(3) = 0.00100167293203; k(3) = 0.00002188016252;
w(4) = 16.19234386986640;zeta(4) = 0.01000472824082; k(4) = 0.00027927762861;

G = tf(0);
for i=1:4
    G = G + k(i)* tf(w(i)^2,[1,2*zeta(i)*w(i),w(i)^2]);
end

%% Modal truncation
Gmt = tf(0);
for i=1:2
    Gmt = Gmt + k(i)* tf(w(i)^2,[1,2*zeta(i)*w(i),w(i)^2]);
end

%% Balanced realization
[Gbal,hankel_sv] = balreal(G); % pouziti funkce z Control System Toolboxu

%% Truncation of the balanced realization
Gbt = modred(Gbal,[5:8],'Trunc');

%% Residualization of the balanced realization
Gbr = modred(Gbal,[5:8],'MatchDC');

%% Balanced truncation using Schur decomposition 
Gst = schurmr(G,4);

%% Minimization of Hankel norm of the error
[Goh,redInfo_hankel] = hankelmr(G,4);

%% Minimizing multiplicative error
[Gbst,redInfo_bst] = bstmr(G,4);

%% Plotting Bode characteristics
figure(1)
bode(G,Gmt,Gbt,Gst,Gbr,Goh,Gbst), legend('G','G_{mt}','G_{bt}','G_{st}','G_{br}','G_{oh}','G_{bst}')
grid on, title('Bodeho diagram'), 

%% Another model, this time with different weights of different modes
% Here we give a different model. The difference is that the modes have
% different weights.

k(1) = 0.01; k(2) = 0.005; k(3) = 0.007; k(4) = 0.004;
G = tf(0);
for i=1:4
    G = G + k(i)* tf(w(i)^2,[1, 2*zeta(i)*w(i), w(i)^2]);
end

%% Modal truncation
Gmt = tf(0);
for i=1:2
    Gmt = Gmt + k(i)* tf(w(i)^2,[1, 2*zeta(i)*w(i), w(i)^2]);
end

%% Balanced realization
[Gbal,hankel_sv] = balreal(G); 
Gbt = modred(Gbal,[5:8],'Trunc');

%% Balanced residualization
Gbr = modred(Gbal,[5:8],'MatchDC');

%% Balanced truncation using Schur decomposition
[Gst,redinfo_schur] = schurmr(G,4);

%% Minimizing Hankel norm of the error of the approximation
[Goh,redInfo_hankel] = hankelmr(G,4);

%% Minimizing the multiplicative error
[Gbst,redInfo_bst] = bstmr(G,4);

%% Plotting the Bode characteristics
figure(2)
bode(G,Gmt,Gbt,Gst,Gbr,Goh,Gbst), legend('G','G_{mt}','G_{bt}','G_{st}','G_{br}','G_{oh}','G_{bst}')
grid on, title('Bodeho diagram'),

%% Plotting the errors
[Gmag,Gphase,Gw] = bode(G);

Emt = G-Gmt; 
Ebt = G-Gbt; 
Est = G-Gst; 
Ebr = G-Gbr; 
Eoh = G-Goh; 
Ebst = G-Gbst; 

figure(3)
bodemag(Emt,Ebt,Est,Ebr,Eoh,Ebst); grid on
legend('E_{mt}','E_{bt}','E_{st}','E_{br}','E_{oh}','E_{bst}')
title('Additive error'), xlabel('Frequency'), ylabel('Magnitude')

[Emtmag,Emtphase,Emtw] = bode(Emt,Gw); 
[Ebtmag,Ebtphase,Ebtw] = bode(Ebt,Gw);
[Estmag,Estphase,Estw] = bode(Est,Gw);
[Ebrmag,Ebrphase,Ebrw] = bode(Ebr,Gw);
[Eohmag,Eohphase,Eohw] = bode(Eoh,Gw);
[Ebstmag,Ebstphase,Ebstw] = bode(Ebst,Gw);

figure(4)
loglog(Emtw,squeeze(Emtmag./Gmag),Ebtw,squeeze(Ebtmag./Gmag),Estw,squeeze(Estmag./Gmag),...
    Ebrw,squeeze(Ebrmag./Gmag),Eohw,squeeze(Eohmag./Gmag),Ebstw,squeeze(Ebstmag./Gmag))
legend('R_{mt}','R_{bt}','R_{st}','R_{br}','R_{oh}','R_{bst}')
title('Relative error'), xlabel('Frequency'), ylabel('Magnitude'), grid on