Consider an \( n \)-DOF robot manipulator with joint coordinates \( \mathbf{q}(t) \in \mathbb{R}^n \) and joint torques (actuator-space inputs) \( \boldsymbol{u}(t) \in \mathbb{R}^n \), often equal to motor torques \( \boldsymbol{u}(t) = \boldsymbol{\tau}(t) \). The standard rigid-body dynamics (which you already know from robotics kinematics/dynamics) can be written as

\[ \mathbf{M}(\mathbf{q}) \, \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\,\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) + \mathbf{f}(\dot{\mathbf{q}}) = \boldsymbol{\tau}_c(t) + \boldsymbol{\tau}_d(t) \]

where

• \( \mathbf{M}(\mathbf{q}) \) is the inertia matrix,
• \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) collects Coriolis and centrifugal terms,
• \( \mathbf{g}(\mathbf{q}) \) is gravity, \( \mathbf{f}(\dot{\mathbf{q}}) \) is friction, and
• \( \boldsymbol{\tau}_c \) is the control torque, \( \boldsymbol{\tau}_d \) an unknown disturbance.

A joint-space controller uses measurements of \( \mathbf{q}(t),\dot{\mathbf{q}}(t) \) and a joint-space reference \( \mathbf{q}_d(t),\dot{\mathbf{q}}_d(t),\ddot{\mathbf{q}}_d(t) \) to generate \( \boldsymbol{\tau}_c(t) \) so that the joints move as desired. The logical information flow is illustrated below.

flowchart TD
  Qd["Joint reference q_d(t)"] --> E["Error e(t) = q_d(t) - q(t)"]
  E --> Ctr["Joint-space controller"]
  Ctr --> Tau["Actuator torques tau(t)"]
  Tau --> Plant["Robot joints and dynamics"]
  Plant --> Meas["Measured q(t), qdot(t)"]
  Meas --> E
        

The central object in joint-space control is the error signal \( \mathbf{e}(t) \). All control objectives are expressed in terms of how this error behaves over time.

We assume the robot configuration is represented by joint angles or displacements \( \mathbf{q}(t) = [q_1(t),\dots,q_n(t)]^\top \). A convenient state vector is

\[ \mathbf{x}(t) = \begin{bmatrix} \mathbf{q}(t) \\ \dot{\mathbf{q}}(t) \end{bmatrix} \in \mathbb{R}^{2n}. \]

From the dynamics, we can write a first-order nonlinear state-space model with input \( \boldsymbol{u}(t) = \boldsymbol{\tau}_c(t) \):

\[ \dot{\mathbf{x}}(t) = \begin{bmatrix} \dot{\mathbf{q}}(t) \\ \mathbf{M}^{-1}(\mathbf{q}(t)) \big( \boldsymbol{\tau}_c(t) - \mathbf{C}(\mathbf{q}(t),\dot{\mathbf{q}}(t))\dot{\mathbf{q}}(t) - \mathbf{g}(\mathbf{q}(t)) - \mathbf{f}(\dot{\mathbf{q}}(t)) + \boldsymbol{\tau}_d(t) \big) \end{bmatrix} = \mathbf{f}\big(\mathbf{x}(t),\boldsymbol{\tau}_c(t),t\big). \]

While later chapters introduce sophisticated nonlinear controllers, in this lesson we only use this model to define what we want the closed-loop behavior to be. The system is also constrained by:

\[ q_i^{\min} \leq q_i(t) \leq q_i^{\max}, \quad |\dot{q}_i(t)| \leq \dot{q}_i^{\max}, \quad |\tau_{c,i}(t)| \leq \tau_i^{\max}, \]

which encode joint limits, velocity limits, and actuator torque limits. These constraints must be respected when specifying realistic joint-space objectives.

Let the desired joint configuration signal be \( \mathbf{q}_d(t) \). We define the tracking error as

