%% Numerical Do's and Don'ts Demo -- High-Order Transfer Functions

%   Copyright 1986-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $  $Date: 2004/08/17 21:33:14 $

%% High-Order Transfer Functions
% When working with medium- to high-order systems (10+ states), avoid the
% transfer function form and use the state-space form whenever possible.
% Computations involving high-order transfer functions often result
% in severe loss of accuracy, and even overflow.  Even a simple
% product of two transfer functions can give surprising results, as shown
% below.
%
% Start with two discrete-time transfer functions Pd and Cd of order 9 and
% 2, respectively.  Their frequency responses are plotted above.

load numdemo               % load Pd,Cd models
bode(Pd,'b',Cd,'r')

h = gcr;
title('Bode diagram of Pd (blue) and Cd (red)');
set(h.AxesGrid,'XUnits','rad/sec','YUnits',{'dB';'deg'});
ax = getaxes(h.AxesGrid);
set(ax,'Xlim',[0.1 1e6]);
set(ax(1),'Ylim',[-100 100]);
set(ax(2),'Ylim',[-540 180]);

%% Pitfalls
% Next compute the product L = Pd*Cd (series connection) three different
% ways and compare the frequency responses: 

Ltf = Pd * Cd;                    % transfer function product
Lzp = zpk(Pd) * Cd;               % zero-pole-gain product
Lss = ss(Pd) * Cd;                % state space product
bode(Ltf,'g',Lzp,'b',Lss,'r-.',{1e-1,1e3})
h = gcr;
title('Bode diagram of Ltf (green), Lzp (blue), and Lss (red dash)');
set(h.AxesGrid,'XUnits','rad/sec','YUnits',{'dB';'deg'});
ax = getaxes(h.AxesGrid);
set(ax(1),'Ylim',[-20 120]);
set(ax(2),'Ylim',[-180 0]);

%% 
% The responses of the state-space and zero-pole-gain products match, but
% the response of the transfer function product Ltf is choppy and erratic
% below 100 rad/sec.
%          
% To explain the loss of accuracy with the transfer function form,
% compare the pole/zero maps of Pd and Cd near z=1:

pzmap(Pd,'b',Cd,'r')
h = gcr;
title('Pole/zero maps of Pd (blue) and Cd (red)');
h.AxesGrid.setxlim([0 1.05]);
h.AxesGrid.setylim([-1 1]);

%%
% Note the multiple roots near z=1.  Because the relative accuracy of
% polynomial values drops near roots, the relative error on the transfer
% function value near z=1 ends up exceeding 100%.  Observing that
% frequencies below 100 rad/s map to |z-1|<1e-3, this explains the erratic
% results below 100 rad/s.


%% Moral
% *Moral:*  Compute with the state-space form and avoid combining transfer
% functions.