In earlier lessons of this chapter we described nonlinear dynamics as \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u}, t) \) and studied equilibria and linearization. Here we focus on a very common and structurally important subclass: a linear plant interconnected with a nonlinear element (often in the actuator or transmission path).

A standard abstraction is a linear state-space plant driven by a scalar “effective input” \( u \) produced by a nonlinearity \( \varphi(\cdot) \):

\[ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\,u, \qquad u = \varphi(v), \qquad v = \mathbf{C}\mathbf{x} + \mathbf{D}\,r(t). \]

The modeling task is to specify \( \varphi(\cdot) \) with correct structure: memoryless (saturation, dead-zone) versus hysteretic (backlash), and then to reason about the resulting trajectories.

flowchart TD
  S["Start: observed nonlinearity in actuator/transmission"] --> Q1["Is there memory (hysteresis)?"]
  Q1 -->|yes| HYS["Use hysteresis operator \n(e.g., backlash/play) \nwith internal state"]
  Q1 -->|no| Q2["Is response bounded (hard limits)?"]
  Q2 -->|yes| SAT["Use saturation model \n(hard clip or smooth approx)"]
  Q2 -->|no| Q3["Is there insensitivity around zero?"]
  Q3 -->|yes| DZ["Use dead-zone \n(deadband) model"]
  Q3 -->|no| PWL["Fit a PWL static map from data \n(lookup / interpolation)"]
  SAT --> PWL
  DZ --> PWL
  HYS --> SIM["Simulate / analyze region-wise dynamics"]
  PWL --> SIM
        

In many systems, these effects appear together: a dead-zone (amplifier deadband) followed by saturation (current/voltage limit) and then backlash (gear mesh). The combination yields non-smooth, sometimes hybrid trajectories, even when the plant is linear.

A PWL approximation replaces a nonlinear static map \( g(u) \) by a function \( \hat{g}(u) \) that is linear on intervals. Choose breakpoints \( u_0 < u_1 < \cdots < u_m \) and values \( g(u_i) \), then define \( \hat{g} \) by linear interpolation:

\[ \hat{g}(u) = \begin{cases} g(u_0), & u < u_0 \\ g(u_i) + \dfrac{g(u_{i+1})-g(u_i)}{u_{i+1}-u_i}\,(u-u_i), & u_i \le u < u_{i+1},\; i=0,\dots,m-1 \\ g(u_m), & u \ge u_m. \end{cases} \]

The flat extrapolation outside \( [u_0,u_m] \) matches many actuator maps (bounded output). On each interval the slope is constant; hence if a nonlinear element is inside a state-space model, the overall dynamics become piecewise-affine.

2.1 From PWL static map to piecewise-affine dynamics

Suppose \( u = \varphi(v) \) and \( \varphi \) is affine in region \( \mathcal{R}_i \): \( \varphi(v) = k_i v + b_i \) for \( v \in \mathcal{R}_i \). Then the plant

\[ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\,(k_i v + b_i), \qquad v = \mathbf{C}\mathbf{x} + \mathbf{D}\,r(t) \]

becomes (in region \( \mathcal{R}_i \)):

\[ \dot{\mathbf{x}} = (\mathbf{A} + \mathbf{B}k_i\mathbf{C})\mathbf{x} + \mathbf{B}k_i\mathbf{D}\,r(t) + \mathbf{B}b_i. \]

Therefore, within each region the dynamics are linear (affine), and switching occurs when the signal \( v \) crosses a breakpoint. Precisely at a breakpoint, the vector field can be non-differentiable; physically, this corresponds to a change in constraint regime.

2.2 Interpolation error bound (university-level result)

If \( g \) is twice continuously differentiable on \( [u_i, u_{i+1}] \) and \( \max_{u \in [u_i,u_{i+1}]} |g''(u)| \le M_2 \), then the error of linear interpolation satisfies:

\[ \max_{u \in [u_i,u_{i+1}]} |g(u) - \hat{g}(u)| \le \frac{M_2}{8}\,(u_{i+1}-u_i)^2. \]

Sketch of proof. Let \( L(u) \) be the line through \( (u_i,g(u_i)) \) and \( (u_{i+1},g(u_{i+1})) \). The interpolation remainder for linear interpolation can be derived from Taylor’s theorem with remainder: for each \( u \) there exists \( \xi \) between \( u_i \) and \( u_{i+1} \) such that \( g(u)-L(u)=\dfrac{g''(\xi)}{2}(u-u_i)(u-u_{i+1}) \). The product \( (u-u_i)(u-u_{i+1}) \) is maximized in magnitude at the midpoint, giving \( \max |(u-u_i)(u-u_{i+1})| = (u_{i+1}-u_i)^2/4 \), and the bound follows.

Saturation represents hard actuator limits. The most common symmetric hard saturation is \( \operatorname{sat}(u;U) \):

\[ \operatorname{sat}(u;U)= \begin{cases} -U, & u < -U \\ u, & -U \le u \le U \\ U, & u > U. \end{cases} \]

3.1 Lipschitz continuity (existence/uniqueness friendly)

