In previous lessons of this chapter you modeled mechanical, electrical, and electromechanical systems, mostly assuming linear constitutive laws (Hooke's law springs, linear resistors, ideal op-amps, etc.). Real actuators, sensors, and robot links, however, often contain nonlinear elements such as saturations, dead zones, Coulomb friction, diode characteristics, and trigonometric terms like \( \sin(\theta) \) in pendulum dynamics.

At the same time, the analysis and design methods in upcoming chapters (time response, stability criteria, root locus, frequency response, PID tuning, etc.) are based on linear models. The key idea of this lesson is:

Replace a nonlinear model by a linear LTI approximation that is accurate for small deviations around a chosen operating point.

A simple illustration is the trigonometric approximation near the origin:

\[ \sin(\theta) = \theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots \approx \theta \quad \text{for small } |\theta|. \]

The linear term \( \theta \) captures the dominant behavior close to \( \theta = 0 \), while higher-order terms are small. Linearization of dynamical systems generalizes this idea: we approximate the nonlinear vector field by its first-order Taylor expansion around a point of interest.

flowchart TD
  P["Physical system (mechanical / electrical / robotic)"]
    --> N["Nonlinear equations of motion"]
  N --> O["Choose operating point (steady state)"]
  O --> L["Linearized LTI model near operating point"]
  L --> A["Use linear analysis and design tools (later chapters)"]
        

We consider time-invariant nonlinear systems written in state-space form (you can think of the state as the minimal set of variables needed to describe the dynamics):

\[ \dot{\mathbf{x}}(t) = \mathbf{f}\big(\mathbf{x}(t),\,\mathbf{u}(t)\big), \qquad \mathbf{y}(t) = \mathbf{g}\big(\mathbf{x}(t),\,\mathbf{u}(t)\big), \]

where \( \mathbf{x}(t) \in \mathbb{R}^n \) is the state vector, \( \mathbf{u}(t) \in \mathbb{R}^m \) the input vector, and \( \mathbf{y}(t) \in \mathbb{R}^p \) the output vector. The functions \( \mathbf{f} \) and \( \mathbf{g} \) are (in general) nonlinear, but assumed to be continuously differentiable so that we can apply Taylor expansions.

2.1 Equilibrium (steady-state) points

A constant input \( \mathbf{u}(t) = \mathbf{u}^\ast \) and constant state \( \mathbf{x}(t) = \mathbf{x}^\ast \) form an equilibrium point (or steady-state operating point) if

\[ \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) = \mathbf{0}. \]

This means that if the system is placed exactly at \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \), it will remain there forever (no motion), because the time derivative of the state is zero.

2.2 Operating point vs. equilibrium

In many control applications, we linearize around an equilibrium point: the operating point is exactly the equilibrium. Later in more advanced courses you may encounter linearization around non-equilibrium trajectories (time-varying operating points), but in this introductory view we focus on the equilibrium case, which already covers many important problems such as robot joint control about a fixed posture or motor speed regulation around a constant velocity.

To understand the general multivariable case, it is helpful to start with a scalar system:

\[ \dot{x}(t) = f\big(x(t),u(t)\big), \]

with scalar state \( x \) and scalar input \( u \). Assume \( (x^\ast,u^\ast) \) is an equilibrium, i.e. \( f(x^\ast,u^\ast) = 0 \). We introduce perturbation variables:

\[ \delta x(t) = x(t) - x^\ast, \qquad \delta u(t) = u(t) - u^\ast. \]

Since \( x^\ast \) is constant, \( \delta \dot{x}(t) = \dot{x}(t) \). We expand \( f(x,u) \) around the operating point using the first-order Taylor series in two variables:

\[ f(x,u) \approx f(x^\ast,u^\ast) + \frac{\partial f}{\partial x}\bigg|_{(x^\ast,u^\ast)} (x - x^\ast) + \frac{\partial f}{\partial u}\bigg|_{(x^\ast,u^\ast)} (u - u^\ast) + \text{higher-order terms}. \]

Because \( f(x^\ast,u^\ast) = 0 \), the constant term vanishes, and in terms of perturbations:

\[ \delta \dot{x}(t) \approx a\, \delta x(t) + b\, \delta u(t), \]

where

