How to implement an objective function in Jump

Specific Domains Optimization (Mathematical)
1

Hello everyone,

I am implementing an optimization problem for a satellite that has to move from a radius 1 orbit to a radius 2 orbit. The objective of my optimization is to minimize the required deltaV. I am trying to define this objective function but I get an error saying that the way I define the deltaV is not supported by Jump. How can I define it in another way so that it accepts the function? I attach the code and the error obtained in case you can help me to solve it.

Thanks in advance

The error is:

ERROR: The objective function deltaV is not supported by JuMP.

The implemented code:


using JuMP
import Ipopt
import Plots
using Plots

# Constants and initial conditions

μ = 3.986e5 #[km^3/s^2]
x0 = 1      #[m]
y0 = 0      #[m]
xf = 2      #[m]
yf = 0      #[m]
m = 100     #[kg]
n = 500     #Number of iterations
const integration_rule = "Trapezoidal integration"

# Create the Model

Hohmann = Model(Ipopt.Optimizer)
set_attribute(Hohmann,"print_level", 4) #Da mas informacion sobre la optimizacion 
set_attribute(Hohmann, "max_iter", 20000) #Define el numero maximo de iteraciones del solucionador
set_attribute(Hohmann, "print_timing_statistics", "yes") #Nos proporciona informacion sobre el tiempo que toma cada fase del proceso de optimizacion 

# Create variables

@variables(Hohmann, begin
    x[1:n] 
    y[1:n]
    u[1:n]
    v[1:n]
end)

# Set initial conditions

fix(x[1], x0, force=true)
fix(x[n], xf, force=true)
fix(y[1], y0, force=true)
fix(y[n], xf, force=true)

#Aditional equations

@NLexpression(Hohmann, Δt, 0.2/n) #Time step

@NLexpressions(Hohmann, begin
    T[j=1:n], sqrt(u[j]^2 + v[j]^2)
    R[j=1:n], sqrt(x[j]^2 + y[j]^2)
    tf[j=1:n], Δt*n
end)

# Dynamic equations

@NLexpressions(Hohmann, begin
    δδx[j=1:n], -μ/R[j]*x[j] + u[j]/m   
    δδy[j=1:n], -μ/R[j]*y[j] + v[j]/m   
end)

# Collocation method (trapezoidal integration function)

function CollocationMethods(x, y, Δt)
    if integration_rule == "Rectangular integration"
        for j in 2:n
            @NLconstraint(Hohmann, x[j] == x[j-1] + Δt*y[j-1])
        end     
    elseif integration_rule == "Trapezoidal integration"
        for j in 2:n
            @NLconstraint(Hohmann, x[j] == x[j-1] + 0.5*Δt*(y[j]+y[j-1]))
        end  
    else
        @error "Unexpected integration rule '$(integration_rule)'"
    end
end

# Objective function 

function deltaV(T,tf)
    term = 0
    for i in 2:n-1
        term = term + 2*(T[i]/m)
    end

    h = tf/n
    deltaV = h/2*(T[0]/m + T[n]/m + term)
end

register(Hohmann, :deltaV, 2, deltaV; autodiff = true)

# Use collocation method

CollocationMethods(x, δδx, Δt)
CollocationMethods(y, δδy, Δt)

@objective(Hohmann, Min, deltaV)

optimize!(Hohmann)

println( "The minimun deltaV is: ", objective_value(Hohmann))
2

Not sure if that’s the only problem, but you are still using the old version of the nonlinear modeling API. Take a look at Nonlinear Modeling · JuMP for the current way of doing things

3

You have a few issues:

  • As @gdalle suggests, you can update to the new JuMP syntax, which basically means dropping the @NL and removing the register
  • Your objective expressions is a function deltaV. You need to actually call it, so @objective(Hohmann, Min, deltaV(T, tf))
  • Your deltaV function calls T[0]. This index does not exist. Did you mean T[1]?
  • You need to think carefully about variable bounds, and where your functions are defined. The initial point has u = 0 and v = 0. The derivative at sqrt(0) is not defined, so Ipopt will throw an error. You need to ensure that all functions are twice differentiable, so ensure that, for example, R[j] cannot be 0.

An improvement, but still not quite correct, is the following. You could improve things by putting tight bounds on x, y, u and v.

