Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
simulation of a variable mass with ode45

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 7 Feb, 2008 11:37:01

Message: 1 of 9

Hello,
I am trying to simulate a simple model of a rocket using ode45.

The problem is that the mass "M" of the rocket is a function
of time "t" and the thrust "Th", which is a function of the
velocity "dx" and position "x" of the rocket --> M(t, Th(dx,x)).

In order to express (in the file containing the equation of
the rocket) the variation of "M" I need:

1.) to know what is the size of the current integration step
taken by the solver;
2.) a way to preserve the current value of "M" in order to
use it during the next integration step;
3.) A way to output the profiles of "M" and "Th" after the
simulation is over.

I have implemented this in Simulink and there I could use
the "Memory block".

Any help is appreciated!
Dimitar

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 8 Feb, 2008 08:24:03

Message: 2 of 9

I wonder weather my question is too trivial!

Just as a follow up, I would like to mention, that with a
simple implementation of any order Runge-Kutta, it is
straightforward to solve this problem, by simply adding the
mass as an output (and input) to the (matlab) function
containing the equations for the rocket. But with ode45 this
method doesn't work or ...?

Cheers,
Dimitar

Subject: simulation of a variable mass with ode45

From: Jacek Kierzenka

Date: 8 Feb, 2008 22:18:02

Message: 3 of 9

Dimitar,

One problem is that in MATLAB ODE solvers, the current step
size is never passed to the derivative function. Can you
somehow work around that?

-Jacek

-------------------------

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in
message <foeqgt$s42$1@fred.mathworks.com>...
> Hello,
> I am trying to simulate a simple model of a rocket using
ode45.
>
> The problem is that the mass "M" of the rocket is a
function
> of time "t" and the thrust "Th", which is a function of
the
> velocity "dx" and position "x" of the rocket --> M(t, Th
(dx,x)).
>
> In order to express (in the file containing the equation
of
> the rocket) the variation of "M" I need:
>
> 1.) to know what is the size of the current integration
step
> taken by the solver;
> 2.) a way to preserve the current value of "M" in order to
> use it during the next integration step;
> 3.) A way to output the profiles of "M" and "Th" after the
> simulation is over.
>
> I have implemented this in Simulink and there I could use
> the "Memory block".
>
> Any help is appreciated!
> Dimitar
>
>

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 9 Feb, 2008 07:43:02

Message: 4 of 9

> One problem is that in MATLAB ODE solvers, the current step
> size is never passed to the derivative function. Can you
> somehow work around that?


Hello Jacek,
Thank you for the post!

In my notation, the variation of the mass depends on the
Thrust, and the Thrust depends on the position and velocity
of the rocket, so since I cannot precompute the Thrust, the
only option I have is to update the mass on-line.

The relation mass --> Thrust is:

M = Fuel_cons*Th*d_time

with "d_time" being the time step, and "Fuel_cons" the fuel
needed to produce 1 N of force for a duration of 1 s.

So I definitely need to have d_time.

If would have been perfect if we could use a "derivative
function" like:

--------------------------------------------------------
function [dy, old_time, M] = EM_rocket(t,y, old_time, M)

delta_t = t - old_time;
old_time = t;

if delta_t > 0
  M = M - Fuel_cons*Th*d_time;
end
--------------------------------------------------------

The above thing can be implemented in Simulink by simply
using the "Memory block". I presume that there is a similar
technique in Matlab?


Regards,
Dimitar

Subject: simulation of a variable mass with ode45

From: Jacek Kierzenka

Date: 11 Feb, 2008 14:53:12

Message: 5 of 9

You could accomplish something similar using and OutputFcn (see solver
options available via ODESET).
In that function, you would store the last accepted time step (old_time). It
would need to be stored in some
sort of a global variable, visible to both the OutputFcn and your derivative
function.
That should do the trick.

Hope that helps,
Jacek

