How to choose value of Ts for converting continuous state space (SS) matrix to discrete SS matrix?

Hi,

I converted Transfer function to state space matrix but I am curious how to choose best Ts?

I am not sure what is the function of this Ts related to plant?

I just used random values of Ts to check the affect of conversion, I found that changing Ts can chane the Discrete SS A,B matrices while it has no affect on C matrix.

I will use c2d distrete state space matrix in MPC controller for designing a controller to deploy on real plant. Please help me in understanding Ts and how to choose the best Ts value based on some facts?

My code is

%Transfer Function By System identification
num = [-1.324 6.6e05 -6.914e09 5.128e13 -1.054e19]
den = [1 1.99e04 1.371e09 1.145e13 4.037e17]
%Conversion to SS 

plant = tf([num],[den], 'OutputDelay',0.001)
% Simplify the model to remove insignificant poles and zeros

SSpzt = ss(plant) %This gives continuous time linear state space model of TFpzt
Ac = SSpzt.A;
Bc = SSpzt.B;
Cc = SSpzt.C;
Dc = SSpzt.D;

Ts = 0.0001 %Sampling time actual  = 3.1309e-05
Td = 0.009; %time Delay baed on frequency data

%Contineous Model to discrete Statespace Matrices
SSpzt_dis = c2dm(Ac,Bc,Cc,Dc,Ts)
       
[Ad,Bd,Cd,Dd]  = c2dm(Ac, Bc, Cc, Dc, Ts) 

sysdis = ss(Ad,Bd,Cd,Dd,Ts)
sysdis.D = 0

[Mathieu NOE]

I simply looked at the Bode plot and my first question is : how reliable is the TF model above 10^5 Hz ?

if the TF is supposed to have a constant 2nd order roll off then the asymptotic constant gain can be neglected

also I suppsoe that you don't want to simplify to the point the resonance at 3.6e3 Hz would disappear

all in all I suppose we consider the TF to be valid up to 10^5 Hz and that would be my lower limit for Fs, therefore Ts = 1e-5 s (max)

beside that , I don't understand why in your code Ts is so low Ts = 0.0001

Answer 1

[Tianyu]

Hi Ayesha,

As Mathieu mentioned, you can eye-ball the bode plot of the continuous time model and select the sampling time so that the system dynamics does not alter too much in the discrete time. A systematic way to address this is to use the bandwidth function. Below is the modified code.

Pay attention to the Select Ts part.

close all
%Transfer Function By System identification
num = [-1.324 6.6e05 -6.914e09 5.128e13 -1.054e19];
den = [1 1.99e04 1.371e09 1.145e13 4.037e17];
%Conversion to SS 
plant = tf([num],[den], 'OutputDelay',0.001);
% Simplify the model to remove insignificant poles and zeros
SSpzt = ss(plant); %This gives continuous time linear state space model of TFpzt
Ac = SSpzt.A;
Bc = SSpzt.B;
Cc = SSpzt.C;
Dc = SSpzt.D;

%% Select Ts
% w (rad/s) is the highest bandwidth given by the model.
% dbstop needs to be selected to cover the model dynamics. 
% For model given by Paul (in comment), we can select dbstop = -30.
dbdrop = -3; 
w = bandwidth(SSpzt,dbdrop); 
% converts unit to (/s), f is the Nyquist frequency.
f = w/(2*pi); 
% The scalar 2 before f comes from Nyquist theorem, one needs to sample 
% the system at least 2 times the Nyquist frequency to avoid aliasing. 
% factor is to ensure that aliasing won't happen. A typical choice 
% of the factor is 2~10. 
factor = 2; 
fs = 2*f*factor; 
Ts = 1/fs; % gives the sample time you want. 

%Contineous Model to discrete Statespace Matrices
SSpzt_dis = c2dm(Ac,Bc,Cc,Dc,Ts);
       
[Ad,Bd,Cd,Dd]  = c2dm(Ac, Bc, Cc, Dc, Ts) ;
sysdis = ss(Ad,Bd,Cd,Dd,Ts);
sysdis.OutputDelay = floor(0.001/Ts);

% Check the plot to observe that discrete time model has 
% similar bode and step plots as the continuous time model. 
subplot(2,1,1)
bode(SSpzt,sysdis)
subplot(2,1,2)
step(SSpzt,sysdis)

[Paul]

Hi Tianyu,

Only using bandwidth with the default dbdrop could be risky, in general, because the bandwidth might not cover all of the relevant dynamics. One really needs to understand all of the dynamics of the system.

For example, consider the following plant

h = tf(1,[1 1])*tf(1,[1/100^2 2*.02/100 1]);

And discretize with a factor of 10

w = bandwidth(h);
f = w/(2*pi);
factor = 10;
fs = 2*f*factor;
Ts = 1/fs; %Sampling time actual  = 3.1309e-05;
hd = c2d(h,Ts);

The phase response might not be (probabably isn't) acceptable.

bode(h,hd)

[Tianyu]

Hi Paul,

Good example! Thanks for pointing this out.

This highlights the importance of comparing the frequency response.