Numerical integration of the missile dynamic model

Hello everyone,

I want to integrate the dynamic model of a 6DOF flying missile (using ode45 and RK4 methods) and plot the trajectory of the missile , this model contains 12 nonlinear ode's and discrete data, external forces and moments are calcualted by a seperate function and when I launch the simulation I get this error :

Warning: Failure at t=5.051155e+00. Unable to meet integration tolerances without reducing the step size below the smallest

value allowed (1.421085e-14) at time t.

> In ode45 (line 360)

In ode45_integration (line 11) .

I would appreciate your help and your suggestions. Best regards

My functions are :

%%%% Dynamic model integration using ode45 %%%%

t0 = 0;

tf = 15;

h = 0.01;

timerange = t0:h:tf;

IC = [0;0;1;0;deg2rad(-18);0;13;0;0;0;0;0];

%% Integration par ode45 %%%

[t,Y] = ode45(@(t,Y) MDD_Missile(t,Y),timerange,IC);

plot(Y(:,1),Y(:,3))

xlabel('range');

ylabel('Altitude(m)');

legend('altitude en fonction de la porteé');

grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dYdt = MDD_Missile(t,Y)

%%% Missile discrete data %%%

time = [0; 0.3; 0.6; 1.2; 1.8; 2.4; 4.2; 5.2];

m = [11.25; 11.16; 11.06; 10.82; 10.58; 10.38; 10.16; 10.15];

TR = [0; 570; 650; 750; 770; 650; 50; 0];

rxA = [0.565; 0.555; 0.544; 0.519; 0.493; 0.471; 0.447; 0.446];

Ixx = [0.0252; 0.025; 0.0248; 0.0244; 0.0239; 0.0235; 0.0231; 0.0231];

Iyy = [0.985; 0.979; 0.973; 0.958; 0.942; 0.929; 0.915; 0.914];

Izz = [0.985; 0.979; 0.973; 0.958; 0.942; 0.929; 0.915; 0.914];

%%% Interpolation data-set with t %%%

m = interp1(time,m,t);

Ixx = interp1(time,Ixx,t);

Iyy = interp1(time,Iyy,t);

Izz = interp1(time,Izz,t);

TR = interp1(time,TR,t);

rxA = interp1(time,rxA,t);

IG = [Ixx 0 0;

0 Iyy 0;

0 0 Izz];

%% Command parameters (Fins orientation angles) %%

sigma_9g = deg2rad (0);

sigma_10g = deg2rad (0);

sigma_11g = deg2rad (0);

sigma_12g = deg2rad (0);

%% Call of fonction calculting external forces and moments acting on the missile %%

ForcesMoments_Aero = F_M (Y,m,TR,rxA,sigma_9g,sigma_10g,sigma_11g,sigma_12g);

Fx = ForcesMoments_Aero(1);

Fy = ForcesMoments_Aero(2);

Fz = ForcesMoments_Aero(3);

L = ForcesMoments_Aero(4);

M = ForcesMoments_Aero(5);

N = ForcesMoments_Aero(6);

%%% 12 différential equation of motion of the missile %%%

VF = [Fx; Fy; Fz];

VM = [L; M; N];

Mat_oRb = [cos(Y(5))*cos(Y(6)) cos(Y(6))*sin(Y(5))*sin(Y(4))-sin(Y(6))*cos(Y(4)) cos(Y(6))*sin(Y(5))*cos(Y(4))+sin(Y(6))*sin(Y(4));

sin(Y(6))*cos(Y(5)) sin(Y(6))*sin(Y(5))*sin(Y(4))+cos(Y(6))*cos(Y(5)) sin(Y(6))*sin(Y(5))*cos(Y(4))-cos(Y(6))*sin(Y(4));

-sin(Y(5)) cos(Y(5))*sin(Y(4)) cos(Y(5))*cos(Y(4))];

Mat_H = [1 sin(Y(4))*tan(Y(5)) cos(Y(4))*tan(Y(5));

0 cos(Y(4)) -sin(Y(4));

0 sin(Y(4))/cos(Y(5)) cos(Y(4))/cos(Y(5))];

