A parallel mechanism is described by a set of holonomic constraint equations linking the end-effector pose vector \( \mathbf{x}\in\mathbb{R}^m \) (task space) and the joint coordinates \( \mathbf{q}\in\mathbb{R}^n \) (usually actuator and passive joint variables). The constraints can be written as

\[ \boldsymbol{\Phi}(\mathbf{x},\mathbf{q}) = \mathbf{0}, \quad \boldsymbol{\Phi} : \mathbb{R}^m\times\mathbb{R}^n \rightarrow \mathbb{R}^k. \]

For a properly designed parallel robot, the number of independent constraints \( k \) is such that the instantaneous mobility of the platform is \( m \) (e.g., \( m=3 \) for planar and \( m=6 \) for spatial mechanisms). Singularities correspond to configurations where this mobility changes abruptly, leading to loss of stiffness, uncontrollable motions, or loss of force transmission.

In parallel mechanisms, three main singularity types are classically distinguished (Gosselin–Angeles classification):

• Type I (serial or inverse-kinematic singularity): related to the loss of rank of the actuation Jacobian.
• Type II (parallel or direct-kinematic singularity): related to the loss of rank of the constraint Jacobian with respect to platform coordinates.
• Type III (combined singularity): both matrices simultaneously lose rank, often associated with assembly-mode changing or self-motions.

At Type II singularities, the platform can move even with actuators locked (loss of stiffness). At Type I singularities, some platform twists cannot be generated even with actuators fully driven. Type III combines both degeneracies.

flowchart TD
  C["Configuration (x, q)"] --> J1["Compute A = dPhi/dx"]
  C --> J2["Compute B = dPhi/dq"]
  J1 --> T2["Check det(A) = 0 ?"]
  J2 --> T1["Check rank(B) drop ?"]
  T1 --> L1["Type I: \nloss of velocity / force control"]
  T2 --> L2["Type II: \nloss of stiffness, extra motion"]
  T1 --> BOTH["Both degenerate ?"]
  T2 --> BOTH
  BOTH --> L3["Type III: \ncombined / complex singularity"]
        

From Lesson 2, the constraint equations \( \boldsymbol{\Phi}(\mathbf{x},\mathbf{q})=\mathbf{0} \) yield the instantaneous kinematics upon differentiation with respect to time:

\[ \frac{\partial\boldsymbol{\Phi}}{\partial\mathbf{x}}(\mathbf{x},\mathbf{q})\, \dot{\mathbf{x}} + \frac{\partial\boldsymbol{\Phi}}{\partial\mathbf{q}}(\mathbf{x},\mathbf{q})\, \dot{\mathbf{q}} = \mathbf{0}. \]

Define the matrices

\[ \mathbf{A}(\mathbf{x},\mathbf{q}) := \frac{\partial\boldsymbol{\Phi}}{\partial\mathbf{x}}(\mathbf{x},\mathbf{q}), \qquad \mathbf{B}(\mathbf{x},\mathbf{q}) := \frac{\partial\boldsymbol{\Phi}}{\partial\mathbf{q}}(\mathbf{x},\mathbf{q}), \]

so that

\[ \mathbf{A}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{x}} + \mathbf{B}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{q}} = \mathbf{0}. \]

For a fully actuated parallel robot with \( m \) platform DoFs and \( m \) independent actuators, we typically have \( \mathbf{A},\mathbf{B}\in\mathbb{R}^{m\times m} \) and \( \operatorname{rank}(\mathbf{A}) = \operatorname{rank}(\mathbf{B}) = m \) in regular configurations.

If \( \mathbf{A} \) is invertible, the mapping from actuator velocities to platform twist is

\[ \dot{\mathbf{x}} = -\mathbf{A}^{-1}(\mathbf{x},\mathbf{q})\, \mathbf{B}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{q}} = \mathbf{J}_{xq}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{q}}, \]

where \( \mathbf{J}_{xq} = -\mathbf{A}^{-1}\mathbf{B} \) is the forward instantaneous kinematic Jacobian. Similarly, if \( \mathbf{B} \) is invertible we can express \( \dot{\mathbf{q}} = -\mathbf{B}^{-1}\mathbf{A}\dot{\mathbf{x}} \).

