Identification of the Stewart Platform using Simscape

Identification
States as the motion of the mobile platform
Simple Model without any sensor
Cartesian Plot
From a force to force sensor
From a force applied in the leg to the displacement of the leg
Transmissibility
Compliance
Inertial

We would like to extract a state space model of the Stewart Platform from the Simscape model.

The inputs are:

Symbol	Meaning
\(\bm{\mathcal{F}}_{d}\)	External forces applied in {B}
\(\bm{\tau}\)	Joint forces
\(\bm{\mathcal{F}}\)	Cartesian forces applied by the Joints
\(\bm{D}_{w}\)	Fixed Based translation and rotations around {A}

The outputs are:

Symbol	Meaning
\(\bm{\mathcal{X}}\)	Relative Motion of {B} with respect to {A}
\(\bm{\mathcal{L}}\)	Joint Displacement
\(\bm{F}_{m}\)	Force Sensors in each strut
\(\bm{v}_{m}\)	Inertial Sensors located at \(b_i\) measuring in the direction of the strut

An important difference from basic Simulink models is that the states in a physical network are not independent in general, because some states have dependencies on other states through constraints.

Identification

Simscape Model

Initialize the Stewart Platform

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = computeJacobian(stewart);
stewart = initializeStewartPose(stewart);

Identification

%% Options for Linearized
options = linearizeOptions;
options.SampleTime = 0;

%% Name of the Simulink File
mdl = 'stewart_platform_identification';

%% Input/Output definition
clear io; io_i = 1;
io(io_i) = linio([mdl, '/tau'],  1, 'openinput');  io_i = io_i + 1;
io(io_i) = linio([mdl, '/Fext'], 1, 'openinput');  io_i = io_i + 1;
io(io_i) = linio([mdl, '/X'],    1, 'openoutput'); io_i = io_i + 1;
io(io_i) = linio([mdl, '/Vm'],   1, 'openoutput'); io_i = io_i + 1;
io(io_i) = linio([mdl, '/Taum'], 1, 'openoutput'); io_i = io_i + 1;
io(io_i) = linio([mdl, '/Lm'],   1, 'openoutput'); io_i = io_i + 1;

%% Run the linearization
G = linearize(mdl, io, options);
G.InputName  = {'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6', ...
                'Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'};

G.OutputName = {'Xdx', 'Xdy', 'Xdz', 'Xrx', 'Xry', 'Xrz', ...
                'Vm1', 'Vm2', 'Vm3', 'Vm4', 'Vm5', 'Vm6', ...
                'taum1', 'taum2', 'taum3', 'taum4', 'taum5', 'taum6', ...
                'Lm1', 'Lm2', 'Lm3', 'Lm4', 'Lm5', 'Lm6'};

States as the motion of the mobile platform

Initialize the Stewart Platform

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = computeJacobian(stewart);
stewart = initializeStewartPose(stewart);

Identification

%% Options for Linearized
options = linearizeOptions;
options.SampleTime = 0;

%% Name of the Simulink File
mdl = 'stewart_platform_identification_simple';

%% Input/Output definition
clear io; io_i = 1;
io(io_i) = linio([mdl, '/tau'],  1, 'openinput');  io_i = io_i + 1;
io(io_i) = linio([mdl, '/X'],    1, 'openoutput'); io_i = io_i + 1;
io(io_i) = linio([mdl, '/Xdot'], 1, 'openoutput'); io_i = io_i + 1;

%% Run the linearization
G = linearize(mdl, io);
% G.InputName  = {'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'};
% G.OutputName = {'Xdx', 'Xdy', 'Xdz', 'Xrx', 'Xry', 'Xrz', 'Vdx', 'Vdy', 'Vdz', 'Vrx', 'Vry', 'Vrz'};

Let’s check the size of G:

size(G)

size(G)
State-space model with 12 outputs, 6 inputs, and 18 states.
'org_babel_eoe'
ans =
    'org_babel_eoe'

We expect to have only 12 states (corresponding to the 6dof of the mobile platform).

Gm = minreal(G);

Gm = minreal(G);
6 states removed.

And indeed, we obtain 12 states.

Coordinate transformation

