Kinematic Study of the Stewart Platform

If we note:

• \({}^A\hat{\bm{s}}_i\) the unit vector representing the direction of the i’th strut and expressed in frame \(\{A\}\)
• \({}^A\bm{b}_i\) the position vector of the i’th joint fixed to the mobile platform and expressed in frame \(\{A\}\)

Then, we can compute the Jacobian with the following equation (the superscript \(A\) is ignored):

The Jacobian matrix \(\bm{J}\) can be computed using the computeJacobian function (accessible here). For instance:

  stewart = computeJacobian(stewart);

This will add three new matrix to the stewart structure:

• J the Jacobian matrix
• K the stiffness matrix
• C the compliance matrix

The Jacobian matrix links the input joint rate \(\dot{\bm{\mathcal{L}}} = [ \dot{l}_1, \dot{l}_2, \dot{l}_3, \dot{l}_4, \dot{l}_5, \dot{l}_6 ]^T\) of each strut to the output twist vector of the mobile platform is denoted by \(\dot{\bm{X}} = [^A\bm{v}_p, {}^A\bm{\omega}]^T\):

The input joint rate \(\dot{\bm{\mathcal{L}}}\) can be measured by taking the derivative of the relative motion sensor in each strut. The output twist vector can be measured with a “Transform Sensor” block measuring the relative velocity and relative angular velocity of frame \(\{B\}\) with respect to frame \(\{A\}\).

If the Jacobian matrix is inversible, we can also compute \(\dot{\bm{\mathcal{X}}}\) from \(\dot{\bm{\mathcal{L}}}\).

The Jacobian matrix can also be used to approximate forward and inverse kinematics for small displacements. This is explained in section 3.3.

If we note:

• \(\bm{\tau} = [\tau_1, \tau_2, \cdots, \tau_6]^T\): vector of actuator forces applied in each strut
• \(\bm{\mathcal{F}} = [\bm{f}, \bm{n}]^T\): external force/torque action on the mobile platform at \(\bm{O}_B\)

We find that the transpose of the Jacobian matrix links the two by the following equation:

If the Jacobian matrix is inversible, we also have the following relation:

As explain in this document, each Actuator is modeled by 3 elements in parallel:

• A spring with a stiffness \(k_{i}\)
• A dashpot with a damping \(c_{i}\)
• A force source \(\tau_i\)

The stiffness of the actuator \(k_i\) links the applied actuator force \(\delta \tau_i\) and the corresponding small deflection \(\delta l_i\):

If we combine these 6 relations:

Substituting \(\bm{\tau} = \bm{J}^{-T} \bm{\mathcal{F}}\) and \(\delta \bm{\mathcal{L}} = \bm{J} \cdot \delta \bm{\mathcal{X}}\) gives

And then we identify the stiffness matrix \(\bm{K}\):

If the stiffness matrix \(\bm{K}\) is inversible, the compliance matrix of the manipulator is defined as

The compliance matrix of a manipulator shows the mapping of the moving platform wrench applied at \(\bm{O}_B\) to its small deflection by

The stiffness and compliance matrices are computed using the computeJacobian function (accessible here).

For inverse kinematic analysis, it is assumed that the position \({}^A\bm{P}\) and orientation of the moving platform \({}^A\bm{R}_B\) are given and the problem is to obtain the joint variables \(\bm{\mathcal{L}} = \left[ l_1, l_2, l_3, l_4, l_5, l_6 \right]^T\).

This problem can be easily solved using the loop closures.

The obtain joint variables are:

If the position and orientation of the platform lie in the feasible workspace, the solution is unique. Otherwise, the solution gives complex numbers.

This inverse kinematic solution can be obtained using the function inverseKinematics (described here).

In forward kinematic analysis, it is assumed that the vector of limb lengths \(\bm{L}\) is given and the problem is to find the position \({}^A\bm{P}\) and the orientation \({}^A\bm{R}_B\).

This is a difficult problem that requires to solve nonlinear equations.

In a next section, an approximate solution of the forward kinematics problem is proposed for small displacements.

For small displacements mobile platform displacement \(\delta \bm{\mathcal{X}} = [\delta x, \delta y, \delta z, \delta \theta_x, \delta \theta_y, \delta \theta_z ]^T\) around \(\bm{\mathcal{X}}_0\), the associated joint displacement can be computed using the Jacobian (approximate solution of the inverse kinematic problem):

Similarly, for small joint displacements \(\delta\bm{\mathcal{L}} = [ \delta l_1,\ \dots,\ \delta l_6 ]^T\) around \(\bm{\mathcal{L}}_0\), it is possible to find the induced small displacement of the mobile platform (approximate solution of the forward kinematic problem):

