Display output of ODE45 after error

Hello Programmers,

I have a model that predicts the altitude (Z) of a UAV.

The model is a sets of ODE and I solved it using ODE45.

The issue is that the UAV altitude (Z) decreases to a negative value and the program terminate as Z cannot be negative.

If Z = 0, it means that the UAV is on the surface of the earth.

Instead of giving error when Z is negative. I will like to display a message that "The object has successfully landed at location, Longitude = .. and Latitude = " and return the output of ODE45. So, I can plot all the output variables.

Please I will appreciate any help.

Thank you.

Answer 1

[Amit Bhowmick]

I dont know your equation but you can use similar trick

tspan = [0 5];
y0 = 20;
[t,y] = ode45(@fun, tspan, y0);
plot(t,y,'-o')


function dydt=fun(t,y)
if(y<=0)
    dydt=0;
else
    dydt=-2*t;
end
end

[Walter Roberson]

No, this will not work.

If the object is starting from ground to move upwards then it will have an initial velocity, but you would cancel that out.

[Amit Bhowmick]

Is your initial velocity -ve ? It doest not mean one should use the similar logic. Logic to be developed to control the output. I think event function does the same thing. Here in the this example y is the height which should not be negetive and negetive y means the object is travelling beyond the ground where if you inplemnt it using implicit or explicit euler, RK same logic to be adopted where equation will look like,

dy/dt=-2*t (t>0,y=>0)

Kindly follow the trajectory of the particle from the plot

[Walter Roberson]

The mathematics of the ode*() functions requires that the first two derivatives of the function be continuous within the range of values evaluated. Using an if rarely gives the required continuity.

An event function set to fire when the height decreases through 0 would not have problems with continuity of derivatives.

[Amit Bhowmick]

In this case the negetive height implies body at rest, with if condtion that is what we can acheive by making velocity zero after plasting colison during landing. In otherway the differential equation is not valid when y becomes negetive as you can see from equation assumption or constarint. Also from the plot of particle trajectory attached it is clear that after the particle touched the ground particle becomes static.

[Walter Roberson]

One of the most dangerous problems with ode45() and related functions is that they cannot reliably detect when they are returning incorrect solutions.

In a case such as this, instead of solving the original equations, ode45() would be solving something similar to multiplying the original equations by an exponential backoff function. The difference compared to the original equations would be in the details right near the landing, with the landing either ending up being too early, or else with there being an overshoot.

[Walter Roberson]

Velocity must be informed by height. If you do not have current height as a boundary condition to the ODE, then the UAV would have to assume that if velocity is negative that the UAV must go full power to slow down to velocity 0 -- and then potentially get stuck hovering in mid-air until fuel ran out.

    dydt=-2*t;

In an UAV situation, you do not decelerate more strongly according to how long you have been flying. There is a maximum lift that the UAV is capable of.

[Amit Bhowmick]

Please follow the quote "I dont know your equation but you can use similar trick" hence i did not write the code for UAV rather it's a free fall case where a particle can not go below the ground i.e. negetive y.

at the zero point after colison the velocity will be governed by two factor,

1. the value of coefficient of restitution(e) and 2. Velocity before colision

And the equation mentioned correctly represents the case of plastic colision e=0. also there will be discontinuity in y as well as dydt (downward and upward velocities ) which suggest modification of function is required (since at y=0 the velocities exist one before colison and other after colision). Hence ode45 will keep producing results coinsidering velocity -2t such that particle moves below the ground. Therefore one need to modify the function with two values at ground so that ode45 can produce correct results.

For UAV landing there must be a deaceleration by some sort jet to neutralize the affect of gravity and accordingly ode should be governed for perfect landing. However, in an actual case the perfect landing never possible there have to be some sort of colision at the ground. If your equation is not fit enought to consider the colision then it will go below the ground.

[Telema Harry]

Thank you for the contribution. But my model is slightly different

Answer 2

[Walter Roberson]

You should use an event function to terminate integration. See the ballode example.

Note that in order to be able to detect height being a particular value, you would probably need to have height be one of the boundary conditions. You might not need its integral on output, but you need it for the steering calculations. For example, the UAV is going to make different decisions about how much lift it needs to generate if the UAV is 50m off the ground than if the UAV is 50cm off the ground.

[Telema Harry]

Thank you for the feedback.

I have attached a graph of my altitude profile and a sample of my code as further clarification.

If I run the simulation for 3000 seconds for example. The height will be less than 0 before the end time and the program will terminates and with an error.

I have used and if statement to track the moment the altitude turns into a negative value.

sample of my code is shown below.

[t,output] = ode45('BalloonODE',tspan,In0);

function dxdt = BalloonODE(t,y)   
      
    Z         = y(7);
    
      if sign(Z) == -1
           error('The Balloon has successfully landed on the surface')
       %  I will like the ODE solver to terminates at this point and return all the previous values of y.    
      end
      
      dZ = Uz;
      
       dxdt = [dlat, dlon dUz, dZ]'
       
end

Thank you for your help.

[Telema Harry]

I have used the code in ballode function

I also saw a nice code from @Jan at https://www.mathworks.com/matlabcentral/answers/397501-how-to-stop-ode45-when-value-reach-certain-value-other-than-zero .

but I am getting an error and I do not know how to resolve the error. Please see a screenshot of the error.

   options = odeset('Events',@stopevents);
    
     In0 = [loc.lat; loc.lon; In.Uz; In.Z]';
     
[t,output] = ode45('BalloonODE',tspan,In0,options);

function [value,isterminal,direction] = stopevents(t,y)

% Locate the time when height passes through zero in a decreasing direction
% and stop integration.
    value = y(7);     % detect height = 0
    isterminal = 1;   % stop the integration
    direction = 0;   % negative direction


end

[Walter Roberson]

What shows up for

which -all lower

[Telema Harry]

Please I do not understand.

[Walter Roberson]

At the MATLAB command line give the command

which -all lower

and show us the output.

[Telema Harry]

Thank you. Please see a screen shot of the output

[Jan]

Use the debugger for investigations:

dbstop if error

Then run your code again. When Matlab stops at the error, check the contents of "eventFcn".

I guess boldly, that the problem is caused as side-effect by providing the function to be integrated by its name instead of a function handle. This notation is outdated for over 20 years, so change:

[t,output] = ode45('BalloonODE',tspan,In0,options);

to

[t,output] = ode45(@BalloonODE,tspan,In0,options);

[Telema Harry]

woooow. thank you so much @JanJan. The code now runs without the error.

However, I have not been able to get my desired result. The Z is still decreasing below 0 and P_air can only be calculated when Z > 0. So the program terminates.

Thank you for the help.

   options = odeset('Events',@stopevents);
    
     In0 = [lat_0; lon_0; T_gas_0; T_film_0; M_gas_0; Uz_0; Z_0]';
      [t,output] = ode45(@BalloonODE,tspan,In0,options);


function [value,isterminal,direction] = stopevents(t,y)

% Locate the time when height passes through zero in a decreasing direction
% and stop integration.
    value = y(7);     % detect height = 0
    isterminal = 1;   % stop the integration
    direction = -1;   % negative direction


end