ECE 481: Digital Control Systems
Yash Pant
Estimated study time: 1 hr 14 min
Table of contents
Sources and References
Primary references — C. L. Phillips, H. T. Nagle, and A. Chakrabortty, Digital Control System Analysis and Design, 4th ed., Pearson, 2015; G. F. Franklin, J. D. Powell, and A. F. Emami-Naeini, Feedback Control of Dynamic Systems, 8th ed., Pearson, 2019. Supplementary texts — K. J. Åström and B. Wittenmark, Computer-Controlled Systems: Theory and Design, 3rd ed., Prentice Hall, 1997; G. F. Franklin, J. D. Powell, and M. L. Workman, Digital Control of Dynamic Systems, 3rd ed., Addison-Wesley, 1998. Online resources — MIT OpenCourseWare 6.302 Feedback Systems; Karl Åström and Richard Murray, Feedback Systems: An Introduction for Scientists and Engineers (free PDF: fbsbook.org).
Chapter 1: Introduction to Digital and Sampled-Data Control
1.1 Why Digital Control?
Modern control systems are implemented on digital hardware: microcontrollers, FPGAs, digital signal processors (DSPs), and general-purpose CPUs running real-time operating systems. Rather than treating the digital implementation as a mere approximation of an analog design, digital control theory treats the discrete-time nature of computation as a first-class design consideration.
The key structural difference between an analog and a digital control loop is the presence of two conversion devices: an analog-to-digital converter (ADC) that samples continuous plant output at discrete time instants, and a digital-to-analog converter (DAC) that holds the computed control signal constant between updates. The digital processor computes the control law between sample instants. This architecture introduces phenomena—sampling, quantization, inter-sample behavior, and computational delay—that have no analog counterpart.
Other canonical examples include: robotics joint controllers (torque updates at 1–4 kHz), flight control computers (sensor fusion and surface actuation at 40–200 Hz), power electronics (PWM duty-cycle updates synchronized to the switching frequency), and network control systems where the “plant” and “controller” communicate over a shared communication channel with variable delay.
1.2 Structure of a Sampled-Data Control System
A sampled-data system contains both continuous-time and discrete-time components in the same feedback loop:
- Plant — a continuous-time dynamical system \(G_p(s)\).
- Sampler (ideal) — produces the sequence \(y[k] = y(kT)\) from the continuous output \(y(t)\).
- Digital controller — a difference-equation algorithm that maps \(\{e[k]\}\) to \(\{u[k]\}\).
- Zero-order hold (ZOH) — reconstructs a piecewise-constant continuous signal \(u(t)\) from \(\{u[k]\}\).
- Anti-aliasing filter — a continuous-time lowpass filter placed before the ADC to prevent frequency aliasing.
The analysis of this hybrid system requires matching the continuous and discrete domains. The approach taken in this course is to derive an exact discrete-time model of the continuous plant with ZOH actuation and then design the controller entirely in discrete time.
Chapter 2: Review of Signals, Systems, and Analog Control
2.1 BIBO and Internal Stability
A linear time-invariant (LTI) system is bounded-input bounded-output (BIBO) stable if every bounded input produces a bounded output. For a system with transfer function \(G(s)\), BIBO stability requires that all poles of \(G(s)\) lie in the open left half of the complex plane (LHP), i.e., \(\text{Re}(s_i) < 0\) for every pole \(s_i\).
Equivalently, all poles of \(G(s) = \mathcal{L}\{g(t)\}\) lie strictly in the open left half-plane.
Feedback stability is a stronger notion. For the unity-feedback system with open-loop transfer function \(L(s) = C(s)G_p(s)\), the closed-loop transfer function is
\[ T(s) = \frac{L(s)}{1 + L(s)}. \]The closed-loop poles are the roots of the characteristic equation \(1 + L(s) = 0\). Even if \(C(s)\) and \(G_p(s)\) individually have unstable pole-zero cancellations, those hidden modes may appear in the closed-loop system. Internal stability requires that all transfer functions in the four-transfer-function matrix of the closed loop have poles only in the open LHP.
2.2 Time-Domain Specifications
For a second-order system
\[ T(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \]the standard time-domain performance metrics for a unit step input are:
- Percent overshoot (PO): \(\text{PO} = 100 \exp\!\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right)\), valid for \(0 < \zeta < 1\).
- Peak time: \(t_p = \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}\).
- Settling time (2%): \(t_s \approx \frac{4}{\zeta\omega_n}\).
- Rise time: approximately \(t_r \approx \frac{1.8}{\omega_n}\) for \(0.3 \leq \zeta \leq 0.8\).
These metrics guide the placement of closed-loop poles in both s-domain and z-domain design.
2.3 Steady-State Performance and System Type
For the unity-feedback system, the steady-state error to polynomial reference inputs is characterized by the system type — the number of open-loop poles at the origin.
| Input | Type 0 | Type 1 | Type 2 |
|---|---|---|---|
| Step | \(\frac{1}{1+K_p}\) | 0 | 0 |
| Ramp | \(\infty\) | \(\frac{1}{K_v}\) | 0 |
| Parabola | \(\infty\) | \(\infty\) | \(\frac{1}{K_a}\) |
where the position, velocity, and acceleration constants are
\[ K_p = \lim_{s\to 0} L(s), \quad K_v = \lim_{s\to 0} s\,L(s), \quad K_a = \lim_{s\to 0} s^2 L(s). \]2.4 Frequency-Domain Analysis
The Nyquist stability criterion counts encirclements of the \(-1 + j0\) point in the Nyquist plot of \(L(j\omega)\) to determine closed-loop stability. For minimum-phase open-loop systems, the Bode stability criterion applies: the system is stable if the phase of \(L(j\omega)\) is greater than \(-180^\circ\) at the gain crossover frequency \(\omega_{gc}\) where \(\lvert L(j\omega_{gc})\rvert = 1\).
Gain margin (GM) and phase margin (PM) are the standard robustness measures:
\[ \text{GM} = \frac{1}{\lvert L(j\omega_{pc})\rvert}, \qquad \text{PM} = 180^\circ + \angle L(j\omega_{gc}), \]where \(\omega_{pc}\) is the phase crossover frequency where \(\angle L = -180^\circ\). Typical design targets are \(\text{GM} \geq 6\,\text{dB}\) and \(\text{PM} \geq 45^\circ\).
Chapter 3: Pole Placement for Continuous-Time Systems
3.1 State-Space Representation
A continuous-time LTI system is written in state-space form as
\[ \dot{x}(t) = Ax(t) + Bu(t), \qquad y(t) = Cx(t) + Du(t), \]where \(x \in \mathbb{R}^n\), \(u \in \mathbb{R}^m\), \(y \in \mathbb{R}^p\). The pair \((A,B)\) is controllable if and only if the controllability matrix
\[ \mathcal{C} = \begin{bmatrix} B & AB & A^2B & \cdots & A^{n-1}B \end{bmatrix} \]has rank \(n\).
3.2 Full-State Feedback and Pole Placement
With full state measurement, the control law \(u = -Kx\) places the closed-loop poles at the eigenvalues of \(A - BK\). If \((A,B)\) is controllable, then for any desired characteristic polynomial \(\alpha_d(s) = \prod_{i=1}^n (s - \lambda_i)\) with \(\lambda_i\) chosen freely, there exists a unique \(K\) achieving this.
where \(e_n^T = [0 \; 0 \;\cdots\; 0 \; 1]\) and \(\alpha_d(A)\) is the desired characteristic polynomial evaluated at the matrix \(A\).
The open-loop poles are at \(s = -1\) and \(s = -2\). Suppose we desire closed-loop poles at \(s = -4 \pm j4\), giving
\[ \alpha_d(s) = (s + 4 - j4)(s + 4 + j4) = s^2 + 8s + 32. \]The controllability matrix is \(\mathcal{C} = [B \; AB] = \begin{bmatrix} 0 & 1 \\ 1 & -3 \end{bmatrix}\), with \(\det(\mathcal{C}) = -1\), so the system is controllable.
Ackermann’s formula gives
\[ \alpha_d(A) = A^2 + 8A + 32I = \begin{bmatrix} -2 & -3 \\ 6 & 7 \end{bmatrix} + \begin{bmatrix} 0 & 8 \\ -16 & -24 \end{bmatrix} + \begin{bmatrix} 32 & 0 \\ 0 & 32 \end{bmatrix} = \begin{bmatrix} 30 & 5 \\ -10 & 15 \end{bmatrix}. \]Then \(K = e_2^T \mathcal{C}^{-1} \alpha_d(A) = [0\;1]\begin{bmatrix} 3 & -1 \\ -1 & 0\end{bmatrix}\begin{bmatrix}30 & 5\\-10 & 15\end{bmatrix}\).
Computing: \(\mathcal{C}^{-1} = \frac{1}{-1}\begin{bmatrix}-3 & -1 \\ -1 & 0\end{bmatrix} = \begin{bmatrix}3 & 1 \\ 1 & 0\end{bmatrix}\).
Row 2 of \(\mathcal{C}^{-1}\) is \([1\;0]\), so \(K = [1\;0]\begin{bmatrix}30 & 5\\-10 & 15\end{bmatrix} = [30\;5]\). One can verify \(\det(sI - A + BK) = s^2 + 8s + 32\). \(\square\)
3.3 Tracking with a Reference Input
Pure state feedback drives the output to zero. To track a nonzero reference \(r\), introduce a feedforward gain \(N\):
\[ u = -Kx + Nr. \]At steady state (\(\dot{x} = 0\)), we need \(y_{ss} = r\), which gives
\[ N = \frac{1}{C(BK - A)^{-1}B}. \]Alternatively, an integral action state (integrator on the tracking error \(e = r - y\)) can be appended to the plant state to achieve zero steady-state error robustly:
\[ \dot{x}_I = r - Cx, \qquad u = -Kx - K_I x_I. \]The augmented system has dimension \(n+1\), and the \(n+1\) poles are placed by Ackermann on the augmented controllability matrix.
Chapter 4: Discretization of Continuous-Time Controllers
4.1 The Discretization Problem
Given a continuous-time controller \(C(s)\) designed by analog methods (e.g., lead-lag, PID, pole placement), the goal is to find a difference equation \(C(z)\) that approximates \(C(s)\) in a computationally implementable form.
The key insight is that \(C(s)\) can be viewed as an integrator (or collection of integrators) in the time domain, and different numerical integration rules yield different \(s\)-to-\(z\) mappings.
4.2 Numerical Integration Methods
Let \(T\) be the sampling period. All methods substitute an approximation for \(s\) to obtain \(z\).
Forward Euler (Forward Difference)
The integral \(\int u\,dt\) is approximated by the left-endpoint rectangle rule:
\[ x[k+1] = x[k] + T \dot{x}[k] \implies \dot{x} \approx \frac{x[k+1] - x[k]}{T}. \]This gives the substitution
\[ s \;\longleftarrow\; \frac{z - 1}{T}. \]The stability region of the s-plane maps to a disk of radius \(T/2\) centered at \((1 + jT\cdot 0, T/2)\). For small \(T\), a stable \(C(s)\) may yield an unstable \(C(z)\), since the entire LHP does not map inside the unit circle.
Backward Euler (Backward Difference)
Using the right-endpoint rule:
\[ s \;\longleftarrow\; \frac{z - 1}{Tz}. \]The entire LHP maps inside the unit circle, so stability is always preserved. However, the frequency response is distorted at high frequencies.
Tustin / Bilinear Transformation
This is the trapezoidal integration rule:
\[ s \;\longleftarrow\; \frac{2}{T}\frac{z - 1}{z + 1}. \]The bilinear map \(z = \frac{1 + sT/2}{1 - sT/2}\) is a Möbius transformation that maps the imaginary axis of the s-plane exactly to the unit circle, and the open LHP to the interior of the unit circle. Therefore, a stable \(C(s)\) always yields a stable \(C(z)\).
The frequency-domain relationship is
\[ \omega_d = \frac{2}{T}\tan\!\left(\frac{\omega_a T}{2}\right), \]where \(\omega_d\) is the digital (analog) frequency corresponding to the analog frequency \(\omega_a\). This nonlinear “frequency warping” can be compensated by pre-warping: before applying the Tustin substitution, replace each design frequency \(\omega_0\) with the pre-warped frequency \(\hat{\omega}_0 = \frac{2}{T}\tan\!\left(\frac{\omega_0 T}{2}\right)\), so that the critical frequency (e.g., crossover or corner frequency) is exactly preserved.
4.3 Step-Invariant (ZOH Equivalent) Discretization
Rather than approximating the controller’s differential equation, the step-invariant method finds \(C(z)\) such that the discrete-time controller produces the same output samples as the continuous-time controller when both are driven by a piecewise-constant (ZOH) input:
\[ C(z) = (1 - z^{-1})\,\mathcal{Z}\!\left\{\frac{C(s)}{s}\right\}. \]This method preserves step-response matching exactly and is most appropriate when designing the controller to work with a ZOH plant. For the plant discretization it is the exact method (see Chapter 6).
- Forward Euler: simple, may destabilize; use only if T is very small relative to the fastest time constant.
- Backward Euler: always stable, but introduces a pole at origin; low-frequency behavior preserved.
- Tustin: best frequency response matching; use pre-warping when a specific frequency must be exactly reproduced.
- ZOH (step-invariant): exact step-response equivalence; preferred for plant discretization and when the design is done entirely in discrete time.
Chapter 5: Nonlinear Systems — Linearization and Compensation
5.1 State-Space Models of Nonlinear Systems
A general continuous-time nonlinear system is described by
\[ \dot{x}(t) = f(x(t), u(t)), \qquad y(t) = h(x(t), u(t)), \]where \(f : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n\) and \(h : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^p\) are smooth (continuously differentiable) functions. The tools of linear control theory do not apply directly, but in a neighborhood of a nominal operating point the system can be well approximated by a linear model.
5.2 Linearization
Let \((x_0, u_0)\) be an equilibrium point, i.e., \(f(x_0, u_0) = 0\). Define small perturbations \(\delta x = x - x_0\), \(\delta u = u - u_0\). Expanding \(f\) in a Taylor series and retaining only first-order terms:
\[ \dot{\delta x} \approx \left.\frac{\partial f}{\partial x}\right|_{(x_0,u_0)} \delta x + \left.\frac{\partial f}{\partial u}\right|_{(x_0,u_0)} \delta u =: A\,\delta x + B\,\delta u. \]Similarly,
\[ \delta y \approx \left.\frac{\partial h}{\partial x}\right|_{(x_0,u_0)} \delta x + \left.\frac{\partial h}{\partial u}\right|_{(x_0,u_0)} \delta u =: C\,\delta x + D\,\delta u. \]The Jacobian matrices \(A, B, C, D\) define the linearized system, which is LTI. This approximation is valid as long as the state remains near \((x_0, u_0)\).
5.3 Inverting Static Nonlinearities (Feedback Linearization)
When the nonlinearity is in the input path and is invertible, it can be canceled exactly. Suppose the plant model is
\[ \dot{x} = f(x) + g(x)\,\phi(v), \]where \(\phi(\cdot)\) is a known static nonlinearity and \(v\) is the new control input. Choosing
\[ v = \phi^{-1}(u), \]with \(u\) the output of a linear controller, cancels \(\phi\) and leaves a linear system for which standard linear design techniques apply.
5.4 Static Friction (Stiction)
Static friction (stiction) occurs when the actuator force is insufficient to overcome the breakaway friction force, causing the system to “stick” near a setpoint. This produces a limit cycle in closed-loop operation: the controller winds up integrator output until the breakaway force is overcome, causing an overshoot, after which the system re-sticks.
Compensation strategies include:
- Dither: a high-frequency low-amplitude signal added to the control input that keeps the mechanism on the verge of motion.
- Impulsive injection: a brief high-gain burst to break free, followed by immediate return to normal gain.
- Dead-zone compensation: feed forward a control increment that directly overcomes the stiction force.
Chapter 6: Linear Discrete-Time Systems
6.1 Difference Equations and the Shift Operator
A discrete-time LTI system is described by the constant-coefficient difference equation
\[ y[k] + a_1 y[k-1] + \cdots + a_n y[k-n] = b_0 u[k] + b_1 u[k-1] + \cdots + b_m u[k-m]. \]The z-transform converts this algebraic recursion into a rational transfer function. Define the forward shift operator \(z\): \(z\{x[k]\} = x[k+1]\), so that \(z^{-1}\{x[k]\} = x[k-1]\).
6.2 The z-Transform
where \(z \in \mathbb{C}\) and the series converges for \(\lvert z \rvert > R\) (the region of convergence, ROC).
Common z-transform pairs:
| Sequence \(x[k]\) | Z-transform \(X(z)\) | ROC |
|---|---|---|
| \(\delta[k]\) (unit impulse) | \(1\) | all \(z\) |
| \(u[k]\) (unit step) | \(\frac{z}{z-1}\) | \(\lvert z \rvert > 1\) |
| \(a^k u[k]\) | \(\frac{z}{z-a}\) | \(\lvert z \rvert > \lvert a \rvert\) |
| \(k\,a^k u[k]\) | \(\frac{az}{(z-a)^2}\) | \(\lvert z \rvert > \lvert a \rvert\) |
| \(\cos(\Omega_0 k)\,u[k]\) | \(\frac{z(z - \cos\Omega_0)}{z^2 - 2z\cos\Omega_0 + 1}\) | \(\lvert z \rvert > 1\) |
| \(\sin(\Omega_0 k)\,u[k]\) | \(\frac{z\sin\Omega_0}{z^2 - 2z\cos\Omega_0 + 1}\) | \(\lvert z \rvert > 1\) |
Key Properties
- Linearity: \(\mathcal{Z}\{\alpha x + \beta y\} = \alpha X(z) + \beta Y(z)\).
- Time delay: \(\mathcal{Z}\{x[k-n]\} = z^{-n}X(z)\) (for causal sequences with zero initial conditions).
- Time advance: \(\mathcal{Z}\{x[k+n]\} = z^n X(z) - z^n x[0] - \cdots - z x[n-1]\).
- Convolution: \(\mathcal{Z}\{x * h\} = X(z)H(z)\).
- Multiplication by \(a^k\): \(\mathcal{Z}\{a^k x[k]\} = X(z/a)\).
Final and Initial Value Theorems
6.3 Sampling and the Relationship Between s and z
The ideal sampler produces an impulse train:
\[ x^*(t) = \sum_{k=0}^{\infty} x(kT)\,\delta(t - kT). \]Taking the Laplace transform of \(x^*(t)\) and substituting \(z = e^{sT}\) gives the z-transform. This substitution,
\[ z = e^{sT}, \]maps the left half-plane to the interior of the unit circle: \(\text{Re}(s) < 0 \iff \lvert z \rvert < 1\). The imaginary axis \(\text{Re}(s) = 0\) maps to the unit circle \(\lvert z \rvert = 1\). Crucially, the mapping is many-to-one: \(s_0\) and \(s_0 + j\frac{2\pi k}{T}\) for integer \(k\) all map to the same \(z\). This is the origin of aliasing.
Nyquist-Shannon Sampling Theorem
where \(\omega_s = 2\pi/T\) is the sampling frequency. The minimum rate \(\omega_s = 2\omega_{max}\) is the Nyquist rate.
When sampling below the Nyquist rate, high-frequency components appear as low-frequency aliases—aliasing. The anti-aliasing (lowpass) filter before the ADC attenuates all components above \(\omega_s/2\).
In practice, sampling rates for control systems are chosen to be 10–30 times the closed-loop bandwidth to ensure inter-sample behavior is well approximated by the sampled model.
6.4 ZOH Equivalent Plant Model
Given a continuous-time plant \(G_p(s)\), the exact discrete-time equivalent model with ZOH actuation is
\[ G_p(z) = (1 - z^{-1})\,\mathcal{Z}\!\left\{\frac{G_p(s)}{s}\right\} = \frac{z-1}{z}\,\mathcal{Z}\!\left\{\frac{G_p(s)}{s}\right\}. \]Taking the z-transform:
\[ \mathcal{Z}\!\left\{\frac{G_p(s)}{s}\right\} = \frac{1}{a}\!\left(\frac{z}{z-1} - \frac{z}{z - e^{-aT}}\right) = \frac{z(1-e^{-aT})}{a(z-1)(z-e^{-aT})}. \]Multiplying by \((z-1)/z\):
\[ G_p(z) = \frac{1 - e^{-aT}}{a}\cdot\frac{1}{z - e^{-aT}}. \]The discrete pole is at \(z = e^{-aT}\), which is the exact image of the continuous pole at \(s = -a\) under \(z = e^{sT}\). \(\square\)
6.5 BIBO Stability of Discrete-Time Systems
Equivalently, all poles of the transfer function \(H(z)\) lie strictly inside the unit circle: \(\lvert z_i \rvert < 1\) for all poles \(z_i\).
Asymptotic stability of the autonomous system \(x[k+1] = Ax[k]\) requires all eigenvalues of \(A\) to lie strictly inside the unit circle.
6.6 Jury Stability Test
The Jury test is the discrete-time analog of the Routh-Hurwitz criterion. Given the characteristic polynomial
\[ p(z) = a_n z^n + a_{n-1}z^{n-1} + \cdots + a_1 z + a_0, \]construct the Jury table with rows built from the coefficients and their reversal. The system is stable if and only if all leading minors (the elements in the first column of the table after a specific construction) are positive. For the first few conditions, a quick check is:
- \(p(1) > 0\)
- \((-1)^n p(-1) > 0\)
- \(\lvert a_0 \rvert < a_n\)
- Proceed with the full Jury array for \(n > 2\).
- \(p(1) = 1 + b + c > 0\)
- \(p(-1) = 1 - b + c > 0\)
- \(\lvert c \rvert < 1\)
6.7 Stability via Bilinear Transform and Routh-Hurwitz
An alternative to the Jury test is to apply the bilinear (Tustin) substitution
\[ z = \frac{1 + w}{1 - w} \quad \Longleftrightarrow \quad w = \frac{z - 1}{z + 1} \]to the characteristic polynomial \(p(z)\), obtaining a polynomial \(\tilde{p}(w)\) in \(w\). The mapping \(z = (1+w)/(1-w)\) sends the interior of the unit circle to the open left half of the \(w\)-plane. Applying the Routh-Hurwitz criterion to \(\tilde{p}(w)\) then determines stability.
6.8 Frequency Response of Discrete-Time Systems
For a stable discrete-time system \(H(z)\), the frequency response is obtained by evaluating \(H(z)\) on the unit circle:
\[ H(e^{j\Omega}) = H(z)\Big|_{z = e^{j\Omega}}, \]where \(\Omega \in [-\pi, \pi]\) is the digital frequency (radians per sample). The relationship to physical frequency \(\omega\) (rad/s) is \(\Omega = \omega T\).
The magnitude and phase plots of \(H(e^{j\Omega})\) constitute the digital Bode plot. Note that \(H(e^{j\Omega})\) is periodic with period \(2\pi\), so only \(\Omega \in [0, \pi]\) (corresponding to \(f \in [0, f_s/2]\)) is physically meaningful.
Gain and Phase Margins in Discrete Time
For the discrete-time unity-feedback loop with open-loop transfer function \(L(z) = C(z)G_p(z)\), the gain margin and phase margin are defined identically to the continuous-time case but evaluated on the unit circle:
\[ \text{GM} = \frac{1}{\lvert L(e^{j\Omega_{pc}})\rvert}, \quad \text{PM} = 180^\circ + \angle L(e^{j\Omega_{gc}}), \]where \(\Omega_{pc}\) satisfies \(\angle L(e^{j\Omega_{pc}}) = -180^\circ\) and \(\Omega_{gc}\) satisfies \(\lvert L(e^{j\Omega_{gc}})\rvert = 1\).
Chapter 7: Control Design in Discrete Time
7.1 Design Workflow
The two main strategies for discrete-time controller design are:
Design-then-discretize: Design \(C(s)\) in continuous time, then discretize to \(C(z)\) using Tustin, ZOH, or another method. Adequate when \(T\) is small relative to the closed-loop bandwidth.
Discretize-then-design (direct digital design): Obtain the ZOH-equivalent plant \(G_p(z)\), then design \(C(z)\) directly using discrete-time frequency response or state-space methods. This is the preferred approach for slower sampling rates.
7.2 Discrete-Time State-Space
A discrete-time state-space model is
\[ x[k+1] = A_d x[k] + B_d u[k], \qquad y[k] = C_d x[k] + D_d u[k]. \]For the ZOH-equivalent of \(\dot{x} = Ax + Bu\), \(y = Cx + Du\):
\[ A_d = e^{AT}, \qquad B_d = \left(\int_0^T e^{A\tau}\,d\tau\right) B = A^{-1}(e^{AT} - I)B \quad (\text{if }A \text{ is invertible}). \]The matrix exponential is computed via the series \(e^{AT} = I + AT + \frac{(AT)^2}{2!} + \cdots\).
Controllability and Observability in Discrete Time
The controllability matrix is identical in form:
\[ \mathcal{C}_d = \begin{bmatrix} B_d & A_d B_d & A_d^2 B_d & \cdots & A_d^{n-1}B_d \end{bmatrix}. \]The pair \((A_d, B_d)\) is controllable if and only if \(\text{rank}(\mathcal{C}_d) = n\). Controllability of the continuous pair \((A, B)\) implies controllability of \((A_d, B_d)\) except at sampling frequencies that are multiples of eigenvalue frequencies (sampling zeros can create issues).
The observability matrix is
\[ \mathcal{O}_d = \begin{bmatrix} C_d \\ C_d A_d \\ \vdots \\ C_d A_d^{n-1} \end{bmatrix}. \]7.3 Pole Placement in Discrete Time
The same Ackermann formula applies for \(A_d\), \(B_d\), and a desired characteristic polynomial \(\alpha_d(z)\) with roots inside the unit circle. The dominant pole placement guidelines for a second-order discrete-time system with natural frequency \(\omega_n\) and damping \(\zeta\) are:
\[ z_{1,2} = e^{(-\zeta \omega_n \pm j\omega_n\sqrt{1-\zeta^2})T}. \]The performance specifications translate from s-domain \(\lambda_i\) to z-domain via \(z_i = e^{\lambda_i T}\).
The ZOH matrices are
\[ A_d = e^{AT} = I + AT = \begin{bmatrix}1 & T \\ 0 & 1\end{bmatrix} = \begin{bmatrix}1 & 0.1\\0 & 1\end{bmatrix}, \quad B_d = \begin{bmatrix}T^2/2\\T\end{bmatrix} = \begin{bmatrix}0.005\\0.1\end{bmatrix}. \]Desired continuous poles at \(s = -5 \pm j5\): \(z_{1,2} = e^{(-5 \pm j5)(0.1)} = e^{-0.5}e^{\pm j0.5} \approx 0.606\angle\pm 28.6^\circ\).
\[ \alpha_d(z) = (z - 0.606 e^{j0.5})(z - 0.606 e^{-j0.5}) = z^2 - 2(0.606)\cos(0.5)\,z + 0.606^2. \]\(\cos(0.5) \approx 0.8776\), so \(\alpha_d(z) = z^2 - 1.0638\,z + 0.3679\).
The controllability matrix \(\mathcal{C}_d = [B_d\;A_dB_d] = \begin{bmatrix}0.005&0.105\\0.1&0.1\end{bmatrix}\), \(\det = -0.0100\). Applying Ackermann’s formula yields the gain vector \(K\) that places the poles at the desired locations.
7.4 Digital PID Controller
The continuous PID controller is
\[ C(s) = K_p + \frac{K_i}{s} + K_d s. \]Applying the Tustin substitution \(s \to \frac{2}{T}\frac{z-1}{z+1}\) gives the digital PID:
\[ C(z) = K_p + K_i \frac{T}{2}\frac{z+1}{z-1} + K_d \frac{2}{T}\frac{z-1}{z+1}. \]In the position form, the output at sample \(k\) is
\[ u[k] = K_p e[k] + K_i \frac{T}{2}\sum_{j=0}^{k}(e[j] + e[j-1]) + K_d \frac{2}{T}(e[k] - e[k-1]), \]where \(e[k] = r[k] - y[k]\). The velocity (incremental) form avoids integrator windup issues:
\[ \Delta u[k] = u[k] - u[k-1] = K_p (e[k] - e[k-1]) + K_i T\,e[k] + K_d \frac{e[k] - 2e[k-1] + e[k-2]}{T}. \]Anti-windup: When the actuator saturates, the integrator should be paused or back-calculated. The back-calculation method feeds the difference between the unsaturated and saturated commands back to the integrator input through a gain \(1/T_t\):
\[ \dot{x}_I = e + \frac{u_{sat} - u}{T_t}, \]where \(T_t\) (the tracking time constant) is typically chosen as \(T_t = \sqrt{T_i T_d}\).
7.5 Lead-Lag Compensation in the z-Domain
Continuous lead and lag compensators have the form
\[ C_{\text{lead}}(s) = K\frac{s + z_c}{s + p_c}, \quad z_c < p_c \quad (\text{lead: zero inside pole}), \]\[ C_{\text{lag}}(s) = K\frac{s + z_c}{s + p_c}, \quad z_c > p_c \quad (\text{lag: zero outside pole}). \]The corresponding discrete-time compensators obtained by Tustin are
\[ C_{\text{lead/lag}}(z) = K' \frac{z - z_d}{z - p_d}, \]where \(z_d\) and \(p_d\) are the Tustin-mapped zeros and poles. Lead compensation adds phase near the crossover frequency, improving phase margin. Lag compensation boosts low-frequency gain to reduce steady-state error with minimal phase degradation at crossover.
7.6 Root Locus in the z-Domain
The root locus for a discrete-time system is constructed with exactly the same rules as in continuous time, but the stability boundary is now the unit circle rather than the imaginary axis. The rules (angles of departure/arrival, breakaway points, asymptote angles) all apply unchanged; the design goal is to keep locus branches inside the unit circle.
Key guidelines:
- Desired closed-loop pole locations in the z-plane correspond to \(z = e^{\lambda T}\), where \(\lambda = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2}\).
- Constant-damping-ratio lines in the z-plane are logarithmic spirals.
- Constant-natural-frequency circles in the z-plane are circles of radius \(e^{-\zeta\omega_nT}\).
- Lines of constant \(\omega_d = \omega_n\sqrt{1-\zeta^2}\) appear as radial lines from the origin at angle \(\pm\omega_d T\).
Deadbeat Control
A special pole placement choice in discrete time is deadbeat control: all closed-loop poles placed at \(z = 0\). A deadbeat controller drives any initial condition to zero in exactly \(n\) steps (for an \(n\)-th order system). While this achieves the fastest possible settling, it typically requires large control effort and produces poor inter-sample behavior.
7.7 Frequency-Domain Design: Loop Shaping in Discrete Time
The discrete-time open-loop transfer function is evaluated on the unit circle:
\[ L(e^{j\Omega}) = C(e^{j\Omega}) G_p(e^{j\Omega}). \]The design objective is to shape \(\lvert L(e^{j\Omega})\rvert\) and \(\angle L(e^{j\Omega})\) to achieve:
- High gain at low frequencies for reference tracking and disturbance rejection.
- Gain crossover at the desired bandwidth frequency \(\Omega_c = \omega_c T\).
- Adequate phase margin (typically \(\geq 45^\circ\)) at crossover.
- Gain rolloff at high frequencies for noise attenuation and robustness to unmodeled dynamics.
The discretization of \(L(s)\) introduces additional phase lag at high frequencies relative to the analog model. This computation delay effect is especially significant: a computational delay of one full sample period \(T\) corresponds to adding \(-\Omega\) radians of phase at digital frequency \(\Omega\), i.e., a phase shift of \(-\omega T\) rad at analog frequency \(\omega\). This must be accounted for in the margin calculations by modeling the delay as \(z^{-1}\) in the loop.
7.8 Observer Design: Luenberger Observer
When not all states are measurable, an observer (state estimator) reconstructs the state from input-output measurements:
\[ \hat{x}[k+1] = A_d \hat{x}[k] + B_d u[k] + L(y[k] - C_d \hat{x}[k]). \]The observer error \(e_o[k] = x[k] - \hat{x}[k]\) obeys
\[ e_o[k+1] = (A_d - LC_d)\,e_o[k]. \]The observer gain \(L\) is chosen so that the eigenvalues of \(A_d - LC_d\) lie inside the unit circle (for stability) and converge significantly faster than the closed-loop poles (typically 2–5 times faster). The separation principle holds for linear systems: the observer poles and controller poles can be placed independently, and the combined observer-based controller achieves the desired closed-loop performance.
The observer gain is designed by duality: the problem is equivalent to pole placement for the pair \((A_d^T, C_d^T)\), so Ackermann’s formula gives
\[ L = \alpha_o(A_d)\,\mathcal{O}_d^{-T}\,e_n, \]where \(\alpha_o(z)\) is the desired observer characteristic polynomial.
Chapter 8: Optimization-Based Control of Constrained Systems
8.1 Introduction to Optimal Control
Standard pole placement and loop-shaping design does not directly account for constraints on states and inputs, nor does it provide a systematic means to trade off performance and control effort. Optimal control formulates the controller design as an optimization problem.
The Linear Quadratic Regulator (LQR) minimizes the infinite-horizon cost:
\[ J = \sum_{k=0}^{\infty} \left( x[k]^T Q\, x[k] + u[k]^T R\, u[k] \right), \]where \(Q = Q^T \succeq 0\) and \(R = R^T \succ 0\) are design weight matrices. The optimal control law is linear state feedback:
\[ u^*[k] = -K_{\text{LQR}}\, x[k], \qquad K_{\text{LQR}} = (R + B_d^T P B_d)^{-1} B_d^T P A_d, \]where \(P = P^T \succ 0\) is the unique solution to the discrete-time algebraic Riccati equation (DARE):
\[ P = Q + A_d^T P A_d - A_d^T P B_d (R + B_d^T P B_d)^{-1} B_d^T P A_d. \]LQR provides guaranteed stability margins (\(\text{GM} \geq 6\,\text{dB}\), \(\text{PM} \geq 60^\circ\) for the continuous-time case) but does not handle hard constraints.
8.2 Effect of Constraints
Physical systems always have constraints:
- Input constraints: actuator saturation, \(u_{\min} \leq u[k] \leq u_{\max}\).
- State constraints: safety bounds, e.g., position limits, temperature limits.
When an LQR controller encounters an active constraint, it cannot achieve the optimal unconstrained behavior, and the closed-loop system may lose its stability guarantees. Constraint handling requires solving an optimization problem at each control step rather than evaluating a fixed gain matrix.
8.3 Model Predictive Control (MPC)
Model Predictive Control (also known as receding horizon control) is the dominant framework for constrained optimal control in industry. At each time step \(k\), MPC solves a finite-horizon optimal control problem:
\[ \min_{\mathbf{u}} \sum_{i=0}^{N-1} \left[ x[k+i]^T Q\, x[k+i] + u[k+i]^T R\, u[k+i] \right] + x[k+N]^T P_f\, x[k+N] \]subject to:
\[ x[k+i+1] = A_d x[k+i] + B_d u[k+i], \quad i = 0,\ldots, N-1, \]\[ u_{\min} \leq u[k+i] \leq u_{\max}, \quad i = 0,\ldots, N-1, \]\[ x_{\min} \leq x[k+i] \leq x_{\max}, \quad i = 1,\ldots, N, \]\[ x[k] = x_{\text{current}}. \]Here \(N\) is the prediction horizon, \(P_f\) is a terminal cost weight (often chosen as the LQR cost-to-go matrix), and \(x[k]\) is the measured (or estimated) current state.
Receding Horizon Principle
Only the first element of the optimal control sequence is applied:
\[ u^*[k] = u^*_{0|k}. \]At the next time step \(k+1\), the optimization is re-solved with the new measured state. This receding horizon (rolling horizon) strategy provides feedback, automatically compensates for model errors and disturbances, and keeps the problem horizon always \(N\) steps long.
8.4 MPC as a Quadratic Program
For linear systems with linear constraints and a quadratic cost, the MPC optimization is a Quadratic Program (QP). Define the stacked input vector
\[ \mathbf{U} = \begin{bmatrix} u[k] \\ u[k+1] \\ \vdots \\ u[k+N-1] \end{bmatrix} \in \mathbb{R}^{mN}. \]The predicted state trajectory satisfies \(\mathbf{X} = \Phi x[k] + \Gamma \mathbf{U}\) for appropriately defined matrices \(\Phi\) and \(\Gamma\). Substituting into the cost yields a QP in \(\mathbf{U}\):
\[ \min_{\mathbf{U}} \;\frac{1}{2}\mathbf{U}^T H\,\mathbf{U} + x[k]^T F\,\mathbf{U} \]subject to linear inequality constraints \(A_{\text{ineq}}\mathbf{U} \leq b_{\text{ineq}}(x[k])\). The QP can be solved in real time using active-set or interior-point methods. For small prediction horizons and fast processors, MPC can be deployed at sampling rates of hundreds of hertz or more.
Stability and Feasibility
MPC without a terminal constraint may not guarantee stability. Two standard approaches to ensure stability:
Terminal constraint: enforce \(x[k+N] = 0\) (or \(x[k+N] \in \mathcal{X}_f\), a positively invariant terminal set). This guarantees recursive feasibility and closed-loop stability.
Terminal cost: choose \(P_f\) as the LQR Riccati solution. This ensures that the cost decreases along the optimal trajectory, yielding a Lyapunov-function argument for stability.
Handling Input/State Constraints in Practice
- Soft constraints on states: replace hard inequality \(x \leq x_{\max}\) with a penalty term \(\rho \max(0, x - x_{\max})^2\) in the objective. This avoids infeasibility at the cost of occasional constraint violation.
- Hard input constraints are always included as they reflect physical limits.
- Disturbance handling: Tube MPC, min-max MPC, and stochastic MPC are extensions that provide robustness guarantees in the presence of bounded or stochastic disturbances.
Chapter 9: Worked Design Examples
9.1 End-to-End Design: DC Motor Position Control
Plant. A DC motor with transfer function
\[ G_p(s) = \frac{K_m}{s(\tau s + 1)} = \frac{10}{s(0.1s + 1)}, \]where \(K_m = 10\,\text{rad/V}\cdot\text{s}\) and \(\tau = 0.1\,\text{s}\). This is a type-1 plant (integrator in the loop), so the closed-loop system will have zero steady-state error to step references without additional integral action.
Sampling. Choose \(T = 0.01\,\text{s}\) (10 ms), giving a sampling frequency of 100 Hz, which is approximately 20 times the open-loop bandwidth.
ZOH Discretization.
\[ G_p(z) = (1 - z^{-1})\,\mathcal{Z}\!\left\{\frac{10}{s^2(0.1s+1)}\right\}. \]Expanding via partial fractions:
\[ \frac{10}{s^2(0.1s+1)} = \frac{10}{s^2} - \frac{1}{s} + \frac{1}{s + 10}. \]Z-transform:
\[ \mathcal{Z}\{\cdots\} = \frac{10Tz}{(z-1)^2} - \frac{z}{z-1} + \frac{z}{z-e^{-10T}}. \]With \(T = 0.01\): \(e^{-0.1} \approx 0.9048\), so
\[ \mathcal{Z}\{\cdots\} = \frac{0.1z}{(z-1)^2} - \frac{z}{z-1} + \frac{z}{z-0.9048}. \]Combining over the common denominator \((z-1)^2(z-0.9048)\) and multiplying by \((z-1)/z\) yields
\[ G_p(z) = \frac{(4.679\times10^{-3})z + (4.321\times10^{-3})}{(z-1)(z-0.9048)}. \]Controller Design (Frequency Domain). A proportional-plus-velocity controller \(C(z) = K_p + K_d(1 - z^{-1})/T\) is tuned to place the open-loop gain crossover at \(\omega_c = 20\,\text{rad/s}\) (\(\Omega_c = 0.2\,\text{rad/sample}\)) with \(\text{PM} \geq 50^\circ\). The Bode plot of \(G_p(e^{j\Omega})\) is evaluated numerically, and the compensator gains are selected iteratively.
9.2 Stability Analysis via Jury Test
For the closed-loop characteristic polynomial of the motor example, suppose
\[ p(z) = z^3 - 2.4z^2 + 2.05z - 0.6. \]Apply the Jury test:
- \(p(1) = 1 - 2.4 + 2.05 - 0.6 = 0.05 > 0\). \(\checkmark\)
- \(p(-1) = -1 - 2.4 - 2.05 - 0.6 = -6.05\), and \((-1)^3 p(-1) = 6.05 > 0\). \(\checkmark\)
- \(\lvert a_0 \rvert = 0.6 < a_3 = 1\). \(\checkmark\)
The Jury array (first-row check) would confirm full stability. All three quick conditions are met, indicating stability.
9.3 MPC for a Constrained Double Integrator
Plant (discrete time):
\[ A_d = \begin{bmatrix}1 & T \\ 0 & 1\end{bmatrix} = \begin{bmatrix}1 & 0.1\\0&1\end{bmatrix}, \quad B_d = \begin{bmatrix}T^2/2\\T\end{bmatrix} = \begin{bmatrix}0.005\\0.1\end{bmatrix}. \]Constraints: \(-1 \leq u[k] \leq 1\), \(-5 \leq x_1[k] \leq 5\).
Cost: \(Q = I_2\), \(R = 0.1\), horizon \(N = 10\).
At each step, the QP (dimension 10 in \(\mathbf{U}\)) is solved using an active-set solver. For initial condition \(x[0] = [3, 0]^T\), the MPC drives the state to the origin while saturating the input for the first several steps and then releasing as the state approaches the origin. Without the input constraint, the LQR would command \(\lvert u \rvert > 1\) initially, so MPC strictly outperforms LQR in terms of constraint satisfaction and prevents actuator windup.
Chapter 10: Synthesis — Key Concepts and Design Checklist
10.1 The s-to-z Mapping Summary
The following table summarizes the key correspondences:
| Continuous-time (s-domain) | Discrete-time (z-domain) |
|---|---|
| Open LHP (\(\text{Re}(s) < 0\)) | Interior of unit circle (\(\lvert z \rvert < 1\)) |
| Imaginary axis | Unit circle |
| Open RHP | Exterior of unit circle |
| Pole at \(s = -a\) | Pole at \(z = e^{-aT}\) |
| \(\omega = 0\) | \(z = 1\) |
| \(\omega = \omega_s/2 = \pi/T\) | \(z = -1\) |
| Underdamped pair \(-\zeta\omega_n \pm j\omega_d\) | \(e^{-\zeta\omega_nT}\angle\pm\omega_dT\) |
10.2 Discretization Method Comparison
| Method | \(s \to\) | Stability preserved? | Frequency matching |
|---|---|---|---|
| Forward Euler | \(\frac{z-1}{T}\) | Not guaranteed | Poor |
| Backward Euler | \(\frac{z-1}{Tz}\) | Always | Moderate |
| Tustin | \(\frac{2}{T}\frac{z-1}{z+1}\) | Always | Good (warped) |
| Tustin + pre-warp | (same with adjusted \(\omega\)) | Always | Exact at \(\omega_0\) |
| ZOH (step-invariant) | (via \(\mathcal{Z}\{G(s)/s\}\)) | Always | Exact step response |
10.3 Design Checklist for a Digital Control System
Specify performance: bandwidth \(\omega_c\), phase margin PM, steady-state error requirements (system type), time-domain specs (\(\zeta\), \(\omega_n\), settling time).
Choose sampling rate: \(\omega_s \geq 20\omega_c\) as a rule of thumb; account for anti-aliasing filter requirements.
Model the plant: identify or derive \(G_p(s)\); perform ZOH discretization to obtain \(G_p(z)\).
Design the controller:
- For high \(\omega_s/\omega_c\): design \(C(s)\) in continuous time, then discretize (Tustin preferred).
- For low \(\omega_s/\omega_c\): design \(C(z)\) directly in the discrete-time domain.
Verify stability: Jury test or bilinear + Routh-Hurwitz on the closed-loop characteristic polynomial.
Check margins: evaluate GM and PM from the Bode plot of \(L(e^{j\Omega})\) on the unit circle.
Simulate: verify time-domain step response; check inter-sample behavior with a continuous-time simulation of the sampled-data loop.
Account for delays: include any computational delay \(\Delta\) as \(z^{-\lceil\Delta/T\rceil}\) in the loop model.
Add constraints (if applicable): implement anti-windup for PID or formulate as MPC.
Test robustness: perturb plant parameters, add noise model, verify stability margins are maintained.
Appendix A: z-Transform Derivation from Laplace Transform
Let \(x(t)\) be a continuous-time signal sampled by an ideal impulse sampler:
\[ x^*(t) = x(t) \cdot \sum_{k=0}^{\infty} \delta(t - kT) = \sum_{k=0}^{\infty} x(kT)\,\delta(t-kT). \]The Laplace transform is
\[ X^*(s) = \mathcal{L}\{x^*(t)\} = \sum_{k=0}^{\infty} x(kT)\,e^{-skT}. \]Substituting \(z = e^{sT}\) (equivalently \(e^{-skT} = z^{-k}\)):
\[ X(z) = X^*(s)\Big|_{e^{sT}=z} = \sum_{k=0}^{\infty} x(kT)\,z^{-k}, \]which is the definition of the z-transform. This derivation makes explicit that:
- The z-transform is the Laplace transform of the sampled signal evaluated on the substitution \(z = e^{sT}\).
- The ROC of \(X(z)\) is determined by the convergence of the power series.
- Every horizontal strip of width \(2\pi/T\) in the s-plane maps to the entire z-plane, explaining the periodicity of the discrete-time frequency response.
Appendix B: The Zero-Order Hold Transfer Function
The ZOH holds its input constant over the interval \([kT, (k+1)T)\). Its continuous-time output is
\[ u(t) = u[k] \quad \text{for } kT \leq t < (k+1)T. \]This is equivalent to convolving the impulse train \(u^*(t)\) with a rectangular pulse of duration \(T\):
\[ p(t) = \begin{cases} 1 & 0 \leq t < T \\ 0 & \text{otherwise} \end{cases}, \quad P(s) = \frac{1 - e^{-sT}}{s}. \]The ZOH transfer function in the s-domain is:
\[ H_{\text{ZOH}}(s) = \frac{1 - e^{-sT}}{s}. \]On the imaginary axis:
\[ H_{\text{ZOH}}(j\omega) = T e^{-j\omega T/2} \text{sinc}\!\left(\frac{\omega T}{2\pi}\right), \]where \(\text{sinc}(x) = \sin(\pi x)/(\pi x)\). The ZOH introduces a magnitude distortion (sinc roll-off) and a phase lag of \(\omega T/2\) radians — equivalent to half a sample period of pure delay. At the design stage, this can be approximated as a delay \(e^{-sT/2}\) for frequency-domain analysis.
Appendix C: Ackermann’s Formula — Full Derivation
For a single-input controllable system \((A, B)\), the characteristic polynomial of \(A\) is \(\phi(s) = \det(sI - A) = s^n + \alpha_{n-1}s^{n-1} + \cdots + \alpha_0\). By the Cayley-Hamilton theorem, \(\phi(A) = 0\).
Transform to controllable canonical form via \(T_c = \mathcal{C}\,W\), where \(W\) is the upper-triangular Toeplitz matrix of coefficients:
\[ W = \begin{bmatrix} \alpha_1 & \alpha_2 & \cdots & \alpha_{n-1} & 1 \\ \alpha_2 & \cdots & 1 & 0 \\ \vdots & \iddots & & \vdots \\ 1 & 0 & \cdots & & 0 \end{bmatrix}. \]In the canonical coordinates, the state-feedback gain that assigns the desired polynomial \(\alpha_d(s)\) is straightforward to write down. Transforming back gives
\[ K = e_n^T \mathcal{C}^{-1}\,\alpha_d(A), \]where \(e_n = [0 \; \cdots \; 0 \; 1]^T\). The formula requires only \(\mathcal{C}^{-1}\) and evaluating the polynomial \(\alpha_d\) at the matrix \(A\), both of which are numerically efficient operations.