Consider a general nonlinear state-space model (continuous time): \( \mathbf{x}(t)\in\mathbb{R}^n \) (state), \( \mathbf{u}(t)\in\mathbb{R}^m \) (input), \( \mathbf{y}(t)\in\mathbb{R}^p \) (output),

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

A central theme in control engineering is that many designs and analyses are performed “locally,” i.e., for trajectories near a constant operating condition. This requires three precise notions:

• \( \mathbf{x}^\star \) is an equilibrium (for a specified constant input) if it makes the state stop changing.
• \( (\mathbf{x}^\star,\mathbf{u}^\star,\mathbf{y}^\star) \) is an operating point if it is a steady state together with consistent constant input/output values.
• \( \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}^\star \), \( \tilde{\mathbf{u}}=\mathbf{u}-\mathbf{u}^\star \), \( \tilde{\mathbf{y}}=\mathbf{y}-\mathbf{y}^\star \) are perturbation variables that shift the operating point to the origin.

flowchart TD
  A["Given nonlinear model: xdot=f(x,u), y=g(x,u)"] --> B["Choose constant input u* (bias/setpoint)"]
  B --> C["Solve equilibrium condition: f(x*, u*) = 0"]
  C --> D["Compute output at operating point: y* = g(x*, u*)"]
  D --> E["Define perturbations: xtilde=x-x*, utilde=u-u*, ytilde=y-y*"]
  E --> F["Rewrite dynamics in perturbation coordinates"]
  F --> G["Local analysis/design is performed around origin in (xtilde, utilde)"]
        

We distinguish continuous-time (CT) and discrete-time (DT) notions. For CT, suppose \( \mathbf{u}(t)=\mathbf{u}^\star \) is constant.

Definition (CT equilibrium). A state \( \mathbf{x}^\star \) is an equilibrium under input \( \mathbf{u}^\star \) if

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

Proposition 1 (equilibrium generates a constant trajectory). If \( \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0} \) and \( \mathbf{u}(t)=\mathbf{u}^\star \), then the trajectory \( \mathbf{x}(t)\equiv \mathbf{x}^\star \) satisfies the dynamics for all \( t \ge 0 \).

Proof. Substitute \( \mathbf{x}(t)=\mathbf{x}^\star \). Then \( \dot{\mathbf{x}}(t)=\mathbf{0} \), while the model gives \( \dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0} \). Hence the differential equation is satisfied identically. □

For DT systems, we assume the model is given by an update map \( \mathbf{x}_{k+1}=\mathbf{F}(\mathbf{x}_k,\mathbf{u}_k) \) and output \( \mathbf{y}_k=\mathbf{G}(\mathbf{x}_k,\mathbf{u}_k) \).

Definition (DT equilibrium / fixed point). Under constant input \( \mathbf{u}_k=\mathbf{u}^\star \), a state \( \mathbf{x}^\star \) is an equilibrium if

\[ \mathbf{x}^\star=\mathbf{F}(\mathbf{x}^\star,\mathbf{u}^\star). \]

Proposition 2 (DT equilibrium is a fixed trajectory). If \( \mathbf{x}^\star=\mathbf{F}(\mathbf{x}^\star,\mathbf{u}^\star) \) and \( \mathbf{u}_k=\mathbf{u}^\star \), then \( \mathbf{x}_k\equiv\mathbf{x}^\star \) for all \( k\ge 0 \) satisfies the update equation.

Proof. If \( \mathbf{x}_k=\mathbf{x}^\star \), then \( \mathbf{x}_{k+1}=\mathbf{F}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{x}^\star \). By induction, the constant sequence is a solution. □

These definitions are purely structural: they do not yet address whether nearby trajectories approach or depart from the equilibrium. That stability question will be formalized later (Chapter 11), but our present aim is to set the coordinate framework needed for local approximations.

In engineering practice, a “steady operating condition” includes not just a constant state, but also the constant inputs (biases, constant loads, setpoints) and the corresponding steady outputs.

Definition (operating point, CT). A triple \( (\mathbf{x}^\star,\mathbf{u}^\star,\mathbf{y}^\star) \) is an operating point if

\[ \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0}, \qquad \mathbf{y}^\star=\mathbf{g}(\mathbf{x}^\star,\mathbf{u}^\star). \]

Definition (operating point, DT). A triple \( (\mathbf{x}^\star,\mathbf{u}^\star,\mathbf{y}^\star) \) is an operating point if

