PSE-823: Advanced Process Dynamics and Control¶

Lecture 15 -- Dead-Time Compensation and Model Predictive Control¶

(Cecil, Ch. 8, Sec. 8.2-8.3)¶

Date: __________

Agenda¶

Part A -- Dead-time compensation: the paper machine (Sec. 8.2)

Part B -- Step-response models and the dynamic matrix (Sec. 8.3)

Part C -- MPC on the jacketed reactor: horizons, move suppression, constraints

Part D -- The model inside MPC, identified from data

Setup (run once)¶

  • The standard imports, plus tools for delays, matrices, optimization and regression
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
import control as ct
from scipy.integrate import solve_ivp
from scipy.optimize import curve_fit
plt.rcParams.update({"font.size": 15})
from collections import deque
from scipy.optimize import brentq
from scipy.linalg import toeplitz
from gekko import GEKKO
from sklearn.linear_model import LinearRegression, RidgeCV

Part A¶

Dead-Time Compensation¶

(Cecil, Ch. 8, Sec. 8.2)

Introduce -- a process dominated by dead time (Fig. 8.8)¶

  • Paper machine: basis weight measured at the reel, about a minute after the stock valve acts.
  • Model: dead time $\theta = 1.2$ min, time constant $\tau = 0.17$ min, gain $K = 0.67$ %/%.
  • Disturbance: actual consistency 2% -> 2.2% (+10% basis weight).
  • PID must wait two dead times before seeing any correction.
  • Minutes matter: at 1600 ft/min, 8 min is over 2 miles of paper.

Derive 1/2 -- PID limits (Fig. 8.9)¶

Crossover: $\tan^{-1}(0.17\omega) + 1.2\omega = \pi$ (Lecture 8):

Cecil: $K_u = 1.6$ %/%, $P_u = 2.7$ min. Quarter decay: $K_C = K_u/2 = 0.8$, $T_I = P_u = 2.7$ min. Lambda tuning with $\tau_{CL} = \theta$: $K_C = \frac{\tau}{\tau_{CL} + \theta}\frac1K = 0.11$, $T_I = \tau = 0.17$ min.

Derive 2/2 -- the Smith predictor (Figs. 8.14-8.15)¶

Run a model beside the process; feed the controller a prediction:

$$\tilde W(t + \theta) = W(t) + [\hat W(t + \theta) - \hat W(t)]$$

$\hat W(t + \theta)$: model output before its dead time; $\hat W(t)$: after it. The bracket is the effect of moves made in the last dead time, not yet measured. With a perfect model the controller sees only the lag: tune it hard.

Code 1/3 -- ultimate gain and period¶

  • Mirrors Derive 1/2 with brentq (Lecture 8)
In [2]:
K, tau, theta = 0.67, 0.17, 1.2
w = brentq(lambda w: np.arctan(tau*w) + theta*w - np.pi, 0.1, 10)
Ku = np.sqrt(1 + (tau * w)**2) / K
print(f"w_co = {w:.3f} rad/min, Ku = {Ku:.2f} %/%, "
      f"Pu = {2*np.pi/w:.2f} min")
print(f"Lambda: Kc = {tau / (2*theta) / K:.3f}, TI = {tau} min")

w_co = 2.307 rad/min, Ku = 1.60 %/%, Pu = 2.72 min
Lambda: Kc = 0.106, TI = 0.17 min

Code walkthrough. The crossover equation is solved as in Lecture 8; Ku = 1/AR there, with AR = K/sqrt(1 + (tau w)^2). Prints w_co = 2.307 rad/min, Ku = 1.60 %/% and Pu = 2.72 min (Cecil: 1.6 and 2.7). Lambda tuning gives Kc = 0.17/(2.4)/0.67 = 0.106 (Cecil rounds to 0.11) and TI = 0.17 min.

Code 2/3 -- the loop, with and without the predictor¶

  • Dead times as queues (Lecture 8); PI with set point 0
  • smith=True adds the bracket of Derive 2/2 to the PV