Mat = [0 -Y(12) Y(11);

Y(12) 0 -Y(10);

-Y(11) Y(10) 0];

V_xyz = Mat_oRb*[Y(7); Y(8); Y(9)]; %[dxdt; dydt; dzdt]

V_euler = Mat_H*[Y(10); Y(11); Y(12)]; %[dalphadt; dbetadt; dgammadt]

V_uvw = (VF/m)-(Mat*[Y(7); Y(8); Y(9)]); %[dudt; dvdt; dwdt]

V_pqr = inv(IG)*(VM-(Mat*IG*[Y(10);Y(11);Y(12)])); % [dpdt; dqdt; drdt]

dYdt = zeros(12,1);

dYdt(1) = V_xyz(1);

dYdt(2) = V_xyz(2);

dYdt(3) = V_xyz(3);

dYdt(4) = V_euler(1);

dYdt(5) = V_euler(2);

dYdt(6) = V_euler(3);

dYdt(7) = V_uvw(1);

dYdt(8) = V_uvw(2);

dYdt(9) = V_uvw(3);

dYdt(10) = V_pqr(1);

dYdt(11) = V_pqr(2);

dYdt(12) = V_pqr(3);

dYdt = dYdt(:);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ForcesMoments_Aero = F_M (Y,m,TR,rxA,sigma_9g,sigma_10g,sigma_11g,sigma_12g)

rxT = - 0.46;

ryT = 0;

rzT = 0;

rT = [rxT;

ryT;

rzT];

rho = 1.225;

g = -9.81;

D = 0.127;

% rxA = 0.12528;

ryA = 0;

rzA = 0;

CX0 = 0.255;

CX2 = 0.484;

CNA = 3.298;

CLP = -0.042;

CMQ = -1.8;

CLa = 1.683;

CD0a = 0.004;

CD2a = 0.268;

CLg = 0.905;

CD0g = 0.004;

CD2g = 0.111;

rgm = 0.0985;

lgm = -0.154;

lgc = -0.465;

rgc = 0.065;

%% Calcul des forces %%

%% Force de gravite %%

FG = g*m*[-sin(Y(5));

cos(Y(5))*sin(Y(4));

cos(Y(4))*cos(Y(5))];

%% Force aerodynamique du corps missile %%

V = sqrt(Y(7)^2+Y(8)^2+Y(9)^2);

FA =-(pi*rho*V^2*D^2/8)*[CX0+CX2*(Y(8)^2+Y(9)^2)/V^2;

CNA*(Y(8)/V);

CNA*(Y(9)/V)];

%% Calcul des moments %%

%% Moments des forces aerodynamique du corps missile %%

rA = [rxA;

ryA;

rzA];

Moment_A = cross(rA,FA);

%% Moments aerodynamiques non reguliers %%

Moment_UA = (pi*rho*V*D^4/16)*[CLP*Y(10);

CMQ*Y(11);

CMQ*Y(12)];

%% Forces est moments des surfaces additionnelles %%

%% Les ailes %%

FW_a = zeros(3,1);

Moment_a = zeros(3,1);

for j = 1 : 8

phi_a = (j - 1) * (pi/4);

rx_a = lgm;

ry_a = rgm * cos(phi_a);

rz_a = rgm * sin(phi_a);

SLi_a = 0.00175;

gamma_a = 0;

sigma_a = 0;

r_a=[rx_a;

ry_a;

rz_a];

phi_i = phi_a;

gamma_i = gamma_a;

K = Ti (phi_i,gamma_i)' * ([Y(7);Y(8);Y(9)]+([0 -Y(12) Y(11);Y(12) 0 -Y(10);-Y(11) Y(10) 0]*r_a));

V_i = sqrt(K(1)^2+K(2)^2+K(3)^2);

alpha_a = sigma_a + atan(K(3)/K(1));

CL_a = alpha_a * CLa;