Singularities occur whenever these mappings cease to be locally invertible or lose rank:

• Loss of rank of \( \mathbf{A} \): unbounded or indeterminate platform motion for bounded actuator motion.
• Loss of rank of \( \mathbf{B} \): actuators cannot span the constraint space, leading to inability to generate some twists or wrenches.

Consider the static equilibrium of a parallel mechanism subjected to an external wrench \( \mathbf{w}\in\mathbb{R}^m \) on the platform and actuator torques/forces \( \boldsymbol{\tau}\in\mathbb{R}^m \). Let \( \boldsymbol{\lambda}\in\mathbb{R}^m \) denote the internal constraint forces. From the principle of virtual work with admissible virtual velocities \( \delta\mathbf{x},\delta\mathbf{q} \) satisfying \( \mathbf{A}\,\delta\mathbf{x} + \mathbf{B}\,\delta\mathbf{q} = \mathbf{0} \), one obtains

\[ \boldsymbol{\tau}^{\mathsf{T}}\delta\mathbf{q} + \mathbf{w}^{\mathsf{T}}\delta\mathbf{x} = \mathbf{0}. \]

Introducing the internal forces \( \boldsymbol{\lambda} \) as Lagrange multipliers for the constraints yields

\[ \mathbf{B}^{\mathsf{T}}(\mathbf{x},\mathbf{q})\,\boldsymbol{\lambda} = \boldsymbol{\tau}, \qquad \mathbf{A}^{\mathsf{T}}(\mathbf{x},\mathbf{q})\,\boldsymbol{\lambda} = -\mathbf{w}. \]

When both \( \mathbf{A} \) and \( \mathbf{B} \) are nonsingular, the mapping

\[ \boldsymbol{\tau} = -\mathbf{B}^{\mathsf{T}}\mathbf{A}^{-\mathsf{T}}\,\mathbf{w} \]

is well-defined and finite. This matrix \( \mathbf{K}_w := -\mathbf{B}^{\mathsf{T}}\mathbf{A}^{-\mathsf{T}} \) is the static Jacobian mapping platform wrenches to actuator torques/forces.

Now we can formalize the three types of singularities.

3.1 Type II (Parallel) Singularities

A Type II singularity occurs when \( \mathbf{A} \) loses rank: \( \operatorname{rank}(\mathbf{A}) < m \). In that case, there exist nonzero platform twists \( \dot{\mathbf{x}}\neq\mathbf{0} \) such that

\[ \mathbf{A}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{x}} = \mathbf{0}, \]

with \( \dot{\mathbf{q}} = \mathbf{0} \). The platform can move despite locked actuators (loss of stiffness, self-motions). From the static perspective, as \( \det(\mathbf{A}) \rightarrow 0 \), the matrix \( \mathbf{A}^{-\mathsf{T}} \) becomes unbounded and the mapping \( \mathbf{w}\mapsto\boldsymbol{\tau} \) may demand arbitrarily large actuator efforts.

Formally, if \( \mathbf{A} \) is singular, there exists \( \mathbf{w}\neq\mathbf{0} \) with \( \mathbf{A}^{\mathsf{T}}\boldsymbol{\lambda}=-\mathbf{w} \) having no solution or infinitely many solutions, which breaks the uniqueness of actuator reactions.

3.2 Type I (Serial) Singularities

A Type I singularity occurs when \( \mathbf{B} \) loses rank: \( \operatorname{rank}(\mathbf{B}) < m \). In this case, there exist nonzero joint velocity vectors \( \dot{\mathbf{q}}\neq\mathbf{0} \) that are compatible with \( \dot{\mathbf{x}} = \mathbf{0} \):

\[ \mathbf{B}(\mathbf{x},\mathbf{q})\,\dot{\mathbf{q}} = \mathbf{0}, \qquad \dot{\mathbf{x}} = \mathbf{0}. \]

From the actuator side, it means the actuators cannot span all required constraint reactions. Certain platform wrenches cannot be balanced using finite \( \boldsymbol{\tau} \). Geometrically, some leg chains are themselves in serial singularity, analogous to the kinematic singularity of a single open chain.

3.3 Type III (Combined) Singularities

A Type III singularity arises when both \( \mathbf{A} \) and \( \mathbf{B} \) lose rank simultaneously. These configurations are typically associated with:

• changes of assembly mode (multiple solutions to the forward kinematics coalesce),
• transition between workspace regions, and
• rich sets of self-motions.

Mathematically, there exist nontrivial pairs \( (\dot{\mathbf{x}},\dot{\mathbf{q}}) \) satisfying \( \mathbf{A}\dot{\mathbf{x}}+\mathbf{B}\dot{\mathbf{q}}=\mathbf{0} \) with both null spaces nontrivial, leading to underconstrained internal motions.

Consider a planar 3-RPR parallel robot: three identical legs connect a fixed base triangle to a moving platform triangle. Each leg consists of a revolute (R) joint at the base, an actuated prismatic (P) joint, and a revolute (R) joint at the platform.

4.1 Geometry and Constraint Equations

Let the base anchor points be \( \mathbf{A}_i = [A_{ix},A_{iy}]^{\mathsf{T}} \) and the platform anchor points in the platform frame be \( \mathbf{b}_i = [r_i\cos\alpha_i,\,r_i\sin\alpha_i]^{\mathsf{T}} \), \( i=1,2,3 \). The platform pose is \( \mathbf{x} = [x,y,\varphi]^{\mathsf{T}} \) (planar position and orientation). The prismatic joint lengths are \( \boldsymbol{\rho}=[\rho_1,\rho_2,\rho_3]^{\mathsf{T}} \).

The position of platform point \( i \) in the base frame is

\[ \mathbf{p}_i(\mathbf{x}) = \begin{bmatrix} x \\ y \end{bmatrix} + \mathbf{R}(\varphi)\,\mathbf{b}_i, \quad \mathbf{R}(\varphi) = \begin{bmatrix} \cos\varphi & -\sin\varphi \\ \sin\varphi & \cos\varphi \end{bmatrix}. \]

The (squared) length of leg \( i \) must equal \( \rho_i^2 \), hence the constraint equations are

\[ \Phi_i(\mathbf{x},\boldsymbol{\rho}) = \big\|\mathbf{p}_i(\mathbf{x})-\mathbf{A}_i\big\|^2 - \rho_i^2 = 0, \quad i=1,2,3. \]

Collecting them gives \( \boldsymbol{\Phi}(\mathbf{x},\boldsymbol{\rho}) = [\Phi_1,\Phi_2,\Phi_3]^{\mathsf{T}} = \mathbf{0} \).

4.2 Jacobians \( \mathbf{A} \) and \( \mathbf{B} \)

Differentiating each constraint \( \Phi_i \) with respect to time:

\[ \frac{\partial\Phi_i}{\partial\mathbf{x}}\dot{\mathbf{x}} + \frac{\partial\Phi_i}{\partial\rho_i}\dot{\rho}_i = 0. \]

Since \( \Phi_i = (\mathbf{p}_i-\mathbf{A}_i)^{\mathsf{T}} (\mathbf{p}_i-\mathbf{A}_i) - \rho_i^2 \), we have

\[ \frac{\partial\Phi_i}{\partial\mathbf{p}_i} = 2(\mathbf{p}_i-\mathbf{A}_i)^{\mathsf{T}},\qquad \frac{\partial\Phi_i}{\partial\rho_i} = -2\rho_i. \]

Using the chain rule, \( \frac{\partial\Phi_i}{\partial\mathbf{x}} = \frac{\partial\Phi_i}{\partial\mathbf{p}_i}\, \frac{\partial\mathbf{p}_i}{\partial\mathbf{x}} \). Writing \( \mathbf{x}=[x,y,\varphi]^{\mathsf{T}} \), we obtain

\[ \frac{\partial\mathbf{p}_i}{\partial\mathbf{x}} = \begin{bmatrix} 1 & 0 & -r_i\sin(\varphi+\alpha_i) \\ 0 & 1 & r_i\cos(\varphi+\alpha_i) \end{bmatrix}. \]

Therefore the \( i \)-th row of \( \mathbf{A} \) is

\[ \mathbf{A}_{i\cdot} = \frac{\partial\Phi_i}{\partial\mathbf{x}} = 2(\mathbf{p}_i-\mathbf{A}_i)^{\mathsf{T}} \frac{\partial\mathbf{p}_i}{\partial\mathbf{x}} \in\mathbb{R}^{1\times 3}. \]