In [3]:
def paper(Kc, TI, smith=False, th_m=1.2, dt=0.005, t_end=12):
    sheet = deque([0.0] * round(theta / dt))    # true dead time
    model = deque([0.0] * round(th_m / dt))     # model dead time
    x = xm = I = 0.0; W = []
    for _ in range(int(t_end / dt)):
        pv = sheet[0] + (xm - model[0] if smith else 0)   # Fig. 8.15
        I += -pv * dt; u = Kc * (-pv + I / TI)            # PI
        x += dt * (K * (u + 15) - x) / tau      # 15% load from t = 0
        xm += dt * (K * u - xm) / tau           # model, no load
        sheet.append(x); model.append(xm); model.popleft()
        W.append(sheet.popleft())
    return np.array(W)

Code walkthrough. sheet carries the true weight through the machine; model carries the model's. Each step: the PV is the measurement sheet[0], plus (Smith) the model output before its delay minus after it. The PI controller acts on -PV. The process lag gets the controller output plus the load: 15% of valve equivalent, i.e. +10% weight (0.67 x 15) if uncontrolled. The model gets the controller output only -- it cannot know the consistency changed. W records the measured weight.

Code 3/3 -- three controllers (Fig. 8.9 and beyond)¶

  • IAE over 12 min in the legend
In [4]:
t = np.arange(int(12 / 0.005)) * 0.005
fig, ax = plt.subplots(figsize=(7, 5.25))
for name, args in [("quarter decay", (0.8, 2.7)),
                   ("Lambda", (0.11, 0.17)),
                   ("Smith + PI", (1.0, 0.17, True))]:
    W = paper(*args)
    ax.plot(t, W, label=f"{name}: IAE {np.trapezoid(abs(W), t):.1f}")
ax.set(xlabel="t (min)", ylabel="basis weight change (%)"); ax.legend()
plt.show()
No description has been provided for this image

Code walkthrough. All three see nothing for one dead time, then the weight jumps toward +10%. Quarter decay: IAE 40.7 and still off by 1.4% at 12 min (TI of 2.7 min is far too long, as Cecil says). Lambda: smooth, IAE 25.6, lined out by about 7 min (Cecil: about 8). Smith predictor with an aggressive PI (Kc = 1.0): IAE 14.6, within 0.5% after 3.4 min -- close to the two-dead-time limit of 2.4 min.

Build on it -- the model's dead time is 20% too long¶

  • Machine speed changed; the predictor was not updated
In [ ]:
W = paper(1.0, 0.17, smith=True, th_m=1.44)
print(f"IAE = {np.trapezoid(abs(W), t):.2f}, "
      f"last time |W| > 0.5: {t[abs(W) > 0.5].max():.2f} min")

Code walkthrough. Same controller, model dead time 1.44 min against the true 1.2. IAE rises from 14.6 to 18.75 and the response is outside +/-0.5% until 8.35 min, worse than Lambda tuning's line-out. The predictor's benefit depends on an accurate dead time (Cecil: on a paper machine it varies with machine speed, so it must be computed from the speed).

Your Turn (6 min)¶

  1. By hand: Lambda tuning with $\tau_{CL} = 2\theta$ instead of $\theta$. New $K_C$?
  2. Read the code: what does this print?
q = deque([0.0, 0.0, 0.0]); out = []
for x in [1.0, 2.0, 3.0, 4.0]:
    q.append(x); out.append(q.popleft())
print(out)

(1) Kc = 0.17/(2.4 + 1.2)/0.67 = 0.070 %/%: slower still. (2) [0.0, 0.0, 0.0, 1.0]: each value comes out three steps later -- a dead time of three samples.

Check -- cold-call¶

The operators beat the PID by: read the scan, turn the stock valve by the "calibrated" amount, then wait a dead time. Which part of the Smith predictor does the waiting correspond to?

The bracket W-hat(t + theta) - W-hat(t): the memory of moves already made but not yet seen. It stops the controller from piling on more correction while the first one is still travelling down the machine.

