Performance Issues with Symbolically Generated Functions for use with DifferentialEquations

Hello anyone reading this, sorry if there are some formatting issues. First post and all. Promise I’ll learn. Anyways, a lot of the work I do involves optimal control for spacecraft using indirect methods, so I would like to solve differential equations fast. This is obviously the major draw of Julia for me. My experience with Julia in this department has been pretty good so far, I have learned a lot from SciMLTutorials.jl and from solving several simple problems on my own. Recently, I wrote the following code.

using LinearAlgebra
using Symbolics
function generate_eoms()
  @variables p f g h k L m           # state variables
  @variables p1 p2 p3 p4 p5 p6 p7    # costates
  @variables ur ut un delta          # controls
  @variables c Thr rho SI2CAN t      # extra parameters
  pd = c, Thr, rho, SI2CAN

  # Form state vectors
  S = [p;f;g;h;k;L;m]
  P = [p1;p2;p3;p4;p5;p6;p7]
  X = [S; P]

  # Get A matrix
  SinL = sin(L)
  CosL = cos(L)
  q = 1+f*CosL+g*SinL
  s = 1+h^2+k^2
  C1 = sqrt(p)
  C2 = 1/q
  C3 = h*SinL-k*CosL
  A = [       0.0                    2*p*C2*C1                   0.0;
          C1*SinL         C1*C2*((q+1)*CosL+f)          -C1*(g/q)*C3;
         -C1*CosL         C1*C2*((q+1)*SinL+g)           C1*(f/q)*C3;
              0.0                          0.0        C1*s*CosL*C2/2;
              0.0                          0.0        C1*s*SinL*C2/2;
              0.0                          0.0              C1*C2*C3;
              0.0                          0.0                   0.0]

  # Get vector optimal control 
  Pv = transpose(A)*P
  normPv = sqrt(Pv[1]*Pv[1] + Pv[2]*Pv[2] + Pv[3]*Pv[3])
  U_op = -Pv/normPv

  # Get scalar optimal control
  SF = c*SI2CAN*normPv/m+p7
  delta_op = 0.5*(1+tanh(SF/rho))
  Delta = Thr/m*SI2CAN*delta_op*U_op

  # Get acceleration terms
  b = [0.0; 0.0; 0.0; 0.0; 0.0; sqrt(p)*(q/p)^2; -Thr/c*delta_op]
  dy = A*Delta+b
  
  # Define Cost
  Lagrangian = 0.0
  
  # Hamiltonian is a scalar
  Hamiltonian = Lagrangian + dot(P,dy)
  
  # Here, we define the costates as the negative gradient of the Hamiltonian 
  costates = -Symbolics.gradient(Hamiltonian, S)
  
  # Full nonlinear system of equations here
  F = [dy;costates]
  
  return eval(build_function(F, X, pd, t)[2])
end

Since the research I do involves solving different problems for lots of complex dynamical systems, it’s nice to be able to symbolically derive the governing equations since half of them are nearly impossible to derive by hand (the “costates,” which are derived identically to the time rate of change of momenta in Hamiltonian mechanics).

The goal of the function that I wrote is to use the physical dynamics to derive the costate dynamics and then generate an optimized function to calculate both that can be integrated with DifferentialEquations.jl. However, the code shown above produces a function that takes an incredibly long time to evaluate just once (on the order of minutes), let alone integrate, so I am obviously doing something wrong here. The strange thing is, I have used this exact process in Julia for different dynamical systems and it has worked very well! But not this time.

I have attempted to take a different approach using the ForwardDiff.jl package like so:

using ForwardDiff
xDot!(dX,ham,X,P,probdata,t) = ForwardDiff.jacobian!(dX, P-> hamiltonian(ham,X,P,probdata), P)
pDot!(dP,ham,X,P,probdata,t) = ForwardDiff.jacobian!(dP, X->-hamiltonian(ham,X,P,probdata), X)

