06. Discrete-Time Control, $z$-Plane Analysis & Discretization¶

Digital Control Systems Design and Analysis with ctrlpy

This tutorial demonstrates discrete-time LTI systems, difference equation simulations, unit-circle stability classification, continuous-to-discrete ($c2d$) conversion methods (ZOH, FOH, Tustin, Pre-warping, Matched), and complex $z$-plane pole-zero visualizations.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np

import ctrlpy as cp

# Ensure inline plotting
%matplotlib inline

import matplotlib.pyplot as plt import numpy as np import ctrlpy as cp # Ensure inline plotting %matplotlib inline

1. Defining Discrete-Time Transfer Functions¶

ctrlpy provides DiscreteTransferFunction (alias cp.dtf) for representing discrete LTI systems with sampling period $T_s > 0$.

Let's define a second-order digital filter/plant in standard variable $z$: $$H(z) = \frac{0.04837 z + 0.04526}{z^2 - 1.718 z + 0.7408}, \quad T_s = 0.1\text{ s}$$

In [ ]:

# Define discrete transfer function with sampling time Ts = 0.1 s
H = cp.dtf([0.04837, 0.04526], [1.0, -1.718, 0.7408], dt=0.1)

print("Discrete Transfer Function H(z):")
print(H)
print(f"Sampling Period Ts: {H.dt} s")
print(f"Discrete Poles: {H.poles()}")
print(f"Discrete Zeros: {H.zeros()}")
print(f"Discrete DC Gain H(1): {H.dcgain():.4f}")

# Define discrete transfer function with sampling time Ts = 0.1 s H = cp.dtf([0.04837, 0.04526], [1.0, -1.718, 0.7408], dt=0.1) print("Discrete Transfer Function H(z):") print(H) print(f"Sampling Period Ts: {H.dt} s") print(f"Discrete Poles: {H.poles()}") print(f"Discrete Zeros: {H.zeros()}") print(f"Discrete DC Gain H(1): {H.dcgain():.4f}")

In Jupyter, discrete systems natively render as formatted LaTeX mathematical expressions with explicit sampling period $T_s$:

In [ ]:

# Render formatted LaTeX representation
H

# Render formatted LaTeX representation H

Formulation in Delay Operator $z^{-1}$¶

Digital signal processing and microcontroller implementations often formulate difference equations in backward shift operator $z^{-1}$: $$H(z^{-1}) = \frac{0.04837 z^{-1} + 0.04526 z^{-2}}{1 - 1.718 z^{-1} + 0.7408 z^{-2}}$$

In [ ]:

H_delay = cp.dtf([0.0, 0.04837, 0.04526], [1.0, -1.718, 0.7408], dt=0.1, var="z^-1")
H_delay

H_delay = cp.dtf([0.0, 0.04837, 0.04526], [1.0, -1.718, 0.7408], dt=0.1, var="z^-1") H_delay

2. Discretization Methods (c2d)¶

Let's discretize a continuous underdamped second-order plant: $$G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad \omega_n = 5\text{ rad/s}, \; \zeta = 0.4$$

We compare all 5 continuous-to-discrete transformation methods supported by ctrlpy:

1. Zero-Order Hold (ZOH)
2. First-Order Hold (FOH)
3. Tustin (Bilinear Transform)
4. Tustin with Frequency Pre-Warping ($\omega_{\text{warp}} = \omega_n$)
5. Matched Pole-Zero Method

In [ ]:

wn = 5.0
zeta = 0.4
G_cont = cp.tf([wn**2], [1.0, 2.0 * zeta * wn, wn**2])
Ts = 0.08  # Sampling time

# Perform continuous-to-discrete conversions
H_zoh = cp.c2d(G_cont, dt=Ts, method="zoh")
H_foh = cp.c2d(G_cont, dt=Ts, method="foh")
H_tustin = cp.c2d(G_cont, dt=Ts, method="tustin")
H_prewarp = cp.c2d(G_cont, dt=Ts, method="tustin", prewarp_frequency=wn)
H_matched = cp.c2d(G_cont, dt=Ts, method="matched")

print("--- Discretized Transfer Functions ---")
print("ZOH:", H_zoh)
print("\nTustin:", H_tustin)
print("\nTustin Pre-warped:", H_prewarp)

wn = 5.0 zeta = 0.4 G_cont = cp.tf([wn**2], [1.0, 2.0 * zeta * wn, wn**2]) Ts = 0.08 # Sampling time # Perform continuous-to-discrete conversions H_zoh = cp.c2d(G_cont, dt=Ts, method="zoh") H_foh = cp.c2d(G_cont, dt=Ts, method="foh") H_tustin = cp.c2d(G_cont, dt=Ts, method="tustin") H_prewarp = cp.c2d(G_cont, dt=Ts, method="tustin", prewarp_frequency=wn) H_matched = cp.c2d(G_cont, dt=Ts, method="matched") print("--- Discretized Transfer Functions ---") print("ZOH:", H_zoh) print("\nTustin:", H_tustin) print("\nTustin Pre-warped:", H_prewarp)