Part B¶

Step-Response Models and the Dynamic Matrix¶

(Cecil, Ch. 8, Sec. 8.3)

Introduce -- a model as a list of numbers (Table 8.3)¶

  • Reactor with a once-through jacket (Fig. 8.21): cooling water $108.8 \to 88.8$ lb/min; temperature $150 \to 154.5$ F.
  • Step response per unit input: $s(k) = c(k)/(-20)$, every $\Delta t = 15$ min; flat after $N = 13$ points.
  • Impulse response $g(k) = s(k) - s(k-1)$.
  • Free-form: no structure assumed. Linear, so superposition holds.

Derive 1/2 -- prediction by superposition¶

Write any input as steps $\Delta m(j)$; add their responses:

$$c(k) = \sum_j s(k - j)\,\Delta m(j)$$

Pulse of +40 lb/min for 30 min: $\Delta m(0) = +40$, $\Delta m(2) = -40$ (Table 8.4).

Derive 2/2 -- the dynamic matrix and QDMC¶

Effect of $L$ future moves on $R$ future outputs:

$$\mathbf x = A\,\Delta\mathbf m,\quad A = \begin{bmatrix} s(1) & 0\\ s(2) & s(1)\\ \vdots & \vdots\\ s(R) & s(R-1)\end{bmatrix}$$

Choose $\Delta\mathbf m$ to cancel the predicted errors $\hat{\mathbf e}$ in least squares: $\Delta\mathbf m = (A^TA)^{-1}A^T\hat{\mathbf e}$. Implement only the first move; repeat next sample.

Code 1/3 -- Table 8.3¶

  • From the step test to $s(k)$ and $g(k)$
In [6]:
c_test = np.array([0, 1.0, 1.9, 2.5, 3.0, 3.4, 3.7, 3.9, 4.1, 4.2,
                   4.3, 4.4, 4.4, 4.5])           # F, every 15 min
s = c_test / -20                                  # F per (lb/min)
g = np.diff(s)
print("s:", s.round(3)); print("g:", g.round(3))

s: [-0.    -0.05  -0.095 -0.125 -0.15  -0.17  -0.185 -0.195 -0.205 -0.21
 -0.215 -0.22  -0.22  -0.225]
g: [-0.05  -0.045 -0.03  -0.025 -0.02  -0.015 -0.01  -0.01  -0.005 -0.005
 -0.005  0.    -0.005]

Code walkthrough. The test changed the flow by -20 lb/min, so dividing by -20 gives the response to +1 lb/min. Prints s = 0, -0.05, -0.095, ..., -0.225 and g = -0.05, -0.045, -0.03, ..., with s(11) = s(12) (g(12) = 0): the 0.1 F resolution of the thermometer shows up as "quantization noise" in the model. np.diff gives the 13 differences g(k).

Code 2/3 -- predict the pulse test (Table 8.4)¶

  • Mirrors Derive 1/2: two shifted, scaled copies of $s$
In [7]:
def s_at(k):                      # 0 before k = 0, flat after 13
    return np.where(k < 0, 0.0, s[np.clip(k, 0, 13)])
k = np.arange(19)
C = 150 + 40 * s_at(k) - 40 * s_at(k - 2)
print(C.round(1))

[150.  148.  146.2 147.  147.8 148.2 148.6 149.  149.2 149.4 149.6 149.6
 149.8 149.8 149.8 150.  150.  150.  150. ]

Code walkthrough. s_at extends the 14 stored values: zero for negative k (before the step), the final value beyond k = 13 (the model is finite). The pulse is +40 at k = 0 and -40 at k = 2 (30 min later). Prints 150.0, 148.0, 146.2, 147.0, 147.8, ... 150.0: Table 8.4's C(k) column. The real test read 148.6 and 147.4 at the first two samples: the reactor is not perfectly linear.

Code 3/3 -- the dynamic matrix, $R = 4$, $L = 2$¶

  • Mirrors Derive 2/2; toeplitz builds the shifted columns
