Continuous vs Discrete linear system simulation Python

Question

Normally, when you discretize continuous system with zoh method, your continuous system output with the corresponding discrete inputs passed through zero-order hold should be the same as discrete system output at sampling instants. That is:

$sys = \frac{s + 1}{s^2+4s+5}$

And the corresponding discretized system with sampling time $T_s = 0.01$ and zero-order hold equivalence is:

$sysd = \frac{0.009851 z - 0.009753}{z^2 - 1.96 z + 0.9608}$.

Now, for a given $u(n) (n = 0, 1, ...)$ I simulate the discrete system inside for loop like:

$y(n+1) = 1.96 y(n) - 0.9608y(n-1)+ 0.009851u(n) - 0.009753u(n-1)$.

Then I want to find continuous system output to $u$ passed through zero order hold, meaning

$u(t) = u(n), for \ nT_s \leq t<(n+1)T_s$.

To this purpose, I decompose $u(t)$ to $u_n(t)$'s, where

$u_n(t) = U(t-nT_s)(u(n)-u(n-1))$,

and $U(t)$ is unit step function. Finally, for superposing the responses to all these $u_n(t)$'s, I find step resonse of continuous system for each $u_n(t)$, weight it by $(u(n) - u(n-1))$, shift it by $nT_s$ and add all the responses. That is,:

$sys(u(t)) = \sum_{0}^{N}{sys(u_i(t))} = \sum_{0}^{N} {sys(U(t-iT_s)(u(i) - u(i-1)))} = \sum_0^N{(u(i) - u(i-1))sys(U(t-iT_s))}$

Expectedly, both outputs (output of continuous system and discrete system with $u(n)$) should be the same at sampling instants, but here in my code it is different.

First of all, do you agree that the theoretical approach is correct? If yes, what is the problem in code, I also tried it in Matlab, but still different.

Just try these u values with the above system and its discrete version, with $T_s = 0.01$:

u =  [10.0, 9.0149, 8.156446201000001, 7.40885031569949, 6.758285462234175, 6.192642593060395, 5.701317169573098, 5.275022326677378, 4.905625237909304, 4.586003800165011, 4.309921114839202, 4.071915555493578, 3.867204486588504, 3.691599938148692, 3.5414347517271514, 3.4134978973869554, 3.3049778228833446, 3.213412837641881, 3.1366476579814817, 3.072795348504841, 3.020203989582603, 2.9774274840643464, 2.9431999892239666, 2.916413523772292, 2.8960983556695097, 2.8814058254283283, 2.871593302477761, 2.866011009711813, 2.864090484238467, 2.8653344711511473, 2.869308072374607, 2.875630994733954, 2.88397076074848, 2.894036762601707, 2.905575054584112, 2.918363792306603, 2.932209238369994, 2.9469422641488374, 2.9624152860827557, 2.9784995825184213, 2.995082943845555, 3.0120676145383865, 3.029368490853509, 3.046911542436347, 3.064632430030726, 3.0824752949388703, 3.100391698903132, 3.118339695729313, 3.136283018291052, 3.1541903665863655, 3.172034784296738, 3.189793112857524, 3.2074455134133073, 3.2249750482272233, 3.2423673141601883, 3.259610121752935, 3.276693214246804, 3.293608021582629, 3.3103474450330097, 3.3269056686628606, 3.343277994285584, 3.359460696996085, 3.375450898724326, 3.391246457570536, 3.406845870961228, 3.4222481909087, 3.437452949869936, 3.452460095887635, 3.4672699358596546, 3.4818830859264702, 3.4963004280917023, 3.510523072300681, 3.52455232329827, 3.5383896516714874, 3.552036668556261, 3.5654951035523497, 3.5787667854471006, 3.5918536253982856, 3.604757602269737, 3.6174807498515147, 3.630025145729702, 3.642392901600066, 3.6545861548454264, 3.666607061218924, 3.6784577884950043, 3.6901405109671046, 3.701657404686073, 3.7130106433465047, 3.7242023947397396, 3.73523481770234, 3.7461100594977537, 3.756830253576574, 3.7673975176676313, 3.777813952158069, 3.788081638725788, 3.798202639192171, 3.808178994567024, 3.818012724261144, 3.8277058254450025, 3.8372602725346994, 3.8466780167887107, 3.855960986000986, 3.865111084277786, 3.8741301918871853, 3.883020165171591, 3.8917828365148144, 3.900420014356283, 3.9089334832459333, 3.917325003934123, 3.9255963134916163]