Proposition. For all real \( u,v \), \( |\operatorname{sat}(u;U) - \operatorname{sat}(v;U)| \le |u-v| \).

Proof. The map is a composition of “clipping” operations \( \min(\max(u,-U),U) \), each of which is 1-Lipschitz. Alternatively, check cases: if both points are inside \( [-U,U] \), the difference is exactly \( |u-v| \); if one or both are clipped, the output moves no faster than the input because the slope is either 0 or 1. Hence the global Lipschitz constant is 1. □

3.2 Sector bound (useful for later stability tools)

For symmetric saturation, the ratio \( \operatorname{sat}(u;U)/u \) lies between 0 and 1 for \( u \ne 0 \), which yields the sector inequality:

\[ 0 \le \operatorname{sat}(u;U)\,u \le u^2 \quad \text{for all } u\in\mathbb{R}. \]

Intuitively: saturation never flips sign (first inequality) and never amplifies magnitude (second inequality). We will reuse this “bounded gain” viewpoint when we study more global nonlinear stability concepts later.

3.3 Smooth saturation (optional modeling)

In some analyses one uses a differentiable approximation such as \( u_\text{eff} = U\tanh(u/U) \). This is still bounded by \( \pm U \), but avoids non-differentiability at \( \pm U \). A PWL fit can approximate such smooth maps with guaranteed error bounds (Section 2.2).

A dead-zone captures insensitivity near zero: small commands do not move the actuator (e.g., valve stiction threshold, amplifier deadband). A symmetric dead-zone of half-width \( d \) is:

\[ \operatorname{dz}(u;d)= \begin{cases} 0, & |u| \le d \\ u - d\,\operatorname{sgn}(u), & |u| > d. \end{cases} \]

4.1 Basic properties

• Lipschitz: \( |\operatorname{dz}(u;d)-\operatorname{dz}(v;d)| \le |u-v| \) (again slopes are 0 or 1).
• Sector: \( 0 \le \operatorname{dz}(u;d)\,u \le u^2 \) for all \( u \).
• Non-differentiability: kinks at \( u=\pm d \) mean local linearization depends on the operating point: about \( u=0 \) the incremental gain is 0, while outside the dead-zone the incremental gain is 1.

4.2 Why dead-zone matters in feedback loops

In a simple proportional loop, small steady-state errors may fail to generate enough control effort to exit the dead-zone, so the error “sticks” at a nonzero value. This is a purely nonlinear effect and is not predicted by linear models.

Backlash is a gap in mechanical transmissions (e.g., gear mesh). Unlike saturation and dead-zone, backlash is not memoryless: the output depends on the past direction of motion. A widely used mathematical model is the play operator (rate-independent hysteresis) of width \( b \).

5.1 Play/backlash as an inequality + update rule

Let \( u(t) \) be the input (e.g., motor shaft angle or commanded displacement) and \( y(t) \) be the output (e.g., load angle). The play operator enforces the band constraint:

\[ y(t) \in [\,u(t) - b/2,\; u(t) + b/2\,]. \]

Inside the band, the output does not change. When \( u(t) \) pushes the boundary, the output follows to keep the constraint active. One compact (engineering) description is the discrete-time update used in simulation:

\[ y_{k+1} = \begin{cases} u_{k+1} - b/2, & u_{k+1} > y_k + b/2 \\ u_{k+1} + b/2, & u_{k+1} < y_k - b/2 \\ y_k, & |u_{k+1} - y_k| \le b/2. \end{cases} \]

flowchart TD
  IN["Input u(k+1), previous output y(k), gap b"] --> C1["Is u(k+1) > y(k) + b/2 ?"]
  C1 -->|yes| UP1["Set y(k+1) = u(k+1) - b/2"]
  C1 -->|no| C2["Is u(k+1) < y(k) - b/2 ?"]
  C2 -->|yes| UP2["Set y(k+1) = u(k+1) + b/2"]
  C2 -->|no| HOLD["Keep y(k+1) = y(k)"]
  UP1 --> OUT["Output y(k+1)"]
  UP2 --> OUT
  HOLD --> OUT
        

5.2 Interpretation

The model creates a hysteresis loop: for increasing \( u \), the output lags by \( b/2 \), and for decreasing \( u \) it lags in the opposite direction. This captures the key physical fact: reversing direction does not immediately reverse the output because the gear gap must be “taken up”.

The implementations below provide: (i) reusable nonlinearity blocks, (ii) a simple closed-loop simulation that shows the effects, and (iii) CSV output for plotting. In later chapters we will formalize simulation, stiffness, and event handling; here the focus is on the models.

6.1 Python (NumPy/SciPy optional; self-contained RK4 here)

File: Chapter14_Lesson4.py

"""
Chapter14_Lesson4.py
System Dynamics (Control Engineering) — Chapter 14, Lesson 4
Piecewise-Linear Approximations, Saturation, Dead-Zone, and Backlash Models

This script:
1) Defines common nonlinearities used in engineering models (static + hysteretic).
2) Builds a simple closed-loop second-order plant.
3) Simulates it with an explicit RK4 integrator (kept self-contained).
4) Writes results to CSV for plotting/analysis.

Note: Numerical integration is treated formally in Chapter 15; here we use RK4 as a practical tool.
"""

