# Introduction to [JuMP](https://github.com/JuliaOpt/JuMP.jl/blob/master/README.md)

February 27, 2020

*Uses JuMP version 0.21*

### What is JuMP?
JuMP stands for **Ju**lia for **M**athematical Programming. It is a "domain-specific modeling language for mathematical optimization embedded in Julia."

What does "domain-specific" mean? It means that some algorithms are useful in some knowledge domains (e.g. supply chain management) that may not be useful in other knowledge domains (e.g. estimating statistical model parameters).

What is great about JuMP is that it allows the user to change algorithms with little or no change to the modeling code. This means that someone without domain-specific knowledge can easily utilize JuMP to accomplish a task that would otherwise be impossible.

### What goes into JuMP?
To use JuMP, a user must specify the following components:

- An objective function (you gotta solve something!)
- Variables
- Constraints

JuMP does the rest! (I'll show you examples shortly of how you can do this.)

### What kinds of problems can JuMP handle?
JuMP provides the ability to solve the following types of objective functions:

- Linear
- Convex Quadratic
- Nonlinear (convex and nonconvex)

It also provides the ability to use constraints of the following types:

- Linear
- Convex Quadratic
- Second-order Conic
- Semidefinite
- Nonlinear (convex and nonconvex)

Finally, it provides support for any of the following types of variables:

- Continuous
- Integer-valued
- Semicontinuous
- Semi-integer

### What examples will we go through today?
Today I'll go through four examples:

1. Solve a very simple constrained system of equations (adapted from JuMP's Github repository [here](https://github.com/JuliaOpt/JuMP.jl/blob/master/examples/basic.jl))
2. Solve a Sudoku puzzle (adapted from JuMP's Github repository [here](https://github.com/JuliaOpt/JuMP.jl/blob/master/examples/sudoku.jl))
3. Solve a professional sports team's salary budget problem
4. Estimate the parameters of a linear regression model

---

## 1. A simple constrained system of equations
Suppose you want to solve the following (constrained) system of equations:

\begin{equation*}
\max\,\, 5x + 3y
\end{equation*}

subject to the following constraints:

\begin{align}
x + 5y & \leq 3 \\
x & \geq 0 \\
y & \geq 0 \\
x & \leq 2 \\
y & \leq 30 \\
\end{align}

In economics, this problem may represent a utility maximization problem, where utility is linear in $x$ and $y$, the individul has an income of 3, and the price of good $x$ is 1, with the price of good $y$ equal to 5. (And there are supply constraints on $x$ and $y$ governed by the market in general.)

If we wanted to, we could solve this using a Lagrangian (although in this case with linear utility, we know we will have a corner solution, so no need to bother with calculus). 

But we want the computer to do this, so let's see how it's done.

### Specifying the model components
As explained above, JuMP needs an objective function, variables, and constraints. Additionally, we need to tell JuMP which optimization algorithm to use.

#### Optimizer
The first step is to declare a model and attach an optimizer. We are going to use the GLPK optimizer:
```julia
model = Model(GLPK.Optimizer)
```

#### Variables
Next, we need to tell it what the variables are, and any constraints on those variables:
```julia
@variable(model, 0 <= x <= 2)
@variable(model, 0 <= y <= 30)
```

#### Objective function
Next, we tell it the objective function (and whether we want to maximize or minimize):
```julia
@objective(model, Max, 5x + 3y)
```

#### Constraints
Finally, we give it the constraints. Note that the single constraints on each variable were incorporated when we declared the variables themselves, but we have one additional constraint that is a function of both variables. We could also impose a constraint that the optimal values of $x$ and $y$ be integers (if these represented indivisible objects, for example).
```julia
@constraint(model, 1x + 5y <= 3)
```

### Optimizing the model
Once we have all of the components declared, we can optimize the model:
```julia
JuMP.optimize!(model)
```

We can then look at the output of the optimization as follows:
```julia
obj_value = JuMP.objective_value(model)
x_value   = JuMP.value(x)
y_value   = JuMP.value(y)
```
We'll put all of the code together below so it can be run on your machine.

Following Julia protocol, we will also wrap everything in a function and then call that function.

### All of the code together


```julia
using JuMP, GLPK


# wrap all of our code inside a function (for better performance)
function example_basic()
    
    # define model and optimizer
    model = Model(GLPK.Optimizer)
    
    # define variables
    @variable(model, 0 <= x <= 2)
    @variable(model, 0 <= y <= 30)

    # define objective function
    @objective(model, Max, 5x + 3y)
    
    # add additional constraints
    @constraint(model, 1x + 5y <= 3.0)

    # display the model
    print(model)
    
    # optimize the model
    JuMP.optimize!(model)

    # return and print objective function and optimal values of variables
    obj_value = JuMP.objective_value(model)
    x_value = JuMP.value(x)
    y_value = JuMP.value(y)
    println("Objective value: ", obj_value)
    println("x = ", x_value)
    println("y = ", y_value)
end


# call the function defined above
example_basic()
```

    Max 5 x + 3 y
    Subject to
     x + 5 y <= 3.0
     x >= 0.0
     y >= 0.0
     x <= 2.0
     y <= 30.0
    Objective value: 10.6
    x = 2.0
    y = 0.2


So our optimal values are:

\begin{align}
x^* & = 2 \\
y^* & = 0.2 \\
\text{Objective} & = 10.6
\end{align}

As I mentioned before we did any programming, this problem would yield a corner solution (where, for one of the goods, the optimal value hits one of the constraints). In this case, it was $x^* = 2$.

#### Other types of constraints
As I mentioned above, we can also add other constraints, for example, that the optimal values be integers (e.g. if $x$ and $y$ are indivisible). In this case, we would add `@constraint(model, x in MOI.Integer())` and similar for $y$.

We could also put a constraint on the objective function itself (e.g. that the objective value be no larger than 10). To do this we would add the following:
```julia
@expression(model, objval, 5x + 3y)
@constraint(model, objval <= 10)
```
---

## 2. Solving a Sudoku puzzle
Now that we are more comfortable with JuMP, we can solve a Sudoku puzzle. All we need to do is appropriately tell JuMP how to understand the puzzle board.

### Sudoku objective
In this case, the objective function is to fill in the puzzle board given a starting grid that has some numbers filled in. We won't have a formal objective function for this; we will just give JuMP a starting grid and tell it to satisfy all of the constraints, where one of the constraints is that the board gets filled. (A blank board has an infinite number of solutions, but a partially completed board should have just one unique solution.)

### Sudoku variables
We will communicate the state of the puzzle board with an array of variables. The variables will be put in a 3-dimensional array, where the first two dimensions tell the "latitude" and "longitude" of the cell on the puzzle board, and the third dimension keeps track of which of the numbers 1-9 will fill that cell.

Mathematically, we have
```julia
x[i, j, k]
```
which, if equal to 1, indicates that cell $(i,j)$ should contain the number $k$. The indices $(i,j,k)$ each must take on integer values from 1 to 9, since the puzzle board has 81 squares. So our `x` array is a 9 x 9 x 9 cube.

### Sudoku constraints
The constraints of Sudoku are as follows:

1. Each cell can only contain one number (duh, but we have to explain this to the computer!)
2. Each row contains each number exactly once
3. Each column contains each number exactly once
4. Each 3x3 subgrid contains each number exactly once

### Code


```julia
using JuMP, GLPK
function example_sudoku()
    
    # input the initial puzzle board (0s mean blanks)
    initial_grid = [
                    3 1 0 0 5 8 0 0 4;
                    0 0 9 3 2 0 0 0 0;
                    0 2 5 1 0 4 0 9 0;
                    0 0 0 0 0 0 3 8 9;
                    0 0 8 0 0 0 5 0 0;
                    5 4 6 0 0 0 0 0 0;
                    0 8 0 2 0 3 6 5 0;
                    0 0 0 0 7 1 4 0 0;
                    7 0 0 4 8 0 0 2 1
                    ]

    # use GLPK Optimizer
    model = Model(GLPK.Optimizer)
    
    # Set up the variables: each one can only take on binary values, so we add "Bin" to the end as a constraint
    @variable(model, x[1:9, 1:9, 1:9], Bin)

    # Add the constraints
    @constraints(model, begin
                     # Constraint 1 - Only one value appears in each cell
                     cell[i in 1:9, j in 1:9], sum(x[i, j, :]) == 1
                     # Constraint 2 - Each value appears in each row once only
                     row[i in 1:9, k in 1:9], sum(x[i, :, k]) == 1
                     # Constraint 3 - Each value appears in each column once only
                     col[j in 1:9, k in 1:9], sum(x[:, j, k]) == 1
                     # Constraint 4 - Each value appears in each 3x3 subgrid once only
                     subgrid[i=1:3:7, j=1:3:7, val=1:9], sum(x[i:i + 2, j:j + 2, val]) == 1
                 end)

    # Add additional constraints that reflect the starting point of the puzzle board
    # (i.e. don't attempt to update the numbers that were given as part of the puzzle)
    for row in 1:9, col in 1:9
        if initial_grid[row, col] != 0
            @constraint(model, x[row, col, initial_grid[row, col]] == 1)
        end
    end

    # Solve it
    JuMP.optimize!(model)

    term_status = JuMP.termination_status(model)
    primal_status = JuMP.primal_status(model)
    is_optimal = term_status == MOI.OPTIMAL

    # Check solution
    if is_optimal
        mip_solution = JuMP.value.(x)
        sol = zeros(Int, 9, 9)
        for row in 1:9, col in 1:9, val in 1:9
            if mip_solution[row, col, val] >= 0.9
                sol[row, col] = val
            end
        end
        return sol
    else
        error("The solver did not find an optimal solution.")
    end
end

function print_sudoku_solution(solution)
    println("Solution:")
    println("[-----------------------]")
    for row in 1:9
        print("[ ")
        for col in 1:9
            print(solution[row, col], " ")
            if col % 3 == 0 && col < 9
                print("| ")
            end
        end
        println("]")
        if row % 3 == 0
            println("[-----------------------]")
        end
    end
end

sol = example_sudoku()
print_sudoku_solution(sol)
```

    Solution:
    [-----------------------]
    [ 3 1 7 | 9 5 8 | 2 6 4 ]
    [ 4 6 9 | 3 2 7 | 8 1 5 ]
    [ 8 2 5 | 1 6 4 | 7 9 3 ]
    [-----------------------]
    [ 2 7 1 | 6 4 5 | 3 8 9 ]
    [ 9 3 8 | 7 1 2 | 5 4 6 ]
    [ 5 4 6 | 8 3 9 | 1 7 2 ]
    [-----------------------]
    [ 1 8 4 | 2 9 3 | 6 5 7 ]
    [ 6 9 2 | 5 7 1 | 4 3 8 ]
    [ 7 5 3 | 4 8 6 | 9 2 1 ]
    [-----------------------]


---
# 3. Solving a salary cap problem
NBA general managers want to build a championship team. What is their objective function? To win a championship. Well, that's difficult to exactly write down, but we can look at other more easily measurable outputs (like points scored, points allowed, etc.)
    
Suppose we want to find the team that will have the best statistics, but at the cheapest price.

### Data
I obtained the data using [this R script](https://github.com/tyleransom/DScourseS20/blob/master/WebData/getNBAplayerStats.R), which makes use of the `nbastatR` package. The data is in CSV format, which we will directly read into Julia using the `HTTP` and `JuliaDB` packages.

### Objective function
It's not clear what objective function we should use, but let's start with a simple one: maximize points scored.

### Constraints
We have at least two constraints: our team can only have 15 players on it, and the total team salary must be below the luxury tax threshold (\\$132.6m). We will talk about other important constraints after trying things out with this most basic constraint.

### Variables
In this case, the variables are the players that we pick. 

### Code


```julia
using HTTP, JuliaDB, JuMP, GLPK

# first function: read in the player data from the class GitHub repository
function read_in_data(url)
    newtable  = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    players   = newtable[1][2]
    salaries  = newtable[1][4]./1000000
    ppg       = newtable[1][11]
    return players,salaries,ppg
end

function SolveModel(players,salary,points)
    N = length(salary) 
    
    m = Model(GLPK.Optimizer)

    # define the variables: they are 0 if the player did not make the team, 1 if the player did make the team
    @variable(m, picked[1:N], Bin)

    # categories: 
    @objective(m, Max, sum( points[i] * picked[i] for i in 1:N)) 

    @constraints m begin
        # Constraint 1 - payroll <= 132.6m
        sum(salary[i] * picked[i] for i in 1:N) <= 132.6
        # Constraint 2 - must have exactly 15 players on roster
        sum(picked[i] for i in 1:N) == 15
    end

    # Solve it
    JuMP.optimize!(m);
    pck     = convert(BitArray,JuMP.value.(picked))
    lineup  = players[pck]
    points  = JuMP.objective_value(m)
    payroll = sum(salary[pck])
    return lineup,points,payroll
end

# call first function (to import data)
players,salaries,pts = read_in_data("https://raw.githubusercontent.com/tyleransom/DScourseS20/master/WebData/playerSalaryStats.csv")

# pass data into second function to get optimal lineup
lineup,total_points,payroll = SolveModel(players,salaries,pts)
println("team: ",lineup)
println("total points scored per game: ",total_points)
println("payroll: ",payroll)
```

    team: ["Trae Young", "Jayson Tatum", "Jaylen Brown", "Zach LaVine", "Collin Sexton", "Domantas Sabonis", "Kendrick Nunn", "Giannis Antetokounmpo", "Brandon Ingram", "Shai Gilgeous-Alexander", "Buddy Hield", "Pascal Siakam", "Donovan Mitchell", "Luka Doncic", "Bradley Beal"]
    total points scored per game: 354.3
    payroll: 132.53548


The output of our optimization tells us that with the team listed above, we should expect them to score more than 354 points per game! Unfortunately, that answer makes no sense. What went wrong?

Two major things that went wrong:

1. We didn't account for the fact that there are only 240 minutes in an NBA game (48 minutes times 5 players on the floor)
2. We didn't account for the fact that teams can typically attempt no more than 80 field goals in an NBA game

Let's adjust our code so that we account for these important constraints and see if we get anything more reasonable. We will need to add minutes and field goal attampts into our data import function, and we will need to add a constraint on total minutes and total field goal attempts.


```julia
# first function: read in the player data from the class GitHub repository
function read_in_data(url)
    newtable = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    players  = newtable[1][2]
    salaries = newtable[1][4]./1000000
    mpg      = newtable[1][5]
    fgaG5    = 1.0.*((newtable[1][9]).>5)
    fga      = newtable[1][9]
    ppg      = newtable[1][11]
    return players,salaries,mpg,fgaG5,fga,ppg
end

function SolveModel(players,salary,minutes,field_goals_over5,field_goals,points)
    N = length(salary) 
    
    m = Model(GLPK.Optimizer)

    # define the variables: they are 0 if the player did not make the team, 1 if the player did make the team
    @variable(m, picked[1:N], Bin)

    # categories: 
    @objective(m, Max, sum( points[i] * picked[i] for i in 1:N)) 

    @constraints m begin
        # Constraint 1 - payroll <= 132.6m
        sum(salary[i] * picked[i] for i in 1:N) <= 132.6
        # Constraint 2 - must have exactly 15 players on roster
        sum(picked[i] for i in 1:N) == 15
        # Constraint 3 - total minutes must not exceed 240
        sum(minutes[i] * picked[i] for i in 1:N) <= 240
        # Constraint 4 - total shot attempts must be lower than 80
        sum(field_goals[i] * picked[i] for i in 1:N) <= 80
    end

    # Solve it
    JuMP.optimize!(m);
    pck      = convert(BitArray,JuMP.value.(picked))
    lineup   = players[pck]
    totmin   = sum(minutes[pck])
    payroll  = sum(salary[pck])
    totshots = sum(field_goals[pck])
    points   = JuMP.objective_value(m)
    return lineup,points,totmin,payroll,totshots
end

# call first function (to import data)
players,salaries,minutes,over5fg,fga,pts = read_in_data("https://raw.githubusercontent.com/tyleransom/DScourseS20/master/WebData/playerSalaryStats.csv")

# pass data into second function to get optimal lineup
lineup,total_points,total_minutes,payroll,totshots = SolveModel(players,salaries,minutes,over5fg,fga,pts)
println("team: ",lineup)
println("total points scored per game: ",total_points)
println("payroll: ",payroll)
println("total shots per game: ",totshots)
println("total minutes per game: ",total_minutes)
```

    team: ["Damian Jones", "Javonte Green", "JaVale McGee", "Dwight Howard", "Goga Bitadze", "Jaxson Hayes", "Nerlens Noel", "Patrick Patterson", "Yogi Ferrell", "Chris Boucher", "Tony Bradley", "Mason Plumlee", "Thon Maker", "Christian Wood", "James Harden"]
    total points scored per game: 121.0
    payroll: 84.22254099999999
    total shots per game: 79.89999999999999
    total minutes per game: 239.90000000000003


The results gives us a much mroe reasonable number of 121 points, 80 field goal attempts, and a super-cheap payroll of \\$84m, which is \\$25m lower than the lowest in the NBA right now.

---

# 4. Estimating Linear Regression Coefficients
We can also use JuMP to estimate linear regression coefficients. In this case, we must use the `Ipopt` (pronounced eye-PEE-opt) optimizer. Why? Because our objective function is nonlinear (we are minimizing the sum of the squared residuals) and the optimizers we have used above are only valid for linear objective functions.

### Objective function
The objective function for OLS is

\begin{equation}
\min_{\beta} \sum_{i} (y_i - \beta_0 - \beta_1 x_1 - \beta_2 x_2 - \cdots - \beta_k x_k)^2
\end{equation}

### Variables
The variables in this case are the parameters we want to estimate---the $\beta$'s.

### Code


```julia
using HTTP, JuliaDB, JuMP, Ipopt
function import_auto(url)
    newtable  = csvread(IOBuffer(HTTP.get(url).body), skiplines_begin=0, header_exists=true)
    depvar    = log.(newtable[1][2]) # log price
    indepvars = cat(ones(size(depvar)),newtable[1][3],newtable[1][5],newtable[1][6]; dims=2) # constant, mpg, headroom, trunk
    return depvar,indepvars
end

Y,X = import_auto("https://tyleransom.github.io/teaching/MetricsLabs/auto.csv")


function jumpOLS(Y,X,startval=zeros(size(X,2),1))
    OLS = Model(Ipopt.Optimizer)
    
    # Declare the variables you are optimizing over
    @variable(OLS, b[i=1:size(X,2)], start = startval[i])
    
    # Write your objective function
    @NLobjective(OLS, Min, sum( (Y[i]-sum( X[i,k]*b[k] for k in 1:size(X,2) ))^2 for i in 1:size(X,1) ) )
    
    # Solve the objective function
    JuMP.optimize!(OLS)
    
    SSR = JuMP.objective_value(OLS)
    b_value = JuMP.value.(b)
    println("Objective value: ", SSR)
    println("beta hat = ", b_value)
    println("RMSE = ", sqrt(SSR/(size(X,1)-size(X,2))))
end

jumpOLS(Y,X)

```

    
    ******************************************************************************
    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.:        0
    Number of nonzeros in Lagrangian Hessian.............:       10
    
    Total number of variables............................:        4
                         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...............:        0
            inequality constraints with only lower bounds:        0
       inequality constraints with lower and upper bounds:        0
            inequality constraints with only upper bounds:        0
    
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
       0  5.5361028e+03 0.00e+00 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
       1  8.0341708e+00 0.00e+00 2.37e-14  -1.0 9.28e+00    -  1.00e+00 1.00e+00f  1
    
    Number of Iterations....: 1
    
                                       (scaled)                 (unscaled)
    Objective...............:   2.9676334739952230e-02    8.0341707763946211e+00
    Dual infeasibility......:   2.3739303738078894e-14    6.4268590449501054e-12
    Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
    Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
    Overall NLP error.......:   2.3739303738078894e-14    6.4268590449501054e-12
    
    
    Number of objective function evaluations             = 2
    Number of objective gradient evaluations             = 2
    Number of equality constraint evaluations            = 0
    Number of inequality constraint evaluations          = 0
    Number of equality constraint Jacobian evaluations   = 0
    Number of inequality constraint Jacobian evaluations = 0
    Number of Lagrangian Hessian evaluations             = 1
    Total CPU secs in IPOPT (w/o function evaluations)   =      7.216
    Total CPU secs in NLP function evaluations           =      3.174
    
    EXIT: Optimal Solution Found.
    Objective value: 8.034170776394621
    beta hat = [9.27849, -0.0296502, -0.115783, 0.0247281]
    RMSE = 0.33878292199313836


We can check our answer in R with the following code:
```r
df <- read.csv("https://tyleransom.github.io/teaching/MetricsLabs/auto.csv") %>% as_tibble %>% 
      mutate(logprice = log(price)) %>% 
      drop_na(foreign)
summary(lm(logprice ~ mpg + headroom + trunk, data=df))
```

which gives us

```
Call:
lm(formula = log(price) ~ mpg + headroom + trunk, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6017 -0.2521 -0.1082  0.2104  1.0445 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.278489   0.314134  29.537  < 2e-16 ***
mpg         -0.029650   0.008434  -3.515 0.000775 ***
headroom    -0.115783   0.062605  -1.849 0.068619 .  
trunk        0.024728   0.013857   1.785 0.078667 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3388 on 70 degrees of freedom
Multiple R-squared:  0.2842,	Adjusted R-squared:  0.2535 
F-statistic: 9.263 on 3 and 70 DF,  p-value: 3.082e-05
```

8d79a01fba363cc360ae872107b9ef7655cf27e9
ipynb
Jupyter Notebook
LectureNotes/14-JuMP/JuMPintro.ipynb
brandyherrera/DScourseS21
05f627c2413a61457b97951915368e26549e398a
MIT
2020-01-09T04:00:38.000Z
2022-03-03T12:22:41.000Z
LectureNotes/14-JuMP/JuMPintro.ipynb
brandyherrera/DScourseS21
05f627c2413a61457b97951915368e26549e398a
MIT
2020-01-08T22:46:56.000Z
2020-11-28T20:39:17.000Z
LectureNotes/14-JuMP/JuMPintro.ipynb
brandyherrera/DScourseS21
05f627c2413a61457b97951915368e26549e398a
MIT
2020-01-09T04:00:45.000Z
2022-01-17T12:30:14.000Z
Qwen/Qwen-72B
1. YES
2. YES