Why does ode45 output "NaN" after some time, with this time depending on the input

Question

Merritt 2023-11-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2045755-why-does-ode45-output-nan-after-some-time-with-this-time-depending-on-the-input

评论： William Rose 2023-11-10

Hi, I am modeling a rocket launching by using ode45 to solve ODEs. I am having a strange problem where for the base problem where the rocket has 1 stage and therefore 1 input for the propellant mass, it is solved with no problem and plotted (rocket height vs. time) as shown below, but when I add another stage to the rocket and model it, ode45 outputs "NaN" starting at seemingly random timesteps.

When I change the script to account for a 2nd stage for the rocket, ode45 keeps running into instances where it can only calculate "NaN" for the solution to the odes, and this makes little sense to me because it also doesn't make sense for the derivative states to suddenly be calculated as Inf. I am able to plot the numbers it does get, and I do not understand why ode45 starts to calculate "NaN" at certain values; it almost seems random where it encounters this problem (depending on the 2 input propellant masses, of course). Shown below are the results from ode45 for the imputs of [5000;5000], [7000;7000], and [10000;0], respectively; strangely the [10000;0] case is calculated as it should be and gives a solution that I would expect. The dotted lines are the time of burnout for each stage of the rocket (so ode45 seems to have a very hard time with the second stage, but not the first). The code including the ODEs to be solved by ode45 is attached to this post as "heightcont_2stage.m". The code I was using to call the ODE function is attached to this post as "testindstudy.m".

I guess the main thing I wanted to ask was if the addition of "if" statements to govern when the equations of motion switch from the 1st stage equations to second stage equations would be what is causing this strange output, because I suspect that is what might be the problem (althought that is a guess at best). The specific code governing which stage the rocket is currently in (and what mass/acceleration to use) is shown below. There are no instances of singularities as far as I can tell either, so the reason behind the "NaN" is beyond me.

if t <= tbo1     % Burning DURING 1st stage
    mt = m0 - mdot1*t;
    vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo1 && t <= tbo2     % Burning DURING 2nd stage
    mt = m0 - mP1 - mS1 - mdot2*t;
    vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo2     % Coasting AFTER 2nd stage
    mt = m0 - mP1 - mS1 - mP2;
    vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
end

I have tried solving these ODEs with several different solvers (ode45, ode15s, ode23, ode89) and they produce the same results, so it has to be something not meshing well with my code. Just wanted some guesses into what might be causing this.

Any guesses as to why this is happening are welcome, thank you.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2023-11-10

编辑：Torsten 2023-11-10

在 MATLAB Online 中打开

clear all

close all

u_2stage = [7000;7000];

u = [10000];

tmax = 3000;

hmax = 2000;

%[heightconstraint] = heightconst(tmax,u,hmax);

[heightconstraint_2stage,xout] = heightconst_2stage(tmax,u_2stage,hmax);

Warning: Failure at t=1.383550e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t.

function [constraint,xout] = heightconst_2stage(tmax,u0,hmax)

% Calculates the height of a rocket through the solution of ODEs

%

% Inputs:

% tmax = Maximum time to run the ODE solver for (s)

% x0 = Initial conditions for height (km), velocity (km/s), theta (rad),

% phi (rad), and psi (rad)

% u0 = Mass of propellant for each stage

%

% Outputs:

% height = Height of rocket (km)

% Solve ODEs for equations of motion

tspan = [0,tmax];

odeoptions = [];

x0 = [0;0;0;0;0];

[tout,xout] = ode15s(@states,tspan,x0,odeoptions,u0);

height = xout(:,1);

maxhreached = max(height);

constraint = hmax - maxhreached;

figure(2)

plot(tout,xout(:,1),'b-')

hold on

grid on

xline(120,'k--')

xline(240,'k--')

ylim([0,5000])

for i = 1:numel(tout)

[~,mt(i)] = states(tout(i),xout(i,:),u0);

end

figure(3)

plot(tout,mt)

% Define state space representation for use with ode45

function [output,mt] = states(t,x,u)

% Calculates the state space representation for the objective

% function for use in the ode45 MATLAB function

%

% Inputs:

% t = Current time step (s)

% x = Height (km), velocity (km/s), theta (rad), phi (rad), and psi (rad)

% for the previous time step

% u = Stage propellant masses mP1 (kg) and mP2 (kg)

%

% Outputs:

% hdot = Time derivative of height (km)

% vdot = Time derivative of velocity (km/s)

% thetadot = Time derivative of theta (rad/s)

% phidot = Time derivative of phi (rad/s)

% psidot = Time derivative of psi (rad/s)

Isp = 353; % Isp of propellant (s)

epsilon = .01; % Rocket structural ratio

mPL = 1000; % Payload mass (kg)

d = 1.65e-3; % Radius of rocket (km)

tbo1 = 120; % Time to burnout for 1st stage (s)

tbo2 = 240; % Time to burnout for 2nd stage (s)

cd = .3; % Drag coefficient

RE = 6378; % Radius of Earth (km)

rho0 = 1.225e9; % Density of air at sea level (kg/km^3)

