Simulate Response of a Gimballed Pendulum with real life data

2 views (last 30 days)
Jun Da Wee on 2 Mar 2022
Answered: Abhimenyu on 31 Jan 2024
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

Abhimenyu on 31 Jan 2024
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.
• 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
% For example:
% 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.
I hope this helps to resolve your query.
Regards,
Abhimenyu

Categories

Find more on Operating on Diagonal Matrices in Help Center and File Exchange

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!