A classical Stewart platform is a 6-DOF parallel manipulator consisting of a fixed base platform and a moving top platform connected by six variable-length legs (typically prismatic actuators) arranged in a 6–6 configuration. Each leg is attached to the base and the moving platform using universal or spherical joints, forming six independent kinematic chains in parallel.

We introduce two right-handed Cartesian frames:

• Base frame \( \{B\} \), fixed to the base platform.
• Moving platform frame \( \{P\} \), fixed to the moving platform, usually at its geometric center.

For leg \( i \in \{1,\dots,6\} \):

• \( \mathbf{b}_i \in \mathbb{R}^3 \): coordinates of the base attachment point in frame \( \{B\} \).
• \( \mathbf{p}_i \in \mathbb{R}^3 \): coordinates of the corresponding platform attachment point in frame \( \{P\} \).
• \( L_i \): actuated leg length (prismatic joint variable).

The platform pose relative to the base is described by:

• \( \mathbf{p} \in \mathbb{R}^3 \): position vector of the origin of \(\{P\}\) in \(\{B\}\).
• \( \mathbf{R} \in SO(3) \): rotation matrix from frame \(\{P\}\) to frame \(\{B\}\), often parameterized by Euler angles or unit quaternions (introduced earlier in the course).

Our objectives in this case study are:

1. Derive geometric constraint equations relating \(\mathbf{p},\mathbf{R}\) and leg lengths \(L_i\).
2. Formulate inverse and forward kinematics for the platform.
3. Obtain the differential kinematics (Jacobian) mapping platform twist to leg length rates.
4. Use virtual work to relate actuator forces to the platform wrench and embed this in a dynamic model.
5. Implement these models in multiple programming environments with robotics-relevant libraries.

flowchart TD
  A["Specify geometry: base points Bi, platform points Pi"] --> B["Choose pose: R, p"]
  B --> C["Compute leg vectors di = p + R*Pi - Bi"]
  C --> D["Compute leg lengths Li = ||di||"]
  D --> E["Check constraints: limits, collisions, workspace"]
  E --> F["Use in IK, dynamics, and control"]
        

The key geometric relationship comes from the fact that each leg connects corresponding attachment points on the base and moving platform. In base coordinates, the position of the \(i\)-th platform attachment point is

\[ \,^{B}\mathbf{x}_i = \mathbf{p} + \mathbf{R}\,\mathbf{p}_i, \quad i = 1,\dots,6. \]

The leg vector from the base to the moving platform is

\[ \mathbf{d}_i = \,^{B}\mathbf{x}_i - \mathbf{b}_i \;=\; \mathbf{p} + \mathbf{R}\,\mathbf{p}_i - \mathbf{b}_i. \]

The corresponding leg length is

\[ L_i = \|\mathbf{d}_i\| \;=\; \sqrt{\mathbf{d}_i^{\top}\mathbf{d}_i} \;=\; \sqrt{\big(\mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i\big)^{\top} \big(\mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i\big)}. \]

Squaring both sides eliminates the square root and yields six scalar constraint equations:

\[ \Phi_i(\mathbf{p},\mathbf{R},L_i) \;=\; L_i^2 - \big(\mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i\big)^{\top} \big(\mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i\big) \;=\; 0,\quad i = 1,\dots,6. \]

These constraints encode the closed-chain geometry of the parallel mechanism. They are the basis for both:

• Inverse kinematics (IK): given \(\mathbf{p},\mathbf{R}\), compute the actuator lengths \(L_i\).
• Forward kinematics (FK): given leg lengths \(L_i\), find the pose \((\mathbf{p},\mathbf{R})\) that satisfies all six constraints.

To obtain a minimal parameterization, we collect the generalized coordinates for the platform pose in

\[ \mathbf{x} = \begin{bmatrix} \mathbf{p} \\ \boldsymbol{\theta} \end{bmatrix} \in \mathbb{R}^6, \quad \boldsymbol{\theta} = [\phi,\theta,\psi]^{\top} \]

where \(\boldsymbol{\theta}\) denotes a set of Euler angles (e.g., ZYX) describing \(\mathbf{R}(\boldsymbol{\theta})\), already introduced in earlier chapters.

3.1 Inverse kinematics (platform pose → leg lengths)

