Chapter 12: Advanced PID Design and Implementation

Lesson 1: PID Design from Time-Domain Specifications

In this lesson we formalize how to design PID controllers directly from time-domain performance specifications such as overshoot, settling time, and steady-state accuracy. Building on earlier chapters (time response of second-order systems, steady-state error, and basic PID control), we derive analytic relationships between time-domain requirements and closed-loop pole locations, then map those to PID gains via polynomial matching in the Laplace domain. Finally, we implement the complete workflow in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

1. Unity-Feedback PID Problem Formulation

We consider a standard single-input single-output unity-feedback loop, where the controller is PID and the plant is a linear time-invariant (LTI) system described by its transfer function. The reference input is \( r(t) \), the output is \( y(t) \), the error is \( e(t) = r(t) - y(t) \), and the control signal is \( u(t) \).

The continuous-time PID controller in parallel form is

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

or, in the Laplace domain (assuming zero initial conditions),

\[ C(s) = \frac{U(s)}{E(s)} = K_p + \frac{K_i}{s} + K_d s = \frac{K_d s^2 + K_p s + K_i}{s}. \]

Let the plant transfer function be \( G(s) \). Under unity feedback, the closed-loop transfer function from reference \( R(s) \) to output \( Y(s) \) is

\[ T(s) = \frac{Y(s)}{R(s)} = \frac{C(s)G(s)}{1 + C(s)G(s)}. \]

Our design objective in this lesson is:

  • Given a plant \( G(s) \) and time-domain specifications (e.g., maximum overshoot, settling time, allowable steady-state error),
  • compute PID gains \( K_p, K_i, K_d \) such that the closed-loop step response satisfies those specifications.
flowchart TD
  R["Reference r(t)"] --> E["Sum: e = r - y"]
  Y["Output y(t)"] --> E
  E --> CBLK["PID C(s)"]
  CBLK --> U["Control u(t)"]
  U --> GBLK["Plant G(s)"]
  GBLK --> Y
        

2. Time-Domain Specifications and Second-Order Prototypes

For servo-type systems, time-domain specifications are usually expressed on the closed-loop response to a unit-step reference. For an underdamped second-order prototype with damping ratio \( \zeta \) and natural frequency \( \omega_n \),

\[ T_{\mathrm{2nd}}(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad 0 < \zeta < 1, \;\; \omega_n > 0. \]

The corresponding unit-step response has the well-known performance metrics:

  • Peak time \( T_p \):

    \[ T_p = \frac{\pi}{\omega_d} = \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}, \]

  • Maximum overshoot \( M_p \) (relative, not percentage):

    \[ M_p = \mathrm{e}^{-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}}. \]

  • Settling time (2% criterion) \( T_s \):

    \[ T_s \approx \frac{4}{\zeta \omega_n}. \]

2.1. Solving for \( \zeta \) From Overshoot

In design, we often know the allowed maximum overshoot \( M_p^\star \) (e.g., 10% → \( M_p^\star = 0.1 \)) and want to compute the corresponding damping ratio \( \zeta \). Starting from

\[ M_p^\star = \mathrm{e}^{-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}}, \]

take natural logarithms:

\[ \ln M_p^\star = -\frac{\pi\zeta}{\sqrt{1-\zeta^2}}. \]

Rearranging,

\[ -\ln M_p^\star = \frac{\pi\zeta}{\sqrt{1-\zeta^2}} \;\Longrightarrow\; \left(-\ln M_p^\star\right)^2 = \frac{\pi^2 \zeta^2}{1-\zeta^2}. \]

Solve for \( \zeta^2 \):

\[ \left(-\ln M_p^\star\right)^2 (1-\zeta^2) = \pi^2 \zeta^2 \;\Longrightarrow\; \left(-\ln M_p^\star\right)^2 = \zeta^2\left[\pi^2 + \left(-\ln M_p^\star\right)^2\right]. \]

\[ \zeta^2 = \frac{\left(-\ln M_p^\star\right)^2}{\pi^2 + \left(-\ln M_p^\star\right)^2}, \quad \zeta = \frac{-\ln M_p^\star}{\sqrt{\pi^2 + \left(\ln M_p^\star\right)^2}}. \]

2.2. Solving for \( \omega_n \) From Settling Time

Given a desired 2% settling time \( T_s^\star \) and a chosen \( \zeta \), we invert the approximate relation