In [8]:
def dyn_matrix(s, R, L):
    return toeplitz(s[1:R + 1], np.zeros(L))
A = dyn_matrix(s, 4, 2)
G = np.linalg.solve(A.T @ A, A.T)                  # (A^T A)^-1 A^T
print(A); print(G.round(3))
print("moves:", (G @ np.full(4, 5.0)).round(1))

[[-0.05   0.   ]
 [-0.095 -0.05 ]
 [-0.125 -0.095]
 [-0.15  -0.125]]
[[-14.052  -8.39   -0.343   3.616]
 [ 18.309   9.09   -3.052  -9.316]]
moves: [-95.8  75.1]

Code walkthrough. toeplitz(col, row) makes a matrix with constant diagonals: first column s(1)..s(R), first row zeros -- exactly A. np.linalg.solve computes (A^T A)^(-1) A^T without forming the inverse. G matches Cecil: [[-14.05, -8.39, -0.343, 3.616], [18.31, 9.09, -3.052, -9.316]]. For a +5 F set-point step every predicted error is 5: moves -95.8 and +75.1 lb/min. Only -95.8 is implemented.

Your Turn (8 min)¶

  1. By hand: square DMC uses only the first row: $\Delta m = \hat e(1)/s(1)$. For the 5 F step with $\Delta t = 15$ min, what move? With $\Delta t = 1$ min, $s(1) = 0$: what happens?
  2. Read the code: what does it print?
print(toeplitz([1, 2, 3], [1, 0]))

(1) 5/(-0.05) = -100 lb/min (flow 108.8 -> 8.8). With s(1) = 0 the move is infinite: the dead time (1.35 min) hides the first sample. QDMC needs R > int(theta/dt) + L. (2) [[1 0] [2 1] [3 2]].

Check -- think-pair-share¶

MPC computes $L$ future moves every sample but implements only the first. Why throw the rest away?

Receding horizon: next sample brings a new measurement (and model error); re-solving uses it. The later moves were based on older information. It is feedback, implemented through repeated optimization.

Part C¶

MPC on the Jacketed Reactor¶

(Cecil, Ch. 8, Sec. 8.3)

Introduce -- making MPC behave¶

  • A perfect-model controller is as aggressive as its math allows: the cooling flow would go negative.
  • Knobs: horizon $R$, number of moves $L$, and move suppression: $\Phi = \|\hat{\mathbf e} - A\Delta\mathbf m\|^2 + k^2\|\Delta\mathbf m\|^2$.
  • Constraints: $0 \le$ flow, i.e. $\Delta M \ge -108.8$ lb/min.
  • Cecil's smoothed model: $K = -0.2312$, $\tau = 55.49$ min, $\theta = 1.353$ min; $\Delta t = 1$ min, $R = 60$, $L = 2$.

Derive 1/1 -- move suppression¶

Set $\partial\Phi/\partial\Delta\mathbf m = 0$:

$$\Delta\mathbf m = (A^TA + k^2I)^{-1}A^T\hat{\mathbf e}$$

Ridge regression on the moves: $k$ trades tracking for smaller moves. $k = 0$ is QDMC; large $k$ gives timid moves.

Code 1/4 -- smooth the model by fitting (Table 8.7)¶

  • FOPDT least squares on the 14 points of $s(k)$ (Lecture 10)
In [9]:
def fopdt(t, K, T, th):
    return K * (1 - np.exp(-np.clip(t - th, 0, None) / T))
p, _ = curve_fit(fopdt, 15.0 * np.arange(14), s, p0=[-0.23, 50, 1])
print("K, tau, theta =", p.round(4))
s1 = fopdt(np.arange(400.0), *p)               # every 1 min
print("s(1), s(2) =", s1[1:3].round(4))

K, tau, theta = [-0.2311 55.4868  1.308 ]
s(1), s(2) = [-0.     -0.0029]