The matrix \( \mathbf{B} \in\mathbb{R}^{3\times 3} \) is diagonal:

\[ \mathbf{B}(\mathbf{x},\boldsymbol{\rho}) = \operatorname{diag}\big(-2\rho_1,\,-2\rho_2,\,-2\rho_3\big). \]

Hence:

• Type I singularities occur when at least one \( \rho_i = 0 \), i.e., \( \det(\mathbf{B})=0 \).
• Type II singularities occur when rows of \( \mathbf{A} \) become linearly dependent, so that \( \det(\mathbf{A}) = 0 \).

4.3 Geometric Interpretation of Type II Singularities

For the 3-RPR robot, the direction of leg \( i \) is along \( \mathbf{u}_i = \frac{\mathbf{p}_i-\mathbf{A}_i}{\|\mathbf{p}_i-\mathbf{A}_i\|} \). Under mild assumptions, it can be shown that the rows of \( \mathbf{A} \) become linearly dependent if and only if the three leg directions \( \mathbf{u}_1,\mathbf{u}_2,\mathbf{u}_3 \) are either concurrent (their lines intersect in a common point) or parallel. In that case, the platform admits a finite instantaneous motion (rotation about the concurrence point or translation) without changing any \( \rho_i \), i.e., a self-motion.

A sketch of the argument is:

1. The \( i \)-th row of \( \mathbf{A} \) is proportional to the screw parameters of the constraint imposed by leg \( i \).
2. Concurrency or parallelism of leg axes implies a linear dependence among the corresponding constraint screws.
3. Linear dependence of constraint screws equals a rank deficiency of \( \mathbf{A} \), which generates extra degrees of freedom.

Thus, for workspace design and path planning, one must avoid configurations where leg axes nearly intersect at a point or become nearly parallel.

In practice, exact singularity conditions \( \det(\mathbf{A}) = 0 \) or \( \operatorname{rank}(\mathbf{B}) < m \) are rarely met due to numerical uncertainties. Instead, we measure proximity to singularities.

5.1 Determinant-Based Measures

For square matrices, the determinant magnitude gives a crude measure:

\[ \eta_A(\mathbf{x},\mathbf{q}) = |\det(\mathbf{A})|,\qquad \eta_B(\mathbf{x},\mathbf{q}) = |\det(\mathbf{B})|. \]

Small values of \( \eta_A \) signal proximity to Type II singularities, whereas small \( \eta_B \) indicate proximity to Type I singularities. However, determinants are scale-dependent and may be numerically ill-conditioned for high-dimensional problems.

5.2 Condition Number and Singular Values

A more robust measure uses the singular value decomposition (SVD). For a matrix \( \mathbf{M}\in\mathbb{R}^{m\times m} \), let \( \sigma_{\max} \) and \( \sigma_{\min} \) denote its largest and smallest singular values. The 2-norm condition number is

\[ \kappa_2(\mathbf{M}) = \frac{\sigma_{\max}(\mathbf{M})}{\sigma_{\min}(\mathbf{M})}. \]

Near singularities, \( \sigma_{\min}(\mathbf{M})\rightarrow 0 \) and \( \kappa_2(\mathbf{M})\rightarrow\infty \). We can define

\[ \kappa_A := \kappa_2(\mathbf{A}), \qquad \kappa_B := \kappa_2(\mathbf{B}), \]

and use thresholds (e.g., \( \kappa_A \) larger than some value) to classify configurations as numerically near-singular.

5.3 Manipulability-Like Measures

For velocity transmission, we may examine the Jacobian \( \mathbf{J}_{xq} = -\mathbf{A}^{-1}\mathbf{B} \) and define a manipulability measure analogous to that for serial manipulators:

\[ w(\mathbf{x},\mathbf{q}) = \sqrt{\det\big(\mathbf{J}_{xq}\mathbf{J}_{xq}^{\mathsf{T}}\big)} = \prod_{i=1}^m \sigma_i(\mathbf{J}_{xq}), \]

where \( \sigma_i \) are singular values of \( \mathbf{J}_{xq} \). When any singular value vanishes, \( w(\mathbf{x},\mathbf{q}) = 0 \), indicating a singularity. A small manipulability value means poor velocity transmission or force transmission in certain directions.

