At the end of Lesson 1, we obtained discrete-time models for linearized robot dynamics by sampling a continuous-time model with sampling period \( T_s \). For a linearized joint or task-space control loop we write

\[ x_{k+1} = A_d x_k + B_d u_k,\quad y_k = C_d x_k + D_d u_k \]

where \( x_k \) is the sampled state at time \( t_k = k T_s \), \( u_k \) is the control input held constant between \( t_k \) and \( t_{k+1} \) by a zero-order hold, and \( y_k \) is the sampled output (e.g., joint position, end-effector position).

In a typical robot joint servo, we use a discrete-time controller

\[ u_k = K_d x_k + K_r r_k \]

where \( K_d \) is a discrete-time feedback gain (e.g., joint-space PD in state form) and \( r_k \) is a reference signal. The autonomous closed-loop dynamics (for constant reference) have the form

\[ x_{k+1} = A_{\text{cl}} x_k,\quad A_{\text{cl}} = A_d + B_d K_d. \]

Discrete-time stability and performance of the robot controller are therefore determined by the properties of \( A_{\text{cl}} \).

flowchart TD
  CT["Linearized robot dynamics (continuous-time)"] --> DISCR["Sample + hold with Ts"]
  DISCR --> DYN["Discrete model (Ad, Bd, Cd, Dd)"]
  DYN --> CTRL["Design discrete controller Kd"]
  CTRL --> CL["Closed-loop Ad_cl = Ad + Bd Kd"]
  CL --> STAB["Check eigenvalues |lambda| < 1"]
  CL --> PERF["Evaluate performance \n(settling, overshoot, bandwidth)"]
        

Consider first the scalar system \( x_{k+1} = a x_k \). By recursion,

\[ x_k = a^k x_0. \]

Then \( x_k \to 0 \) as \( k \to \infty \) if and only if \( |a| < 1 \). If \( |a| > 1 \), the trajectories diverge; if \( |a| = 1 \), the system is marginally stable at best (oscillatory for \( a = -1 \)).

For the multi-dimensional LTI system \( x_{k+1} = A_d x_k \), we similarly have \( x_k = A_d^k x_0 \). The spectral radius \( \rho(A_d) \) is defined as

\[ \rho(A_d) = \max_i |\lambda_i(A_d)| \]

where \( \lambda_i(A_d) \) are the eigenvalues of \( A_d \).

Theorem (Schur stability). The discrete-time system \( x_{k+1} = A_d x_k \) is asymptotically stable if and only if

\[ \rho(A_d) < 1, \]

i.e., all eigenvalues of \( A_d \) lie strictly inside the unit circle in the complex plane.

For robot joint control, we often start from a continuous-time closed-loop matrix \( A_c \) and obtain \( A_d = e^{A_c T_s} \) by exact discretization. Let \( \lambda_c \) be an eigenvalue of \( A_c \). Its discrete-time counterpart is

\[ z = e^{\lambda_c T_s}. \]

If \( \Re(\lambda_c) < 0 \), then \( |z| = e^{\Re(\lambda_c) T_s} < 1 \), so continuous-time asymptotic stability implies discrete-time asymptotic stability under exact discretization. However, aggressive sampling or approximate discretization schemes can distort performance, even if stability is preserved.

Many robot controllers are designed in state space; Lyapunov methods provide constructive tools to certify stability, and they extend naturally to the discrete-time case.

Consider \( x_{k+1} = A_d x_k \). A quadratic Lyapunov function candidate is

\[ V(x_k) = x_k^\top P x_k,\quad P \succ 0, \]

where \( P \) is a symmetric positive definite matrix. The Lyapunov difference is

\[ \Delta V = V(x_{k+1}) - V(x_k) = x_k^\top \big(A_d^\top P A_d - P\big) x_k. \]

If we can find \( P \succ 0 \) and \( Q \succ 0 \) such that

\[ A_d^\top P A_d - P = -Q, \]