Code walkthrough. The same curve_fit as Lecture 10, on Table 8.3 at 15-min spacing. Prints K = -0.2311, tau = 55.49 min, theta = 1.31 min (Cecil: -0.2312, 55.49, 1.353). The fitted curve, sampled every minute, is a smooth step-response model with no quantization noise. s(1) = 0 because the dead time exceeds one sample; s(2) = -0.003.

Code 2/4 -- a QDMC controller in 12 lines¶

  • Model and plant both as step responses (superposition)
  • Model error added to every prediction; flow constraint by clipping
In [10]:
def qdmc(sm, sp, R, L, k, n=120, r=5.0, lo=-108.8):
    A = dyn_matrix(sm, R, L)
    gain = np.linalg.solve(A.T @ A + k**2 * np.eye(L), A.T)[0]
    cm, cp = np.zeros(n + R + 1), np.zeros(n + R + 1)   # model, plant
    u, y, M = 0.0, [], []
    for i in range(n):
        e = r - (cm[i+1:i+R+1] + cp[i] - cm[i])   # predicted errors
        dm = max(gain @ e, lo - u); u += dm       # first move, clipped
        cm[i:] += dm * sm[:n+R+1-i]; cp[i:] += dm * sp[:n+R+1-i]
        y.append(cp[i]); M.append(u)
    return np.array(y), np.array(M)

Code walkthrough. gain is the first row of (A^T A + k^2 I)^(-1) A^T: only the first move is ever used. cm and cp hold the future trajectories implied by all moves so far, for the model and the plant. Each sample: model error = measured cp[i] minus predicted cm[i] (Cecil's EM), added to the R predictions; errors against the set point; first move, clipped so the flow stays >= 0; both trajectories updated by superposition.

Code 3/4 -- set point 150 -> 155 F¶

  • Perfect model; $k = 0$ against $k = 0.17$ (Figs. 8.31, 8.33)
In [11]:
fig, (a1, a2) = plt.subplots(2, 1, figsize=(7, 5.6), sharex=True)
for k in [0.0, 0.17]:
    y, M = qdmc(s1, s1, 60, 2, k)
    a1.plot(150 + y, label=f"k = {k}: IAE {np.abs(5 - y).sum():.0f}")
    a2.plot(108.8 + M)
a1.set(ylabel="T (F)"); a1.legend()
a2.set(xlabel="t (min)", ylabel="flow (lb/min)"); plt.show()
No description has been provided for this image

Code walkthrough. k = 0: the flow drops straight to zero (the constraint) and stays there until the temperature is nearly at 155; tracking is fast, IAE 39. k = 0.17: the first move is only -26.5 lb/min, the flow never goes below about 60, and the temperature approaches 155 more slowly (IAE 122). With a perfect model, suppression only costs speed...

Build on it -- ...but the model is never perfect¶

  • Plant gain 50% higher, dead time 5 min; same two controllers
In [ ]:
sp_bad = fopdt(np.arange(400.0), 1.5 * p[0], p[1], 5.0)
for k in [0.0, 0.17]:
    y, M = qdmc(s1, sp_bad, 60, 2, k)
    print(f"k = {k}: peak {150 + y.max():.2f} F, "
          f"IAE {np.abs(5 - y).sum():.0f}")

Code walkthrough. Same controllers; the plant is now more sensitive and slower to respond than the model. k = 0: peak 157.58 F and a sustained oscillation (IAE 322). k = 0.17: peak 155.0 F, IAE 82 -- better than with the perfect model, because the extra plant gain speeds up the cautious controller. Move suppression buys robustness to model error, Cecil's reason for using it.

Code 4/4 -- the same problem in GEKKO¶

  • Model as an ODE, flow bounds built in, moves penalized (DCOST)
  • Dead time (1.35 min against $\tau = 55$ min) neglected here
In [13]:
m = GEKKO(remote=False); m.time = np.linspace(0, 60, 61)
Mv = m.MV(value=0, lb=-108.8, ub=91.2); Mv.STATUS = 1; Mv.DCOST = 0.3
Cv = m.CV(value=0); Cv.STATUS = 1; Cv.SP = 5; Cv.TR_INIT = 0
m.Equation(p[1] * Cv.dt() == -Cv + p[0] * Mv)   # FOPDT, no dead time
m.options.IMODE = 6; m.options.CV_TYPE = 2       # MPC, squared error
m.solve(disp=False)
T = 150 + np.array(list(Cv.value))
print("planned:", np.round(Mv.value[1:5], 1), "T(10), T(60):",
      T[[10, 60]].round(2))

planned: [-38.3 -64.7 -82.  -92.7] T(10), T(60): [153.25 154.96]

Code walkthrough. GEKKO builds and solves the optimization: MV is the manipulated variable with bounds, STATUS = 1 lets the optimizer move it, DCOST penalizes each move (like k); CV is the controlled variable with set point 5; the FOPDT is an ODE in Cv.dt(). IMODE 6 is MPC; CV_TYPE 2 squares the error. remote=False solves on this machine. Planned flow changes (from 108.8) -38.3, -64.7, -82.0, -92.7: a gradual descent toward the lower bound; T reaches 153.25 F at 10 min, 154.96 at 60. Constraints are part of the problem, not clipped after.

Your Turn (6 min)¶

  1. By hand: at steady state, how much must the cooling flow change for +5 F? Compare with the sum of the two QDMC moves of Part B, $-95.8 + 75.1$.
  2. Read the code: in qdmc, what would happen without the term cp[i] - cm[i]?

(1) 5/(-0.225) = -22.2 lb/min; the two moves sum to -20.7: close, and closer as R grows (Cecil). (2) No model-error feedback: the controller would act on the model alone. With a mismatched plant (or a disturbance) it would settle with an offset -- the term is MPC's integral action.

Check -- poll¶

In GEKKO, raising DCOST from 0.3 to 30 would make the first planned flow change:

A) larger B) smaller C) unchanged, only later moves change

