Forward Diff Trying to Cast Dual Type to Float64

Specific Domains Optimization (Mathematical)
, , , ,
1

Hello,

I’m new to JuMP and trying to implement a yaw-plane vehicle model. The inputs to the function are arrays of the steering and throttle inputs (for every discrete time step). The model runs fine on its own, but adding JuMP has introduced dual numbers. It appears as though every time I try allowing the Dual type to propagate throughout my function, it still tries to cast it as a Float64 value. As I understand, this is a relatively common issue with JuMP newbies.

My question: how can I pass the steering angles and throttle inputs to my function without casting them as Float64?

Code:

#From what I've found online, it helps to define the arguments as a splat of #type Real
function opt_function(control::Type...) where {Type :< Real}
#INITIAL PARAMETERS
    x0 = 0 #metres
    y0 = -3 #metres
    theta0 = 0 #rad
    u0 = 20 #m/s
    v0 = 0 #m/s
    r0 = 0 #m/s
    X0 = [x0, y0, theta0, u0, v0, r0]

########################################
#PROBLEM SECTION
#Trying to store the dual number values passed in the control array in "slip" and "steer"
    slip = Vector{Type}(undef, Int(ceil(T/delT)))
    steer = Vector{Type}(undef, Int(ceil(T/delT)))
   
#Control[] is formatted as a vcat(slip, steer)
    for i in 1:Int(ceil(length(control)/2))
        slip[i] = control[i]
        steer[i] = control[i+Int(length(control)/2)]
    end
#########################################
    #vehicle parameters
    m = 1914 #kg
    Izz = 2600 #kg*m^2
    a = 1.35 #metres
    b = 1.5 #metres
    kx = 100000 #N/m
    cf = 80000#2*1437*dpr #N/m
    cr = 80000#2*1507*dpr #N/m
    params = [m, Izz, a, b, kx, cf, cr]
    
#Calls the vehicle model
    return bicycle_discrete_opt([slip,steer], X0, params, T, delT)
end

###############################################
#Leaving this in case it helps
using JuMP
using Ipopt


global T = 5 # end time {seconds}
global delT = 2.5 #time step {seconds}

model = Model(Ipopt.Optimizer)

#X = x (global), y (global), theta, u ,v ,r (1x6 vector)
#@variable(model, slip_coeffs[1:3])
#@variable(model, steer_coeffs[1:3])
#@variable(model, slip[1:Int(ceil(T/delT))])
#@variable(model, steer[1:Int(ceil(T/delT))])
@variable(model, control[1:2*Int(ceil(T/delT))])

@operator(model, bicycle_op, 2, (control...) -> opt_function(collect(control)))
@objective(model, Min, bicycle(control...))
#optimize!(model)

#############################################
#yaw-plane model, I hope this formats better than it shows in this window

function bicycle_discrete_opt(slip, steer, x0, p, T, delT)
    m, Izz, a, b, kx, cf, cr = p

    x = Array{Real, 2}(undef, Int(ceil(T/delT))+1, 6)
    x[1,:] .= x0

    for i in 1:Int(ceil(T/delT))
        t = (i-1)*delT

        x[i+1,1] = x[i,1] + (x[i,4]*cos(x[i,3])-x[i,5]*sin(x[i,3]))*delT #x_dot (global coordinate frame) #Note the negative sign (x aligns with initial vehicle dir, y is left positive)
        x[i+1,2] = x[i,2] + (x[i,4]*sin(x[i,3])+x[i,5]*cos(x[i,3]))*delT #y_dot (global coordinate frame)
        x[i+1,3] = x[i,3] + (x[i,6]*pi/180)*delT #theta_dot (radians)
        x[i+1,4] = x[i,4] + (2*kx*slip_ratio[i]/m)*delT #u_dot
        x[i+1,5] = x[i,5] + ((cf*steer[i] - (x[i,5]*(cf+cr) + (x[i,6]*pi/180)*(a*cf - b*cr + m*x[i,4]^2))/x[i,4])/m)*delT #v_dot
        x[i+1,6] = x[i,6] + ((180/pi)*(cf*a*steer[i] - (x[i,5]*(a*cf - b*cr) +(x[i,6]*pi/180)*(a^2 * cf + b^2 *cr))/x[i,4])/Izz)*delT #r_dot (deg/s)
    end