then \( \Delta V = - x_k^\top Q x_k \leq - \lambda_{\min}(Q) \|x_k\|^2 \), which proves global asymptotic stability of the origin. This is the discrete-time Lyapunov equation. Equivalently, the Lyapunov inequality

\[ A_d^\top P A_d - P \prec 0 \]

also certifies asymptotic stability. For joint-space controllers expressed as \( A_{\text{cl}} = A_d + B_d K_d \), we simply replace \( A_d \) by \( A_{\text{cl}} \) in the Lyapunov equation.

In practice, \( P \) is found numerically by solving a discrete Lyapunov equation (e.g., in MATLAB or Mathematica). This is widely used in robot control design to verify that the sampled closed loop is stable for a given sampling period \( T_s \) and gain matrix \( K_d \).

For robot control, stability is necessary but not sufficient: we also require fast, well-damped tracking with bounded overshoot and limited control effort. For linear time-invariant discrete-time systems with a real or complex pole structure, time-domain performance is tightly linked to the location of poles in the \( z \)-plane.

Consider a second-order closed-loop system with complex conjugate poles \( z_{1,2} = r e^{\pm j\theta} \), where \( 0 < r < 1 \) and \( 0 < \theta < \pi \). Suppose these arise from sampling a continuous-time second-order system with natural frequency \( \omega_n \) and damping ratio \( \zeta \). Let \( \omega_d = \omega_n \sqrt{1 - \zeta^2} \) be the damped natural frequency. Then the mapping

\[ z = e^{(-\zeta \omega_n \pm j \omega_d) T_s} \quad \Rightarrow \quad r = e^{-\zeta \omega_n T_s},\;\; \theta = \omega_d T_s \]

connects continuous and discrete-time specifications. Roughly:

• Poles closer to the origin (\( r \) small) correspond to faster decay (shorter settling time).
• Poles closer to the unit circle (\( r \approx 1 \)) produce slow decay and longer settling time.
• The angle \( \theta \) controls the oscillation frequency in samples.

For a second-order underdamped system, the continuous-time step response overshoot \( M_p \) is approximately

\[ M_p \approx \exp\!\left( -\frac{\pi \zeta}{\sqrt{1 - \zeta^2}} \right). \]

The widely used 2 % settling time approximation is

\[ T_s \approx \frac{4}{\zeta \omega_n}. \]

Given \( T_s \) and \( T_s \approx N_s T_s \) in samples, we get an approximate number of samples to settle:

\[ N_s \approx \frac{4}{\zeta \omega_n T_s}. \]

For a first-order discrete-time system with pole \( z \), \( 0 < z < 1 \), the step response behaves like \( 1 - z^{k+1} \). A useful rule of thumb is that the response is essentially settled after about \( k \approx \frac{4}{1 - z} \) samples.

flowchart TD
  ZP["Closed-loop poles (z-plane)"] --> R["Radius r = |z| (speed of decay)"]
  ZP --> TH["Angle theta (oscillation per sample)"]
  R --> TS["Samples to settle Ns ~ 4 / ((1 - r) or zeta*wn*Ts)"]
  TH --> OV["Overshoot Mp ~ exp(-pi*zeta/sqrt(1 - zeta^2))"]
  TS --> SPEC["Check Ts, bandwidth vs robot task"]
  OV --> SPEC
        

Consider a single revolute joint with linearized dynamics around an operating point

\[ J \ddot q(t) + b \dot q(t) = \tau(t), \]

where \( J > 0 \) is the equivalent inertia, \( b \geq 0 \) is viscous friction, and \( \tau(t) \) is the actuator torque. Let \( q_{\text{ref}} \) be a constant joint reference. Define the tracking error and its derivative

\[ e(t) = q_{\text{ref}} - q(t), \quad \dot e(t) = - \dot q(t). \]

A continuous-time PD torque is

\[ \tau(t) = K_p e(t) + K_d \dot e(t). \]