\[ \mathbf{e}(t) = \mathbf{q}_d(t) - \mathbf{q}(t), \quad \dot{\mathbf{e}}(t) = \dot{\mathbf{q}}_d(t) - \dot{\mathbf{q}}(t). \]

Two standard classes of joint-space objectives are:

• Regulation: The desired configuration is constant, \( \mathbf{q}_d(t) \equiv \mathbf{q}^\star \) and the objective is

  \[ \lim_{t \to \infty} \mathbf{q}(t) = \mathbf{q}^\star, \quad \lim_{t \to \infty} \dot{\mathbf{q}}(t) = \mathbf{0}. \]

• Tracking: The desired configuration is time-varying, and the objective is

  \[ \lim_{t \to \infty} \|\mathbf{e}(t)\| = 0, \quad \lim_{t \to \infty} \|\dot{\mathbf{e}}(t)\| = 0. \]

In practice, exact convergence is impossible due to noise, unmodeled dynamics, and limited actuator bandwidth. Therefore we specify performance envelopes, for example:

\[ \|\mathbf{e}(t)\| \leq \varepsilon_{\infty} \quad \text{for all } t \geq T_s, \]

where \( T_s \) is a settling time and \( \varepsilon_{\infty} \) is the steady-state tracking accuracy.

A convenient way to describe joint-space objectives is to rewrite the dynamics in error coordinates. Define the error state vector

\[ \mathbf{x}_e(t) = \begin{bmatrix} \mathbf{e}(t) \\ \dot{\mathbf{e}}(t) \end{bmatrix} \in \mathbb{R}^{2n}. \]

Using the definitions of \( \mathbf{e}(t) \) and \( \dot{\mathbf{e}}(t) \), we can compute

\[ \ddot{\mathbf{e}}(t) = \ddot{\mathbf{q}}_d(t) - \ddot{\mathbf{q}}(t). \]

From the robot dynamics, assuming disturbances are negligible for now (\( \boldsymbol{\tau}_d(t) \approx \mathbf{0} \)), the joint accelerations satisfy

\[ \ddot{\mathbf{q}}(t) = \mathbf{M}^{-1}(\mathbf{q}(t)) \big( \boldsymbol{\tau}_c(t) - \mathbf{C}(\mathbf{q}(t),\dot{\mathbf{q}}(t))\dot{\mathbf{q}}(t) - \mathbf{g}(\mathbf{q}(t)) - \mathbf{f}(\dot{\mathbf{q}}(t)) \big). \]

Hence the error dynamics can always be written as

\[ \ddot{\mathbf{e}}(t) = \ddot{\mathbf{q}}_d(t) - \mathbf{M}^{-1}(\mathbf{q}(t)) \big( \boldsymbol{\tau}_c(t) - \mathbf{C}(\mathbf{q}(t),\dot{\mathbf{q}}(t))\dot{\mathbf{q}}(t) - \mathbf{g}(\mathbf{q}(t)) - \mathbf{f}(\dot{\mathbf{q}}(t)) \big), \]

which shows that any joint-space controller is essentially a mapping from \( (\mathbf{x}_e, \mathbf{q}_d,\dot{\mathbf{q}}_d,\ddot{\mathbf{q}}_d) \) to \( \boldsymbol{\tau}_c \) that shapes \( \mathbf{x}_e(t) \).

flowchart TD
  Meas["Measured q, qdot"] --> Err["Compute e = q_d - q, edot = qd_dot - qdot"]
  Ref["Reference q_d, qd_dot, qd_ddot"] --> Err
  Err --> Xe["Error state x_e = [e; edot]"]
  Xe --> Spec["Desired error dynamics + performance specs"]
        

In later lessons we will propose specific control laws whose closed-loop error system can be made to approximate a desired target form.

A common joint-space control objective is to make each joint error behave approximately like a decoupled linear second-order system:

\[ \ddot{e}_i(t) + 2 \zeta_i \omega_{n,i} \dot{e}_i(t) + \omega_{n,i}^2 e_i(t) = 0, \quad i = 1,\dots,n, \]

where \( \omega_{n,i} > 0 \) is the natural frequency and \( \zeta_i > 0 \) is the damping ratio for joint \( i \).

For a single joint (dropping index \( i \)), the associated characteristic polynomial is

\[ s^2 + 2 \zeta \omega_n s + \omega_n^2 = 0, \]

with roots

\[ s_{1,2} = -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1}. \]

For \( 0 < \zeta < 1 \) the system is underdamped, and the solution has oscillatory form

\[ e(t) = \mathrm{e}^{-\zeta \omega_n t} \big( A \cos(\omega_d t) + B \sin(\omega_d t) \big), \quad \omega_d = \omega_n \sqrt{1 - \zeta^2}, \]

for some constants \( A,B \) determined by initial conditions. From this standard form we obtain classical performance indices:

• Peak time:

  \[ T_p = \frac{\pi}{\omega_d}. \]

• Settling time (2 percent criterion):

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

• Maximum overshoot:

  \[ M_p = \exp\!\Bigg( - \frac{\pi \zeta}{\sqrt{1 - \zeta^2}} \Bigg). \]

Defining a joint-space objective often means: pick \( \zeta_i,\omega_{n,i} \) for each joint so that \( T_s,M_p \) satisfy the application requirements, then design the controller (next lesson) such that the closed-loop error dynamics approximate this target template.

Beyond peak time and overshoot, we can express joint-space control objectives as cost functionals. A common quadratic performance index on a finite horizon \( [0,T] \) is

\[ J_T = \int_0^T \Big( \mathbf{e}^\top(t) \mathbf{Q} \mathbf{e}(t) + \dot{\mathbf{e}}^\top(t) \mathbf{R} \dot{\mathbf{e}}(t) + \boldsymbol{\tau}_c^\top(t) \mathbf{S} \boldsymbol{\tau}_c(t) \Big) \, \mathrm{d}t, \]

where \( \mathbf{Q},\mathbf{R},\mathbf{S} \) are positive semidefinite weighting matrices. Special cases include:

• Integral of squared error (ISE):

  \[ \mathrm{ISE} = \int_0^\infty \mathbf{e}^\top(t)\mathbf{e}(t) \, \mathrm{d}t. \]

• Integral of absolute error (IAE):

  \[ \mathrm{IAE} = \int_0^\infty \|\mathbf{e}(t)\|_1 \, \mathrm{d}t. \]

In optimal control (later chapters) we explicitly minimize such indices, but even in purely feedback design they provide a precise language for what "good joint tracking" means.

Real controllers are digital. Let the sampling period be \( T_s \). The controller observes \( \mathbf{q}[k] \approx \mathbf{q}(kT_s) \) and \( \dot{\mathbf{q}}[k] \approx \dot{\mathbf{q}}(kT_s) \), and outputs a torque command \( \boldsymbol{\tau}_c[k] \) held approximately constant over \( [kT_s,(k+1)T_s) \).

The discrete-time tracking error is

\[ \mathbf{e}[k] = \mathbf{q}_d[k] - \mathbf{q}[k], \quad \dot{\mathbf{e}}[k] = \dot{\mathbf{q}}_d[k] - \dot{\mathbf{q}}[k]. \]

The joint-space objectives translate directly:

• Exponential-like decay of \( \|\mathbf{e}[k]\| \) at a specified rate.
• Bounds on discretized overshoot and settling time in samples.
• Respecting discrete-time torque limits \( |\tau_{c,i}[k]| \leq \tau_i^{\max} \).

Later digital-control chapters will refine these notions, but for now the key idea is that all objectives defined in continuous time can be discretized consistently once \( T_s \) is fixed.