For computation, we follow a systematic procedure: construct the constraint Jacobians \( \mathbf{A},\mathbf{B} \), evaluate rank and condition numbers, and classify the configuration.

flowchart TD
  S["Start: pose x, joints q"] --> C["Build Phi(x, q)"]
  C --> J["Differentiate to get A = dPhi/dx, B = dPhi/dq"]
  J --> R1["Compute rank(A), cond(A)"]
  J --> R2["Compute rank(B), cond(B)"]
  R1 --> T2["Is rank(A) < m or cond(A) big ?"]
  R2 --> T1["Is rank(B) < m or cond(B) big ?"]
  T2 --> L2["Type II candidate"]
  T1 --> L1["Type I candidate"]
  T1 --> BOTH["Both conditions met ?"]
  T2 --> BOTH
  BOTH --> L3["Type III candidate / re-check geometry"]
        

In practice, the rank tests use singular values with a numerical tolerance. The geometry of the mechanism helps confirm whether a numerically near-singular configuration is structurally meaningful or just a numerical artifact.

We now implement the 3-RPR model in Python using numpy. The script constructs \( \mathbf{A} \) and \( \mathbf{B} \) for a given configuration, computes determinants and condition numbers, and classifies the singularity type.

import numpy as np

# ---------- Geometry of a symmetric 3-RPR robot ----------
# Base anchor points (equilateral triangle of radius Rb)
Rb = 1.0
Ab = np.array([
    [Rb, 0.0],
    [Rb * np.cos(2.0 * np.pi / 3.0), Rb * np.sin(2.0 * np.pi / 3.0)],
    [Rb * np.cos(4.0 * np.pi / 3.0), Rb * np.sin(4.0 * np.pi / 3.0)]
])

# Platform anchor points (equilateral triangle of radius Rp) in platform frame
Rp = 0.5
alpha = np.array([0.0, 2.0 * np.pi / 3.0, 4.0 * np.pi / 3.0])
b_local = np.stack([Rp * np.cos(alpha), Rp * np.sin(alpha)], axis=1)

def forward_points(x_vec):
    """
    Compute platform points p_i in base frame.
    x_vec = [x, y, phi].
    """
    x, y, phi = x_vec
    R = np.array([[np.cos(phi), -np.sin(phi)],
                  [np.sin(phi),  np.cos(phi)]])
    p = np.empty_like(b_local)
    for i in range(3):
        p[i, :] = np.array([x, y]) + R @ b_local[i, :]
    return p

def compute_AB(x_vec, rho):
    """
    Build A and B matrices for the 3-RPR mechanism.
    A: 3x3, B: 3x3 diagonal.
    """
    p = forward_points(x_vec)
    A = np.zeros((3, 3))
    B = np.zeros((3, 3))

    x, y, phi = x_vec
    for i in range(3):
        # leg vector and its norm
        d = p[i, :] - Ab[i, :]
        # partial of p_i w.r.t. x, y, phi
        ri = Rp
        ang = phi + alpha[i]
        dp_dxi = np.array([[1.0, 0.0, -ri * np.sin(ang)],
                           [0.0, 1.0,  ri * np.cos(ang)]])
        # A row: 2 (p_i - A_i)^T * dp_dxi
        A[i, :] = 2.0 * d @ dp_dxi

        # B diagonal: dPhi_i/drho_i = -2 rho_i
        B[i, i] = -2.0 * rho[i]

    return A, B

def classify_singularity(A, B, tol=1e-6):
    """
    Classify configuration based on A and B.
    """
    # determinants (may be near zero)
    detA = np.linalg.det(A)
    detB = np.linalg.det(B)

    # singular values for rank and condition number
    sA = np.linalg.svd(A, compute_uv=False)
    sB = np.linalg.svd(B, compute_uv=False)

    rankA = np.sum(sA > tol)
    rankB = np.sum(sB > tol)

    kappaA = (sA[0] / sA[-1]) if sA[-1] > tol else np.inf
    kappaB = (sB[0] / sB[-1]) if sB[-1] > tol else np.inf

    m = A.shape[0]
    type_I = rankB < m
    type_II = rankA < m

    if type_I and type_II:
        singular_type = "Type III (combined)"
    elif type_I:
        singular_type = "Type I (serial)"
    elif type_II:
        singular_type = "Type II (parallel)"
    else:
        singular_type = "Regular"

    return {
        "type": singular_type,
        "detA": detA,
        "detB": detB,
        "rankA": int(rankA),
        "rankB": int(rankB),
        "kappaA": float(kappaA),
        "kappaB": float(kappaB),
    }

