In state-feedback design, the implementable design variable is the gain matrix \( K \) in \( u=-Kx+r \). For a single-input LTI system

\[ \dot{x}=Ax+bu,\qquad u=-Kx+r,\qquad \dot{x}=(A-bK)x+br. \]

The free closed-loop motion is determined by the eigenvalues of \( A-bK \). Therefore specifications such as “less than 10% overshoot” or “settle within 2 seconds” must first be mapped into a desired characteristic polynomial \( \Delta_d(s) \) or desired closed-loop poles \( \{p_1,\dots,p_n\} \).

\[ \det(sI-(A-bK))=\Delta_d(s)=\prod_{i=1}^{n}(s-p_i). \]

Lesson 1 formulated this assignment problem, Lesson 2 solved it in controllable canonical form, and Lesson 3 introduced Ackermann’s formula. The present lesson answers the design question: which poles should we ask Ackermann’s formula to assign?

flowchart TD
  A["Time-domain requirements"] --> B["Choose dominant prototype"]
  B --> C["Compute damping ratio and natural frequency"]
  C --> D["Form dominant complex pole pair"]
  D --> E["Add remaining non-dominant poles"]
  E --> F["Build desired characteristic polynomial"]
  F --> G["Compute state-feedback gain K"]
  G --> H["Check closed-loop response and control effort"]
        

Most classical transient specifications are easiest to interpret through the normalized second-order closed-loop model

\[ G_2(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_n s+ \omega_n^2},\qquad 0<\zeta<1. \]

Its poles are

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

Define the exponential decay rate and damped natural frequency by

\[ \sigma=\zeta\omega_n, \qquad \omega_d=\omega_n\sqrt{1-\zeta^2}, \qquad p_{1,2}=-\sigma\pm j\omega_d. \]

The real part controls decay speed, while the imaginary part controls oscillation. If the full system has order \( n>2 \), we often select this pair as the dominant pair and place the remaining poles farther left so that they decay faster and contribute less to the visible transient.

For the underdamped second-order unit-step response, the maximum overshoot fraction \( M_p \) satisfies

\[ M_p=e^{-\frac{\pi\zeta}{\sqrt{1-\zeta^2} } }, \qquad 0<M_p<1. \]

If overshoot is specified as a percentage \( PO=100M_p \), solve the preceding equation for \( \zeta \):

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

Squaring both sides and isolating \( \zeta^2 \) gives

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

Thus a smaller allowed overshoot corresponds to a larger damping ratio, which geometrically rotates the desired poles closer to the negative real axis.

The transient envelope of the second-order response decays essentially like \( e^{-\zeta\omega_n t} \). For the 2% settling criterion, the standard engineering approximation is

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

For the 5% criterion, a common approximation is

\[ T_s\approx \frac{3}{\zeta\omega_n}=\frac{3}{\sigma}. \]

Therefore, if a 2% settling time specification says \( T_s\le T_s^{\max} \), then the dominant real part should approximately satisfy

\[ \sigma\ge \frac{4}{T_s^{\max} }, \qquad \operatorname{Re}(p_{1,2})=-\sigma\le -\frac{4}{T_s^{\max} }. \]

Faster settling requires more negative real parts, but this usually increases feedback gain magnitude and control effort.

The first peak occurs at approximately

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

Hence a peak-time requirement fixes the imaginary part:

\[ \omega_d=\frac{\pi}{T_p}, \qquad p_{1,2}=-\frac{\zeta\omega_d}{\sqrt{1-\zeta^2} } \pm j\omega_d. \]

Rise-time formulae depend on the definition of rise time. A frequently used rough approximation for underdamped systems is

\[ T_r\approx \frac{1.8}{\omega_n}. \]

Since pole placement is an exact eigenvalue assignment method but the classical time-domain formulae are approximations, the design process is normally iterative: map specifications to initial poles, compute \( K \), simulate the closed-loop state-space system, and adjust if the high-order dynamics or zeros distort the response.

Time-domain requirements can be interpreted as regions in the \( s \)-plane. A minimum damping ratio places poles inside a cone around the negative real axis; a maximum settling time places poles to the left of a vertical line; a limit on oscillation speed constrains the imaginary part.

\[ \zeta=\frac{-\operatorname{Re}(p)}{|p|}, \qquad |\operatorname{Im}(p)|=\omega_d, \qquad -\operatorname{Re}(p)=\sigma. \]

