Simulate Response of a Gimballed Pendulum with real life data

I have coded a gimballed pendulum model and I would like to modify the code so that I would be able to use real life data for the velocity(X_dot, Y_dot, Z_dot) and acceleration(X_ddot, Y_ddot, Z_ddot).

close all
clear
clc   %CLC clears the command window and homes the cursor

global  mt lt It tt mp lp Ip tp g H T k d t omega i

%System parameters .

mt=1.24; lt=0.28; It=0.1; tt=0.03;
mp=2.24; lp=0.15; Ip=0.1; tp=0.06;
g=9.81;

H=2; T=1; omega=2*pi/T; lambda=2.5; k=2*pi/lambda; t=[0:0.12:30]; x=0;
d=lambda/10;  

xs0=[0  0  0  0];   stsz=0.001;
i=1;
%Solver parameters setting


options=odeset('AbsTol',1e-6,'RelTol',1e-6);
[t x]=ode45(@Two_Dof_Air,t,[xs0(1) xs0(2) xs0(3) xs0(4)],options);
%ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb ,options

%%
function xs=Two_Dof_Air(t,x)
 
global  mt lt It tt mp lp Ip tp g H T k d t omega i X_dot Y_dot Z_dot X_ddot Y_ddot Z_ddot
 
H=2; T=1; omega=2*pi/T; lambda=2.5; k=2*pi/lambda; t=[0:0.5:50]; 
d=lambda/10; 
%state space vector
x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4);   

%     X_dot  = pi*H/T*cosh(k*(d))/sinh(k*d)*cos(-omega);              % velocity in x-axes
%     Y_dot  = pi*H/T*cosh(k*(d))/sinh(k*d)*cos(-omega);              % velocity in y-axes
%     Z_dot  = pi*H/T*sinh(k*(d))/sinh(k*d)*sin(-omega);              % velocity in z-axes
%     X_ddot = 2*(pi^2)*H/(T^2)*cosh(k*(d))/sinh(k*d)*cos(-omega);    % acceration in x-axes
%     Y_ddot = 2*(pi^2)*H/(T^2)*cosh(k*(d))/sinh(k*d)*cos(-omega);    % acceration in y-axes
%     Z_ddot = -2*(pi^2)*H/(T^2)*sinh(k*(d))/sinh(k*d)*sin(-omega);   % acceration in z-axes

%State space model equations
     xs1=x2;
     xs2=(((-mt*g*lt*sin(x1))/(It*cos(x3))) ...
         +(x4*sin(x3))*((2*x2/cos(x3))+(mt*lt*(omega*((cosh(k)*d)*(cos(-omega))/(sinh(k*d)))*cos(x1)+omega*((sinh(k)*d)*(sin(-omega))/(sinh(k*d)))*sin(x1)))/(It*(cos(x3).^2))) ...
         -((mt*lt*((-(omega^2)*H/(T^2))*((cosh(k)*d)*(sin(-omega))/(sinh(k*d)))*cos(x1)+(-(omega^2)*H/(T^2))*((sinh(k)*d)*(cos(-omega))/(sinh(k*d)))*sin(x1)))/(It*cos(x3))) ...
         +(-tt*sign(x2))/(It*(cos(x3).^2)));
     xs3=x4;
     xs4=(((-mp*g*lp*sin(x3))/(Ip*cos(x1))) ...
         +(x2*sin(x1))*((2*x4/cos(x1))+(mp*lp*(omega*((cosh(k)*d)*(cos(-omega))/(sinh(k*d)))*cos(x3)+omega*((sinh(k)*d)*(sin(-omega))/(sinh(k*d)))*sin(x3)))/(Ip*(cos(x1).^2))) ...
         -((mp*lp*((-(omega^2)*H/(T^2))*((cosh(k)*d)*(sin(-omega))/(sinh(k*d)))*cos(x3)+(-(omega^2)*H/(T^2))*((sinh(k)*d)*(cos(-omega))/(sinh(k*d)))*sin(x3)))/(Ip*cos(x1))) ...
         +(-tp*sign(x4))/(Ip*(cos(x1).^2)));