These two relations solve the forward and inverse kinematic problems for small displacement in a approximate way. As the inverse kinematic can be easily solved exactly this is not much useful, however, as the forward kinematic problem is difficult to solve, this approximation can be very useful for small displacements.

The function forwardKinematicsApprox (described here) can be used to solve the forward kinematic problem using the Jacobian matrix.

We first define some general Stewart architecture.

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart, 'H', 90e-3, 'MO_B', 45e-3);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStewartPose(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeJointDynamics(stewart);
stewart = computeJacobian(stewart);

Let’s first compare the perfect and approximate solution of the inverse for pure \(x\) translations.

We compute the approximate and exact required strut stroke to have the wanted mobile platform \(x\) displacement. The estimate required strut stroke for both the approximate and exact solutions are shown in Figure 1. The relative strut length displacement is shown in Figure 2.

Xrs = logspace(-6, -1, 100); % Wanted X translation of the mobile platform [m]

Ls_approx = zeros(6, length(Xrs));
Ls_exact = zeros(6, length(Xrs));

for i = 1:length(Xrs)
    Xr = Xrs(i);
    L_approx(:, i) = stewart.kinematics.J*[Xr; 0; 0; 0; 0; 0;];
    [~, L_exact(:, i)] = inverseKinematics(stewart, 'AP', [Xr; 0; 0]);
end

Let’s first define the Stewart platform architecture that we want to study.

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart, 'H', 90e-3, 'MO_B', 45e-3);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStewartPose(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeJointDynamics(stewart);
stewart = computeJacobian(stewart);

Let’s now define the wanted extreme translations and rotations.

Tx_max = 50e-6; % Translation [m]
Ty_max = 50e-6; % Translation [m]
Tz_max = 50e-6; % Translation [m]
Rx_max = 30e-6; % Rotation [rad]
Ry_max = 30e-6; % Rotation [rad]
Rz_max = 0;     % Rotation [rad]

As a first estimation, we estimate the needed actuator stroke for “pure” rotations and translation. We do that using either the Inverse Kinematic solution or the Jacobian matrix as an approximation.

LTx = stewart.kinematics.J*[Tx_max 0 0 0 0 0]';
LTy = stewart.kinematics.J*[0 Ty_max 0 0 0 0]';
LTz = stewart.kinematics.J*[0 0 Tz_max 0 0 0]';
LRx = stewart.kinematics.J*[0 0 0 Rx_max 0 0]';
LRy = stewart.kinematics.J*[0 0 0 0 Ry_max 0]';
LRz = stewart.kinematics.J*[0 0 0 0 0 Rz_max]';

The obtain required stroke is:

From -3.8e-05[m] to 3.8e-05[m]: Total stroke = 76.1[um]

This is surely a low estimation of the required stroke.

We know would like to have a more precise estimation.

To do so, we may estimate the required actuator stroke for all possible combination of translation and rotation.

Let’s first generate all the possible combination of maximum translation and rotations.

Ps = [2*(dec2bin(0:5^2-1,5)-'0')-1, zeros(5^2, 1)].*[Tx_max Ty_max Tz_max Rx_max Ry_max Rz_max];

For all possible combination, we compute the required actuator stroke using the inverse kinematic solution.

L_min = 0;
L_max = 0;

for i = 1:size(Ps,1)
    Rx = [1 0        0;
          0 cos(Ps(i, 4)) -sin(Ps(i, 4));
          0 sin(Ps(i, 4))  cos(Ps(i, 4))];

    Ry = [ cos(Ps(i, 5)) 0 sin(Ps(i, 5));
           0        1 0;
           -sin(Ps(i, 5)) 0 cos(Ps(i, 5))];

    Rz = [cos(Ps(i, 6)) -sin(Ps(i, 6)) 0;
          sin(Ps(i, 6))  cos(Ps(i, 6)) 0;
          0        0       1];

    ARB = Rz*Ry*Rx;
    [~, Ls] = inverseKinematics(stewart, 'AP', Ps(i, 1:3)', 'ARB', ARB);

    if min(Ls) < L_min
        L_min = min(Ls)
    end
    if max(Ls) > L_max
        L_max = max(Ls)
    end
end

We obtain the required actuator stroke:

From -8.9e-05[m] to 8.9e-05[m]: Total stroke = 177.2[um]

This is probably a much realistic estimation of the required actuator stroke.

Let’s first define the Stewart platform architecture that we want to study.

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart, 'H', 90e-3, 'MO_B', 45e-3);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStewartPose(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeJointDynamics(stewart);
stewart = computeJacobian(stewart);

Let’s now define the actuator stroke.

L_min = -50e-6; % [m]
L_max =  50e-6; % [m]

Let’s first estimate the mobility in translation when the orientation of the Stewart platform stays the same.

As shown previously, for such small stroke, we can use the approximate Forward Dynamics solution using the Jacobian matrix:

To obtain the mobility “volume” attainable by the Stewart platform when it’s orientation is set to zero, we use the spherical coordinate \((r, \theta, \phi)\).

For each possible value of \((\theta, \phi)\), we compute the maximum radius \(r\) attainable with the constraint that the stroke of each actuator should be between L_min and L_max.

thetas = linspace(0, pi, 50);
phis = linspace(0, 2*pi, 50);
rs = zeros(length(thetas), length(phis));

for i = 1:length(thetas)
    for j = 1:length(phis)
        Tx = sin(thetas(i))*cos(phis(j));
        Ty = sin(thetas(i))*sin(phis(j));
        Tz = cos(thetas(i));

        dL = stewart.kinematics.J*[Tx; Ty; Tz; 0; 0; 0;]; % dL required for 1m displacement in theta/phi direction

        rs(i, j) = max([dL(dL<0)*L_min; dL(dL>0)*L_max]);
    end
end

Now that we have found the corresponding radius \(r\), we plot the obtained mobility. We can also approximate the mobility by a sphere with a radius equal to the minimum obtained value of \(r\), this is however a pessimistic estimation of the mobility.

Let’s first define the Stewart Platform Geometry.

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart, 'H', 90e-3, 'MO_B', 150e-3);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart, 'AP', [0; 0; 0]);

