A Script for Symbolic and Numerical Solution of System Time Response for a Simple Closed Loop Discrete Time System

by Nasser Abbasi
May 11, 2005.

This is a simple Mathematica script which takes as input the reference signal R(s), the plant transfer function G(s), and the sample time T, and generate the sampled output c[kT]. The intermediate computation of G(z) and the closed loop transfer function T(z) and C(z) and c[KT] are all shown both in symbolic format and in numerical format.

This script solves only for this system configuration, for arbitrary G(s), T, and R(s).



To illustrate, the script is run to solve examples 6.1 and example 6.4 in the text book Digital Control System Analysis and Design, 3rd edition by Phillips and Nagle. The result generated matches the text book result with the additional advantage of being able to see the analytical steps taken during the solution steps.
Next I need to add the ability to generate the solution c(t) to the sampled input.

Data Input section

Here type in the input for the system (this could also be done via an Input statement if needed) for now, simply edit these lines to change the input as needed.

inputS = (1/s);    

plantS61 = 
      (4/((s + 2)));   
T61 = 
      0.1;  (*in seconds .  
          This is T for example 6.1*) 

plantS64 = 
      (1/(s (s + 1)));  (*this is G (s) for example 6.1 in text*) 
T64  = 
      1;  (*in seconds .  
          This is T for example 6.4*)

Example 6.1 output

process[T61, inputS, 
        plantS61];

Example 6.4 output

process[T64, inputS, 
      plantS64];

Code section

\(\(\(\ \)\(
\)\(\(Clear["Global`*"];\)

    \(normalize[expr_] := 
        Module[{d, c, l}, 
d = 
            Denominator[expr]; 
l = 
            CoefficientList[d, \ z]; 
c = 
            Coefficient[\ d\ , z, 
              Length[l] - 
                1]; 
\(\(Numerator[expr]\/c\)\(\ \)\)\/Expand[Denominator[expr]\/c]
];\)\)\)\)

\(\(process[sampleTime_, inputS_, plantS_] := 
        Module[{inputZ, zoh, plantZ, plantZ2, openLoopZ, closedLoopZ, c, cc, 
            tab}, 
Clear[T, k, 
            z]; 
inputZ = \ 
            ZTransform[
              InverseLaplaceTransform[inputS, s, t] /. t -> k\ T, k, 
              z]\ \ ; 
zoh = \ \(z - 1\)\/z; 
plantZ = \ 
            ZTransform[
              InverseLaplaceTransform[plantS, s, t] /. t -> k\ T, k, 
              z]\ \ ; 
plantZ2 = \ 
            ZTransform[
              InverseLaplaceTransform[plantS\/s, s, t] /. t -> k\ T, k, 
              z]\ \ ; 
openLoopZ = 
            zoh*\ plantZ2\  // Simplify\ ; 
closedLoopZ = 
            openLoopZ\/\(1 + openLoopZ\) // 
              Simplify; 
outputZ = 
            inputZ*\ closedLoopZ // Simplify; 
c = 
            Together[
              InverseZTransform[outputZ, z, 
                k]]\ ; 

Print["********** Symbolic solution is *************"]; 
Print["Input R(z) = ", 
            inputZ]; 
Print[\*""Plant (Subscript[G, p])(z) = "", normalize[
              plantZ]]; 
 (*\(Print[\*""modified Plant (Subscript[G, p2])(z) = "", 
                normalize[
                  plantZ2]];\)*) 
Print["G(z) = "\ , \ 
            normalize[
              openLoopZ]]; 
Print["Closed loop transfer function T(z) = ", 
            normalize[
              closedLoopZ]]; 
Print["Output C(z) = ", 
            normalize[
              outputZ]]; 
Print["Output c[k] = ", 
            c // FullSimplify]; 
Print[]; 
Print["********* Numerical solution is  *************"]; 
T = 
            sampleTime; 
Print["Input R(z) = ", 
            inputZ // 
              N]; 
Print[\*""Plant (Subscript[G, p])(z) = "", normalize[
              plantZ // 
                N]]; 
 (*\(Print[\*""modified Plant (Subscript[G, p2])(z) = "", 
                normalize[
                  plantZ2 // 
                    N]];\)*) 
Print["G(z) = "\ , \ 
            normalize[
              openLoopZ // 
                N]]; 
Print["Closed loop transfer function T(z) = ", 
            normalize[
              closedLoopZ // 
                N]]; 
Print["Output C(z) = ", 
            normalize[
              outputZ // 
                N]]; 
Print["Output c[k] = ", \(c // 
                N\) // FullSimplify]; 
Print[]; 

cc[i_] := 
            c /. k -> i; 
tab = 
            Table[{i, cc[i]}, {i, 0, 8}]; 
 (*\(Print[
                tab];\)*) 
p1 = 
            ListPlot[tab, PlotJoined -> True, 
              DisplayFunction -> Identity]; 
p2 = 
            ListPlot[tab, PlotStyle -> PointSize[0.02], 
              DisplayFunction -> Identity]; 
Show[{p1, 
              p2}, DisplayFunction -> $DisplayFunction, 
            AxesLabel -> {"Time (sec)", "c(t)"}, 
            PlotLabel -> "Discrete time response"];
\ ];\)\)