%% Simplifying Higher-Order Plant Models
% This demo shows how you can approximate high-order plant models by simpler,
% low-order models using the tools available in the Robust Control Toolbox.

%% Motivation
% There are situations where you have to deal with models that are quite
% large:
%
% * *High-Order Plants:* Detailed first-principles or finite-element models
% of your plant tend to have high order. It is often desirable to simplify
% such models for simulation or control design purposes.
%
% * *High-Order Controllers:* Robust control techniques often yield
% high-order controllers. This is common, e.g., 
% when frequency weighting functions are used for shaping the open-loop
% response. It is often desirable to simplify such controllers for
% implementation purposes.
%
% For control purposes, it is generally enough to have an accurate model near
% the crossover frequency. For simulation purposes, it is enough to capture
% the essential dynamics in the frequency range of the excitation signals.
% This means that it is often possible to find low-order approximations of
% high-order models.  The Robust Control Toolbox offers a variety of model
% reduction algorithms to best suit your requirements and your model
% characteristics.

%% The Model Reduction Process
% A model reduction task typically involves the following steps:
% 
% * Analyze the important characteristics of the model from its time- or 
% frequency-domain responses using, e.g., |step| or |bode|
% * Determine an appropriate reduced order by plotting the  Hankel singular
% values of the original model (|hankelsv|) to determine which modes
% (states) can be discarded without sacrificing the key characteristics.
% * Choose a reduction algorithm. The various reduction methods
% available in the toolbox are:  
%    |balancmr|, |bstmr|, |schurmr|, |hankelmr|, and |ncfmr|.
% You can easily access these algorithms through the top-level interface
% |reduce|.  The methods employ different measures of "closeness"
% between the original and reduced models. The choice is application dependent.
% You may try all of them to investigate their relative merits. 
% * Validation: Validate your results by comparing the dynamics of the
% reduced model to the original. You may need to adjust your reduction
% parameters (choice of model order, algorithm, error bounds etc.) if the
% results are not satisfactory.  

%% Example: A Model for Rigid Body Motion of a Building 
% Let us see how to apply the reduction methods on a model of a building.
% The building (the Los Angeles University Hospital) has 8 floors, each
% considered as having 3 degrees of freedom - two displacements and one 
% rotation. We can represent the input-output relationship for any one of
% these displacements using a 48-state model, where each state represents
% a displacement or its rate of change (velocity).

load build.mat

%% Examining the Plant Dynamics
% Let us begin by analyzing the frequency response of the model.
bode(G)
grid on

%%
% As observed from the frequency response of the model, the essential
% dynamics of the system lie in the frequency range of 3 to 50 rad/s. The
% magnitude drops in both the very low as well has the high frequency
% ranges. Your objective is to find a 
% low-order model that preserves the information content in this
% frequency range to an acceptable level of accuracy.

%% Computing Hankel Singular Values
% In order to gain insight into what states of the model can be safely
% discarded, look at the Hankel singular values of the model.
%
hsv_add = hankelsv(G);
bar(hsv_add)
title('Hankel Singular Values of the Model (G)');
xlabel('Number of States')
ylabel('Singular Values (\sigma_i)')
line([10.5 10.5],[0 1.5e-3],'Color','r','linestyle','--','linewidth',1)
text(6,1.6e-3,'10 dominant states.')
%%
% The Hankel singular value plot suggests that there are 4
% dominant modes in this system. However, the contribution of the remaining
% modes is still significant. Let us draw the line at 10 states and discard
% the remaining ones. In other words, let's try to find a 10th-order reduced model
% |Gr| that best approximates the original system |G|.

%% Performing Model Reduction using an Additive Error Bound
% The function |reduce| is the gateway to all model reduction routines
% available in the toolbox. Let us use the default, Square-root Balance
% Truncation ('balancmr') option of |reduce| as the first step. This method
% uses an "additive" error bound for reduction, meaning that it tries to
% keep the absolute approximation error uniformly small across frequencies.

% Compute 10th-order reduced model (reduce uses balancmr method by default)
[Gr_add,info_add] = reduce(G,10);

% Now compare the original model G to the reduced model Gr_add
bode(G,'b',Gr_add,'r'), grid on
title('Comparing Original (G) to the Reduced model (Gr\_add)')
legend('G - 48-state original ','Gr\_add - 10-state reduced','location','northeast')

%% Performing Model Reduction using a Multiplicative Error Bound
% As seen from the Bode plot above, the reduced model captures the
% resonances below 30 rad/s quite well, but the match in the low
% frequency region (<2 rad/s) is poor. Also, the reduced model does not
% fully capture the dynamics in the 30-50 rad/s frequency range. A possible
% explanation for large errors at low frequencies is the relatively low
% gain of the model at these frequencies. Consequently, even large errors
% at these frequencies contribute little to the overall
% error. 
%
% To get around this problem, try a multiplicative-error method like |bstmr|.
% This algorithm emphasizes relative errors rather than absolute ones.
% Because relative comparisons do not work when gains are close to zero,
% you need to add a minimum-gain threshold, e.g., by adding a feedthrough 
% gain |d| to your original model.  Assuming you are not concerned about
% errors at gains below -100 dB, you can set the feedthrough to 1e-5.

