# JuMP operators for \> and mod2pi?
**URL:** https://discourse.julialang.org/t/jump-operators-for-and-mod2pi/122871
**Category:** Optimization (Mathematical)
**Tags:** jump
**Created:** [November 20, 2024, 8:41pm UTC](https://discourse.julialang.org/t/jump-operators-for-and-mod2pi/122871 "2024-11-20T20:41:44Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@lamont](https://discourse.julialang.org/u/lamont)
#### Post date: [November 20, 2024, 8:41pm UTC](https://discourse.julialang.org/t/jump-operators-for-and-mod2pi/122871/1 "2024-11-20T20:41:44Z")
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:
```julia
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:
```julia
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:
```julia
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:
```julia
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).
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [November 20, 2024, 9:00pm UTC](https://discourse.julialang.org/t/jump-operators-for-and-mod2pi/122871/2 "2024-11-20T21:00:07Z")
You’re hitting this limitation: [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Limitations)
So you cannot use the function directly. You’ll need to modify it.
The `>` you could write as:
```julia
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:
```Julia
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.
---
### Author:  [@lamont](https://discourse.julialang.org/u/lamont)
#### Post date: [November 21, 2024, 1:50am UTC](https://discourse.julialang.org/t/jump-operators-for-and-mod2pi/122871/3 "2024-11-21T01:50:56Z")
> [@odow](#):
>
> You’re hitting this limitation: [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Limitations)
Yep, that’s what I thought.
> [@](#):
>
> The `mod2pi` is muuuch harder, since it is inherently non-convex and discontinuous.
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.