\[ T_s^\star \approx \frac{4}{\zeta \omega_n} \;\Longrightarrow\; \omega_n \approx \frac{4}{\zeta T_s^\star}. \]

Thus, once \( M_p^\star \) and \( T_s^\star \) are specified, we can compute \( \zeta \) and \( \omega_n \), and thereby the desired dominant closed-loop poles

\[ s_{1,2} = -\zeta\omega_n \pm j \omega_n\sqrt{1-\zeta^2}. \]

The essence of time-domain PID design is to choose \( K_p, K_i, K_d \) so that the closed-loop characteristic polynomial has these complex conjugate poles as its dominant roots.

3. PID Design via Polynomial Matching for a Simple Plant

To make the design ideas concrete, we focus on a common second-order plant with one pole at the origin (type-1 behavior), e.g., a simplified DC motor speed or position model:

\[ G(s) = \frac{1}{s(s+1)}. \]

The PID controller has transfer function

\[ C(s) = K_p + \frac{K_i}{s} + K_d s = \frac{K_d s^2 + K_p s + K_i}{s}. \]

The open-loop transfer function is

\[ L(s) = C(s)G(s) = \frac{K_d s^2 + K_p s + K_i}{s} \cdot \frac{1}{s(s+1)} = \frac{K_d s^2 + K_p s + K_i}{s^2 (s+1)}. \]

The closed-loop characteristic equation is obtained from \( 1 + L(s) = 0 \):

\[ 1 + \frac{K_d s^2 + K_p s + K_i}{s^2 (s+1)} = 0 \;\Longrightarrow\; s^2 (s+1) + K_d s^2 + K_p s + K_i = 0. \]

Expanding:

\[ s^3 + s^2 + K_d s^2 + K_p s + K_i = 0 \;\Longrightarrow\; s^3 + (1 + K_d) s^2 + K_p s + K_i = 0. \]

Hence the closed-loop characteristic polynomial is

\[ \Delta_{\text{cl}}(s) = s^3 + (1 + K_d) s^2 + K_p s + K_i. \]

3.1. Desired Third-Order Dominant Pole Polynomial

Our closed-loop system is third order. A standard design trick is:

  • Enforce that the two dominant poles match a desired second-order prototype determined by \( \zeta, \omega_n \),
  • Place the third (real) pole sufficiently far to the left in the complex plane so that its effect is small.

Let the third pole be at \( s = -p_3 \) with \( p_3 > 0 \). Then the desired characteristic polynomial is

\[ \Delta_{\text{des}}(s) = \bigl(s^2 + 2\zeta\omega_n s + \omega_n^2\bigr)\,(s + p_3). \]

Multiply out:

\[ \begin{aligned} \Delta_{\text{des}}(s) &= s^3 + (2\zeta\omega_n + p_3) s^2 + (\omega_n^2 + 2\zeta\omega_n p_3) s + \omega_n^2 p_3. \end{aligned} \]

3.2. Coefficient Matching and Explicit PID Formulas

We now equate the coefficients of \( \Delta_{\text{cl}}(s) \) and \( \Delta_{\text{des}}(s) \). Specifically,

\[ s^3 + (1 + K_d) s^2 + K_p s + K_i \equiv s^3 + (2\zeta\omega_n + p_3) s^2 + (\omega_n^2 + 2\zeta\omega_n p_3) s + \omega_n^2 p_3. \]

Thus we obtain the linear system

\[ \begin{cases} 1 + K_d = 2\zeta\omega_n + p_3, \\ K_p = \omega_n^2 + 2\zeta\omega_n p_3, \\ K_i = \omega_n^2 p_3. \end{cases} \]

Solving:

\[ \boxed{ \begin{aligned} K_d &= 2\zeta\omega_n + p_3 - 1, \\ K_p &= \omega_n^2 + 2\zeta\omega_n p_3, \\ K_i &= \omega_n^2 p_3. \end{aligned} } \]

3.3. Choice of the Third Pole \( p_3 \)

To preserve the second-order behavior, the real pole at \( -p_3 \) must be sufficiently faster (more negative) than the dominant pair. A common heuristic is

\[ p_3 \geq 3\zeta\omega_n \quad\text{or even}\quad p_3 \geq 5\zeta\omega_n, \]

which places the third pole at least 3 to 5 times further left than the real part \( -\zeta\omega_n \) of the dominant poles. For this lesson, we will adopt

