In previous lessons we have treated PID controllers as continuous-time operators acting on an error signal \( e(t) = r(t) - y(t) \). A common parallel PID structure (with filtered derivative from Lesson 4) is

\[ u(t) = K_p e(t) + K_i \int_0^t e(\tau)\,\mathrm{d}\tau + K_d \dot{e}_f(t), \]

where \( \dot{e}_f(t) \) is a filtered derivative such as

\[ \dot{e}_f(t) = \frac{N}{1 + \frac{1}{N}\frac{\mathrm{d}}{\mathrm{d}t}}\,e(t), \]

with \( N > 0 \) a filter parameter from Lesson 4 that limits high-frequency amplification.

In an embedded realization, signals are sampled, quantized, and processed by a microcontroller. A typical single-loop position or speed controller for a robot joint has the following structure.

flowchart TD
  R["Reference r(k)"] --> E["Error e(k) = r(k) - y(k)"]
  E --> PID["Discrete PID algorithm"]
  PID --> SAT["Actuator limits sat(u(k))"]
  SAT --> PWM["DAC / PWM drive"]
  PWM --> PLANT["Motor / robot joint plant"]
  PLANT --> SENSOR["Encoder / sensor"]
  SENSOR --> ADC["Sampling & quantization"]
  ADC --> E
  

Practical implementation issues arise at almost every block in this chain: sampling and computational delay, numerical approximation of the integral and derivative, actuator saturation, sensor noise and quantization, and finite-precision arithmetic.

Let the sampling period be \( T_s > 0 \), and define discrete-time instants \( t_k = kT_s \) with integer \( k \ge 0 \). The sampled signals are \( e_k = e(t_k) \), \( u_k = u(t_k) \), and \( y_k = y(t_k) \).

A positional discrete-time PID controller can be written as

\[ u_k = K_p e_k + I_k + D_k, \]

where \( I_k \) approximates the integral term and \( D_k \) approximates the derivative term.

2.1 Integral Approximation

A simple (rectangle) rule approximation of the integral is

\[ \int_0^{t_k} e(\tau)\,\mathrm{d}\tau \approx T_s \sum_{j=0}^{k} e_j. \]

Therefore, if we define \( I_k = K_i T_s \sum_{j=0}^{k} e_j \), we obtain the recursion

\[ I_k = I_{k-1} + K_i T_s e_k,\quad k \ge 1,\quad I_0 = 0. \]

Proof of the recursion: By definition,

\[ I_k = K_i T_s \sum_{j=0}^{k} e_j = K_i T_s \sum_{j=0}^{k-1} e_j + K_i T_s e_k = I_{k-1} + K_i T_s e_k. \]

A higher-accuracy trapezoidal rule uses

\[ I_k = I_{k-1} + \frac{K_i T_s}{2}\bigl(e_k + e_{k-1}\bigr), \]

which better approximates the continuous integral but increases implementation complexity slightly (requires access to both \( e_k \) and \( e_{k-1} \)).

2.2 Derivative Approximation and Filter

A backward-difference approximation of the derivative is

\[ \dot{e}(t_k) \approx \frac{e_k - e_{k-1}}{T_s}. \]

Directly using this approximation leads to very strong amplification of high-frequency sensor noise. In Lesson 4 we introduced filtered derivatives; a discrete-time first-order low-pass filter for the derivative can be written as

\[ D_k = \alpha D_{k-1} + (1-\alpha)\frac{e_k - e_{k-1}}{T_s},\quad 0 < \alpha < 1. \]

The parameter \( \alpha \) controls the effective derivative bandwidth. Substituting \( I_k \) and \( D_k \) into the PID law yields a fully discrete-time implementable controller.

Real actuators (motors, valves, thrusters) cannot generate arbitrarily large control signals. Let the actuator limits be \( u_{\min} \le u_k \le u_{\max} \). Define the saturation function

\[ \operatorname{sat}(v) = \begin{cases} u_{\min}, & v < u_{\min},\\[3pt] v, & u_{\min} \le v \le u_{\max},\\[3pt] u_{\max}, & v > u_{\max}. \end{cases} \]

