How to implement hermite-simpson collocation method?

Hello everyone.

I am trying to do a trajectory optimization using Hermite-Simpson collocation method. I want to know if in the package jump or in other package exist a function call hermite-simpson or if i have to create by myself.

Thank you so much

1 Like

Hi @Coral, welcome to the forum.

Take a look at InfiniteOpt.jl

Otherwise there are a couple of optimal-control related tutorials in the documentation:


1 Like

Hi @odow again! I have been looking in the documentation that you sent and any of them reply to my question. In the examples that you provide they are doing the discretization for every function. My idea is to make a function that you can use in every case to do the discretization, i mean, create a function to which I pass a variable and it returns it to me discretized. Do you have any idea of how to do that? I was trying to do it using the trapezoidal integration method

Hi…here is a script using JuMP for hermite-simpson. You can modify it for your example.

using JuMP
import Ipopt
using GLMakie

dubins = Model(Ipopt.Optimizer)
vm = 1.0   # Max forward speed
um = 1.0   # Max turning speed
ns = 3     # Number of states
nu = 1     # Number of inputs
n = 200    # Time steps

# System dynamics
function f(x, u)
    return [vm*cos(x[3]), vm*sin(x[3]), u[1]]

# Objective Function
function obj_f(x, u, Δt)
    return Δt

# Decision variables
@variables(dubins, begin
    0 <= Δt <= 0.05  # Time step
    # State variables
    x[1:3,1:n]               # Velocity
    # Control variables
    -um ≤ u[1:1,1:n] ≤ um    # Thrust

# Objective
@objective(dubins, Min, obj_f(x[1:3,1:n], u[1,1:n], Δt))
#Initial conditions
@constraint(dubins, x[1:3, 1] .== [0, 0, -π/2 ] )
# Final conditions
@constraint(dubins, x[1:3, n] == [5.0, 5.0, π/2 + π/4])

# Dynamic constraints: hermite_simpson
for j in 1:n-1
    ub = (u[1:nu, j] + u[1:nu, j+1]) / 2

    f1 = f(x[1:ns, j], u[1:nu, j])
    f2 = f(x[1:ns, j+1], u[1:nu, j+1])
    xb = (x[1:ns, j] + x[1:ns, j+1]) / 2 + Δt * (f1 - f2) / 8
    fb = f(xb, ub)

    @constraint(dubins, x[1:ns, j+1] == x[1:ns, j] + Δt * (f1 + 4 * fb + f2) / 6)

# Solve for the control and state

# Display results
println("Min time: ", objective_value(dubins)*n)

fig = Figure()
ax = Axis(fig[1,1], autolimitaspect = 1)
lines!(ax, value.(x[1,1:n]), value.(x[2, 1:n]))

I have an experimental package DirectOptimalControl.jl which allows you to define the problem similar to GPOPS. However, the documentation is limited and there are some limitations with large problems.

1 Like

Thank you so much, i would try to use it for my example!


Hi again! How can you implement hermite-simpson as a function? My goal in my project is to make a function to which I pass the state space and return it directly discretized using hermite-simpson or trapezoidal integration. Basically a function using direct placement methods to be able to do the discretization.

Thank you so much for your help

Also if you can add in a message the information about the problem like the dynamics equations, a definition of the problem it would be really helpful to understand the code

Thank you