After completing this lesson, students should be able to:

• derive the quarter-car vertical-dynamics equations and identify uncertain parameters;
• form a linear-in-the-parameters regressor for adaptive active suspension;
• construct a filtered tracking error and a Lyapunov-based parameter update law;
• distinguish passive, active, and semi-active suspension actuation constraints;
• evaluate ride comfort, suspension travel, road holding, and actuator effort together;
• extend parameter adaptation to a bicycle-model chassis and yaw-control problem.

Suspension control is intrinsically multi-objective. Reducing body acceleration improves ride comfort, but aggressive isolation may enlarge suspension travel or dynamic tire load. Adaptation is useful because vehicle mass changes with passengers and payload, damping is temperature dependent, tire and suspension parameters age, and road excitation is neither stationary nor known in advance. The controller must therefore learn useful combinations of uncertain parameters while respecting actuator force, stroke, damping, and tire-contact constraints.

flowchart TD
  R["Road displacement z_r"] --> T["Tire stiffness k_t"]
  T --> U["Unsprung mass m_u"]
  U --> SD["Suspension spring k_s and damper c_s"]
  SD --> B["Sprung mass m_s"]
  SENS["Body and wheel measurements"] --> ADAPT["Adaptive law and reference model"]
  ADAPT --> ACT["Active force or commanded damping"]
  ACT --> SD
  B --> SENS
  U --> SENS
        

Let \( z_s \) be the sprung-mass displacement, \( z_u \) the unsprung-mass displacement, and \( z_r \) the road displacement. The suspension deflection and relative velocity are \( q=z_s-z_u \) and \( \dot q=\dot z_s-\dot z_u \). For an ideal active actuator that applies force \( u \) between the two masses, Newton's equations are

\[ m_s\ddot z_s=-k_s(z_s-z_u)-c_s(\dot z_s-\dot z_u)+u, \]

\[ m_u\ddot z_u=k_s(z_s-z_u)+c_s(\dot z_s-\dot z_u) -k_t(z_u-z_r)-u. \]

A common state vector is \( \mathbf{x}=[z_s,\dot z_s,z_u,\dot z_u]^T \). With known parameters, the model has the linear state-space form

\[ \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}_u u+ \mathbf{B}_r z_r, \]

\[ \mathbf{A}=\begin{bmatrix} 0&1&0&0\\ -k_s/m_s&-c_s/m_s&k_s/m_s&c_s/m_s\\ 0&0&0&1\\ k_s/m_u&c_s/m_u&-(k_s+k_t)/m_u&-c_s/m_u \end{bmatrix}, \]

\[ \mathbf{B}_u=\begin{bmatrix}0&1/m_s&0&-1/m_u\end{bmatrix}^T, \qquad \mathbf{B}_r=\begin{bmatrix}0&0&0&k_t/m_u\end{bmatrix}^T. \]

The adaptive derivation below does not require all parameters to be unknown. We select \(m_s\), \(k_s\), and \(c_s\) because payload changes directly alter \(m_s\), while effective suspension stiffness and damping may vary with operating condition. The unsprung mass and tire stiffness can be included in a larger indirect estimator, but they are not needed for the body-side direct adaptive law developed here.

The principal suspension performance variables are:

\[ y_1=\ddot z_s \quad \text{(ride comfort)},\qquad y_2=z_s-z_u \quad \text{(suspension travel)}, \]

\[ y_3=k_t(z_u-z_r) \quad \text{(dynamic tire force)},\qquad y_4=u \quad \text{(actuator effort)}. \]

A finite-horizon engineering objective can be expressed as

\[ J=\int_0^T\left( w_a\ddot z_s^2+w_s(z_s-z_u)^2+w_t(z_u-z_r)^2+w_u u^2 \right)dt, \]

where the weights represent the comfort, stroke, road-holding, and energy priorities. Adaptive control does not remove this trade-off. It preserves a chosen closed-loop behavior when the plant changes. Hard constraints must still be checked:

\[ |z_s-z_u|\le z_{\max},\qquad |u|\le u_{\max},\qquad k_t(z_u-z_r)+F_{t0}>0. \]

The last inequality requires the total tire normal load to remain positive. This is a safety condition, not merely a comfort metric.

Choose a stable second-order reference model for the desired body motion:

\[ \ddot z_m+2\zeta_m\omega_m\dot z_m+\omega_m^2 z_m =\omega_m^2 z_r, \qquad \zeta_m>0,\;\omega_m>0. \]

Define the tracking error and filtered error

\[ e=z_s-z_m,\qquad s=\dot e+\lambda e,\qquad \lambda>0. \]

Introduce a virtual acceleration command

\[ v=\ddot z_m-\lambda\dot e-\kappa s, \qquad \kappa>0. \]

The body-side dynamics can be written as

\[ m_s\ddot z_s=-k_s q-c_s\dot q+u. \]

Select the direct adaptive control law

\[ u=\hat m_s v+\hat k_s q+\hat c_s\dot q =\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}, \]

\[ \hat{\boldsymbol{\theta} }= \begin{bmatrix}\hat m_s&\hat k_s&\hat c_s\end{bmatrix}^T, \qquad \boldsymbol{\phi}=\begin{bmatrix}v&q&\dot q\end{bmatrix}^T. \]

This is linear in the unknown physical parameters even though the closed-loop plant is dynamic. Let \(\tilde{\boldsymbol{\theta} }= \hat{\boldsymbol{\theta} }-\boldsymbol{\theta}\). Substitution gives the scalar filtered-error model

\[ m_s\dot s=-m_s\kappa s+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

This equation is the key adaptive-control representation: a stable nominal error dynamic plus a parameter-error term that is linear in the regressor.

Choose the Lyapunov candidate

\[ V=\frac{1}{2}m_s s^2+ \frac{1}{2}\tilde{\boldsymbol{\theta} }^T \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }, \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^T>0. \]

Since the true parameters are constant in the nominal derivation, \(\dot{\tilde{\boldsymbol{\theta} } }= \dot{\hat{\boldsymbol{\theta} } }\). Differentiating gives

\[ \dot V=m_s s\dot s+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \]

