Source-linked AI summary
JuMP: A Modeling Language for Mathematical Optimization
Iain Dunning, Joey Huchette, Miles Lubin
TL;DR
JuMP addresses the need for a lightweight, extensible modeling language that fits modern scientific workflows without sacrificing performance. It embeds algebraic modeling in Julia, adds solver interaction and automatic differentiation, and supports many optimization classes. The paper reports competitive derivative-evaluation performance and faster model generation on larger instances, while identifying limits for forward-mode differentiation and parameter changes.
Problem
Existing commercial AMLs were relatively standalone, while open-source AMLs embedded in high-level languages motivated concerns about workflow integration and performance.
Method
JuMP embeds algebraic modeling in Julia and uses Julia features to support optimization classes, solver callbacks, automatic differentiation, and extensible modeling workflows.
Results
JuMP becomes significantly faster than Pyomo and YALMIP as instances grow, while derivative evaluation remains within 2.2x of AMPL and is up to 3x faster than GAMS.
Takeaways & Limitations
JuMP combines broad optimization modeling, high-level solver interaction, and integration with interactive visualization and scientific-computing workflows.
Takeaways & Limitations
Forward-mode automatic differentiation grows linearly with input dimension, and complex model-data changes may require rebuilding a JuMP model from scratch.
Abstract
from arXiv · showhide
JuMP is an open-source modeling language that allows users to express a wide range of optimization problems (linear, mixed-integer, quadratic, conic-quadratic, semidefinite, and nonlinear) in a high-level, algebraic syntax. JuMP takes advantage of advanced features of the Julia programming language to offer unique functionality while achieving performance on par with commercial modeling tools for standard tasks. In this work we will provide benchmarks, present the novel aspects of the implementation, and discuss how JuMP can be extended to new problem classes and composed with state-of-the-art tools for visualization and interactivity.
1. Introduction.
JuMP is an open-source algebraic modeling language embedded in Julia, designed to combine natural optimization modeling with performance, extensibility, and integration into modern scientific workflows.
- JuMP is an algebraic modeling language embedded in Julia and presented as an open-source alternative to commercial systems.
- The project targets lightweight integration with simulations, interactive visualization, programmatic model construction, and solver interaction during execution.
- Existing open-source AMLs address workflow integration but can suffer performance limitations from high-level host languages such as MATLAB and Python.
- JuMP extends modeling beyond linear and mixed-integer optimization to quadratic, conic-quadratic, semidefinite, and derivative-based nonlinear problems.
- The paper presents JuMP’s technical and usability features, including its Julia implementation, extensions, and interactive visualization applications.
2. The role of a modeling language.
Modeling languages let users express optimization problems algebraically while generating solver input and handling communication details, avoiding low-level matrix construction. JuMP adds high-level access to solver interaction and emphasizes natural expression of models.
- AMLs translate closed-form algebraic expressions into solver input structures, removing the need to write matrix generators by hand.
- JuMP’s modeling syntax is evaluated for how naturally mathematical statements translate into code, illustrated through a minimum cost flow formulation across four AMLs.
- A modeling language loads user input into memory and generates solver-specific input according to the problem class.
- JuMP provides a high-level in-memory interface for branch-and-cut callbacks, enabling bidirectional communication with solvers during solution.
3. Syntactic macros: parsing without a parser.
JuMP uses Julia’s syntactic macros to capture algebraic expressions as syntax and generate efficient modeling code without a custom text-based parser. This approach avoids key performance drawbacks of operator overloading.
- Operator overloading can make repeated quadratic-expression additions costly because intermediate expressions may contain O(n^2) terms.
- Figure 2 compares how JuMP, AMPL, Pyomo, and GAMS express the same minimum cost flow model.
- The naive Python accumulation example can require O(d^4) operations and excessive memory allocations.
- JuMP macros construct expression objects from Julia’s parsed syntax rather than treating algebraic model input as plain text.
- Syntactic macros provide natural algebraic modeling syntax without a custom text-based parser or operator-overloading drawbacks.
4. Code generation for linear and conic-quadratic models.
JuMP translates algebraic optimization expressions into solver-ready models using Julia’s macro and code-generation capabilities. This approach reduces expression-building overhead while supporting large linear, quadratic, and conic-quadratic models.
- Scope and evaluation: JuMP targets linear and conic-quadratic optimization, whose structural properties support efficient processing of large-scale models.The paper evaluates model-generation time before the solver begins and uses linear-quadratic control and facility-location instances as stress tests.
- Code generation: JuMP uses macros to transform algebraic input expressions into compiled code that generates solver input structures.Macros receive code, construct an expression data structure, and substitute the generated expression before compilation.
- Code generation: The generated code pre-allocates coefficient and index storage and appends terms directly to linear and quadratic expressions.This avoids creating a new partial expression for every addition during nested summation.
- Code generation: O(d^2) operations replace the O(d^4) naive operator-overloading approach for the illustrated quadratic expression.JuMP’s generated code is similar to a hand-written matrix generator.
JuMP
JuMP supports repeated solver interaction and model modification in memory, enabling iterative optimization workflows. Its design trades automatic parametric-data propagation for this solver-access model and incurs compilation startup cost.
- Evaluation: JuMP’s computational experiments use linear-quadratic control and facility-location problems implemented across seven modeling languages.The authors describe these models as stress tests rather than representative of all conic-quadratic problems.
- Trade-offs: A few seconds of compilation startup cost affects even the smallest JuMP instances, but repeated solves amortize this cost within one session.The compilation cost is paid only the first time an instance is solved in a sequence.
- Repeated solves: JuMP provides efficient in-memory solver access so related optimization problems can preserve solver state across solves when possible.This supports branch-and-bound, Benders decomposition, and cutting-plane algorithms without regenerating every model from scratch.
- Repeated solves: A cutting-plane implementation in GAMS was 5.8x slower overall than an implementation of the same algorithm in C++.The comparison illustrates overhead associated with traditional AMLs that lack direct in-memory solver access.
- Trade-offs: JuMP does not support automatic propagation of parametric values through an existing LP model as parameters change.For more complex data or structural changes, the stated idiom is to construct a new model, potentially inside a parameterized function.
5. Computing derivatives for nonlinear models.
JuMP computes derivatives for nonlinear models by representing algebraic expressions as graphs and applying automatic differentiation, including to user-defined functions. Its derivative implementation supports Hessian computation and performs competitively in benchmark comparisons with established AMLs.
- 5.1. Expression graphs and reverse-mode AD: JuMP represents nonlinear expressions as directed acyclic graphs encoding operation order and dependencies, enabling expression evaluation and derivative computation.Macros generate these graph objects from JuMP’s algebraic syntax.
- 5.1. Expression graphs and reverse-mode AD: Because JuMP restricts input to a specific syntax, it can guarantee efficient differentiation for valid expressions, subject to expression graphs fitting in memory.This contrasts with general-purpose AD tools that accept arbitrary code and require more specialized preparation.
- 5.1. Expression graphs and reverse-mode AD: JuMP uses reverse-mode automatic differentiation for nonlinear expressions and can compute Hessian-vector products before recovering sparse Hessian entries with fewer evaluations than dense methods.The illustrated Hessian is symmetric, and exploiting sparsity reduces the number of Hessian-vector products required.
- 5.2. User-defined functions: JuMP lets users register functions with hand-coded or automatically differentiated derivatives, extending nonlinear modeling to functions outside its built-in library.The example defines a squareroot function in generic Julia code, registers it, and enables automatic differentiation with autodiff=true.
- 5.2. User-defined functions: Forward-mode automatic differentiation is used for user-defined functions, but its cost grows linearly with input dimension, so callbacks may be preferable for high-dimensional functions.Users must write function code generic with respect to the numerical input type.
- 5.4. Benchmarks: Excluding the smallest instances, JuMP’s derivative evaluation stays within 2.2x of AMPL, is up to 3x faster than GAMS, and is at worst 25% slower.For model generation, JuMP has a large compiler-dominated startup cost but becomes significantly faster than Pyomo and YALMIP as instances grow.
JuMP
Tables 2 and 3 benchmark JuMP’s model-generation and derivative-evaluation performance against existing modeling languages for nonlinear optimization.
- Table 2 measures time to generate each model and pass it to the solver across derivative-based nonlinear optimization languages.A dash indicates that a model class is not implemented.
- Table 3 measures derivative-evaluation time, including gradients, Jacobians, and Hessians, during three Ipopt iterations.Pyomo uses AMPL’s solver library, while YALMIP does not provide second-order derivatives.
6. Extensions.
JuMP is designed to support extensible solver interfaces and higher-level modeling extensions, including parallel stochastic programming, robust optimization, and chance constraints.
- Extensions: JuMP supports extending its solver set and syntax to new problem classes, beyond the more typical AML focus on existing classes.Higher-level syntax can represent model structures more naturally than standard-form formulations.
- 6.1. Extension for parallel multistage stochastic programming: StructJuMP targets block-structured optimization arising in stochastic programming, contingency analysis, and multicommodity flow.The structure exposes matrices and vectors required by specialized methods such as Benders decomposition.
- 6.1. Extension for parallel multistage stochastic programming: StructJuMP scaled to 2048 cores, while model-generation overhead remained a small fraction of total solution time.Its implementation comprised fewer than 500 lines of code, including C++ solver and MPI interfaces.
- 6.2. Extension for robust optimization: JuMPeR models robust optimization directly through an Uncertain primitive and supports either reformulation or cutting-plane solution methods.Users can switch uncertainty sets and solution techniques with minimal code changes.
- 6.3. Extension for chance constraints: JuMPChance provides algebraic syntax for chance constraints, including Gaussian random variables and a probability parameter.In optimal power flow under wind-generation uncertainty, a more conservative distributionally robust model was associated with realized cost savings under inaccurate Gaussian assumptions.
7. Interactivity and visualization.
JuMP integrates optimization with interactive scientific workflows so users can manipulate inputs, re-solve models, and visualize resulting solution behavior in notebooks.
- 7. Interactivity and visualization: Interactive optimization and solution visualization are presented as important complements to communicating models to solvers.The paper contrasts this goal with the limited interactivity of standalone commercial systems and the commercial, Windows-only requirements of some alternatives.
- 7. Interactivity and visualization: Jupyter notebooks combine code, rich text, equations, visualizations, and interactive widgets in one shareable document.IJulia provides this notebook environment for Julia and JuMP.
- 7.1. Example: Portfolio Optimization: A Markowitz portfolio problem minimizes variance subject to a minimum mean-return requirement.The model is a quadratic optimization problem with linear constraints.
- 7.1. Example: Portfolio Optimization: The portfolio notebook combines formulation text, JuMP code, an interactive widget, and a visualization.Moving the rmin slider re-solves the model and updates the historical return-distribution plot.
- 7.2. Example: Goddard Rocket: The Goddard Rocket notebook varies maximum thrust and an altitude-drag coefficient through sliders, re-solving and plotting state and control trajectories.The implementation uses JuMP in IJulia and Interact.jl for exploration.
8. Appendix: Benchmark models for Section 4.1.
The appendix defines benchmark models used to evaluate JuMP across quadratic, constrained, mixed-integer, and conic-quadratic formulations.
- 8.1. lqcp: The linear-quadratic control benchmark scales by increasing the two-dimensional discretization parameters m and n.Benchmarking fixes m = n and varies n over {500, 1000, 1500, 2000}.
- 8.1. lqcp: The benchmark includes quadratic objectives, linear constraints, discretized index sets, control bounds, and state bounds.The displayed constraints use index sets I and J together with reduced sets I′ and J′.
- 8.2. fac: The facility-location benchmark minimizes the maximum distance from each customer to its nearest facility.It is formulated naturally as a mixed-integer second-order cone problem with deterministically generated data.
- 8.2. fac: Binary variables zc,f indicate whether facility f is closer to customer c than any other facility.The conic-quadratic constraint may be translated into an equivalent quadratic form depending on the modeling system.
9. Benchmark models for Section 5.4.
The benchmark models cover nonlinear beam control, nonlinear AC power flow, and facility location, with instances scaled to test model-generation performance. The AC power-flow benchmark is enlarged by duplicating its network, producing proportional increases in problem dimensions.
- Benchmark models: The facility-location example places three facilities among a four-by-four customer grid and minimizes the maximum customer-to-closest-facility distance.Customers are represented by rectangles, facilities by circles, and dotted circles show the objective distance.
- Benchmark models: clnlbeam is a nonlinear beam control model scaled with n ∈ {5000, 50000, 500000}.Each instance has 3n variables, 2n constraints, and diagonal Hessians.
- Benchmark models: The AC power-flow model minimizes active power losses while balancing active and reactive loads under nonlinear Kirchhoff constraints.Its parameters include nodal active power demand, voltage magnitude, phase angle, and complex-valued admittance.
- Benchmark models: The base AC network has 662 nodes, 1017 edges, 1489 decision variables, and 1324 constraints.Its Lagrangian Hessian has 8121 nonzero elements.
- Benchmark models: 10-fold and 100-fold network duplication produces proportional increases in the AC power-flow problem dimensions.The enlarged instances are generated from the base network.