The unsaturated PID output is \( v_k = K_p e_k + I_k + D_k \), while the physical actuator receives \( u_k = \operatorname{sat}(v_k) \). When saturation is active, \( u_k \ne v_k \) and the integral term may continue to grow, a phenomenon known as integral windup.

3.1 Mathematical View of Windup

Assume a constant positive error \( e_k = \bar{e} > 0 \) and an integral update \( I_k = I_{k-1} + K_i T_s \bar{e} \). Then by recursion

\[ I_k = I_0 + k K_i T_s \bar{e}, \]

which grows unbounded as \( k \to \infty \). If the actuator is saturated at \( u_{\max} \), the plant output may not move enough to reduce the error, so the integrator continues to increase. Once the saturation is removed (e.g., when the reference changes), the accumulated \( I_k \) can drive the plant far beyond the target, causing excessive overshoot or even instability.

3.2 Back-Calculation Anti-Windup

A widely used back-calculation anti-windup modifies the integral update to

\[ I_k = I_{k-1} + K_i T_s e_k + K_{aw}\bigl(u_k - v_k\bigr), \]

where \( K_{aw} \ge 0 \) is an anti-windup gain. Note that \( u_k - v_k = 0 \) when the actuator is not saturated, so the extra term only acts under saturation.

Property: When the actuator remains saturated at \( u_k = u_{\max} \) and the unsaturated command \( v_k > u_{\max} \), the anti-windup term tends to reduce the integrator value:

\[ u_k - v_k = u_{\max} - v_k < 0 \quad\Rightarrow\quad K_{aw}(u_k - v_k) < 0, \]

so the integral contribution is driven back toward a value compatible with saturation. This prevents long recoveries and large overshoots once the system comes out of saturation.

3.3 Bumpless Transfer

Another practical issue is the transition between manual and automatic modes (for example, a human operator directly commanding motor voltage, then handing control to the PID). Bumpless transfer means that the automatic controller engages without a large jump in \( u_k \).

For a positional PID, bumpless transfer can be achieved by initializing the internal integral state so that the PID reproduces the current actuator value:

\[ I_0 = u_{\text{manual}} - K_p e_0 - D_0. \]

This guarantees that at the moment of switching, \( u_0 = K_p e_0 + I_0 + D_0 = u_{\text{manual}} \), so the control signal is continuous.

Let the measured output be corrupted by additive noise: \( y_m(t) = y(t) + n(t) \), where \( n(t) \) is measurement noise. Sampling gives \( y_{m,k} = y_k + n_k \).

The derivative term is the most sensitive to noise. Approximating \( \dot{e}(t_k) \) by a finite difference,

\[ \frac{e_k - e_{k-1}}{T_s} = \frac{(r_k - y_k) - (r_{k-1} - y_{k-1})}{T_s} = \frac{(r_k - r_{k-1}) - (y_k - y_{k-1})}{T_s}, \]

but if \( y_k \) is noisy, then even when the true plant output changes slowly, \( y_k - y_{k-1} \) can be dominated by noise. This motivates derivative filtering (Section 2.2) and sometimes measurement path filters (Chapter 26).

4.1 Quantization

Suppose the encoder and ADC together implement a uniform quantizer with step size \( \Delta > 0 \). Then the quantization error \( q_k = y_{m,k} - y_k \) satisfies

\[ -\frac{\Delta}{2} \le q_k \le \frac{\Delta}{2}. \]

This bounded but non-vanishing error can cause small limit cycles where the control signal and output oscillate within one or a few quantization levels. From a design perspective, reducing \( \Delta \) (higher encoder resolution or higher ADC bits) and limiting the derivative gain \( K_d \) alleviate the effect.

4.2 Finite Precision and Overflow

Microcontrollers use finite-precision arithmetic. If \( I_k \) is stored in a fixed-width integer, then windup can cause overflow. A simple numerical safeguard is to clamp the internal integral state:

\[ I_k \in [I_{\min}, I_{\max}], \]

with \( I_{\min}, I_{\max} \) chosen so that \( K_p e_k + I_k + D_k \) remains within the actuator and numeric limits under realistic conditions.

Consider a DC motor driving a single robot joint with angular position \( \theta(t) \). The mechanical dynamics are

\[ J \ddot{\theta}(t) + b \dot{\theta}(t) = K_t i(t), \]

