Consider a translational mass-spring-damper plant with actuator force \( u(t) \), position \( q(t) \), and output \( y(t)=q(t) \):

\[ m\ddot{q}(t)+b\dot{q}(t)+kq(t)=u(t),\qquad m>0,\; b>0,\; k>0. \]

The design target is a full-state feedback controller \( u=-K\mathbf{x}+\bar{N}r \) that makes the position track a constant reference \( r \). The state is selected as

\[ \mathbf{x}=\begin{bmatrix}x_1\\x_2\end{bmatrix} =\begin{bmatrix}q\\\dot{q}\end{bmatrix}. \]

flowchart TD
  A["Physical parameters: m, b, k"] --> B["Choose states: position and velocity"]
  B --> C["Build state-space matrices"]
  C --> D["Check controllability"]
  D --> E["Choose desired closed-loop poles"]
  E --> F["Compute feedback gain K"]
  F --> G["Compute reference gain Nbar"]
  G --> H["Simulate and verify"]
        

From \( x_1=q \) and \( x_2=\dot{q} \),

\[ \dot{x}_1=x_2,\\ \dot{x}_2=-\frac{k}{m}x_1-\frac{b}{m}x_2+\frac{1}{m}u. \]

Therefore

\[ \dot{\mathbf{x} }=A\mathbf{x}+Bu,\qquad y=C\mathbf{x}+Du, \]

\[ A=\begin{bmatrix}0&1\\-\frac{k}{m}&-\frac{b}{m}\end{bmatrix},\quad B=\begin{bmatrix}0\\\frac{1}{m}\end{bmatrix},\quad C=\begin{bmatrix}1&0\end{bmatrix},\quad D=0. \]

For constant input \( u_e \), the equilibrium satisfies \( \mathbf{0}=A\mathbf{x}_e+Bu_e \). Hence

\[ x_{2,e}=0,\qquad -\frac{k}{m}x_{1,e}+\frac{1}{m}u_e=0, \qquad y_e=x_{1,e}=\frac{u_e}{k}. \]

The open-loop characteristic polynomial is

\[ p_A(s)=\det(sI-A)=s^2+\frac{b}{m}s+\frac{k}{m}. \]

Since all coefficients are positive, the second-order Routh-Hurwitz test gives open-loop asymptotic stability. But state-feedback pole placement additionally requires controllability. The Kalman controllability matrix is

\[ \mathcal{C}=\begin{bmatrix}B&AB\end{bmatrix} =\begin{bmatrix}0&\frac{1}{m}\\\frac{1}{m}&-\frac{b}{m^2}\end{bmatrix}. \]

\[ \det(\mathcal{C})=-\frac{1}{m^2}\neq 0. \]

Thus \( \operatorname{rank}(\mathcal{C})=2 \) for \( m>0 \), so the pair \( (A,B) \) is controllable and both poles can be assigned.

Choose a standard second-order target polynomial

\[ s^2+2\zeta\omega_n s+ \omega_n^2,\qquad 0<\zeta<1,\quad \omega_n>0. \]

The corresponding poles are

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

With the approximate two-percent settling-time rule \( T_s\approx 4/(\zeta\omega_n) \), the natural frequency can be selected as

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

For \( \zeta=0.70 \) and \( T_s=2\;\mathrm{s} \), \( \omega_n=2.8571\;\mathrm{rad/s} \) and \( \lambda_{1,2}\approx -2.0000\pm j2.0408 \).

Let \( K=\begin{bmatrix}k_1&k_2\end{bmatrix} \) and \( u=-K\mathbf{x} \). Then

\[ A_{\mathrm{cl} }=A-BK =\begin{bmatrix}0&1\\-\frac{k+k_1}{m}&-\frac{b+k_2}{m}\end{bmatrix}. \]

Therefore

\[ p_{\mathrm{cl} }(s)=s^2+\frac{b+k_2}{m}s+\frac{k+k_1}{m}. \]

Matching coefficients with the target polynomial gives

\[ k_1=m\omega_n^2-k, \qquad k_2=2\zeta\omega_n m-b. \]

