```julia
using JuMP, Cbc, Plots
```

# MS-E2121 - Linear optimization
## Exercise session 1

First, we warm up with two small optimization problems. After that, you are expected to solve a small transportation problem.
#' 
### Demo exercise a)
\begin{align}
    \max \ & x_1 + 2x_2 + 5x_3 \\
    \text{s.t. } & x_1 - x_2 -3x_3 && \geq 5 \\
    & x_1 + 3x_2 - 7x_3 && \leq 10 \\
    & x_1 && \leq 10 \\
    & x_i \geq 0, \ \forall i \in \{1,2,3\}
\end{align}


```julia
model_a = Model(Cbc.Optimizer);
```


```julia
@variable(model_a, x[1:3] >= 0)
```


```julia
@objective(model_a, Max, x[1] + 2*x[2] + 5*x[3])
```


```julia
@constraint(model_a, x[1] - x[2] - 3*x[3] >= 5)
@constraint(model_a, x[1] + 3*x[2] - 7*x[3] <= 10)
@constraint(model_a, x[1] <= 10)
```


```julia
println(model_a)
```


```julia
optimize!(model_a)
```

#### Extracting the solutions:
- termination_status(model) gives the reason for stopping the optimization process.
- objective_value(model) gives the final objective value
- value(model[:variable]) gives the final value of the given variable in the model. We could also do value(x) instead, since we only have one model. If we have multiple models each with a variable named x, value(x) would return the values corresponding to the variable that was defined last.


```julia
println(termination_status(model_a))
println(objective_value(model_a))
println(value.(model_a[:x]))        # the dot notation value.() means the function value() element-wise to collection type input.
```


```julia
if termination_status(model_a) == MOI.OPTIMAL
    println("Found optimal solution $(objective_value(model_a)) at x=$(value.(model_a[:x]))")
else
    println("Optimization terminated with status $(termination_status(model_a))")
end
```

### Demo exercise b)
\begin{align}
    \max \ & 2x_1 + 4x_2 \\
    \text{s.t. } & x_1 + x_2 && \leq 5 \\
    & -x_1 + 3x_2 && \leq 1 \\
    & x_1 && \leq 5 \\
    & x_2 && \leq 5 \\
    & x_i \geq 0, \ \forall i \in \{1,2\}
\end{align}

$~$
$~$

While it is possible to write out the constraints individually, it is sometimes preferable to use the matrix form:
\begin{align}
    \max \ & c^\top x \\
    \text{s.t. } & Ax && \leq b \\
    & x_i \geq 0, \ \forall i \in \{1,2\},
\end{align}

$~$
$~$

where for this problem, $c = \begin{bmatrix}
2 \\
4 
\end{bmatrix}$, $A = \begin{bmatrix}
1 & 1 \\
-1 & 3 \\
1 & 0 \\
0 & 1 
\end{bmatrix}$ and $b = \begin{bmatrix}
5 \\
1 \\
5 \\
5 
\end{bmatrix}$


```julia
c = [2 4]
A = [1 1; -1 3; 1 0; 0 1]
b = [5 1 5 5]'
model_b = Model(Cbc.Optimizer);
```


```julia
@variable(model_b, x[i=1:2] >= 0)
```


```julia
@objective(model_b, Max, sum(c[i]*x[i] for i in 1:2))
```


```julia
@constraint(model_b, [j in 1:4], sum(A[j,i]*x[i] for i in 1:2) .<= b[j])
```


```julia
println(model_b)
```


```julia
optimize!(model_b)
```


```julia
println(termination_status(model_b))
println(objective_value(model_b))
println(value.(model_b[:x]))
```


```julia
if termination_status(model_b) == MOI.OPTIMAL
    println("Found optimal solution $(objective_value(model_b)) at x=$(value.(model_b[:x]))")
else
    println("Optimization terminated with status $(termination_status(model_b))")
end
```

We can also plot the constraints, the corresponding feasible region and the objective function at the optimal point to verify our result. We can see that the objective function indeed touches a vertex of the feasible region at x = (3.5, 1.5).