\[ p_3 = \alpha \zeta\omega_n, \quad \alpha \in [3,5]. \]

This completes the analytic mapping \( (M_p^\star, T_s^\star, \alpha) \mapsto (\zeta, \omega_n, p_3) \mapsto (K_p,K_i,K_d) \) for our example plant \( G(s) = 1/(s(s+1)) \).

4. Step-by-Step Design Algorithm

We summarize the PID design from time-domain specifications for \( G(s) = 1/(s(s+1)) \) as an algorithm. The same idea extends to other plants when their closed-loop characteristic polynomial can be expressed as a linear function of \( K_p,K_i,K_d \).

  1. Specify desired maximum overshoot \( M_p^\star \) and settling time \( T_s^\star \) (2% criterion).
  2. Compute damping ratio \( \zeta \) from \( M_p^\star \) using the formula in Section 2.1.
  3. Compute \( \omega_n \) from \( \zeta \) and \( T_s^\star \) using Section 2.2.
  4. Choose third pole factor \( \alpha \) (e.g. 4) and set \( p_3 = \alpha \zeta\omega_n \).
  5. Compute \( K_p,K_i,K_d \) from the closed-loop coefficient matching formulas in Section 3.2.
  6. Validate the design by simulating the step response and verifying that the actual \( M_p, T_s \) are within specification.
flowchart TD
  A["Specs: Mp*, Ts*"] --> B["Compute zeta from Mp*"]
  B --> C["Compute omega_n from Ts* and zeta"]
  C --> D["Choose alpha, set p3 = alpha * zeta * omega_n"]
  D --> E["Match coefficients to get Kp, Ki, Kd"]
  E --> F["Simulate step response"]
  F --> G{"Specs satisfied?"}
  G -->|yes| H["Accept PID gains"]
  G -->|no| C2["Adjust Mp*, Ts*, or alpha, repeat"]
        

5. Numerical Design Example

Consider again \( G(s) = 1/(s(s+1)) \). Suppose we require:

  • Maximum overshoot \( M_p^\star = 0.1 \) (10%),
  • Settling time \( T_s^\star = 2\,\mathrm{s} \) (2% criterion).

5.1. Compute \( \zeta \) and \( \omega_n \)

First compute \( \zeta \). Let \( M_p^\star = 0.1 \). Then

\[ \ln M_p^\star = \ln(0.1) \approx -2.3026, \quad -\ln M_p^\star \approx 2.3026. \]

\[ \zeta = \frac{-\ln M_p^\star}{\sqrt{\pi^2 + (\ln M_p^\star)^2}} \approx \frac{2.3026}{\sqrt{9.8696 + 5.3019}} \approx 0.591. \]

With \( T_s^\star = 2 \),

\[ \omega_n \approx \frac{4}{\zeta T_s^\star} \approx \frac{4}{0.591 \cdot 2} \approx 3.386 \;\mathrm{rad/s}. \]

5.2. Choose Third Pole and Compute PID Gains

Choose \( \alpha = 4 \), so

\[ p_3 = \alpha \zeta\omega_n \approx 4 \cdot 0.591 \cdot 3.386 \approx 8.0. \]

Now compute the desired PID parameters:

\[ \begin{aligned} K_d &= 2\zeta\omega_n + p_3 - 1, \\ K_p &= \omega_n^2 + 2\zeta\omega_n p_3, \\ K_i &= \omega_n^2 p_3. \end{aligned} \]

Numerically,

\[ \omega_n^2 \approx 11.47, \quad 2\zeta\omega_n \approx 4.0. \]

\[ \begin{aligned} K_d &\approx 4.0 + 8.0 - 1 = 11.0, \\ K_p &\approx 11.47 + 2 \cdot 0.591 \cdot 3.386 \cdot 8.0 \approx 11.47 + 32.0 \approx 43.47, \\ K_i &\approx 11.47 \cdot 8.0 \approx 91.76. \end{aligned} \]

In practice, one might round these to \( K_d = 11 \), \( K_p = 43.5 \), \( K_i = 92 \) and then validate and refine using simulation or root-locus-based fine-tuning.

6. Python Implementation for Time-Domain PID Design

We now implement the full design workflow in Python using the control library (widely used in robotics and control applications). The code:

  1. Computes \( \zeta, \omega_n, p_3 \) from specs,
  2. Computes \( K_p,K_i,K_d \),
  3. Builds the closed-loop system and plots the step response.

