Digital Control Design

Digital Control Design (gnc-autonomy/control/digital-control-design)

Use when the task is designing a sampled-data control loop in the z-domain: the plant is continuous and must be discretized with a zero-order hold, the compensator is continuous and must be emulated in discrete time, or the loop is implemented as a discrete PID. This leaf implements the sampled-data control methods in pure Python, stdlib only, deterministic and offline: ZOH step-invariant discretization with exact closed-form coefficients, Tustin bilinear emulation with frequency prewarping, discrete PID coefficient forms, unit-circle stability and the sample-rate rule. It pairs with the continuous s-domain design leaves (pid-control-design, lead-lag-compensation) that produce the continuous controllers this leaf samples, and with state-space-analysis for the continuous state-space context.

Domain quick reference

• ZOH first-order plant G(s) = a/(s + a): the step-invariant map gives A = exp(-a*T) and B = 1 - A, so A + B == 1.0 exactly and the sampled DC gain is the continuous one (1.0). zoh_first_order(a, T).
• ZOH second-order plant G(s) = wn^2/(s^2 + 2zetawns + wn^2), underdamped (0 < zeta < 1), in phase-variable companion form with C = [1, 0]: with sigma = exp(-zetawnT), omega_d = wnsqrt(1 - zeta^2) the discrete state matrix is the exact closed-form matrix exponential A11 = sigma*(cos(wdT) + (zetawn/wd)sin(wdT)), A12 = sigmasin(wdT)/wd, A21 = -sigmawn^2sin(wdT)/wd, A22 = sigma(cos(wdT) - (zetawn/wd)sin(wdT)), and B = [1 - A11, sigmawn^2sin(wdT)/wd]. The step response settles to 1.0 and the sampled natural frequency equals the continuous omega_d. zoh_second_order(zeta, wn, T). Assumption recorded: the spec hint "A11 = sigmacos(omega_d*T)" describes the scaled companion variant; the phase-variable form above is the standard textbook realization and satisfies every stated validation identity.
• Tustin (bilinear) emulation substitutes s = (2/T)(z - 1)/(z + 1), or with frequency prewarping at wc the constant c = wc/tan(wcT/2) replaces 2/T; the prewarped map sends s = jwc exactly to z = exp(jwcT), so the emulated phase matches the continuous phase at wc. A continuous pole s_p maps to z = (c + s_p)/(c - s_p). Prewarping requires wcT < pi (wc below the folding frequency). tustin_emulate(cont_coeffs, T, wc = None) -> {num_z, den_z} with a monic denominator; tustin_frequency_check(cont_coeffs, z_coeffs, wc, T) reports the phase error in degrees at wc.
• Position-form discrete PID: u(k) = Kpe(k) + KiTsum(e) + Kd(e(k) - e(k-1))/T, so the difference-equation coefficients are {kp: Kp, ki: Ki*T, kd: Kd/T}. discrete_pid_position(Kp, Ki, Kd, T).
• Velocity-form discrete PID: delta_u(k) = u(k) - u(k-1) = b0e(k) + b1e(k-1) + b2e(k-2) with b0 = Kp + KiT + Kd/T, b1 = -Kp - 2Kd/T, b2 = Kd/T and a1 = -1 on u(k-1) (the 1 - z^-1 delta denominator). Identity b0 + b1 + b2 = KiT: a constant error injects Ki*T of control per step. discrete_pid_velocity(Kp, Ki, Kd, T).
• Unit-circle stability: a sampled pole is stable only when its modulus is strictly below 1; a pole exactly on the unit circle is unstable. unit_circle_poles(den_z) -> {poles, stable} with exact closed forms for degree 1 and 2 denominators (a 1e-12 boundary epsilon keeps numerically unit-modulus poles from being mislabeled stable; higher degrees raise ValueError, numeric root finding is out of scope).
• Sample-rate rule (minimum-rate rule, not a band): sample 10-20 times per closed-loop cycle, w_s_min = SAMPLING_RULE_LOWwb with SAMPLING_RULE_LOW = 10 and T_max = 2pi/w_s_min. Sampling at or faster than the rule minimum is acceptable: verdict "ok" when T <= T_max, "too-slow" when T > T_max. Assumption recorded: the spec's model line lists the signature as sample_rate_rule(wb) but the verdict and validation list require the candidate sample period, so the implemented signature is sample_rate_rule(wb, T) -> {w_s_min_rad_s, t_max_s, verdict}.

