Classical linear control theory provides analytic tools (Laplace transforms, transfer functions, root locus, Bode plots, Nyquist criterion) for predicting and shaping the behavior of feedback systems. Software tools implement these constructions algorithmically. For a single-input single-output (SISO) LTI system with input \( u(t) \) and output \( y(t) \), the dynamics can be written as a linear differential equation

\[ a_n y^{(n)}(t) + a_{n-1} y^{(n-1)}(t) + \cdots + a_0 y(t) = b_m u^{(m)}(t) + b_{m-1} u^{(m-1)}(t) + \cdots + b_0 u(t), \]

where \( a_n \neq 0 \) and typically \( m \leq n \). Taking Laplace transforms under zero initial conditions yields

\[ (a_n s^n + \cdots + a_0)Y(s) = (b_m s^m + \cdots + b_0)U(s), \quad G(s) = \frac{Y(s)}{U(s)} = \frac{b_m s^m + \cdots + b_0}{a_n s^n + \cdots + a_0}. \]

Software tools store the numerator and denominator polynomials of \( G(s) \) as arrays of real numbers and then:

• compute time responses (step, impulse, arbitrary inputs) by solving the differential equation;
• compute frequency response by evaluating \( G(j\omega) \) on a frequency grid;
• analyze stability through pole locations, gain/phase margins, and root locus;
• assist controller design (PID tuning, compensator placement, loop shaping, etc.).

The typical workflow in computer-aided linear control design is summarized below.

flowchart TD
  A["Physical system & assumptions"] --> B["Derive ODE model"]
  B --> C["Form transfer function G(s)"]
  C --> D["Implement G(s) in tool (MATLAB / Python / C++ / Java / Mathematica)"]
  D --> E["Analyze: time response, frequency response, stability"]
  E --> F["Design controller (e.g. PID, lead-lag)"]
  F --> G["Simulate closed loop (possibly with robot model)"]
  G --> H["Validate and iterate design"]
        

Any proper rational SISO transfer function with real coefficients can be written as

\[ G(s) = \frac{b_m s^m + b_{m-1} s^{m-1} + \cdots + b_0} {a_n s^n + a_{n-1} s^{n-1} + \cdots + a_0}, \quad a_n \neq 0. \]

Software packages represent this object as two vectors of coefficients:

\[ \mathbf{b} = [b_m,\dots,b_0], \qquad \mathbf{a} = [a_n,\dots,a_0]. \]

For example, a standard second-order plant \( G(s) = \dfrac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2} \) corresponds to \( \mathbf{b} = [\omega_n^2] \) and \( \mathbf{a} = [1,\,2\zeta \omega_n,\,\omega_n^2] \).

In time-domain simulations, the software solves the ODE

\[ a_n y^{(n)}(t) + \cdots + a_0 y(t) = b_m u^{(m)}(t) + \cdots + b_0 u(t), \]

whereas in frequency-domain analysis it evaluates \( G(j\omega) \) on a grid \( \{\omega_k\}_{k=0}^{N-1} \), typically logarithmically spaced:

\[ \omega_k = \omega_{\min} \left( \frac{\omega_{\max}}{\omega_{\min}} \right)^{k/(N-1)}, \quad G(j\omega_k) = \frac{\sum_{i=0}^{m} b_i (j\omega_k)^i} {\sum_{i=0}^{n} a_i (j\omega_k)^i}. \]

The Bode magnitude and phase are then given by

\[ |G(j\omega_k)| = \sqrt{\Re\{G(j\omega_k)\}^2 + \Im\{G(j\omega_k)\}^2},\quad \angle G(j\omega_k) = \operatorname{atan2}(\Im\{G(j\omega_k)\},\Re\{G(j\omega_k)\}). \]

Thus, all major tools share the same mathematical core: polynomial evaluation in the complex plane and numerical solution of linear ODEs.

MATLAB with the Control System Toolbox is a de facto standard for classical linear control. Models are specified using tf (transfer functions) or related constructs. Simulink provides block-diagram simulation, which closely matches the block-diagram representations developed in earlier chapters.