where \( J \) is the equivalent inertia, \( b \) viscous friction, \( K_t \) torque constant, and \( i(t) \) armature current. Neglecting armature inductance for simplicity and modeling the electrical equation as

\[ R i(t) + K_e \dot{\theta}(t) = u(t), \]

with \( R \) armature resistance and \( K_e \) back-EMF constant, we obtain

\[ i(t) = \frac{u(t) - K_e \dot{\theta}(t)}{R}. \]

Substituting into the mechanical equation gives

\[ J \ddot{\theta}(t) + \left(b + \frac{K_t K_e}{R}\right)\dot{\theta}(t) = \frac{K_t}{R}u(t). \]

Taking Laplace transforms with zero initial conditions leads to the transfer function

\[ \frac{\Theta(s)}{U(s)} = \frac{\frac{K_t}{R}}{J s^2 + \left(b + \frac{K_t K_e}{R}\right)s} = \frac{K}{s(\tau_m s + 1)}, \]

where \( K \) and \( \tau_m \) are positive constants that can be identified from motor parameters or experiments.

We implement a discrete-time PID with sampling period \( T_s \), integral and derivative recursions from Section 2, and anti-windup from Section 3. This same structure will be realized in each programming language in the following sections.

In Python, basic PID experiments can be performed using numerical libraries (numpy, scipy) and control-oriented packages such as control or robotics toolboxes. The code below shows:

• Definition of a motor transfer function.
• Discrete-time PID implementation with anti-windup.
• Simulation of step response for a robot joint position command.

import numpy as np
from math import copysign

# Example motor parameters (normalized)
K = 10.0       # dc gain
tau_m = 0.05   # mechanical time constant
Ts = 0.001     # sampling period (1 ms)

# PID gains
Kp = 2.0
Ki = 50.0
Kd = 0.001
alpha = 0.9    # derivative filter coefficient
K_aw = 20.0    # anti-windup gain

u_min, u_max = -12.0, 12.0  # actuator voltage limits

# Simulation horizon
T = 0.5
N = int(T / Ts)

# Allocate arrays
t = np.linspace(0.0, T, N+1)
theta = np.zeros(N+1)      # position
theta_dot = np.zeros(N+1)  # velocity (state)
u = np.zeros(N+1)          # control
r = np.ones(N+1) * 1.0     # 1 radian step

# PID internal states
I = 0.0
D = 0.0
e_prev = 0.0

for k in range(N):
    # Error
    e = r[k] - theta[k]

    # Derivative filter
    D = alpha * D + (1.0 - alpha) * (e - e_prev) / Ts

    # Unsaturated PID output
    v = Kp * e + I + Kd * D

    # Saturation
    u[k] = max(u_min, min(u_max, v))

    # Anti-windup back-calculation
    I = I + Ki * Ts * e + K_aw * (u[k] - v)

    # Motor dynamics: theta_dot_dot = -(1/tau_m)*theta_dot + (K/tau_m)*u
    theta_ddot = -(1.0 / tau_m) * theta_dot[k] + (K / tau_m) * u[k]

    # Integrate dynamics with simple Euler
    theta_dot[k+1] = theta_dot[k] + Ts * theta_ddot
    theta[k+1] = theta[k] + Ts * theta_dot[k+1]

    e_prev = e

# The arrays t, theta, u now contain the simulated response.
# For robotics contexts, one may integrate this PID loop into a ROS node:
# - Use sensor_msgs/JointState for theta
# - Publish u as an effort command through control_msgs/JointJog or similar.
      

In robotics, Python-based frameworks such as ROS with ros_control or MoveIt often implement their own PID loops; however, the internal logic closely follows the structure of the above discrete-time PID.

C++ is common in embedded controllers and in high-performance robotics stacks (e.g., ros_control and real-time joint controllers). The following snippet shows a PID class tailored for a robot joint, suitable for integration into a control loop at period Ts.

class DiscretePID {
public:
    DiscretePID(double Kp_, double Ki_, double Kd_,
                double Ts_, double alpha_, double Kaw_,
                double u_min_, double u_max_)
        : Kp(Kp_), Ki(Ki_), Kd(Kd_), Ts(Ts_),
          alpha(alpha_), Kaw(Kaw_),
          u_min(u_min_), u_max(u_max_),
          I(0.0), D(0.0), e_prev(0.0)
    {}