\[ \mathbf{x}^\star=\mathbf{F}(\mathbf{x}^\star,\mathbf{u}^\star), \qquad \mathbf{y}^\star=\mathbf{G}(\mathbf{x}^\star,\mathbf{u}^\star). \]

Note the asymmetry: equilibrium is a state condition (under an assumed constant input), whereas an operating point bundles state, input, and output values. This bundling is essential for “small-signal” modeling because deviations are defined relative to the chosen operating point.

Let \( (\mathbf{x}^\star,\mathbf{u}^\star,\mathbf{y}^\star) \) be a CT operating point. Define perturbation variables: \( \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}^\star \), \( \tilde{\mathbf{u}}=\mathbf{u}-\mathbf{u}^\star \), \( \tilde{\mathbf{y}}=\mathbf{y}-\mathbf{y}^\star \).

Proposition 3 (exact perturbation dynamics). The perturbation state satisfies

\[ \dot{\tilde{\mathbf{x}}}(t) = \mathbf{f}\!\big(\mathbf{x}^\star+\tilde{\mathbf{x}}(t),\mathbf{u}^\star+\tilde{\mathbf{u}}(t)\big) - \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star). \]

Since \( \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0} \) at an operating point, this simplifies to

\[ \dot{\tilde{\mathbf{x}}}(t) = \mathbf{f}\!\big(\mathbf{x}^\star+\tilde{\mathbf{x}}(t),\mathbf{u}^\star+\tilde{\mathbf{u}}(t)\big). \]

Proof. By definition, \( \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}^\star \). Differentiate both sides: \( \dot{\tilde{\mathbf{x}}}=\dot{\mathbf{x}}-\mathbf{0}=\dot{\mathbf{x}} \) because \( \mathbf{x}^\star \) is constant. Substitute the original dynamics \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \), then replace \( \mathbf{x}=\mathbf{x}^\star+\tilde{\mathbf{x}} \) and \( \mathbf{u}=\mathbf{u}^\star+\tilde{\mathbf{u}} \). □

Proposition 4 (origin is an equilibrium in perturbation coordinates). If \( \tilde{\mathbf{u}}(t)\equiv\mathbf{0} \), then \( \tilde{\mathbf{x}}(t)\equiv\mathbf{0} \) is a solution of the perturbation system.

Proof. Substitute \( \tilde{\mathbf{x}}=\mathbf{0} \), \( \tilde{\mathbf{u}}=\mathbf{0} \): \( \dot{\tilde{\mathbf{x}}}=\mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0} \). Hence the origin is an equilibrium. □

The output perturbation satisfies the exact relation \( \tilde{\mathbf{y}}=\mathbf{g}(\mathbf{x}^\star+\tilde{\mathbf{x}},\mathbf{u}^\star+\tilde{\mathbf{u}})-\mathbf{g}(\mathbf{x}^\star,\mathbf{u}^\star) \).

flowchart TD
  X["Physical coordinates (x,u,y)"] --> S["Shift by operating point (x*,u*,y*)"]
  S --> T["Perturbation coordinates (xtilde, utilde, ytilde)"]
  T --> O["Origin corresponds to operating point"]
  O --> L["Local behavior studied for small |xtilde| and small |utilde|"]
        

This coordinate shift is exact and introduces no approximation. The approximation step (linearization) will be introduced next lesson by replacing the nonlinear dependence on \( (\tilde{\mathbf{x}},\tilde{\mathbf{u}}) \) with a first-order model.

Consider a simple pendulum (angle \( \theta \), angular rate \( \omega \)) with viscous damping and an applied torque input \( u \). A common nonlinear model is:

\[ \dot{\theta}=\omega,\qquad \dot{\omega}= -\frac{g}{\ell}\sin(\theta) - \frac{b}{m\ell^2}\,\omega + \frac{1}{m\ell^2}u. \]

Let \( u(t)=u^\star \) be constant. An equilibrium must satisfy \( \dot{\theta}=0 \) and \( \dot{\omega}=0 \), hence \( \omega^\star=0 \) and

\[ -\frac{g}{\ell}\sin(\theta^\star) + \frac{1}{m\ell^2}u^\star = 0 \quad\Longleftrightarrow\quad \sin(\theta^\star)=\frac{u^\star}{mg\ell}. \]