Consider the benchmark second-order plant \( G(s) = \dfrac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} \). The following MATLAB script constructs the transfer function and computes its step and Bode responses.

% Second-order benchmark plant
zeta = 0.5;
wn   = 4.0;

num = [wn^2];              % numerator coefficients
den = [1  2*zeta*wn  wn^2];% denominator coefficients

G = tf(num, den);

% Time-domain analysis
figure;
step(G);
title('Step response of second-order plant');

% Frequency-domain analysis
figure;
bode(G);
grid on;
title('Bode plot of second-order plant');
      

For robotics applications, MATLAB offers Robotics System Toolbox and Robotics Toolbox (third-party), which allow:

• building rigid-body models of manipulators and mobile robots;
• linearizing nonlinear robot dynamics about an operating point to obtain an LTI model;
• connecting this LTI model to standard controllers (e.g., PID position loops) in Simulink.

A typical Simulink diagram for a robot joint position loop uses blocks corresponding to: the actuator dynamics \( G(s) \), a position sensor, a PID block, and summing junctions that represent the error computation \( e(t) = r(t) - y(t) \).

Python provides open-source alternatives to MATLAB via numpy, scipy, and the control (python-control) package, which mimics many MATLAB Control System Toolbox functions. In robotics, Python is widely used in ROS/ROS 2 nodes and for high-level planning using libraries such as roboticstoolbox and ROS control interfaces.

The code below recreates the same second-order plant and computes its step and Bode responses. We also show a basic “from scratch” implementation of the frequency response using polynomial evaluation.

import numpy as np
import matplotlib.pyplot as plt
import control as ctrl  # pip install control

# Second-order plant parameters
zeta = 0.5
wn   = 4.0

num = [wn**2]
den = [1.0, 2.0*zeta*wn, wn**2]

G = ctrl.TransferFunction(num, den)

# Time response: step
t, y = ctrl.step_response(G)

plt.figure()
plt.plot(t, y)
plt.xlabel("t [s]")
plt.ylabel("y(t)")
plt.title("Step response of second-order plant")
plt.grid(True)

# Frequency response using python-control
w = np.logspace(-1, 2, 400)
mag, phase, omega = ctrl.bode(G, w, Plot=False)

plt.figure()
plt.subplot(2, 1, 1)
plt.semilogx(omega, 20*np.log10(mag))
plt.ylabel("Magnitude [dB]")
plt.grid(True)

plt.subplot(2, 1, 2)
plt.semilogx(omega, phase * 180.0/np.pi)
plt.xlabel("Frequency [rad/s]")
plt.ylabel("Phase [deg]")
plt.grid(True)
plt.tight_layout()

# From-scratch frequency response evaluation
def freq_response(num, den, omega):
    """
    Evaluate G(j*omega) for a SISO transfer function defined by
    numerator and denominator coefficient lists.
    """
    s = 1j * omega
    num_val = np.polyval(num, s)
    den_val = np.polyval(den, s)
    return num_val / den_val

omega_test = np.logspace(-1, 2, 5)
G_vals = freq_response(num, den, omega_test)
print("Sampled G(j*omega):", G_vals)

plt.show()
      

In a robotics context, the same G could represent the linear approximation of a robot joint servo. The step response then corresponds to a joint angle command, and Bode analysis gives insight into bandwidth and noise sensitivity of that joint loop.

C++ is commonly used for real-time control on embedded hardware and in robotics middleware such as ROS/ROS 2. While C++ does not ship with a dedicated control toolbox, it provides:

• <complex> for complex arithmetic;
• <vector> for dynamic arrays;
• linear algebra libraries such as Eigen for matrix operations;
• robotics-oriented frameworks like ROS control and Orocos for integrating controllers with robots.

The following C++ code illustrates a small utility to evaluate \( G(j\omega) \) from numerator and denominator coefficients, suitable for computing points of a Bode plot. This is a basic “from scratch” implementation of the rational function evaluation.

#include <iostream>
#include <vector>
#include <complex>