And this is my code:

import numpy as np
import matplotlib.pyplot as plt
import control.matlab as cn

sys = cn.tf([1, 1], [1, 4, 5])

Ts = 0.01

sysd = cn.sample_system(sys, Ts, method = 'zoh')

u = np.array([10.0, 9.0149, 8.156446201000001, 7.40885031569949, 6.758285462234175, 6.192642593060395, 5.701317169573098, 5.275022326677378, 4.905625237909304, 4.586003800165011, 4.309921114839202, 4.071915555493578, 3.867204486588504, 3.691599938148692, 3.5414347517271514, 3.4134978973869554, 3.3049778228833446, 3.213412837641881, 3.1366476579814817, 3.072795348504841, 3.020203989582603, 2.9774274840643464, 2.9431999892239666, 2.916413523772292, 2.8960983556695097, 2.8814058254283283, 2.871593302477761, 2.866011009711813, 2.864090484238467, 2.8653344711511473, 2.869308072374607, 2.875630994733954, 2.88397076074848, 2.894036762601707, 2.905575054584112, 2.918363792306603, 2.932209238369994, 2.9469422641488374, 2.9624152860827557, 2.9784995825184213, 2.995082943845555, 3.0120676145383865, 3.029368490853509, 3.046911542436347, 3.064632430030726, 3.0824752949388703, 3.100391698903132, 3.118339695729313, 3.136283018291052, 3.1541903665863655, 3.172034784296738, 3.189793112857524, 3.2074455134133073, 3.2249750482272233, 3.2423673141601883, 3.259610121752935, 3.276693214246804, 3.293608021582629, 3.3103474450330097, 3.3269056686628606, 3.343277994285584, 3.359460696996085, 3.375450898724326, 3.391246457570536, 3.406845870961228, 3.4222481909087, 3.437452949869936, 3.452460095887635, 3.4672699358596546, 3.4818830859264702, 3.4963004280917023, 3.510523072300681, 3.52455232329827, 3.5383896516714874, 3.552036668556261, 3.5654951035523497, 3.5787667854471006, 3.5918536253982856, 3.604757602269737, 3.6174807498515147, 3.630025145729702, 3.642392901600066, 3.6545861548454264, 3.666607061218924, 3.6784577884950043, 3.6901405109671046, 3.701657404686073, 3.7130106433465047, 3.7242023947397396, 3.73523481770234, 3.7461100594977537, 3.756830253576574, 3.7673975176676313, 3.777813952158069, 3.788081638725788, 3.798202639192171, 3.808178994567024, 3.818012724261144, 3.8277058254450025, 3.8372602725346994, 3.8466780167887107, 3.855960986000986, 3.865111084277786, 3.8741301918871853, 3.883020165171591, 3.8917828365148144, 3.900420014356283, 3.9089334832459333, 3.917325003934123, 3.9255963134916163])

y = []
y.append(0.)
y.append(0.009851*u[0])

for i in np.arange(1, 110):   # 110 is arbitrary
    y.append(y[i]*1.96 - y[i-1]*0.9608 + u[i]*0.009851 - u[i-1]*0.009753)

T = np.arange(0, 0.8, Ts)
yc = cn.step(sys, T)[0]
yc = yc*(u[0])

for i in np.arange(1, T.size -5):
    yc[i:] = yc[i:] +  (u[i] - u[i-1]) * cn.step(sys, T[0:(T.size - i)])[0]

y = np.array(y)

plt.plot(np.arange(0, y.size * 0.01 - 0.0000001, 0.01), y)
plt.plot(T , yc)



plt.show()