Convert Transfer Function from Discrete to Continuous Time

Due to the solid establishment of theories and methods, control system design is typically performed in continuous-time domain. But when it comes to implementation phase, nobody would want to construct an analog controller nowadays. So the continuous-time controller needs to be converted to one in the discrete-time domain. Scilab has mainly two commands for the task: dscr and cls2dls, which use zero-order hold and bilinear transformation methods, respectively. Both commands work with data in state-space form.

From time to time, a control engineer faces a problem in opposite direction. For example, she obtains a discrete-time plant from identification process and wants to convert it to continuous-time domain for analysis and design. Surprisingly enough, I could not find a Scilab command for such purpose. Maybe … I don’t know, it might appear in future version. But at this time I have to get around this problem by writing my own dls2cls function*. Bilinear transform method is used.

* Actually I prefer a shorter command name like d2c, but well, dls2cls makes it recognized as an exact inverse of cls2dls.

Before coding anything, I google to find if there exists a solution already. Could not find one other than some suggestion to use horner command in Scilab to convert polynomial of s to z using the Bilinear transform relationship. Conceptually, that should work. However, I tested it with a stable 5-th order transfer function and found that the result came out having higher order and, worse, an unstable continuous-time system. So this approach is not reliable due to some Scilab internal computational error beyond my comprehension.

I have no choice but to write my own algorithm from scratch. So here it goes.

Start from the Bilinear transform relationship

(1)   \begin{equation*}  s= \frac{2(z-1)}{T(z+1)}   \end{equation*}

where T is sampling period. From (1), we derive for z in terms of s as

(2)   \begin{equation*}  z= \frac{-(sT+2)}{(sT-2)}   \end{equation*}

To make the notation concise, define

(3)   \begin{equation*}  f(s) = -(sT+2) \end{equation*}

(4)   \begin{equation*}  g(s) = (sT-2) \end{equation*}

From now on we simply write them as f and g. So the relationship (2) becomes

(5)   \begin{equation*}  z = \frac{f}{g} \end{equation*}

To demonstrate the process, let us consider a third order discrete-time transfer function of the form

(6)   \begin{equation*}  P(z) = \frac{b_1+b_2z+b_3z^2+b_4z^3}{a_1+a_2z+a_3z^2+a_4z^3} \end{equation*}

Then its continuous-time counterpart can be found using (5) as

(7)   \begin{equation*}  P(s) = \frac{b_1+b_2\frac{f}{g}+b_3\frac{f^2}{g^2}+b_4\frac{f^3}{g^3}}{a_1+a_2\frac{f}{g}+a_3\frac{f^2}{g^2}+a_4\frac{f^3}{g^3}} \end{equation*}

Multiplying both numerator and denominator by g^3 yields

(8)   \begin{equation*}  P(s) = \frac{b_1g^3+b_2g^2f+b_3gf^2+b_4f^3}{a_1g^3+a_2g^2f+a_3gf^2+a_4f^3} \end{equation*}

This derivation can be extended to an n^{th} order discrete-time transfer function

(9)   \begin{equation*}  P(z) = \frac{b_1+b_2z+ \ldots+b_nz^{n-1}+b_{n+1}z^n}{a_1+a_2z+ \ldots+a_nz^{n-1}+a_{n+1}z^n} \end{equation*}

One can show that the corresponding continuous-time transfer function is

(10)   \begin{equation*}  P(s) = \frac{b_1g^n+b_2g^{n-1}f+\ldots+b_ngf^{n-1}+b_{n+1}f^n}{a_1g^n+a_2g^{n-1}f+\ldots+a_ngf^{n-1}+a_{n+1}f^n} \end{equation*}

The equation (10) can be easily written as Scilab algorithm. Indeed, this is a good exercise on writing a Scilab function. We show it step by step below.

Create a new file in Scinotes. Start with Scilab function placeholder.

function Sysc=dls2cls(Sysd, T)
 
endfunction

Save the file as dls2cls.sci in your chosen directory. The algorithm will sit between the function and endfunction keyword. The arguments passed to dls2cls are a discrete-time transfer function in z domain and sampling time T. The function returns a continuous-time transfer function Sysc.

At the beginning of function, define s and z variable

s=poly(0,'s');
z=poly(0,'z');

While dls2cls should work with any valid discrete-time system, we restrict it only to a proper (which includes strictly-proper) rational function; i.e., the denominator order must be greater or equal to the numerator order. This represents a system that is realizable. Use the following code to perform this check.

numdeg = degree(Sysd.num);
dendeg = degree(Sysd.den);
if numdeg>dendeg
    error('The discrete-time system is improper');
end

The error command displays the quoted message and exits the function in case an improper transfer function is passed to dls2cls.

To compute the algorithm (10), define f and g

f = -(s*T+2); 
g = (s*T-2);

and extract the coefficients from numerator and denominator

bn = coeff(Sysd.num);
an = coeff(Sysd.den);

Now, if the transfer function is strictly proper, the length of bn is shorter than an. Then we have to pad bn with 0 so that they have equal length.

npad = dendeg - numdeg;  // number of zero to be padded, if any
if npad > 0
    for i=0:npad
       bn = [bn 0]; // pad numerator coeff vector with zero 
    end 
end

The essence of this algorithm is the for loop that implements (10)

nums = 0;
dens = 0;
for i=0:dendeg,
    nums = nums + bn(i+1)*(f^i)*(g^(dendeg-i));
    dens = dens + an(i+1)*(f^i)*(g^(dendeg-i));
end

The computation is straightforward. After we derive the numerator and denominator of (10), a continuous-time transfer function is formed using syslin command

Sysc = syslin('c',nums,dens);

and returned to the caller.


Download the complete file –> dls2cls.sci


To test the function, load it into Scilab workspace

-->exec('dls2cls.sci',-1);

then create an arbitrary discrete-time model, say

-->numd = -2*z^4-z^3+z^2+2*z+3; 
-->dend = 3*z^4+2*z^3+z^2 -z -2; 
-->Sysd = syslin('d',numd/dend)
 Sysd  = 
              2   3    4  
    3 + 2z + z - z - 2z   
    -------------------   
             2    3    4  
  - 2 - z + z + 2z + 3z

Then pass to dls2cls with some specified sampling time T, say, 0.01

-->Sysc=dls2cls(Sysd,0.01)
 Sysc  =
                        2           3            4  
    48 - 2.08s + 0.0016s - 0.000028s + 1.000D-08s   
    ---------------------------------------------   
                        2           3            4  
    48 + 2.08s + 0.0016s + 0.000028s + 1.000D-08s

The function returns the continuous-time conversion. To verify the result, convert it back to discrete-time using Scilab cls2dls command. Note that this command requires a system with state-space representation.

-->Syscss=tf2ss(Sysc);
-->Sysdss = cls2dls(Syscss,0.01);

Then converts back to transfer function

-->Sysd1 = ss2tf(Sysdss)
 Sysd1  =
                               2            3            4  
    1 + 0.6666667z + 0.3333333z - 0.3333333z - 0.6666667z   
    -----------------------------------------------------   
                                       2            3   4   
  - 0.6666667 - 0.3333333z + 0.3333333z + 0.6666667z + z

We get the same discrete-time transfer function back (only normalized by dividing the numerator and denominator by 3).

I am happy with the result, though I have not yet tested this function extensively. If you find a bug kindly let me know.


Download the complete file –> dls2cls.sci


Comments

comments

Comments are closed.