For \( m=1.2 \), \( b=0.8 \), \( k=3.0 \), \( \zeta=0.70 \), and \( T_s=2 \),

\[ K=\begin{bmatrix}6.7959&4.0000\end{bmatrix}. \]

This closed-form derivation is the two-state version of the same coefficient-matching idea behind companion-form and Ackermann-type pole placement.

With \( u=-K\mathbf{x}+\bar{N}r \), the closed-loop forced dynamics are

\[ \dot{\mathbf{x} }=A_{\mathrm{cl} }\mathbf{x}+B\bar{N}r, \qquad y=C\mathbf{x}. \]

For constant \( r \), the steady-state vector is

\[ \mathbf{x}_{ss}=-A_{\mathrm{cl} }^{-1}B\bar{N}r. \]

Therefore

\[ y_{ss}=-CA_{\mathrm{cl} }^{-1}B\bar{N}r. \]

Enforcing \( y_{ss}=r \) gives the standard reference scaling formula

\[ \bar{N}=-\frac{1}{C(A-BK)^{-1}B}, \qquad C(A-BK)^{-1}B\neq 0. \]

For this plant, \( \bar{N}=k+k_1=m\omega_n^2=9.7959 \).

flowchart TD
  R["reference r"] --> N["gain Nbar"]
  N --> S["sum"]
  X["state x"] --> K["gain K"]
  K --> M["negative feedback"]
  M --> S
  S --> U["force u"]
  U --> P["plant"]
  P --> Y["position y"]
  P --> X
        

The analytical verification conditions are

\[ \operatorname{rank}\begin{bmatrix}B&AB\end{bmatrix}=2, \qquad \sigma(A-BK)=\{\lambda_1,\lambda_2\}, \qquad -C(A-BK)^{-1}B\bar{N}=1. \]

Numerical simulation then checks transient response, settling time, overshoot, and maximum actuator effort.

Proposition: If all eigenvalues of \( A-BK \) satisfy \( \operatorname{Re}(\lambda_i)<0 \) and \( C(A-BK)^{-1}B\neq 0 \), then \( u=-K\mathbf{x}+\bar{N}r \) with \( \bar{N}=-1/[C(A-BK)^{-1}B] \) gives nominal asymptotic tracking of a constant reference.

Proof: Let \( A_{\mathrm{cl} }=A-BK \). Since \( A_{\mathrm{cl} } \) is Hurwitz, every solution of \( \dot{\tilde{\mathbf{x} } }=A_{\mathrm{cl} }\tilde{\mathbf{x} } \) converges to zero. For constant \( r \), define

\[ \mathbf{x}_{ss}=-A_{\mathrm{cl} }^{-1}B\bar{N}r. \]

With \( \tilde{\mathbf{x} }=\mathbf{x}-\mathbf{x}_{ss} \),

\[ \dot{\tilde{\mathbf{x} } }=A_{\mathrm{cl} }\tilde{\mathbf{x} }. \]

Thus \( \mathbf{x} \) converges to \( \mathbf{x}_{ss} \). The final output is

\[ y_{ss}=C\mathbf{x}_{ss} =-CA_{\mathrm{cl} }^{-1}B\bar{N}r =-CA_{\mathrm{cl} }^{-1}B \left(-\frac{1}{CA_{\mathrm{cl} }^{-1}B}\right)r=r. \]

The nominal constant-reference tracking proof is complete.

Using the numerical values in the implementations,

\[ A=\begin{bmatrix}0&1\\-2.5&-0.6667\end{bmatrix},\quad B=\begin{bmatrix}0\\0.8333\end{bmatrix},\quad K=\begin{bmatrix}6.7959&4.0000\end{bmatrix},\quad \bar{N}=9.7959. \]

\[ A-BK=\begin{bmatrix}0&1\\-8.1633&-4.0000\end{bmatrix}, \\ \sigma(A-BK)=\{-2.0000+j2.0408,\;-2.0000-j2.0408\}. \]

Faster pole placement increases \( \omega_n \), which increases \( k_1 \), \( k_2 \), and the likely peak force \( u(t)=-k_1q(t)-k_2\dot{q}(t)+\bar{N}r \).