from __future__ import annotations
from dataclasses import dataclass
from typing import List, Tuple
import math
import csv


def sat(u: float, u_max: float) -> float:
    """Hard saturation: clip to [-u_max, u_max]."""
    if u > u_max:
        return u_max
    if u < -u_max:
        return -u_max
    return u


def dead_zone(u: float, d: float) -> float:
    """
    Symmetric dead-zone of half-width d:
      y = 0                if |u| <= d
      y = u - d*sign(u)    if |u| >  d
    """
    if abs(u) <= d:
        return 0.0
    return u - d * (1.0 if u > 0.0 else -1.0)


@dataclass
class Backlash:
    """
    Rate-independent backlash (play operator) of width b.
    The output y stays within a band around the input u:
        y ∈ [u - b/2, u + b/2]
    and only changes when u pushes y to the band boundary.
    """
    b: float
    y: float = 0.0

    def step(self, u: float) -> float:
        half = 0.5 * self.b
        if u > self.y + half:
            self.y = u - half
        elif u < self.y - half:
            self.y = u + half
        return self.y


@dataclass
class PWL:
    """
    1D piecewise-linear interpolant defined by breakpoints (x_i, y_i), x_i strictly increasing.
    """
    xs: List[float]
    ys: List[float]

    def __post_init__(self) -> None:
        if len(self.xs) != len(self.ys) or len(self.xs) < 2:
            raise ValueError("PWL requires >=2 points with equal-length xs and ys.")
        for i in range(len(self.xs) - 1):
            if not (self.xs[i] < self.xs[i+1]):
                raise ValueError("xs must be strictly increasing.")

    def eval(self, x: float) -> float:
        # Extrapolate flat outside range (common in actuator maps).
        if x <= self.xs[0]:
            return self.ys[0]
        if x >= self.xs[-1]:
            return self.ys[-1]

        # Find interval (linear search; for large lists use bisect).
        for i in range(len(self.xs) - 1):
            if self.xs[i] <= x <= self.xs[i+1]:
                x0, x1 = self.xs[i], self.xs[i+1]
                y0, y1 = self.ys[i], self.ys[i+1]
                t = (x - x0) / (x1 - x0)
                return (1.0 - t) * y0 + t * y1

        # Should not reach here
        return self.ys[-1]