h0 = 7.5; % Reference height for calculation of rho (km)

g0 = 9.81e-3; % Gravitational acceleration at sea level (km/s^2)

mu = 398600; % Gravitational constant for Earth (km^3/s^2)

mP1 = u(1,1);

mP2 = u(2,1);

rho = rho0*exp(-x(1)/h0); % Air density for current height (kg/km^3)

mS1 = mP1*(epsilon/(1 - epsilon)); % Structural mass of 1st stage (kg)

mS2 = mP2*(epsilon/(1 - epsilon)); % Structural mass of 2nd stage (kg)

m0 = mPL + mP1 + mS1 + mP2 + mS2; % Initial mass of rocket (kg)

A = (pi/4)*d^2; % Cross-sectional area of rocket (km^2)

mdot1 = mP1/tbo1; % Rate of propellant use for 1st stage (kg/s)

mdot2 = mP2/(tbo2 - tbo1); % Rate of propellant use for 2nd stage (kg/s)

if x(1) < 0

x(1) = 0;

end

if t <= tbo1

mt = m0 - mdot1*t;

vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));

elseif t > tbo1 && t <= tbo2

mt = m0 - mP1 - mS1 - mdot2*t;

vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));

elseif t > tbo2

mt = m0 - mP1 - mS1 - mP2;

vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));

end

if x(2) == 0

phidot = 0;

else

phidot = (mu*sin(x(5)))/(((RE + x(1))^2)*x(2));

end

hdot = x(2)*cos(x(5));

thetadot = (x(2)*sin(x(5)))/(RE + x(1));

if x(1) > .5 && x(5) < .1*(pi/180)

psidot = .1*(pi/180);

else

psidot = phidot - thetadot;

end

output = [hdot;vdot;thetadot;phidot;psidot];

end

Torsten 2023-11-10

编辑：Torsten 2023-11-10

It seems I found the reason for failure - you divide by mt which becomes 0 in the course of integration (see above).

Merritt 2023-11-10

This seems to have been the issue! Thank you, you are a lifesaver.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2023-11-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2045755-why-does-ode45-output-nan-after-some-time-with-this-time-depending-on-the-input#answer_1350835

在 MATLAB Online 中打开

@Merritt,

There was an error when I ran testinstudy, because it calls heightconst(), which you did not provide. SO I commented out that line. It runs fine with u_2stage=[10000;0]. It runs fine with u_2stage=[5000;5000]: it doesn;t get very high, but it runs without NaNs. With u_2stage=[7000;7000], I get the result

>> testindstudy

Warning: Failure at t=2.406072e+02. Unable to meet integration tolerances without reducing

the step size below the smallest value allowed (4.547474e-13) at time t.

> In ode45 (line 353)

In heightconst_2stage (line 17)

In testindstudy (line 12)

ode45() is giving an error because it keeps trying to reduce the stepsize to get a smooth solution, and it gets to its minimum step size without finding a smooth solution. So it bails out. I have had similar problems when modeling the cardiovascular system, with one-way valves. My solution was to replace a statement like

if pressureDiff<0, 
    flow=0; 
else 
    flow=pressureDiff/resistance; 
end

with a smooth function such as

flow=pressureDiff/(resistance*(1+exp(-k*pressureDiff));

where k is large. This avoided a discontinuity in the equations, and stopped the crashing of ode45() in my simulation. I recommend you do something equivalent at the times of stage separation.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

William Rose 2023-11-10

在 MATLAB Online 中打开

@Merritt,

Elaborating on what I said above:

Replace

if t <= tbo1
    mt = m0 - mdot1*t;
    vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo1 && t <= tbo2
    mt = m0 - mP1 - mS1 - mdot2*t;
    vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo2
    mt = m0 - mP1 - mS1 - mP2;
    vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
end

with

tau=0.5;  % transition duration (s)
mt = m0...
    - mdot1*t/(1+exp((t-tbo1)/tau))...      % turn off -mdot1*t when t>tbo1
    - (mP1+mS1)/(1+exp(-(t-tbo1)/tau))...   % turn on -mP1-mS1 when t>tbo1
    - mdot2*t/((1+exp(-(t-tbo1)/tau))*(1+exp((t-tbo2)/tau)))...  % turn on mdot2*t when tbo1<t<tbo2
    - mP2/(1+exp(-(t-tbo2)/tau));           % turn on -mP2 when t>tbo2
vdot = - 0.5*rho*cd*A*x(2)^2/mt - mu*cos(x(5))/(RE + x(1))^2...
       + (1/mt)*Isp*g0*mdot1/(1+exp((t-tbo1)/tau))...   % turn off +mdot1 when t>tbo1
       + (1/mt)*Isp*g0*mdot2/((1+exp(-(t-tbo1)/tau))*(1+exp((t-tbo2)/tau))); % turn on +mdot2 when tbo1<t<tbo2

Check my work, because I might have made some typos.

I don;t promise it will work, but this is the kindof think I had to do.

I am using sigmoidal funcitons that go from 0 to 1, or 1 to 0, to turn on or turn off different terms for different values of t. This gives the general idea.

请先登录，再进行评论。