    double update(double r, double y) {
        double e = r - y;

        // Derivative filter
        D = alpha * D + (1.0 - alpha) * (e - e_prev) / Ts;

        // Unsaturated output
        double v = Kp * e + I + Kd * D;

        // Saturation
        double u = v;
        if (u > u_max) u = u_max;
        if (u < u_min) u = u_min;

        // Anti-windup
        I += Ki * Ts * e + Kaw * (u - v);

        e_prev = e;
        return u;
    }

private:
    double Kp, Ki, Kd, Ts, alpha, Kaw;
    double u_min, u_max;
    double I, D, e_prev;
};

// Example usage in a loop with Ts sampling time:
// DiscretePID pid(2.0, 50.0, 0.001, 0.001, 0.9, 20.0, -12.0, 12.0);
// while (running) {
//     double theta = readEncoderJoint();
//     double r = getReferencePosition();
//     double u = pid.update(r, theta);
//     writeMotorVoltage(u);
// }
      

Within a ros_control controller, the update() method would be called from the real-time control loop, reading joint states from hardware_interface::JointStateHandle and writing efforts or commands to hardware_interface::JointHandle.

Java is used in some robotics platforms (for example, educational mobile robots or competition environments that offer Java APIs). The following code illustrates a similar discrete PID structure in Java.

public class DiscretePID {
    private final double Kp, Ki, Kd, Ts, alpha, Kaw;
    private final double uMin, uMax;
    private double I, D, ePrev;

    public DiscretePID(double kp, double ki, double kd,
                       double ts, double alpha, double kaw,
                       double uMin, double uMax) {
        this.Kp = kp;
        this.Ki = ki;
        this.Kd = kd;
        this.Ts = ts;
        this.alpha = alpha;
        this.Kaw = kaw;
        this.uMin = uMin;
        this.uMax = uMax;
        this.I = 0.0;
        this.D = 0.0;
        this.ePrev = 0.0;
    }

    public double update(double r, double y) {
        double e = r - y;

        // Derivative filter
        D = alpha * D + (1.0 - alpha) * (e - ePrev) / Ts;

        double v = Kp * e + I + Kd * D;

        // Saturation
        double u = Math.max(uMin, Math.min(uMax, v));

        // Anti-windup
        I += Ki * Ts * e + Kaw * (u - v);

        ePrev = e;
        return u;
    }
}

// In a differential-drive robot, one might have two PID instances
// for left and right wheel velocities, using encoders and motor drivers
// provided by the platform's Java robotics library.
      

Java-based robotics frameworks often provide sensor and actuator classes (for example, wheel encoders, motor controllers) that can be wrapped by such PID classes to form closed-loop speed or position controllers.

MATLAB and Simulink are standard tools for control and robotics, with dedicated toolboxes for modeling robot dynamics and designing controllers. The following MATLAB script sets up the DC motor transfer function, designs a PID with fixed gains, and simulates a discretized implementation with saturation and anti-windup.

% Motor parameters (normalized)
K = 10; tau_m = 0.05;
s = tf('s');
P = K / (s*(tau_m*s + 1));

% PID gains (obtained from previous design steps)
Kp = 2.0; Ki = 50.0; Kd = 0.001;
Ts = 0.001; alpha = 0.9; Kaw = 20;

u_min = -12; u_max = 12;

T = 0.5;
N = round(T/Ts);

theta = zeros(1,N+1);
theta_dot = zeros(1,N+1);
u = zeros(1,N+1);
r = ones(1,N+1); % 1 rad step

I = 0; D = 0; e_prev = 0;

for k = 1:N
    e = r(k) - theta(k);

    % Derivative filter
    D = alpha*D + (1-alpha)*(e - e_prev)/Ts;

    v = Kp*e + I + Kd*D;
    u(k) = min(max(v, u_min), u_max);

    % Anti-windup
    I = I + Ki*Ts*e + Kaw*(u(k) - v);

    % Motor dynamics: theta_ddot = -(1/tau_m)*theta_dot + (K/tau_m)*u
    theta_ddot = -(1/tau_m)*theta_dot(k) + (K/tau_m)*u(k);

    theta_dot(k+1) = theta_dot(k) + Ts*theta_ddot;
    theta(k+1) = theta(k) + Ts*theta_dot(k+1);

    e_prev = e;