\[ \zeta\ge \zeta_{\min},\qquad \operatorname{Re}(p)\le -\frac{4}{T_s^{\max} }, \qquad |\operatorname{Im}(p)|\le \omega_d^{\max}. \]

flowchart TD
  A["Candidate pole p"] --> B["Check damping: \nzeta above limit"]
  A --> C["Check settling: \nreal part far left enough"]
  A --> D["Check oscillation: \nimaginary part acceptable"]
  B --> E["Accept feasible pole region"]
  C --> E
  D --> E
  E --> F["Select dominant pair and extra poles"]
        

A pole is a good candidate only if it satisfies all active constraints. However, excessively fast poles may demand large actuator effort and can amplify sensor noise or unmodeled dynamics.

For an \( n \)-state single-input system with a desired dominant pair \( p_{1,2}=-\sigma_d\pm j\omega_d \), the remaining \( n-2 \) poles are often chosen real and farther left:

\[ p_3,p_4,\dots,p_n\in\mathbb{R}, \qquad \operatorname{Re}(p_k)\le -\eta\sigma_d, \qquad \eta\in[4,10]. \]

The larger \( \eta \) is, the faster the non-dominant modes decay. But a very large \( \eta \) can produce a large feedback gain \( K \), which may be unacceptable under actuator saturation or model uncertainty.

\[ \Delta_d(s)=\left(s^2+2\zeta\omega_n s+ \omega_n^2\right)\prod_{k=3}^{n}(s-p_k). \]

\[ K=e_n^T\mathcal{C}^{-1}\Delta_d(A), \qquad \mathcal{C}=\begin{bmatrix}b&Ab&\cdots&A^{n-1}b\end{bmatrix}. \]

Suppose a controllable third-order single-input system must satisfy \( PO=10\% \) and 2% settling time \( T_s\le 2 \) seconds. The overshoot fraction is \( M_p=0.10 \), so

\[ \zeta=\frac{-\ln(0.10)}{\sqrt{\pi^2+(\ln(0.10))^2} } \approx 0.5912. \]

\[ \sigma=\frac{4}{T_s}=2, \qquad \omega_n=\frac{\sigma}{\zeta}\approx 3.3832. \]

\[ \omega_d=\omega_n\sqrt{1-\zeta^2}\approx 2.7288, \qquad p_{1,2}=-2\pm j2.7288. \]

For a third-order system, select one extra non-dominant pole, for example \( p_3=-12 \). Then

\[ \Delta_d(s)=(s^2+4s+11.4461)(s+12) =s^3+16s^2+59.4461s+137.3531. \]

For the plant

\[ A=\begin{bmatrix}0&1&0\\0&0&1\\-2&-3&-1\end{bmatrix}, \qquad b=\begin{bmatrix}0\\0\\1\end{bmatrix}, \]

Ackermann’s formula gives approximately

\[ K=\begin{bmatrix}135.3531&56.4461&15.0000\end{bmatrix}. \]

The following implementations perform the same workflow: convert overshoot and settling-time requirements into desired poles, append a faster real pole for a third-order example, and compute \( K \) by Ackermann’s formula. Python uses NumPy; MATLAB optionally compares with the Control System Toolbox function acker; C++ and Java implement the required linear-algebra steps from scratch for teaching transparency.

Chapter23_Lesson4.py

"""
Chapter23_Lesson4.py
Mapping time-domain specifications to desired closed-loop poles for SISO pole placement.
Self-contained Ackermann implementation using NumPy.
"""

import numpy as np


def damping_ratio_from_overshoot(percent_overshoot: float) -> float:
    """Return zeta from percent overshoot Mp in percent, valid for 0 < Mp < 100."""
    if not (0.0 < percent_overshoot < 100.0):
        raise ValueError("percent_overshoot must be between 0 and 100")
    m = percent_overshoot / 100.0
    L = np.log(m)
    return -L / np.sqrt(np.pi**2 + L**2)