In Python, joint-space variables are often represented using numpy arrays. Robotics libraries such as roboticstoolbox (Corke) or modern_robotics (Lynch & Park) provide convenient interfaces to actual robot models. Here we focus on computing errors and a simple integral-of-squared-error metric over a discrete trajectory.

import numpy as np

# Example: n-DOF joint trajectory, sampled at N steps
# q_d[k], q[k] are arrays in R^n

def joint_space_error(q_d, q, qd_d=None, qd=None, dt=0.001):
    """
    Compute joint-space position and velocity tracking error over a trajectory.

    Parameters
    ----------
    q_d : array_like, shape (N, n)
        Desired joint positions.
    q : array_like, shape (N, n)
        Measured joint positions.
    qd_d : array_like, shape (N, n) or None
        Desired joint velocities (optional).
    qd : array_like, shape (N, n) or None
        Measured joint velocities (optional).
    dt : float
        Sample time.

    Returns
    -------
    e, edot : np.ndarray
        Position and velocity error trajectories.
    ise : float
        Integral of squared position error.
    """
    q_d = np.asarray(q_d)
    q = np.asarray(q)
    assert q_d.shape == q.shape

    e = q_d - q

    if qd_d is not None and qd is not None:
        qd_d = np.asarray(qd_d)
        qd = np.asarray(qd)
        assert qd_d.shape == qd.shape
        edot = qd_d - qd
    else:
        # crude finite-difference approximation
        edot = np.zeros_like(e)
        edot[1:] = (e[1:] - e[:-1]) / dt

    # ISE over all joints and time
    ise = float(np.sum(np.sum(e**2, axis=1)) * dt)
    return e, edot, ise


# Minimal example: 2-DOF arm tracking a sinusoidal reference
if __name__ == "__main__":
    dt = 0.001
    T = 2.0
    t = np.arange(0.0, T, dt)
    N = t.size
    n = 2

    # desired joint trajectories
    qd1 = 0.5 * np.sin(2.0 * np.pi * t)
    qd2 = 0.3 * np.sin(4.0 * np.pi * t)
    q_d = np.stack([qd1, qd2], axis=1)

    # measured: assume a small phase lag and amplitude error
    q = np.stack([
        0.48 * np.sin(2.0 * np.pi * (t - 0.01)),
        0.29 * np.sin(4.0 * np.pi * (t - 0.01)),
    ], axis=1)

    e, edot, ise = joint_space_error(q_d, q, dt=dt)
    print("ISE over horizon T = ", T, " seconds:", ise)
      

This kind of function is a building block for evaluating whether a particular controller meets the joint-space objectives specified in terms of error norms and integral criteria.

9.1 C++ (Eigen + robotics middlewares)

In C++, joint-space computations are typically performed with Eigen for linear algebra, possibly integrated with middleware such as ROS Control or Orocos. The snippet below assumes Eigen::VectorXd for joint vectors.

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

struct JointTrajectorySample {
    Eigen::VectorXd q_d;   // desired positions
    Eigen::VectorXd q;     // measured positions
};

double computeISE(const std::vector<JointTrajectorySample>& traj,
                  double dt)
{
    if (traj.empty()) return 0.0;
    const int n = traj.front().q_d.size();
    double ise = 0.0;

    for (const auto& s : traj) {
        Eigen::VectorXd e = s.q_d - s.q;
        ise += e.squaredNorm() * dt;
    }
    return ise;
}

int main() {
    const int n = 2;
    double dt = 0.001;
    std::vector<JointTrajectorySample> traj;

    // Fill with some fake data (for illustration)
    for (int k = 0; k < 2000; ++k) {
        double t = k * dt;
        JointTrajectorySample s;
        s.q_d = Eigen::VectorXd(n);
        s.q = Eigen::VectorXd(n);

        s.q_d(0) = 0.5 * std::sin(2.0 * M_PI * t);
        s.q_d(1) = 0.3 * std::sin(4.0 * M_PI * t);

        s.q(0) = 0.48 * std::sin(2.0 * M_PI * (t - 0.01));
        s.q(1) = 0.29 * std::sin(4.0 * M_PI * (t - 0.01));

        traj.push_back(s);
    }

    double ise = computeISE(traj, dt);
    std::cout << "ISE = " << ise << std::endl;
    return 0;
}
      