```julia
function plot_model(A, b, c; sol=nothing, xlim=5, ylim=5)
    
    # Check problem dimensions and throw an error if not 2D
    m,n = size(A)
    if n != 2
        throw(ArgumentError("Model is not 2-dimensional"))
    end
    
    # Some ranges to help with plotting, these are similar to what you may have used in Matlab
    x1_range = collect(0:0.1:5)
    x2_range = collect(0:0.1:5)
    
    # Create empty plot
    plot(xlim=(-1,6), ylim=(-1,6), xlabel="x1", ylabel="x2")
    
    # Add lines Ax=b to the plot 
    for con in 1:m
        a_x1 = A[con,1]
        a_x2 = A[con,2]
        if a_x1!=0 && a_x2!=0
            plot!(x1_range, (-a_x1*x1_range.+b[con])/a_x2, label="$(a_x1==1 ? "" : a_x1)x1+$(a_x2==1 ? "" : a_x2)x2=$(b[con])")
        elseif a_x1!=0 && a_x2==0
            vline!([b[con]], label="$(a_x1==1 ? "" : a_x1)x1=$(b[con])")
        elseif a_x1==0 && a_x2!=0
            hline!([b[con]], label="$(a_x2==1 ? "" : a_x2)x2=$(b[con])")
        end
    end
    vline!([0], label="x1=0")
    hline!([0], label="x2=0")

    # Collect the points in the feasible region using a grid
    feasible_region = []
    for x1 in x1_range
        for x2 in x2_range
            if all(A*[x1;x2].<=b)
                push!(feasible_region, [x1, x2])
            end
        end
    end

    # Turn the array of arrays into a 2D-array and plot the points
    feasible_region = hcat(feasible_region...)
    scatter!(feasible_region[1,:], feasible_region[2,:], ms=1, label=false)
    
    # If objective function value was given as an argument, plot the objective function with that value
    if !isnothing(sol)
        plot!(x1_range, (-c[1]*x1_range.+sol)/c[2], linestyle=:dash, linewidth=2, color=:black, label="Objective function")
    end
    
    # This is how you would add a title to your plot
    title!("")
end

plot_model(A,b,c; sol=13)
```

### Student exercise

_Source: Julia Programming for Operations Research 2$^{nd}$ edition (Online)._
_Available at: http://www.chkwon.net/julia - published by Changhyun Kwon_

$~$
$~$

Consider the following network where the combined supplies from the Austin and Buffalo nodes need to meet the demands coming from Chicago, Denver, and Erie. The supplies available are represented by $S$ and the demands by $D$, the costs of transporting in each arc connecting supply and demand nodes are shown as $c$.


Solve the transportation problem to the minimum cost. The data is given, you will need to:

- create decision variables (What are the decisions in this problem? The amount transported from city i to city j.)
- formulate the objective function (What are we minimizing? The total transportation cost.)
- write constraints (What constraints need to be satisfied? Suppliers cannot send out more than their supply capacity and the demand in each node must be fulfilled.)

The correct optimal cost is 251.

#### Mathematical Fomularion of Transportation Problem
##### Algebaric (Symbolic) Form
For sets $I = \{ \text{Austin}, \text{Buffalo} \}$ and $J = \{ \text{Chicago}, \text{Denver}, \text{Erie} \}$ of the supply and demand cities, respectively,
\begin{align}
    \min. \ \ & z = \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} \\
    \text{s.t. } & \sum_{j \in J} x_{ij} \leq C_i, \forall i \\
    & \sum_{i \in I} x_{ij} \geq D_j, \forall j \\
    & x_{ij} \geq 0, \forall i \in I, \forall j \in J
\end{align}

$~$
$~$

where for this particular problem, 
$c = \begin{bmatrix}
10 & 7 & 9 \\
4 & 9 & 8 
\end{bmatrix}$, 
$C = \begin{bmatrix}
15 \\
25
\end{bmatrix}$, and 
$D = \begin{bmatrix}
15 \\
12 \\
13 
\end{bmatrix}$.

#### Loose Format


```julia
using JuMP, Cbc

# Supply nodes and their capacities
supply_nodes = ["Austin" "Buffalo"]
S = [15 25]
# Demand nodes and corresponding demands
demand_nodes = ["Chicago" "Denver" "Erie"]
D = [15 12 13]
# Transportation unit costs
c = [10 7 9; 4 9 8]

# Converting arrays to dictionaries
s_dict = Dict(supply_nodes .=> S)
d_dict = Dict(demand_nodes .=> D)
c_dict = Dict( (supply_nodes[i], demand_nodes[j]) => c[i,j]
          for i in 1:length(supply_nodes), j in 1:length(demand_nodes) );
```


```julia
# Set the model
model_transport = Model(Cbc.Optimizer);
```


```julia
# Variables
@variable(model_transport, x[supply_nodes, demand_nodes] ≥ 0);
```


```julia
# Formulate objective function
@objective(model_transport, Min, sum(c_dict[i,j]*x[i,j] for i in supply_nodes for j in demand_nodes));
```


```julia
# Formulate constraints
@constraint(model_transport, supply_cap[i in supply_nodes], sum(x[i,j] for j in demand_nodes) ≤ s_dict[i]);
@constraint(model_transport, demand_limit[j in demand_nodes], sum(x[i,j] for i in supply_nodes) ≥ d_dict[j]);
```


```julia
# Print the model and optimize it
println(model_transport)
optimize!(model_transport)

x_star = value.(x)
for s in supply_nodes, d in demand_nodes
    println("from $s to $d: $(x_star[s,d]) units")
end
```

584925be87996a4149b4c32bfa9e4a26ca99ad61
ipynb
Jupyter Notebook
exercises/ex 01/exercise1_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
exercises/ex 01/exercise1_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
exercises/ex 01/exercise1_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
Qwen/Qwen-72B
1. YES
2. YES