2

Edit: bicycle_discrete_opt() returns a scalar from a custom cost function.

3

Hi @Bocconocini, welcome to the forum :smile:

You don’t need to use a custom operator for this. As you have found, it can be a bit difficult.

I would start to write your code like this:

using JuMP
using Ipopt
T = 5        # end time {seconds}
delT = 2.5   # time step {seconds} 
x0 = 0       # metres
y0 = -3      # metres
theta0 = 0   # rad
u0 = 20      # m/s
v0 = 0       # m/s
r0 = 0       # m/s
m = 1914     # kg
Izz = 2600   # kg*m^2
a = 1.35     # metres
b = 1.5      # metres
kx = 100_000 # N/m
cf = 80_000  # 2*1437*dpr # N/m
cr = 80_000  # 2*1507*dpr # N/m
# slip_ratio = ???
N = ceil(Int, T / delT)
model = Model(Ipopt.Optimizer)
@variable(model, slip[1:N])
@variable(model, steer[1:N])
@variable(model, x[1:N + 1, 1:6])
fix.(x[1,:], [x0, y0, theta0, u0, v0, r0])
for i in 1:N
    t = (i - 1) * delT
    @constraints(model, begin
        x[i+1,1] == x[i,1] + (x[i,4]*cos(x[i,3])-x[i,5]*sin(x[i,3]))*delT 
        x[i+1,2] == x[i,2] + (x[i,4]*sin(x[i,3])+x[i,5]*cos(x[i,3]))*delT
        x[i+1,3] == x[i,3] + (x[i,6]*pi/180)*delT #theta_dot (radians)
        x[i+1,4] == x[i,4] + (2*kx*slip_ratio[i]/m)*delT #u_dot
        x[i+1,5] == x[i,5] + ((cf*steer[i] - (x[i,5]*(cf+cr) + (x[i,6]*pi/180)*(a*cf - b*cr + m*x[i,4]^2))/x[i,4])/m)*delT #v_dot
        x[i+1,6] == x[i,6] + ((180/pi)*(cf*a*steer[i] - (x[i,5]*(a*cf - b*cr) +(x[i,6]*pi/180)*(a^2 * cf + b^2 *cr))/x[i,4])/Izz)*delT #r_dot (deg/s)
    end)
end

There are a few bits missing, but it should point you in the right direction.

There are a few optimal control tutorials in the documentation. Check out: Example: nonlinear optimal control of a rocket · JuMP, Example: optimal control for a Space Shuttle reentry trajectory · JuMP

1 Like
4

Hi @odow,

I really appreciate your response. I’m working up from this problem to importing a discretized state-space of a full-vehicle model with separate a non-linear tire model :sweat_smile:. Because I’m not certain that I’ll be able to import that model easily with constraint equations, I’m trying to work with the @operator call.

In the meantime, I’ll definitely give this solution a go and see if I can get it going! If you’ve got any more information on getting @operator to work, I’d greatly appreciate it!

1 Like
5

See Nonlinear Modeling · JuMP

There are a few things in your code to fix.

Instead of

function opt_function(control::Type...) where {Type :< Real}

do

function opt_function(control::Type...) where {Type<:Real}

You’ve called bicycle_discrete_opt([slip,steer], X0, params, T, delT) but defined bicycle_discrete_opt(slip, steer, x0, p, T, delT).

Ignore JuMP for now, and make sure that opt_function works if you call it with some Float64 input arguments.

6

Hi @Odow,

Thank you for pointing out some of my hasty mistakes. I had gone through several iterations while frustrated and tried to quickly revert my code to what I had believed to be the “closest functional iteration”. Obviously, when posting to forums, I should be much more careful in the future. Apologies.

Having taken the weekend to step away and come back, I’ve managed to make some more progress on using the @operator command and I think I’m starting to see some more silly mistakes.

Thank you for your time, I’m going to mark this as solved for the time being.

1 Like
7

No need to be “careful” in future. We don’t mind pointing out issues like this if you have been struggling!

1 Like