import numpy as np
import control as ct
import matplotlib.pyplot as plt

# --- Specifications ---
Mp_star = 0.10      # 10% overshoot (relative)
Ts_star = 2.0       # 2% settling time [s]
alpha   = 4.0       # third pole factor p3 = alpha * zeta * omega_n

# --- Step 1: Compute zeta from Mp* ---
lnMp = np.log(Mp_star)
zeta = -lnMp / np.sqrt(np.pi**2 + lnMp**2)

# --- Step 2: Compute omega_n from Ts* ---
omega_n = 4.0 / (zeta * Ts_star)

# --- Step 3: Choose p3 and compute PID gains ---
p3 = alpha * zeta * omega_n

Kd = 2*zeta*omega_n + p3 - 1.0
Kp = omega_n**2 + 2*zeta*omega_n*p3
Ki = omega_n**2 * p3

print(f"Computed gains: Kp={Kp:.3f}, Ki={Ki:.3f}, Kd={Kd:.3f}")

# --- Step 4: Build plant and controller ---
# Plant G(s) = 1/(s*(s+1))
s = ct.TransferFunction.s
G = 1 / (s*(s + 1))

# PID controller in parallel form: C(s) = Kp + Ki/s + Kd*s
C = Kp + Ki/s + Kd*s

# Closed-loop system under unity feedback
T = ct.feedback(C*G, 1)

# --- Step 5: Simulate step response ---
t = np.linspace(0, 10, 1000)
t_out, y_out = ct.step_response(T, t)

# Compute overshoot and settling time empirically
y_final = y_out[-1]
Mp_actual = (np.max(y_out) - y_final) / y_final
# 2% settling time
within = np.where(np.abs(y_out - y_final) < 0.02 * np.abs(y_final))[0]
Ts_actual = t_out[within[0]] if len(within) > 0 else np.nan

print(f"Approx. actual Mp={Mp_actual*100:.1f}%")
print(f"Approx. actual Ts={Ts_actual:.2f} s")

plt.figure()
plt.plot(t_out, y_out)
plt.axhline(y_final, linestyle="--")
plt.xlabel("t [s]")
plt.ylabel("y(t)")
plt.title("Closed-loop Step Response with Designed PID")
plt.grid(True)
plt.show()
      

In a robotics context, similar computations are used when designing joint or actuator controllers for specified settling time and overshoot, before embedding the resulting PID gains into lower-level firmware or ROS-based control nodes.

7. C++ Implementation: Discrete-Time PID and Simulation

In embedded control for robotics, PID is almost always implemented in discrete time. A simple forward-Euler discretization of the PID law is

\[ \begin{aligned} e[k] &= r[k] - y[k], \\ I[k] &= I[k-1] + e[k]\Delta t, \\ D[k] &= \frac{e[k] - e[k-1]}{\Delta t}, \\ u[k] &= K_p e[k] + K_i I[k] + K_d D[k], \end{aligned} \]

where \( \Delta t \) is the sampling period. For simulation, we can also integrate the plant dynamics numerically. Below is a minimal C++ example (without any external control libraries) that:

  • uses pre-computed \( K_p,K_i,K_d \),
  • approximates the plant \( G(s) = 1/(s(s+1)) \) via a second-order ODE,
  • applies discrete-time PID control.

#include <iostream>
#include <vector>

int main() {
    // PID gains from the design example
    double Kp = 43.47;
    double Ki = 91.76;
    double Kd = 11.0;

    // Simulation parameters
    double dt = 0.001;
    double T_end = 10.0;
    int N = static_cast<int>(T_end / dt);

    // Reference (unit step)
    auto r = [](double /*t*/) { return 1.0; };

    // State of plant G(s) = 1/(s*(s+1))
    // State-space realization:
    // x1_dot = x2
    // x2_dot = -x2 - 0*x1 + u
    double x1 = 0.0;  // output position y
    double x2 = 0.0;  // velocity
    double y  = 0.0;

    // PID states
    double e_prev = 0.0;
    double I      = 0.0;

    std::vector<double> t_vec, y_vec;

    for (int k = 0; k < N; ++k) {
        double t = k * dt;
        double ref = r(t);
        y = x1;

        double e = ref - y;
        I += e * dt;
        double D = (e - e_prev) / dt;
        double u = Kp * e + Ki * I + Kd * D;
        e_prev = e;

        // Simple Euler integration of plant
        double x1_dot = x2;
        double x2_dot = -x2 + u; // corresponds roughly to 1/(s*(s+1))
        x1 += x1_dot * dt;
        x2 += x2_dot * dt;

        t_vec.push_back(t);
        y_vec.push_back(y);
    }

    // Print a few samples
    for (size_t i = 0; i < t_vec.size(); i += 1000) {
        std::cout << "t=" << t_vec[i]
                  << " y=" << y_vec[i] << std::endl;
    }

    return 0;
}
      