The Hamiltonian function is nearly identical to generate_eoms(), but just calculates the value of the Hamiltonian. This approach works better (and is a lot nicer to look at), but the resulting functions are still pretty slow. The only thing I have been able to do to speed up the evaluation of the generated function is to avoid some of the substitutions, write the code to its own file, and then manually add the lines that I avoided substituting to the new, written file. This is obviously not ideal, but it is very fast. The code written by build_function() is extremely lengthy otherwise.

My main questions: What is the best way to programmatically generate the differential equations of a Hamiltonian system in an efficient form and then integrate them? Where did I go wrong here?

Please leave any and all criticisms of my code/questions/person. Please leave them as soon as you have a chance, my advisor is breathing down my neck.
Thanks! - Sam

Which version of Julia are you using? On my Windows laptop, running Julia v1.7.2, it executes in 3.6 seconds. If you are using an old version of Julia maybe upgrade and also update all your packages.

1 Like

I also see a runtime of a few seconds so I’m confused. Maybe I’m missing something. Example of the code that takes longer?

Apologies, I should have been more clear. Evaluating the function produced by this code with some argument does not execute in a reasonable amount of time. I am currently running Julia 1.7.2 on a Windows Surface Book 2, and I have been using the package manager to keep my packages up to date almost every day. For example, the following code,

f = generate_eoms()
X = rand(14,1)
dX = rand(14,1)
pd = 1.0, 2.0, 3.0, 4.0
f(dX,X,pd,0)

takes a significant amount of time to run. The generate_eoms() function itself takes only a few seconds to produce f . Please let me know if you encounter the same issue on your machines running this code.

Your function is wild, Wild function generation · GitHub, hundred thousand characters or so. That’s why it takes so long to compile. Is this expected?

Note that A/b is a linear solve. If those are supposed to be element-wise, then they need a dot: A ./ b.

julia> b = [0.5 0.2]
1×2 Matrix{Float64}:
 0.5  0.2

julia> A = [1 1
            1 1]
2×2 Matrix{Int64}:
 1  1
 1  1

julia> A / b
2×1 Matrix{Float64}:
 2.413793103448275
 2.413793103448275

julia> A ./ b
2×2 Matrix{Float64}:
 2.0  5.0
 2.0  5.0

If you’re doing A / b and that’s what you mean, then of course it needs to calculate the lu-factorization or inverse of A, which is what grows the equation size and the symbolic expression.

If you’re looking to do this over inversions, then you will need to use automatic differentiation because of the expression growth.

2 Likes

Ah, that makes sense. It is indeed supposed to be element-wise. I updated my code to reflect this, but the size of the output does not change, still around a hundred thousand characters. I may still be approaching this wrong. This size output was not expected. I have previously used this work flow in MATLAB (I know), which produces much more compact output, but the rules in Julia are very different for a lot of good reasons. I decided to move on to the Automatic Differentiation route, as it ends up being easier and evaluating much faster, but thanks so much for your advice!

I ran into another problem (of course). If you would like me to open up another thread about it I can, but I will show the code that I am curious about below.

using ForwardDiff
using DifferentialEquations

function hamiltonian(X,P,probdata)
  # This function defines the Hamiltonian for the system
  # Extract problem data
  c, Thr, rho, SI2CAN = probdata

  # Extract states (MEEs)
  p = X[1]
  f = X[2]
  g = X[3]
  h = X[4]
  k = X[5] 
  L = X[6]
  m = X[7]

  # Get A matrix
  SinL = @. sin(L)
  CosL = @. cos(L)
  q = @. 1+f*CosL+g*SinL
  s = @. 1+h^2+k^2
  C1 = @. sqrt(p)
  C2 = @. 1/q
  C3 = @. h*SinL-k*CosL

  A = @. [       0            2*p*C2*C1               0;
           C1*SinL C1*C2*((q+1)*CosL+f)    -C1*(g/q)*C3;
          -C1*CosL C1*C2*((q+1)*SinL+g)     C1*(f/q)*C3;
                 0                    0  C1*s*CosL*C2/2;
                 0                    0  C1*s*SinL*C2/2;
                 0                    0        C1*C2*C3;
                 0                    0               0]

  # Get vector optimal control
  Pv = transpose(A)*P
  normPv = sqrt(dot(Pv,Pv))
  U_op = @. -Pv/normPv

  # Get scalar optimal control
  SF = @. c*SI2CAN*normPv/m + P[7]
  delta_op = @. 0.5*(1+tanh(SF/rho))
  Delta = @. Thr/m*SI2CAN*delta_op*U_op

  # Define b matrix
  b = @. [0.0; 0.0; 0.0; 0.0; 0.0; sqrt(p)*(q/p)^2; -Thr/c*delta_op]

  # Return Hamiltonian, which is just H = Lagrangian + transpose(P)*dy
  return transpose(P)*(A*Delta + b)