% The original equation for the pendulum     
%    xs2=(((-mt*g*lt*sin(x1))/(It*cos(x3))) ...                       
%        +(x4*sin(x3))*((2*x2/cos(x3))+(mt*lt/(It*(cos(x3).^2))*(X_dot*cos(x1)+Z_dot*sin(x1)))) ...
%        -((mt*lt*(X_ddot*cos(x1)+Z_ddot*sin(x1)))/(It*cos(x3))) ...
%        +(-tt*sign(x2))/(It*(cos(x3).^2)));
%    xs4=(((-mp*g*lp*sin(x3))/(Ip*cos(x1)))+ ...                      
%        (x2*sin(x1))*((2*x4/cos(x1))+(mp*lp*(Y_dot*cos(x3)+Z_dot*sin(x3)))/(Ip*(cos(x1).^2))) ...
%        -((mp*lp*(Y_ddot*cos(x3)+Z_ddot*sin(x3)))/(Ip*cos(x1))) ...
%        +(-tp*sign(x4))/(Ip*(cos(x1).^2)));     
   xs=[xs1; xs2;  xs3;  xs4];
    
 end

Answer 1

[Abhimenyu]

Hi Jun,

I understand that you want to modify your code so that it can use real-life data for velocity and acceleration. This can be done by loading the data from a file inside your function.

Please follow the below mentioned steps to use the real-life data into your code:

• Load the real-life data: Please load the real-life data into your MATLAB environment using the "load" function if the data is saved in a .mat file. This data could be in the form of arrays that correspond to the "t" array. Please follow the example MATLAB code below:

close all 
clear 
clc 

% Load your data here or pass it into the function 
% For example: 
% data = load('datafile.mat'); 
% X_dot_data = data.X_dot; 
% Y_dot_data = data.Y_dot; 
% Z_dot_data = data.Z_dot; 
% X_ddot_data = data.X_ddot; 
% Y_ddot_data = data.Y_ddot; 
% Z_ddot_data = data.Z_ddot; 
% time_data = data.time; 

%define the global variables along with other variables 
global mt lt It tt mp lp Ip tp g H T k d t omega i 
global X_dot_data Y_dot_data Z_dot_data X_ddot_data Y_ddot_data Z_ddot_data time_data 

• Interpolate the data and modify the "Two_Dof_Air" function: Since the ODE solver may use time points that are not exactly the same as the ones in the real-life data, the data should be interpolated using the "interp1" MATLAB function to find the values of velocity and acceleration at the exact time points needed by the solver. Please follow the example MATLAB code below:

function xs = Two_Dof_Air(t, x) 

    global mt lt It tt mp lp Ip tp g H T k d t omega i 
    global X_dot_data Y_dot_data Z_dot_data X_ddot_data Y_ddot_data Z_ddot_data time_data 
    
    % Interpolate the data to get the values at the current time t 
    
    %other interpolation options can be included as well 
    
    X_dot = interp1(time_data, X_dot_data, t, 'linear', 'extrap'); 
    Y_dot = interp1(time_data, Y_dot_data, t, 'linear', 'extrap'); 
    Z_dot = interp1(time_data, Z_dot_data, t, 'linear', 'extrap'); 
    X_ddot = interp1(time_data, X_ddot_data, t, 'linear', 'extrap'); 
    Y_ddot = interp1(time_data, Y_ddot_data, t, 'linear', 'extrap'); 
    Z_ddot = interp1(time_data, Z_ddot_data, t, 'linear', 'extrap'); 
    
    % state space vector 
    x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); 
    
    %rest of the code remains same 
    %Use the real data (X_dot, X_ddot etc) instead of the theoretical 
    %expressions to compute xs1, xs2 etc. 
    xs1 = x2; 
    xs2 = (((-mt*g*lt*sin(x1))/(It*cos(x3))) ... 
           +(x4*sin(x3))*((2*x2/cos(x3))+(mt*lt/(It*(cos(x3).^2))*(X_dot*cos(x1)+Z_dot*sin(x1)))) ... 
           -((mt*lt*(X_ddot*cos(x1)+Z_ddot*sin(x1)))/(It*cos(x3))) ... 
           +(-tt*sign(x2))/(It*(cos(x3).^2)));
    
    %compute values for xs3 and xs4 similarly
    xs=[xs1; xs2;  xs3;  xs4];
    
 end

The above-mentioned steps will help to use real-life values for velocity and acceleration in your code.

Please follow the below mentioned MATLAB R2023b documentation links to get more information on "load" and "interp1" functions respectively:

• https://www.mathworks.com/help/matlab/ref/load.html
• https://www.mathworks.com/help/matlab/ref/interp1.html

I hope this helps to resolve your query.

Regards,

Abhimenyu