For IK, the pose \((\mathbf{p},\mathbf{R})\) is known from the task-space specification, and we seek the six leg lengths \(L_i\). Using the geometric relationships from the previous section, the IK is straightforward:

\[ L_i(\mathbf{p},\mathbf{R}) = \big\|\mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i\big\|, \quad i = 1,\dots,6. \]

This closed-form computation is one of the main advantages of many parallel mechanisms: IK is simple and computationally cheap.

3.2 Forward kinematics (leg lengths → platform pose)

In FK, we know the actuator lengths \(L_i\) and must solve for the pose \(\mathbf{x}\). We rewrite the six constraint equations as a nonlinear system:

\[ F(\mathbf{x}) = \begin{bmatrix} \Phi_1(\mathbf{x}) \\ \vdots \\ \Phi_6(\mathbf{x}) \end{bmatrix} = \mathbf{0}, \quad \Phi_i(\mathbf{x}) = L_i^2 - \big\|\mathbf{p} + \mathbf{R}(\boldsymbol{\theta})\mathbf{p}_i - \mathbf{b}_i\big\|^2. \]

The FK problem is significantly more challenging:

• It is a system of six coupled, nonlinear equations in six unknowns.
• There can be multiple valid solutions (assembly modes) for the same set of leg lengths.
• Closed-form algebraic solutions exist only for specially structured geometries; generic platforms require numerical methods.

A standard numerical solution uses Newton–Raphson iterations:

\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \big(\mathbf{J}_F(\mathbf{x}_k)\big)^{-1} F(\mathbf{x}_k), \quad \mathbf{J}_F(\mathbf{x}) = \frac{\partial F}{\partial \mathbf{x}}(\mathbf{x}) \in \mathbb{R}^{6 \times 6}. \]

Here \(\mathbf{J}_F\) is the Jacobian of the constraint vector with respect to the pose parameters \(\mathbf{x}\). In practice, one may compute \(\mathbf{J}_F\) analytically, or approximate it numerically using finite differences.

To use the platform in motion control, we need a relationship between leg length rates \( \dot{L}_i \) and the platform twist (linear and angular velocity). Define the platform twist

\[ \mathbf{t} = \begin{bmatrix} \mathbf{v} \\ \boldsymbol{\omega} \end{bmatrix} \in \mathbb{R}^6, \]

where \(\mathbf{v}\) is the translational velocity of the origin of \(\{P\}\) expressed in \(\{B\}\), and \(\boldsymbol{\omega}\) is the angular velocity of the platform.

From the geometry, recall

\[ \mathbf{d}_i = \mathbf{p} + \mathbf{R}\mathbf{p}_i - \mathbf{b}_i, \quad L_i = \|\mathbf{d}_i\|, \quad \hat{\mathbf{u}}_i = \frac{\mathbf{d}_i}{L_i} \]

Differentiating \(\mathbf{d}_i\) with respect to time:

\[ \dot{\mathbf{d}}_i = \dot{\mathbf{p}} + \dot{\mathbf{R}}\mathbf{p}_i = \mathbf{v} + \boldsymbol{\omega} \times \big(\mathbf{R}\mathbf{p}_i\big), \]

where we used the standard rigid-body relationship \( \dot{\mathbf{R}}\mathbf{p}_i = \boldsymbol{\omega} \times (\mathbf{R}\mathbf{p}_i) \).

Differentiating the leg length gives

\[ \dot{L}_i = \frac{\mathbf{d}_i^{\top}\dot{\mathbf{d}}_i}{\|\mathbf{d}_i\|} = \hat{\mathbf{u}}_i^{\top} \dot{\mathbf{d}}_i = \hat{\mathbf{u}}_i^{\top} \mathbf{v} + \hat{\mathbf{u}}_i^{\top} \big( \boldsymbol{\omega} \times (\mathbf{R}\mathbf{p}_i) \big). \]

Using the vector triple product identity \( \mathbf{a}^{\top} (\boldsymbol{\omega} \times \mathbf{b}) = (\mathbf{b} \times \mathbf{a})^{\top} \boldsymbol{\omega} \), we obtain

\[ \dot{L}_i = \hat{\mathbf{u}}_i^{\top} \mathbf{v} + \big(\mathbf{R}\mathbf{p}_i \times \hat{\mathbf{u}}_i\big)^{\top}\boldsymbol{\omega}. \]