def desired_poles_from_specs(percent_overshoot: float, settling_time: float, settling_fraction: float = 0.02):
    """
    Compute dominant second-order pole pair from percent overshoot and settling time.
    settling_fraction=0.02 uses the 2 percent criterion, approximately Ts = 4/(zeta*wn).
    settling_fraction=0.05 uses the 5 percent criterion, approximately Ts = 3/(zeta*wn).
    """
    zeta = damping_ratio_from_overshoot(percent_overshoot)
    c = 4.0 if abs(settling_fraction - 0.02) < 1e-12 else 3.0
    omega_n = c / (zeta * settling_time)
    sigma = zeta * omega_n
    omega_d = omega_n * np.sqrt(max(0.0, 1.0 - zeta**2))
    return zeta, omega_n, np.array([-sigma + 1j * omega_d, -sigma - 1j * omega_d], dtype=complex)


def companion_extra_poles(dominant_poles, n: int, multiplier: float = 5.0):
    """Place remaining real poles multiplier times farther left than dominant real part."""
    if n < len(dominant_poles):
        raise ValueError("system order n must be at least the number of dominant poles")
    poles = list(dominant_poles)
    sigma = min(-p.real for p in dominant_poles if p.real < 0)
    for i in range(n - len(dominant_poles)):
        poles.append(-(multiplier + i) * sigma)
    return np.array(poles, dtype=complex)


def controllability_matrix(A, b):
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float).reshape(-1, 1)
    n = A.shape[0]
    return np.hstack([np.linalg.matrix_power(A, k) @ b for k in range(n)])


def ackermann_gain(A, b, poles):
    """Compute K such that eig(A - bK) equals poles for a controllable SISO pair."""
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float).reshape(-1, 1)
    n = A.shape[0]
    Wc = controllability_matrix(A, b)
    if np.linalg.matrix_rank(Wc) != n:
        raise ValueError("(A,b) is not controllable")

    coeff = np.poly(poles)  # s^n + a_{n-1}s^{n-1}+...+a0
    phi_A = np.linalg.matrix_power(A, n)
    for i in range(n):
        power = n - 1 - i
        phi_A += coeff[i + 1] * np.linalg.matrix_power(A, power)

    e_nT = np.zeros((1, n))
    e_nT[0, -1] = 1.0
    K = e_nT @ np.linalg.inv(Wc) @ phi_A
    return np.real_if_close(K)


if __name__ == "__main__":
    # Example: third-order SISO plant. Design dominant poles from Mp and Ts.
    A = np.array([[0.0, 1.0, 0.0],
                  [0.0, 0.0, 1.0],
                  [-2.0, -3.0, -1.0]])
    b = np.array([[0.0], [0.0], [1.0]])

    Mp = 10.0       # percent overshoot
    Ts = 2.0        # seconds, 2 percent settling time
    zeta, omega_n, pair = desired_poles_from_specs(Mp, Ts)
    desired = companion_extra_poles(pair, n=A.shape[0], multiplier=6.0)
    K = ackermann_gain(A, b, desired)

    print("zeta =", zeta)
    print("omega_n =", omega_n)
    print("desired poles =", desired)
    print("K =", K)
    print("closed-loop eigenvalues =", np.linalg.eigvals(A - b @ K))

Chapter23_Lesson4.cpp

// Chapter23_Lesson4.cpp
// Mapping time-domain specifications to desired poles and Ackermann pole placement.
// Compile: g++ -std=c++17 Chapter23_Lesson4.cpp -O2 -o Chapter23_Lesson4

#include <cmath>
#include <complex>
#include <iostream>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
using Vec = std::vector<double>;
using Cx = std::complex<double>;

