# The JuMP ecosystem for mathematical optimization

## 2019 INFORMS-ALIO meeting


## Minimum # of Passports to Visit all Countries?
[](https://www.passportindex.org)

199 passports = $10^{33}$ times the age of the universe to enumerate at $10^{17}$ flops!


```julia
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
```

    [32m[1m  Updating[22m[39m registry at `~/.julia/registries/General`
    [32m[1m  Updating[22m[39m git-repo `https://github.com/JuliaRegistries/General.git`
    [?25l[2K[?25h

Download data from https://github.com/ilyankou/passport-index-dataset


```julia
;git clone https://github.com/ilyankou/passport-index-dataset.git
```

    Cloning into 'passport-index-dataset'...



```julia
using DelimitedFiles
data = readdlm(joinpath("passport-index-dataset","passport-index-matrix.csv"),',')
cntr = data[2:end,1]
vf = (x ->  x == -1 || x == 3 ? 1 : 0).(data[2:end,2:end]);
```

## (Constrained) Mathematical Optimization and JuMP

$$
\begin{alignedat}{3} 
&\min_{x,y} &\quad \sum_{\mathit{cntr} \;\in\; \mathit{World}} \mathit{pass}_{\,\mathit{cntr}} \\
&\text{s.t.}&\quad  \sum_{\mathit{cntr} \;\in\; \mathit{World}}\mathit{vf}(\mathit{cntr},\mathit{dst}) \cdot \mathit{pass}_{\,\mathit{cntr}} &\geq 1\quad  &\quad& \forall \; \mathit{dst} \;\in \; \mathit{World}\\
 &&\mathit{pass}_{\,\mathit{cntr}}  &\in \{0,1\}&\quad& \forall \; \mathit{cntr}\in \; \mathit{World}.
\end{alignedat}
$$


```julia
using JuMP, GLPK
model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, pass[1:length(cntr)], Bin)
@constraint(model, [j=1:length(cntr)], sum( vf[i,j]*pass[i] for i in 1:length(cntr)) >= 1)
@objective(model, Min, sum(pass))
JuMP.optimize!(model)
print(JuMP.objective_value(model)," passports: ",join(cntr[findall(JuMP.value.(pass) .== 1)],", "))
```

    23.0 passports: Afghanistan, Angola, Australia, Austria, Comoros, Congo, Eritrea, Gambia, Georgia, Hong Kong, India, Iraq, Kenya, Madagascar, Maldives, North Korea, Papua New Guinea, Singapore, Somalia, Sri Lanka, Tunisia, United Arab Emirates, United States

# A step by step example
Let's see how we translate a simple, 2 variable LP to JuMP code.

$$
\begin{align*}
&\max_{x,y}& \quad x + 2y \\
&\text{s.t.}&\quad x + y &\leq 1 \\
&&0\leq x, y &\leq 1
\end{align*}
$$



Load JuMP, MathOptInterface (MOI), and GLPK (GNU LP/MIP solver):


```julia
using JuMP  
using MathOptInterface # Replaces MathProgBase
# shortcuts
const MOI = MathOptInterface

using GLPK # Loading the GLPK module for using its solver
```

Construct a model object (a container for variables, constraints, solver options, etc.):


```julia
model = Model(with_optimizer(GLPK.Optimizer, msg_lev = 4));  
# Old syntax: model = Model(solver=GLPKSolverLP(msg_lev = 4)))
```

Define variables $0\leq x, y \leq 1$:


```julia
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
```




$$ y $$



Add constraint $x + y \leq 1$:


```julia
@constraint(model, x + y <= 1)
```




$ x + y \leq 1.0 $



Add objective $\max x + 2y$:


```julia
@objective(model, Max, x + 2y)
```




$$ x + 2 y $$



To solve the optimization problem, call the `optimize` function.


```julia
JuMP.optimize!(model) # Old syntax: status = JuMP.solve(model)
```

    GLPK Simplex Optimizer, v4.64
    1 row, 2 columns, 2 non-zeros
    *     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
    *     2: obj =   2.000000000e+00 inf =   0.000e+00 (0)
    OPTIMAL LP SOLUTION FOUND


We can then check the status of the optimization call.


```julia
@show JuMP.has_values(model)
@show JuMP.termination_status(model) == MOI.OPTIMAL
@show JuMP.primal_status(model) == MOI.FEASIBLE_POINT
@show JuMP.dual_status(model) == MOI.FEASIBLE_POINT
```

    JuMP.has_values(model) = true
    JuMP.termination_status(model) == MOI.OPTIMAL = true
    JuMP.primal_status(model) == MOI.FEASIBLE_POINT = true
    JuMP.dual_status(model) == MOI.FEASIBLE_POINT = true





    true



# New Solver Status

Much more details than old `:Optimal, :Unbounded, :Infeasible, :UserLimit, :Error, :NotSolved`

```julia 
@show JuMP.termination_status(model) == MOI.OPTIMAL
```


```julia
display(typeof(MOI.OPTIMAL))
```


    Enum MathOptInterface.TerminationStatusCode:
    OPTIMIZE_NOT_CALLED = 0
    OPTIMAL = 1
    INFEASIBLE = 2
    DUAL_INFEASIBLE = 3
    LOCALLY_SOLVED = 4
    LOCALLY_INFEASIBLE = 5
    INFEASIBLE_OR_UNBOUNDED = 6
    ALMOST_OPTIMAL = 7
    ALMOST_INFEASIBLE = 8
    ALMOST_DUAL_INFEASIBLE = 9
    ALMOST_LOCALLY_SOLVED = 10
    ITERATION_LIMIT = 11
    TIME_LIMIT = 12
    NODE_LIMIT = 13
    SOLUTION_LIMIT = 14
    MEMORY_LIMIT = 15
    OBJECTIVE_LIMIT = 16
    NORM_LIMIT = 17
    OTHER_LIMIT = 18
    SLOW_PROGRESS = 19
    NUMERICAL_ERROR = 20
    INVALID_MODEL = 21
    INVALID_OPTION = 22
    INTERRUPTED = 23
    OTHER_ERROR = 24


```julia
@show JuMP.primal_status(model) == MOI.FEASIBLE_POINT
```


```julia
display(typeof(MOI.FEASIBLE_POINT))
```


    Enum MathOptInterface.ResultStatusCode:
    NO_SOLUTION = 0
    FEASIBLE_POINT = 1
    NEARLY_FEASIBLE_POINT = 2
    INFEASIBLE_POINT = 3
    INFEASIBILITY_CERTIFICATE = 4
    NEARLY_INFEASIBILITY_CERTIFICATE = 5
    UNKNOWN_RESULT_STATUS = 6
    OTHER_RESULT_STATUS = 7


We can also inspect the solution values and optimal cost:


```julia
@show JuMP.value(x)              # Old syntax: getvalue(x)
@show JuMP.value(y)              # Old syntax: getvalue(y)
@show JuMP.objective_value(model)       # Old syntax: getobjectivevalue(model)
```

    JuMP.value(x) = 0.0
    JuMP.value(y) = 1.0
    JuMP.objective_value(model) = 2.0





    2.0



I can also "name" constraints for later reference.


```julia
model = Model(with_optimizer(GLPK.Optimizer, msg_lev = 0))
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@constraint(model, inequality, x + y <= 1)     # <=============== constraint can be referenced later as "inequality"
@objective(model, Max, x + 2y)
JuMP.optimize!(model)
@show JuMP.termination_status(model) == MOI.OPTIMAL
```

    JuMP.termination_status(model) == MOI.OPTIMAL = true





    true



constraint references  can by used to delete them


```julia
JuMP.delete(model, inequality)
JuMP.optimize!(model)
@show JuMP.termination_status(model) == MOI.OPTIMAL
@show JuMP.objective_value(model)
```

    JuMP.termination_status(model) == MOI.OPTIMAL = true
    JuMP.objective_value(model) = 3.0





    3.0



Constraint references can be used to modify problem (see [MOI](MOI.ipynb#Model-modifications)) and to get duals (see [Topics notebook](Topics.ipynb#Duality)).

## Collections of variables/constraints and summation in JuMP

You can also create collections of variables like $x_i \geq 0 \quad \forall \; i\in\{1,\ldots,10\}$


```julia
model = Model()
@variable(model, x[1:10] >= 0)
```




    10-element Array{VariableRef,1}:
     x[1] 
     x[2] 
     x[3] 
     x[4] 
     x[5] 
     x[6] 
     x[7] 
     x[8] 
     x[9] 
     x[10]



Also multidimensional indexing, separated by commas:


```julia
@variable(model, y[1:10,["red","blue"]] <= 1)
```




    2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
        Dimension 1, 1:10
        Dimension 2, ["red", "blue"]
    And data, a 10×2 Array{VariableRef,2}:
     y[1,red]   y[1,blue] 
     y[2,red]   y[2,blue] 
     y[3,red]   y[3,blue] 
     y[4,red]   y[4,blue] 
     y[5,red]   y[5,blue] 
     y[6,red]   y[6,blue] 
     y[7,red]   y[7,blue] 
     y[8,red]   y[8,blue] 
     y[9,red]   y[9,blue] 
     y[10,red]  y[10,blue]



and more complicated expressions like $\quad
i \leq z_{ij} \leq u_j \;\;\; \forall i \in \{1,...,10\}, j \in \{i+1, ..., 10\}
$:


```julia
u = rand(10)
@variable(model, i <= z[i=1:10,j=(i+1):10] <= u[j])
```




    JuMP.Containers.SparseAxisArray{VariableRef,2,Tuple{Any,Any}} with 45 entries:
      [8, 10]  =  z[8,10]
      [3, 6 ]  =  z[3,6]
      [6, 9 ]  =  z[6,9]
      [8, 9 ]  =  z[8,9]
      [1, 10]  =  z[1,10]
      [4, 5 ]  =  z[4,5]
      [2, 4 ]  =  z[2,4]
      [4, 9 ]  =  z[4,9]
      [1, 2 ]  =  z[1,2]
      [3, 4 ]  =  z[3,4]
               ⋮
      [4, 6 ]  =  z[4,6]
      [2, 10]  =  z[2,10]
      [1, 9 ]  =  z[1,9]
      [6, 10]  =  z[6,10]
      [1, 8 ]  =  z[1,8]
      [7, 10]  =  z[7,10]
      [6, 7 ]  =  z[6,7]
      [3, 10]  =  z[3,10]
      [3, 8 ]  =  z[3,8]
      [1, 5 ]  =  z[1,5]
      [3, 5 ]  =  z[3,5]



To specify conditions on the indexing, you can add conditionals inside the ``[...]`` block, separated by a semicolon:


```julia
@variable(model, w[i=1:10, c=["red","blue"]; iseven(i) || c == "red"] >= 0)
```




    JuMP.Containers.SparseAxisArray{VariableRef,2,Tuple{Any,Any}} with 15 entries:
      [10, blue]  =  w[10,blue]
      [7, red  ]  =  w[7,red]
      [8, red  ]  =  w[8,red]
      [10, red ]  =  w[10,red]
      [6, blue ]  =  w[6,blue]
      [1, red  ]  =  w[1,red]
      [5, red  ]  =  w[5,red]
      [2, blue ]  =  w[2,blue]
      [6, red  ]  =  w[6,red]
      [3, red  ]  =  w[3,red]
      [8, blue ]  =  w[8,blue]
      [4, red  ]  =  w[4,red]
      [4, blue ]  =  w[4,blue]
      [2, red  ]  =  w[2,red]
      [9, red  ]  =  w[9,red]



Also easy to create constrainta like $ \sum _{i = 1}^{10} x_i \leq 1$:


```julia
@constraint(model, sum(x[i] for i in 1:10) <= 1)
```




$ x_{1} + x_{2} + x_{3} + x_{4} + x_{5} + x_{6} + x_{7} + x_{8} + x_{9} + x_{10} \leq 1.0 $



Or more complicated ones like 
$
\sum_{\substack{i\in\{1,...,10\}\\
                c\in\{"red","blue"\}}}
       coef(c) \cdot y_{ic} = 1
$


```julia
coef = Dict("red" => 2, "blue" => 3)
@constraint(model, sum(coef[c]*y[i,c] for i in 1:10, c in ["red","blue"]) == 1)
```




$ 2 y_{1,red} + 3 y_{1,blue} + 2 y_{2,red} + 3 y_{2,blue} + 2 y_{3,red} + 3 y_{3,blue} + 2 y_{4,red} + 3 y_{4,blue} + 2 y_{5,red} + 3 y_{5,blue} + 2 y_{6,red} + 3 y_{6,blue} + 2 y_{7,red} + 3 y_{7,blue} + 2 y_{8,red} + 3 y_{8,blue} + 2 y_{9,red} + 3 y_{9,blue} + 2 y_{10,red} + 3 y_{10,blue} = 1.0 $



or $
\sum_{i = 1}^{10} \sum_{j = i+1}^{10} 
       i \cdot j \cdot z_{ij} \leq
\sum_{\substack{i\in\{1,...,10\},
                c\in\{"red","blue"\} \\
                \text{s.t. } iseven(i) \text{ or } c = "red"}}
       i^2 \cdot w_{ic} 
$:


```julia
@constraint(model, sum(i*j*z[i,j] for i in 1:10, j in (i+1):10) <=
               sum(i^2*w[i,c] for i in 1:10, c in ["red","blue"] if iseven(i) || c == "red"))
```




$ 2 z_{1,2} + 3 z_{1,3} + 4 z_{1,4} + 5 z_{1,5} + 6 z_{1,6} + 7 z_{1,7} + 8 z_{1,8} + 9 z_{1,9} + 10 z_{1,10} + 6 z_{2,3} + 8 z_{2,4} + 10 z_{2,5} + 12 z_{2,6} + 14 z_{2,7} + 16 z_{2,8} + 18 z_{2,9} + 20 z_{2,10} + 12 z_{3,4} + 15 z_{3,5} + 18 z_{3,6} + 21 z_{3,7} + 24 z_{3,8} + 27 z_{3,9} + 30 z_{3,10} + 20 z_{4,5} + 24 z_{4,6} + 28 z_{4,7} + 32 z_{4,8} + 36 z_{4,9} + 40 z_{4,10} + 30 z_{5,6} + 35 z_{5,7} + 40 z_{5,8} + 45 z_{5,9} + 50 z_{5,10} + 42 z_{6,7} + 48 z_{6,8} + 54 z_{6,9} + 60 z_{6,10} + 56 z_{7,8} + 63 z_{7,9} + 70 z_{7,10} + 72 z_{8,9} + 80 z_{8,10} + 90 z_{9,10} - w_{1,red} - 4 w_{2,red} - 4 w_{2,blue} - 9 w_{3,red} - 16 w_{4,red} - 16 w_{4,blue} - 25 w_{5,red} - 36 w_{6,red} - 36 w_{6,blue} - 49 w_{7,red} - 64 w_{8,red} - 64 w_{8,blue} - 81 w_{9,red} - 100 w_{10,red} - 100 w_{10,blue} \leq 0.0 $



Can also do collections of constraints (named or unamed):
    $$ 
\begin{align}
x_i &\leq 0.9 \quad \forall i \in \{1,2,3\} \quad\text{ (large bounds)}\\
x_i &\leq 0.5 \quad \forall i \in \{4,5,6\} 
\end{align}
$$


```julia
@constraint(model,largebounds[i=1:3], x[i] <= 0.9)
@constraint(model,[i=4:6], x[i] <= 0.5)
```




    1-dimensional DenseAxisArray{ConstraintRef{Model,C,Shape} where Shape<:AbstractShape where C,1,...} with index sets:
        Dimension 1, 4:6
    And data, a 3-element Array{ConstraintRef{Model,C,Shape} where Shape<:AbstractShape where C,1}:
     x[4] ≤ 0.5
     x[5] ≤ 0.5
     x[6] ≤ 0.5



# Classes of Constraints Beyond Linear Inequalities

Broadcasted and two sided linear inequalities:


```julia
A = [1.0 2.0; 3.0 4.0]
l = [4.0, 5.0]
u = [2.0, 3.0]
model = Model()
@variable(model, x[1:2])
@constraint(model, l .<= A*x .<= u)
```




    2-element Array{ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.Interval{Float64}},ScalarShape},1}:
     x[1] + 2 x[2] ∈ [4.0, 2.0]  
     3 x[1] + 4 x[2] ∈ [5.0, 3.0]



## Quadratic Inequalities:

Both convex:


```julia
model = Model()
@variable(model, x[1:2])
@constraint(model, x[1]^2 + x[2] <= 1)
```




$ x_{1}^2 + x_{2} \leq 1.0 $



and non-convex:


```julia
@constraint(model, x[1]*x[2] - 1.0 == 0.0)
```




$ x_{1}\times x_{2} = 1.0 $



## Conic constraints including...

Semidefinite constraints:


```julia
model = Model()                         # using CSDP; model = Model(with_optimizer(CSDP.CSDPOptimizer))
@variable(model, y[1:2,1:2], Symmetric)
@constraint(model, y in PSDCone())  
@variable(model, t)
@variable(model, w)
@SDconstraint(model,  [t 1; 1 -w] ⪰ [1 t; t -2])
```




$ \begin{bmatrix}
t - 1 & -t + 1\\
-t + 1 & -w + 2\\
\end{bmatrix} \in PSDCone() $



Second order cone constraints:
$$ 
\begin{equation}
\left\| Ax+u \right\|_2 \leq t
\end{equation}
$$


```julia
model = Model()
@variable(model, x[1:2])
@variable(model, t)
@constraint(model, [t;A*x+u] in SecondOrderCone())
```




$ [t, x_{1} + 2 x_{2} + 2, 3 x_{1} + 4 x_{2} + 3] \in MathOptInterface.SecondOrderCone(3) $



Rotated second order cone constraints:
$$ 
\begin{equation}
\left\| Ax+u \right\|_2 \leq t \cdot w,\quad w\geq 0
\end{equation}
$$


```julia
model = Model()
@variable(model, x[1:2])
@variable(model, t)
@variable(model, w)
@constraint(model, [t;w;A*x+u] in RotatedSecondOrderCone())
```




$ [t, w, x_{1} + 2 x_{2} + 2, 3 x_{1} + 4 x_{2} + 3] \in MathOptInterface.RotatedSecondOrderCone(4) $



# Also "derivative based" nonlinear constraints

Remains unchanged:


```julia
using Ipopt

model = Model(with_optimizer(Ipopt.Optimizer))
@variable(model, x)
@variable(model, y)

@NLobjective(model, Min, (1-x)^2 + 100(y-x^2)^2)
@NLconstraint(model, exp(x)+sin(x) <=0)

JuMP.optimize!(model)
```

    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    This is Ipopt version 3.12.10, running with linear solver mumps.
    NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
    
    Number of nonzeros in equality constraint Jacobian...:        0
    Number of nonzeros in inequality constraint Jacobian.:        1
    Number of nonzeros in Lagrangian Hessian.............:        4
    
    Total number of variables............................:        2
                         variables with only lower bounds:        0
                    variables with lower and upper bounds:        0
                         variables with only upper bounds:        0
    Total number of equality constraints.................:        0
    Total number of inequality constraints...............:        1
            inequality constraints with only lower bounds:        0
       inequality constraints with lower and upper bounds:        0
            inequality constraints with only upper bounds:        1
    
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
       0  1.0000000e+00 1.00e+00 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
       1  1.1271001e+01 6.01e-02 6.56e+01  -1.0 5.46e-01    -  1.00e+00 1.00e+00h  1
       2  3.1638662e+00 4.28e-02 1.79e+01  -1.0 3.17e-01    -  1.00e+00 7.13e-01f  1
       3  2.5851172e+00 0.00e+00 1.08e+00  -1.0 1.41e-01    -  1.00e+00 1.00e+00f  1
       4  2.6099717e+00 0.00e+00 1.25e-02  -1.0 1.21e-02    -  1.00e+00 1.00e+00h  1
       5  2.5289573e+00 0.00e+00 1.55e-01  -2.5 3.43e-02    -  1.00e+00 1.00e+00f  1
       6  2.5262625e+00 0.00e+00 3.78e-04  -2.5 1.51e-03    -  1.00e+00 1.00e+00h  1
       7  2.5235888e+00 0.00e+00 1.70e-04  -3.8 1.17e-03    -  1.00e+00 1.00e+00f  1
       8  2.5234381e+00 0.00e+00 5.52e-07  -5.7 6.62e-05    -  1.00e+00 1.00e+00h  1
       9  2.5234363e+00 0.00e+00 8.16e-11  -8.6 8.09e-07    -  1.00e+00 1.00e+00h  1
    
    Number of Iterations....: 9
    
                                       (scaled)                 (unscaled)
    Objective...............:   2.5234362583018233e+00    2.5234362583018233e+00
    Dual infeasibility......:   8.1634699000687760e-11    8.1634699000687760e-11
    Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
    Complementarity.........:   2.5074330278896033e-09    2.5074330278896033e-09
    Overall NLP error.......:   2.5074330278896033e-09    2.5074330278896033e-09
    
    
    Number of objective function evaluations             = 10
    Number of objective gradient evaluations             = 10
    Number of equality constraint evaluations            = 0
    Number of inequality constraint evaluations          = 10
    Number of equality constraint Jacobian evaluations   = 0
    Number of inequality constraint Jacobian evaluations = 10
    Number of Lagrangian Hessian evaluations             = 9
    Total CPU secs in IPOPT (w/o function evaluations)   =      2.360
    Total CPU secs in NLP function evaluations           =      1.220
    
    EXIT: Optimal Solution Found.



```julia

```

136219a468d2923b69283a60e3ab5022dff45b18
ipynb
Jupyter Notebook
JuMP.ipynb
juan-pablo-vielma/INFORMS-ALIO-2019-JuMP-Tutorial
30f8e37c35f833d967c39179b11c391d7ae43e6a
MIT
2019-06-03T02:33:02.000Z
2020-10-11T21:39:25.000Z
JuMP.ipynb
juan-pablo-vielma/INFORMS-ALIO-2019-JuMP-Tutorial
30f8e37c35f833d967c39179b11c391d7ae43e6a
MIT
JuMP.ipynb
juan-pablo-vielma/INFORMS-ALIO-2019-JuMP-Tutorial
30f8e37c35f833d967c39179b11c391d7ae43e6a
MIT
2019-08-06T06:05:01.000Z
2019-08-06T06:05:01.000Z
Qwen/Qwen-72B
1. YES
2. YES