CD_a = CD0a + CD2a * alpha_a^2;

Ma = (0.5*rho*SLi_a*V_i^2)*Ti (phi_i,gamma_i)*[CL_a*sin(alpha_a-sigma_a)-CD_a*cos(alpha_a-sigma_a);

0;

-CL_a*sin(alpha_a-sigma_a)-CD_a*cos(alpha_a-sigma_a)];

FW_a = FW_a + Ma;

Moment_a = Moment_a + cross(r_a,Ma);

end

%% Les ailerons %%

V_orientation = [sigma_9g;

sigma_10g;

sigma_11g;

sigma_12g];

FW_g = zeros(3,1);

Moment_g = zeros(3,1);

F_TVC = zeros(3,1);

Moment_TVC = zeros(3,1);

for j = 9 : 12

phi_g = (pi/4) + (j-8) * (pi/2);

rx_g = lgc;

ry_g = rgc * cos(phi_g);

rz_g = rgc * sin(phi_g);

SLi_g = 0.0035;

gamma_g = pi/6;

sigma_g = V_orientation(j-8);

r_g=[rx_g;

ry_g;

rz_g];

phi_i= phi_g;

gamma_i = gamma_g;

K = Ti (phi_i,gamma_i)'*([Y(7);Y(8);Y(9)]+[0 -Y(12) Y(11);Y(12) 0 -Y(10);-Y(11) Y(10) 0]*r_g);

V_i = sqrt(K(1)^2+K(2)^2+K(3)^2);

alpha_g = sigma_g + atan(K(3)/K(1));

CL_g = alpha_g * CLg;

CD_g = CD0g + CD2g * alpha_g^2;

Mg = (0.5*rho* SLi_g *V_i^2)*Ti (phi_i,gamma_i)*[CL_g*sin(alpha_g-sigma_g)-CD_g*cos(alpha_g-sigma_g);

0;

-CL_g*sin(alpha_g-sigma_g)-CD_g*cos(alpha_g-sigma_g);];

FW_g = FW_g + Mg;

Moment_g = Moment_g + cross(r_g,Mg);

F_TVC = F_TVC + (TR/4)*[cos(sigma_g);

sin(phi_g)*sin(sigma_g);

-cos(phi_g)*sin(sigma_g)];

Moment_TVC = Moment_TVC + cross(rT,F_TVC);

end

FW = FW_a + FW_g;

Moment_W = Moment_a + Moment_g;

%% Forces et moments totales %%

F_aero = FG + FA + FW+ F_TVC;

F_aero = F_aero(:);

Fx = F_aero(1);

Fy = F_aero(2);

Fz = F_aero(3);

M_aero = Moment_A + Moment_TVC + Moment_W + Moment_UA;

M_aero = M_aero(:);

L = M_aero(1);

M = M_aero(2);

N = M_aero(3);

ForcesMoments_Aero = [Fx; Fy; Fz; L; M; N];

end

[Alan Stevens]

What is function Ti?

K = Ti(phi_i,gamma_i)' * ([Y(7);Y(8);Y(9)]+....etc.

It doesn't seem to be defined in your listing above.

[MOHANDAREZKI AKTOUF]

Hello Alan and thanks for your answer, Ti is a function that i defined so I can call it everytime I need it.

function Matrice_passage = Ti (phi_i,gamma_i)

Sgi = sin (gamma_i);

Cgi = cos (gamma_i);

Spi = sin (phi_i);

Cpi = cos (phi_i);

Matrice_passage = [Cgi -Sgi 0;

Cpi*Sgi Cpi*Cgi -Spi;

Spi*Sgi Spi*Cgi Cpi];

end

[MOHANDAREZKI AKTOUF]

And it works very well when I test it

Ti(pi/4,0)

ans =

1.0000 0 0

0 0.7071 -0.7071

0 0.7071 0.7071

Answer 1

[Alan Stevens]

I don't know the actual problem, but if you look at your values of Y against time you will see that some of them are clearly diverging towards +/- infinity. That is what is causing Matab to complain. You need to investigate that in some detail.