\[ a = \frac{\partial f}{\partial x}\bigg|_{(x^\ast,u^\ast)}, \qquad b = \frac{\partial f}{\partial u}\bigg|_{(x^\ast,u^\ast)}. \]

The omitted higher-order terms are of order \( \mathcal{O}\big(\delta x^2,\,\delta x\,\delta u,\,\delta u^2\big) \), so if \( \delta x \) and \( \delta u \) are sufficiently small, the linear approximation is accurate. This is exactly the regime in which classical linear control is most effective.

Consider now the general nonlinear system

\[ \dot{\mathbf{x}}(t) = \mathbf{f}\big(\mathbf{x}(t),\mathbf{u}(t)\big), \qquad \mathbf{y}(t) = \mathbf{g}\big(\mathbf{x}(t),\mathbf{u}(t)\big), \]

with \( \mathbf{x}(t) \in \mathbb{R}^n \), \( \mathbf{u}(t) \in \mathbb{R}^m \), \( \mathbf{y}(t) \in \mathbb{R}^p \). Let \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) be an equilibrium:

\[ \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) = \mathbf{0}. \]

Define perturbations \( \delta \mathbf{x}(t) = \mathbf{x}(t) - \mathbf{x}^\ast \), \( \delta \mathbf{u}(t) = \mathbf{u}(t) - \mathbf{u}^\ast \), and similarly \( \delta \mathbf{y}(t) = \mathbf{y}(t) - \mathbf{y}^\ast \) with \( \mathbf{y}^\ast = \mathbf{g}(\mathbf{x}^\ast,\mathbf{u}^\ast) \).

The multivariable Taylor expansion (keeping only first-order terms) gives

\[ \begin{aligned} \dot{\mathbf{x}}(t) &= \mathbf{f}\big(\mathbf{x}^\ast + \delta \mathbf{x}(t), \mathbf{u}^\ast + \delta \mathbf{u}(t)\big) \\ &\approx \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \underbrace{\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\bigg|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}}_{\mathbf{A}} \delta \mathbf{x}(t) + \underbrace{\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\bigg|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}}_{\mathbf{B}} \delta \mathbf{u}(t). \end{aligned} \]

Using \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)=\mathbf{0} \) and \( \delta \dot{\mathbf{x}}(t) = \dot{\mathbf{x}}(t) \), we obtain the linearized state equation:

\[ \delta \dot{\mathbf{x}}(t) = \mathbf{A}\,\delta \mathbf{x}(t) + \mathbf{B}\,\delta \mathbf{u}(t), \]

where

\[ \mathbf{A} = \left[ \frac{\partial f_i}{\partial x_j}(\mathbf{x}^\ast,\mathbf{u}^\ast) \right]_{i=1,\dots,n}^{j=1,\dots,n}, \qquad \mathbf{B} = \left[ \frac{\partial f_i}{\partial u_k}(\mathbf{x}^\ast,\mathbf{u}^\ast) \right]_{i=1,\dots,n}^{k=1,\dots,m}. \]

Similarly, linearizing the output equation:

\[ \begin{aligned} \mathbf{y}(t) &= \mathbf{g}\big(\mathbf{x}^\ast + \delta \mathbf{x}(t), \mathbf{u}^\ast + \delta \mathbf{u}(t)\big) \\ &\approx \mathbf{g}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \underbrace{\frac{\partial \mathbf{g}}{\partial \mathbf{x}}\bigg|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}}_{\mathbf{C}} \delta \mathbf{x}(t) + \underbrace{\frac{\partial \mathbf{g}}{\partial \mathbf{u}}\bigg|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}}_{\mathbf{D}} \delta \mathbf{u}(t), \end{aligned} \]

yields the linearized output equation:

\[ \delta \mathbf{y}(t) = \mathbf{C}\,\delta \mathbf{x}(t) + \mathbf{D}\,\delta \mathbf{u}(t). \]

The matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) are constants (time-invariant) as long as we linearize around a constant equilibrium of a time-invariant nonlinear system. Thus the linearized model is an LTI system in the perturbation variables.

Consider a simple pendulum of mass \( m \) and length \( L \), with viscous damping \( b \) and torque input \( u(t) \) applied at the pivot. Let \( \theta(t) \) be the angle measured from the downward vertical. The equation of motion is