std::complex<double> eval_poly(const std::vector<double> &coeffs,
                               std::complex<double> s) {
    // Horner's rule
    std::complex<double> acc(0.0, 0.0);
    for (double c : coeffs) {
        acc = acc * s + c;
    }
    return acc;
}

std::complex<double> freq_response(const std::vector<double> &num,
                                    const std::vector<double> &den,
                                    double omega) {
    std::complex<double> s(0.0, omega); // s = j*omega
    std::complex<double> num_val = eval_poly(num, s);
    std::complex<double> den_val = eval_poly(den, s);
    return num_val / den_val;
}

int main() {
    // Second-order plant G(s) = wn^2 / (s^2 + 2*zeta*wn*s + wn^2)
    double zeta = 0.5;
    double wn   = 4.0;

    std::vector<double> num{wn*wn};                         // [wn^2]
    std::vector<double> den{1.0, 2.0*zeta*wn, wn*wn};       // [1, 2*zeta*wn, wn^2]

    std::vector<double> omegas{0.1, 1.0, 10.0};
    for (double w : omegas) {
        std::complex<double> G = freq_response(num, den, w);
        double mag   = std::abs(G);
        double phase = std::arg(G);
        std::cout << "omega = " << w
                  << ", |G(j*omega)| = " << mag
                  << ", angle = " << phase << std::endl;
    }
    return 0;
}
      

In a robotic joint controller implemented in ROS 2, such C++ code can be integrated into a node that periodically computes control actions based on the linear model, while still relying on the analytic design performed offline in MATLAB or Python.

Java is less common in low-level control loops but is used in higher-level supervisory control and simulation frameworks (for example in some industrial SCADA systems or educational robotics platforms such as FIRST Robotics with WPILib). Numerical libraries such as Apache Commons Math provide ODE solvers and linear algebra.

Below we define a minimal Java class to simulate the step response of a second-order plant using an explicit Euler discretization. For the ODE

\[ \ddot{y}(t) + 2\zeta\omega_n \dot{y}(t) + \omega_n^2 y(t) = \omega_n^2 u(t), \]

introduce state variables \( x_1 = y \), \( x_2 = \dot{y} \). Then

\[ \dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = -2\zeta\omega_n x_2(t) - \omega_n^2 x_1(t) + \omega_n^2 u(t), \quad y(t) = x_1(t). \]

Using a fixed step \( h \), explicit Euler gives \( x_{k+1} = x_k + h f(x_k,u_k) \). The Java code is:

public class SecondOrderStepEuler {
    private final double zeta;
    private final double wn;

    public SecondOrderStepEuler(double zeta, double wn) {
        this.zeta = zeta;
        this.wn = wn;
    }

    public double[] simulateStep(double tFinal, double h) {
        int nSteps = (int) Math.round(tFinal / h);
        double[] y = new double[nSteps + 1];

        double x1 = 0.0; // position
        double x2 = 0.0; // velocity
        double u = 1.0;  // unit step

        for (int k = 0; k <= nSteps; ++k) {
            y[k] = x1;

            double dx1 = x2;
            double dx2 = -2.0 * zeta * wn * x2 - wn * wn * x1 + wn * wn * u;

            x1 += h * dx1;
            x2 += h * dx2;
        }
        return y;
    }

    public static void main(String[] args) {
        SecondOrderStepEuler sim = new SecondOrderStepEuler(0.5, 4.0);
        double[] y = sim.simulateStep(5.0, 0.001);
        System.out.println("Final value y(T) = " + y[y.length - 1]);
    }
}
      

In a Java-based robotics stack, such a class could be wrapped by higher-level components that compare measured joint positions to the simulated linear model, enabling model-based diagnostic or supervisory control.

Wolfram Mathematica offers powerful symbolic manipulation, which complements numerical tools by enabling closed-form derivations of transfer functions, partial fraction expansions, and analytic step responses. Mathematica also has functions for transfer function models and time/frequency response.

The following Mathematica code defines the second-order plant, computes its step response, and derives the analytic expression symbolically.

zeta = 0.5;
wn   = 4.0;

G = TransferFunctionModel[ wn^2 / (s^2 + 2 zeta wn s + wn^2), s ];