def rk4_step(f, t: float, x: Tuple[float, float], dt: float, u: float) -> Tuple[float, float]:
    """One step of RK4 for a 2-state system x' = f(t, x, u)."""
    x1, x2 = x

    k1 = f(t, (x1, x2), u)
    k2 = f(t + 0.5 * dt, (x1 + 0.5 * dt * k1[0], x2 + 0.5 * dt * k1[1]), u)
    k3 = f(t + 0.5 * dt, (x1 + 0.5 * dt * k2[0], x2 + 0.5 * dt * k2[1]), u)
    k4 = f(t + dt, (x1 + dt * k3[0], x2 + dt * k3[1]), u)

    xn1 = x1 + (dt / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0])
    xn2 = x2 + (dt / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
    return (xn1, xn2)


def plant_dynamics(t: float, x: Tuple[float, float], u: float,
                   wn: float = 5.0, zeta: float = 0.25, b: float = 1.0) -> Tuple[float, float]:
    """
    Second-order plant:
        x1' = x2
        x2' = -2*zeta*wn*x2 - wn^2*x1 + b*u
    """
    x1, x2 = x
    dx1 = x2
    dx2 = -2.0 * zeta * wn * x2 - (wn * wn) * x1 + b * u
    return (dx1, dx2)


def simulate(T: float = 6.0, dt: float = 1e-3) -> None:
    # Reference (step)
    def r(t: float) -> float:
        return 1.0

    # PD controller (simple)
    Kp, Kd = 18.0, 4.5

    # Nonlinearities
    u_max = 2.0     # saturation limit
    d = 0.25        # dead-zone half-width
    backlash = Backlash(b=0.20, y=0.0)

    # Optional: approximate a smooth nonlinearity g(u)=tanh(u) using PWL
    # (Illustrates PWL approximation; not required for saturation/deadzone/backlash.)
    xs = [-3.0, -1.5, -0.5, 0.0, 0.5, 1.5, 3.0]
    ys = [math.tanh(x) for x in xs]
    pwl_tanh = PWL(xs=xs, ys=ys)

    use_pwl_map = False  # if True, apply u <- pwl_tanh.eval(u) as a static map

    x = (0.0, 0.0)
    t = 0.0

    rows = []
    while t <= T + 1e-12:
        # Control law
        e = r(t) - x[0]
        u_cmd = Kp * e - Kd * x[1]

        # Apply dead-zone (e.g., valve stiction / quantization / amplifier deadband)
        u_dz = dead_zone(u_cmd, d=d)

        # Apply saturation (actuator limit)
        u_sat = sat(u_dz, u_max=u_max)

        # Optional static PWL map (e.g., linearizing an empirically measured nonlinearity)
        u_map = pwl_tanh.eval(u_sat) if use_pwl_map else u_sat

        # Apply backlash (gear train); introduces hysteresis (memory)
        u_bl = backlash.step(u_map)

        # Integrate plant
        x = rk4_step(lambda tt, xx, uu: plant_dynamics(tt, xx, uu), t, x, dt, u_bl)

        rows.append((t, x[0], x[1], u_cmd, u_dz, u_sat, u_bl))
        t += dt

    out_csv = "Chapter14_Lesson4_output.csv"
    with open(out_csv, "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["t", "x1", "x2", "u_cmd", "u_deadzone", "u_sat", "u_backlash"])
        w.writerows(rows)

    print(f"Wrote {out_csv} with {len(rows)} samples.")
    print("Tip: plot x1(t) and compare u_cmd vs u_backlash to see dead-zone + saturation + hysteresis effects.")


if __name__ == "__main__":
    simulate()
      

6.2 C++ (standard library; optional Eigen for larger models)

File: Chapter14_Lesson4.cpp

/*
Chapter14_Lesson4.cpp
System Dynamics (Control Engineering) — Chapter 14, Lesson 4
Piecewise-Linear Approximations, Saturation, Dead-Zone, and Backlash Models

Build (example, MSYS2/MinGW64):
  g++ -O2 -std=c++17 Chapter14_Lesson4.cpp -o Chapter14_Lesson4.exe

Run:
  ./Chapter14_Lesson4.exe

It writes: Chapter14_Lesson4_output.csv
*/

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <stdexcept>
#include <iomanip>
#include <string>

static double sat(double u, double umax) {
    if (u > umax) return umax;
    if (u < -umax) return -umax;
    return u;
}

static double dead_zone(double u, double d) {
    if (std::abs(u) <= d) return 0.0;
    return u - d * ((u > 0.0) ? 1.0 : -1.0);
}

struct Backlash {
    double b;
    double y;
    explicit Backlash(double width, double y0=0.0) : b(width), y(y0) {}

    double step(double u) {
        const double half = 0.5 * b;
        if (u > y + half) y = u - half;
        else if (u < y - half) y = u + half;
        return y;
    }
};

struct PWL {
    std::vector<double> xs;
    std::vector<double> ys;

    PWL(std::vector<double> x, std::vector<double> y) : xs(std::move(x)), ys(std::move(y)) {
        if (xs.size() != ys.size() || xs.size() < 2) throw std::runtime_error("PWL needs >=2 points.");
        for (size_t i=0; i+1<xs.size(); ++i) {
            if (!(xs[i] < xs[i+1])) throw std::runtime_error("PWL xs must be strictly increasing.");
        }
    }

    double eval(double x) const {
        if (x <= xs.front()) return ys.front();
        if (x >= xs.back())  return ys.back();
        for (size_t i=0; i+1<xs.size(); ++i) {
            if (xs[i] <= x && x <= xs[i+1]) {
                const double x0 = xs[i], x1 = xs[i+1];
                const double y0 = ys[i], y1 = ys[i+1];
                const double t = (x - x0) / (x1 - x0);
                return (1.0 - t)*y0 + t*y1;
            }
        }
        return ys.back();
    }
};

static void plant(double /*t*/, const double x1, const double x2, const double u,
                  double& dx1, double& dx2, double wn=5.0, double zeta=0.25, double b=1.0) {
    dx1 = x2;
    dx2 = -2.0*zeta*wn*x2 - (wn*wn)*x1 + b*u;
}

static void rk4_step(double t, double& x1, double& x2, double dt, double u) {
    double k1_1, k1_2; plant(t, x1, x2, u, k1_1, k1_2);
    double k2_1, k2_2; plant(t+0.5*dt, x1+0.5*dt*k1_1, x2+0.5*dt*k1_2, u, k2_1, k2_2);
    double k3_1, k3_2; plant(t+0.5*dt, x1+0.5*dt*k2_1, x2+0.5*dt*k2_2, u, k3_1, k3_2);
    double k4_1, k4_2; plant(t+dt, x1+dt*k3_1, x2+dt*k3_2, u, k4_1, k4_2);

    x1 += (dt/6.0)*(k1_1 + 2.0*k2_1 + 2.0*k3_1 + k4_1);
    x2 += (dt/6.0)*(k1_2 + 2.0*k2_2 + 2.0*k3_2 + k4_2);
}

int main() {
    const double T = 6.0;
    const double dt = 1e-3;

    // PD controller
    const double Kp = 18.0, Kd = 4.5;

    // Nonlinearities
    const double u_max = 2.0;
    const double d = 0.25;
    Backlash backlash(0.20, 0.0);

    // Optional PWL map approximating tanh(u)
    std::vector<double> xs = {-3.0,-1.5,-0.5,0.0,0.5,1.5,3.0};
    std::vector<double> ys; ys.reserve(xs.size());
    for (double x : xs) ys.push_back(std::tanh(x));
    PWL pwl_tanh(xs, ys);

    const bool use_pwl_map = false;

    double x1 = 0.0, x2 = 0.0;

    std::ofstream out("Chapter14_Lesson4_output.csv");
    out << "t,x1,x2,u_cmd,u_deadzone,u_sat,u_backlash\n";
    out << std::fixed << std::setprecision(9);

    for (double t=0.0; t <= T + 1e-12; t += dt) {
        const double r = 1.0;               // step reference
        const double e = r - x1;
        const double u_cmd = Kp*e - Kd*x2;

        const double u_dz  = dead_zone(u_cmd, d);
        const double u_sat = sat(u_dz, u_max);
        const double u_map = use_pwl_map ? pwl_tanh.eval(u_sat) : u_sat;
        const double u_bl  = backlash.step(u_map);

        rk4_step(t, x1, x2, dt, u_bl);

        out << t << "," << x1 << "," << x2 << ","
            << u_cmd << "," << u_dz << "," << u_sat << "," << u_bl << "\n";
    }

    std::cout << "Wrote Chapter14_Lesson4_output.csv\n";
    std::cout << "Plot x1(t) and compare u_cmd vs u_backlash to see nonlinear effects.\n";
    return 0;
}
      

6.3 Java (standard library; optional Apache Commons Math)

File: Chapter14_Lesson4.java

/*
Chapter14_Lesson4.java
System Dynamics (Control Engineering) — Chapter 14, Lesson 4
Piecewise-Linear Approximations, Saturation, Dead-Zone, and Backlash Models

Compile:
  javac Chapter14_Lesson4.java

Run:
  java Chapter14_Lesson4

Output:
  Chapter14_Lesson4_output.csv
*/

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

public class Chapter14_Lesson4 {

    static double sat(double u, double umax) {
        if (u > umax) return umax;
        if (u < -umax) return -umax;
        return u;
    }

    static double deadZone(double u, double d) {
        if (Math.abs(u) <= d) return 0.0;
        return u - d * ((u > 0.0) ? 1.0 : -1.0);
    }

    static class Backlash {
        final double b;
        double y;
        Backlash(double width, double y0) { b = width; y = y0; }
        double step(double u) {
            final double half = 0.5 * b;
            if (u > y + half) y = u - half;
            else if (u < y - half) y = u + half;
            return y;
        }
    }

    static class PWL {
        final double[] xs;
        final double[] ys;

        PWL(double[] xs, double[] ys) {
            if (xs.length != ys.length || xs.length < 2) {
                throw new IllegalArgumentException("PWL needs >=2 points with equal lengths.");
            }
            for (int i=0; i<xs.length-1; i++) {
                if (!(xs[i] < xs[i+1])) throw new IllegalArgumentException("xs must be strictly increasing.");
            }
            this.xs = Arrays.copyOf(xs, xs.length);
            this.ys = Arrays.copyOf(ys, ys.length);
        }

        double eval(double x) {
            if (x <= xs[0]) return ys[0];
            if (x >= xs[xs.length-1]) return ys[ys.length-1];
            for (int i=0; i<xs.length-1; i++) {
                if (xs[i] <= x && x <= xs[i+1]) {
                    double x0 = xs[i], x1 = xs[i+1];
                    double y0 = ys[i], y1 = ys[i+1];
                    double t = (x - x0) / (x1 - x0);
                    return (1.0 - t)*y0 + t*y1;
                }
            }
            return ys[ys.length-1];
        }
    }

    static double[] plant(double x1, double x2, double u, double wn, double zeta, double b) {
        double dx1 = x2;
        double dx2 = -2.0*zeta*wn*x2 - (wn*wn)*x1 + b*u;
        return new double[]{dx1, dx2};
    }

    static double[] rk4Step(double x1, double x2, double u, double dt) {
        double wn = 5.0, zeta = 0.25, b = 1.0;

        double[] k1 = plant(x1, x2, u, wn, zeta, b);
        double[] k2 = plant(x1 + 0.5*dt*k1[0], x2 + 0.5*dt*k1[1], u, wn, zeta, b);
        double[] k3 = plant(x1 + 0.5*dt*k2[0], x2 + 0.5*dt*k2[1], u, wn, zeta, b);
        double[] k4 = plant(x1 + dt*k3[0], x2 + dt*k3[1], u, wn, zeta, b);

        double xn1 = x1 + (dt/6.0)*(k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]);
        double xn2 = x2 + (dt/6.0)*(k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1]);

        return new double[]{xn1, xn2};
    }

    public static void main(String[] args) throws IOException {
        final double T = 6.0;
        final double dt = 1e-3;

        // PD controller
        final double Kp = 18.0, Kd = 4.5;

        // Nonlinearities
        final double uMax = 2.0;
        final double d = 0.25;
        Backlash backlash = new Backlash(0.20, 0.0);

        // Optional PWL map approximating tanh(u)
        double[] xs = new double[]{-3.0,-1.5,-0.5,0.0,0.5,1.5,3.0};
        double[] ys = new double[xs.length];
        for (int i=0; i<xs.length; i++) ys[i] = Math.tanh(xs[i]);
        PWL pwlTanh = new PWL(xs, ys);

        final boolean usePwlMap = false;

        double x1 = 0.0, x2 = 0.0;

        try (FileWriter fw = new FileWriter("Chapter14_Lesson4_output.csv")) {
            fw.write("t,x1,x2,u_cmd,u_deadzone,u_sat,u_backlash\n");
            for (double t=0.0; t <= T + 1e-12; t += dt) {
                double r = 1.0;
                double e = r - x1;
                double uCmd = Kp*e - Kd*x2;

                double uDz  = deadZone(uCmd, d);
                double uSat = sat(uDz, uMax);
                double uMap = usePwlMap ? pwlTanh.eval(uSat) : uSat;
                double uBl  = backlash.step(uMap);

                double[] xn = rk4Step(x1, x2, uBl, dt);
                x1 = xn[0]; x2 = xn[1];

                fw.write(String.format(java.util.Locale.US,
                        "%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f\n",
                        t, x1, x2, uCmd, uDz, uSat, uBl));
            }
        }

        System.out.println("Wrote Chapter14_Lesson4_output.csv");
        System.out.println("Plot x1(t) and compare u_cmd vs u_backlash to see nonlinear effects.");
    }
}
      

