JuMP operators for > and mod2pi?

Specific Domains Optimization (Mathematical)
1

I’m trying to solve an NLP with JuMP and the problem that I’m trying to solve canonically converts orbital state vectors into keplerian elements in order to match terminal conditions. This is horrible and very nonlinear, and there are better approaches, but for this problem everyone does it this way. I’m also trying to reuse a function that I have already written, which AFAIK handles all the edge cases correctly. Remarkably, most of the function works, until it gets down to some edge cases for circular equatorial orbits and handling the mod2pi wrapping of the angles.

The function is this:

function rv2oe(μ, r, v)
    rmag = sqrt(sum(i^2 for i in r))         # Magnitude of position vector
    vmag = sqrt(sum(i^2 for i in v))         # Magnitude of velocity vector
    rhat = r ./ rmag                         # Unit vector in the direction of r
    hv = cross(r, v)                         # Specific angular momentum vector
    hhat = hv ./ sqrt(sum(i^2 for i in hv))  # Unit vector of angular momentum
    eccvec = cross(v / μ, hv) - rhat         # Eccentricity vector
    sma = 1.0 / (2.0 / rmag - vmag^2 / μ)    # Semi-major axis
    l = sum(i^2 for i in hv) / μ             # Semi-latus rectum

    # Parameters for frame transformation
    d = 1.0 + hhat[3]
    p = d == 0 ? 0 : hhat[1] / d
    q = d == 0 ? 0 : -hhat[2] / d
    const1 = 1.0 / (1.0 + p^2 + q^2)

    fhat = [
        const1 * (1.0 - p^2 + q^2),
        const1 * 2.0 * p * q,
        -const1 * 2.0 * p
    ]

    ghat = [
        const1 * 2.0 * p * q,
        const1 * (1.0 + p^2 - q^2),
        const1 * 2.0 * q
    ]

    # Calculate Keplerian elements
    h = dot(eccvec, ghat)
    xk = dot(eccvec, fhat)
    x1 = dot(r, fhat)
    y1 = dot(r, ghat)
    xlambdot = atan(y1, x1)                       # True longitude
    ecc = sqrt(h^2 + xk^2)                        # Eccentricity
    inc = 2.0 * atan(sqrt(p^2 + q^2))             # Inclination
    lan = inc > eps() ? atan(p, q) : 0.0          # Longitude of ascending node
    argp = ecc > eps() ? atan(h, xk) - lan : 0.0  # Argument of periapsis
    nu = xlambdot - lan - argp                    # True anomaly

    # Normalize angles to [0, 2π]
    lan = mod2pi(lan)
    argp = mod2pi(argp)
    nu = mod2pi(nu)

    return sma, ecc, inc, lan, argp, nu, l
end

The way I’m trying to use it just to get it to “compile” is:

include("rv2oe.jl")
using LinearAlgebra
using JuMP
model = Model();
@variable(model, vf[1:3]);
@variable(model, rf[1:3]);
rv2oe(1.0, rf, vf)

The error on the inc > eps() conditional is:

ERROR: MethodError: no method matching isless(::Float64, ::NonlinearExpr)
The function `isless` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  isless(::Missing, ::Any)
   @ Base missing.jl:87
  isless(::Any, ::Missing)
   @ Base missing.jl:88
  isless(::AbstractFloat, ::ForwardDiff.Dual{Ty}) where Ty
   @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:145
  ...

Stacktrace:
 [1] <(x::Float64, y::NonlinearExpr)
   @ Base ./operators.jl:353
 [2] >(x::NonlinearExpr, y::Float64)
   @ Base ./operators.jl:379
 [3] rv2oe(μ::Float64, r::Vector{VariableRef}, v::Vector{VariableRef})
   @ Main ~/julia/jump-pseudospectral-examples.jl/rv2oe.jl:37
 [4] top-level scope
   @ REPL[7]:1

The error on the mod2pi() is:

ERROR: MethodError: no method matching rem2pi(::NonlinearExpr, ::RoundingMode{:Down})
The function `rem2pi` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  rem2pi(::Float32, ::RoundingMode)
   @ Base math.jl:1414
  rem2pi(::Float64, ::RoundingMode{:Down})
   @ Base math.jl:1349
  rem2pi(::Int32, ::RoundingMode)
   @ Base math.jl:1416
  ...

Stacktrace:
 [1] mod2pi(x::NonlinearExpr)
   @ Base.Math ./math.jl:1446
 [2] rv2oe(μ::Float64, r::Vector{VariableRef}, v::Vector{VariableRef})
   @ Main ~/julia/jump-pseudospectral-examples.jl/rv2oe.jl:42
 [3] top-level scope
   @ REPL[9]:1

Is there a way to fix this, while preserving the reusability of the rv2oe() function for other applications outside JuMP? I can obviously fork the function just for this application, but ideally I’d like to keep everything resuable, with the treatment of the edge conditions maintained, and not have to maintain a separate version just for this problem.

And I can promise for this application that the NLP solver stays away from these edge conditions and doesn’t need to e.g. calculate a derivative at the points where the functions become horrible (which is why I could fork the function if I really need to).

2

You’re hitting this limitation: Nonlinear Modeling · JuMP

So you cannot use the function directly. You’ll need to modify it.

The > you could write as:

lan = op_ifelse(
        op_strictly_greater_than(
            inc,
            eps(),
        ),
        atan(p, q),
        0.0,
    )          # Longitude of ascending node
    argp = op_ifelse(
        op_strictly_greater_than(
            ecc,
            eps(),
        ),
        atan(h, xk) - lan,
        0.0,
      )  # Argument of periapsis

But I would strongly consider changing variable bounds etc so that these cases cannot happen.

The mod2pi is muuuch harder, since it is inherently non-convex and discontinuous.

The usual approach would be something like this:

function my_mod2pi(model, x)
    div = @variable(model, integer = true)
    rem = @variable(model, lower_bound = 0, upper_bound = 2pi)
    @constraint(model, 2pi * div + rem == x)
    return rem
end

but now you’ve added integer variables to the problem and it becomes a MINLP.

Someone else might have suggestions for what to do here.

3

Yep, that’s what I thought.

Yeah, I was guessing it wouldn’t be possible, but I had to try to ask.

It isn’t actually that big of a deal, I’ll just have to fork it and work around it.

2 Likes