Therefore the \(i\)-th row of the length-rate Jacobian \(\mathbf{J}_L(\mathbf{p},\mathbf{R}) \in \mathbb{R}^{6 \times 6}\) is

\[ \mathbf{J}_L^{(i,:)} = \begin{bmatrix} \hat{\mathbf{u}}_i^{\top} & \big(\mathbf{R}\mathbf{p}_i \times \hat{\mathbf{u}}_i\big)^{\top} \end{bmatrix}. \]

Stacking all six legs gives

\[ \dot{\mathbf{L}} = \mathbf{J}_L(\mathbf{p},\mathbf{R}) \,\mathbf{t}, \quad \dot{\mathbf{L}} = \begin{bmatrix} \dot{L}_1 \\ \vdots \\ \dot{L}_6 \end{bmatrix}. \]

If \(\mathbf{J}_L\) is nonsingular, we can invert it locally:

\[ \mathbf{t} = \mathbf{J}_L^{-1}(\mathbf{p},\mathbf{R})\,\dot{\mathbf{L}}, \]

providing the mapping from leg length rates to platform twist. Singularities of the parallel mechanism correspond to configurations where \(\mathbf{J}_L\) loses rank, and either:

• a platform motion exists with \(\dot{\mathbf{L}} = \mathbf{0}\) (uncontrollable internal motion), or
• not all platform twists can be generated, even with arbitrary \(\dot{\mathbf{L}}\) (loss of mobility).

Consider the moving platform as a rigid body of mass \( m \) and inertia tensor \( \mathbf{I}_P \) expressed at the origin of \(\{P\}\). Let the generalized platform coordinates be \(\mathbf{x}\) and the twist \(\mathbf{t}\), as defined previously.

The rigid-body equations of motion of the platform (ignoring leg masses for simplicity) can be written in generalized form as

\[ \mathbf{M}(\mathbf{x})\dot{\mathbf{t}} + \mathbf{C}(\mathbf{x},\mathbf{t})\mathbf{t} + \mathbf{g}(\mathbf{x}) = \mathbf{w}_{\text{ext}} + \mathbf{w}_{\text{act}}, \]

where:

• \( \mathbf{M}(\mathbf{x}) \): 6 × 6 inertia matrix of the platform.
• \( \mathbf{C}(\mathbf{x},\mathbf{t}) \): Coriolis/centrifugal terms.
• \( \mathbf{g}(\mathbf{x}) \): gravitational wrench expressed in \(\{B\}\).
• \( \mathbf{w}_{\text{ext}} \): external wrench (e.g., payload, environment).
• \( \mathbf{w}_{\text{act}} \): wrench applied by the six actuators.

5.1 Mapping actuator forces to platform wrench

Each actuator applies a scalar force \( f_i \) along the direction of leg \(i\). Collect the forces into \(\mathbf{f} = [f_1,\dots,f_6]^{\top}\). The virtual work of the actuators for a virtual leg length variation \(\delta \mathbf{L}\) is

\[ \delta W_{\text{act}} = \mathbf{f}^{\top} \delta \mathbf{L}. \]

The corresponding virtual work on the platform (generalized virtual displacement \(\delta \mathbf{x}\)) is

\[ \delta W_{\text{plat}} = \mathbf{w}_{\text{act}}^{\top} \delta \mathbf{x}. \]

Ideal constraint transmission implies \( \delta W_{\text{act}} = \delta W_{\text{plat}} \) for all virtual displacements compatible with the constraints. Using the differential kinematics \(\delta \mathbf{L} = \mathbf{J}_L(\mathbf{x}) \delta \mathbf{x}\), we obtain

\[ \mathbf{f}^{\top} \delta \mathbf{L} = \mathbf{f}^{\top} \mathbf{J}_L(\mathbf{x}) \delta \mathbf{x} = \mathbf{w}_{\text{act}}^{\top} \delta \mathbf{x}. \]

Since this must hold for all \(\delta \mathbf{x}\), we conclude

\[ \mathbf{w}_{\text{act}} = \mathbf{J}_L^{\top}(\mathbf{x})\mathbf{f}. \]

The dynamic model of the platform driven by leg forces can therefore be written as