In a real robotic actuator controller, the above logic would be embedded inside a real-time loop, often combined with sensor filtering and safety checks (e.g., saturation, anti-windup, which will be treated in later lessons).

8. Java Implementation for PID in a Robotics Context

Java-based robotics frameworks (e.g., many educational robot platforms) often expose PID controllers as reusable classes. Below is a simple continuous-time inspired, discrete implementation of a PID class, which you could integrate with a simulation loop or real sensors/actuators.


public class PIDController {
    private double Kp;
    private double Ki;
    private double Kd;

    private double integral;
    private double prevError;
    private boolean firstCall;

    public PIDController(double Kp, double Ki, double Kd) {
        this.Kp = Kp;
        this.Ki = Ki;
        this.Kd = Kd;
        this.integral = 0.0;
        this.prevError = 0.0;
        this.firstCall = true;
    }

    public double update(double ref, double y, double dt) {
        double e = ref - y;
        integral += e * dt;

        double derivative = 0.0;
        if (!firstCall) {
            derivative = (e - prevError) / dt;
        } else {
            firstCall = false;
        }
        prevError = e;

        return Kp * e + Ki * integral + Kd * derivative;
    }

    public void reset() {
        integral = 0.0;
        prevError = 0.0;
        firstCall = true;
    }
}

// Example usage in a main program (simulation skeleton)
public class Example {
    public static void main(String[] args) {
        double Kp = 43.47;
        double Ki = 91.76;
        double Kd = 11.0;

        PIDController pid = new PIDController(Kp, Ki, Kd);

        double dt = 0.001;
        double T_end = 10.0;
        int N = (int)(T_end / dt);

        double x1 = 0.0;
        double x2 = 0.0;
        double y  = 0.0;

        for (int k = 0; k < N; ++k) {
            double t = k * dt;
            double ref = 1.0; // step

            double u = pid.update(ref, y, dt);

            // simple plant approximation: G(s) ~ 1/(s*(s+1))
            double x1_dot = x2;
            double x2_dot = -x2 + u;
            x1 += x1_dot * dt;
            x2 += x2_dot * dt;
            y = x1;

            if (k % 1000 == 0) {
                System.out.println("t=" + t + " y=" + y);
            }
        }
    }
}
      

The gains are obtained analytically as in Sections 2–5, while the Java code focuses on reliable real-time implementation within a robotics or mechatronics software stack.

9. MATLAB/Simulink Implementation

MATLAB and Simulink are standard tools for control engineering and robotics. We first compute the gains symbolically or numerically in MATLAB, then build a Simulink model using basic blocks.


% Specifications
Mp_star = 0.10;
Ts_star = 2.0;
alpha   = 4.0;

% Step 1: zeta from Mp*
lnMp = log(Mp_star);
zeta = -lnMp / sqrt(pi^2 + lnMp^2);

% Step 2: omega_n from Ts*
omega_n = 4 / (zeta * Ts_star);

% Step 3: choose p3 and compute gains
p3 = alpha * zeta * omega_n;
Kd = 2*zeta*omega_n + p3 - 1;
Kp = omega_n^2 + 2*zeta*omega_n*p3;
Ki = omega_n^2 * p3;

fprintf('Kp = %.3f, Ki = %.3f, Kd = %.3f\n', Kp, Ki, Kd);

% Plant G(s) = 1/(s*(s+1))
s = tf('s');
G = 1/(s*(s+1));

% PID controller
C = Kp + Ki/s + Kd*s;

% Closed-loop transfer function
T = feedback(C*G, 1);

% Step response and time-domain characteristics
figure;
step(T);
grid on;
title('Closed-loop Step Response with Designed PID');

info = stepinfo(T);
disp(info);
      

To create an equivalent Simulink model programmatically:


% Create new Simulink model
modelName = 'PID_from_time_specs';
new_system(modelName);
open_system(modelName);

% Add blocks: Step, Sum, PID, Transfer Fcn, Scope
add_block('simulink/Sources/Step', [modelName '/Step']);
add_block('simulink/Math Operations/Sum', [modelName '/Sum'], ...
          'Inputs', '+-');
add_block('simulink/Continuous/PID Controller', [modelName '/PID']);
add_block('simulink/Continuous/Transfer Fcn', [modelName '/Plant']);
add_block('simulink/Sinks/Scope', [modelName '/Scope']);

% Set PID parameters
set_param([modelName '/PID'], 'P', num2str(Kp), ...
                             'I', num2str(Ki), ...
                             'D', num2str(Kd));

% Set plant transfer function 1/(s*(s+1)) => num=1, den=[1 1 0]
set_param([modelName '/Plant'], 'Numerator', '[1]', ...
                                 'Denominator', '[1 1 0]');

% Connect blocks (positions omitted for brevity)
add_line(modelName, 'Step/1', 'Sum/1');
add_line(modelName, 'Sum/1',  'PID/1');
add_line(modelName, 'PID/1',  'Plant/1');
add_line(modelName, 'Plant/1','Scope/1');
add_line(modelName, 'Plant/1','Sum/2');

save_system(modelName);
      

Running the Simulink model yields a closed-loop step response that should closely match the time-domain design targets.

10. Wolfram Mathematica Implementation

Wolfram Mathematica offers symbolic and numeric tools for control design. The following notebook-style code demonstrates the same design workflow.


(* Specifications *)
MpStar  = 0.10;
TsStar  = 2.0;
alpha   = 4.0;

(* Step 1: zeta from Mp* *)
lnMp    = Log[MpStar];
zeta    = -lnMp / Sqrt[Pi^2 + lnMp^2] // N;

(* Step 2: omega_n from Ts* *)
omegaN  = 4.0/(zeta*TsStar);

(* Step 3: p3 and PID gains *)
p3      = alpha*zeta*omegaN;
Kd      = 2*zeta*omegaN + p3 - 1;
Kp      = omegaN^2 + 2*zeta*omegaN*p3;
Ki      = omegaN^2 * p3;

{Kp, Ki, Kd} // N

(* Plant and controller *)
s  = LaplaceTransformVariable;
G  = TransferFunctionModel[1/(s*(s + 1)), s];
C  = TransferFunctionModel[Kp + Ki/s + Kd*s, s];

T  = SystemsModelFeedbackConnect[C*G];

(* Step response *)
step = OutputResponse[T, UnitStep[t], {t, 0, 10}];
Plot[step, {t, 0, 10}, GridLines -> Automatic,
     AxesLabel -> {"t", "y(t)"},
     PlotLabel -> "Closed-loop Step Response"];

(* Time-domain metrics: overshoot and settling-time estimates *)
yFinal = step /. t -> 10.0;
yMax   = Maximize[{step, 0 <= t <= 10}, t][[1]];
MpActual = (yMax - yFinal)/yFinal;
MpActual

(* Numerical settling time using 2% band *)
tsol = t /. FindFirst[
   Table[{Abs[step /. t -> tau - yFinal], tau}, {tau, 0, 10, 0.01}],
   #[[1]] < 0.02*Abs[yFinal] &][[2]];
tsol
      

Mathematica is particularly convenient if you want to symbolically derive PID gains for more general plant structures while still enforcing time-domain pole-placement constraints.

11. Problems and Solutions

Problem 1 (Deriving \( \zeta \) from \( M_p \)): Starting from the overshoot formula for an underdamped second-order system, \( M_p = \mathrm{e}^{-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}} \), derive the explicit expression \( \zeta = \frac{-\ln M_p}{\sqrt{\pi^2 + (\ln M_p)^2}} \).

Solution:

Take natural logarithms on both sides:

\[ \ln M_p = -\frac{\pi\zeta}{\sqrt{1-\zeta^2}}. \]

Multiply both sides by \( -1 \):

\[ -\ln M_p = \frac{\pi\zeta}{\sqrt{1-\zeta^2}}. \]

Square both sides:

\[ (-\ln M_p)^2 = \frac{\pi^2 \zeta^2}{1-\zeta^2}. \]

Rearranging,