Step Response Comparison: Continuous vs. Discretization Methods¶

Let's compare the continuous response $y(t)$ with the discrete step responses $y[k]$:

In [ ]:

# Simulate continuous response
t_cont = np.linspace(0.0, 3.0, 500)
res_cont = G_cont.step(T=t_cont)

# Simulate discrete responses
res_zoh = H_zoh.step(T=3.0)
res_foh = H_foh.step(T=3.0)
res_tustin = H_tustin.step(T=3.0)
res_prewarp = H_prewarp.step(T=3.0)
res_matched = H_matched.step(T=3.0)

plt.figure(figsize=(10, 6))
plt.plot(res_cont.t, res_cont.y, "k-", linewidth=2.0, label="Continuous G(s)")
plt.step(res_zoh.t, res_zoh.y, "b--", where="post", label="ZOH")
plt.step(res_foh.t, res_foh.y, "m-.", where="post", label="FOH")
plt.step(res_tustin.t, res_tustin.y, "g:", where="post", linewidth=2.0, label="Tustin")
plt.step(res_prewarp.t, res_prewarp.y, "r--", where="post", label="Tustin Pre-warped")
plt.step(res_matched.t, res_matched.y, "c-.", where="post", label="Matched PZ")

plt.title("Step Response Comparison: Continuous vs. Discretization Methods (Ts = 0.08 s)")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True, linestyle="--", alpha=0.6)
plt.legend(loc="lower right")
plt.tight_layout()
plt.show()

# Simulate continuous response t_cont = np.linspace(0.0, 3.0, 500) res_cont = G_cont.step(T=t_cont) # Simulate discrete responses res_zoh = H_zoh.step(T=3.0) res_foh = H_foh.step(T=3.0) res_tustin = H_tustin.step(T=3.0) res_prewarp = H_prewarp.step(T=3.0) res_matched = H_matched.step(T=3.0) plt.figure(figsize=(10, 6)) plt.plot(res_cont.t, res_cont.y, "k-", linewidth=2.0, label="Continuous G(s)") plt.step(res_zoh.t, res_zoh.y, "b--", where="post", label="ZOH") plt.step(res_foh.t, res_foh.y, "m-.", where="post", label="FOH") plt.step(res_tustin.t, res_tustin.y, "g:", where="post", linewidth=2.0, label="Tustin") plt.step(res_prewarp.t, res_prewarp.y, "r--", where="post", label="Tustin Pre-warped") plt.step(res_matched.t, res_matched.y, "c-.", where="post", label="Matched PZ") plt.title("Step Response Comparison: Continuous vs. Discretization Methods (Ts = 0.08 s)") plt.xlabel("Time [s]") plt.ylabel("Amplitude") plt.grid(True, linestyle="--", alpha=0.6) plt.legend(loc="lower right") plt.tight_layout() plt.show()

3. Complex $z$-Plane Pole-Zero Maps & Unit-Circle Stability¶

In the $z$-plane, stability is determined strictly by the Unit Circle $|z| = 1$:

• Inside the circle ($|z| < 1$): Stable (exponentially decaying modes).
• On the circle ($|z| = 1$): Marginally stable (non-decaying oscillations if simple).
• Outside the circle ($|z| > 1$): Unstable (exponentially diverging modes).

Let's render both static Matplotlib and interactive Plotly pole-zero maps:

In [ ]:

# Static Matplotlib Pole-Zero Map with Unit Circle
fig, ax = cp.plot_pzmap(H_zoh)
plt.show()

# Static Matplotlib Pole-Zero Map with Unit Circle fig, ax = cp.plot_pzmap(H_zoh) plt.show()

In [ ]:

# Interactive Plotly Pole-Zero Map with hover tooltips
fig_plotly = cp.iplot_pzmap(H_zoh)
fig_plotly.show()

# Interactive Plotly Pole-Zero Map with hover tooltips fig_plotly = cp.iplot_pzmap(H_zoh) fig_plotly.show()

4. Discrete Frequency Response (Bode Analysis)¶

In discrete systems, frequency response is defined over $\omega \in [0, \pi / T_s]$ up to the Nyquist folding frequency $\omega_N = \pi / T_s$:

In [ ]:

# Compute frequency response up to Nyquist limit
bdata_disc = H_zoh.bode(n_points=300)

fig, (ax_mag, ax_phase) = plt.subplots(2, 1, sharex=True, figsize=(8, 6))
ax_mag.semilogx(bdata_disc.w, bdata_disc.mag_db, "b-", linewidth=1.5)
ax_mag.set_ylabel("Magnitude (dB)")
ax_mag.set_title(
    f"Discrete Frequency Response (Nyquist Limit $\\omega_N = {np.pi / Ts:.2f}$ rad/s)"
)
ax_mag.grid(True, which="both", linestyle="--", alpha=0.6)

ax_phase.semilogx(bdata_disc.w, bdata_disc.phase, "b-", linewidth=1.5)
ax_phase.set_ylabel("Phase (deg)")
ax_phase.set_xlabel("Frequency (rad/s)")
ax_phase.grid(True, which="both", linestyle="--", alpha=0.6)