Set state variables \( x_1(t) = e(t) \), \( x_2(t) = \dot e(t) \). Using \( \ddot e(t) = - \ddot q(t) \) and the dynamics, we obtain

\[ \ddot e(t) = -\frac{1}{J}\big( K_p e(t) + K_d \dot e(t) + b \dot e(t) \big). \]

Therefore

\[ \dot x_1(t) = x_2(t), \quad \dot x_2(t) = -\frac{K_p}{J} x_1(t) -\frac{K_d + b}{J} x_2(t). \]

The continuous-time closed-loop matrix is

\[ A_c = \begin{bmatrix} 0 & 1 \\ -K_p/J & -(K_d + b)/J \end{bmatrix}. \]

For a sampling period \( T_s \), the exact discrete-time closed-loop matrix (under zero-order hold and constant reference) is

\[ A_d = e^{A_c T_s}. \]

The eigenvalues of \( A_d \) are \( z_i = e^{\lambda_i(A_c) T_s} \). Hence:

• Choose \( K_p, K_d \) so that the continuous-time eigenvalues \( \lambda_i(A_c) \) have negative real parts satisfying the desired \( \zeta, \omega_n \).
• Select \( T_s \) such that \( \omega_n T_s \ll \pi \), typically \( T_s \approx \frac{1}{20} \) of the dominant time constant or faster, to preserve performance.
• Verify discrete-time stability by checking \( |z_i| < 1 \) and performance by checking \( r, \theta \) as in Section 4.

This workflow is standard in joint-space PD control implemented on a digital controller: tuning is done in continuous-time, but the final certification is discrete-time.

In practice, discrete-time stability and performance analysis for robot control is implemented via numerical libraries. Below we assume given \( A_d \) and illustrate how to compute eigenvalues and settling times in several languages. We also mention relevant robotics-related libraries in each ecosystem that interact with the control layer.

6.1 Python (NumPy, SciPy, robotics-toolbox, python-control)

In Python, robot control is often prototyped with numpy, scipy, python-control, and higher-level packages such as roboticstoolbox-python.

import numpy as np
from numpy.linalg import eig
from math import log, sqrt

# Discrete-time closed-loop matrix Ad (example for one joint)
Ad = np.array([[0.95, 0.01],
               [-0.2,  0.90]])

Ts = 0.005  # sampling period [s]

# Eigenvalues (discrete poles)
z, _ = eig(Ad)

print("Discrete poles:", z)

# Check stability: all |z| < 1
stable = np.all(np.abs(z) < 1.0)
print("Schur-stable?", stable)

# Approximate dominant pole and settling time
r = max(np.abs(z))
if r < 1.0:
    # approximate decay time constant tau by matching r = exp(-Ts/tau)
    tau = -Ts / log(r)
    Ts_settle = 4.0 * tau   # 2 percent criterion
    print("Approx continuous settling time [s]:", Ts_settle)
    print("Approx samples to settle:", Ts_settle / Ts)

# Using python-control to construct a discrete system
import control as ct

sysd = ct.ss(Ad, np.zeros((2,1)), np.eye(2), np.zeros((2,1)), Ts)
print("System is stable?", ct.isdtime(sysd) and ct.pole(sysd))
      

6.2 C++ (Eigen, ROS control, Orocos KDL)

In C++, real-time robot controllers typically run within ROS-based or custom real-time frameworks. Libraries such as Eigen (linear algebra), ros_control, and Orocos KDL handle kinematics and dynamics, while the discrete-time control loop is implemented explicitly.

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

int main() {
  Eigen::Matrix2d Ad;
  Ad << 0.95, 0.01,
         -0.2, 0.90;

  Eigen::EigenSolver<Eigen::Matrix2d> es(Ad);
  Eigen::VectorXcd z = es.eigenvalues();

  std::cout << "Discrete poles: " << z.transpose() << std::endl;

  bool stable = true;
  for (int i = 0; i < z.size(); ++i) {
    if (std::abs(z[i]) >= 1.0) {
      stable = false;
    }
  }
  std::cout << "Schur-stable? " << (stable ? "yes" : "no") << std::endl;

  return 0;
}

