Method to speedup loop with ode45 inside it

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Method to speedup loop with ode45 inside it

octavecontrib
Simulation with a loop and ode 45 inside it executed much slower in Octave (5.2 and 6.x) compared to Matlab. Before optimising the differential equation I was integrating, I decided to profile the code to see if my differential equation code or the discrete-time-controller code was the bottleneck.

Profiling indicated that most of the execution time is spent inside functions other than the differential equation functions specific to my problem.
@inputParser_xxxxx and runge_kutta_45_dorpri are consuming most of the time.

I am simulating a system which has both discrete time sub-system (digital control) and continuous time system (plant being controlled). The fixed-time-step part is implemented by the while loop and the continuous time part is implemented by the ode45.

Is there a way to speed up my simulation ? If the functions at the top of the profile were user written, I could have attempted to convert them to C++ mkoct files.

Below is the code and the profiling result. The code is a minimum working example to illustrate the problem. Even for much more complicated (non linear, time varying) differential equations, I am still getting similar results.

CODE
=========================
clc;
clear;

fprintf('example profiled running...\n');

endtime = 1000;
stepsize = 0.5; % the fixed step used by the controller

t = 0;
ii = 1;
% reserve space so that it doesn't slow down the program
tl = zeros(ceil(endtime/stepsize), 1);
yl = zeros(ceil(endtime/stepsize), 2);

ic = [1, 0];

control = 0.1;

profile on;

while(t < endtime)
% run the non linear time varying differential equation
% using a variable time step method.
% example shown here is a linear system.
% but profiling results are the same.
[tt, yy] = ode45(@(t,x)[0,1; -1, control]*x/2, [0, stepsize], ic);

% collect results for plotting later
yl(ii, :) = yy(end, :);
tl(ii) = t + tt(end, :);

ic = yy(end, :);

if(norm(ic) > 5)
  % this is where the controller logic would have been written in the real application
  control = -0.05;
end

ii = ii+1;

t = t+stepsize;
end

profile off;
profshow;

plot(tl, yl);

fprintf('example profiled ran...\n');

Profiling result
==============================
   #                  Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------------------
  75     runge_kutta_45_dorpri            14.695      15.45        67522
  23        @inputParser/parse            14.132      14.86        11200
  73        integrate_adaptive            13.456      14.15         5600
  29  @inputParser/add_missing             7.984       8.40        22400
  55   @inputParser/is_argname             6.316       6.64       128800

  76                     feval             4.414       4.64       410732
  57            anonymous@:0:0             3.895       4.10       550732
  56 @inputParser/validate_arg             3.865       4.06       128800
  63               AbsRel_norm             3.263       3.43        84322
  35                    unique             2.574       2.71        28000
  30                   setdiff             2.478       2.61        22400
   2                     ode45             2.419       2.54         5600
  24                fieldnames             1.296       1.36        84000
  60              odemergeopts             1.146       1.20         5600
  62         starting_stepsize             1.064       1.12         5600
   5                    odeset             1.040       1.09        11200
  31              validsetargs             1.032       1.08        22400
  27            __fieldnames__             0.934       0.98        84000
  52                  binary *             0.827       0.87      1225452
  41                      sort             0.782       0.82        22400