using JuMP
import Ipopt
function main()
    μ = 3.986e5
    x0 = 1     
    y0 = 0     
    xf = 2     
    yf = 0     
    m = 100    
    n = 500    
    Δt = 0.2 / n
    tf = Δt * n
    integration_rule = "Trapezoidal integration"
    Hohmann = Model(Ipopt.Optimizer)
    @variables(Hohmann, begin
        x[1:n]
        y[1:n]
        u[1:n]
        v[1:n]
        uv[1:n] >= 0.00001
        xy[1:n] >= 0.00001
        T[1:n] >= 0.00001
        R[1:n] >= 0.00001
        δδx[1:n]
        δδy[1:n]
    end)
    fix(x[1], x0)
    fix(x[n], xf)
    fix(y[1], y0)
    fix(y[n], xf)
    @constraints(Hohmann, begin
        [j=1:n], uv[j] == u[j]^2 + v[j]^2
        [j=1:n], xy[j] == x[j]^2 + y[j]^2
        [j=1:n], T[j] == sqrt(uv[j])
        [j=1:n], R[j] == sqrt(xy[j])
        [j=1:n], δδx[j] == -μ / R[j] * x[j] + u[j] / m   
        [j=1:n], δδy[j] == -μ / R[j] * y[j] + v[j] / m   
    end)
    function CollocationMethods(x, y, Δt)
        if integration_rule == "Rectangular integration"
            @constraint(Hohmann, [j in 2:n], x[j] == x[j-1] + Δt*y[j-1])
        elseif integration_rule == "Trapezoidal integration"
            @constraint(Hohmann, [j in 2:n], x[j] == x[j-1] + 0.5*Δt*(y[j]+y[j-1]))
        else
            @error "Unexpected integration rule '$(integration_rule)'"
        end
    end
    function deltaV(T, tf)
        term = 0
        for i in 2:n-1
            term = term + 2 * (T[i] / m)
        end
        h = tf / n
        return h / 2 * (T[1] / m + T[n] / m + term)
    end
    CollocationMethods(x, δδx, Δt)
    CollocationMethods(y, δδy, Δt)
    @objective(Hohmann, Min, deltaV(T, tf))
    optimize!(Hohmann)
    println( "The minimun deltaV is: ", objective_value(Hohmann))
end
2 Likes
4

Thank you so much for your reply! It was really helpful

1 Like
5

Hello again! I tried to do the following implementation reducing my equation to a first order equations. I did it following the same example that you send to me but i am not able to make it work. Do you see why is not working??


using JuMP
import Ipopt

μ = 3.986e5 #[km^3/s^2]
x0 = 6858   #[km]
y0 = 0      #[km]
xf = 27178  #[km]
yf = 0      #[m]
m = 100     #[kg]
n = 500     #Number of iterations
Δt = 0.2/n  #Step time
tf = 2925.5078 # Time of flight[s]
const integration_rule = "Trapezoidal integration"

Hohmann = Model(Ipopt.Optimizer)
@variables(Hohmann, begin
    x[1:n]
    y[1:n]
    Tx[1:n]
    Ty[1:n]
    uv[1:n] >= 0.00001
    xy[1:n] >= 0.00001
    T[1:n] >= 0.00001
    R[1:n] >= 0.00001
    δx[1:n]
    δy[1:n]
    u[1:n]               #For reducction to first order equations
    v[1:n]               #For reducction to first order equations
    δu[1:n]
    δv[1:n]
end)

fix(x[1], x0)
fix(x[n], xf)
fix(y[1], y0)
fix(y[n], xf)

@constraints(Hohmann, begin
    [j=1:n], T[j] = sqrt(u[j]^2 + v[j]^2)
    [j=1:n], R[j] = sqrt(x[j]^2 + y[j]^2)
    [j=1:n], δx = u[j]
    [j=1:n], δy = s[j]
    [j=1:n], δu = -μ/R[j]*x[j] + Tx[j]/m
    [j=1:n], δv = -μ/R[j]*y[j] + Ty[j]/m
end)

@constraints(Hohmann, begin
    [j=1:n], uv[j] == Tx[j]^2 + Ty[j]^2
    [j=1:n], xy[j] == x[j]^2 + y[j]^2
    [j=1:n], T[j] == sqrt(uv[j])
    [j=1:n], R[j] == sqrt(xy[j])
    [j=1:n], δx = u[j]
    [j=1:n], δy = s[j]
    [j=1:n], δu[j] == -μ / R[j] * x[j] + u[j] / m   
    [j=1:n], δv[j] == -μ / R[j] * y[j] + v[j] / m   
end)

function CollocationMethods(x, y, Δt)
    if integration_rule == "Rectangular integration"
        @constraint(Hohmann, [j in 2:n], x[j] == x[j-1] + Δt*y[j-1])
    elseif integration_rule == "Trapezoidal integration"
        @constraint(Hohmann, [j in 2:n], x[j] == x[j-1] + 0.5*Δt*(y[j]+y[j-1]))
    else
        @error "Unexpected integration rule '$(integration_rule)'"
    end
end

function deltaV(T, tf)
    term = 0
    for i in 2:n-1
        term = term + 2 * (T[i] / m)
    end
    h = tf / n
    return h / 2 * (T[1] / m + T[n] / m + term)
end
CollocationMethods(x, δx, Δt)
CollocationMethods(y, δy, Δt)
@objective(Hohmann, Min, deltaV(T, tf))
optimize!(Hohmann)
println( "The minimun deltaV is: ", objective_value(Hohmann))