%%%% Dynamic model integration using ode45 %%%%
t0 = 0; 
tf = 5.09;
h = 0.1;
timerange = [t0 tf]; %t0:h:tf;
IC = [0;0;1;0;deg2rad(-18);0;13;0;0;0;0;0];
                            %% Integration par ode45 %%%
[t,Y] = ode15s(@(t,Y) MDD_Missile(t,Y),timerange,IC);

plot(t,Y)
xlabel('time')
ylabel('Y')
% xlabel('range');
% ylabel('Altitude(m)');
% legend('altitude en fonction de la porteé');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dYdt = MDD_Missile(t,Y)
%%% Missile discrete data %%%
time = [0; 0.3; 0.6; 1.2; 1.8; 2.4; 4.2; 5.2; 15];
mm = [11.25; 11.16; 11.06; 10.82; 10.58; 10.38; 10.16; 10.15; 10.15];
TRR = [0; 570; 650; 750; 770; 650; 50; 0; 0];
rxAA = [0.565; 0.555; 0.544; 0.519; 0.493; 0.471; 0.447; 0.446; 0.446];
Ixxx = [0.0252; 0.025; 0.0248; 0.0244; 0.0239; 0.0235; 0.0231; 0.0231; 0.231];
Iyyy = [0.985; 0.979; 0.973; 0.958; 0.942; 0.929; 0.915; 0.914; 0.914];
Izzz = [0.985; 0.979; 0.973; 0.958; 0.942; 0.929; 0.915; 0.914; 0.914];
%%% Interpolation data-set with t %%%
m = interp1(time,mm,t);
Ixx = interp1(time,Ixxx,t);
Iyy = interp1(time,Iyyy,t);
Izz = interp1(time,Izzz,t);
TR = interp1(time,TRR,t);
rxA = interp1(time,rxAA,t);
IG = [Ixx     0      0;
       0     Iyy     0;
       0      0     Izz];

%% Command parameters (Fins orientation angles) %%
sigma_9g = deg2rad (0);      
sigma_10g = deg2rad (0);
sigma_11g = deg2rad (0);
sigma_12g = deg2rad (0);
%% Call of fonction calculting external forces and moments acting on the missile %%
ForcesMoments_Aero = F_M (Y,m,TR,rxA,sigma_9g,sigma_10g,sigma_11g,sigma_12g);
Fx = ForcesMoments_Aero(1);
Fy = ForcesMoments_Aero(2);
Fz = ForcesMoments_Aero(3);
L = ForcesMoments_Aero(4);
M = ForcesMoments_Aero(5);
N = ForcesMoments_Aero(6);
%%% 12 différential equation of motion of the missile %%%
VF = [Fx; Fy; Fz];
VM = [L; M; N];
Mat_oRb = [cos(Y(5))*cos(Y(6))   cos(Y(6))*sin(Y(5))*sin(Y(4))-sin(Y(6))*cos(Y(4))   cos(Y(6))*sin(Y(5))*cos(Y(4))+sin(Y(6))*sin(Y(4));
           sin(Y(6))*cos(Y(5))   sin(Y(6))*sin(Y(5))*sin(Y(4))+cos(Y(6))*cos(Y(5))   sin(Y(6))*sin(Y(5))*cos(Y(4))-cos(Y(6))*sin(Y(4));
          -sin(Y(5))             cos(Y(5))*sin(Y(4))                                 cos(Y(5))*cos(Y(4))];
Mat_H = [1     sin(Y(4))*tan(Y(5))      cos(Y(4))*tan(Y(5));
         0     cos(Y(4))               -sin(Y(4));
         0     sin(Y(4))/cos(Y(5))      cos(Y(4))/cos(Y(5))];
Mat = [0       -Y(12)     Y(11); 
       Y(12)     0       -Y(10);  
      -Y(11)    Y(10)       0];