B: a much smaller first change than -38.3 lb/min. DCOST penalizes every move, the first included; the controller spreads the change over time and reaches the set point more slowly.

Part D¶

The Model Inside MPC, Identified from Data¶

(Cecil, Ch. 8, Sec. 8.3)

Introduce -- a short, noisy plant test¶

  • Cecil: the step-response coefficients can be fitted by least squares to any test data (finite impulse response, FIR).
  • 60 impulse coefficients at $\Delta t = 5$ min; only 150 samples of a random +/-10 lb/min test; 0.1 F noise and 0.1 F resolution.
  • Plain least squares fits the noise. Ridge regression (Lecture 10's ML tools, plus a penalty) shrinks the coefficients.
  • Question: which model makes the better controller?

Derive 1/1 -- FIR regression and its penalty¶

$c_k = \sum_{j=1}^{60} g_j\,m_{k-j}$: linear in $g$:

Least squares: $\min \|c - Xg\|^2$. Ridge: $\min \|c - Xg\|^2 + \alpha\|g\|^2$ -- the same form as move suppression, now on the model. RidgeCV picks $\alpha$ by cross-validation.

Code 1/3 -- the test data and the regressors¶

  • Column $j$ of $X$ is the input delayed by $j$ samples
In [14]:
rng = np.random.default_rng(3)
n, N, dt5 = 150, 60, 5.0
m_test = np.repeat(rng.choice([-10.0, 10.0], n // 6 + 1), 6)[:n]
g_true = np.diff(fopdt(dt5 * np.arange(N + 1), *p))
c_meas = np.convolve(m_test, np.r_[0, g_true])[:n]
c_meas = np.round(c_meas + 0.1 * rng.standard_normal(n), 1)
X = np.column_stack([np.r_[np.zeros(j), m_test[:n - j]]
                     for j in range(1, N + 1)])
print(X.shape, c_meas[:5])

(150, 60) [-0.1 -0.3 -0.4 -0.5 -0.7]

Code walkthrough. The input holds each random +/-10 for 6 samples (30 min). g_true is the true impulse response (the fitted reactor at 5-min spacing); np.convolve applies the FIR sum; noise is added and the result rounded to 0.1 F like the real thermometer. np.r_ concatenates; column j of X is the input delayed j samples, so X g is the FIR prediction. 150 rows, 60 unknowns.

Code 2/3 -- least squares versus ridge¶

  • Two scikit-learn regressors; cumulative sum turns $g$ into $s$
In [15]:
fits = {"least squares": LinearRegression(fit_intercept=False),
        "ridge": RidgeCV(alphas=np.logspace(-2, 5, 30),
                         fit_intercept=False)}
s_id = {}
for name, mdl in fits.items():
    g_hat = mdl.fit(X, c_meas).coef_
    s_id[name] = np.r_[0, np.cumsum(g_hat)]
    print(f"{name}: gain {s_id[name][-1]:.3f}")
print("alpha chosen:", round(fits["ridge"].alpha_))

least squares: gain -0.230
ridge: gain -0.227
alpha chosen: 1172

Code walkthrough. fit_intercept=False: deviation variables, no offset term. coef_ holds the 60 impulse coefficients; their running sum is the step-response model s(k). Both give a gain near the true -0.231 (-0.230 and -0.227), but the least-squares s(k) is jagged and the ridge one is smooth (the next plot). Cross-validation chose alpha = 1172.

Code 3/3 -- put each model in the controller¶

  • QDMC, $\Delta t = 5$ min, $R = 12$, $L = 2$, $k = 0$; true plant
In [16]:
s_plant = fopdt(dt5 * np.arange(400), *p)
fig, ax = plt.subplots(figsize=(7, 5.25))
for name, sm in [("true", s_plant), *s_id.items()]:
    sm = np.r_[sm, np.full(400 - len(sm), sm[-1])][:400]
    y, M = qdmc(sm, s_plant, 12, 2, 0.0, n=40)
    tv = np.abs(np.diff(np.r_[0, M])).sum()
    ax.plot(dt5 * np.arange(40), 150 + y,
            label=f"{name}: IAE {np.abs(5-y).sum():.0f}, "
                  f"moves {tv:.0f}")
ax.set(xlabel="t (min)", ylabel="T (F)"); ax.legend(); plt.show()
No description has been provided for this image

Code walkthrough. Each identified s(k) is padded with its final value to the length qdmc needs, then used as the controller's model against the true plant. "moves" is the total variation, the sum of |change in flow|: how hard the valve works. True model: IAE 10, moves 203. Least squares: IAE 41, moves 2887, and the temperature swings grow -- the noise in s(k) becomes valve thrashing. Ridge: IAE 15, moves 587. Regularizing the model did what smoothing did for Cecil.

Your Turn (6 min)¶

Read the code: what does each print?

print(np.r_[0, np.cumsum([1.0, 2.0, 3.0])])
print(np.convolve([1.0, 0.0, 0.0], [0.0, 0.5, 0.3])[:3])

Which line turns an impulse response into a step response?

[0. 1. 3. 6.] and [0. 0.5 0.3]. The first: s(k) is the running sum of g(k). The second shows a unit pulse at k = 0 producing the impulse response itself.

Check -- cold-call¶

Both identified models had nearly the right gain. Why did one make a much worse controller?

MPC uses the whole s(k) shape, especially the early coefficients that set the first moves; noise there is amplified by an aggressive controller (Cecil's quantization-noise example). A good model for control needs the right dynamics, not just the right gain: fit quality and control quality are different tests (Lectures 10-11).

Wrap-up¶

Derived, then coded: the dead-time limit on PID and the Smith predictor; step and impulse response models and superposition; the dynamic matrix and QDMC; move suppression and constraints; GEKKO MPC; FIR identification by least squares and by ridge regression, judged inside the controller.

Tools: deque as a dead time, toeplitz, np.linalg.solve, receding-horizon loops, GEKKO (MV, CV, IMODE = 6), RidgeCV

Course end: from first-order lags to model predictive control.