\[ \mathbf{M}(\mathbf{x})\dot{\mathbf{t}} + \mathbf{C}(\mathbf{x},\mathbf{t})\mathbf{t} + \mathbf{g}(\mathbf{x}) = \mathbf{w}_{\text{ext}} + \mathbf{J}_L^{\top}(\mathbf{x})\mathbf{f}. \]

This model directly connects actuator forces to platform accelerations. Conversely, a desired platform wrench can be transformed to actuator forces by inverting (or pseudo-inverting) \(\mathbf{J}_L^{\top}\) under appropriate conditions.

flowchart TD
  Fcmd["Actuator forces f (6x1)"] --> Wrench["Platform wrench w_act = J_L^T(x) * f"]
  Wrench --> Dyn["Dynamics: M(x)*t_dot + C(x,t)*t + g(x) = w_ext + w_act"]
  Dyn --> Int["Integrate t, x over time"]
  Int --> Update["Update geometry and J_L(x) for next step"]
        

For Python, we rely on numpy for linear algebra. For more advanced work, one can combine this with robotics libraries such as spatialmath-python or roboticstoolbox-python, but here we implement the core Stewart platform kinematics from scratch.

import numpy as np

class StewartPlatform:
    """
    Simple geometric model of a 6-6 Stewart platform.

    base_points:    array-like, shape (6, 3) with Bi in base frame
    platform_points: array-like, shape (6, 3) with Pi in platform frame
    """
    def __init__(self, base_points, platform_points):
        self.B = np.asarray(base_points, dtype=float).reshape(6, 3)
        self.P = np.asarray(platform_points, dtype=float).reshape(6, 3)

    @staticmethod
    def rot_zyx(phi, theta, psi):
        """Rotation matrix R = Rz(psi) * Ry(theta) * Rx(phi)."""
        cphi, sphi = np.cos(phi), np.sin(phi)
        cth, sth = np.cos(theta), np.sin(theta)
        cpsi, spsi = np.cos(psi), np.sin(psi)

        Rx = np.array([[1.0, 0.0, 0.0],
                       [0.0, cphi, -sphi],
                       [0.0, sphi, cphi]])
        Ry = np.array([[cth, 0.0, sth],
                       [0.0, 1.0, 0.0],
                       [-sth, 0.0, cth]])
        Rz = np.array([[cpsi, -spsi, 0.0],
                       [spsi, cpsi, 0.0],
                       [0.0, 0.0, 1.0]])
        return Rz @ Ry @ Rx

    def leg_vectors(self, p, R):
        """
        Compute leg vectors di, lengths Li, and unit vectors ui.

        Returns:
            d (6,3), L (6,), u (6,3)
        """
        p = np.asarray(p).reshape(3)
        R = np.asarray(R).reshape(3, 3)
        d = p + (R @ self.P.T).T - self.B  # di = p + R Pi - Bi
        L = np.linalg.norm(d, axis=1)
        u = d / L[:, None]
        return d, L, u

    def inverse_kinematics(self, pose):
        """
        IK: given pose x = [px, py, pz, phi, theta, psi],
        compute leg lengths L.
        """
        px, py, pz, phi, theta, psi = pose
        p = np.array([px, py, pz])
        R = self.rot_zyx(phi, theta, psi)
        _, L, _ = self.leg_vectors(p, R)
        return L

    def jacobian_length_rate(self, pose):
        """
        Compute J_L mapping platform twist t = [v; omega] to leg rates Ldot.
        """
        px, py, pz, phi, theta, psi = pose
        p = np.array([px, py, pz])
        R = self.rot_zyx(phi, theta, psi)
        d, L, u = self.leg_vectors(p, R)

        JL = np.zeros((6, 6))
        for i in range(6):
            ui = u[i, :]
            Rp = R @ self.P[i, :]
            # row: [ui^T, (Rp x ui)^T]
            JL[i, 0:3] = ui
            JL[i, 3:6] = np.cross(Rp, ui)
        return JL

    def forward_kinematics_newton(self, L_target, x0, max_iter=20, tol=1e-8):
        """
        Simple Newton-Raphson forward kinematics using constraint equations.
        pose x = [px, py, pz, phi, theta, psi].
        """
        x = np.array(x0, dtype=float).reshape(6)
        L_target = np.asarray(L_target, dtype=float).reshape(6)

        def constraint(xvec):
            px, py, pz, phi, theta, psi = xvec
            p = np.array([px, py, pz])
            R = self.rot_zyx(phi, theta, psi)
            _, L, _ = self.leg_vectors(p, R)
            # Phi_i = Li^2 - L_target_i^2
            return L**2 - L_target**2

        for k in range(max_iter):
            Phi = constraint(x)
            normPhi = np.linalg.norm(Phi)
            if normPhi < tol:
                break

            # Numerical Jacobian J_F via finite differences
            JF = np.zeros((6, 6))
            eps = 1e-6
            for j in range(6):
                dx = np.zeros(6)
                dx[j] = eps
                Phi_plus = constraint(x + dx)
                JF[:, j] = (Phi_plus - Phi) / eps

            # Newton step: x_{k+1} = x_k - JF^{-1} Phi
            delta = np.linalg.solve(JF, Phi)
            x = x - delta

        return x, normPhi
      

