Performance problem with deeply nested nonlinear expressions
Les mainteneurs répondent en général sous 1 jour
Évaluation
- Difficulté
- 4/5
- Temps estimé
- 3-5 jours
- Accessibilité débutants
- 45/100
- Type d'issue
- Bug
- Clarté
- Plutôt claire
- Activité
- Calme
- Stack technique
- julia
- Domaine
- performance, tooling
Piste de recherche
Commencez par les reproductions Julia fournies, en particulier le modèle RK4 multistep avec N = 2 et la comparaison en un seul pas avec N = 50. Mesurez séparément par benchmark la construction du modèle et l’optimisation, et examinez comment les expressions imbriquées croissent au cours de l’expansion répétée de RK4. Le travail sera considéré comme terminé une fois qu’il aura été déterminé si le ralentissement se situe dans JuMP/MathOptInterface ou dans le schéma de modélisation, avec un benchmark minimal et une orientation de codage documentée.
Rédigé par le modèle d'indexation à partir du texte de l'issue.
Description
Hello, I am interested in switching from CasADi tu JuMP for optimal control (the system is quaternion kinematics, so not a lot of variables). My MATLAB CasADi model does multistep RK4 integration for multiple shooting (4 steps/interval, 10 intervals) and builds and solves the problem in < 2s. The equivalent JuMP code builds the model for many minutes, using up to 10GB of RAM, and is basically unusable. This is the JuMP code I wrote:
using ReferenceFrameRotations, DifferentialEquations, Plots
##
function qdot(X, u)
collect(dquat(Quaternion(X), u))
end
function RK4(f, X, u, dt)
k1 = dt*f(X , u)
k2 = dt*f(X+k1/2, u)
k3 = dt*f(X+k2/2, u)
k4 = dt*f(X+k3 , u)
X + (k1+2k2+2k3+k4)/6
end
model = Model(Ipopt.Optimizer)
DCMf = [
-1 0 0
0 0 1
0 1 0
]
# two quaternions represent the same orientation
qf = -convert(Quaternion, DCM(DCMf)) |> collect
N = 2
T = 1
dt = T / (N-1)
tabq = @variable(model, [1:4, 1:N], base_name="q")
tabu = @variable(model, [1:3, 1:N-1], base_name="ω_b")
cost = 0
for i = 1:N-1
F = [tabq[:, i]; cost]
for M = 1:4
F = RK4((X, u) -> [qdot(X[1:4], u); 1/2 * u' * u], F, tabu[:, i], dt)
end
@constraint(model, tabq[:, i+1] == F[1:4])
cost = F[end]
end
#traj constraints
@constraint(model, tabq[:, 1] == q0)
@constraint(model, tabq[:, N] == qf)
# @constraint(model, c[i=1:N], tabq[:, i]' * tabq[:, i] == 1)
@objective(model, Min, cost)
#start values
set_start_value.(tabq, q0)
I set N = 2 to prove it's slow even for a single step (similar problems happen with single-step high order RK schemes, such as RK8, which leads me to believe JuMP can't handle big expressions really well). For comparison, if I build a single-step model with 50 steps (approx. 4steps/interval x 10 intervals for the working CasADi model) like this:
model = Model(Ipopt.Optimizer)
DCMf = [
-1 0 0
0 0 1
0 1 0
]
# two quaternions represent the same orientation
qf = -convert(Quaternion, DCM(DCMf)) |> collect
N = 50
T = 1
dt = T / (N-1)
tabq = @variable(model, [1:4, 1:N], base_name="q")
tabu = @variable(model, [1:3, 1:N-1], base_name="ω_b")
cost = 0
for i = 1:N-1
#can't do multistep integration
F = RK4((X, u) -> [qdot(X[1:4], u); 1/2 * u' * u], [tabq[:, i]; cost], tabu[:, i], dt)
@constraint(model, tabq[:, i+1] == F[1:4])
cost = F[end]
end
#traj constraints
@constraint(model, tabq[:, 1] == q0)
@constraint(model, tabq[:, N] == qf)
# @constraint(model, c[i=1:N], tabq[:, i]' * tabq[:, i] == 1)
@objective(model, Min, cost)
#start values
set_start_value.(tabq, q0)
model
##
optimize!(model)
It then works reasonably fast. I don't understand why JuMP can't handle the first code; in my view, this is a big limitation. Is there a particular way of coding in JuMP I'm missing, or is this truly JuMP's fault?
- Langage dominant
- Julia
- Étoiles
- 506
- Forks
- 103
- Merge moyen
- 1 j 10 min
- PR mergées (30 j)
- 12
Préparer son environnement
Ce projet ne fournit ni conteneur de développement, ni Dockerfile, ni guide de contribution : l'installation est à votre charge. Commencez par son README, et consultez notre guide de la première contribution pour les étapes générales.
Par où commencer
- Lisez l'issue en entier, puis le guide de contribution du projet.
- Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
- Forkez le dépôt et travaillez sur une branche.
- Ouvrez une pull request qui référence le numéro de l'issue.
Autres issues de jump-dev/MathOptInterface.jl
-
[Bridges] MultirowChange in set_map.jl drops zero coefficientsPeut-être pris @odow l’a pris il y a 5 jours. Ouverte
jump-dev/MathOptInterface.jl#3106 · 2 commentaires · 1 personne assignée ·
Les mainteneurs répondent en général sous 1 jour
-
Precision loss in Hessian computationPeut-être pris @blegat l’a pris il y a 6 jours. Ouverte
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 15/100
jump-dev/MathOptInterface.jl#3104 · 11 commentaires ·
Les mainteneurs répondent en général sous 1 jour
-
Layer instantiatorOuverte
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 35/100
jump-dev/MathOptInterface.jl#2986 · 3 commentaires ·
Les mainteneurs répondent en général sous 1 jour
-
Differentiable Optimization APIOuverte
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 25/100
jump-dev/MathOptInterface.jl#2969 · 2 commentaires · 1 réaction ·
Les mainteneurs répondent en général sous 1 jour
-
Add Batched{S} setOuverteType: Set Request
Difficulté 5/5 Plus d'une semaine Accessibilité débutants 35/100
jump-dev/MathOptInterface.jl#2904 · 15 commentaires · 7 réactions ·
Les mainteneurs répondent en général sous 1 jour
Toutes les issues de jump-dev/MathOptInterface.jl
Issues similaires
-
Difficulté 2/5 1-3 heures Accessibilité débutants 76/100
Les mainteneurs répondent en général sous 1 jour
-
Difficulté 2/5 1-3 heures Accessibilité débutants 72/100
EnzymeAD/Reactant.jl#3549 ·
Les mainteneurs répondent en général sous 1 jour
-
bug
Difficulté 2/5 1-3 heures Accessibilité débutants 75/100
-
bug help wanted
Difficulté 2/5 1-3 heures Accessibilité débutants 68/100
JuliaGeometry/Meshes.jl#1469 ·
Les mainteneurs répondent en général sous 1 jour
-
distributed 🕸️
Difficulté 2/5 1-3 heures Accessibilité débutants 62/100
CliMA/Oceananigans.jl#6295 · 1 commentaire ·
Les mainteneurs répondent en général sous 1 jour