GG = G;
GG.d = 1e-5;

%%
% Now look at the singular values for multiplicative (relative) errors
% (use 'mult' option of |hankelsv|)

hsv_mult = hankelsv(GG,'mult');
bar(hsv_mult)
title('Multiplicative-Error Singular Values of the Model (G)');
xlabel('Number of States')
ylabel('Singular Values (\sigma_i)')
%%
% A 26th-order model looks promising, but for sake of comparison with
% the previous result, we stick to a 10th order reduction. 

% Use bstmr algorithm option for model reduction 
[Gr_mult,info_mult] = reduce(GG,10,'algorithm','bst');

%now compare the original model G to the reduced model Gr_mult
bode(G,Gr_add,Gr_mult,{1e-2,1e4}), grid on
title('Comparing Original (G) to the Reduced models (Gr\_add and Gr\_mult)')
legend('G - 48-state original ','Gr\_add (balancmr)','Gr\_mult (bstmr)','location','northeast')

%%
% The fit between the original and the reduced model is 
% much better with the multiplicative-error approach, even at low frequencies. 
% This is confirmed by comparing the step responses:

step(G,Gr_add,Gr_mult,15) %step response until 15 seconds
legend('G: 48-state original ','Gr\_add: 10-state (balancmr)','Gr\_mult: 10-state (bstmr)')

%% Validating the Results 
% All algorithms provide bounds on the approximation error. For
% additive-error methods like |balancmr|, the approximation error
% is measured by the peak (maximum) gain of the error model
% |G-Greduced| across all frequencies. This peak gain is also known
% as the H-infinity norm of the error model. The error bound for 
% additive-error algorithms looks like
%
% $$\| G-Gr_{add} \|_\infty \leq  2 \sum_{i=11}^{48} \sigma_i := ErrorBound$$
% 
% where the sum is over all discarded Hankel singular values of |G| (entries 
% 11 through 48 of |hsv_add|). You can verify that this bound is satisfied 
% by comparing the two sides of this inequality:

norm(G-Gr_add,inf) % actual error

% theoretical bound (stored in the "ErrorBound" field of the "INFO" 
% struct returned by |reduce|)
info_add.ErrorBound 

%%
% For multiplicative-error methods like |bstmr|, the approximation error
% is measured by the peak gain across frequency of the relative error model
% |G\(G-Greduced)|, and the error bound looks like
%
% $$\|G^{-1}(G-Gr_{mult}) \|_\infty \leq  \prod_{i=11}^{48}(1+2\sigma_i(\sigma_i+\sqrt{1+\sigma^2_i}))-1 := ErrorBound$$
%
% where the sum is over the discarded *multiplicative* Hankel singular
% values computed by |hankelsv(G,'mult')|.  Again compare these bounds
% for the reduced model |Gr_mult|

norm(GG\(GG-Gr_mult),inf) % actual error

% Theoretical bound 
info_mult.ErrorBound 

%%
% Plot the relative error for confirmation
bodemag(GG\(GG-Gr_mult),{1e-2,1e3})
grid on
text(0.1,-50,'Peak Gain: -4.6 dB (59%) at 17.2 rad/s')
title('Relative error between original model (G) and reduced model (Gr\_mult)')

%%
% From the relative error plot above, there is up to 59% relative error 
% at 17.2 rad/s, which may be more than you are willing to accept.

%% Picking the Lowest Order Compatible with a Desired Accuracy Level    
% To improve the accuracy of |Gr_mult|, you need to increase the 
% order. If you want at least 5% relative accuracy, what is the 
% lowest order you can get? The function |reduce| can automatically
% select the lowest-order model compatible with your desired
% level of accuracy.

% Specify a maximum of 5% approximation error
[Gred,info] = reduce(GG,'ErrorType','mult','MaxError',0.05);
size(Gred)

%%
% The algorithm has picked a 34-state reduced model |Gred|.
% Compare the actual error with the theoretical bound:
norm(GG\(GG-Gred),inf)
info.ErrorBound

%%
% Look at the relative error magnitude as a function of frequency.
% Higher accuracy has been achieved at the expense of a larger model
% order (from 10 to 34). Note that the actual maximum error is 0.6%,
% much less than the 5% target. This discrepancy comes from the fact that
% |bstmr| uses the error bound rather than the actual error to select
% the order.

bodemag(GG\(GG-Gred),{1,1e3})
grid on
text(5,-75,'Peak Gain: -43.3 dB (0.6%) at 73.8 rad/s')
title('Relative error between original model (G) and reduced model (Gred)')

%%
% Compare the Bode responses
bode(G,Gred,{1e-2,1e4}), grid on
legend('G - 48-state original','Gred - 34-state reduced')

%%
% Finally, compare the step responses of the original and reduced models.
% They are virtually indistinguishable.

step(G,'b',Gred,'r--',15) %step response until 15 seconds
legend('G: 48-state original ','Gred: 34-state (bstmr)')
text(5,-4e-4,'Maximum relative error <0.05')

%% References
% The example discussed here was taken from the "Benchmark
% examples for model reduction of linear time invariant dynamical systems"
% available at: http://www.win.tue.nl/niconet/NIC2/benchmodred.html.