Workflow

1. Discretize the plant: zoh_first_order(a, T) for the first-order plant a/(s + a), or zoh_second_order(zeta, wn, T) for the underdamped second-order plant; the returned (A, B) drive the difference equation x(k+1) = Ax(k) + Bu(k).
2. Emulate the continuous compensator with tustin_emulate, passing the real num and den coefficient lists in descending powers of s and the sample time T; give wc when the emulation must match phase at the crossover target (frequency prewarping).
3. Check the emulation: tustin_frequency_check reports the phase error in degrees at wc between the continuous and the emulated compensator; the spec bound is under 1 deg for a prewarped emulation.
4. For a discrete PID, take discrete_pid_position or discrete_pid_velocity coefficients and write the difference equations from the forms above (velocity form increments u each step, so it is the natural form for incremental actuators).
5. Check stability: unit_circle_poles on the sampled denominator (plant loop denominator or emulated compensator) returns the stability verdict against the unit circle.
6. Select the sample time: sample_rate_rule(wb, T) with the closed-loop bandwidth wb returns T_max and the ok / too-slow verdict.
7. Confirm the deterministic checks with the contract test scripts/test_digital_control_design.py.

Worked example

Lead compensator D(s) = 20*(s + 20)/(s + 200) (DC gain 2) on plant G(s) = 10/(s + 10), crossover target wc = 10 rad/s, T = 0.01 s. Real module outputs:

• zoh_first_order(10, 0.01) -> A = 0.9048374180, B = 0.0951625820; A is exp(-0.1) and A + B == 1.0 exactly (sampled DC gain 1).
• zoh_second_order(0.5, 10, 0.01) -> A = [[0.9951665847, 0.0095004083], [-0.9500408335, 0.9001625014]], B = [0.0048334153, 0.9500408335]; the sampled step response settles to 1.0000000000 and the sampled natural frequency is 8.660254 rad/s, matching the continuous omega_d = 10*sqrt(0.75) exactly.
• tustin_emulate with prewarp at wc = 10: num_z = [10.9962478112, -8.9954139914], den_z = [1.0, 0.0004169099]; the DC gain at z = 1 is 2.0000000000 and tustin_frequency_check returns 6.36e-15 deg, far under the 1 deg bound. Plain Tustin gives den_z = [1.0, 0.0], the continuous pole at -200 rad/s mapped to z = 0.
• discrete_pid_position(2.0, 1.0, 0.1, 0.01) -> {kp: 2.0, ki: 0.01 (Ki*T), kd: 10.0 (Kd/T)}.
• discrete_pid_velocity(2.0, 1.0, 0.1, 0.01) -> {b0: 12.01, b1: -22.0, b2: 10.0, a1: -1.0}; b0 + b1 + b2 = 0.01 = Ki*T.
• unit_circle_poles: pole z = 0.5 -> stable True; z = 1.05 -> stable False; z = exp(j*0.3) (unit modulus) -> stable False under the strict rule.
• sample_rate_rule(10, 0.01) -> {w_s_min_rad_s: 100.0, t_max_s: 0.0628318531, verdict: "ok"}; at T = 0.1 s (slower than T_max) the verdict is "too-slow".

Verification