Using the filtered-error model:

\[ \dot V=-m_s\kappa s^2+s\tilde{\boldsymbol{\theta} }^T \boldsymbol{\phi}+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \]

Select the gradient update law

\[ \dot{\hat{\boldsymbol{\theta} } }=- \boldsymbol{\Gamma}\boldsymbol{\phi}s. \]

The cross term cancels exactly:

\[ \dot V=-m_s\kappa s^2\le 0. \]

Therefore \(s\) and \(\tilde{\boldsymbol{\theta} }\) are bounded. Under the standard bounded-signal conditions for the quarter-car model, \(s\in L_2\cap L_\infty\) and \(\dot s\in L_\infty\); Barbalat's lemma gives \(s\to0\). Since \(\dot e+\lambda e=s\) is a stable first-order filter, \(e\to0\) and \(\dot e\to0\).

Parameter convergence is a separate question. Tracking convergence does not imply \(\hat{\boldsymbol{\theta} }\to \boldsymbol{\theta}\). Exact parameter recovery additionally requires persistent excitation of the regressor:

\[ \exists\,T_0>0,\alpha>0:\quad \int_t^{t+T_0}\boldsymbol{\phi}(\tau) \boldsymbol{\phi}^T(\tau)d\tau\ge\alpha\mathbf{I} \quad \forall t. \]

Real suspension controllers face measurement noise, unmodeled actuator dynamics, force saturation, road-estimation errors, and slowly varying parameters. A practical update law is

\[ \dot{\hat{\boldsymbol{\theta} } }= \operatorname{Proj}_{\Omega}\left[ -\boldsymbol{\Gamma}\boldsymbol{\phi}s -\sigma\boldsymbol{\Gamma} (\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}_0) \right], \qquad \sigma>0. \]

The projection operator keeps estimates in a physically meaningful set \(\Omega\), while the leakage term pulls poorly excited estimates toward a nominal vector \(\boldsymbol{\theta}_0\). The actuator command is

\[ u=\operatorname{sat}_{[-u_{\max},u_{\max}]} (\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}). \]

Let \(\Delta_u=u-u_d\) denote the difference between the realized force and the ideal adaptive force. The filtered-error dynamics become

\[ m_s\dot s=-m_s\kappa s+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}+\Delta_u+d(t), \]

where \(d(t)\) aggregates bounded modeling error. With projection and leakage, the Lyapunov derivative admits a bound of the form

\[ \dot V\le-c_1 s^2-c_2\|\tilde{\boldsymbol{\theta} }\|^2+c_3, \qquad c_1,c_2>0, \]

so the tracking and parameter errors are uniformly ultimately bounded. Increasing leakage may reduce drift but introduces bias; increasing adaptation gains accelerates learning but amplifies noise. These are implementation trade-offs, not purely mathematical tuning choices.

A semi-active damper cannot generate an arbitrary force. In the adopted sign convention its force is

\[ u_{sa}=-c_d(t)(\dot z_s-\dot z_u), \qquad c_{\min}\le c_d(t)\le c_{\max}. \]

Consequently, the damper is dissipative because

\[ u_{sa}(\dot z_s-\dot z_u)=-c_d(t)(\dot z_s-\dot z_u)^2\le0. \]

Given a desired active force \(u_d\), a regularized clipped-damping realization is

\[ c_d^*=\operatorname{sat}_{[c_{\min},c_{\max}]} \left( -\frac{u_d\dot q}{\dot q^2+\varepsilon} \right), \qquad u_{sa}=-c_d^*\dot q, \]

where \(\varepsilon>0\) prevents division by a small relative velocity. The semi-active force generally differs from the desired active force; the mismatch must be treated as a bounded input in the robust stability analysis. For magnetorheological dampers, hysteresis states and current-to-force dynamics create an additional inner-loop control problem.

Adaptive chassis control also addresses uncertain mass, yaw inertia, tire cornering stiffness, and friction. For the linear bicycle model at longitudinal speed \(v_x>0\), define lateral velocity \(v_y\), yaw rate \(r\), steering angle \(\delta\), and direct yaw moment \(M_z\). The equations are

\[ m(\dot v_y+v_x r)=F_{yf}+F_{yr}, \qquad I_z\dot r=aF_{yf}-bF_{yr}+M_z. \]

Under the small-slip linear tire approximation,

\[ F_{yf}=C_f\left(\delta-\frac{v_y+ar}{v_x}\right), \qquad F_{yr}=C_r\left(-\frac{v_y-br}{v_x}\right). \]

Unknown cornering stiffnesses enter linearly. For example, the yaw equation can be rearranged as

\[ I_z\dot r-M_z= \begin{bmatrix} a\left(\delta-\frac{v_y+ar}{v_x}\right)& -b\left(\frac{v_y-br}{v_x}\right) \end{bmatrix} \begin{bmatrix}C_f\\C_r\end{bmatrix}. \]

This permits gradient or RLS estimation of \(C_f,C_r\), followed by adaptation of a yaw-moment or rear-steering controller. Practical systems must suspend or normalize adaptation when \(v_x\) is small, because the regressor contains \(1/v_x\). In an integrated chassis controller, vertical-load changes from suspension motion also modify available tire forces, so suspension and yaw objectives should share saturation and tire-force allocation logic.

flowchart TD
  A["Select quarter-car or half-car model"] --> B["Choose comfort, stroke, road-holding, and force metrics"]
  B --> C["Identify uncertain but linearly parameterized terms"]
  C --> D["Choose stable reference model and filtered error"]
  D --> E["Derive nominal Lyapunov adaptive law"]
  E --> F["Add projection, normalization, leakage, and force limits"]
  F --> G["Simulate bumps, random roads, payload changes, and sensor noise"]
  G --> H["Check RMS acceleration, travel, tire load, and actuator energy"]
  H --> I["Run software-in-the-loop and hardware-in-the-loop tests"]
  I --> J["Enable supervised vehicle testing with fallback controller"]
        

A fixed passive or robust fallback controller should remain available during development. Adaptation should be inhibited during sensor faults, implausible states, excessive actuator temperature, low-voltage events, loss of tire contact, or violation of the model's operating region.

