How to interface symbolic solutions with numeric ode solvers

Hello Community,

My question is of an optimisation nature, and as the title suggests, I have been having difficulty getting results in acceptable time when finding the solution to a dynamical system using the symbolic toolbox and then applying it to an ode solver.

I have copied and pasted the equations of motion, in symbolic form, then substituted parameters, and then solved for the second derivatives and stored this as 'solution', i.e. solution is a symbolic structure which contains the solutions to the second derivatives for the physical system.

When I use the function 'subs' on my solution variable within my derivative function that gets passed to ode45 the issue is that the software spends all it's time within the environment of the symbolic toolbox and getting a time series for 50 seconds takes >20 minutes, which is ridiculous (confirmed with the profiler).

As a work around I have defined a structure consisting of functions that have been converted from the solutions into function handles, and I pass this structure to the ode solver so that the second derivatives can be found outside of the symbolic toolbox and this has made the simulation at least feasible (in the order of 3mins). However when I look at the profiler I can see that the functions themselves are using the symbolic toolbox, even though they contain no symbols - specifically they are using symengine. They have no reason to use the symengine since they have no symbolic variables (all their symbolic variables are converted to function variables with matlabFunction)

eom = subs(euler_equations+generalised_dissipation-generalised_drag,parameters);

sol = solve(eom, [x_dd phi_dd theta_dd psi_dd]);

solution.x_dd = matlabFunction(sol.x_dd); % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)
solution.phi_dd = matlabFunction(sol.phi_dd); % @(Qdrag_phi,phi_d,theta_d)Qdrag_phi-phi_d+theta_d
solution.theta_dd = matlabFunction(sol.theta_dd); % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)
solution.psi_dd = matlabFunction(sol.psi_dd); % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)

Where euler_equations, generalised_dissipation and generalised_drag are symbolic variables that describe the system and parameters is a structure containing the values of the parameters.

The symbolic toolbox is still eating up ~1/3 of the simulation time, and I would like to know how to define my functions so that they wont under any circumstance use the symbolic toolbox. I understand that I am have most likely used bad MatLab practise when interfacing between symbolic toolbox and solutions that are applicable to ode45 and I am open to suggestions on how to more optimally achieve my aims.

The code that leads to the simulation is as follows:

X0 = [0;30;0;0;0;0;pi;0];
opts = odeset('RelTol',1e-3,'AbsTol',1e-6,'MaxStep',1e-3);
[t, X] = ode45(@(t,X) unicycleODEs(t, X, solution, inf_drag, parameters), [0,50], X0, opts);

function Xdot = unicycleODEs(t, X, X_dd, inf_drag, parameters)
    % Extract state variables from X
    x = X(1);
    x_d = X(2);
    phi = X(3);
    phi_d = X(4);
    theta = X(5);
    theta_d = X(6);
    psi_ = X(7);
    psi_d = X(8);

    % find generalised force of drag
    inf_drag_funct = @(chi) inf_drag(chi,psi_,psi_d,theta,theta_d,x_d);
    generalised_drag = integral(inf_drag_funct,0,parameters.l,'ArrayValued',true);
    % grab components
    Qdrag_x = generalised_drag(1);
    Qdrag_phi = generalised_drag(2);
    Qdrag_theta = generalised_drag(3);
    Qdrag_psi = generalised_drag(4);

    % get Xdd

    % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)
    x_dd = X_dd.x_dd(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d); 
    % @(Qdrag_phi,phi_d,theta_d)
    phi_dd = X_dd.phi_dd(Qdrag_phi,phi_d,theta_d);
     % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)
    theta_dd = X_dd.theta_dd(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d);
    % @(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d)
    psi_dd = X_dd.psi_dd(Qdrag_x,Qdrag_psi,Qdrag_theta,phi_d,psi_,psi_d,theta,theta_d); 

    % Assign the derivatives of each state variable to Xdot
    Xdot = [x_d; x_dd; phi_d; phi_dd; theta_d; theta_dd; psi_d; psi_dd];
end

For completness, the most time is spent doing a numeric integration of inf_drag at every time step of the derivative function.

In summary, I have gotten the simulation to work but I know that the way I've done it would go against wider wisdom and convention, and this shows in the simulation times. How can I optimise the process, why does symbolic toolbox completely kill simulation times whenever it is found within the derivative function loop?

The system itself is not that complicated and should not take this long to simulate 50 seconds - once this project is done I will likely want to simulate this system over a time span of ~1000, which is not feasible as it stands.

Thank you in advance,

Should you require the full code that includes the equations of motion I would be happy to provide it.

Answer 1

[Star Strider]

With an initially symbolic set of differential equations, the symbolic versions should not be part of a numerical integration. The usual approach to creating linear or nonlinear differential equations using the Symbolic Math Toolbox to set them up and then using that result to numerically integrate them is something like this —

syms x(t) t Y 
Dx  = diff(x);
D2x = diff(Dx);
DE = D2x + Dx + x == exp(0.01*x) * sin(2*pi*x)

DE(t) = 

[VF,Sbs] = odeToVectorField(DE)

VF = 

Sbs = 