6.4 MATLAB / Simulink

File: Chapter14_Lesson4.m

% Chapter14_Lesson4.m
% System Dynamics (Control Engineering) — Chapter 14, Lesson 4
% Piecewise-Linear Approximations, Saturation, Dead-Zone, and Backlash Models
%
% This script implements the nonlinearities and runs a simple RK4 simulation
% of a 2nd-order plant under a PD controller. Results are written to CSV.
%
% Note: Numerical simulation methods are covered systematically in Chapter 15.
% Here RK4 is used as a practical, self-contained integrator.

clear; clc;

T  = 6.0;
dt = 1e-3;
t  = 0:dt:T;

% PD ←controller
Kp = 18.0; Kd = 4.5;

% Nonlinearities
u_max = 2.0;      % saturation
d     = 0.25;     % dead-zone half width
b_gap = 0.20;     % backlash width

% Plant parameters
wn = 5.0; zeta = 0.25; b = 1.0;

% Backlash internal state
y_bl = 0.0;

% Optional PWL map approximating tanh(u)
xs = [-3.0 -1.5 -0.5 0.0 0.5 1.5 3.0];
ys = tanh(xs);
use_pwl_map = false;

% Storage
x1 = zeros(size(t));  % position
x2 = zeros(size(t));  % velocity
u_cmd = zeros(size(t));
u_dz  = zeros(size(t));
u_sat = zeros(size(t));
u_bl  = zeros(size(t));