Matrix identity(int n) {
    Matrix I(n, Vec(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}

Matrix matmul(const Matrix& A, const Matrix& B) {
    int n = (int)A.size(), m = (int)B[0].size(), r = (int)B.size();
    Matrix C(n, Vec(m, 0.0));
    for (int i = 0; i < n; ++i)
        for (int k = 0; k < r; ++k)
            for (int j = 0; j < m; ++j)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Vec matvec(const Matrix& A, const Vec& x) {
    Vec y(A.size(), 0.0);
    for (size_t i = 0; i < A.size(); ++i)
        for (size_t j = 0; j < x.size(); ++j) y[i] += A[i][j] * x[j];
    return y;
}

Matrix matadd(const Matrix& A, const Matrix& B, double scaleB = 1.0) {
    Matrix C = A;
    for (size_t i = 0; i < A.size(); ++i)
        for (size_t j = 0; j < A[0].size(); ++j) C[i][j] += scaleB * B[i][j];
    return C;
}

Matrix matscale(const Matrix& A, double s) {
    Matrix B = A;
    for (auto& row : B) for (double& v : row) v *= s;
    return B;
}

Matrix matpow(Matrix A, int p) {
    int n = (int)A.size();
    Matrix R = identity(n);
    while (p > 0) {
        if (p & 1) R = matmul(R, A);
        A = matmul(A, A);
        p >>= 1;
    }
    return R;
}

Matrix inverse(Matrix A) {
    int n = (int)A.size();
    Matrix aug(n, Vec(2 * n, 0.0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) aug[i][j] = A[i][j];
        aug[i][n + i] = 1.0;
    }
    for (int col = 0; col < n; ++col) {
        int piv = col;
        for (int r = col + 1; r < n; ++r)
            if (std::abs(aug[r][col]) > std::abs(aug[piv][col])) piv = r;
        if (std::abs(aug[piv][col]) < 1e-12) throw std::runtime_error("singular matrix");
        std::swap(aug[piv], aug[col]);
        double div = aug[col][col];
        for (double& v : aug[col]) v /= div;
        for (int r = 0; r < n; ++r) if (r != col) {
            double f = aug[r][col];
            for (int j = 0; j < 2 * n; ++j) aug[r][j] -= f * aug[col][j];
        }
    }
    Matrix inv(n, Vec(n, 0.0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j) inv[i][j] = aug[i][n + j];
    return inv;
}

std::vector<Cx> polyFromRoots(const std::vector<Cx>& roots) {
    std::vector<Cx> c = {Cx(1.0, 0.0)};
    for (Cx r : roots) {
        std::vector<Cx> next(c.size() + 1, Cx(0.0, 0.0));
        for (size_t i = 0; i < c.size(); ++i) {
            next[i] += c[i];
            next[i + 1] += -r * c[i];
        }
        c = next;
    }
    return c;
}

double dampingRatioFromOvershoot(double percentOvershoot) {
    double m = percentOvershoot / 100.0;
    double L = std::log(m);
    return -L / std::sqrt(M_PI * M_PI + L * L);
}

std::vector<Cx> dominantPoles(double percentOvershoot, double settlingTime) {
    double zeta = dampingRatioFromOvershoot(percentOvershoot);
    double wn = 4.0 / (zeta * settlingTime);
    double sigma = zeta * wn;
    double wd = wn * std::sqrt(std::max(0.0, 1.0 - zeta * zeta));
    return {Cx(-sigma, wd), Cx(-sigma, -wd)};
}

Matrix controllabilityMatrix(const Matrix& A, const Vec& b) {
    int n = (int)A.size();
    Matrix W(n, Vec(n, 0.0));
    Vec col = b;
    for (int k = 0; k < n; ++k) {
        for (int i = 0; i < n; ++i) W[i][k] = col[i];
        col = matvec(A, col);
    }
    return W;
}

Vec ackermannGain(const Matrix& A, const Vec& b, const std::vector<Cx>& poles) {
    int n = (int)A.size();
    Matrix W = controllabilityMatrix(A, b);
    Matrix Winv = inverse(W);
    auto coeffCx = polyFromRoots(poles); // s^n + a_{n-1}s^{n-1}+...+a0

    Matrix phi = matpow(A, n);
    for (int i = 0; i < n; ++i) {
        int power = n - 1 - i;
        phi = matadd(phi, matscale(matpow(A, power), coeffCx[i + 1].real()));
    }
    Matrix temp = matmul(Winv, phi);
    Vec K(n, 0.0); // e_n^T * temp
    for (int j = 0; j < n; ++j) K[j] = temp[n - 1][j];
    return K;
}

int main() {
    Matrix A = { {0, 1, 0}, {0, 0, 1}, {-2, -3, -1} };
    Vec b = {0, 0, 1};

    auto poles = dominantPoles(10.0, 2.0);
    double sigma = -poles[0].real();
    poles.push_back(Cx(-6.0 * sigma, 0.0));

    Vec K = ackermannGain(A, b, poles);
    std::cout << "Desired poles:\n";
    for (auto p : poles) std::cout << p << "\n";
    std::cout << "K = [ ";
    for (double k : K) std::cout << k << " ";
    std::cout << "]\n";
    return 0;
}

Chapter23_Lesson4.java

// Chapter23_Lesson4.java
// Mapping time-domain specifications to desired poles and Ackermann pole placement.
// Compile: javac Chapter23_Lesson4.java && java Chapter23_Lesson4

import java.util.Arrays;

public class Chapter23_Lesson4 {
    static double[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] multiply(double[][] A, double[][] B) {
        int n = A.length, m = B[0].length, r = B.length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int k = 0; k < r; k++)
                for (int j = 0; j < m; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    static double[] multiply(double[][] A, double[] x) {
        double[] y = new double[A.length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < x.length; j++) y[i] += A[i][j] * x[j];
        return y;
    }

    static double[][] add(double[][] A, double[][] B, double scaleB) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++) C[i][j] = A[i][j] + scaleB * B[i][j];
        return C;
    }

    static double[][] scale(double[][] A, double s) {
        double[][] B = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++) B[i][j] = s * A[i][j];
        return B;
    }

    static double[][] power(double[][] A, int p) {
        double[][] R = eye(A.length), B = A;
        while (p > 0) {
            if ((p & 1) == 1) R = multiply(R, B);
            B = multiply(B, B);
            p >>= 1;
        }
        return R;
    }

    static double[][] inverse(double[][] A) {
        int n = A.length;
        double[][] aug = new double[n][2 * n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(A[i], 0, aug[i], 0, n);
            aug[i][n + i] = 1.0;
        }
        for (int col = 0; col < n; col++) {
            int piv = col;
            for (int r = col + 1; r < n; r++)
                if (Math.abs(aug[r][col]) > Math.abs(aug[piv][col])) piv = r;
            if (Math.abs(aug[piv][col]) < 1e-12) throw new RuntimeException("singular matrix");
            double[] tmp = aug[piv]; aug[piv] = aug[col]; aug[col] = tmp;
            double div = aug[col][col];
            for (int j = 0; j < 2 * n; j++) aug[col][j] /= div;
            for (int r = 0; r < n; r++) if (r != col) {
                double f = aug[r][col];
                for (int j = 0; j < 2 * n; j++) aug[r][j] -= f * aug[col][j];
            }
        }
        double[][] inv = new double[n][n];
        for (int i = 0; i < n; i++) System.arraycopy(aug[i], n, inv[i], 0, n);
        return inv;
    }

    static double dampingRatioFromOvershoot(double percentOvershoot) {
        double m = percentOvershoot / 100.0;
        double L = Math.log(m);
        return -L / Math.sqrt(Math.PI * Math.PI + L * L);
    }

    static double[] secondOrderPolynomialFromSpecs(double percentOvershoot, double settlingTime) {
        double zeta = dampingRatioFromOvershoot(percentOvershoot);
        double wn = 4.0 / (zeta * settlingTime);
        return new double[]{1.0, 2.0 * zeta * wn, wn * wn};
    }

    static double[] convolve(double[] a, double[] b) {
        double[] c = new double[a.length + b.length - 1];
        for (int i = 0; i < a.length; i++)
            for (int j = 0; j < b.length; j++) c[i + j] += a[i] * b[j];
        return c;
    }

    static double[][] controllabilityMatrix(double[][] A, double[] b) {
        int n = A.length;
        double[][] W = new double[n][n];
        double[] col = Arrays.copyOf(b, n);
        for (int k = 0; k < n; k++) {
            for (int i = 0; i < n; i++) W[i][k] = col[i];
            col = multiply(A, col);
        }
        return W;
    }

    static double[] ackermannGain(double[][] A, double[] b, double[] desiredPoly) {
        int n = A.length;
        double[][] W = controllabilityMatrix(A, b);
        double[][] Winv = inverse(W);
        double[][] phi = power(A, n);
        for (int i = 0; i < n; i++) {
            int pow = n - 1 - i;
            phi = add(phi, scale(power(A, pow), desiredPoly[i + 1]), 1.0);
        }
        double[][] temp = multiply(Winv, phi);
        return Arrays.copyOf(temp[n - 1], n);
    }

    public static void main(String[] args) {
        double[][] A = { {0, 1, 0}, {0, 0, 1}, {-2, -3, -1} };
        double[] b = {0, 0, 1};

        double Mp = 10.0, Ts = 2.0;
        double zeta = dampingRatioFromOvershoot(Mp);
        double wn = 4.0 / (zeta * Ts);
        double sigma = zeta * wn;

        double[] secondOrder = secondOrderPolynomialFromSpecs(Mp, Ts);
        double[] extraPole = {1.0, 6.0 * sigma}; // factor s + 6*sigma
        double[] desiredPoly = convolve(secondOrder, extraPole);
        double[] K = ackermannGain(A, b, desiredPoly);

        System.out.println("zeta = " + zeta);
        System.out.println("omega_n = " + wn);
        System.out.println("desired polynomial = " + Arrays.toString(desiredPoly));
        System.out.println("K = " + Arrays.toString(K));
    }
}

Chapter23_Lesson4.m

% Chapter23_Lesson4.m
% Mapping time-domain specifications to desired poles for SISO pole placement.
% Includes from-scratch Ackermann implementation and MATLAB Control System Toolbox comparison.

clear; clc;

A = [0 1 0;
     0 0 1;
    -2 -3 -1];
b = [0; 0; 1];

Mp = 10;        % percent overshoot
Ts = 2.0;       % seconds, 2 percent settling time

[zeta, wn, pair] = desired_poles_from_specs(Mp, Ts);
sigma = -real(pair(1));
poles = [pair, -6*sigma];

K_ack = ackermann_gain(A, b, poles);
disp('zeta ='); disp(zeta);
disp('omega_n ='); disp(wn);
disp('desired poles ='); disp(poles);
disp('K from Ackermann ='); disp(K_ack);
disp('closed-loop eigenvalues ='); disp(eig(A - b*K_ack));

% If Control System Toolbox is installed, compare with acker/place:
if exist('acker', 'file') == 2
    K_toolbox = acker(A, b, poles);
    disp('K from acker ='); disp(K_toolbox);
end

function zeta = damping_ratio_from_overshoot(percent_overshoot)
    if percent_overshoot <= 0 || percent_overshoot >= 100
        error('percent_overshoot must be between 0 and 100');
    end
    m = percent_overshoot / 100;
    L = log(m);
    zeta = -L / sqrt(pi^2 + L^2);
end

function [zeta, wn, poles] = desired_poles_from_specs(percent_overshoot, settling_time)
    zeta = damping_ratio_from_overshoot(percent_overshoot);
    wn = 4 / (zeta * settling_time);
    sigma = zeta * wn;
    wd = wn * sqrt(max(0, 1 - zeta^2));
    poles = [-sigma + 1i*wd, -sigma - 1i*wd];
end

function Wc = controllability_matrix(A, b)
    n = size(A, 1);
    Wc = zeros(n, n);
    for k = 1:n
        Wc(:, k) = A^(k-1) * b;
    end
end

function K = ackermann_gain(A, b, poles)
    n = size(A, 1);
    Wc = controllability_matrix(A, b);
    if rank(Wc) < n
        error('(A,b) is not controllable');
    end
    coeff = poly(poles); % [1 a_{n-1} ... a0]
    phiA = A^n;
    for i = 1:n
        power = n - i;
        phiA = phiA + coeff(i+1) * A^power;
    end
    eT = zeros(1, n); eT(n) = 1;
    K = real(eT / Wc * phiA);
end

Chapter23_Lesson4.nb

Notebook[{
Cell["Chapter23_Lesson4.nb", "Title"],
Cell["Mapping time-domain specifications to desired poles for SISO pole placement.", "Text"],
Cell[BoxData["
ClearAll[dampingRatioFromOvershoot, desiredPolesFromSpecs, controllabilityMatrix, ackermannGain];

dampingRatioFromOvershoot[percentOvershoot_] := Module[{m, L},
  If[percentOvershoot <= 0 || percentOvershoot >= 100,
    Print[\"percentOvershoot must be between 0 and 100\"]; Abort[]
  ];
  m = percentOvershoot/100;
  L = Log[m];
  -L/Sqrt[Pi^2 + L^2]
];

desiredPolesFromSpecs[percentOvershoot_, settlingTime_] := Module[{zeta, wn, sigma, wd},
  zeta = dampingRatioFromOvershoot[percentOvershoot];
  wn = 4/(zeta settlingTime);
  sigma = zeta wn;
  wd = wn Sqrt[Max[0, 1 - zeta^2]];
  {zeta, wn, {-sigma + I wd, -sigma - I wd} }
];

controllabilityMatrix[A_, b_] := Module[{n = Length[A]},
  Transpose[Table[MatrixPower[A, k].b, {k, 0, n - 1}]]
];

ackermannGain[A_, b_, poles_] := Module[{n, Wc, coeff, phiA, eT},
  n = Length[A];
  Wc = controllabilityMatrix[A, b];
  If[MatrixRank[Wc] < n, Print[\"(A,b) is not controllable\"]; Abort[]];
  coeff = CoefficientList[Expand[Times @@ (s - # & /@ poles)], s] // Reverse;
  phiA = MatrixPower[A, n] + Sum[coeff[[i + 1]] MatrixPower[A, n - i], {i, 1, n}];
  eT = UnitVector[n, n];
  Chop[eT.Inverse[Wc].phiA]
];

A = { {0, 1, 0}, {0, 0, 1}, {-2, -3, -1} };
b = {0, 0, 1};
{zeta, wn, pair} = desiredPolesFromSpecs[10, 2.0];
sigma = -Re[pair[[1]]];
poles = Join[pair, {-6 sigma}];
K = ackermannGain[A, b, poles];
{zeta, wn, poles, K, Eigenvalues[A - Outer[Times, b, K]]}
"], "Input"]
}]

For a Simulink verification model, use a State-Space block for \( (A,b,C,D) \), a Gain block for \( -K \), a summing junction for \( u=-Kx+r \), and a Scope block for the output. The MATLAB script above provides the numerical gain used in the Simulink feedback loop.

Problem 1 (Overshoot Mapping): A closed-loop system must have no more than 5% overshoot. Compute the minimum damping ratio for a dominant second-order design.

Solution: Here \( M_p=0.05 \). Thus

\[ \zeta=\frac{-\ln(0.05)}{\sqrt{\pi^2+(\ln(0.05))^2} } \approx 0.6901. \]

Therefore choose dominant poles on or inside the damping-ratio cone corresponding to \( \zeta\ge 0.6901 \).

Problem 2 (Settling-Time Line): A 2% settling-time requirement is \( T_s\le 1.5 \) seconds. What real-part bound should the dominant poles satisfy?

Solution: Using \( T_s\approx 4/\sigma \),

\[ \sigma\ge \frac{4}{1.5}=2.6667, \qquad \operatorname{Re}(p)\le -2.6667. \]

Problem 3 (Pole Pair from Two Specifications): Design a dominant pole pair for \( PO=20\% \) and \( T_s=4 \) seconds using the 2% settling criterion.

Solution:

\[ \zeta=\frac{-\ln(0.20)}{\sqrt{\pi^2+(\ln(0.20))^2} } \approx 0.4559, \qquad \sigma=\frac{4}{4}=1. \]

\[ \omega_n=\frac{1}{0.4559}\approx 2.193, \qquad \omega_d\approx 1.952, \qquad p_{1,2}=-1\pm j1.952. \]

Problem 4 (Desired Polynomial): A fourth-order system has desired dominant poles \( -2\pm j3 \). Choose two non-dominant real poles at \( -10 \) and \( -12 \). Find the desired characteristic polynomial.

Solution:

\[ \Delta_d(s)=((s+2)^2+3^2)(s+10)(s+12) =(s^2+4s+13)(s^2+22s+120). \]

\[ \Delta_d(s)=s^4+26s^3+221s^2+766s+1560. \]

Problem 5 (Ackermann Substitution): For a controllable third-order pair \( (A,b) \) and desired polynomial \( \Delta_d(s)=s^3+\alpha_2s^2+\alpha_1s+ \alpha_0 \), write the matrix polynomial used in Ackermann’s formula.

Solution:

\[ \Delta_d(A)=A^3+\alpha_2A^2+\alpha_1A+ \alpha_0I. \]

\[ K=\begin{bmatrix}0&0&1\end{bmatrix} \mathcal{C}^{-1}\Delta_d(A), \qquad \mathcal{C}=\begin{bmatrix}b&Ab&A^2b\end{bmatrix}. \]

Time-domain pole placement begins by converting performance requirements into geometric constraints in the complex plane. Percent overshoot maps primarily to damping ratio; settling time maps primarily to real-part decay rate; peak time maps to damped natural frequency. For higher-order systems, a dominant second-order pair is supplemented by faster non-dominant poles. Once the desired polynomial is constructed, single-input pole placement proceeds exactly through the controllability matrix and Ackermann’s formula.