odefcn = matlabFunction(VF, 'Vars',{t,Y})

odefcn = function_handle with value:
    @(t,Y)[Y(2);exp(Y(1)./1.0e+2).*sin(pi.*Y(1).*2.0)-Y(1)-Y(2)]

[t,y] = ode45(odefcn, [0 15], [0 1]);

figure
plot(t, y)
grid
xlabel('t')
ylabel('Value')
legend(string(Sbs), 'Location','best')

.

[Matthew Ediz Beadman]

Thank you for your answer; I will keep in mind the procedure, although in my case it is not too helpful. Put simply there is a term within the differential equation which is an integration that is too complicated to do symbolically, and so instead I replace it's vector values with dummy variables ([Qdrag_x; Qdrag_y] for example), and so the solution of the differential equations is a function of the state variables and [Qdrag_x; Qdrag_y].

From here I find what Q_dragx and Qdrag_y is at each time step of the numerical solver via a numerical integration. Thinking about it, I could just define what the numerical integration is in terms of the state variables (i.e. a numerical integration of symbolic variables) beforehand and that should optimise simulation time.

edit: clarification; the D.E are not known as a function of the state variables a priori, and only become known at each step of the ode solver after a numerical integration is conducted which uses the state variables of the current time step.

[Star Strider]

If you are integrating a system of differential equations to fit to data, that is relatively straightforward, however it can be difficult with complicated nonlinear systems and noisy data.

You must be able to get your code in a form that the MATLAB ODE integrators can work with (systems of first-order differential equations) to do this, and the integrated result must produce a column vector or column-major matrix that you can compare with your data.

[Matthew Ediz Beadman]

It is an integration of an equation. At the expense of getting technical, I will explain. It is the 3 dimensional drag equation which I am integrating over the surface of a moving body (assuming drag acts at the CoM is inaccurate).

The equation depends on the state variables, which are being solved for by ode45, and so each iteration I simply run a numerical integration with the new values of the state variables.

[Star Strider]

The ‘unicycleODEs’ function, when integrated, will produce an (Nx8) matrix (where ‘N’ is the length of the time or other independent variable vector the deriviatives are with respect to). It is relatively straightforward to select any 3 columns to compare with your data (that I assume describe the dimensions).

If your posted code already runs correctly, and you have data you want to fit your system of differential equaations to, I woulld need for you to attach the data, decide which three columns of the ode45 integrated result ‘X’ you want to fit to it, and what parameters you want to estimate. The rest of what is necessary to estimate the parameters I have already done a few times, so adapting that approach to your code and data should be relatively straightforward.

[Matthew Ediz Beadman]

There is no data to fit nor to compare, it is a simulation of dynamics which I can only guess to be correct based on how the dynamics evolve and plots of the variables in phase space and time series.

% find generalised force of drag
    inf_drag_funct = @(chi) inf_drag(chi,psi_,psi_d,theta,theta_d,x_d);
    generalised_drag = integral(inf_drag_funct,0,parameters.l,'ArrayValued',true);
    % grab components
    Qdrag_x = generalised_drag(1);
    Qdrag_phi = generalised_drag(2);
    Qdrag_theta = generalised_drag(3);
    Qdrag_psi = generalised_drag(4);

The above is the fundamental issue which causes the slow simulation.

The inf_drag_func is a function of the state variables that are evolving as the function is being integrated by ode45, and in turn affects the evolution of the state variables, so its a feedback loop - I can't see a way to let the dynamics run without drag to then calculate drag afterwards and apply them back to the dynamics.

[Star Strider]

I needed to clarify if you’re fitting data, since what you want to do it isn’t obvious to me.

You’re running the simulation with known parameters to simulate the sytem, then? You’re not optimising the parameters?

[Matthew Ediz Beadman]

No, I'm not.

[Star Strider]

With respect to this and your other question, have we resolved this, or is there more to do?

[Matthew Ediz Beadman]

You solved this, I only replied as I thought you had something to add. Thank you :)

[Star Strider]

As always, my pleasure!

Answer 2

[Walter Roberson]

When you symfun(), an anonymous function is generated that does not use the symbolic toolbox.

The traceback for the anonymous function shows it to have been generated inside the symbolic toolbox, but that does not mean that it uses the symbolic toolbox in any way.

syms x
f(x) = x^2

f(x) = 

h = matlabFunction(f)

h = function_handle with value:
    @(x)x.^2

functions(h)

ans = struct with fields:
            function: '@(x)x.^2'
                type: 'anonymous'
                file: '/MATLAB/toolbox/symbolic/symbolic/symengine.p'
           workspace: {[1×1 struct]}
    within_file_path: ''

Generated inside symengine.p does not mean that it relies on the symbolic toolbox to execute.

[Matthew Ediz Beadman]

Hi, thanks for your answer. I'd just like to clarify to other readers that I have tested this by defining all the functions that are defined within the program as their own functions (without defining symbolic variables) in a seperate script and timing the execution of both scripts.

Both scripts take the same time so you are indeed correct, thank you.

Both replies to this post have answered different questions I had, however Star Strider's response was more in line the question title and so I will accept his answer. I am grateful for your answer regardless :)