\[ (-\ln M_p)^2(1-\zeta^2) = \pi^2\zeta^2 \;\Longrightarrow\; (-\ln M_p)^2 = \zeta^2\bigl[\pi^2 + (-\ln M_p)^2\bigr]. \]

Therefore

\[ \zeta^2 = \frac{(-\ln M_p)^2}{\pi^2 + (-\ln M_p)^2} \quad\Rightarrow\quad \zeta = \frac{-\ln M_p}{\sqrt{\pi^2 + (\ln M_p)^2}}. \]

This matches the formula used in the design algorithm.

Problem 2 (Polynomial Matching for a Different Plant): Let the plant be \( G(s) = \frac{1}{s(s+2)} \) and the controller still be \( C(s) = K_p + \frac{K_i}{s} + K_d s \). Derive the closed-loop characteristic polynomial under unity feedback and identify the coefficient relations between \( K_p,K_i,K_d \) and a desired polynomial \( s^3 + a_2 s^2 + a_1 s + a_0 \).

Solution:

First compute the open-loop transfer function:

\[ L(s) = C(s)G(s) = \left(K_p + \frac{K_i}{s} + K_d s\right)\frac{1}{s(s+2)} = \frac{K_d s^2 + K_p s + K_i}{s^2(s+2)}. \]

The closed-loop characteristic equation is \( 1 + L(s) = 0 \), which yields

\[ s^2(s+2) + K_d s^2 + K_p s + K_i = 0. \]

Expand the plant term:

\[ s^2(s+2) = s^3 + 2s^2, \]

so

\[ s^3 + 2s^2 + K_d s^2 + K_p s + K_i = 0 \;\Longrightarrow\; s^3 + (2 + K_d) s^2 + K_p s + K_i = 0. \]

Matching with the desired polynomial \( s^3 + a_2 s^2 + a_1 s + a_0 \) yields

\[ \boxed{ \begin{aligned} 2 + K_d &= a_2, \\ K_p &= a_1, \\ K_i &= a_0. \end{aligned} } \]

Once \( a_2, a_1, a_0 \) are chosen based on desired dominant poles and third pole, the PID gains follow immediately.

Problem 3 (Third Pole Placement Heuristic): Consider a desired dominant pole pair with real part \( -\sigma_d = -\zeta\omega_n \). Show that if the third pole is located at \( -p_3 \) with \( p_3 \geq 5\sigma_d \), then the contribution of the third mode to the step response decays at least \( \mathrm{e}^{-5} \) times faster than the envelope of the dominant mode.

Solution:

The dominant envelope of the underdamped second-order response decays approximately as \( \mathrm{e}^{-\sigma_d t} \), whereas the third mode decays as \( \mathrm{e}^{-p_3 t} \). At time \( t \),

\[ \frac{\mathrm{e}^{-p_3 t}}{\mathrm{e}^{-\sigma_d t}} = \mathrm{e}^{-(p_3 - \sigma_d) t}. \]

If \( p_3 \geq 5\sigma_d \), then \( p_3 - \sigma_d \geq 4\sigma_d \) and for \( t \approx 1/\sigma_d \) (a characteristic time of the dominant mode),

\[ \mathrm{e}^{-(p_3 - \sigma_d)t} \leq \mathrm{e}^{-4\sigma_d (1/\sigma_d)} = \mathrm{e}^{-4}. \]

If we increase to \( p_3 \geq 6\sigma_d \), we obtain \( \mathrm{e}^{-5} \). Thus, choosing \( p_3 \geq 5\sigma_d \) ensures the third mode decays several e-foldings faster than the dominant envelope, justifying the second-order approximation.

Problem 4 (Effect of Integral Action on Steady-State Error): For the plant \( G(s) = 1/(s(s+1)) \) and PID controller \( C(s) = K_p + K_i/s + K_d s \), show that the closed-loop system has zero steady-state error to a unit step for any positive \( K_i \), assuming closed-loop stability.

Solution:

The closed-loop transfer function is \( T(s) = \frac{C(s)G(s)}{1 + C(s)G(s)} \). The steady-state error to a unit step is

\[ e_{\mathrm{ss}} = \lim_{t \to \infty} e(t) = \lim_{s \to 0} s\left(\frac{1}{s} - T(s)\frac{1}{s}\right) = 1 - \lim_{s \to 0} T(s). \]

But