Therefore equilibria exist only if \( |u^\star| \le mg\ell \). In typical regulation around the downward position, one selects the equilibrium near \( \theta^\star \approx 0 \) when \( u^\star \approx 0 \). The corresponding steady output (if \( y=\theta \)) is \( y^\star=\theta^\star \).

Define perturbations: \( \tilde{\theta}=\theta-\theta^\star \), \( \tilde{\omega}=\omega-\omega^\star=\omega \), \( \tilde{u}=u-u^\star \). The exact perturbation dynamics become

\[ \dot{\tilde{\theta}}=\tilde{\omega},\qquad \dot{\tilde{\omega}}= -\frac{g}{\ell}\sin(\theta^\star+\tilde{\theta}) -\frac{b}{m\ell^2}\tilde{\omega} +\frac{1}{m\ell^2}(u^\star+\tilde{u}). \]

Using the equilibrium condition \( -\frac{g}{\ell}\sin(\theta^\star)+\frac{1}{m\ell^2}u^\star=0 \), you can equivalently write

\[ \dot{\tilde{\omega}}= -\frac{g}{\ell}\Big(\sin(\theta^\star+\tilde{\theta})-\sin(\theta^\star)\Big) -\frac{b}{m\ell^2}\tilde{\omega} +\frac{1}{m\ell^2}\tilde{u}. \]

This form isolates how deviations drive deviations; it is the standard starting point for local approximation in the next lesson.

For a CT model with constant input \( \mathbf{u}^\star \), computing an equilibrium is a root-finding problem:

\[ \text{Find }\mathbf{x}^\star\in\mathbb{R}^n \text{ such that } \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0}. \]

For DT models, it is a fixed-point problem:

\[ \text{Find }\mathbf{x}^\star \text{ such that } \mathbf{x}^\star-\mathbf{F}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0}. \]

A widely used numerical method is Newton’s method. In the scalar case, to solve \( \phi(z)=0 \), Newton updates