for k = 1:(numel(t)-1)
    r = 1.0; % step reference

    e = r - x1(k);
    u_cmd(k) = Kp*e - Kd*x2(k);

    u_dz(k)  = dead_zone(u_cmd(k), d);
    u_sat(k) = sat(u_dz(k), u_max);

    u_map = u_sat(k);
    if use_pwl_map
        u_map = pwl_eval(xs, ys, u_sat(k));
    end

    [u_bl(k), y_bl] = backlash_step(u_map, y_bl, b_gap);

    % RK4 integration for x' = f(x,u)
    x = [x1(k); x2(k)];
    f = @(xx,uu) plant(xx, uu, wn, zeta, b);

    k1 = f(x, u_bl(k));
    k2 = f(x + 0.5*dt*k1, u_bl(k));
    k3 = f(x + 0.5*dt*k2, u_bl(k));
    k4 = f(x + dt*k3, u_bl(k));

    x_next = x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);
    x1(k+1) = x_next(1);
    x2(k+1) = x_next(2);
end

% Last control sample (for completeness)
u_cmd(end) = u_cmd(end-1);
u_dz(end)  = u_dz(end-1);
u_sat(end) = u_sat(end-1);
u_bl(end)  = u_bl(end-1);

tbl = table(t(:), x1(:), x2(:), u_cmd(:), u_dz(:), u_sat(:), u_bl(:), ...
    'VariableNames', {'t','x1','x2','u_cmd','u_deadzone','u_sat','u_backlash'});

writetable(tbl, 'Chapter14_Lesson4_output.csv');
disp('Wrote Chapter14_Lesson4_output.csv');

% Optional quick plots (uncomment if desired)
% figure; plot(t, x1); grid on; xlabel('t'); ylabel('x1');
% figure; plot(t, u_cmd, t, u_bl); grid on; xlabel('t'); ylabel('u'); legend('u\_cmd','u\_backlash');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonlinearity functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = sat(u, umax)
    y = min(max(u, -umax), umax);
end

function y = dead_zone(u, d)
    if abs(u) <= d
        y = 0.0;
    else
        y = u - d*sign(u);
    end
end

function [y, y_state] = backlash_step(u, y_state, b_gap)
    half = 0.5*b_gap;
    if u > y_state + half
        y_state = u - half;
    elseif u < y_state - half
        y_state = u + half;
    end
    y = y_state;