\[ mL^2 \ddot{\theta}(t) + b \dot{\theta}(t) + mgL \sin\big(\theta(t)\big) = u(t). \]

This model is directly analogous to a single-link robot joint with gravity and viscous friction. We define the state variables \( x_1(t) = \theta(t) \) and \( x_2(t) = \dot{\theta}(t) \). Writing the system as first-order equations:

\[ \begin{aligned} \dot{x}_1(t) &= x_2(t), \\ \dot{x}_2(t) &= -\frac{b}{mL^2} x_2(t) - \frac{g}{L} \sin\big(x_1(t)\big) + \frac{1}{mL^2} u(t). \end{aligned} \]

We choose the equilibrium corresponding to the pendulum hanging straight down with zero torque: \( x_1^\ast = 0,\, x_2^\ast = 0,\, u^\ast = 0 \).

5.1 Computing the Jacobians

Let \( \mathbf{x} = [x_1\; x_2]^\top \). Define

\[ \mathbf{f}(\mathbf{x},u) = \begin{bmatrix} f_1(\mathbf{x},u) \\[4pt] f_2(\mathbf{x},u) \end{bmatrix} = \begin{bmatrix} x_2 \\ -\dfrac{b}{mL^2} x_2 - \dfrac{g}{L} \sin(x_1) + \dfrac{1}{mL^2} u \end{bmatrix}. \]

The state Jacobian \( \mathbf{A} = \partial \mathbf{f}/\partial \mathbf{x} \) is

\[ \mathbf{A}(\mathbf{x},u) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2} \\[4pt] \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\dfrac{g}{L}\cos(x_1) & -\dfrac{b}{mL^2} \end{bmatrix}. \]

The input Jacobian \( \mathbf{B} = \partial \mathbf{f}/\partial u \) is

\[ \mathbf{B}(\mathbf{x},u) = \begin{bmatrix} \dfrac{\partial f_1}{\partial u} \\[4pt] \dfrac{\partial f_2}{\partial u} \end{bmatrix} = \begin{bmatrix} 0 \\[4pt] \dfrac{1}{mL^2} \end{bmatrix}. \]

Evaluating at the equilibrium \( (x_1^\ast,x_2^\ast,u^\ast) = (0,0,0) \):

\[ \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -\dfrac{g}{L} & -\dfrac{b}{mL^2} \end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} 0 \\[4pt] \dfrac{1}{mL^2} \end{bmatrix}. \]

The linearized model in perturbation coordinates \( \delta \mathbf{x}(t) = \mathbf{x}(t) - \mathbf{x}^\ast \), \( \delta u(t) = u(t) - u^\ast \) is

\[ \delta \dot{\mathbf{x}}(t) = \begin{bmatrix} 0 & 1 \\ -\dfrac{g}{L} & -\dfrac{b}{mL^2} \end{bmatrix} \delta \mathbf{x}(t) + \begin{bmatrix} 0 \\[4pt] \dfrac{1}{mL^2} \end{bmatrix} \delta u(t). \]

This is an LTI model of a linearized pendulum, which behaves like a damped second-order oscillator near the vertical configuration. In robotics, multi-link manipulators are often linearized in the same way around a nominal joint configuration.

Given a nonlinear model of the form \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},\mathbf{u}) \), the procedure to obtain an LTI approximation around an equilibrium point is:

1. Choose state variables \( \mathbf{x} \) and inputs \( \mathbf{u} \).
2. Find a constant pair \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) satisfying \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) = \mathbf{0} \).
3. Compute the Jacobian matrices \( \mathbf{A} = \partial \mathbf{f}/\partial \mathbf{x} \) and \( \mathbf{B} = \partial \mathbf{f}/\partial \mathbf{u} \) at \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \).
4. (Optionally) linearize the output equation to obtain \( \mathbf{C} \) and \( \mathbf{D} \).
5. Form the local LTI model \( \delta \dot{\mathbf{x}} = \mathbf{A}\,\delta \mathbf{x} + \mathbf{B}\,\delta \mathbf{u} \), \( \delta \mathbf{y} = \mathbf{C}\,\delta \mathbf{x} + \mathbf{D}\,\delta \mathbf{u} \).

