Controller discretization: sample, compute, and hold

Explore controller discretization with exact held-input dynamics. Compare immediate and delayed commands, calculate discrete poles, and connect sample timing to control stability.

By 11 min read

What you will learn

  • State when a controller samples, computes, and applies its command.
  • Derive the exact discrete model of a first-order plant with held input.
  • Calculate how one sample of delay changes closed-loop poles.
  • Distinguish sampling limits, controller approximations, and numerical integration.
  • Check discrete stability before interpreting a finite response plot.

Before you start

A digital controller acts on measurements that arrive at distinct times. The actuator often holds its last command while the controller waits for the next measurement. That timing can turn a stable continuous design into an unstable sampled loop.

You will derive the exact update for a small speed-control model. Then you will add one sample of delay and check its poles. MathWorks documents zero-order hold as an exact plant conversion for staircase inputs; the timing of the feedback connection still needs its own model.

Define controller discretization and timing

Use a horizontal joint with speed ω, commanded torque u, inertia J = 1 kg·m², and viscous damping B = 1 N·m·s/rad:

J ω̇ = u − Bω

The target steps from zero to r = 1 rad/s at t = 0. The joint starts at rest. This teaching model assumes exact speed measurements and ideal torque actuation, with no saturation, friction beyond B, or mechanical limits.

At sample k, time is t = kh, where h is the sample period in seconds. The controller reads ω[k] = ω(kh) and computes c[k] = K(r[k] − ω[k]). K has units N·m·s/rad.

We compare two schedules:

  • Immediate update: apply c[k] during [kh, (k+1)h).
  • One-sample delay: apply c[k−1] during that same interval. Initialize c[−1] = 0.

The second schedule models a full sample of latency. Real hardware can have a fractional delay or varying latency; those cases need a different timing model.

Follow a command between samples

A zero-order hold, or ZOH, keeps u constant until the next update. Speed continues changing between measurements. With elapsed time θ after the latest sample:

ω(kh + θ) = exp(−Bθ/J) ω[k] + [1 − exp(−Bθ/J)] u[k]/B
0 ≤ θ < h

For K = 2 and h = 0.2 s, immediate control starts with c[0] = u[0] = 2 N·m. Halfway through that hold, speed is 2(1 − exp(−0.1)) = 0.190325 rad/s. The controller still uses its sample of zero and keeps applying 2 N·m.

At t = 0.2 s, speed reaches 0.362538 rad/s. The controller then computes 2(1 − 0.362538) = 1.274923 N·m. Speed stays continuous when the command changes.

Derive the exact held-input plant

Evaluate the hold solution at θ = h:

ω[k+1] = aω[k] + bu[k]
a = exp(−Bh/J), b = (1 − a)/B

Here a is dimensionless. The coefficient b converts torque into speed, so its units are rad/(N·m·s). With J = B = 1 in the stated units and h = 0.2 s, a = 0.818731 and b = 0.181269.

This update exactly matches the plant at sample times under the hold assumption. A finer numerical integration step cannot remove the effects of h: the controller still receives one measurement per sample period.

Write the sampled plant in z

For a sequence x[k], the one-sided z-transform is X(z) = Σ x[k]z⁻ᵏ, starting at k = 0. A one-sample delay with zero stored history contributes z⁻¹. Taking the transform of the recurrence with ω[0] = 0 gives:

zΩ(z) = aΩ(z) + bU(z)
Gd(z) = Ω(z)/U(z) = b/(z − a) = bz⁻¹/(1 − az⁻¹)

The numerator z⁻¹ records that u[k] first affects ω[k+1]. Adding computation delay creates another delay in the loop. The plant already includes the hold interval.

A discrete mode evolves like zᵏ. Its magnitude decays when |z| < 1. MathWorks states this open-unit-disk condition for stable discrete-time models; continuous-time poles use the open left half-plane.

Close the loop with immediate updates

Substitute u[k] = K(r[k] − ω[k]):

ω[k+1] = (a − bK)ω[k] + bKr[k]
Ω(z)/R(z) = bK/(z − a + bK)
p = a − bK

For K ≥ 0, asymptotic stability requires K < (1+a)/b. At equality, p = −1: perturbations alternate without decaying. Beyond that boundary, their magnitude grows.

The default values give p = 0.456192. If the loop is stable, its final speed is K/(B+K) × r = 2/3 rad/s. Proportional control leaves a steady error of 1/3 rad/s because damping needs a nonzero torque command.

For comparison, continuously updated proportional feedback gives ω(t) = [K/(1+K)] [1 − exp(−(1+K)t)] with our numerical J and B. That controller changes its torque at every instant. Sampling its already-closed response would describe a different input from the staircase command in this experiment.

Add one sample of computation delay

With delayed actuation, u[k] = c[k−1]. The first interval applies zero torque. For k ≥ 1:

ω[k+1] = aω[k] + bK(r[k−1] − ω[k−1])
Ω(z)/R(z) = bK/(z² − az + bK)
p₁,₂ = [a ± √(a² − 4bK)]/2

For 0 < a < 1 and K ≥ 0, both roots lie inside the unit circle exactly when bK < 1. The real quadratic's three strict stability inequalities reduce to bK < 1, since 1 − a + bK and 1 + a + bK are already positive.

Use h = 1 s and K = 2. Immediate control has pole −0.896362, which is stable. Delayed control has poles 0.183940 ± 1.109237i, each with magnitude 1.124385, so it is unstable.

The exact gain bounds are 2.163953 for immediate control and 1.581977 for delayed control. At the delayed boundary bK = 1, the distinct conjugate modes persist without decaying. A unit-circle boundary does not establish bounded-input, bounded-output stability.

The algebraic equilibrium still exists beyond the boundary. The lab reports a steady limit only when the poles establish asymptotic stability. It labels poles within floating-point roundoff of the unit circle as near the boundary and withholds that limit.

Choose a sample period with a model

Start with the plant dynamics and required control bandwidth. Include sensing, computation, communication, and actuator timing. Then check the sampled loop's poles and its response between samples.

For the delayed K = 2 example, stability needs 2(1 − exp(−h)) < 1, or h < ln 2 ≈ 0.693147 s. This exact bound belongs to this plant and timing rule. Acceptable tracking can demand a shorter period than the stability bound.

Sampling also limits what the sensor can distinguish. At 10 samples per second, cosine signals at 8 Hz and 2 Hz give the same samples:

cos(2π × 8k/10) = cos(2π × 2k/10)

Analog Devices explains anti-alias filtering before conversion. The Nyquist limit concerns reconstructing signals under stated spectral assumptions. It does not guarantee acceptable feedback delay, tracking, or stability. This lab assumes a clean speed sensor and does not simulate aliasing.

Separate model conversion from integration

Three choices enter a digital-control implementation:

  • Plant discretization: derive a sample-to-sample model for the assumed input hold. We use an exact ZOH model here.
  • Controller discretization: convert a controller with internal dynamics into difference equations. A filtered derivative needs an update rule; the static proportional controller here needs no numerical integration.
  • Simulation integration: approximate the continuous plant between events when an exact solution is unavailable. The integrator must respect command-change times.

For a dynamic controller, common substitutions include forward Euler s ≈ (z−1)/h, backward Euler s ≈ (z−1)/(hz), and Tustin s ≈ 2(z−1)/(h(z+1)). They produce different discrete dynamics. MathWorks compares ZOH and Tustin conversion assumptions, including their time- and frequency-domain goals.

Check the final connected loop after choosing the rule and sample period. A controller's discrete update, the input hold, and any delay all contribute to that loop.

Inspect the sampled response

At the default t = 2 s, the current speed is 0.666406 rad/s and the immediate command is 0.667187 N·m. Move the cursor between samples to compare current speed with the latest measurement. The held input stays fixed until the next event.

Choose Slow sampling, then Unstable delayed loop. Both use K = 2 and h = 1 s. The changed timing alone moves the poles outside the unit circle.

Inspect the held command

What changes when the command arrives later?

Compare the latest measurement, the new command, and the torque currently reaching the joint.

2.00
2.00

J = 1 kg·m², B = 1 N·m·s/rad, target = 1 rad/s. Initial speed and stored command are zero. The plant follows its exact exponential solution under each held input. The model includes no actuator limit or sensor noise.

Sample, compute, select, and hold one control commandAt time kh, read speed omega k and compute c k equal to K times reference minus sampled speed. Apply the newly computed command c k immediately. Hold that selected torque until the next sample at k plus one times h. Speed evolves continuously during the hold.At time khRead ω[k]Compute c[k]K(r[k] − ω[k])Apply c[k]Hold u[k]kh ≤ t < (k+1)hNext sample: ω[k+1]
The measurement comes first. The new command reaches the actuator at the same sample time. The readouts show the command selected just after an event.
Exact sampled-loop speed and continuous-feedback comparisonAt 2.000000 seconds, joint speed is 0.666406 radians per second. Blue follows exact motion between samples, with dots at sample events. Orange dashes show the distinct continuously updated proportional controller. The horizontal axis runs from zero to eight seconds. Values are not clipped; the vertical scale includes the complete plotted response. Speed curves join exact values at twenty subdivisions per hold interval.01202468Speed (rad/s)Time (s)
Blue: sampled control with exact motion between measurement dots. Orange dashes: continuously updated proportional control. Gray dashes: the 1 rad/s target. These controllers apply different torque signals.
Computed commands and the held actuator inputAt 2.000000 seconds, held torque is 0.667187 newton meters and the latest computed command is 0.667187 newton meters. Blue stairs show the held input; orange crosses mark newly computed commands. The horizontal axis runs from zero to eight seconds. Values are not clipped; the vertical scale includes the complete plotted response. Speed curves join exact values at twenty subdivisions per hold interval.01.5302468Torque (N·m)Time (s)
Blue stairs: held actuator torque. Orange crosses: freshly computed commands. With delay, each new command reaches the actuator one sample later. At the final event, the new hold starts at 8 s and extends beyond the plot.
Sample index k
10
Latest sampled speed (rad/s)
0.666406
Current speed (rad/s)
0.666406
Computed command (N·m)
0.667187
Held command (N·m)
0.667187
Plant coefficient a
0.818731
Plant coefficient b
0.181269
Closed-loop poles
0.456192
Largest pole magnitude
0.456192
Stability
Asymptotically stable
Stable steady speed (rad/s)
0.666667
Stable steady error (rad/s)
0.333333