\[ T(s) = \frac{C(s)G(s)}{1 + C(s)G(s)} = \frac{K_d s^2 + K_p s + K_i}{s^2(s+1) + (K_d s^2 + K_p s + K_i)}. \]

Evaluate at \( s \to 0 \):

\[ \lim_{s \to 0} T(s) = \frac{K_i}{0 + K_i} = 1 \]

whenever \( K_i \neq 0 \). Therefore

\[ e_{\mathrm{ss}} = 1 - 1 = 0. \]

Thus, any PID controller with positive integral gain \( K_i > 0 \) yields zero steady-state error to a unit step, provided the closed-loop is stable.

Problem 5 (Design Check for Different Specifications): For the same plant \( G(s) = 1/(s(s+1)) \), suppose we desire \( M_p^\star = 0.05 \) (5%) and \( T_s^\star = 3\,\mathrm{s} \). Using the design formulas in this lesson and \( \alpha = 4 \), compute approximate PID gains. (You may leave intermediate values in symbolic form, but give numerical values for \( K_p,K_i,K_d \).)

Solution (outline):

  1. Compute \( \zeta \) from \( M_p^\star = 0.05 \) using the formula derived in Problem 1.
  2. Compute \( \omega_n = 4/(\zeta T_s^\star) \).
  3. Set \( p_3 = 4\zeta\omega_n \).
  4. Use \( K_d = 2\zeta\omega_n + p_3 - 1 \), \( K_p = \omega_n^2 + 2\zeta\omega_n p_3 \), \( K_i = \omega_n^2 p_3 \).

Numerically, you obtain slightly larger \( \zeta \) and smaller \( \omega_n \) than in the worked example, resulting in smaller overshoot and slower response. (Students are encouraged to perform actual computations in Python or MATLAB as an exercise.)

12. Summary

In this lesson, we connected classical time-domain specifications (maximum overshoot, settling time) to the parameters of a standard second-order prototype and then used polynomial matching to design PID gains for a simple but representative plant. The key steps are:

  • Translate \( M_p^\star, T_s^\star \) into \( \zeta, \omega_n \) and desired dominant poles.
  • Augment with a suitably fast third pole \( p_3 \) to accommodate the higher-order closed-loop.
  • Match the closed-loop characteristic polynomial to obtain explicit formulas for \( K_p, K_i, K_d \).
  • Implement the design in multiple software environments and verify the time response numerically.

This time-domain PID design approach is a bridge between intuitive specifications (settling time, overshoot) and more systematic design methods. Subsequent lessons will incorporate frequency-domain viewpoints and address robustness, anti-windup, and noise-aware derivative action.

13. References

  1. Åström, K.J., & Hägglund, T. (1995). PID Controllers: Theory, Design, and Tuning. ISA (monograph; various journal papers derived from this work).
  2. Ho, W.K., Datta, A., & Bhattacharyya, S.P. (1997). Analytical expressions for PID controller design satisfying mixed time and frequency domain specifications. IEEE Transactions on Automatic Control, 42(6), 815–819.
  3. Åström, K.J., & Hägglund, T. (2000). Benchmark systems for PID control. IFAC Journal of Control Engineering Practice, 8(6), 673–682.
  4. Hang, C.C., Åström, K.J., & Ho, W.K. (1991). Refinements of the Ziegler–Nichols tuning formula. IEE Proceedings D (Control Theory and Applications), 138(2), 111–118.
  5. Hägglund, T. (1999). An industrial dead-time compensating PI controller. Control Engineering Practice, 7(4), 389–399.
  6. Ho, W.K., Hang, C.C., & Cao, L.S. (1994). Tuning of PID controllers based on gain and phase margin specifications. Automatica, 30(3), 497–502.
  7. Silva, G.J., Datta, A., & Bhattacharyya, S.P. (2002). New results on the synthesis of PID controllers. IEEE Transactions on Automatic Control, 47(2), 241–252.
  8. O’Dwyer, A. (2003). PI and PID controller tuning rules: an overview and personal perspective. IFAC Conference on Digital Control: Past, Present and Future of PID Control, 1–24.
  9. Middleton, R.H., & Goodwin, G.C. (1988). Improved finite word length characteristics in digital control using delta operators. IEEE Transactions on Automatic Control, 33(10), 1046–1058.
  10. Visioli, A. (2001). Optimal tuning of PID controllers for integral and unstable processes. IEE Proceedings – Control Theory and Applications, 148(2), 180–188.