flowchart TD
  N1["Start from nonlinear model: dx/dt = f(x,u)"]
    --> N2["Choose candidate operating point (x_op,u_op)"]
  N2 --> N3["Solve f(x_op,u_op) = 0 (equilibrium check)"]
  N3 --> N4["Define perturbations: delta_x = x - x_op, delta_u = u - u_op"]
  N4 --> N5["Compute Jacobians A = df/dx, B = df/du at (x_op,u_op)"]
  N5 --> N6["Obtain linear model: delta_x_dot = A delta_x + B delta_u"]
        

In Python, symbolic linearization can be efficiently performed with sympy. In robotics, such symbolic models are often embedded into toolboxes like sympy.physics.mechanics or roboticstoolbox-python to derive small-signal models of manipulators and mobile robots. Below we compute the matrices \( \mathbf{A} \) and \( \mathbf{B} \) for the pendulum example.

import sympy as sp

# Symbols
theta, omega, u = sp.symbols('theta omega u')
m, L, b, g = sp.symbols('m L b g', positive=True)

# State and dynamics
x = sp.Matrix([theta, omega])

f1 = omega
f2 = -(b/(m*L**2))*omega - (g/L)*sp.sin(theta) + u/(m*L**2)
f = sp.Matrix([f1, f2])

# Jacobians: A = df/dx, B = df/du
A = f.jacobian(x)
B = f.jacobian(sp.Matrix([u]))

print("A(theta,omega,u) =")
sp.pprint(A)
print("\nB(theta,omega,u) =")
sp.pprint(B)

# Operating point (pendulum hanging down with zero torque)
theta_op = 0
omega_op = 0
u_op = 0

subs_dict = {theta: theta_op, omega: omega_op, u: u_op}
A_op = sp.simplify(A.subs(subs_dict))
B_op = sp.simplify(B.subs(subs_dict))

print("\nLinearized A matrix at (0,0,0) =")
sp.pprint(A_op)
print("\nLinearized B matrix at (0,0,0) =")
sp.pprint(B_op)

# Example numerical parameters (typical small pendulum)
A_num = A_op.subs({m: 1.0, L: 1.0, b: 0.1, g: 9.81})
B_num = B_op.subs({m: 1.0, L: 1.0})

print("\nA_num =")
sp.pprint(A_num)
print("\nB_num =")
sp.pprint(B_num)
      

In more complex robotic systems, symbolic differentiation may be replaced or complemented by automatic differentiation libraries and specialized dynamics engines (e.g., for articulated robots), but the mathematical principle is always the same: compute Jacobians of the vector field at a nominal operating point.

In C++-based robotic software stacks (e.g. ROS, simulation engines, whole-body controllers), the Eigen library is a standard choice for matrix operations. Below we implement the analytic linearization of the pendulum and compute the matrices \( \mathbf{A} \) and \( \mathbf{B} \) at a given operating point.

#include <iostream>
#include <Eigen/Dense>

using Eigen::Matrix2d;
using Eigen::Vector2d;

struct PendulumParams {
    double m;
    double L;
    double b;
    double g;
};

void linearizePendulum(const PendulumParams& p,
                       double theta_op,
                       double omega_op,
                       double u_op,
                       Matrix2d& A,
                       Vector2d& B)
{
    // For the simple pendulum, the Jacobians are analytically known.
    (void)u_op; // unused: equilibrium torque is zero here

    A.setZero();
    A(0,0) = 0.0;
    A(0,1) = 1.0;
    A(1,0) = -p.g / p.L * std::cos(theta_op);
    A(1,1) = -p.b / (p.m * p.L * p.L);

    B(0) = 0.0;
    B(1) = 1.0 / (p.m * p.L * p.L);
}

int main() {
    PendulumParams p{1.0, 1.0, 0.1, 9.81};

    Matrix2d A;
    Vector2d B;

    double theta_op = 0.0;
    double omega_op = 0.0;
    double u_op     = 0.0;

    linearizePendulum(p, theta_op, omega_op, u_op, A, B);

    std::cout << "A =\n" << A << "\n\n";
    std::cout << "B =\n" << B << "\n";

    return 0;
}
      

In larger robotic systems, similar functions are used to linearize around nominal configurations or trajectories. For example, model-based controllers may repeatedly call such routines online to update local LTI approximations.