We can perform the following transformation using the ss2ss command.

Gt = ss2ss(Gm, Gm.C);

Then, the C matrix of Gt is the unity matrix which means that the states of the state space model are equal to the measurements \(\bm{Y}\).

The measurements are the 6 displacement and 6 velocities of mobile platform with respect to \(\{B\}\).

We could perform the transformation by hand:

At = Gm.C*Gm.A*pinv(Gm.C);

Bt = Gm.C*Gm.B;

Ct = eye(12);
Dt = zeros(12, 6);

Gt = ss(At, Bt, Ct, Dt);

Analysis

[V,D] = eig(Gt.A);

Mode Number	Resonance Frequency [Hz]	Damping Ratio [%]
1.0	174.5	0.9
2.0	174.5	0.7
3.0	202.1	0.7
4.0	237.3	0.6
5.0	237.3	0.5
6.0	283.8	0.5

Visualizing the modes

To visualize the i’th mode, we may excite the system using the inputs \(U_i\) such that \(B U_i\) is co-linear to \(\xi_i\) (the mode we want to excite).

\[ U(t) = e^{\alpha t} ( ) \]

Let’s first sort the modes and just take the modes corresponding to a eigenvalue with a positive imaginary part.

ws = imag(diag(D));
[ws,I] = sort(ws)
ws = ws(7:end); I = I(7:end);

for i = 1:length(ws)

i_mode = I(i); % the argument is the i'th mode

lambda_i = D(i_mode, i_mode);
xi_i = V(:,i_mode);

a_i = real(lambda_i);
b_i = imag(lambda_i);

Let do 10 periods of the mode.

t = linspace(0, 10/(imag(lambda_i)/2/pi), 1000);
U_i = pinv(Gt.B) * real(xi_i * lambda_i * (cos(b_i * t) + 1i*sin(b_i * t)));

U = timeseries(U_i, t);

Simulation:

load('mat/conf_simscape.mat');
set_param(conf_simscape, 'StopTime', num2str(t(end)));
sim(mdl);

Save the movie of the mode shape.

smwritevideo(mdl, sprintf('figs/mode%i', i), ...
             'PlaybackSpeedRatio', 1/(b_i/2/pi), ...
             'FrameRate', 30, ...
             'FrameSize', [800, 400]);

end

Figure 1: Identified mode - 1

Figure 2: Identified mode - 3

Figure 3: Identified mode - 5

Identification

%% Options for Linearized
options = linearizeOptions;
options.SampleTime = 0;

%% Name of the Simulink File
mdl = 'stewart_platform_identification';

%% Input/Output definition
clear io; io_i = 1;
io(io_i) = linio([mdl, '/tau'],  1, 'openinput');  io_i = io_i + 1;
io(io_i) = linio([mdl, '/Lm'],    1, 'openoutput'); io_i = io_i + 1;

%% Run the linearization
G = linearize(mdl, io, options);
% G.InputName  = {'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'};
% G.OutputName = {'Xdx', 'Xdy', 'Xdz', 'Xrx', 'Xry', 'Xrz', 'Vdx', 'Vdy', 'Vdz', 'Vrx', 'Vry', 'Vrz'};

size(G)

Change of states

At = G.C*G.A*pinv(G.C);

Bt = G.C*G.B;

Ct = eye(12);
Dt = zeros(12, 6);

Gt = ss(At, Bt, Ct, Dt);

size(Gt)

Simple Model without any sensor

Simscape Model

open 'stewart_identification_simple.slx'

Initialize the Stewart Platform

stewart = initializeStewartPlatform();
stewart = initializeFramesPositions(stewart);
stewart = generateGeneralConfiguration(stewart);
stewart = computeJointsPose(stewart);
stewart = initializeStrutDynamics(stewart);
stewart = initializeCylindricalPlatforms(stewart);
stewart = initializeCylindricalStruts(stewart);
stewart = computeJacobian(stewart);
stewart = initializeStewartPose(stewart);

Identification