• zoh_first_order: A + B == 1.0 exactly for any valid (a, T); A equals exp(-a*T); ValueError for a <= 0 or T <= 0.
• zoh_second_order: step response settles to 1.0 (DC-gain identity); sampled natural frequency within 1% of continuous omega_d; ValueError for wn <= 0, zeta outside (0, 1), T <= 0.
• tustin_emulate: monic denominator; DC gain preserved; continuous pole mapping z = (c + s_p)/(c - s_p); prewarped phase error under 1 deg at wc; ValueError for T <= 0, empty coefficient lists, wc <= 0 when given, and wc*T >= pi.
• discrete PID forms: worked-example coefficient identities hold; the velocity-form identity b0 + b1 + b2 = Ki*T holds; ValueError for T <= 0.
• unit_circle_poles: truth table for z = 0.5 (stable), z = 1.05 (unstable), z = exp(j*0.3) (boundary, unstable); ValueError for empty denominators and for degree 3 and higher.
• sample_rate_rule: verdict "ok" at T = 0.01 s for wb = 10 (T < T_max) and exactly at T = T_max; "too-slow" at T = 0.1 s; ValueError for wb <= 0 or T <= 0.
• All functions deterministic (no RNG, identical floats run to run); the contract test re-runs every function and compares exact results.

Pitfalls

• Prewarping without the folding constraint: prewarping requires wc*T < pi (wc below the folding frequency); wc*T >= pi (and wc <= 0 or T <= 0) raises ValueError, and plain Tustin at a fast pole can collapse the denominator (the -200 rad/s pole maps to z = 0 without prewarp).
• Skipping the phase check: tustin_frequency_check reports the phase error at wc and the spec bound is under 1 deg for a prewarped emulation; verify it before trusting the emulated phase at crossover.
• Calling a unit-circle pole stable: a sampled pole is stable only strictly inside the unit circle - a pole exactly on it (z = exp(j*0.3), or unit modulus within the 1e-12 epsilon) is labeled unstable by the strict rule.
• Velocity-form PID arithmetic: the delta coefficients obey b0 + b1 + b2 = KiT with a1 = -1 on u(k-1); use the velocity form for incremental actuators and the position form (ki = KiT, kd = Kd/T) otherwise.
• Sampling slower than the rule: sample 10-20 times per closed-loop cycle (w_s_min = 10wb, T_max = 2pi/w_s_min); the verdict is too-slow when T > T_max.
• The signature is sample_rate_rule(wb, T) (the candidate sample period is required for the verdict), and unit_circle_poles handles degree 1-2 denominators exactly while higher degrees raise ValueError.

Related leaves

• gnc-autonomy/control/pid-control-design: continuous PID tuning sibling; its s-domain gains become the Kp, Ki, Kd sampled here.
• gnc-autonomy/control/lead-lag-compensation: continuous compensator sibling; its D(s) is the input to tustin_emulate.
• cross-cutting/numerics/digital-filter-design: signal-filter sibling (IIR Butterworth design); distinct scope, it shapes signals, not control loops.
• gnc-autonomy/control/state-space-analysis: continuous state-space context for the plants discretized by the zero-order hold.

Behavior contract (gate 3)

Run the deterministic contract test (stdlib unittest, offline):

python3 scripts/test_digital_control_design.py

(35 test methods, deterministic, under a second). It covers the ZOH first- and second-order discretization against the worked example with the exact A + B == 1 identity and the settle-to-1 DC-gain identity, the 1% sampled-frequency match, plain and prewarped Tustin emulation with the closed-form pole mapping and the under-1-deg phase-error bound at wc, the discrete PID position and velocity coefficient identities (b0 + b1 + b2 = Ki*T), the unit-circle truth table including the unit-modulus boundary case, the sample-rate rule verdicts at T = 0.01 and T = 0.1 s, exact dict key sets, ValueError rejection of every non-physical input, and run-to-run determinism.

Compliance

• Standards referenced, not reproduced: ARP4754A (development assurance of aircraft functions) supplies the assurance context; every equation above is standard sampled-data control methodology, summary-only per standards-map.yaml.
• compliance: STANDARDS-REF, gated: false.