The Python program uses NumPy for vector operations, a from-scratch RK4 integrator, componentwise projection, actuator saturation, and Matplotlib for visualization. The estimated parameter vector is initialized away from the true values to demonstrate bounded online adjustment.

Chapter26_Lesson4.py

"""Chapter 26, Lesson 4: adaptive active-suspension simulation.

Dependencies:
    pip install numpy matplotlib

The controller uses a linearly parameterized sprung-mass model, a filtered
tracking error, projection, sigma modification, and actuator saturation.
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Plant:
    ms: float = 350.0
    mu: float = 45.0
    ks: float = 15000.0
    cs: float = 1200.0
    kt: float = 190000.0


@dataclass(frozen=True)
class Controller:
    lambda_e: float = 8.0
    kappa: float = 12.0
    sigma: float = 0.03
    u_max: float = 4500.0
    gamma_m: float = 8.0
    gamma_k: float = 1200.0
    gamma_c: float = 250.0
    wn_ref: float = 1.2
    zeta_ref: float = 0.9


THETA_INITIAL = np.array([280.0, 11000.0, 800.0], dtype=float)
THETA_LOWER = np.array([180.0, 7000.0, 300.0], dtype=float)
THETA_UPPER = np.array([650.0, 26000.0, 3000.0], dtype=float)


def road_height(t: float) -> float:
    """A smooth 50 mm cosine bump."""
    height, start, duration = 0.05, 1.0, 0.35
    if start <= t <= start + duration:
        tau = (t - start) / duration
        return 0.5 * height * (1.0 - math.cos(2.0 * math.pi * tau))
    return 0.0


def projected_rate(theta: np.ndarray, raw_rate: np.ndarray) -> np.ndarray:
    """Componentwise projection onto known physically plausible bounds."""
    rate = raw_rate.copy()
    for i in range(theta.size):
        moving_below = theta[i] <= THETA_LOWER[i] and rate[i] < 0.0
        moving_above = theta[i] >= THETA_UPPER[i] and rate[i] > 0.0
        if moving_below or moving_above:
            rate[i] = 0.0
    return rate


def rhs(t: float, y: np.ndarray, plant: Plant, ctrl: Controller) -> np.ndarray:
    """Continuous dynamics for plant, reference model, and parameter estimates."""
    zs, zsd, zu, zud, zm, zmd, hat_m, hat_k, hat_c = y
    theta_hat = np.array([hat_m, hat_k, hat_c], dtype=float)

    zr = road_height(t)
    zmdd = (
        -2.0 * ctrl.zeta_ref * ctrl.wn_ref * zmd
        - ctrl.wn_ref**2 * zm
        + ctrl.wn_ref**2 * zr
    )

    error = zs - zm
    error_rate = zsd - zmd
    sliding_error = error_rate + ctrl.lambda_e * error
    suspension_deflection = zs - zu
    relative_velocity = zsd - zud

    virtual_acceleration = (
        zmdd
        - ctrl.lambda_e * error_rate
        - ctrl.kappa * sliding_error
    )
    regressor = np.array(
        [virtual_acceleration, suspension_deflection, relative_velocity],
        dtype=float,
    )

    commanded_force = float(theta_hat @ regressor)
    actuator_force = float(np.clip(commanded_force, -ctrl.u_max, ctrl.u_max))

    gamma = np.array([ctrl.gamma_m, ctrl.gamma_k, ctrl.gamma_c], dtype=float)
    raw_adaptation = (
        -gamma * regressor * sliding_error
        - ctrl.sigma * (theta_hat - THETA_INITIAL)
    )
    theta_rate = projected_rate(theta_hat, raw_adaptation)

    sprung_acceleration = (
        -plant.ks * suspension_deflection
        - plant.cs * relative_velocity
        + actuator_force
    ) / plant.ms
    unsprung_acceleration = (
        plant.ks * suspension_deflection
        + plant.cs * relative_velocity
        - plant.kt * (zu - zr)
        - actuator_force
    ) / plant.mu

    return np.array(
        [
            zsd,
            sprung_acceleration,
            zud,
            unsprung_acceleration,
            zmd,
            zmdd,
            theta_rate[0],
            theta_rate[1],
            theta_rate[2],
        ],
        dtype=float,
    )


def rk4_step(
    t: float, y: np.ndarray, step: float, plant: Plant, ctrl: Controller
) -> np.ndarray:
    k1 = rhs(t, y, plant, ctrl)
    k2 = rhs(t + 0.5 * step, y + 0.5 * step * k1, plant, ctrl)
    k3 = rhs(t + 0.5 * step, y + 0.5 * step * k2, plant, ctrl)
    k4 = rhs(t + step, y + step * k3, plant, ctrl)
    next_y = y + (step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    next_y[6:9] = np.clip(next_y[6:9], THETA_LOWER, THETA_UPPER)
    return next_y


def simulate(duration: float = 8.0, step: float = 0.001):
    plant, ctrl = Plant(), Controller()
    time = np.arange(0.0, duration + step, step)
    state = np.zeros((time.size, 9), dtype=float)
    state[0, 6:9] = THETA_INITIAL

    for i in range(time.size - 1):
        state[i + 1] = rk4_step(time[i], state[i], step, plant, ctrl)

    acceleration = np.zeros(time.size)
    force = np.zeros(time.size)
    road = np.zeros(time.size)
    for i, (t, y) in enumerate(zip(time, state)):
        derivative = rhs(t, y, plant, ctrl)
        acceleration[i] = derivative[1]
        road[i] = road_height(t)

        zs, zsd, zu, zud, zm, zmd, hat_m, hat_k, hat_c = y
        zmdd = (
            -2.0 * ctrl.zeta_ref * ctrl.wn_ref * zmd
            - ctrl.wn_ref**2 * zm
            + ctrl.wn_ref**2 * road[i]
        )
        error = zs - zm
        error_rate = zsd - zmd
        sliding_error = error_rate + ctrl.lambda_e * error
        virtual_acceleration = (
            zmdd - ctrl.lambda_e * error_rate - ctrl.kappa * sliding_error
        )
        regressor = np.array(
            [virtual_acceleration, zs - zu, zsd - zud], dtype=float
        )
        force[i] = np.clip(
            np.array([hat_m, hat_k, hat_c]) @ regressor,
            -ctrl.u_max,
            ctrl.u_max,
        )

    return time, state, road, acceleration, force


def main() -> None:
    time, state, road, acceleration, force = simulate()
    suspension_travel = state[:, 0] - state[:, 2]
    tire_deflection = state[:, 2] - road

    print(f"RMS body acceleration: {np.sqrt(np.mean(acceleration**2)):.6f} m/s^2")
    print(f"Maximum suspension travel: {np.max(np.abs(suspension_travel)):.6f} m")
    print(f"Maximum tire deflection: {np.max(np.abs(tire_deflection)):.6f} m")
    print(f"Maximum actuator force: {np.max(np.abs(force)):.3f} N")
    print("Final estimates [m_s, k_s, c_s]:", state[-1, 6:9])

    fig, axes = plt.subplots(4, 1, figsize=(10, 11), sharex=True)
    axes[0].plot(time, road, label="road")
    axes[0].plot(time, state[:, 0], label="sprung mass")
    axes[0].plot(time, state[:, 4], "--", label="reference model")
    axes[0].set_ylabel("displacement [m]")
    axes[0].legend()
    axes[0].grid(True)

    axes[1].plot(time, acceleration)
    axes[1].set_ylabel("body acceleration [m/s^2]")
    axes[1].grid(True)

    axes[2].plot(time, force)
    axes[2].set_ylabel("actuator force [N]")
    axes[2].grid(True)

    axes[3].plot(time, state[:, 6], label="hat m_s")
    axes[3].plot(time, state[:, 7] / 50.0, label="hat k_s / 50")
    axes[3].plot(time, state[:, 8], label="hat c_s")
    axes[3].set_xlabel("time [s]")
    axes[3].set_ylabel("scaled estimates")
    axes[3].legend()
    axes[3].grid(True)

    fig.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

The C++17 implementation is self-contained and uses fixed-size std::array objects. For a full-car MIMO controller, Eigen is a suitable replacement for the small hand-written vector operations.

Chapter26_Lesson4.cpp

// Chapter 26, Lesson 4: self-contained adaptive active-suspension simulation.
// Build: g++ -std=c++17 -O2 Chapter26_Lesson4.cpp -o Chapter26_Lesson4

#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>

namespace {
constexpr double PI = 3.14159265358979323846;
using State = std::array<double, 9>;

struct Plant {
    double ms = 350.0;
    double mu = 45.0;
    double ks = 15000.0;
    double cs = 1200.0;
    double kt = 190000.0;
};

struct Controller {
    double lambda_e = 8.0;
    double kappa = 12.0;
    double sigma = 0.03;
    double u_max = 4500.0;
    std::array<double, 3> gamma{8.0, 1200.0, 250.0};
    double wn_ref = 1.2;
    double zeta_ref = 0.9;
};

constexpr std::array<double, 3> THETA_INITIAL{280.0, 11000.0, 800.0};
constexpr std::array<double, 3> THETA_LOWER{180.0, 7000.0, 300.0};
constexpr std::array<double, 3> THETA_UPPER{650.0, 26000.0, 3000.0};

double roadHeight(double t) {
    constexpr double height = 0.05;
    constexpr double start = 1.0;
    constexpr double duration = 0.35;
    if (t >= start && t <= start + duration) {
        const double tau = (t - start) / duration;
        return 0.5 * height * (1.0 - std::cos(2.0 * PI * tau));
    }
    return 0.0;
}

State addScaled(const State& a, const State& b, double scale) {
    State out{};
    for (std::size_t i = 0; i < out.size(); ++i) {
        out[i] = a[i] + scale * b[i];
    }
    return out;
}

State derivative(double t, const State& y, const Plant& p, const Controller& c) {
    const double zs = y[0], zsd = y[1], zu = y[2], zud = y[3];
    const double zm = y[4], zmd = y[5];
    const std::array<double, 3> theta{y[6], y[7], y[8]};

    const double zr = roadHeight(t);
    const double zmdd = -2.0 * c.zeta_ref * c.wn_ref * zmd
                      - c.wn_ref * c.wn_ref * zm
                      + c.wn_ref * c.wn_ref * zr;
    const double e = zs - zm;
    const double ed = zsd - zmd;
    const double s = ed + c.lambda_e * e;
    const double q = zs - zu;
    const double qd = zsd - zud;
    const double v = zmdd - c.lambda_e * ed - c.kappa * s;
    const std::array<double, 3> phi{v, q, qd};

    double uCommand = 0.0;
    for (int i = 0; i < 3; ++i) {
        uCommand += theta[i] * phi[i];
    }
    const double u = std::clamp(uCommand, -c.u_max, c.u_max);

    std::array<double, 3> thetaRate{};
    for (int i = 0; i < 3; ++i) {
        thetaRate[i] = -c.gamma[i] * phi[i] * s
                     - c.sigma * (theta[i] - THETA_INITIAL[i]);
        const bool below = theta[i] <= THETA_LOWER[i] && thetaRate[i] < 0.0;
        const bool above = theta[i] >= THETA_UPPER[i] && thetaRate[i] > 0.0;
        if (below || above) {
            thetaRate[i] = 0.0;
        }
    }

    const double zsdd = (-p.ks * q - p.cs * qd + u) / p.ms;
    const double zudd = (p.ks * q + p.cs * qd - p.kt * (zu - zr) - u) / p.mu;

    return State{zsd, zsdd, zud, zudd, zmd, zmdd,
                 thetaRate[0], thetaRate[1], thetaRate[2]};
}

State rk4(double t, const State& y, double h, const Plant& p, const Controller& c) {
    const State k1 = derivative(t, y, p, c);
    const State k2 = derivative(t + 0.5 * h, addScaled(y, k1, 0.5 * h), p, c);
    const State k3 = derivative(t + 0.5 * h, addScaled(y, k2, 0.5 * h), p, c);
    const State k4 = derivative(t + h, addScaled(y, k3, h), p, c);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = y[i] + (h / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    for (int i = 0; i < 3; ++i) {
        next[6 + i] = std::clamp(next[6 + i], THETA_LOWER[i], THETA_UPPER[i]);
    }
    return next;
}
}  // namespace

int main() {
    const Plant plant;
    const Controller controller;
    constexpr double step = 0.001;
    constexpr double duration = 8.0;

    State state{};
    state[6] = THETA_INITIAL[0];
    state[7] = THETA_INITIAL[1];
    state[8] = THETA_INITIAL[2];

    double accelerationSquareIntegral = 0.0;
    double maxSuspensionTravel = 0.0;
    double maxTireDeflection = 0.0;

    const int sampleCount = static_cast<int>(duration / step) + 1;
    for (int i = 0; i < sampleCount - 1; ++i) {
        const double t = i * step;
        const State dy = derivative(t, state, plant, controller);
        accelerationSquareIntegral += dy[1] * dy[1] * step;
        maxSuspensionTravel = std::max(maxSuspensionTravel, std::abs(state[0] - state[2]));
        maxTireDeflection = std::max(maxTireDeflection, std::abs(state[2] - roadHeight(t)));
        state = rk4(t, state, step, plant, controller);
    }

    const double rmsAcceleration = std::sqrt(accelerationSquareIntegral / duration);
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "RMS body acceleration: " << rmsAcceleration << " m/s^2\n";
    std::cout << "Maximum suspension travel: " << maxSuspensionTravel << " m\n";
    std::cout << "Maximum tire deflection: " << maxTireDeflection << " m\n";
    std::cout << "Final estimates [m_s, k_s, c_s]: "
              << state[6] << ", " << state[7] << ", " << state[8] << "\n";
    return 0;
}

The Java version mirrors the C++ equations and requires only the JDK. Apache Commons Math can be introduced later for matrix-based half-car or full-car models.

Chapter26_Lesson4.java

// Chapter 26, Lesson 4: self-contained adaptive active-suspension simulation.
// Build and run: javac Chapter26_Lesson4.java && java Chapter26_Lesson4

import java.util.Arrays;

public final class Chapter26_Lesson4 {
    private static final double PI = Math.PI;
    private static final double[] THETA_INITIAL = {280.0, 11000.0, 800.0};
    private static final double[] THETA_LOWER = {180.0, 7000.0, 300.0};
    private static final double[] THETA_UPPER = {650.0, 26000.0, 3000.0};

    private static final class Plant {
        double ms = 350.0;
        double mu = 45.0;
        double ks = 15000.0;
        double cs = 1200.0;
        double kt = 190000.0;
    }

    private static final class Controller {
        double lambdaE = 8.0;
        double kappa = 12.0;
        double sigma = 0.03;
        double uMax = 4500.0;
        double[] gamma = {8.0, 1200.0, 250.0};
        double wnRef = 1.2;
        double zetaRef = 0.9;
    }

    private static double clamp(double value, double lower, double upper) {
        return Math.max(lower, Math.min(upper, value));
    }

    private static double roadHeight(double t) {
        final double height = 0.05;
        final double start = 1.0;
        final double duration = 0.35;
        if (t >= start && t <= start + duration) {
            double tau = (t - start) / duration;
            return 0.5 * height * (1.0 - Math.cos(2.0 * PI * tau));
        }
        return 0.0;
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        double[] out = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            out[i] = a[i] + scale * b[i];
        }
        return out;
    }

    private static double[] derivative(
            double t, double[] y, Plant p, Controller c) {
        double zs = y[0], zsd = y[1], zu = y[2], zud = y[3];
        double zm = y[4], zmd = y[5];
        double[] theta = {y[6], y[7], y[8]};

        double zr = roadHeight(t);
        double zmdd = -2.0 * c.zetaRef * c.wnRef * zmd
                - c.wnRef * c.wnRef * zm
                + c.wnRef * c.wnRef * zr;
        double error = zs - zm;
        double errorRate = zsd - zmd;
        double slidingError = errorRate + c.lambdaE * error;
        double deflection = zs - zu;
        double relativeVelocity = zsd - zud;
        double virtualAcceleration = zmdd
                - c.lambdaE * errorRate
                - c.kappa * slidingError;
        double[] phi = {virtualAcceleration, deflection, relativeVelocity};

        double command = theta[0] * phi[0] + theta[1] * phi[1] + theta[2] * phi[2];
        double force = clamp(command, -c.uMax, c.uMax);

        double[] thetaRate = new double[3];
        for (int i = 0; i < 3; i++) {
            thetaRate[i] = -c.gamma[i] * phi[i] * slidingError
                    - c.sigma * (theta[i] - THETA_INITIAL[i]);
            boolean below = theta[i] <= THETA_LOWER[i] && thetaRate[i] < 0.0;
            boolean above = theta[i] >= THETA_UPPER[i] && thetaRate[i] > 0.0;
            if (below || above) {
                thetaRate[i] = 0.0;
            }
        }

        double sprungAcceleration = (-p.ks * deflection - p.cs * relativeVelocity + force) / p.ms;
        double unsprungAcceleration = (p.ks * deflection + p.cs * relativeVelocity
                - p.kt * (zu - zr) - force) / p.mu;

        return new double[] {
            zsd, sprungAcceleration, zud, unsprungAcceleration,
            zmd, zmdd, thetaRate[0], thetaRate[1], thetaRate[2]
        };
    }

    private static double[] rk4(
            double t, double[] y, double h, Plant p, Controller c) {
        double[] k1 = derivative(t, y, p, c);
        double[] k2 = derivative(t + 0.5 * h, addScaled(y, k1, 0.5 * h), p, c);
        double[] k3 = derivative(t + 0.5 * h, addScaled(y, k2, 0.5 * h), p, c);
        double[] k4 = derivative(t + h, addScaled(y, k3, h), p, c);

        double[] next = new double[y.length];
        for (int i = 0; i < y.length; i++) {
            next[i] = y[i] + (h / 6.0)
                    * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        for (int i = 0; i < 3; i++) {
            next[6 + i] = clamp(next[6 + i], THETA_LOWER[i], THETA_UPPER[i]);
        }
        return next;
    }

    public static void main(String[] args) {
        Plant plant = new Plant();
        Controller controller = new Controller();
        double step = 0.001;
        double duration = 8.0;
        int samples = (int) (duration / step) + 1;

        double[] state = new double[9];
        System.arraycopy(THETA_INITIAL, 0, state, 6, 3);

        double accelerationSquareIntegral = 0.0;
        double maxSuspensionTravel = 0.0;
        double maxTireDeflection = 0.0;

        for (int i = 0; i < samples - 1; i++) {
            double t = i * step;
            double[] dy = derivative(t, state, plant, controller);
            accelerationSquareIntegral += dy[1] * dy[1] * step;
            maxSuspensionTravel = Math.max(
                    maxSuspensionTravel, Math.abs(state[0] - state[2]));
            maxTireDeflection = Math.max(
                    maxTireDeflection, Math.abs(state[2] - roadHeight(t)));
            state = rk4(t, state, step, plant, controller);
        }

        double rmsAcceleration = Math.sqrt(accelerationSquareIntegral / duration);
        System.out.printf("RMS body acceleration: %.6f m/s^2%n", rmsAcceleration);
        System.out.printf("Maximum suspension travel: %.6f m%n", maxSuspensionTravel);
        System.out.printf("Maximum tire deflection: %.6f m%n", maxTireDeflection);
        System.out.println("Final estimates [m_s, k_s, c_s]: "
                + Arrays.toString(Arrays.copyOfRange(state, 6, 9)));
    }
}

The MATLAB script uses a local-function structure and an explicit RK4 solver. A direct Simulink realization uses two second-order mechanical channels or four Integrator blocks for the plant, two Integrator blocks for the reference model, a MATLAB Function block for the adaptive law, Saturation blocks for actuator and parameter limits, and To Workspace or Simulation Data Inspector logging for the performance outputs.

Chapter26_Lesson4.m

% Chapter 26, Lesson 4: adaptive active-suspension simulation.
% Run this script directly in MATLAB. It uses only core MATLAB functions.
% For Simulink, map the same equations to Integrator, Sum, Gain,
% Saturation, and MATLAB Function blocks as described in the lesson.

clear; clc; close all;

p.ms = 350; p.mu = 45; p.ks = 15000; p.cs = 1200; p.kt = 190000;
c.lambda = 8; c.kappa = 12; c.sigma = 0.03; c.umax = 4500;
c.gamma = [8; 1200; 250]; c.wn = 1.2; c.zeta = 0.9;

theta0 = [280; 11000; 800];
thetaLower = [180; 7000; 300];
thetaUpper = [650; 26000; 3000];

dt = 1e-3; tf = 8; t = (0:dt:tf)';
y = zeros(numel(t), 9); y(1, 7:9) = theta0';

for k = 1:numel(t)-1
    k1 = dynamics(t(k), y(k,:)', p, c, theta0, thetaLower, thetaUpper);
    k2 = dynamics(t(k)+dt/2, y(k,:)'+dt*k1/2, p, c, theta0, thetaLower, thetaUpper);
    k3 = dynamics(t(k)+dt/2, y(k,:)'+dt*k2/2, p, c, theta0, thetaLower, thetaUpper);
    k4 = dynamics(t(k)+dt, y(k,:)'+dt*k3, p, c, theta0, thetaLower, thetaUpper);
    y(k+1,:) = (y(k,:)' + dt*(k1 + 2*k2 + 2*k3 + k4)/6)';
    y(k+1,7:9) = min(max(y(k+1,7:9), thetaLower'), thetaUpper');
end

road = arrayfun(@roadHeight, t);
acceleration = zeros(size(t));
force = zeros(size(t));
for k = 1:numel(t)
    dy = dynamics(t(k), y(k,:)', p, c, theta0, thetaLower, thetaUpper);
    acceleration(k) = dy(2);
    force(k) = controllerForce(t(k), y(k,:)', c);
end

suspensionTravel = y(:,1) - y(:,3);
tireDeflection = y(:,3) - road;
fprintf('RMS body acceleration: %.6f m/s^2\n', rms(acceleration));
fprintf('Maximum suspension travel: %.6f m\n', max(abs(suspensionTravel)));
fprintf('Maximum tire deflection: %.6f m\n', max(abs(tireDeflection)));
fprintf('Maximum actuator force: %.3f N\n', max(abs(force)));
fprintf('Final estimates [m_s, k_s, c_s]: %.3f, %.3f, %.3f\n', y(end,7:9));

figure('Name','Adaptive Suspension');
tiledlayout(4,1);
nexttile; plot(t, road, t, y(:,1), t, y(:,5),'--'); grid on;
ylabel('Displacement [m]'); legend('Road','Sprung mass','Reference');
nexttile; plot(t, acceleration); grid on; ylabel('Body acc. [m/s^2]');
nexttile; plot(t, force); grid on; ylabel('Force [N]');
nexttile; plot(t, y(:,7), t, y(:,8)/50, t, y(:,9)); grid on;
ylabel('Scaled estimates'); xlabel('Time [s]');
legend('hat m_s','hat k_s / 50','hat c_s');

function dy = dynamics(t, y, p, c, theta0, thetaLower, thetaUpper)
    zs=y(1); zsd=y(2); zu=y(3); zud=y(4); zm=y(5); zmd=y(6);
    theta=y(7:9); zr=roadHeight(t);
    zmdd=-2*c.zeta*c.wn*zmd-c.wn^2*zm+c.wn^2*zr;
    e=zs-zm; ed=zsd-zmd; s=ed+c.lambda*e;
    q=zs-zu; qd=zsd-zud;
    v=zmdd-c.lambda*ed-c.kappa*s;
    phi=[v; q; qd];
    u=min(max(theta'*phi,-c.umax),c.umax);
    thetaDot=-c.gamma.*phi*s-c.sigma*(theta-theta0);
    for i=1:3
        if (theta(i)<=thetaLower(i) && thetaDot(i)<0) || ...
           (theta(i)>=thetaUpper(i) && thetaDot(i)>0)
            thetaDot(i)=0;
        end
    end
    zsdd=(-p.ks*q-p.cs*qd+u)/p.ms;
    zudd=(p.ks*q+p.cs*qd-p.kt*(zu-zr)-u)/p.mu;
    dy=[zsd; zsdd; zud; zudd; zmd; zmdd; thetaDot];
end

function u = controllerForce(t, y, c)
    zs=y(1); zsd=y(2); zu=y(3); zud=y(4); zm=y(5); zmd=y(6);
    theta=y(7:9); zr=roadHeight(t);
    zmdd=-2*c.zeta*c.wn*zmd-c.wn^2*zm+c.wn^2*zr;
    e=zs-zm; ed=zsd-zmd; s=ed+c.lambda*e;
    phi=[zmdd-c.lambda*ed-c.kappa*s; zs-zu; zsd-zud];
    u=min(max(theta'*phi,-c.umax),c.umax);
end

function zr = roadHeight(t)
    h=0.05; t0=1.0; duration=0.35;
    if t>=t0 && t<=t0+duration
        tau=(t-t0)/duration;
        zr=0.5*h*(1-cos(2*pi*tau));
    else
        zr=0;
    end
end

The notebook expression defines the model, projection logic, RK4 map, simulation, RMS acceleration calculation, and displacement plot. The implementation is intentionally numerical so that it follows exactly the same saturation and projection logic as the other languages.

Chapter26_Lesson4.nb

Notebook[{
  Cell["Chapter 26, Lesson 4: Adaptive Suspension and Chassis Systems", "Title"],
  Cell["This notebook simulates the same projected, sigma-modified adaptive active-suspension controller used in the Python, C++, Java, and MATLAB files.", "Text"],
  Cell[BoxData[
"ClearAll[\"Global`*\"];
ms = 350.; mu = 45.; ks = 15000.; cs = 1200.; kt = 190000.;
lambda = 8.; kappa = 12.; sigma = 0.03; umax = 4500.;
gamma = {8., 1200., 250.}; wn = 1.2; zeta = 0.9;
theta0 = {280., 11000., 800.};
thetaLower = {180., 7000., 300.};
thetaUpper = {650., 26000., 3000.};
road[t_] := Piecewise[{ {0.025 (1 - Cos[2 Pi (t - 1.)/0.35]),
  1. <= t <= 1.35} }, 0.];
clip[x_, lo_, hi_] := Min[Max[x, lo], hi];
projectedRate[theta_, raw_] := MapThread[
  If[(#1 <= #3 && #2 < 0) || (#1 >= #4 && #2 > 0), 0., #2] &,
  {theta, raw, thetaLower, thetaUpper}];
rhs[t_, y_] := Module[
  {zs, zsd, zu, zud, zm, zmd, theta, zr, zmdd, e, ed, s, q, qd,
   v, phi, u, raw, td, zsdd, zudd},
  {zs, zsd, zu, zud, zm, zmd} = y[[1 ;; 6]];
  theta = y[[7 ;; 9]]; zr = road[t];
  zmdd = -2 zeta wn zmd - wn^2 zm + wn^2 zr;
  e = zs - zm; ed = zsd - zmd; s = ed + lambda e;
  q = zs - zu; qd = zsd - zud;
  v = zmdd - lambda ed - kappa s; phi = {v, q, qd};
  u = clip[theta.phi, -umax, umax];
  raw = -gamma phi s - sigma (theta - theta0);
  td = projectedRate[theta, raw];
  zsdd = (-ks q - cs qd + u)/ms;
  zudd = (ks q + cs qd - kt (zu - zr) - u)/mu;
  Join[{zsd, zsdd, zud, zudd, zmd, zmdd}, td]
];
rk4[t_, y_, h_] := Module[{k1, k2, k3, k4, yn},
  k1 = rhs[t, y]; k2 = rhs[t + h/2, y + h k1/2];
  k3 = rhs[t + h/2, y + h k2/2];
  k4 = rhs[t + h, y + h k3];
  yn = y + h (k1 + 2 k2 + 2 k3 + k4)/6;
  Join[yn[[1 ;; 6]],
   MapThread[clip, {yn[[7 ;; 9]], thetaLower, thetaUpper}]]
];

dt = 0.001; times = Range[0., 8., dt];
y0 = Join[ConstantArray[0., 6], theta0];
states = Rest@FoldList[rk4[#2 - dt, #1, dt] &, y0, Rest[times]];
states = Prepend[states, y0];
accelerations = MapThread[rhs[#1, #2][[2]] &, {times, states}];
Print[\"RMS body acceleration: \", Sqrt[Mean[accelerations^2]],
  \" m/s^2\"];
Print[\"Final estimates: \", Last[states][[7 ;; 9]]];

ListLinePlot[
 {Transpose[{times, road /@ times}],
  Transpose[{times, states[[All, 1]]}],
  Transpose[{times, states[[All, 5]]}]},
 PlotLegends -> {\"Road\", \"Sprung mass\", \"Reference\"},
 Frame -> True,
 FrameLabel -> {\"Time [s]\", \"Displacement [m]\"}]
"], "Input"]
}, WindowTitle -> "Chapter26_Lesson4"]        

The example applies a smooth 50 mm road bump. The controller estimates the sprung mass, suspension stiffness, and damping while forcing the body motion toward the stable reference model. Because a single transient is not persistently exciting in all parameter directions, the estimates are not expected to converge exactly to the physical values. The correct interpretation is:

• tracking and all internal signals should remain bounded;
• the filtered error should decay after the bump;
• the projection bounds should never be violated;
• the actuator command should remain inside its force limit;
• parameter estimates may settle at values that reproduce the required closed-loop force without identifying the plant uniquely.

To test parameter convergence, replace the isolated bump with a bounded broadband road profile or a sequence of bumps with sufficiently rich frequency content, while retaining safe stroke and tire-load limits.

Problem 1 — State-Space Derivation: Starting from the two quarter-car equations, derive the matrices \(\mathbf{A}\), \(\mathbf{B}_u\), and \(\mathbf{B}_r\) for \(\mathbf{x}=[z_s,\dot z_s,z_u,\dot z_u]^T\).

Solution: Set

\[ x_1=z_s,\quad x_2=\dot z_s,\quad x_3=z_u,\quad x_4=\dot z_u. \]

Then \(\dot x_1=x_2\) and \(\dot x_3=x_4\). Dividing the force equations by the respective masses gives

\[ \dot x_2=-\frac{k_s}{m_s}x_1-\frac{c_s}{m_s}x_2+ \frac{k_s}{m_s}x_3+\frac{c_s}{m_s}x_4+\frac{1}{m_s}u, \]

\[ \dot x_4=\frac{k_s}{m_u}x_1+\frac{c_s}{m_u}x_2- \frac{k_s+k_t}{m_u}x_3-\frac{c_s}{m_u}x_4- \frac{1}{m_u}u+\frac{k_t}{m_u}z_r. \]

Collecting coefficients yields the matrices in Section 2.

Problem 2 — Filtered-Error Dynamics: Show that the control law \(u=\hat m_s v+\hat k_s q+\hat c_s\dot q\) with \(v=\ddot z_m-\lambda\dot e-\kappa s\) produces \(m_s\dot s=-m_s\kappa s+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}\).

Solution: Since \(s=\dot e+\lambda e\),

\[ \dot s=\ddot z_s-\ddot z_m+\lambda\dot e. \]

From the sprung-mass equation and the control law,

\[ m_s\ddot z_s=-k_s q-c_s\dot q+ \hat m_s v+\hat k_s q+\hat c_s\dot q =m_s v+\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Substitute \(v\) and the expression for \(\dot s\):

\[ m_s\dot s=m_s(v-\ddot z_m+\lambda\dot e)+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi} =-m_s\kappa s+ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Problem 3 — Lyapunov Cancellation: Derive the gradient adaptation law that cancels the parameter cross term in \(\dot V\).

Solution: With

\[ V=\frac{1}{2}m_s s^2+\frac{1}{2} \tilde{\boldsymbol{\theta} }^T\boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }, \]

differentiation and substitution give

\[ \dot V=-m_s\kappa s^2+ \tilde{\boldsymbol{\theta} }^T \left(\boldsymbol{\phi}s+ \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta} } }\right). \]

Set the term in parentheses to zero:

\[ \dot{\hat{\boldsymbol{\theta} } }=- \boldsymbol{\Gamma}\boldsymbol{\phi}s. \]

This gives \(\dot V=-m_s\kappa s^2\le0\).

Problem 4 — Semi-Active Force Projection: For a desired active force \(u_d\) and relative velocity \(\dot q\), derive the unconstrained damping coefficient that minimizes \((u_d+c_d\dot q)^2\). Then impose \(c_{\min}\le c_d\le c_{\max}\).

Solution: Define

\[ L(c_d)=(u_d+c_d\dot q)^2. \]

For \(\dot q\ne0\), the stationary point satisfies

\[ \frac{dL}{dc_d}=2(u_d+c_d\dot q)\dot q=0, \qquad c_d=-\frac{u_d}{\dot q}. \]

A numerically regularized and physically admissible command is

\[ c_d^*=\operatorname{sat}_{[c_{\min},c_{\max}]} \left(-\frac{u_d\dot q}{\dot q^2+\varepsilon}\right). \]

The realized force is \(u_{sa}=-c_d^*\dot q\) and always satisfies dissipativity.

Problem 5 — Ultimate-Bound Estimate: Suppose the robust closed-loop analysis gives \(\dot V\le-\alpha\|\mathbf{z}\|^2+\beta\), where \(\mathbf{z}=[s,\tilde{\boldsymbol{\theta} }^T]^T\), \(\alpha>0\), and \(\beta>0\). Determine a sufficient ultimate-bound radius.

Solution: Whenever

\[ \|\mathbf{z}\|^2>\frac{\beta}{\alpha}, \]

the derivative is negative. Therefore trajectories enter and remain in the set

\[ \|\mathbf{z}\|\le\sqrt{\frac{\beta}{\alpha} }. \]

Tighter bounds require retaining the exact quadratic coefficients rather than replacing them by the common lower bound \(\alpha\).

Problem 6 — Chassis Regressor: Derive a regression for estimating \(C_f\) and \(C_r\) from the bicycle-model yaw equation when \(I_z\) is known.

Solution: Substitute the linear tire forces into

\[ I_z\dot r-M_z=aF_{yf}-bF_{yr}. \]

This yields

\[ y_c=\boldsymbol{\phi}_c^T\boldsymbol{\theta}_c, \]

\[ y_c=I_z\dot r-M_z, \quad \boldsymbol{\theta}_c=\begin{bmatrix}C_f&C_r\end{bmatrix}^T, \]

\[ \boldsymbol{\phi}_c= \begin{bmatrix} a\left(\delta-\frac{v_y+ar}{v_x}\right)\\ -b\left(\frac{v_y-br}{v_x}\right) \end{bmatrix}. \]

A normalized gradient estimator can use \(\dot{\hat{\boldsymbol{\theta} } }_c= \boldsymbol{\Gamma}_c\boldsymbol{\phi}_c (y_c-\boldsymbol{\phi}_c^T\hat{\boldsymbol{\theta} }_c)/ (1+\boldsymbol{\phi}_c^T\boldsymbol{\phi}_c)\), with projection onto positive stiffness bounds.

An adaptive suspension controller can be derived directly from the sprung-mass equation by selecting a stable reference model, defining the filtered error \(s=\dot e+\lambda e\), and expressing the actuator force as a linear function of estimated mass, stiffness, and damping. A gradient update law cancels the parameter-error cross term in the Lyapunov derivative and guarantees nominal tracking convergence. Projection, leakage, normalization, and saturation are essential for practical boundedness. Semi-active hardware requires an additional dissipative force projection, while lateral chassis adaptation uses the same regressor concept to estimate tire and vehicle parameters. The resulting design must always be evaluated against comfort, suspension travel, tire contact, force, energy, and fault-management constraints.