This implementation provides:

• A compact representation of the geometry.
• Closed-form inverse kinematics for leg lengths.
• The length-rate Jacobian \(\mathbf{J}_L\).
• A basic Newton–Raphson solver for forward kinematics (sensitive to initial guess and singularities).

In C++, we typically use the Eigen library for linear algebra and can interface with robotics libraries such as Orocos KDL and ROS. Below is a minimal Stewart platform geometry and Jacobian implementation using Eigen.

#include <Eigen/Dense>
#include <array>

class StewartPlatform {
public:
    using Vec3 = Eigen::Vector3d;
    using Mat3 = Eigen::Matrix3d;
    using Mat6 = Eigen::Matrix<double, 6, 6>;

    StewartPlatform(const std::array<Vec3, 6>& basePoints,
                    const std::array<Vec3, 6>& platformPoints)
        : B_(basePoints), P_(platformPoints) {}

    static Mat3 rotZYX(double phi, double theta, double psi) {
        double cphi = std::cos(phi), sphi = std::sin(phi);
        double cth  = std::cos(theta), sth  = std::sin(theta);
        double cpsi = std::cos(psi), spsi = std::sin(psi);

        Mat3 Rx, Ry, Rz;
        Rx << 1.0, 0.0, 0.0,
               0.0, cphi, -sphi,
               0.0, sphi, cphi;
        Ry << cth, 0.0, sth,
               0.0, 1.0, 0.0,
               -sth, 0.0, cth;
        Rz << cpsi, -spsi, 0.0,
               spsi, cpsi, 0.0,
               0.0, 0.0, 1.0;
        return Rz * Ry * Rx;
    }

    // IK: pose = [px, py, pz, phi, theta, psi]
    Eigen::Matrix<double, 6, 1> inverseKinematics(const Eigen::Matrix<double, 6, 1>& pose) const {
        double px = pose(0), py = pose(1), pz = pose(2);
        double phi = pose(3), theta = pose(4), psi = pose(5);
        Vec3 p(px, py, pz);
        Mat3 R = rotZYX(phi, theta, psi);

        Eigen::Matrix<double, 6, 1> L;
        for (int i = 0; i < 6; ++i) {
            Vec3 d = p + R * P_[i] - B_[i];
            L(i) = d.norm();
        }
        return L;
    }

    // Length-rate Jacobian JL
    Mat6 jacobianLengthRate(const Eigen::Matrix<double, 6, 1>& pose) const {
        double px = pose(0), py = pose(1), pz = pose(2);
        double phi = pose(3), theta = pose(4), psi = pose(5);
        Vec3 p(px, py, pz);
        Mat3 R = rotZYX(phi, theta, psi);

        Mat6 JL;
        JL.setZero();

        for (int i = 0; i < 6; ++i) {
            Vec3 d = p + R * P_[i] - B_[i];
            double L = d.norm();
            Vec3 u = d / L;
            Vec3 Rp = R * P_[i];
            Vec3 cross = Rp.cross(u);

            JL.block<1,3>(i, 0) = u.transpose();
            JL.block<1,3>(i, 3) = cross.transpose();
        }
        return JL;
    }

private:
    std::array<Vec3, 6> B_;  // base points
    std::array<Vec3, 6> P_;  // platform points
};
      

This C++ class can be integrated into a ROS node or used with Orocos KDL-based dynamics, where actuator forces and platform dynamics are handled by higher-level components.