# ---------- Example usage ----------
if __name__ == "__main__":
    # Example configuration: center pose with moderate orientation
    x_vec = np.array([0.0, 0.0, 0.3])
    rho = np.array([1.0, 1.0, 1.0])

    A, B = compute_AB(x_vec, rho)
    info = classify_singularity(A, B)

    print("A =\n", A)
    print("B =\n", B)
    print("Classification:", info["type"])
    print("det(A) =", info["detA"], "det(B) =", info["detB"])
    print("rank(A) =", info["rankA"], "rank(B) =", info["rankB"])
    print("kappa(A) =", info["kappaA"], "kappa(B) =", info["kappaB"])
      

This code uses numpy for matrix computations. A more advanced toolbox could integrate symbolic differentiation (sympy) to generate \( \mathbf{A},\mathbf{B} \) automatically from the constraint equations \( \boldsymbol{\Phi}(\mathbf{x},\mathbf{q}) \).

In C++, the Eigen library is a standard choice for linear algebra in robotics. The snippet below assumes that the matrices \( \mathbf{A},\mathbf{B} \) have been computed (e.g., using the same formulas as in the Python code).

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

using Eigen::Matrix3d;
using Eigen::JacobiSVD;

struct SingularityInfo {
    std::string type;
    double detA;
    double detB;
    int rankA;
    int rankB;
    double kappaA;
    double kappaB;
};

SingularityInfo classifySingularity(const Matrix3d& A,
                                    const Matrix3d& B,
                                    double tol = 1e-6)
{
    SingularityInfo info;

    info.detA = A.determinant();
    info.detB = B.determinant();

    JacobiSVD<Matrix3d> svdA(A);
    JacobiSVD<Matrix3d> svdB(B);

    Eigen::Vector3d sA = svdA.singularValues();
    Eigen::Vector3d sB = svdB.singularValues();

    info.rankA = (sA.array() > tol).count();
    info.rankB = (sB.array() > tol).count();

    info.kappaA = (sA(2) > tol) ? sA(0) / sA(2) : std::numeric_limits<double>::infinity();
    info.kappaB = (sB(2) > tol) ? sB(0) / sB(2) : std::numeric_limits<double>::infinity();

    int m = 3;
    bool typeI  = info.rankB < m;
    bool typeII = info.rankA < m;

    if (typeI && typeII) {
        info.type = "Type III (combined)";
    } else if (typeI) {
        info.type = "Type I (serial)";
    } else if (typeII) {
        info.type = "Type II (parallel)";
    } else {
        info.type = "Regular";
    }

    return info;
}

int main() {
    Matrix3d A, B;
    // Fill A and B here (e.g., using a function similar to compute_AB in Python)

    // Dummy example:
    A.setIdentity();
    B.setIdentity();

    SingularityInfo info = classifySingularity(A, B);

    std::cout << "Type: " << info.type << std::endl;
    std::cout << "det(A) = " << info.detA
              << ", det(B) = " << info.detB << std::endl;
    std::cout << "rank(A) = " << info.rankA
              << ", rank(B) = " << info.rankB << std::endl;
    std::cout << "kappa(A) = " << info.kappaA
              << ", kappa(B) = " << info.kappaB << std::endl;

    return 0;
}
      

In a full robotics codebase, this classification would be embedded within a kinematics module, with A and B generated by an analytic kinematics engine or symbolic preprocessing step.

For Java, the Apache Commons Math library provides matrix operations and SVD. The code below assumes that the entries of \( \mathbf{A} \) and \( \mathbf{B} \) are stored in a two-dimensional double array.

import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularValueDecomposition;

public class ParallelSingularity {

    public static class SingularityInfo {
        public String type;
        public double detA;
        public double detB;
        public int rankA;
        public int rankB;
        public double kappaA;
        public double kappaB;
    }