(* Numeric step response *)
stepPlot = OutputResponse[G, UnitStep[t], {t, 0, 5}];

(* Symbolic step response *)
symStep = InverseLaplaceTransform[
  (wn^2 / (s^2 + 2 zeta wn s + wn^2)) * (1/s),
  s, t
] // FullSimplify;

Print["Symbolic step response y(t) = ", symStep];

Plot[Evaluate[stepPlot], {t, 0, 5},
  AxesLabel -> {"t", "y(t)"},
  PlotLabel -> "Step response of second-order plant"
]
      

For robotics, Mathematica can model multibody dynamics and then linearize them around a nominal configuration to produce a transfer function \( G(s) \) that is exported to other environments (e.g., MATLAB or C++ code generation) for implementation.

In robotic systems, linear control tools rarely operate in isolation. Instead, there is a toolchain that connects:

• modeling and linearization (MATLAB/Simulink, Mathematica, or Python);
• controller design (loop shaping, PID tuning, compensators) in MATLAB or Python;
• deployment to embedded controllers (C++ in ROS control, Java in educational robotics, or auto-generated C from Simulink);
• online logging and analysis (Python notebooks, MATLAB scripts).

Robotics-specific libraries include, among others:

• MATLAB/Simulink: Robotics System Toolbox, Simscape Multibody, which provide linearization functions and ready-made joint controllers.
• Python: roboticstoolbox, ROS/ROS 2 Python APIs for control nodes, and simulation environments such as Gazebo.
• C++: ROS control, Orocos, and custom libraries structured around Eigen for real-time kinematics and dynamics.
• Java: WPILib and related libraries in educational robotics, providing PID controllers and model-based utilities.
• Mathematica: multibody and control functions for off-line symbolic computation and export.

All of these environments ultimately implement the same underlying mathematical operations on LTI models that you derive analytically in the rest of the course.

flowchart TD
  M["Derive LTI model G(s) from robot dynamics"] --> O["Offline design (MATLAB / Python / Mathematica)"]
  O --> C["Controller parameters (e.g. PID gains, compensator coefficients)"]
  C --> D["Deployment to C++ / Java / embedded code"]
  D --> R["Robot hardware / simulator"]
  R --> L["Logged data (trajectories, errors)"]
  L --> A["Analysis and redesign in high-level tools"]
        

Problem 1 (Polynomial Representation and Uniqueness): Let \( G(s) \) be a proper rational transfer function with real coefficients. Show that, up to a common nonzero scalar factor, the representation of \( G(s) \) by numerator and denominator coefficient vectors \( (\mathbf{b},\mathbf{a}) \) is unique.

Solution: Suppose we have two representations \( G(s) = \dfrac{B_1(s)}{A_1(s)} = \dfrac{B_2(s)}{A_2(s)} \) with real polynomials \( A_i(s), B_i(s) \) and \( A_i(0) \neq 0 \). Then

\[ B_1(s) A_2(s) = B_2(s) A_1(s). \]

Since the ring of real polynomials is an integral domain, if \( A_1(s) \) and \( A_2(s) \) are both chosen monic (leading coefficient equal to 1), then \( A_1(s) = A_2(s) \) and hence \( B_1(s) = B_2(s) \). Therefore the coefficient vectors \( \mathbf{a}, \mathbf{b} \) are unique up to a common nonzero scalar multiple. In software implementations it is standard to normalize by making the leading denominator coefficient \( a_n = 1 \), which fixes this ambiguity.

Problem 2 (Frequency Grid for Bode Plots): Consider the logarithmic frequency grid \( \omega_k = \omega_{\min} \left( \dfrac{\omega_{\max}}{\omega_{\min}} \right)^{k/(N-1)} \), \( k = 0,\dots,N-1 \). Show that the grid is geometrically spaced and compute the constant ratio \( r \) such that \( \omega_{k+1} = r \,\omega_k \).

Solution: By definition,

\[ \omega_k = \omega_{\min} \left( \frac{\omega_{\max}}{\omega_{\min}} \right)^{k/(N-1)}. \]