end

% Plot step response
t = (0:N)*Ts;
figure; 
subplot(2,1,1); plot(t, r, '--', t, theta); grid on;
ylabel('theta (rad)');
subplot(2,1,2); plot(t, u); grid on;
ylabel('u (V)'); xlabel('time (s)');

% In Simulink, one can instead use:
% - A PID Controller block in discrete-time mode
% - Saturation block for actuator limits
% - Anti-windup option in the PID block (back-calculation)
% connected to a motor plant block or Transfer Fcn block P.
      

In robotics, the Robotics System Toolbox and Simscape Multibody can be used to simulate multi-joint manipulators driven by similar joint-level PID controllers.

Wolfram Mathematica offers symbolic and numeric tools for control. We can model the motor, discretize the PID controller, and simulate the closed-loop system using difference equations or state-space representations.

(* Parameters *)
K = 10.0;
tauM = 0.05;
Ts = 0.001;
Kp = 2.0; Ki = 50.0; Kd = 0.001;
alpha = 0.9; Kaw = 20.0;
uMin = -12.0; uMax = 12.0;

tFinal = 0.5;
nSteps = Round[tFinal/Ts];

(* Difference equations for motor *)
theta = ConstantArray[0.0, nSteps + 1];
thetaDot = ConstantArray[0.0, nSteps + 1];
u = ConstantArray[0.0, nSteps + 1];
r = ConstantArray[1.0, nSteps + 1];

I = 0.0; D = 0.0; ePrev = 0.0;

Do[
  e = r[[k]] - theta[[k]];

  D = alpha*D + (1 - alpha)*(e - ePrev)/Ts;

  v = Kp*e + I + Kd*D;
  u[[k]] = Min[Max[v, uMin], uMax];

  I = I + Ki*Ts*e + Kaw*(u[[k]] - v);

  thetaDDot = -(1.0/tauM)*thetaDot[[k]] + (K/tauM)*u[[k]];
  thetaDot[[k + 1]] = thetaDot[[k]] + Ts*thetaDDot;
  theta[[k + 1]] = theta[[k]] + Ts*thetaDot[[k + 1]];

  ePrev = e;
, {k, 1, nSteps}];

ListLinePlot[
  {
    Transpose[{Range[0, nSteps]*Ts, r}],
    Transpose[{Range[0, nSteps]*Ts, theta}]
  },
  PlotLegends -> {"reference", "theta"},
  Frame -> True, FrameLabel -> {"time (s)", "position (rad)"}
]
      

In more complex robotic systems, symbolic derivation of dynamics (e.g., via Lagrange equations) can be combined with numeric PID control implemented as above.

Problem 1 (Integral Windup Quantification): Consider a discrete-time PI controller without anti-windup:

\[ u_k = K_p e_k + I_k,\quad I_k = I_{k-1} + K_i T_s e_k. \]

The actuator saturates at \( u_{\max} \), and the error remains constant at \( e_k = \bar{e} > 0 \) for \( k = 0,1,\dots,N \). Assume that the controller output is saturated for all these steps, i.e., \( u_k = u_{\max} \) for \( 0 \le k \le N \). Derive an explicit expression for \( I_N \) in terms of \( I_0, K_i, T_s, \bar{e}, N \), and compute the extra control effort \( u_N - u_{\max} \) that would occur if saturation were suddenly removed at step \( N \).

Solution: The integral recursion can be solved explicitly:

\[ I_k = I_{k-1} + K_i T_s \bar{e} \quad\Rightarrow\quad I_k = I_0 + k K_i T_s \bar{e}. \]

In particular,

\[ I_N = I_0 + N K_i T_s \bar{e}. \]

If saturation were removed at step \( N \), the unsaturated output would be

\[ u_N^{\text{unsat}} = K_p \bar{e} + I_N = K_p \bar{e} + I_0 + N K_i T_s \bar{e}. \]

Hence the extra control effort relative to the actuator limit is