Java is used in several robotics-related environments (for example competition robots or Java-based control frameworks). Libraries such as EJML or Apache Commons Math provide matrix utilities, but for a 2-dimensional example we can use plain arrays:

public class PendulumLinearization {

    public static double[][] computeA(double m, double L, double b, double g,
                                      double thetaOp, double omegaOp) {

        // omegaOp is not needed explicitly here but included for completeness
        double[][] A = new double[2][2];
        A[0][0] = 0.0;
        A[0][1] = 1.0;

        A[1][0] = -g / L * Math.cos(thetaOp);
        A[1][1] = -b / (m * L * L);
        return A;
    }

    public static double[] computeB(double m, double L) {
        double[] B = new double[2];
        B[0] = 0.0;
        B[1] = 1.0 / (m * L * L);
        return B;
    }

    public static void main(String[] args) {
        double m = 1.0, L = 1.0, b = 0.1, g = 9.81;
        double thetaOp = 0.0, omegaOp = 0.0;

        double[][] A = computeA(m, L, b, g, thetaOp, omegaOp);
        double[] B = computeB(m, L);

        System.out.println("A matrix:");
        for (int i = 0; i < 2; ++i) {
            System.out.println(A[i][0] + "  " + A[i][1]);
        }

        System.out.println("\nB vector:");
        System.out.println(B[0] + "  " + B[1]);
    }
}
      

In practical Java-based robotics projects, you would typically wrap such linearization routines inside higher-level classes that interact with sensor data and actuator commands, and store matrices in dedicated linear algebra types.

MATLAB is a standard tool in control engineering. For symbolic linearization we can use the Symbolic Math Toolbox; for Simulink-based models, the linearize or linmod functions compute linear models around operating points, heavily used in robotics and mechatronics workflows.

10.1 Symbolic linearization in MATLAB

syms theta omega u m L b g real

x = [theta; omega];

f1 = omega;
f2 = -(b/(m*L^2))*omega - (g/L)*sin(theta) + u/(m*L^2);
f  = [f1; f2];

A = jacobian(f, x);
B = jacobian(f, u);

% Operating point (pendulum hanging down)
theta_op = 0;
omega_op = 0;
u_op     = 0;

A_op = subs(A, {theta, omega, u}, {theta_op, omega_op, u_op});
B_op = subs(B, {theta, omega, u}, {theta_op, omega_op, u_op});

% Example numerical parameters
A_num = double(subs(A_op, {m, L, b, g}, {1, 1, 0.1, 9.81}))
B_num = double(subs(B_op, {m, L}, {1, 1}))
      

10.2 Simulink linearization around an operating point

Suppose you have a Simulink model pendulum_model with input block u and output theta. You can obtain an LTI state-space approximation via:

% Define linear analysis points
io(1) = linio('pendulum_model/u',     1, 'input');
io(2) = linio('pendulum_model/theta', 1, 'output');

% Compute an operating point (e.g. steady-state)
op = operpoint('pendulum_model');

% Linearize about operating point
sys_lin = linearize('pendulum_model', op, io);

A = sys_lin.A;
B = sys_lin.B;
C = sys_lin.C;
D = sys_lin.D;
      

This workflow is standard in industrial control and robotics: nonlinear Simulink models of actuators and mechanisms are linearized automatically for use in controller design and analysis.

Wolfram Mathematica has strong symbolic capabilities useful for automatic linearization of nonlinear models. With functions like D and Grad, Jacobians can be computed directly.

Clear[theta, omega, u, m, L, b, g];

x = {theta, omega};

f = {
   omega,
   -(b/(m L^2)) omega - (g/L) Sin[theta] + u/(m L^2)
   };

(* Jacobians *)
A = D[f, {x}];
B = D[f, {u}];

(* Operating point: hanging down *)
thetaOp = 0;
omegaOp = 0;
uOp     = 0;

Aop = A /. {theta -> thetaOp, omega -> omegaOp, u -> uOp};
Bop = B /. {theta -> thetaOp, omega -> omegaOp, u -> uOp};

MatrixForm[Aop]
MatrixForm[Bop]
      

In more advanced projects, Mathematica can be used to symbolically linearize high-order robotic manipulators and then export the resulting matrices to C++, Python, or MATLAB for real-time implementation.

