2

I have a MATLAB code that solves a large scale ODE system of following type

     function [F] = myfun(t,C,u)

u here is a time dependent vector used as a input in the function myfun. My current code reads as:-

 u    = zeros(4,1);
 u(1) = 0.1; 
 u(2) = 0.1;
 u(3) = 5.01/36;
 u(4) = 0.1;

C0            = zeros(15*4*20+12,1);
options = odeset('Mass',Mf,'RelTol',1e-5,'AbsTol',1e-5);
[T,C]   = ode15s(@OCFDonecolumnA, [0,5000] ,C0,options,u);

I would like to have the entire time domain split into 5 different sections and have the elements of u take different values at different times(something like a multiple step function). what would be the best way to code it considering that I later want to solve a optimization problem to determine values of u for each time interval in the domain [0,5000].

Thanks

anad
  • 63
  • 5
  • Changing parameters of an ODE discontinuously is usually not a good idea (and violates the mathematical smoothness requirements for ODEs) if you're looking for correct and accurate results – even using a stiff solve like `ode15s`. The way to do this is with separate piecewise integrations in time, stopping and starting the integrator for each. See [my answer here](http://stackoverflow.com/a/17370733/2278029) and let me know if it helps. Also note how I pass parameters via an anonymous function rather than the deprecated method of adding them after the options input. – horchler Feb 03 '15 at 00:51

1 Answers1

3

ode15s expects your model function to accept a parameter t for time and a vector x (or C as you named it) of differential state variables. Since ode15s will not let us pass in another vector u of control inputs, I usually wrap the ODE function ffcn inside a model function:

function [t, x] = model(x0, t_seg, u)
    idx_seg = 0;
    function y = ffcn(t, x)
      % simple example of exponential growth
      y(1) = u(idx_seg) * x(1)
    end
    ...
end

In the above, the parameters of model are the initial state values x0, the vector of switching times t_seg, and the vector u of control input values in different segments. You can see that idx_seg is visible from within ffcn. This allows us to integrate the model over all segments (replace ... in the above with the following):

t_start = 0;
t = t_start;
x = x0;
while idx_seg < length(t_seg)
    idx_seg = idx_seg + 1;
    t_end = t_seg(idx_seg);
    [t_sol, x_sol] = ode15s(@ffcn, [t_start, t_end], x(end, :));
    t = [t; t_sol(2 : end)];
    x = [x; x_sol(2 : end, :)];
    t_start = t_end;
end

In the first iteration of the loop, t_start is 0 and t_end is the first switching point. We use ode15s now to only integrate over this interval, and our ffcn evaluates the ODE with u(1). The solution y_sol and the corresponding time points t_sol are appended to our overall solution (t, x).

For the next iteration, we set t_start to the end of the current segment and set the new t_end to the next switching point. You can also see that the last element of t_seg must be the time at which the simulation ends. Importantly, we pass to ode15s the current tail y(end, :) of the simulated trajectory as the initial state vector of the next segment.

In summary, the function model will use ode15s to simulate the model segment by segment and return the overall trajectory y and its time points t. You could convince yourself with an example like

[t, x] = model1(1, [4, 6, 12], [0.4, -0.7, 0.3]);
plot(t, x);

which for my exponential growth example should yield

Simulation Result

For your optimization runs, you will need to write an objective function. This objective function can pass u to model, and then calculate the merit associated with u by looking at x.

s.bandara
  • 5,636
  • 1
  • 21
  • 36