end

xDot!(dX,X,P,probdata,t) = ForwardDiff.jacobian!(dX, P-> hamiltonian(X,P,probdata), P)
pDot!(dP,X,P,probdata,t) = ForwardDiff.jacobian!(dP, X->-hamiltonian(X,P,probdata), X)

# Constants
muS = 132712440018 # Gravitational constant for Earth (km^3/s^2)
AU = 149.6e6      # Astronomical unit (km)
g0 = 9.80665     # m/s^2
dtr = pi/180

DU = AU
TU = DU^1.5/sqrt(muS)
VU = DU/TU
SI2CAN = TU^2/(DU*1000)

TOF_TU = 3534*86400/TU # sec;
tspan = [0 TOF_TU]
N_rev = 5

m_i = 4000  # initial mass (kg)
Isp = 3000  # specific impulse (sec)

Thr = 0.34 # N
c = Isp*g0/TU

# Initial and final states
S_i = 1e3.*[0.000999694371880; -0.000003766790161; 0.000016286899122; -0.000000007702061; 0.000000000618817; 0.014161892533821; 4.000000000000000]

S_f = [1.553697085372486; 0.153029069608839; -0.519948174200711; 0.016183102238719; 0.118139527451067; 46.330240341961556]

P0 = [-1388.8509555559442; -229.94430250015753; 673.2433383721474; -884.7938842474259; -2169.4582779313832; 1.5104219965599175; -0.6098385196774864]

# Define initial ode problem
pd = c, Thr, 0.5, SI2CAN
_prob = DynamicalODEProblem(xDot!, pDot!, reshape(S_i,7,1), zeros(7,1), (0,TOF_TU), pd)

# Remake ode problem with new 
prob = remake(_prob; u0 = ArrayPartition((reshape(S_i, 7,1), reshape(P0,7,1))), p = pd)
sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)

Now I am using AD to generate the dynamics and costate dynamics. However, to solve the problem, I need to feed this to a nonlinear solver to constrain the optimization. This involves redefinition of some of the initial conditions (the costates) between each iteration of the nonlinear solver. To do this efficiently, I am using the remake function provided by DifferentialEquations.jl. This is usually a straightforward operation with a regular ODEProblem, just use remake and set your new u0 = newIC. However, this is the only way I can redefine a DynamicalODEProblem and still have solve work.

# Initial ODE problem with initial costates = zeros(7,1)
_prob = DynamicalODEProblem(xDot!, pDot!, reshape(S_i,7,1), zeros(7,1), (0,TOF_TU), pd)
# New ODE problem with initial costates = P0
prob = remake(_prob; u0 = ArrayPartition((reshape(S_i, 7,1), reshape(P0,7,1))), p = pd)

Is this intended? If so, it is pretty confusing for newcomers like me, and there seems to be a lot of type inflexibility when remaking a DynamicalODEProblem compared to remaking an ODEProblem. If it is intended, then I will just make sure to leave comments in my code detailing this process. Thanks again for your help!

Oh this is worth an issue. We need to specialize remake for the DynamicalODEProblem case. That’s worth an issue in DifferentialEquations.jl.