9.2 Java (using EJML for linear algebra)

In Java, matrix libraries such as EJML provide efficient vector operations. Here, we sketch a simple error computation.

import org.ejml.simple.SimpleMatrix;
import java.util.List;

class JointSample {
    public SimpleMatrix qd;  // n x 1
    public SimpleMatrix q;   // n x 1
}

public class JointSpaceMetrics {

    public static double computeISE(List<JointSample> traj, double dt) {
        double ise = 0.0;
        for (JointSample s : traj) {
            SimpleMatrix e = s.qd.minus(s.q);
            double norm2 = e.dot(e);  // squared Euclidean norm
            ise += norm2 * dt;
        }
        return ise;
    }

    public static void main(String[] args) {
        // In practice: acquire samples from your robot middleware
        // and feed them into computeISE.
    }
}
      

9.3 MATLAB / Simulink

MATLAB (with Robotics System Toolbox) provides high-level interfaces to joint trajectories. The following script computes joint-space error and a simple ISE metric.

function [e, edot, ISE] = joint_space_error_matlab(qd, q, dt)
% qd, q : N-by-n arrays (desired and measured joint positions)
% dt    : sample time

arguments
    qd double
    q  double
    dt (1,1) double {mustBePositive}
end

assert(isequal(size(qd), size(q)), 'Size mismatch between qd and q');

e = qd - q;
edot           = zeros(size(e));
edot(2:end,:)  = diff(e) / dt;
edot(1,:)      = edot(2,:);

ISE = sum(sum(e.^2, 2)) * dt;
end

% Example usage:
% t  = 0:0.001:2;
% qd = [0.5*sin(2*pi*t') 0.3*sin(4*pi*t')];
% q  = [0.48*sin(2*pi*(t'-0.01)) 0.29*sin(4*pi*(t'-0.01))];
% [e, edot, ISE] = joint_space_error_matlab(qd, q, 0.001);
      

In Simulink, the same computations would be implemented using Sum, Gain, and Integrator blocks to construct \( e(t) \) and \( \dot{e}(t) \), followed by Math Function blocks to square and integrate the error.

9.4 Wolfram Mathematica

Mathematica is well suited for symbolic manipulation and continuous-time simulation of error dynamics. The snippet below simulates the second-order error template from Section 5.

(* Parameters for a single joint error dynamics *)
zeta   = 0.7;
omegaN = 8.0;

eq = e''[t] + 2 zeta omegaN e'[t] + omegaN^2 e[t] == 0;
ic = {e[0] == 0.1, e'[0] == 0};

sol = NDSolve[{eq, ic}, e, {t, 0, 2}][[1]];

Plot[Evaluate[e[t] /. sol], {t, 0, 2},
  AxesLabel -> {"t", "e(t)"},
  PlotLabel -> "Joint error dynamics"
]
      

Here, the objective is encoded entirely in the differential equation; a real controller must be designed so that the actual error trajectory matches this idealized solution as closely as possible.

In this lesson we rigorously defined joint-space control objectives for robot manipulators. Starting from the standard rigid-body dynamics, we introduced joint-space state variables, regulation and tracking goals, and error coordinates. We then specified desired second-order error dynamics in terms of natural frequency and damping ratio, connected these to classical performance indices (settling time, overshoot), and introduced quadratic performance functionals such as ISE. Finally, we illustrated how to compute joint-space error trajectories and associated metrics in Python, C++, Java, MATLAB/Simulink, and Mathematica. Subsequent lessons will design specific joint-space controllers (e.g., PD, PID with gravity compensation) to realize the objectives formalized here.