In Java, one can use libraries such as Apache Commons Math or EJML for matrix computations. The following skeleton uses plain arrays for clarity and can be wrapped around a matrix library for efficiency and robustness.

public class StewartPlatform {
    private final double[][] B; // base points [6][3]
    private final double[][] P; // platform points [6][3]

    public StewartPlatform(double[][] basePoints, double[][] platformPoints) {
        if (basePoints.length != 6 || platformPoints.length != 6) {
            throw new IllegalArgumentException("Need 6 base and 6 platform points.");
        }
        this.B = basePoints;
        this.P = platformPoints;
    }

    public static double[][] rotZYX(double phi, double theta, double psi) {
        double cphi = Math.cos(phi), sphi = Math.sin(phi);
        double cth  = Math.cos(theta), sth  = Math.sin(theta);
        double cpsi = Math.cos(psi), spsi = Math.sin(psi);

        double[][] Rx = {
            {1.0, 0.0, 0.0},
            {0.0, cphi, -sphi},
            {0.0, sphi, cphi}
        };
        double[][] Ry = {
            {cth, 0.0, sth},
            {0.0, 1.0, 0.0},
            {-sth, 0.0, cth}
        };
        double[][] Rz = {
            {cpsi, -spsi, 0.0},
            {spsi, cpsi, 0.0},
            {0.0, 0.0, 1.0}
        };
        return matMul(Rz, matMul(Ry, Rx));
    }

    private static double[] matVec(double[][] M, double[] v) {
        double[] r = new double[3];
        for (int i = 0; i < 3; ++i) {
            r[i] = M[i][0] * v[0] + M[i][1] * v[1] + M[i][2] * v[2];
        }
        return r;
    }