Useful libraries for larger models include numpy, scipy.linalg, scipy.signal, and python-control. This file keeps the design explicit.

Chapter30_Lesson4.py

# Chapter30_Lesson4.py
# Integrated case study: mass-spring-damper model to state-feedback controller
# Libraries: NumPy for linear algebra. SciPy is optional for comparison with place_poles.

import math
import numpy as np


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return [B, AB] for a two-state SISO model."""
    return np.column_stack((B, A @ B))


def state_feedback_from_specs(m: float, b: float, k: float, zeta: float, settling_time: float):
    """Compute K = [k1, k2] by matching the desired second-order polynomial."""
    omega_n = 4.0 / (zeta * settling_time)
    k1 = m * omega_n**2 - k
    k2 = 2.0 * zeta * omega_n * m - b
    return np.array([[k1, k2]], dtype=float), omega_n


def feedforward_gain(A: np.ndarray, B: np.ndarray, C: np.ndarray, K: np.ndarray) -> float:
    """Nbar = -1/(C (A-BK)^(-1) B), valid when the scalar denominator is nonzero."""
    Acl = A - B @ K
    denom = C @ np.linalg.solve(Acl, B)
    return float(-1.0 / denom[0, 0])


def rk4_step(f, t, x, h):
    k1 = f(t, x)
    k2 = f(t + 0.5 * h, x + 0.5 * h * k1)
    k3 = f(t + 0.5 * h, x + 0.5 * h * k2)
    k4 = f(t + h, x + h * k3)
    return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate_closed_loop(A, B, C, K, nbar, r=1.0, t_final=6.0, dt=0.002):
    Acl = A - B @ K
    Bcl = B * nbar

    def f(t, x):
        return Acl @ x + Bcl.flatten() * r

    steps = int(t_final / dt)
    t = np.zeros(steps + 1)
    x = np.zeros((steps + 1, 2))
    y = np.zeros(steps + 1)
    u = np.zeros(steps + 1)
    for i in range(steps):
        t[i + 1] = t[i] + dt
        y[i] = float(C @ x[i])
        u[i] = float(-K @ x[i].reshape(-1, 1) + nbar * r)
        x[i + 1] = rk4_step(f, t[i], x[i], dt)
    y[-1] = float(C @ x[-1])
    u[-1] = float(-K @ x[-1].reshape(-1, 1) + nbar * r)
    return t, x, y, u


def main():
    m = 1.2       # kg
    b = 0.8       # N*s/m
    k = 3.0       # N/m
    zeta = 0.70
    settling_time = 2.0

    A = np.array([[0.0, 1.0], [-k / m, -b / m]])
    B = np.array([[0.0], [1.0 / m]])
    C = np.array([[1.0, 0.0]])
    D = np.array([[0.0]])

    Wc = controllability_matrix(A, B)
    print("A =\n", A)
    print("B =\n", B)
    print("Controllability matrix =\n", Wc)
    print("rank(Wc) =", np.linalg.matrix_rank(Wc))

    K, omega_n = state_feedback_from_specs(m, b, k, zeta, settling_time)
    nbar = feedforward_gain(A, B, C, K)
    Acl = A - B @ K

    print("omega_n =", omega_n)
    print("K =", K)
    print("Nbar =", nbar)
    print("closed-loop poles =", np.linalg.eigvals(Acl))

    t, x, y, u = simulate_closed_loop(A, B, C, K, nbar, r=1.0)
    print("final output y(t_final) =", y[-1])
    print("maximum absolute control input =", np.max(np.abs(u)))

    # Save a small CSV for plotting in any external tool.
    data = np.column_stack((t, x[:, 0], x[:, 1], y, u))
    np.savetxt(
        "Chapter30_Lesson4_response.csv",
        data,
        delimiter=",",
        header="t,position,velocity,output,control",
        comments="",
    )


if __name__ == "__main__":
    main()

For larger C++ state-space projects, use Eigen, Armadillo, or Blaze. This version uses standard-library arrays for transparency.

Chapter30_Lesson4.cpp

// Chapter30_Lesson4.cpp
// Integrated case study: mass-spring-damper model to state-feedback controller
// This implementation uses only the C++ standard library. For larger systems, use Eigen or Armadillo.

#include <array>
#include <cmath>
#include <fstream>
#include <iostream>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

Vec2 add(Vec2 a, Vec2 b) { return {a[0] + b[0], a[1] + b[1]}; }
Vec2 scale(double c, Vec2 a) { return {c * a[0], c * a[1]}; }
Vec2 matvec(const Mat2& A, Vec2 x) {
    return {A[0][0] * x[0] + A[0][1] * x[1], A[1][0] * x[0] + A[1][1] * x[1]};
}

double dot(Vec2 a, Vec2 b) { return a[0] * b[0] + a[1] * b[1]; }

Vec2 dynamics(const Mat2& Acl, Vec2 Bcl, double reference, Vec2 x) {
    Vec2 ax = matvec(Acl, x);
    return add(ax, scale(reference, Bcl));
}

Vec2 rk4Step(const Mat2& Acl, Vec2 Bcl, double reference, Vec2 x, double h) {
    Vec2 k1 = dynamics(Acl, Bcl, reference, x);
    Vec2 k2 = dynamics(Acl, Bcl, reference, add(x, scale(0.5 * h, k1)));
    Vec2 k3 = dynamics(Acl, Bcl, reference, add(x, scale(0.5 * h, k2)));
    Vec2 k4 = dynamics(Acl, Bcl, reference, add(x, scale(h, k3)));
    return add(x, scale(h / 6.0, add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))));
}

int main() {
    const double m = 1.2;
    const double b = 0.8;
    const double k = 3.0;
    const double zeta = 0.70;
    const double settlingTime = 2.0;
    const double omegaN = 4.0 / (zeta * settlingTime);

    const Mat2 A = { { {0.0, 1.0}, {-k / m, -b / m} } };
    const Vec2 B = {0.0, 1.0 / m};
    const Vec2 C = {1.0, 0.0};

    const double k1 = m * omegaN * omegaN - k;
    const double k2 = 2.0 * zeta * omegaN * m - b;
    const Vec2 K = {k1, k2};
    const double nbar = k + k1;

    const Mat2 Acl = { { {0.0, 1.0}, {-(k + k1) / m, -(b + k2) / m} } };
    const Vec2 Bcl = {0.0, nbar / m};

    const double trace = Acl[0][0] + Acl[1][1];
    const double determinant = Acl[0][0] * Acl[1][1] - Acl[0][1] * Acl[1][0];

    std::cout << "K = [" << K[0] << ", " << K[1] << "]\n";
    std::cout << "Nbar = " << nbar << "\n";
    std::cout << "Closed-loop characteristic: s^2 - (" << trace << ") s + " << determinant << "\n";

    const double reference = 1.0;
    const double dt = 0.002;
    const double tFinal = 6.0;
    const int steps = static_cast<int>(tFinal / dt);
    Vec2 x = {0.0, 0.0};

    std::ofstream csv("Chapter30_Lesson4_response_cpp.csv");
    csv << "t,position,velocity,output,control\n";

    double maxAbsU = 0.0;
    for (int i = 0; i <= steps; ++i) {
        double t = i * dt;
        double y = dot(C, x);
        double u = -dot(K, x) + nbar * reference;
        if (std::abs(u) > maxAbsU) maxAbsU = std::abs(u);
        csv << t << "," << x[0] << "," << x[1] << "," << y << "," << u << "\n";
        if (i < steps) x = rk4Step(Acl, Bcl, reference, x, dt);
    }

    std::cout << "Final output = " << dot(C, x) << "\n";
    std::cout << "Maximum absolute control input = " << maxAbsU << "\n";
    return 0;
}

For larger Java state-space calculations, use Apache Commons Math or EJML. This version implements the two-state case from scratch.

Chapter30_Lesson4.java

// Chapter30_Lesson4.java
// Integrated case study: mass-spring-damper model to state-feedback controller
// This implementation uses only core Java. For larger systems, use Apache Commons Math or EJML.

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;

public class Chapter30_Lesson4 {
    static double[] add(double[] a, double[] b) {
        return new double[] {a[0] + b[0], a[1] + b[1]};
    }

    static double[] scale(double c, double[] a) {
        return new double[] {c * a[0], c * a[1]};
    }

    static double[] matVec(double[][] A, double[] x) {
        return new double[] {
            A[0][0] * x[0] + A[0][1] * x[1],
            A[1][0] * x[0] + A[1][1] * x[1]
        };
    }

    static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1];
    }

    static double[] dynamics(double[][] Acl, double[] Bcl, double reference, double[] x) {
        double[] ax = matVec(Acl, x);
        return add(ax, scale(reference, Bcl));
    }

    static double[] rk4Step(double[][] Acl, double[] Bcl, double reference, double[] x, double h) {
        double[] k1 = dynamics(Acl, Bcl, reference, x);
        double[] k2 = dynamics(Acl, Bcl, reference, add(x, scale(0.5 * h, k1)));
        double[] k3 = dynamics(Acl, Bcl, reference, add(x, scale(0.5 * h, k2)));
        double[] k4 = dynamics(Acl, Bcl, reference, add(x, scale(h, k3)));
        return add(x, scale(h / 6.0, add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))));
    }

    public static void main(String[] args) throws IOException {
        double m = 1.2;
        double b = 0.8;
        double k = 3.0;
        double zeta = 0.70;
        double settlingTime = 2.0;
        double omegaN = 4.0 / (zeta * settlingTime);

        double[][] A = { {0.0, 1.0}, {-k / m, -b / m} };
        double[] B = {0.0, 1.0 / m};
        double[] C = {1.0, 0.0};

        double k1 = m * omegaN * omegaN - k;
        double k2 = 2.0 * zeta * omegaN * m - b;
        double[] K = {k1, k2};
        double nbar = k + k1;

        double[][] Acl = { {0.0, 1.0}, {-(k + k1) / m, -(b + k2) / m} };
        double[] Bcl = {0.0, nbar / m};

        double trace = Acl[0][0] + Acl[1][1];
        double determinant = Acl[0][0] * Acl[1][1] - Acl[0][1] * Acl[1][0];

        System.out.printf("K = [%.8f, %.8f]%n", K[0], K[1]);
        System.out.printf("Nbar = %.8f%n", nbar);
        System.out.printf("Closed-loop characteristic: s^2 - (%.8f) s + %.8f%n", trace, determinant);

        double reference = 1.0;
        double dt = 0.002;
        double tFinal = 6.0;
        int steps = (int) (tFinal / dt);
        double[] x = {0.0, 0.0};
        double maxAbsU = 0.0;

        try (PrintWriter csv = new PrintWriter(new FileWriter("Chapter30_Lesson4_response_java.csv"))) {
            csv.println("t,position,velocity,output,control");
            for (int i = 0; i <= steps; i++) {
                double t = i * dt;
                double y = dot(C, x);
                double u = -dot(K, x) + nbar * reference;
                maxAbsU = Math.max(maxAbsU, Math.abs(u));
                csv.printf("%.6f,%.10f,%.10f,%.10f,%.10f%n", t, x[0], x[1], y, u);
                if (i < steps) x = rk4Step(Acl, Bcl, reference, x, dt);
            }
        }

        System.out.printf("Final output = %.8f%n", dot(C, x));
        System.out.printf("Maximum absolute control input = %.8f%n", maxAbsU);
    }
}

MATLAB Control System Toolbox functions such as ss, place, eig, and lsim are useful here. The script also attempts to create a simple closed-loop Simulink State-Space model when Simulink is installed.

Chapter30_Lesson4.m

% Chapter30_Lesson4.m
% Integrated case study: mass-spring-damper model to state-feedback controller
% Libraries: base MATLAB for matrix work; Control System Toolbox for ss/place/lsim;
% Simulink section is optional and creates a closed-loop State-Space block.

clear; clc; close all;

m = 1.2;      % kg
b = 0.8;      % N*s/m
k = 3.0;      % N/m
zeta = 0.70;
Ts = 2.0;
omega_n = 4/(zeta*Ts);

A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
D = 0;

Wc = [B A*B];
fprintf('rank(Wc) = %d\n', rank(Wc));

% Manual pole-matching formula for this second-order plant.
k1 = m*omega_n^2 - k;
k2 = 2*zeta*omega_n*m - b;
K_manual = [k1 k2];

% Optional toolbox verification using place.
desiredPoles = roots([1 2*zeta*omega_n omega_n^2]);
if exist('place', 'file') == 2
    K_place = place(A, B, desiredPoles);
else
    K_place = K_manual;
end

K = K_manual;
Acl = A - B*K;
Nbar = -1/(C*(Acl\B));
Bcl = B*Nbar;

fprintf('K_manual = [%g %g]\n', K_manual(1), K_manual(2));
fprintf('K_place  = [%g %g]\n', K_place(1), K_place(2));
fprintf('Nbar = %g\n', Nbar);
disp('closed-loop poles ='); disp(eig(Acl));

% Simulate the closed-loop reference response.
t = 0:0.002:6;
r = ones(size(t));
if exist('ss', 'file') == 2 && exist('lsim', 'file') == 2
    sys_cl = ss(Acl, Bcl, C, D);
    y = lsim(sys_cl, r, t);
else
    f = @(tt, x) Acl*x + Bcl;
    [t_ode, x_ode] = ode45(f, t, [0; 0]);
    y = x_ode*C.';
    t = t_ode;
end

figure;
plot(t, y, 'LineWidth', 1.5); grid on;
xlabel('time (s)'); ylabel('position output');
title('Chapter 30 Lesson 4: closed-loop step response');

% Optional Simulink model for closed-loop state-space response.
% It uses the already-designed closed-loop matrices Acl and Bcl.
try
    if exist('new_system', 'file') == 4
        model = 'Chapter30_Lesson4_Simulink';
        if bdIsLoaded(model)
            close_system(model, 0);
        end
        new_system(model);
        add_block('simulink/Sources/Step', [model '/reference_step'], 'Position', [50 80 100 110]);
        add_block('simulink/Continuous/State-Space', [model '/closed_loop_state_space'], ...
            'A', mat2str(Acl), 'B', mat2str(Bcl), 'C', mat2str(C), 'D', mat2str(D), ...
            'Position', [170 68 310 122]);
        add_block('simulink/Sinks/Scope', [model '/scope'], 'Position', [380 75 430 115]);
        add_line(model, 'reference_step/1', 'closed_loop_state_space/1');
        add_line(model, 'closed_loop_state_space/1', 'scope/1');
        set_param(model, 'StopTime', '6');
        save_system(model);
        fprintf('Simulink model saved: %s.slx\n', model);
    end
catch ME
    fprintf('Simulink model creation skipped: %s\n', ME.message);
end

Mathematica is useful for symbolic verification, exact eigenvalue checks, and direct simulation using NDSolve.

Chapter30_Lesson4.nb

(* Chapter30_Lesson4.nb *)
(* Integrated case study: mass-spring-damper model to state-feedback controller *)

ClearAll["Global`*"];

m = 1.2;
b = 0.8;
k = 3.0;
zeta = 0.70;
settlingTime = 2.0;
omegaN = 4/(zeta*settlingTime);

A = { {0, 1}, {-k/m, -b/m} };
B = { {0}, {1/m} };
Cmat = { {1, 0} };

Wc = ArrayFlatten[{ {B, A.B} }];
MatrixRank[Wc]

k1 = m*omegaN^2 - k;
k2 = 2*zeta*omegaN*m - b;
K = { {k1, k2} };
Acl = A - B.K;
nbar = -1/First[First[Cmat.Inverse[Acl].B]];
Bcl = B*nbar;

Eigenvalues[Acl]
K
nbar

reference = 1;
solution = NDSolve[
  {
    x1'[t] == x2[t],
    x2'[t] == Acl[[2, 1]]*x1[t] + Acl[[2, 2]]*x2[t] + Bcl[[2, 1]]*reference,
    x1[0] == 0,
    x2[0] == 0
  },
  {x1, x2},
  {t, 0, 6}
];

Plot[Evaluate[x1[t] /. solution], {t, 0, 6}, PlotRange -> All,
  AxesLabel -> {"time (s)", "position output"},
  PlotLabel -> "Chapter 30 Lesson 4: closed-loop step response"]

control[t_] := Evaluate[(-K[[1, 1]]*x1[t] - K[[1, 2]]*x2[t] + nbar*reference) /. solution[[1]]];
NMaximize[{Abs[control[t]], 0 <= t <= 6}, t]

Mathematical controllability does not guarantee numerically robust gain computation. A basic conditioning diagnostic is

\[ \kappa(\mathcal{C})=\|\mathcal{C}\|\,\|\mathcal{C}^{-1}\|. \]

A similarity scaling \( \mathbf{x}=T\mathbf{z} \) changes the numerical representation:

\[ \dot{\mathbf{z} }=T^{-1}AT\mathbf{z}+T^{-1}Bu, \qquad y=CT\mathbf{z}. \]

If \( K_z \) is designed in scaled coordinates, then

\[ u=-K_z\mathbf{z}+\bar{N}r=-K_zT^{-1}\mathbf{x}+\bar{N}r, \qquad K=K_zT^{-1}. \]

The closed-loop matrices are similar, so the eigenvalues are unchanged.

Problem 1: Derive the state-space matrices for \( m\ddot{q}+b\dot{q}+kq=u \) using \( x_1=q \), \( x_2=\dot{q} \), and \( y=q \).

Solution:

\[ \dot{x}_1=x_2, \qquad \dot{x}_2=-\frac{k}{m}x_1-\frac{b}{m}x_2+\frac{1}{m}u. \]

\[ A=\begin{bmatrix}0&1\\-\frac{k}{m}&-\frac{b}{m}\end{bmatrix}, \quad B=\begin{bmatrix}0\\\frac{1}{m}\end{bmatrix}, \quad C=\begin{bmatrix}1&0\end{bmatrix},\quad D=0. \]

Problem 2: Prove controllability for \( m>0 \).

Solution: The controllability matrix is

\[ \mathcal{C}=\begin{bmatrix}0&\frac{1}{m}\\\frac{1}{m}&-\frac{b}{m^2}\end{bmatrix}. \]

Its determinant is

\[ \det(\mathcal{C})=-\frac{1}{m^2}\neq 0. \]

Hence the rank is two, so the two-state system is controllable.

Problem 3: Derive the gains for desired polynomial \( s^2+2\zeta\omega_n s+\omega_n^2 \).

Solution: Since

\[ p_{\mathrm{cl} }(s)=s^2+\frac{b+k_2}{m}s+ rac{k+k_1}{m}, \]

coefficient matching gives

\[ k_1=m\omega_n^2-k, \qquad k_2=2\zeta\omega_n m-b. \]

Problem 4: Prove the feedforward formula \( \bar{N}=-1/[C(A-BK)^{-1}B] \).

Solution: At steady state,

\[ \mathbf{x}_{ss}=-(A-BK)^{-1}B\bar{N}r. \]

Thus

\[ y_{ss}=-C(A-BK)^{-1}B\bar{N}r. \]

Choosing \( \bar{N}=-1/[C(A-BK)^{-1}B] \) gives \( y_{ss}=r \).

Problem 5: Explain why reducing \( T_s \) usually increases actuator effort when \( \zeta \) is fixed.

Solution: Since

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

reducing \( T_s \) increases \( \omega_n \). The gain \( k_1=m\omega_n^2-k \) grows quadratically, and \( k_2=2\zeta\omega_n m-b \) grows linearly. Larger gains generally produce larger peak control input and more sensitivity to noise or modeling error.

This integrated case study showed the complete path from a physical equation to a verified state-feedback controller. The central checks are controllability, closed-loop eigenvalue placement, reference DC gain, numerical conditioning, and simulated actuator effort.