    public static SingularityInfo classify(double[][] Adata,
                                           double[][] Bdata,
                                           double tol) {
        SingularityInfo info = new SingularityInfo();

        RealMatrix A = MatrixUtils.createRealMatrix(Adata);
        RealMatrix B = MatrixUtils.createRealMatrix(Bdata);

        info.detA = new org.apache.commons.math3.linear.LUDecomposition(A).getDeterminant();
        info.detB = new org.apache.commons.math3.linear.LUDecomposition(B).getDeterminant();

        SingularValueDecomposition svdA = new SingularValueDecomposition(A);
        SingularValueDecomposition svdB = new SingularValueDecomposition(B);

        double[] sA = svdA.getSingularValues();
        double[] sB = svdB.getSingularValues();

        info.rankA = 0;
        info.rankB = 0;
        for (double v : sA) {
            if (v > tol) info.rankA++;
        }
        for (double v : sB) {
            if (v > tol) info.rankB++;
        }

        info.kappaA = (sA[sA.length - 1] > tol)
                      ? sA[0] / sA[sA.length - 1]
                      : Double.POSITIVE_INFINITY;
        info.kappaB = (sB[sB.length - 1] > tol)
                      ? sB[0] / sB[sB.length - 1]
                      : Double.POSITIVE_INFINITY;

        int m = A.getRowDimension();
        boolean typeI  = info.rankB < m;
        boolean typeII = info.rankA < m;

        if (typeI && typeII) {
            info.type = "Type III (combined)";
        } else if (typeI) {
            info.type = "Type I (serial)";
        } else if (typeII) {
            info.type = "Type II (parallel)";
        } else {
            info.type = "Regular";
        }

        return info;
    }