V_xyz = Mat_oRb*[Y(7); Y(8); Y(9)];                    %[dxdt; dydt; dzdt]
V_euler = Mat_H*[Y(10); Y(11); Y(12)];                 %[dalphadt; dbetadt; dgammadt]
V_uvw = (VF/m)-(Mat*[Y(7); Y(8); Y(9)]);               %[dudt; dvdt; dwdt]
V_pqr = IG\(VM-(Mat*IG*[Y(10);Y(11);Y(12)]));     % [dpdt; dqdt; drdt]
dYdt = zeros(12,1);                    
dYdt(1) = V_xyz(1);
dYdt(2) = V_xyz(2);
dYdt(3) = V_xyz(3);
dYdt(4) = V_euler(1);
dYdt(5) = V_euler(2);
dYdt(6) = V_euler(3);
dYdt(7) = V_uvw(1);
dYdt(8) = V_uvw(2);
dYdt(9) = V_uvw(3);
dYdt(10) = V_pqr(1);
dYdt(11) = V_pqr(2);
dYdt(12) = V_pqr(3);
dYdt = dYdt(:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ForcesMoments_Aero = F_M (Y,m,TR,rxA,sigma_9g,sigma_10g,sigma_11g,sigma_12g)  
rxT = - 0.46;
ryT = 0;
rzT = 0;
rT = [rxT;
      ryT;
      rzT];
rho = 1.225;
g = -9.81;
D = 0.127;
% rxA = 0.12528;
ryA = 0;
rzA = 0;
CX0 = 0.255;
CX2 = 0.484;
CNA = 3.298;
CLP = -0.042;
CMQ = -1.8;
CLa = 1.683;
CD0a = 0.004;
CD2a = 0.268;
CLg = 0.905;
CD0g = 0.004;
CD2g = 0.111;
rgm = 0.0985;
lgm = -0.154;
lgc = -0.465;
rgc = 0.065;
                        %% Calcul des forces %%
%% Force de gravite %% 
FG = g*m*[-sin(Y(5)); 
           cos(Y(5))*sin(Y(4)); 
           cos(Y(4))*cos(Y(5))];
%% Force aerodynamique du corps missile %%
V = sqrt(Y(7)^2+Y(8)^2+Y(9)^2);
FA =-(pi*rho*V^2*D^2/8)*[CX0+CX2*(Y(8)^2+Y(9)^2)/V^2;
                         CNA*(Y(8)/V);
                         CNA*(Y(9)/V)];
                         %% Calcul des moments %%
%% Moments des forces aerodynamique du corps missile %%
rA = [rxA;
      ryA;
      rzA];
Moment_A = cross(rA,FA);
%% Moments aerodynamiques non reguliers %%
Moment_UA = (pi*rho*V*D^4/16)*[CLP*Y(10);
                               CMQ*Y(11);
                               CMQ*Y(12)];
%% Forces est moments des surfaces additionnelles %%
%% Les ailes %%
FW_a = zeros(3,1);
Moment_a = zeros(3,1);
for j = 1 : 8
    phi_a = (j - 1) * (pi/4);
    rx_a = lgm;
    ry_a = rgm * cos(phi_a);
    rz_a = rgm * sin(phi_a);
    SLi_a = 0.00175;
    gamma_a = 0;
    sigma_a = 0;
    r_a=[rx_a;
          ry_a;
          rz_a];
    phi_i = phi_a;
    gamma_i = gamma_a; 
    K = Ti(phi_i,gamma_i)' * ([Y(7);Y(8);Y(9)]+([0 -Y(12) Y(11);Y(12) 0 -Y(10);-Y(11) Y(10) 0]*r_a));
    V_i = sqrt(K(1)^2+K(2)^2+K(3)^2);
    alpha_a = sigma_a + atan(K(3)/K(1));
    CL_a = alpha_a * CLa;
    CD_a = CD0a + CD2a * alpha_a^2;
    Ma = (0.5*rho*SLi_a*V_i^2)*Ti (phi_i,gamma_i)*[CL_a*sin(alpha_a-sigma_a)-CD_a*cos(alpha_a-sigma_a);
                                                   0;    
                                                  -CL_a*sin(alpha_a-sigma_a)-CD_a*cos(alpha_a-sigma_a)];
    FW_a = FW_a + Ma;
    Moment_a = Moment_a + cross(r_a,Ma);                                                               
end
                              %% Les ailerons %%
V_orientation = [sigma_9g;
                 sigma_10g; 
                 sigma_11g;
                 sigma_12g];
FW_g = zeros(3,1);
Moment_g = zeros(3,1);
F_TVC = zeros(3,1);
Moment_TVC = zeros(3,1);
for j = 9 : 12
    phi_g = (pi/4) + (j-8) * (pi/2);
    rx_g = lgc;
    ry_g = rgc * cos(phi_g);
    rz_g = rgc * sin(phi_g);
    SLi_g = 0.0035;
    gamma_g = pi/6;
    sigma_g = V_orientation(j-8);
    r_g=[rx_g;
         ry_g;
         rz_g];
    phi_i= phi_g;
    gamma_i = gamma_g;
    K = Ti(phi_i,gamma_i)'*([Y(7);Y(8);Y(9)]+[0 -Y(12) Y(11);Y(12) 0 -Y(10);-Y(11) Y(10) 0]*r_g);
    V_i = sqrt(K(1)^2+K(2)^2+K(3)^2);
    alpha_g = sigma_g + atan(K(3)/K(1));
    CL_g = alpha_g * CLg;
    CD_g = CD0g + CD2g * alpha_g^2;
    Mg = (0.5*rho* SLi_g *V_i^2)*Ti (phi_i,gamma_i)*[CL_g*sin(alpha_g-sigma_g)-CD_g*cos(alpha_g-sigma_g);
                                                      0;    
                                                     -CL_g*sin(alpha_g-sigma_g)-CD_g*cos(alpha_g-sigma_g);];                                                                 
   
    FW_g = FW_g + Mg;
    Moment_g = Moment_g + cross(r_g,Mg);
    F_TVC = F_TVC + (TR/4)*[cos(sigma_g); 
                            sin(phi_g)*sin(sigma_g); 
                           -cos(phi_g)*sin(sigma_g)];
    Moment_TVC = Moment_TVC + cross(rT,F_TVC);
end
FW = FW_a + FW_g;
Moment_W = Moment_a + Moment_g;
%% Forces et moments totales %%
F_aero = FG + FA + FW+ F_TVC;
F_aero = F_aero(:);
Fx = F_aero(1);
Fy = F_aero(2);
Fz = F_aero(3);
M_aero = Moment_A + Moment_TVC + Moment_W + Moment_UA;
M_aero = M_aero(:);
L = M_aero(1);
M = M_aero(2);
N = M_aero(3);
ForcesMoments_Aero = [Fx; Fy; Fz; L; M; N];
end
function Matrice_passage = Ti (phi_i,gamma_i)
Sgi = sin (gamma_i);
Cgi = cos (gamma_i);
Spi = sin (phi_i);
Cpi = cos (phi_i);
Matrice_passage = [Cgi       -Sgi      0;
                   Cpi*Sgi   Cpi*Cgi  -Spi;
                   Spi*Sgi   Spi*Cgi   Cpi];
end

[MOHANDAREZKI AKTOUF]

Thanks again Alan, I am checking everyday but I can't really find the source of the problem, I will re-check again maybe there's a problem in one of my functions, Can I ask why did you add a value to each data vector?

[MOHANDAREZKI AKTOUF]

Hello Alan, after rechecking, I got this graph but always with the same error :

Warning: Failure at t=4.802987e+00. Unable to meet integration tolerances without reducing the step size below the smallest value

allowed (1.421085e-14) at time t.

[Alan Stevens]

Your original end time went beyond the limits of your interpolation vectors. I thought this might be a cause of the problem - it wasn't!

Your program is too complicated or me to work out what is going on from the listing. If you were able to upload a mathematical model of the system I might take a look (though I promise nothing!).