How to implement hermite-simpson collocation method?

Specific Domains Optimization (Mathematical)
, ,
1

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
2

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:

3

Thanks!!

1 Like
4

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

5

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

Hermit-Simpson 
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]]
end

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

# 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
end)

# 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)
end

# Solve for the control and state
println("Solving...")
JuMP.optimize!(dubins)
solution_summary(dubins)

# 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]))
fig

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
6

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

2 Likes
7

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

8

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