// In a ROS control loop, Ad would come from offline design, and the loop would
// update x_{k+1} = Ad * x_k at each timer callback with a fixed period Ts.
      

6.3 Java (Apache Commons Math, Java-based robotics frameworks)

Java-based robotic platforms (e.g., some industrial controllers, educational robots, FRC's WPILib) can use Apache Commons Math for linear algebra and eigen computations. The control loop runs in a periodic task.

import org.apache.commons.math3.linear.*;

public class DiscreteStabilityDemo {
    public static void main(String[] args) {
        double[][] adData = {
            {0.95, 0.01},
            {-0.2, 0.90}
        };
        RealMatrix Ad = MatrixUtils.createRealMatrix(adData);

        EigenDecomposition ed = new EigenDecomposition(Ad);
        double[] real = ed.getRealEigenvalues();
        double[] imag = ed.getImagEigenvalues();

        boolean stable = true;
        for (int i = 0; i < real.length; ++i) {
            double mag = Math.hypot(real[i], imag[i]);
            System.out.println("Pole " + i + ": " + real[i] + " + j" + imag[i]
                               + " |z| = " + mag);
            if (mag >= 1.0) stable = false;
        }
        System.out.println("Schur-stable? " + stable);
    }
}
      

6.4 MATLAB/Simulink (Robotics System Toolbox, Control System Toolbox)

In MATLAB, the combination of Robotics System Toolbox and Control System Toolbox is standard. One designs a controller in continuous-time, discretizes it with c2d, and analyzes discrete-time poles.

Ad = [0.95  0.01;
      -0.2  0.90];
Ts = 0.005;

sysd = ss(Ad, [], eye(2), [], Ts);
z = eig(Ad);

disp('Discrete poles:');
disp(z);

if all(abs(z) < 1)
    disp('System is Schur-stable.');
else
    disp('System is NOT Schur-stable.');
end

% Example: from continuous-time PD design for a joint
J = 0.5; b = 0.1;
Kp = 25; Kd = 4;

Ac = [0  1;
      -Kp/J  -(Kd + b)/J];

Tc = 0.005;
sysc = ss(Ac, [], eye(2), []);
sysd_pd = c2d(sysc, Tc, 'zoh');
eig(sysd_pd)   % discrete-time poles of PD controlled joint

% In Simulink, one uses Discrete-Time blocks with Ts and scopes to
% visualize settling time and overshoot of joint position.
      

6.5 Wolfram Mathematica (StateSpaceModel, Discrete-time Analysis)

Mathematica supports symbolic and numeric analysis of discrete-time state-space models.

Ad = { {0.95, 0.01},
      {-0.2,  0.90} };
Ts = 0.005;

sysd = StateSpaceModel[{Ad, { {0}, {0} }, IdentityMatrix[2], { {0}, {0} } },
                      SamplingPeriod -> Ts];

(* Eigenvalues (poles) *)
z = Eigenvalues[Ad];
Print["Discrete poles: ", z];

stableQ = Max[Abs[z]] < 1.0;
Print["Schur-stable? ", stableQ];

(* Approximate settling time from dominant pole *)
r = Max[Abs[z]];
tau = -Ts/Log[r];
TsSettle = 4.0 * tau;
Print["Approx. settling time [s]: ", TsSettle];
      

Problem 1 (Scalar Discrete Stability and Settling Samples). A scalar discrete-time error dynamics for a robot joint is \( e_{k+1} = 0.92\, e_k \) with sampling period \( T_s = 1 \,\text{ms} \).
(a) Is the system asymptotically stable?
(b) Approximate how many samples are required for the error magnitude to decrease to at most 2 % of its initial value.

Solution.
(a) The pole is \( z = 0.92 \). Since \( |z| = 0.92 < 1 \), the system is asymptotically stable.
(b) The error is \( e_k = 0.92^k e_0 \). We require \( |e_k| \leq 0.02 |e_0| \), i.e.

\[ 0.92^k \leq 0.02 \quad \Rightarrow \quad k \geq \frac{\ln(0.02)}{\ln(0.92)} \approx \frac{-3.912}{-0.0834} \approx 46.9. \]

So about 47 samples are needed. With \( T_s = 1 \,\text{ms} \), the settling time is roughly \( 47 \,\text{ms} \).

Problem 2 (Schur Stability of a 2x2 Closed Loop). Consider the discrete-time closed-loop matrix of a robot joint controller

\[ A_d = \begin{bmatrix} 0.8 & 0.01 \\ -0.3 & 0.85 \end{bmatrix}. \]

(a) Derive the characteristic polynomial of \( A_d \).
(b) Determine whether the system is asymptotically stable by checking the modulus of the eigenvalues (you may work symbolically or numerically).

Solution.
(a) The characteristic polynomial is

\[ \det(z I - A_d) = \det \begin{bmatrix} z - 0.8 & -0.01 \\ 0.3 & z - 0.85 \end{bmatrix} = (z - 0.8)(z - 0.85) + 0.003. \]

Expanding,

\[ p(z) = z^2 - 1.65 z + (0.68 + 0.003) = z^2 - 1.65 z + 0.683. \]

(b) The eigenvalues are the roots of \( z^2 - 1.65 z + 0.683 = 0 \):

\[ z = \frac{1.65 \pm \sqrt{1.65^2 - 4 \cdot 0.683}}{2} = \frac{1.65 \pm \sqrt{2.7225 - 2.732}}{2} \approx \frac{1.65 \pm j 0.0987}{2}. \]

The magnitude of each eigenvalue is approximately \( |z| \approx \sqrt{(0.825)^2 + (0.04935)^2} \approx 0.826 \), which is less than 1. Hence the system is Schur-stable.

Problem 3 (Range of Gains for Scalar Discrete PD). A scalar discrete-time error dynamics for a linearized joint is modeled as

\[ e_{k+1} = e_k + T_s \big(-k_p e_k\big) = (1 - T_s k_p) e_k, \]

where \( T_s \) is fixed, and \( k_p > 0 \) is the proportional gain in a discrete controller.
(a) For a given \( T_s \), determine the range of \( k_p \) for which the system is asymptotically stable.
(b) For \( T_s = 0.01 \,\text{s} \), give the admissible range of \( k_p \).

Solution.
The discrete-time pole is \( z = 1 - T_s k_p \). Asymptotic stability requires \( |z| < 1 \), i.e.

\[ |1 - T_s k_p| < 1 \quad \Rightarrow \quad -1 < 1 - T_s k_p < 1. \]

This yields

\[ -1 < 1 - T_s k_p \quad \Rightarrow \quad T_s k_p < 2 \quad \Rightarrow \quad k_p < \frac{2}{T_s}, \]

\[ 1 - T_s k_p < 1 \quad \Rightarrow \quad 0 < T_s k_p \quad \Rightarrow \quad k_p > 0. \]

Therefore the range is \( 0 < k_p < \frac{2}{T_s} \). For \( T_s = 0.01 \), we get \( 0 < k_p < 200 \).

Problem 4 (Discrete Lyapunov Equation). Let

\[ A_d = \begin{bmatrix} 0.9 & 0 \\ 0 & 0.7 \end{bmatrix}, \quad Q = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. \]

(a) Solve the discrete Lyapunov equation \( A_d^\top P A_d - P = -Q \) for diagonal \( P = \mathrm{diag}(p_1, p_2) \).
(b) Show that the resulting \( P \) is positive definite and conclude that the system is asymptotically stable.

Solution.
(a) With \( P = \mathrm{diag}(p_1, p_2) \), we have

\[ A_d^\top P A_d = \begin{bmatrix} 0.9 & 0 \\ 0 & 0.7 \end{bmatrix}^\top \begin{bmatrix} p_1 & 0 \\ 0 & p_2 \end{bmatrix} \begin{bmatrix} 0.9 & 0 \\ 0 & 0.7 \end{bmatrix} = \begin{bmatrix} 0.9^2 p_1 & 0 \\ 0 & 0.7^2 p_2 \end{bmatrix}. \]

The Lyapunov equation becomes

\[ \begin{bmatrix} 0.9^2 p_1 - p_1 & 0 \\ 0 & 0.7^2 p_2 - p_2 \end{bmatrix} = -\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. \]

Hence

\[ (0.9^2 - 1) p_1 = -1 \quad \Rightarrow \quad p_1 = \frac{1}{1 - 0.9^2} = \frac{1}{1 - 0.81} = \frac{1}{0.19} > 0, \]

\[ (0.7^2 - 1) p_2 = -1 \quad \Rightarrow \quad p_2 = \frac{1}{1 - 0.7^2} = \frac{1}{1 - 0.49} = \frac{1}{0.51} > 0. \]

(b) Since \( p_1 > 0 \) and \( p_2 > 0 \), \( P \) is positive definite. The Lyapunov equation \( A_d^\top P A_d - P = -Q \) with \( Q \succ 0 \) implies asymptotic stability of \( x_{k+1} = A_d x_k \).

Problem 5 (Sampling Time and Performance Degradation). A continuous-time second-order closed loop for a robot joint has \( \zeta = 0.7 \) and \( \omega_n = 10 \,\text{rad/s} \). You implement the controller digitally with sampling period \( T_s \).
(a) Compute the desired continuous-time 2 % settling time \( T_s^{(c)} \).
(b) For \( T_s = 0.01 \,\text{s} \), approximate the number of samples required to settle.
(c) Qualitatively discuss what happens if \( T_s \) is increased to \( 0.05 \,\text{s} \).

Solution.
(a) The continuous-time approximation gives

\[ T_s^{(c)} \approx \frac{4}{\zeta \omega_n} = \frac{4}{0.7 \cdot 10} \approx 0.571 \,\text{s}. \]

(b) For \( T_s = 0.01 \,\text{s} \), the number of samples to settle is

\[ N_s \approx \frac{T_s^{(c)}}{T_s} \approx \frac{0.571}{0.01} \approx 57 \,\text{samples}. \]

(c) If \( T_s \) is increased to \( 0.05 \,\text{s} \), we would have only about 11 samples over the settling time. The discrete-time poles move closer to the unit circle, and the sampled response becomes coarser; performance degrades, potentially with increased overshoot and poor tracking. For too large \( T_s \), aliasing and discretization error may destabilize the closed-loop implementation, even if the continuous-time design is stable.

In this lesson we formalized stability and performance for discrete-time robot controllers. For linearized joint and task-space loops with sampled dynamics \( x_{k+1} = A_{\text{cl}} x_k \), Schur stability (\( \rho(A_{\text{cl}}) < 1 \)) guarantees asymptotic stability. Discrete-time Lyapunov equations provide a powerful matrix-based test and are widely used in robot control verification.

Performance metrics such as overshoot, settling time, and bandwidth were linked to the location of poles in the \( z \)-plane via the exponential mapping from continuous-time eigenvalues. For typical PD-controlled robot joints, tuning is often done in continuous-time and then translated to discrete-time, with sampling rate chosen to be sufficiently fast relative to the dominant dynamics.

Finally, we illustrated how to compute discrete-time poles and approximate settling times in Python, C++, Java, MATLAB/Simulink, and Mathematica, and connected these tools to robotics-specific libraries that integrate kinematics and dynamics with control loops. These foundations will be crucial when we consider multi-rate loops and real-time scheduling in the next lessons.