\[ z_{k+1}=z_k-\frac{\phi(z_k)}{\phi'(z_k)}. \]

In the vector case, Newton’s method solves linear systems involving the derivative matrix of \( \mathbf{f} \) with respect to \( \mathbf{x} \). The full Jacobian construction is deferred to Lesson 2; here, we focus on the operating-point framing and provide both library-based and from-scratch implementations.

We use the pendulum equilibrium equation (near the downward equilibrium) with parameters \( m=1 \), \( \ell=1 \), \( g=9.81 \), \( b=0.2 \), and a small bias torque \( u^\star=1.0 \). We then simulate a small perturbation response with \( \tilde{\theta}(0) \), \( \tilde{\omega}(0) \) small and \( \tilde{u}(t) \) as a small step.

7.1 Python (NumPy/SciPy)

import numpy as np
from math import sin, cos
from scipy.optimize import root
from scipy.integrate import solve_ivp

# Parameters
m, ell, g, b = 1.0, 1.0, 9.81, 0.2
u_star = 1.0  # bias torque

# --- Equilibrium solve: sin(theta*) = u*/(m g ell)
def phi(theta):
    return (g/ell) * np.sin(theta) - (1.0/(m*ell**2)) * u_star

def dphi(theta):
    return (g/ell) * np.cos(theta)

# Newton (from scratch) near theta=0
theta = 0.0
for _ in range(10):
    theta = theta - phi(theta)/dphi(theta)
theta_star = theta
omega_star = 0.0
print("theta_star (Newton):", theta_star)

# Library root finder (as a cross-check)
sol = root(lambda th: (g/ell)*np.sin(th) - u_star/(m*ell**2), x0=np.array([0.0]))
print("theta_star (root):", sol.x[0])

# --- Perturbation simulation
def pendulum_state(t, x, u_func):
    theta, omega = x
    u = u_func(t)
    return [
        omega,
        -(g/ell)*np.sin(theta) - (b/(m*ell**2))*omega + (1.0/(m*ell**2))*u
    ]

# small step perturbation on input: u(t) = u* + utilde(t)
def u_of_t(t):
    utilde = 0.05 if t >= 0.5 else 0.0
    return u_star + utilde

# initial condition = operating point + small perturbation
x0 = np.array([theta_star + 0.02, omega_star + 0.0])

t_span = (0.0, 5.0)
ts = np.linspace(t_span[0], t_span[1], 1000)

sol_ivp = solve_ivp(lambda t, x: pendulum_state(t, x, u_of_t),
                    t_span, x0, t_eval=ts, rtol=1e-8, atol=1e-10)

theta = sol_ivp.y[0]
omega = sol_ivp.y[1]

theta_tilde = theta - theta_star
omega_tilde = omega - omega_star
print("max |theta_tilde|:", np.max(np.abs(theta_tilde)))
      

7.2 C++ (from-scratch Newton + RK4; optional Eigen)

#include <iostream>
#include <cmath>
#include <vector>

struct Params {
  double m, ell, g, b, u_star;
};

double phi(double theta, const Params& p) {
  return (p.g/p.ell)*std::sin(theta) - (1.0/(p.m*p.ell*p.ell))*p.u_star;
}
double dphi(double theta, const Params& p) {
  return (p.g/p.ell)*std::cos(theta);
}

// Newton (scalar) to find theta_star near 0
double newton_theta_star(const Params& p, double theta0, int iters=12) {
  double th = theta0;
  for (int k=0; k<iters; ++k) {
    th = th - phi(th, p)/dphi(th, p);
  }
  return th;
}

// Pendulum dynamics
void f(double t, const std::vector<double>& x, std::vector<double>& xdot, const Params& p) {
  double theta = x[0];
  double omega = x[1];
  double utilde = (t >= 0.5) ? 0.05 : 0.0;
  double u = p.u_star + utilde;

  xdot[0] = omega;
  xdot[1] = -(p.g/p.ell)*std::sin(theta) - (p.b/(p.m*p.ell*p.ell))*omega + (1.0/(p.m*p.ell*p.ell))*u;
}

// RK4 step
std::vector<double> rk4_step(double t, const std::vector<double>& x, double h, const Params& p) {
  std::vector<double> k1(2), k2(2), k3(2), k4(2), xtmp(2), out(2);
  f(t, x, k1, p);

  xtmp[0] = x[0] + 0.5*h*k1[0]; xtmp[1] = x[1] + 0.5*h*k1[1];
  f(t + 0.5*h, xtmp, k2, p);

  xtmp[0] = x[0] + 0.5*h*k2[0]; xtmp[1] = x[1] + 0.5*h*k2[1];
  f(t + 0.5*h, xtmp, k3, p);

  xtmp[0] = x[0] + h*k3[0]; xtmp[1] = x[1] + h*k3[1];
  f(t + h, xtmp, k4, p);

  out[0] = x[0] + (h/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
  out[1] = x[1] + (h/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);
  return out;
}

int main() {
  Params p{1.0, 1.0, 9.81, 0.2, 1.0};
  double theta_star = newton_theta_star(p, 0.0);
  double omega_star = 0.0;

  std::cout << "theta_star=" << theta_star << std::endl;

  // simulate with small perturbation
  double t0=0.0, tf=5.0, h=0.001;
  std::vector<double> x{theta_star + 0.02, omega_star};

  double max_abs_theta_tilde = 0.0;
  for (double t=t0; t<=tf; t+=h) {
    double theta_tilde = x[0] - theta_star;
    if (std::fabs(theta_tilde) > max_abs_theta_tilde) max_abs_theta_tilde = std::fabs(theta_tilde);
    x = rk4_step(t, x, h, p);
  }
  std::cout << "max |theta_tilde|=" << max_abs_theta_tilde << std::endl;
  return 0;
}
      

Practical library note: for larger systems, C++ workflows often use Eigen for vectors/matrices and Boost.Odeint for ODE integration. The equilibrium computation generalizes to multidimensional root finding (e.g., damped Newton) once derivatives are available.

7.3 Java (Apache Commons Math or from-scratch)

import java.util.function.DoubleUnaryOperator;

public class PendulumOperatingPoint {

  static class Params {
    double m=1.0, ell=1.0, g=9.81, b=0.2, uStar=1.0;
  }

  // Scalar Newton method
  static double newton(DoubleUnaryOperator phi, DoubleUnaryOperator dphi, double x0, int iters) {
    double x = x0;
    for (int k=0; k<iters; k++) {
      x = x - phi.applyAsDouble(x) / dphi.applyAsDouble(x);
    }
    return x;
  }

  // RK4 integrator for 2-state pendulum
  static double[] rk4Step(double t, double[] x, double h, Params p) {
    double[] k1 = f(t, x, p);
    double[] x2 = new double[]{ x[0] + 0.5*h*k1[0], x[1] + 0.5*h*k1[1] };
    double[] k2 = f(t + 0.5*h, x2, p);
    double[] x3 = new double[]{ x[0] + 0.5*h*k2[0], x[1] + 0.5*h*k2[1] };
    double[] k3 = f(t + 0.5*h, x3, p);
    double[] x4 = new double[]{ x[0] + h*k3[0], x[1] + h*k3[1] };
    double[] k4 = f(t + h, x4, p);

    return new double[]{
      x[0] + (h/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]),
      x[1] + (h/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
    };
  }

  static double[] f(double t, double[] x, Params p) {
    double theta = x[0];
    double omega = x[1];
    double uTilde = (t >= 0.5) ? 0.05 : 0.0;
    double u = p.uStar + uTilde;

    double dtheta = omega;
    double domega = -(p.g/p.ell)*Math.sin(theta) - (p.b/(p.m*p.ell*p.ell))*omega + (1.0/(p.m*p.ell*p.ell))*u;
    return new double[]{dtheta, domega};
  }

  public static void main(String[] args) {
    Params p = new Params();

    DoubleUnaryOperator phi = (th) -> (p.g/p.ell)*Math.sin(th) - (1.0/(p.m*p.ell*p.ell))*p.uStar;
    DoubleUnaryOperator dphi = (th) -> (p.g/p.ell)*Math.cos(th);

    double thetaStar = newton(phi, dphi, 0.0, 12);
    double omegaStar = 0.0;

    System.out.println("thetaStar=" + thetaStar);

    double[] x = new double[]{thetaStar + 0.02, omegaStar};
    double t0=0.0, tf=5.0, h=1e-3;
    double maxAbsThetaTilde = 0.0;

    for (double t=t0; t<=tf; t+=h) {
      double thetaTilde = x[0] - thetaStar;
      maxAbsThetaTilde = Math.max(maxAbsThetaTilde, Math.abs(thetaTilde));
      x = rk4Step(t, x, h, p);
    }
    System.out.println("max |thetaTilde|=" + maxAbsThetaTilde);
  }
}
      

Practical library note: Apache Commons Math provides ODE integrators and root solvers for more advanced workflows, but the above illustrates the operating-point and perturbation concepts with minimal dependencies.

7.4 MATLAB (ODE + equilibrium) and Simulink operating point workflow

% Parameters
m = 1; ell = 1; g = 9.81; b = 0.2;
u_star = 1.0;

% Equilibrium equation: (g/ell)*sin(theta_star) - u_star/(m*ell^2) = 0
phi  = @(th) (g/ell)*sin(th) - u_star/(m*ell^2);
dphi = @(th) (g/ell)*cos(th);

% Newton (from scratch)
th = 0.0;
for k = 1:12
  th = th - phi(th)/dphi(th);
end
theta_star = th;
omega_star = 0.0;

disp(theta_star);

% Simulate full nonlinear system with a small input perturbation
u_of_t = @(t) u_star + (t >= 0.5)*0.05;

f = @(t,x) [ x(2);
            -(g/ell)*sin(x(1)) - (b/(m*ell^2))*x(2) + (1/(m*ell^2))*u_of_t(t) ];

x0 = [theta_star + 0.02; omega_star];
tspan = [0 5];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,x] = ode45(f, tspan, x0, opts);

theta_tilde = x(:,1) - theta_star;
omega_tilde = x(:,2) - omega_star;
disp(max(abs(theta_tilde)));
      

Simulink note (operating point concept). If you model the same dynamics in Simulink, the operating point corresponds to constant states and inputs. Typical workflows include:

• Set constant torque source to \( u^\star \) and solve for steady-state initial conditions.
• If available, use operating-point tools (e.g., findop, trim) to compute a consistent steady state for nonlinear models.
• Once \( (\mathbf{x}^\star,\mathbf{u}^\star) \) is computed, define perturbations by subtracting bias values in signals (or by using IC blocks with offsets).

7.5 Wolfram Mathematica (FindRoot + NDSolve)

(* Parameters *)
m = 1; ell = 1; g = 9.81; b = 0.2;
uStar = 1.0;

(* Equilibrium solve near 0 *)
phi[th_] := (g/ell) Sin[th] - uStar/(m ell^2);
thetaStar = th /. FindRoot[phi[th] == 0, {th, 0.0}];
omegaStar = 0.0;

thetaStar

(* Input with small step perturbation *)
u[t_] := uStar + If[t >= 0.5, 0.05, 0.0];

(* Nonlinear simulation *)
sol = NDSolve[
  {
    theta'[t] == omega[t],
    omega'[t] == -(g/ell) Sin[theta[t]] - (b/(m ell^2)) omega[t] + (1/(m ell^2)) u[t],
    theta[0] == thetaStar + 0.02,
    omega[0] == omegaStar
  },
  {theta, omega},
  {t, 0, 5}
];

(* Perturbations *)
thetaTilde[t_] := (theta[t] /. sol[[1]]) - thetaStar;
maxAbsThetaTilde = NMaxValue[{Abs[thetaTilde[t]], 0 <= t <= 5}, t];
maxAbsThetaTilde
      

Problem 1 (CT equilibrium verification): Consider \( \dot{x}=f(x,u) \) with constant input \( u(t)=u^\star \). Show that if \( f(x^\star,u^\star)=0 \), then \( x(t)\equiv x^\star \) is a solution.

Solution: Since \( x(t)=x^\star \) is constant, \( \dot{x}(t)=0 \). The dynamics require \( \dot{x}(t)=f(x^\star,u^\star)=0 \). Hence the differential equation holds for all \( t \ge 0 \). □

Problem 2 (DT fixed-point characterization): For \( x_{k+1}=F(x_k,u_k) \) with \( u_k=u^\star \), prove that \( x^\star=F(x^\star,u^\star) \) implies \( x_k\equiv x^\star \).

Solution: If \( x_k=x^\star \), then \( x_{k+1}=F(x^\star,u^\star)=x^\star \). By induction from \( k=0 \), the constant sequence satisfies the update. □

Problem 3 (perturbation dynamics derivation): Let \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \). Given an operating point \( (\mathbf{x}^\star,\mathbf{u}^\star) \) with \( \mathbf{f}(\mathbf{x}^\star,\mathbf{u}^\star)=\mathbf{0} \), define \( \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}^\star \), \( \tilde{\mathbf{u}}=\mathbf{u}-\mathbf{u}^\star \). Derive the exact equation for \( \dot{\tilde{\mathbf{x}}} \).