end

function y = pwl_eval(xs, ys, x)
    % Flat extrapolation outside range
    if x <= xs(1), y = ys(1); return; end
    if x >= xs(end), y = ys(end); return; end
    % Find interval
    idx = find(xs <= x, 1, 'last');
    x0 = xs(idx); x1 = xs(idx+1);
    y0 = ys(idx); y1 = ys(idx+1);
    t = (x - x0)/(x1 - x0);
    y = (1-t)*y0 + t*y1;
end

function dx = plant(x, u, wn, zeta, b)
    % x = [x1; x2]
    dx = zeros(2,1);
    dx(1) = x(2);
    dx(2) = -2*zeta*wn*x(2) - (wn^2)*x(1) + b*u;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulink note (blocks):
% - Use "Saturation" block for sat()
% - Use "Dead Zone" block for dead_zone()
% - Use "Backlash" block for backlash_step()
% - Use "1-D Lookup Table" for PWL maps (xs, ys)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      

6.5 Wolfram Mathematica

File: Chapter14_Lesson4.nb

(* ::Package:: *)

(* Chapter14_Lesson4.nb
   System Dynamics (Control Engineering) — Chapter 14, Lesson 4
   Piecewise-Linear Approximations, Saturation, Dead-Zone, and Backlash Models *)

(* This notebook-style file is a plain-text Wolfram Notebook expression.
   You can open it in Wolfram Mathematica. *)