stateorder = {...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_1_1_1',...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_2_1_1',...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_3_1_1',...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_4_1_1',...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_5_1_1',...
    'stewart_platform_identification_simple/Solver Configuration/EVAL_KEY/INPUT_6_1_1',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.cylindrical_joint.Rz.q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.cylindrical_joint.Pz.p',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.cylindrical_joint.Rz.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.cylindrical_joint.Pz.v',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.spherical_joint_F.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_2.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_3.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_4.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_5.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_6.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.spherical_joint_F.S.w',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.spherical_joint_M.S.Q',...
    'stewart_platform_identification_simple.Stewart_Platform.Strut_1.Subsystem.spherical_joint_M.S.w'};

%% Options for Linearized
options = linearizeOptions;
options.SampleTime = 0;

%% Name of the Simulink File
mdl = 'stewart_platform_identification_simple';

%% Input/Output definition
clear io; io_i = 1;
io(io_i) = linio([mdl, '/tau'],  1, 'openinput');  io_i = io_i + 1;
io(io_i) = linio([mdl, '/X'],     1, 'openoutput'); io_i = io_i + 1;
io(io_i) = linio([mdl, '/Xdot'],  1, 'openoutput'); io_i = io_i + 1;

%% Run the linearization
G = linearize(mdl, io, options);
G.InputName  = {'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'};

G.OutputName = {'Xdx', 'Xdy', 'Xdz', 'Xrx', 'Xry', 'Xrz', 'Vdx', 'Vdy', 'Vdz', 'Vrx', 'Vry', 'Vrz'};

size(G)

G.StateName

Cartesian Plot

From a force applied in the Cartesian frame to a displacement in the Cartesian frame.

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_cart(1, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_cart(2, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_cart(3, 1), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude');

figure;
bode(G.G_cart, freqs);

From a force to force sensor

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_forc(1, 1), freqs, 'Hz'))), 'k-', 'DisplayName', '$F_{m_i}/F_{i}$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [N/N]');
legend('location', 'southeast');

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_forc(1, 1), freqs, 'Hz'))), 'k-', 'DisplayName', '$F_{m_i}/F_{i}$');
plot(freqs, abs(squeeze(freqresp(G.G_forc(2, 1), freqs, 'Hz'))), 'k--', 'DisplayName', '$F_{m_j}/F_{i}$');
plot(freqs, abs(squeeze(freqresp(G.G_forc(3, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_forc(4, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_forc(5, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_forc(6, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [N/N]');
legend('location', 'southeast');

From a force applied in the leg to the displacement of the leg

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_legs(1, 1), freqs, 'Hz'))), 'k-', 'DisplayName', '$D_{i}/F_{i}$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_legs(1, 1), freqs, 'Hz'))), 'k-', 'DisplayName', '$D_{i}/F_{i}$');
plot(freqs, abs(squeeze(freqresp(G.G_legs(2, 1), freqs, 'Hz'))), 'k--', 'DisplayName', '$D_{j}/F_{i}$');
plot(freqs, abs(squeeze(freqresp(G.G_legs(3, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_legs(4, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_legs(5, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G.G_legs(6, 1), freqs, 'Hz'))), 'k--', 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
legend('location', 'northeast');

Transmissibility

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_tran(1, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(2, 2), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(3, 3), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/m]');

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_tran(4, 4), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(5, 5), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(6, 6), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [$\frac{rad/s}{rad/s}$]');

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_tran(1, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(1, 2), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_tran(1, 3), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/m]');

Compliance

From a force applied in the Cartesian frame to a relative displacement of the mobile platform with respect to the base.

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_comp(1, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_comp(2, 2), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_comp(3, 3), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');

Inertial

From a force applied on the Cartesian frame to the absolute displacement of the mobile platform.

figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G.G_iner(1, 1), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_iner(2, 2), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G.G_iner(3, 3), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');

Identification of the Stewart Platform using Simscape

Table of Contents

Identification

Simscape Model

Initialize the Stewart Platform

Identification

States as the motion of the mobile platform

Initialize the Stewart Platform

Identification

Coordinate transformation

Analysis

Visualizing the modes

Identification

Change of states

Simple Model without any sensor

Simscape Model

Initialize the Stewart Platform

Identification

Cartesian Plot

From a force to force sensor

From a force applied in the leg to the displacement of the leg

Transmissibility

Compliance

Inertial