A pole magnitude above one means the model is unstable, even if this short trace looks modest. Values are not clipped; large responses rescale the axes. Steady limits require all poles strictly inside the unit circle. Within numerical roundoff of that boundary, the limit is not reported.

The blue speed curve follows exact held-input motion, with dots at measurement times. The orange speed curve uses continuously updated proportional feedback. In the command plot, blue steps show applied torque and orange marks show newly computed commands.

The plots include every computed value without saturation or clipping. Unstable settings can produce large speeds and commands, so the axes rescale. A short response window can look modest even when a pole already proves instability.

Reproduce the timing in Python

This standard-library script implements the sample, compute, apply, and hold sequence. Each stored row describes the start of an interval, after selecting its command.

import math

def simulate(gain, period, delayed, intervals):
    a = math.exp(-period)
    b = 1 - a
    speed, previous_command = 0.0, 0.0
    rows = []
    for k in range(intervals + 1):
        command = gain * (1 - speed)
        applied = previous_command if delayed else command
        rows.append((speed, command, applied))
        speed = a * speed + b * applied
        previous_command = command
    return rows

a, b = math.exp(-0.2), 1 - math.exp(-0.2)
print(f"Plant: a={a:.6f}, b={b:.6f}")
print(f"Immediate pole: {a - 2*b:.6f}; stable limit: {2/3:.6f}")
speed, command, applied = simulate(2, 0.2, False, 10)[10]
print(f"At 2 s: speed={speed:.6f}, command={command:.6f}")
a, b = math.exp(-1), 1 - math.exp(-1)
print(f"At h=1: immediate pole={a-2*b:.6f}, delayed radius={math.sqrt(2*b):.6f}")
speed, command, applied = simulate(2, 1, True, 2)[2]
print(f"Delayed at 2 s: speed={speed:.6f}, computed={command:.6f}, held={applied:.6f}")
print(f"Unstable K=4 at 8 s: {simulate(4, 1, False, 8)[8][0]:.6f}")
print(f"Aliased samples: {math.cos(2*math.pi*8/10):.6f}, {math.cos(2*math.pi*2/10):.6f}")

Expected output:

Plant: a=0.818731, b=0.181269
Immediate pole: 0.456192; stable limit: 0.666667
At 2 s: speed=0.666406, command=0.667187
At h=1: immediate pole=-0.896362, delayed radius=1.124385
Delayed at 2 s: speed=1.264241, computed=-0.528482, held=2.000000
Unstable K=4 at 8 s: -379.117636
Aliased samples: 0.309017, 0.309017

Try it yourself

Exercise 1. Use K = 2 and h = ln 2 s with the same plant. Find a and b, the immediate pole, and the delayed pole magnitudes. Which loop has decaying perturbations?

Check the exact stability boundary

Here a = 1/2 and b = 1/2 in the stated units. The immediate pole is 1/2 − (1/2)2 = −1/2, so perturbations decay while alternating sign.

The delayed polynomial is z² − 0.5z + 1. Its roots are 0.25 ± (√15/4)i, each with magnitude 1. The delayed loop sits exactly on the boundary and its free modes do not decay.

Exercise 2. At h = 1 s and K = 2, delayed control starts with c[−1] = 0 and ω[0] = 0. Compute ω[1], c[1], ω[2], c[2], and the held command just after t = 2 s.

Check which command reaches the actuator

At t = 0, c[0] = 2 while u[0] = 0, so ω[1] = 0. At t = 1, c[1] = 2 and the actuator applies c[0] = 2 for one second.

Thus ω[2] = 2(1 − exp(−1)) = 1.264241 and c[2] = 2(1 − ω[2]) = −0.528482. Just after t = 2, the actuator still applies c[1] = 2 N·m. The new negative command reaches it at t = 3.

Sources and further study

Next, apply an explicit sampled recurrence to derivative filtering. Include its extra state and timing when you recheck the full control loop.