-------------------------------

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message
news:fojli6$57b$1@fred.mathworks.com...
>> One problem is that in MATLAB ODE solvers, the current step
>> size is never passed to the derivative function. Can you
>> somehow work around that?
>
>
> Hello Jacek,
> Thank you for the post!
>
> In my notation, the variation of the mass depends on the
> Thrust, and the Thrust depends on the position and velocity
> of the rocket, so since I cannot precompute the Thrust, the
> only option I have is to update the mass on-line.
>
> The relation mass --> Thrust is:
>
> M = Fuel_cons*Th*d_time
>
> with "d_time" being the time step, and "Fuel_cons" the fuel
> needed to produce 1 N of force for a duration of 1 s.
>
> So I definitely need to have d_time.
>
> If would have been perfect if we could use a "derivative
> function" like:
>
> --------------------------------------------------------
> function [dy, old_time, M] = EM_rocket(t,y, old_time, M)
>
> delta_t = t - old_time;
> old_time = t;
>
> if delta_t > 0
> M = M - Fuel_cons*Th*d_time;
> end
> --------------------------------------------------------
>
> The above thing can be implemented in Simulink by simply
> using the "Memory block". I presume that there is a similar
> technique in Matlab?
>
>
> Regards,
> Dimitar
>
>
>
>

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 12 Feb, 2008 16:29:03

Message: 6 of 9

Hello Jacek,

After thinking about your suggestion,
I came to the conclusion that since the OutputFcn is called
after each successful iteration I could directly compute the
variable mass in the OutputFcn. As for the time step I can
compute it like d_time = t_o(4) - t_o(1) (assuming I use
ode45). Of course, in order to have a reliable simulation I
would have to limit the MAX step size.

Using the above method I will be able to change the mass
every four steps (using the "accumulated" Thrust during
d_time), which can be an "OK" approximation if d_time is
"sufficiently" small.

I have two question on how to implement the above thing:
1. How can I use the history of the Thrust during d_time
(computed in the derivative function) in my Output function.
Do I have to define it as a global variable, or there is a
way to directly output it?

2. The mass should definitely be a global variable, however,
I am a bit stuck with how to define this in a simple and
intuitive way.

Of course if you have an idea about an alternative way of
implementing the problem I would appreciate if you share it
with me!

Best regards,
Dimitar

Subject: simulation of a variable mass with ode45

From: Jacek Kierzenka

Date: 12 Feb, 2008 17:16:10

Message: 7 of 9

The OutputFcn only returns the 'status' flag, and the ODEFcn only returns
the derivative vector.
I guess, you are stuck with globals (of some sort) if you need to pass info
between them...

To accumulate info within one function, you could use a persistent variable
(in case you need to
store the time of two consecutive steps.)

If your OutputFcn and ODEFcn happen to be nested functions, they both have
access to local
variables defined in the outer function -- you could use those to exchange
information, shielding
them from being visible to the whole MATLAB session.

Hope that helps,
Jacek
------------------

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message
news:foshgf$5sg$1@fred.mathworks.com...
> Hello Jacek,
>
> After thinking about your suggestion,
> I came to the conclusion that since the OutputFcn is called
> after each successful iteration I could directly compute the
> variable mass in the OutputFcn. As for the time step I can
> compute it like d_time = t_o(4) - t_o(1) (assuming I use
> ode45). Of course, in order to have a reliable simulation I
> would have to limit the MAX step size.
>
> Using the above method I will be able to change the mass
> every four steps (using the "accumulated" Thrust during
> d_time), which can be an "OK" approximation if d_time is
> "sufficiently" small.
>
> I have two question on how to implement the above thing:
> 1. How can I use the history of the Thrust during d_time
> (computed in the derivative function) in my Output function.
> Do I have to define it as a global variable, or there is a
> way to directly output it?
>
> 2. The mass should definitely be a global variable, however,
> I am a bit stuck with how to define this in a simple and
> intuitive way.
>
> Of course if you have an idea about an alternative way of
> implementing the problem I would appreciate if you share it
> with me!
>
> Best regards,
> Dimitar
>

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 12 Feb, 2008 19:37:01

Message: 8 of 9

Now I understand what you meant by "sort of a global
variable" (thank you). I implemented a simplified version,
and (for the moment) I can see only one problem:

At Iteration 1, I get (in the OutputFunction) for example:
-------------
t(1) = 0.000005
t(2) = 0.000010
t(3) = 0.000015
t(4) = 0.000020

At Iteration 2 I get:
-------------
t(1) = 0.000029
t(2) = 0.000037
t(3) = 0.000046
t(4) = 0.000054

so if I form d_time = t(4) - t(1)

at each Iteration, I will lose the time t(4) (iteration 1)
to t(1) (iteration 2), and so on.

In order to illustrate the problem I attached an example
(if you copy and past it should run). I solve it by saving
t(4) in a temporary file (see variable "tmp") and then
reading it (from this file) every time in order to form the
step as:

d_time = t(4) - tmp,


The example:
-----------

function main1
%
% Simplified rocket example using ode45 with nested functions

% Parameters (global)
%----------------
g = 9.81; % acceleration due to gravity
k = 0.00021525; %(Air drag)
Fuel_use = 0.0005; % in order to produce 1N force for 1 s we
need 0.0005 kg fuel

mr = 0.25; % mass of the rocket
mf = 0.25; % mass of the fues
m = mr+mf; % initial mass of the rocket with fuel

Th = 0.0; % initial Thrust

GAIN = 2.1525; % control gain

MAX_ST = 0.05; % maximum step size (for ode45)
%----------------

% IMPORTANT ???:
tmp = 0; % t = 0

fid=fopen('tmp.dat','w');
fprintf(fid, '%f ', tmp);
fclose(fid);

%%%% CALL TO ODE45
options = odeset('OutputFcn',@OutFun, 'MaxStep', MAX_ST);
[T, Y] = ode45(@EM_rocket, [0 250], [0 ; 0], options);

color = 'r';

figure(1)
subplot(2,1,1);hold on;grid on;title('Position')
plot(T,Y(:,1), color);
subplot(2,1,2);hold on;grid on;title('Velocity')
plot(T,Y(:,2), color);

%---------------------------------------------------
    function [dy]= EM_rocket(t,y)

        % Rocket EM
        %----------------
        dy = zeros(2,1);
        dy(1) = y(2);
        dy(2) = Th/m - k*sign(y(2))*y(2)^2/m - g;

    end
%---------------------------------------------------
    function status = OutFun(t,y,flag)

        n = length(t);
        if n ~= 2 && n ~= 0

            tmp = load('tmp.dat'); % t(4) during the
previous iteration

            d_time = t(n) - tmp;

            % d_time = t(n) - t(1); % this is wrong because
the time between
            % t(4) (iteration 1) and t(1) (iteration 2) is
"lost" and so on
            %
            % t(1) t(2) t(3) t(4)
            % 0.000005 0.000010 0.000015 0.000020
iteration 1
            % 0.000029 0.000037 0.000046 0.000054
iteration 2
            % 0.000063 0.000071 0.000079 0.000088
iteration 3
            % ...

            % simple control
            if m > mr
                Th = m*g*GAIN;
            else
                Th = 0;
            end

            m = m - Th*d_time*Fuel_use;

            % IMPORTANT ???:
            fid=fopen('To.dat','a');
            fprintf(fid, '%f \t %f \t %f \t %f \t %f \n', t,
d_time);
            fclose(fid);

            fid=fopen('tmp.dat','w');
            fprintf(fid, '%f ', t(end));
            fclose(fid);


        end
        status = 0;

    end
%---------------------------------------------------
end


Regards,
Dimitar

Subject: simulation of a variable mass with ode45

From: Dimitar Dimitrov

Date: 12 Feb, 2008 20:26:01

Message: 9 of 9

Thats why people have to take a rest!!!

Please disregard the previous post. The "tmp" variable could
be defined as a "global" variable as well :). And everything
works.

So something as a conclusion:
-------------------
For the example in my previous post the system is NOT stiff,
and an approximation to update the mass and the Thrust every
four "steps" is OK.
However, for a stiff closed loop system, such method will
not lead to good results, so my guess is that for a stiff
system with a variable mass (as in my case) the ODE
integrators can not be used, or ? (quite a bold statement)

Probably there is a different approach to the problem?
Any ideas?

Best regards,
Dimitar

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us