Notebook[{
  Cell["Chapter 14 — Lesson 4: PWL, Saturation, Dead-Zone, Backlash", "Title"],
  Cell["1) Nonlinearities (static + backlash/play operator)", "Section"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    sat[u_, umax_] := Clip[u, {-umax, umax}]
  ], "Input"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    deadZone[u_, d_] := Piecewise[{
      {0, Abs[u] <= d},
      {u - d Sign[u], Abs[u] > d}
    }]
  ], "Input"],

  Cell[TextData[{
    "Backlash as a play operator: output y stays within [u - b/2, u + b/2] ",
    "and only updates when u pushes the band boundary."
  }], "Text"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    backlashStep[u_, y_, b_] := Module[{half = b/2},
      Which[
        u > y + half, {u - half, u - half},
        u < y - half, {u + half, u + half},
        True,         {y, y}
      ]
    ]
  ], "Input"],

  Cell["2) Piecewise-linear (PWL) map via linear interpolation", "Section"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    pwlEval[xs_List, ys_List, x_] := Module[{n = Length[xs], i, t},
      If[x <= xs[[1]], Return[ys[[1]]]];
      If[x >= xs[[n]], Return[ys[[n]]]];
      i = First@FirstPosition[UnitStep[xs - x], 1] - 1; (* xs[[i]] <= x < xs[[i+1]] *)
      t = (x - xs[[i]])/(xs[[i + 1]] - xs[[i]]);
      (1 - t) ys[[i]] + t ys[[i + 1]]
    ]
  ], "Input"],

  Cell["3) Closed-loop simulation (RK4, self-contained)", "Section"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    plant[{x1_, x2_}, u_, wn_, zeta_, b_] := {x2, -2 zeta wn x2 - wn^2 x1 + b u}
  ], "Input"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    rk4Step[f_, x_, u_, dt_] := Module[{k1, k2, k3, k4},
      k1 = f[x, u];
      k2 = f[x + (dt/2) k1, u];
      k3 = f[x + (dt/2) k2, u];
      k4 = f[x + dt k3, u];
      x + (dt/6) (k1 + 2 k2 + 2 k3 + k4)
    ]
  ], "Input"],

  Cell[BoxData @ ToBoxes @ HoldForm[
    Module[{
      T = 6.0, dt = 0.001, steps, wn = 5.0, zeta = 0.25, b = 1.0,
      Kp = 18.0, Kd = 4.5, umax = 2.0, d = 0.25, bgap = 0.20,
      xs, ys, usePwl = False,
      x = {0.0, 0.0}, ybl = 0.0,
      data = {}
    },
      steps = Round[T/dt];
      xs = {-3.0, -1.5, -0.5, 0.0, 0.5, 1.5, 3.0};
      ys = Tanh /@ xs;

      Do[
        Module[{t = k dt, r = 1.0, e, ucmd, udz, usat, umap, ubl, tmp},
          e = r - x[[1]];
          ucmd = Kp e - Kd x[[2]];
          udz = deadZone[ucmd, d];
          usat = sat[udz, umax];
          umap = If[usePwl, pwlEval[xs, ys, usat], usat];
          tmp = backlashStep[umap, ybl, bgap];
          ubl = tmp[[1]]; ybl = tmp[[2]];

          x = rk4Step[(plant[#, #2, wn, zeta, b] &), x, ubl, dt];
          AppendTo[data, {t, x[[1]], x[[2]], ucmd, udz, usat, ubl}];
        ],
        {k, 0, steps}
      ];

      Export["Chapter14_Lesson4_output.csv",
        Prepend[data, {"t","x1","x2","u_cmd","u_deadzone","u_sat","u_backlash"}]
      ];
      "Wrote Chapter14_Lesson4_output.csv"
    ]
  ], "Input"]
}]
      

Problem 1 (Interpolation error of PWL): Let \( g \in C^2([a,b]) \) and assume \( \max_{u\in[a,b]} |g''(u)| \le M_2 \). Show that the linear interpolant \( \hat{g} \) on \( [a,b] \) satisfies \( \max_{u\in[a,b]} |g(u)-\hat{g}(u)| \le M_2(b-a)^2/8 \).

Solution: Use the interpolation remainder formula \( g(u)-\hat{g}(u)=\dfrac{g''(\xi)}{2}(u-a)(u-b) \) for some \( \xi \in (a,b) \). Since \( |g''(\xi)| \le M_2 \) and \( \max_{u\in[a,b]} |(u-a)(u-b)|=(b-a)^2/4 \), the bound follows.

Problem 2 (Saturation is 1-Lipschitz): Prove that for all \( u,v \), \( |\operatorname{sat}(u;U)-\operatorname{sat}(v;U)| \le |u-v| \).

Solution: Write \( \operatorname{sat}(u;U)=\min(\max(u,-U),U) \). The map \( u \mapsto \max(u,-U) \) is nonexpansive (it is a projection onto a half-line), and \( u \mapsto \min(u,U) \) is also nonexpansive. The composition of nonexpansive maps is nonexpansive, hence 1-Lipschitz.

Problem 3 (Equilibria with saturation): Consider the scalar system \( \dot{x} = -a x + \operatorname{sat}(k x;U) \) with \( a > 0 \), \( k > 0 \). (a) Find equilibria. (b) Determine stability using piecewise linearization.

Solution: Equilibria solve \( a x^* = \operatorname{sat}(k x^*;U) \). If \( |k x^*| \le U \) then \( a x^* = k x^* \Rightarrow x^*=0 \) unless \( a=k \) (non-isolated case). If \( k x^* > U \), then \( a x^* = U \Rightarrow x^* = U/a \) provided \( k(U/a) > U \Rightarrow k/a > 1 \). Similarly, if \( k x^* < -U \), then \( x^*=-U/a \) provided \( k/a > 1 \). For stability: around \( x=0 \) (unsaturated region), linearization is \( \dot{x} \approx -(a-k)x \). Hence \( x=0 \) is locally asymptotically stable if \( a-k > 0 \), unstable if \( a-k < 0 \). Around \( x^*=\pm U/a \) the saturation is flat (incremental gain 0), so linearization is \( \dot{\delta x}= -a\,\delta x \), which is stable. Thus when \( k/a > 1 \), stable saturated equilibria appear at \( \pm U/a \).

Problem 4 (Dead-zone steady-state error intuition): Consider a proportional control law \( u = K(r-x) \) applied to the stable first-order plant \( \dot{x} = -a x + u_\text{eff} \), where \( u_\text{eff} = \operatorname{dz}(u;d) \). Explain why a constant reference \( r \) may lead to a nonzero steady-state tracking error.

Solution: If the loop settles such that \( |u|=|K(r-x)| \le d \), then \( u_\text{eff}=0 \) and the plant converges to \( x \to 0 \), not \( r \). To maintain \( x \approx r \), the controller must produce a command outside the dead-zone. Therefore, the equilibrium condition becomes piecewise: either \( x^*=0 \) (inside dead-zone) or \( a x^* = K(r-x^*) - d\,\operatorname{sgn}(K(r-x^*)) \) (outside dead-zone). For small references, the first regime is consistent and yields nonzero error.

Problem 5 (Backlash band invariance): Using the update rule in Section 5.1, prove that if \( y_k \in [u_k-b/2,\;u_k+b/2] \), then the next output satisfies \( y_{k+1} \in [u_{k+1}-b/2,\;u_{k+1}+b/2] \) for any input sequence \( u_k \).

Solution: Consider three cases. (i) If \( u_{k+1} > y_k + b/2 \), the rule sets \( y_{k+1}=u_{k+1}-b/2 \), which lies exactly on the lower band boundary. (ii) If \( u_{k+1} < y_k - b/2 \), then \( y_{k+1}=u_{k+1}+b/2 \), which lies exactly on the upper boundary. (iii) Otherwise \( |u_{k+1}-y_k| \le b/2 \), and we keep \( y_{k+1}=y_k \), which immediately implies \( y_{k+1} \in [u_{k+1}-b/2,\;u_{k+1}+b/2] \). Hence the band constraint is invariant. □

We formalized four ubiquitous engineering nonlinearities: PWL static maps (region-wise linearization), saturation (hard limits), dead-zone (insensitivity), and backlash (hysteresis via an internal-state operator). We derived key mathematical properties (Lipschitz continuity, sector-type bounds, interpolation error bounds) that support well-posed simulation and later stability reasoning. The accompanying implementations provide reusable building blocks and a closed-loop example to observe their qualitative impact on trajectories.