As_init = stewart.geometry.As;

Tx_max = 60e-6; % Translation [m]
Ty_max = 60e-6; % Translation [m]
Tz_max = 60e-6; % Translation [m]
Rx_max = 30e-6; % Rotation [rad]
Ry_max = 30e-6; % Rotation [rad]
Rz_max = 0; % Rotation [rad]

Ps = [2*(dec2bin(0:5^2-1,5)-'0')-1, zeros(5^2, 1)].*[Tx_max Ty_max Tz_max Rx_max Ry_max Rz_max];

flex_ang = zeros(size(Ps, 1), 6);
Rxs = zeros(size(Ps, 1), 6)
Rys = zeros(size(Ps, 1), 6)
Rzs = zeros(size(Ps, 1), 6)

for Ps_i = 1:size(Ps, 1)
    Rx = [1 0        0;
          0 cos(Ps(Ps_i, 4)) -sin(Ps(Ps_i, 4));
          0 sin(Ps(Ps_i, 4))  cos(Ps(Ps_i, 4))];

    Ry = [ cos(Ps(Ps_i, 5)) 0 sin(Ps(Ps_i, 5));
           0        1 0;
           -sin(Ps(Ps_i, 5)) 0 cos(Ps(Ps_i, 5))];

    Rz = [cos(Ps(Ps_i, 6)) -sin(Ps(Ps_i, 6)) 0;
          sin(Ps(Ps_i, 6))  cos(Ps(Ps_i, 6)) 0;
          0        0       1];

    ARB = Rz*Ry*Rx;

    stewart = computeJointsPose(stewart, 'AP', Ps(Ps_i, 1:3)', 'ARB', ARB);

    flex_ang(Ps_i, :) = acos(sum(As_init.*stewart.geometry.As));

    for l_i = 1:6
        MRf = stewart.platform_M.MRb(:,:,l_i)*(ARB')*(stewart.platform_F.FRa(:,:,l_i)');

        Rys(Ps_i, l_i) = atan2(MRf(1,3), sqrt(MRf(1,1)^2 + MRf(2,2)^2));
        Rxs(Ps_i, l_i) = atan2(-MRf(2,3)/cos(Rys(Ps_i, l_i)), MRf(3,3)/cos(Rys(Ps_i, l_i)));
        Rzs(Ps_i, l_i) = atan2(-MRf(1,2)/cos(Rys(Ps_i, l_i)), MRf(1,1)/cos(Rys(Ps_i, l_i)));
    end
end

And the maximum bending of the flexible joints is: (in [mrad])

1e3*max(max(abs(flex_ang)))

1.1704

1e3*max(max(sqrt(Rxs.^2 + Rys.^2)))

0.075063

1e3*max(max(Rzs))

0.04666

This Matlab function is accessible here.

This Matlab function is accessible here.

This Matlab function is accessible here.