# Copyright (c) 2019 Arpit Bhatia and contributors #src
# #src
# Permission is hereby granted, free of charge, to any person obtaining a copy #src
# of this software and associated documentation files (the "Software"), to deal #src
# in the Software without restriction, including without limitation the rights #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell #src
# copies of the Software, and to permit persons to whom the Software is #src
# furnished to do so, subject to the following conditions: #src
# #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software. #src
# #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE #src
# SOFTWARE. #src
# # Getting started with JuMP
# This tutorial is aimed at providing a quick introduction to writing and
# solving optimization models with JuMP.
# If you're new to Julia, start by reading [Getting started with Julia](@ref).
# ## What is JuMP?
# JuMP ("Julia for Mathematical Programming") is an open-source modeling
# language that is embedded in Julia. It allows users to formulate various
# classes of optimization problems (linear, mixed-integer, quadratic, conic
# quadratic, semidefinite, and nonlinear) with easy-to-read code. These problems
# can then be solved using state-of-the-art open-source and commercial solvers.
# JuMP also makes advanced optimization techniques easily accessible from a
# high-level language.
# ## What is a solver?
# A solver is a software package that incorporates algorithms for finding
# solutions to one or more classes of problem.
# For example, HiGHS is a solver for linear programming (LP) and mixed integer
# programming (MIP) problems. It incorporates algorithms such as the simplex
# method and the interior-point method.
# The [Supported-solvers](@ref) table lists the open-source and commercial
# solvers that JuMP currently supports.
# ## What is MathOptInterface?
# Each solver has its own concepts and data structures for representing
# optimization models and obtaining results.
# [MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl) (MOI) is
# an abstraction layer that JuMP uses to convert from the problem written in
# JuMP to the solver-specific data structures for each solver.
# MOI can be used directly, or through a higher-level modeling interface like
# JuMP.
# Because JuMP is built on top of MOI, you'll often see the `MathOptInterface.`
# prefix displayed when JuMP types are printed. However, you'll only need to
# understand and interact with MOI to accomplish advanced tasks such as creating
# [solver-independent callbacks](@ref callbacks_manual).
# ## Installation
# JuMP is a package for Julia. From Julia, JuMP is installed by using the
# built-in package manager.
# ```julia
# import Pkg
# Pkg.add("JuMP")
# ```
# You also need to include a Julia package which provides an appropriate solver.
# One such solver is `HiGHS.Optimizer`, which is provided by the
# [HiGHS.jl package](https://github.com/jump-dev/HiGHS.jl).
# ```julia
# import Pkg
# Pkg.add("HiGHS")
# ```
# See [Installation Guide](@ref) for a list of other solvers you can use.
# ## An example
# Let's solve the following linear programming problem using JuMP and HiGHS.
# We will first look at the complete code to solve the problem and then go
# through it step by step.
# Here's the problem:
# ```math
# \begin{aligned}
# & \min & 12x + 20y \\
# & \;\;\text{s.t.} & 6x + 8y \geq 100 \\
# & & 7x + 12y \geq 120 \\
# & & x \geq 0 \\
# & & y \in [0, 3] \\
# \end{aligned}
# ```
# And here's the code to solve this problem:
using JuMP
using HiGHS
model = Model(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@objective(model, Min, 12x + 20y)
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
print(model)
optimize!(model)
termination_status(model)
primal_status(model)
dual_status(model)
objective_value(model)
value(x)
value(y)
shadow_price(c1)
shadow_price(c2)
# ## Step-by-step
# Once JuMP is installed, to use JuMP in your programs write:
using JuMP
# We also need to include a Julia package which provides an appropriate solver.
# We want to use `HiGHS.Optimizer` here which is provided by the `HiGHS.jl`
# package:
using HiGHS
# JuMP builds problems incrementally in a `Model` object. Create a model by
# passing an optimizer to the [`Model`](@ref) function:
model = Model(HiGHS.Optimizer)
# Variables are modeled using [`@variable`](@ref):
@variable(model, x >= 0)
# !!! info
# The macro creates a new Julia object, `x`, in the current scope. We could
# have made this more explicit by writing:
# ```julia
# x = @variable(model, x >= 0)
# ```
# but the macro does this automatically for us to save writing `x` twice.
# Variables can have lower and upper bounds:
@variable(model, 0 <= y <= 30)
# The objective is set using [`@objective`](@ref):
@objective(model, Min, 12x + 20y)
# Constraints are modeled using [`@constraint`](@ref). Here, `c1` and `c2` are
# the names of our constraint:
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
# Call `print` to display the model:
print(model)
# To solve the optimization problem, call the [`optimize!`](@ref) function:
optimize!(model)
# !!! info
# The `!` after optimize is part of the name. It's nothing special.
# Julia has a convention that functions which mutate their arguments should
# end in `!`. A common example is `push!`.
# Now let's see what information we can query about the solution,
# starting with [`is_solved_and_feasible`](@ref):
is_solved_and_feasible(model)
# We can get more information about the solution by querying the three types of
# statuses.
# [`termination_status`](@ref) tells us why the solver stopped:
termination_status(model)
# In this case, the solver found an optimal solution.
# Check [`primal_status`](@ref) to see if the solver found a primal feasible
# point:
primal_status(model)
# and [`dual_status`](@ref) to see if the solver found a dual feasible point:
dual_status(model)
# Now we know that our solver found an optimal solution, and that it has a
# primal and a dual solution to query.
# Query the objective value using [`objective_value`](@ref):
objective_value(model)
# the primal solution using [`value`](@ref):
value(x)
value(y)
# and the dual solution using [`shadow_price`](@ref):
shadow_price(c1)
shadow_price(c2)
# !!! warning
# You should always check whether the solver found a solution before calling
# solution functions like [`value`](@ref) or [`objective_value`](@ref). A
# common workflow is:
# ```julia
# optimize!(model)
# if !is_solved_and_feasible(model)
# error("Solver did not find an optimal solution")
# end
# ```
# That's it for our simple model. In the rest of this tutorial, we expand on
# some of the basic JuMP operations.
# ## Model basics
# Create a model by passing an optimizer:
model = Model(HiGHS.Optimizer)
# Alternatively, call [`set_optimizer`](@ref) at any point before calling
# [`optimize!`](@ref):
model = Model()
set_optimizer(model, HiGHS.Optimizer)
# For some solvers, you can also use [`direct_model`](@ref), which offers a more
# efficient connection to the underlying solver:
model = direct_model(HiGHS.Optimizer())
# !!! warning
# Some solvers do not support [`direct_model`](@ref)!
# ### Solver Options
# Pass options to solvers with [`optimizer_with_attributes`](@ref):
model =
Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false))
# !!! note
# These options are solver-specific. To find out the various options
# available, see the GitHub README of the individual solver packages. The
# link to each solver's GitHub page is in the [Supported solvers](@ref)
# table.
# You can also pass options with [`set_attribute`](@ref):
model = Model(HiGHS.Optimizer)
set_attribute(model, "output_flag", false)
# ## Solution basics
# We saw above how to use [`termination_status`](@ref) and
# [`primal_status`](@ref) to understand the solution returned by the solver.
# However, only query solution attributes like [`value`](@ref) and
# [`objective_value`](@ref) if there is an available solution. Here's a
# recommended way to check:
function solve_infeasible()
model = Model(HiGHS.Optimizer)
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@constraint(model, x + y >= 3)
@objective(model, Max, x + 2y)
optimize!(model)
if !is_solved_and_feasible(model)
@warn("The model was not solved correctly.")
return
end
return value(x), value(y)
end
solve_infeasible()
# ## Variable basics
# Let's create a new empty model to explain some of the variable syntax:
model = Model()
# ### Variable bounds
# All of the variables we have created till now have had a bound. We can also
# create a free variable.
@variable(model, free_x)
# While creating a variable, instead of using the <= and >= syntax, we can also
# use the `lower_bound` and `upper_bound` keyword arguments.
@variable(model, keyword_x, lower_bound = 1, upper_bound = 2)
# We can query whether a variable has a bound using the `has_lower_bound` and
# `has_upper_bound` functions. The values of the bound can be obtained using the
# `lower_bound` and `upper_bound` functions.
has_upper_bound(keyword_x)
upper_bound(keyword_x)
# Note querying the value of a bound that does not exist will result in an error.
try #hide
lower_bound(free_x)
catch err #hide
showerror(stderr, err) #hide
end #hide
# ### [Containers](@id tutorial_variable_container)
# We have already seen how to add a single variable to a model using the
# [`@variable`](@ref) macro. Now let's look at ways to add multiple variables to
# a model.
# JuMP provides data structures for adding collections of variables to a model.
# These data structures are referred to as _containers_ and are of three types:
# `Arrays`, `DenseAxisArrays`, and `SparseAxisArrays`.
# #### Arrays
# JuMP arrays are created when you have integer indices that start at `1`:
@variable(model, a[1:2, 1:2])
# Index elements in `a` as follows:
a[1, 1]
#-
a[2, :]
# Create an n-dimensional variable $x \in {R}^n$ with bounds $l \le x \le u$
# ($l, u \in {R}^n$) as follows:
n = 10
l = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
u = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
@variable(model, l[i] <= x[i = 1:n] <= u[i])
# We can also create variable bounds that depend upon the indices:
@variable(model, y[i = 1:2, j = 1:2] >= 2i + j)
# #### DenseAxisArrays
# `DenseAxisArrays` are used when the indices are not one-based integer ranges.
# The syntax is similar except with an arbitrary vector as an index as opposed
# to a one-based range:
@variable(model, z[i = 2:3, j = 1:2:3] >= 0)
# Indices do not have to be integers. They can be any Julia type:
@variable(model, w[1:5, ["red", "blue"]] <= 1)
# Index elements in a `DenseAxisArray` as follows:
z[2, 1]
#-
w[2:3, ["red", "blue"]]
# See [Forcing the container type](@ref variable_forcing) for more details.
# #### SparseAxisArrays
# `SparseAxisArrays` are created when the indices do not form a Cartesian product.
# For example, this applies when indices have a dependence upon previous indices
# (called triangular indexing):
@variable(model, u[i = 1:2, j = i:3])
# We can also conditionally create variables by adding a comparison check that
# depends upon the named indices and is separated from the indices by a
# semi-colon `;`:
@variable(model, v[i = 1:9; mod(i, 3) == 0])
# Index elements in a `DenseAxisArray` as follows:
u[1, 2]
#-
v[[3, 6]]
# ### Integrality
# JuMP can create binary and integer variables. Binary variables are constrained
# to the set ``\{0, 1\}``, and integer variables are constrained to the set
# ``\mathbb{Z}``.
# #### Integer variables
# Create an integer variable by passing `Int`:
@variable(model, integer_x, Int)
# or setting the `integer` keyword to `true`:
@variable(model, integer_z, integer = true)
# #### Binary variables
# Create a binary variable by passing `Bin`:
@variable(model, binary_x, Bin)
# or setting the `binary` keyword to `true`:
@variable(model, binary_z, binary = true)
# ## Constraint basics
# We'll need a new model to explain some of the constraint basics:
model = Model()
@variable(model, x)
@variable(model, y)
@variable(model, z[1:10]);
# ### [Containers](@id tutorial_constraint_container)
# Just as we had containers for variables, JuMP also provides `Arrays`,
# `DenseAxisArrays`, and `SparseAxisArrays` for storing collections of
# constraints. Examples for each container type are given below.
# #### Arrays
# Create an `Array` of constraints:
@constraint(model, [i = 1:3], i * x <= i + 1)
# #### DenseAxisArrays
# Create an `DenseAxisArray` of constraints:
@constraint(model, [i = 1:2, j = 2:3], i * x <= j + 1)
# #### SparseAxisArrays
# Create an `SparseAxisArray` of constraints:
@constraint(model, [i = 1:2, j = 1:2; i != j], i * x <= j + 1)
# ### Constraints in a loop
# We can add constraints using regular Julia loops:
for i in 1:3
@constraint(model, 6x + 4y >= 5i)
end
# or use for each loops inside the `@constraint` macro:
@constraint(model, [i in 1:3], 6x + 4y >= 5i)
# We can also create constraints such as $\sum _{i = 1}^{10} z_i \leq 1$:
@constraint(model, sum(z[i] for i in 1:10) <= 1)
# ## Objective functions
# Set an objective function with [`@objective`](@ref):
model = Model(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
@objective(model, Min, 2x + y)
# Create a maximization objective using `Max`:
@objective(model, Max, 2x + y)
# !!! tip
# Calling [`@objective`](@ref) multiple times will over-write the previous
# objective. This can be useful when you want to solve the same problem with
# different objectives.
# ## Vectorized syntax
# We can also add constraints and an objective to JuMP using vectorized linear
# algebra. We'll illustrate this by solving an LP in standard form that is,
# ```math
# \begin{aligned}
# & \min & c^T x \\
# & \;\;\text{s.t.} & A x = b \\
# & & x \ge 0
# \end{aligned}
# ```
vector_model = Model(HiGHS.Optimizer)
A = [1 1 9 5; 3 5 0 8; 2 0 6 13]
b = [7, 3, 5]
c = [1, 3, 5, 2]
@variable(vector_model, x[1:4] >= 0)
@constraint(vector_model, A * x .== b)
@objective(vector_model, Min, c' * x)
optimize!(vector_model)
@assert is_solved_and_feasible(vector_model)
objective_value(vector_model)