plt.tight_layout()
plt.show()

# Compute frequency response up to Nyquist limit bdata_disc = H_zoh.bode(n_points=300) fig, (ax_mag, ax_phase) = plt.subplots(2, 1, sharex=True, figsize=(8, 6)) ax_mag.semilogx(bdata_disc.w, bdata_disc.mag_db, "b-", linewidth=1.5) ax_mag.set_ylabel("Magnitude (dB)") ax_mag.set_title( f"Discrete Frequency Response (Nyquist Limit $\\omega_N = {np.pi / Ts:.2f}$ rad/s)" ) ax_mag.grid(True, which="both", linestyle="--", alpha=0.6) ax_phase.semilogx(bdata_disc.w, bdata_disc.phase, "b-", linewidth=1.5) ax_phase.set_ylabel("Phase (deg)") ax_phase.set_xlabel("Frequency (rad/s)") ax_phase.grid(True, which="both", linestyle="--", alpha=0.6) plt.tight_layout() plt.show()

5. Closed-Loop Digital Feedback Control Design¶

Let's design a discrete PID controller for our plant $G(s)$:

1. Define continuous PID controller $C(s) = K_p + \frac{K_i}{s} + \frac{K_d s}{1 + T_f s}$.
2. Discretize $C(s)$ using Tustin transform.
3. Form closed-loop digital feedback system $T(z) = \frac{C(z) G(z)}{1 + C(z) G(z)}$.
4. Simulate closed-loop response and extract transient specifications.

In [ ]:

# Continuous PID controller
Kp = 1.2
Ki = 2.0
Kd = 0.1
C_cont = cp.pid(Kp=Kp, Ki=Ki, Kd=Kd, Tf=0.01)

# Discretize controller via Tustin
C_disc = cp.c2d(C_cont, dt=Ts, method="tustin")

# Discretize plant via ZOH (physical D/A hold)
G_disc = cp.c2d(G_cont, dt=Ts, method="zoh")

# Form open-loop and closed-loop discrete systems
L_disc = C_disc * G_disc
T_closed = cp.feedback(L_disc, 1.0)

print("Closed-Loop Digital System T(z):")
print(T_closed)
print(f"Closed-Loop Stability: {T_closed.stability()}")
print(f"Closed-Loop Poles: {np.abs(T_closed.poles())}")

# Continuous PID controller Kp = 1.2 Ki = 2.0 Kd = 0.1 C_cont = cp.pid(Kp=Kp, Ki=Ki, Kd=Kd, Tf=0.01) # Discretize controller via Tustin C_disc = cp.c2d(C_cont, dt=Ts, method="tustin") # Discretize plant via ZOH (physical D/A hold) G_disc = cp.c2d(G_cont, dt=Ts, method="zoh") # Form open-loop and closed-loop discrete systems L_disc = C_disc * G_disc T_closed = cp.feedback(L_disc, 1.0) print("Closed-Loop Digital System T(z):") print(T_closed) print(f"Closed-Loop Stability: {T_closed.stability()}") print(f"Closed-Loop Poles: {np.abs(T_closed.poles())}")

In [ ]:

# Simulate closed-loop step response
res_cl = T_closed.step(T=3.0)

print("--- Digital Control Performance Metrics ---")
print(f"Steady-State Value: {res_cl.steady_state_value():.4f}")
print(f"Rise Time (10%-90%): {res_cl.rise_time():.4f} s")
print(f"Settling Time (2%):  {res_cl.settling_time(tolerance=0.02):.4f} s")
print(f"Percent Overshoot:   {res_cl.overshoot():.2f} %")
print(f"Peak Time:           {res_cl.peak_time():.4f} s")

# Plot closed-loop step response
fig, ax = T_closed.plot_step(T=3.0)
ax.set_title("Closed-Loop Digital PID Control Step Response")
plt.show()

# Simulate closed-loop step response res_cl = T_closed.step(T=3.0) print("--- Digital Control Performance Metrics ---") print(f"Steady-State Value: {res_cl.steady_state_value():.4f}") print(f"Rise Time (10%-90%): {res_cl.rise_time():.4f} s") print(f"Settling Time (2%): {res_cl.settling_time(tolerance=0.02):.4f} s") print(f"Percent Overshoot: {res_cl.overshoot():.2f} %") print(f"Peak Time: {res_cl.peak_time():.4f} s") # Plot closed-loop step response fig, ax = T_closed.plot_step(T=3.0) ax.set_title("Closed-Loop Digital PID Control Step Response") plt.show()

---

Summary¶

You have learned how to:

1. Construct and manipulate DiscreteTransferFunction models in $z$ and $z^{-1}$.
2. Perform continuous-to-discrete ($c2d$) transformations via ZOH, FOH, Tustin, Pre-warping, and Matched pole-zero methods.
3. Classify stability relative to the unit circle $|z| = 1$.
4. Simulate discrete time-domain step, impulse, and forced difference equations.
5. Design and verify closed-loop digital control systems with PID controllers.