Problem 1 (Scalar Nonlinear System): Consider the scalar nonlinear system \( \dot{x} = -x + x^3 + u \).
(a) Show that \( (x^\ast,u^\ast) = (0,0) \) is an equilibrium.
(b) Derive the linearized model around \( (0,0) \) in terms of \( \delta x \) and \( \delta u \).

Solution:

(a) At an equilibrium, \( 0 = f(x^\ast,u^\ast) = -x^\ast + (x^\ast)^3 + u^\ast \). Substituting \( x^\ast = 0 \), \( u^\ast = 0 \) gives \( 0 = 0 + 0 + 0 \), so \( (0,0) \) is indeed an equilibrium.

(b) We compute the partial derivatives:

\[ \frac{\partial f}{\partial x} = -1 + 3x^2, \qquad \frac{\partial f}{\partial u} = 1. \]

Evaluating at \( (x^\ast,u^\ast) = (0,0) \):

\[ a = \frac{\partial f}{\partial x}\bigg|_{(0,0)} = -1, \qquad b = \frac{\partial f}{\partial u}\bigg|_{(0,0)} = 1. \]

Therefore the linearized model in perturbation coordinates is

\[ \delta \dot{x}(t) \approx -\,\delta x(t) + \delta u(t). \]

Problem 2 (Mass-Spring-Damper with Cubic Stiffness): A 1-DOF mechanical system satisfies \( m\ddot{y} + c\dot{y} + ky + \alpha y^3 = u \).
(a) Choose state variables and write the system as \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},u) \).
(b) Find the equilibrium corresponding to \( y^\ast = 0 \) and \( u^\ast = 0 \).
(c) Derive the linearized state matrix \( \mathbf{A} \) and input matrix \( \mathbf{B} \) at that equilibrium.

Solution:

(a) Define \( x_1 = y \), \( x_2 = \dot{y} \). Then

\[ \begin{aligned} \dot{x}_1 &= x_2, \\ \dot{x}_2 &= -\frac{c}{m} x_2 - \frac{k}{m} x_1 - \frac{\alpha}{m} x_1^3 + \frac{1}{m} u. \end{aligned} \]

So the vector field is

\[ \mathbf{f}(\mathbf{x},u) = \begin{bmatrix} x_2 \\ -\dfrac{c}{m} x_2 - \dfrac{k}{m} x_1 - \dfrac{\alpha}{m} x_1^3 + \dfrac{1}{m} u \end{bmatrix}. \]

(b) With \( y^\ast = 0 \), \( \dot{y}^\ast = 0 \), and \( u^\ast = 0 \), we have \( x_1^\ast = 0 \), \( x_2^\ast = 0 \). Substituting into the dynamics:

\[ m\ddot{y}^\ast + c\dot{y}^\ast + ky^\ast + \alpha (y^\ast)^3 = 0 + 0 + 0 + 0 = 0 = u^\ast, \]

so \( (\mathbf{x}^\ast,u^\ast) = (\mathbf{0},0) \) is an equilibrium.

(c) Compute the Jacobian matrices. First,

\[ \mathbf{A}(\mathbf{x},u) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}} = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2} \\[4pt] \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\dfrac{k}{m} - \dfrac{3\alpha}{m} x_1^2 & -\dfrac{c}{m} \end{bmatrix}. \]

At the equilibrium \( x_1^\ast = 0 \),

\[ \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -\dfrac{k}{m} & -\dfrac{c}{m} \end{bmatrix}. \]

For \( \mathbf{B} \) we have

\[ \mathbf{B}(\mathbf{x},u) = \frac{\partial \mathbf{f}}{\partial u} = \begin{bmatrix} \dfrac{\partial f_1}{\partial u} \\[4pt] \dfrac{\partial f_2}{\partial u} \end{bmatrix} = \begin{bmatrix} 0 \\[4pt] \dfrac{1}{m} \end{bmatrix}, \]

which is independent of the state and remains the same at the equilibrium.

Problem 3 (Pendulum Linearized Around the Upright Position): Using the pendulum dynamics from Section 5, linearize the system around the upright equilibrium \( \theta^\ast = \pi \), \( \dot{\theta}^\ast = 0 \). Find the corresponding \( \mathbf{A} \) and \( \mathbf{B} \) matrices.