Solution: Differentiate \( \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}^\star \) to get \( \dot{\tilde{\mathbf{x}}}=\dot{\mathbf{x}} \). Substitute \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \), then replace \( \mathbf{x}=\mathbf{x}^\star+\tilde{\mathbf{x}} \), \( \mathbf{u}=\mathbf{u}^\star+\tilde{\mathbf{u}} \): \( \dot{\tilde{\mathbf{x}}}=\mathbf{f}(\mathbf{x}^\star+\tilde{\mathbf{x}},\mathbf{u}^\star+\tilde{\mathbf{u}}) \). □

Problem 4 (pendulum operating point existence): For the pendulum model in Section 5, show that equilibria exist only if \( |u^\star| \le mg\ell \).

Solution: At equilibrium, \( \omega^\star=0 \) and \( \sin(\theta^\star)=u^\star/(mg\ell) \). Since \( \sin(\cdot) \) satisfies \( -1 \le \sin(\cdot) \le 1 \), we must have \( -1 \le u^\star/(mg\ell) \le 1 \), i.e. \( |u^\star| \le mg\ell \). □

Problem 5 (numerical equilibrium computation): Let \( \phi(\theta)=(g/\ell)\sin(\theta)-u^\star/(m\ell^2) \). Write one Newton step from \( \theta_k \) to \( \theta_{k+1} \) and specify \( \phi'(\theta) \).

Solution: We have \( \phi'(\theta)=(g/\ell)\cos(\theta) \). Newton’s update is \( \theta_{k+1}=\theta_k-\phi(\theta_k)/\phi'(\theta_k) \), i.e. \( \theta_{k+1}=\theta_k-\dfrac{(g/\ell)\sin(\theta_k)-u^\star/(m\ell^2)}{(g/\ell)\cos(\theta_k)} \). □

We defined equilibria (CT and DT) and operating points as consistent steady state tuples \( (\mathbf{x}^\star,\mathbf{u}^\star,\mathbf{y}^\star) \). We then introduced perturbation variables \( (\tilde{\mathbf{x}},\tilde{\mathbf{u}},\tilde{\mathbf{y}}) \) and proved that the coordinate shift is exact, with the operating point mapped to the origin. Finally, we framed operating point computation as a root/fixed-point problem and implemented equilibrium computation plus small-perturbation simulations across Python, C++, Java, MATLAB/Simulink, and Mathematica. These constructions are the necessary setup for Jacobian linearization in Lesson 2.