Therefore

\[ \frac{\omega_{k+1}}{\omega_k} = \frac{ \omega_{\min} \left( \dfrac{\omega_{\max}}{\omega_{\min}} \right)^{(k+1)/(N-1)} }{ \omega_{\min} \left( \dfrac{\omega_{\max}}{\omega_{\min}} \right)^{k/(N-1)} } = \left( \frac{\omega_{\max}}{\omega_{\min}} \right)^{1/(N-1)} =: r, \]

which is independent of \( k \). Hence the grid is geometric with constant ratio \( r = \left( \dfrac{\omega_{\max}}{\omega_{\min}} \right)^{1/(N-1)} \). This geometric spacing is what Bode plots implement in practice.

Problem 3 (Local Truncation Error of Explicit Euler): For the scalar ODE \( \dot{x}(t) = f(t,x(t)) \), explicit Euler with step size \( h \) is \( x_{k+1} = x_k + h f(t_k, x_k) \). Show that the local truncation error is \( \mathcal{O}(h^2) \).

Solution: Let \( x(t) \) be the exact solution. Expand \( x(t_k + h) \) in a Taylor series around \( t_k \):

\[ x(t_k + h) = x(t_k) + h \dot{x}(t_k) + \frac{h^2}{2} \ddot{x}(\xi_k), \]

for some \( \xi_k \) between \( t_k \) and \( t_k + h \) (Taylor's theorem). Using the ODE, \( \dot{x}(t_k) = f(t_k, x(t_k)) \), so

\[ x(t_k + h) = x(t_k) + h f(t_k, x(t_k)) + \frac{h^2}{2} \ddot{x}(\xi_k). \]

The Euler update with starting value \( x(t_k) \) gives \( \tilde{x}_{k+1} = x(t_k) + h f(t_k, x(t_k)) \). The local truncation error is

\[ \tau_{k+1} := x(t_k + h) - \tilde{x}_{k+1} = \frac{h^2}{2} \ddot{x}(\xi_k), \]

so \( |\tau_{k+1}| \leq C h^2 \) for some constant \( C \) if \( \ddot{x}(t) \) is bounded on the interval of interest. Hence the local truncation error is of order \( h^2 \). This estimate underlies the accuracy of the Java explicit Euler simulation shown earlier.

Problem 4 (Toolchain Choice for Robot Joint Control): You are tasked with designing a position controller for a single revolute robot joint, modeled as a second-order LTI system. Explain a reasonable division of labor between MATLAB/Simulink, Python, and C++ in a typical workflow.

Solution: A typical workflow is:

flowchart TD
  P["Plant identification or derivation of joint model"] --> M["Design in MATLAB/Simulink (PID / compensators)"]
  M --> PY["Rapid prototyping & data analysis in Python"]
  PY --> CXX["Real-time implementation in C++ on robot controller"]
  CXX --> LOG["Log data for validation"]
  LOG --> PY
        

Specifically:

• Use MATLAB/Simulink to derive \( G(s) \), tune PID gains, and verify time and frequency responses on accurate models.
• Use Python for flexible offline analysis, scripting parameter sweeps, and post-processing experimental data from the robot.
• Implement the final controller in C++ within ROS control or a similar real-time framework to meet timing guarantees and interface directly with robot sensors and actuators.

This division respects the strengths of each tool while preserving the core linear control mathematics learned in the course.

In this lesson we connected the mathematical representation of LTI systems via transfer functions to their implementations in major software ecosystems. The core idea is that all tools manipulate the same rational model \( G(s) \), represented through numerator and denominator polynomials, to compute time responses, frequency responses, and closed-loop behavior. MATLAB and Simulink provide a mature, integrated environment; Python offers open-source flexibility; C++ and Java enable deployment in real-time and large-scale software systems; Mathematica supports symbolic derivations that can be exported to other tools. Robotics-oriented libraries in each ecosystem embed these linear control capabilities into complete robot software stacks.

Subsequent lessons in this chapter will build on these tools to automate time-response and frequency-response analysis, root-locus and Bode design, and systematic robustness studies.