Solution:

We use the same state definition: \( x_1 = \theta \), \( x_2 = \dot{\theta} \). The vector field was

\[ \mathbf{f}(\mathbf{x},u) = \begin{bmatrix} x_2 \\ -\dfrac{b}{mL^2} x_2 - \dfrac{g}{L}\sin(x_1) + \dfrac{1}{mL^2} u \end{bmatrix}. \]

The Jacobians are unchanged:

\[ \mathbf{A}(\mathbf{x},u) = \begin{bmatrix} 0 & 1 \\ -\dfrac{g}{L}\cos(x_1) & -\dfrac{b}{mL^2} \end{bmatrix}, \qquad \mathbf{B}(\mathbf{x},u) = \begin{bmatrix} 0 \\[4pt] \dfrac{1}{mL^2} \end{bmatrix}. \]

At the upright position, \( x_1^\ast = \pi \), \( x_2^\ast = 0 \). We use \( \cos(\pi) = -1 \), so

\[ \mathbf{A} = \begin{bmatrix} 0 & 1 \\ \dfrac{g}{L} & -\dfrac{b}{mL^2} \end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} 0 \\[4pt] \dfrac{1}{mL^2} \end{bmatrix}. \]

Note that the sign of the \( g/L \) term changed relative to the downward equilibrium, reflecting the fact that the upright position is unstable: small perturbations grow rather than decay.

Problem 4 (Removal of Constant Terms in Perturbation Coordinates): Let \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},\mathbf{u}) \) and let \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) be an equilibrium. Show that when you express the first-order Taylor expansion in terms of \( \delta \mathbf{x} = \mathbf{x} - \mathbf{x}^\ast \), \( \delta \mathbf{u} = \mathbf{u} - \mathbf{u}^\ast \), the constant term disappears automatically.

Solution:

The Taylor expansion of \( \mathbf{f} \) around \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) is

\[ \begin{aligned} \mathbf{f}(\mathbf{x},\mathbf{u}) &\approx \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \mathbf{A}\,(\mathbf{x}-\mathbf{x}^\ast) + \mathbf{B}\,(\mathbf{u}-\mathbf{u}^\ast) \\ &= \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \mathbf{A}\,\delta \mathbf{x} + \mathbf{B}\,\delta \mathbf{u}. \end{aligned} \]

Because \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) is an equilibrium, we have \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) = \mathbf{0} \). Therefore,

\[ \mathbf{f}(\mathbf{x},\mathbf{u}) \approx \mathbf{A}\,\delta \mathbf{x} + \mathbf{B}\,\delta \mathbf{u}, \]

and hence the linearized perturbation dynamics are

\[ \delta \dot{\mathbf{x}}(t) = \mathbf{A}\,\delta \mathbf{x}(t) + \mathbf{B}\,\delta \mathbf{u}(t), \]

with no additional constant term. This is exactly the standard LTI form used throughout linear control.

• Real-world control systems, especially in robotics and mechatronics, are often nonlinear, but many analysis tools require linear time-invariant models.
• Linearization around an equilibrium uses first-order Taylor expansions of the nonlinear vector field \( \mathbf{f}(\mathbf{x},\mathbf{u}) \) and output map \( \mathbf{g}(\mathbf{x},\mathbf{u}) \).
• The Jacobian matrices \( \mathbf{A} = \partial \mathbf{f}/\partial \mathbf{x} \), \( \mathbf{B} = \partial \mathbf{f}/\partial \mathbf{u} \), \( \mathbf{C} = \partial \mathbf{g}/\partial \mathbf{x} \), \( \mathbf{D} = \partial \mathbf{g}/\partial \mathbf{u} \) define an LTI perturbation model \( \delta \dot{\mathbf{x}} = \mathbf{A}\delta \mathbf{x} + \mathbf{B}\delta \mathbf{u} \), \( \delta \mathbf{y} = \mathbf{C}\delta \mathbf{x} + \mathbf{D}\delta \mathbf{u} \).
• For mechanical examples like pendulums and mass-spring-damper systems, linearization recovers familiar second-order linear models when operating near small perturbations.
• Symbolic and numeric tools in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica allow control engineers to automate the process of linearization, which is foundational for the linear control techniques developed in the next chapters.