    public static void main(String[] args) {
        double[][] A = { {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
        double[][] B = { {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };

        SingularityInfo info = classify(A, B, 1e-6);
        System.out.println("Type: " + info.type);
        System.out.println("det(A) = " + info.detA + ", det(B) = " + info.detB);
        System.out.println("rank(A) = " + info.rankA + ", rank(B) = " + info.rankB);
        System.out.println("kappa(A) = " + info.kappaA + ", kappa(B) = " + info.kappaB);
    }
}
      

Within a Java-based robotics framework, the constraint equations for a specific parallel robot would generate the matrices A and B, which are then passed to such a classifier.

MATLAB provides native support for matrix operations and visualization. Below is a script that scans a grid of platform poses for the 3-RPR robot and records where \( \det(\mathbf{A}) \) is close to zero.

function singularity_map_3RPR()
    % Geometry
    Rb = 1.0;
    Ab = [ Rb, 0.0;
           Rb*cos(2*pi/3), Rb*sin(2*pi/3);
           Rb*cos(4*pi/3), Rb*sin(4*pi/3) ];

    Rp    = 0.5;
    alpha = [0.0; 2*pi/3; 4*pi/3];
    b_loc = [Rp*cos(alpha), Rp*sin(alpha)];

    % Grid in x-y plane at fixed orientation
    phi  = 0.3;
    xs   = linspace(-0.3, 0.3, 40);
    ys   = linspace(-0.3, 0.3, 40);
    detA = zeros(length(ys), length(xs));

    % Nominal leg lengths (for example)
    rho_nom = [1.0; 1.0; 1.0];

    for ix = 1:length(xs)
        for iy = 1:length(ys)
            x_vec = [xs(ix); ys(iy); phi];
            [A, ~] = computeAB_3RPR(x_vec, rho_nom, Ab, b_loc);
            detA(iy, ix) = det(A);
        end
    end

    figure;
    imagesc(xs, ys, log10(abs(detA)));
    set(gca, "YDir", "normal");
    colorbar;
    title("log_{10}(|det(A)|) for 3-RPR");
    xlabel("x");
    ylabel("y");
end

function [A, B] = computeAB_3RPR(x_vec, rho, Ab, b_loc)
    x   = x_vec(1);
    y   = x_vec(2);
    phi = x_vec(3);

    R = [cos(phi), -sin(phi);
         sin(phi),  cos(phi)];

    p = zeros(3, 2);
    for i = 1:3
        p(i, :) = [x, y] + (R * b_loc(i, :).').';
    end

    A = zeros(3, 3);
    B = zeros(3, 3);

    Rp = norm(b_loc(1, :));  % same for each
    for i = 1:3
        d = p(i, :) - Ab(i, :);
        ang = phi + alpha(i); %#ok<NASGU> % can be precomputed

        dp_dxi = [1.0, 0.0, -Rp * sin(phi + alpha(i));
                  0.0, 1.0,  Rp * cos(phi + alpha(i))];

        A(i, :) = 2.0 * d * dp_dxi;
        B(i, i) = -2.0 * rho(i);
    end
end
      

In Simulink, a MATLAB Function block can call computeAB_3RPR at each simulation step using the current state vector x_vec and rho. A separate block can compute det(A) and feed it to a Scope or Display block, allowing online monitoring of proximity to Type II singularities along a simulated trajectory.

Wolfram Mathematica offers symbolic differentiation, making it convenient to derive \( \mathbf{A},\mathbf{B} \) directly from the constraints \( \boldsymbol{\Phi}(\mathbf{x},\mathbf{q}) \).

(* Variables *)
Clear[x, y, phi, rho1, rho2, rho3, Rb, Rp, alpha1, alpha2, alpha3];
xVec   = {x, y, phi};
rhoVec = {rho1, rho2, rho3};

(* Base anchors (parameters) *)
A1 = {Rb, 0};
A2 = {Rb*Cos[2*Pi/3], Rb*Sin[2*Pi/3]};
A3 = {Rb*Cos[4*Pi/3], Rb*Sin[4*Pi/3]};

(* Platform anchors in platform frame *)
b1 = {Rp*Cos[alpha1], Rp*Sin[alpha1]};
b2 = {Rp*Cos[alpha2], Rp*Sin[alpha2]};
b3 = {Rp*Cos[alpha3], Rp*Sin[alpha3]};

R[theta_] := { {Cos[theta], -Sin[theta]},
              {Sin[theta],  Cos[theta]} };

p1 = {x, y} + R[phi].b1;
p2 = {x, y} + R[phi].b2;
p3 = {x, y} + R[phi].b3;

Phi1 = (p1 - A1).(p1 - A1) - rho1^2;
Phi2 = (p2 - A2).(p2 - A2) - rho2^2;
Phi3 = (p3 - A3).(p3 - A3) - rho3^2;

PhiVec = {Phi1, Phi2, Phi3};

(* Jacobians A and B *)
A = D[PhiVec, {{x, y, phi}}];   (* dPhi/dx *)
B = D[PhiVec, {rhoVec}];        (* dPhi/drho *)

DetA = Simplify[Det[A]];

(* Example numerical substitution *)
params = {Rb -> 1.0, Rp -> 0.5,
          alpha1 -> 0, alpha2 -> 2*Pi/3, alpha3 -> 4*Pi/3,
          x -> 0.0, y -> 0.0, phi -> 0.3,
          rho1 -> 1.0, rho2 -> 1.0, rho3 -> 1.0};

A_num   = A /. params;
B_num   = B /. params;
detAVal = Det[A_num];
detBVal = Det[B_num];
      

With symbolic expressions for \( \det(\mathbf{A}) \), one can analytically derive the geometry of the Type II singularity curves in the workspace, for instance by solving \( \det(\mathbf{A}) = 0 \) for \( x,y \) at fixed orientation \( \varphi \).

In this lesson we formalized singularities in parallel mechanisms through the constraint Jacobians \( \mathbf{A} = \partial\boldsymbol{\Phi}/\partial\mathbf{x} \) and \( \mathbf{B} = \partial\boldsymbol{\Phi}/\partial\mathbf{q} \). Type II singularities arise when \( \mathbf{A} \) loses rank, causing self-motions and loss of stiffness; Type I singularities correspond to rank deficiency of \( \mathbf{B} \), limiting the ability of actuators to control platform motion; Type III singularities combine both phenomena. A planar 3-RPR parallel manipulator served as a concrete example, illustrating how determinant and condition number metrics reveal proximity to singularity sets. Finally, we implemented singularity detection and classification using Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, bridging the theoretical definitions with practical computation that will be essential for workspace design and safe trajectory generation in subsequent lessons.