    private static double[][] matMul(double[][] A, double[][] B) {
        double[][] R = new double[3][3];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                R[i][j] = 0.0;
                for (int k = 0; k < 3; ++k) {
                    R[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return R;
    }

    private static double[] vecSub(double[] a, double[] b) {
        return new double[]{a[0] - b[0], a[1] - b[1], a[2] - b[2]};
    }

    private static double norm(double[] v) {
        return Math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
    }

    // IK: pose = [px, py, pz, phi, theta, psi]
    public double[] inverseKinematics(double[] pose) {
        double px = pose[0], py = pose[1], pz = pose[2];
        double phi = pose[3], theta = pose[4], psi = pose[5];
        double[] p = new double[]{px, py, pz};
        double[][] R = rotZYX(phi, theta, psi);

        double[] L = new double[6];
        for (int i = 0; i < 6; ++i) {
            double[] Rp = matVec(R, P[i]);
            double[] di = vecSub(new double[]{p[0] + Rp[0], p[1] + Rp[1], p[2] + Rp[2]}, B[i]);
            L[i] = norm(di);
        }
        return L;
    }
}
      

This Java skeleton covers inverse kinematics and can be extended to compute the Jacobian and to interface with simulation frameworks or real-time controllers written in Java.

MATLAB is widely used in control engineering and robotics. The following script computes leg lengths and the Jacobian. In Simulink, these functions can be wrapped in MATLAB Function blocks and combined with a rigid-body dynamics block for simulation.

function stewart_demo()
    % Base and platform geometry (example values)
    B = [ 0.5, -0.5, 0.0;
          0.5,  0.5, 0.0;
          0.0,  0.6, 0.0;
         -0.5,  0.5, 0.0;
         -0.5, -0.5, 0.0;
          0.0, -0.6, 0.0 ];  % 6x3

    P = [ 0.3, -0.3, 0.0;
          0.3,  0.3, 0.0;
          0.0,  0.4, 0.0;
         -0.3,  0.3, 0.0;
         -0.3, -0.3, 0.0;
          0.0, -0.4, 0.0 ];  % 6x3

    % pose = [px, py, pz, phi, theta, psi]
    pose = [0.0, 0.0, 0.8, 0.0, 0.0, 0.0];

    L = stewart_inverse_kin(pose, B, P);
    JL = stewart_jacobian(pose, B, P);

    disp('Leg lengths:');
    disp(L.');
    disp('Length-rate Jacobian JL:');
    disp(JL);
end

function R = rotZYX(phi, theta, psi)
    cphi = cos(phi); sphi = sin(phi);
    cth  = cos(theta); sth  = sin(theta);
    cpsi = cos(psi); spsi = sin(psi);

    Rx = [1, 0, 0;
          0, cphi, -sphi;
          0, sphi, cphi];
    Ry = [cth, 0, sth;
          0, 1, 0;
          -sth, 0, cth];
    Rz = [cpsi, -spsi, 0;
          spsi,  cpsi, 0;
          0, 0, 1];
    R = Rz * Ry * Rx;
end

function L = stewart_inverse_kin(pose, B, P)
    px = pose(1); py = pose(2); pz = pose(3);
    phi = pose(4); theta = pose(5); psi = pose(6);

    p = [px; py; pz];
    R = rotZYX(phi, theta, psi);

    L = zeros(6, 1);
    for i = 1:6
        d = p + R * P(i, :).';
        d = d - B(i, :).';
        L(i) = norm(d);
    end
end

function JL = stewart_jacobian(pose, B, P)
    px = pose(1); py = pose(2); pz = pose(3);
    phi = pose(4); theta = pose(5); psi = pose(6);

    p = [px; py; pz];
    R = rotZYX(phi, theta, psi);

    JL = zeros(6, 6);
    for i = 1:6
        d = p + R * P(i, :).';
        d = d - B(i, :).';
        L = norm(d);
        u = d / L;
        Rp = R * P(i, :).';
        JL(i, 1:3) = u.';
        JL(i, 4:6) = cross(Rp, u).';
    end
end
      

In Simulink, one typical configuration is:

• A block computing IK for a desired platform pose, outputting leg lengths \(L_i\).
• A block integrating actuator dynamics and applying forces \(f_i\).
• A dynamics block computing the platform motion using \(\mathbf{M},\mathbf{C},\mathbf{g}\) and \(\mathbf{J}_L^{\top}\mathbf{f}\).

Wolfram Mathematica is very convenient for symbolic derivations and automatic code generation of the Stewart platform model. The snippet below sets up symbolic geometry and derives the constraint equations and their Jacobian.

(* Define symbolic variables *)
ClearAll["Global`*"];
(* Base and platform points (example) *)
B = Array[Subscript[b, #] &, {6}];
P = Array[Subscript[p, #] &, {6}];

(* Pose variables: position and ZYX Euler angles *)
px, py, pz, phi, theta, psi = Sequence @@ Array[x, 6];
p = {px, py, pz};

Rx[phi_] := { {1, 0, 0},
             {0, Cos[phi], -Sin[phi]},
             {0, Sin[phi], Cos[phi]} };

Ry[theta_] := { {Cos[theta], 0, Sin[theta]},
               {0, 1, 0},
               {-Sin[theta], 0, Cos[theta]} };

Rz[psi_] := { {Cos[psi], -Sin[psi], 0},
             {Sin[psi], Cos[psi], 0},
             {0, 0, 1} };

R[phi_, theta_, psi_] := Rz[psi].Ry[theta].Rx[phi];

(* Leg lengths and constraints *)
L = Array[Subscript[L, #] &, 6];
di[i_] := p + R[phi, theta, psi].P[[i]] - B[[i]];
len[i_] := Sqrt[di[i].di[i]];

Phi = Table[L[[i]]^2 - (di[i].di[i]), {i, 1, 6}];

(* Jacobian of constraints wrt pose x = {px,py,pz,phi,theta,psi} *)
xvars = {px, py, pz, phi, theta, psi};
JF = D[Phi, {xvars}];

(* Export symbolic Jacobian or generate C code *)
CForm[JF]
      

From this symbolic description, one can generate optimized C or C++ code for use in embedded controllers or real-time simulations of the Stewart platform.

In this case study, we modeled a Stewart platform as a parallel manipulator with six prismatic actuators. Starting from its geometry, we derived the constraint equations relating leg lengths and platform pose, emphasized the ease of inverse kinematics versus the nonlinear complexity of forward kinematics, and constructed the differential kinematics Jacobian mapping platform twist to leg length rates.

Using the virtual work principle, we related actuator forces to the platform wrench via the transpose of the Jacobian and embedded this mapping into a rigid-body dynamic model of the moving platform. Finally, we illustrated multi-language implementations (Python, C++, Java, MATLAB/Simulink, and Mathematica), demonstrating how the theoretical model can be transferred into practical software for simulation and control of Stewart platforms and other parallel robots.