\[ u_N^{\text{unsat}} - u_{\max} = \bigl(K_p \bar{e} + I_0 - u_{\max}\bigr) + N K_i T_s \bar{e}, \]

which grows linearly with \( N \) and can be very large if windup is not mitigated.

Problem 2 (Effect of Sampling Period on Derivative Noise): Let the noise component of the error signal be a zero-mean white sequence \( n_k \) with variance \( \sigma_n^2 \). The derivative estimate is

\[ d_k = \frac{n_k - n_{k-1}}{T_s}. \]

Assuming \( n_k \) are independent, compute the variance of \( d_k \) and interpret its dependence on \( T_s \).

Solution: Using independence and zero mean,

\[ \begin{aligned} \operatorname{Var}(d_k) &= \operatorname{Var}\!\left(\frac{n_k - n_{k-1}}{T_s}\right) = \frac{1}{T_s^2}\operatorname{Var}(n_k - n_{k-1}) \\[3pt] &= \frac{1}{T_s^2}\bigl(\operatorname{Var}(n_k) + \operatorname{Var}(n_{k-1})\bigr) = \frac{2\sigma_n^2}{T_s^2}. \end{aligned} \]

Thus the derivative noise variance grows like \( 1/T_s^2 \). Very small \( T_s \) (high sampling frequency) can dramatically increase the noise level in the derivative estimate unless appropriate filtering or gain reduction is used.

Problem 3 (Back-Calculation Anti-Windup Equilibrium): For a constant reference \( r \) and constant plant output \( y \), suppose that the actuator is saturated at \( u_k = u_{\max} \) and the unsaturated output is \( v_k = K_p e + I_k \) with constant error \( e = r - y \). The integral update uses back-calculation:

\[ I_k = I_{k-1} + K_i T_s e + K_{aw}(u_{\max} - v_{k-1}). \]

Show that if \( e \) remains constant and the controller reaches a steady state under saturation, then \( v_k \) converges to \( u_{\max} \).

Solution: At steady state, we require

\[ I_k = I_{k-1} = I^{\star},\quad v_k = v_{k-1} = v^{\star}. \]

The integral update equation becomes

\[ 0 = K_i T_s e + K_{aw}(u_{\max} - v^{\star}), \]

so

\[ v^{\star} = u_{\max} + \frac{K_i T_s}{K_{aw}} e. \]

On the other hand, by definition \( v^{\star} = K_p e + I^{\star} \) and \( u_k = u_{\max} \) regardless of \( v^{\star} \). If \( K_{aw} \) is chosen sufficiently large, the term \( \frac{K_i T_s}{K_{aw}} e \) becomes small, so \( v^{\star} \approx u_{\max} \). In the exact limit \( K_{aw} \to \infty \), we obtain \( v^{\star} = u_{\max} \) exactly, meaning that the internal unsaturated command is kept equal to the saturated actuator output at equilibrium.

Problem 4 (Implementation Flow for a Robot Joint): Sketch a logical flow for a joint-level PID loop that includes sampling, PID computation, saturation, and anti-windup. Indicate the order of operations required in each sampling period.

Solution (flow):

flowchart TD
  A["Start of sample k"] --> B["Read sensor theta(k)"]
  B --> C["Compute error e(k) = r(k) - theta(k)"]
  C --> D["Update derivative and integral states"]
  D --> E["Compute unsaturated command v(k)"]
  E --> F["Apply saturation to get u(k)"]
  F --> G["Apply anti-windup using u(k) - v(k)"]
  G --> H["Send u(k) to actuator (PWM/DAC)"]
  H --> I["Wait until next sampling instant"]
  

The key point is that anti-windup uses the difference between saturated and unsaturated commands, so saturation must be applied before updating the integral state.

In this lesson we bridged the gap between ideal continuous-time PID control and actual embedded implementations for mechatronic and robotic systems. We derived discrete-time recursions for the integral and derivative terms, highlighted how actuator saturation leads to integral windup, and introduced back-calculation anti-windup and bumpless transfer as essential remedies. We also analyzed the impact of sampling, noise, quantization, and finite precision, and implemented a joint-level PID controller for a DC motor case study in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. These practical tools prepare us for more advanced analysis of frequency-domain behavior, robustness, and linear controller architectures in subsequent chapters.