# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Julia 1.5.4
#     language: julia
#     name: julia-1.5
# ---

# ### Problem 5.1 - Capacitated transportation problem revisited
# *(Marginal costs via dual calculation)*

# Adding necessary packages:

# +
# using Pkg
# Pkg.add("GLPK") 
# -

## Packages
using JuMP, Cbc, GLPK # We need GLPK for obtaining duals

# We start by defining the problem data

# +
nS = 3 # Number of suppliers
nD = 3 # Number of demand points
nP = 2 # Number of products

S = 1:nS
D = 1:nD
P = 1:nP

A = [] # Set of arcs, we include arcs from all suppliers to all demand points
for s in S, d in D
    push!(A, (s,d))
end

# +
costs = zeros(nP,nS,nD) # Cost of transporting one unit of product d from supplier s to demand point d 
costs[1,:,:] = [5 5 Inf; 8 9 7; Inf 10 8]
costs[2,:,:] = [Inf 18 Inf; 15 12 14; Inf 20 Inf]

sup = [80 400; 200 1500; 200 300]
dem = [60 300; 100 1000; 200 500]

cap = [Inf 300 0; 300 700 600; 0 Inf Inf];

# +
## Define model using GLPK
model = Model(GLPK.Optimizer);

## Variables
@variable(model, x[a in A, p in P; costs[p,a[1],a[2]] < Inf] >= 0);

## OF
@objective(model, Min, sum(costs[p,a[1],a[2]]*x[a,p] for p in P, a in A if costs[p,a[1],a[2]] < Inf));

## Constraints
@constraint(model, sup_con[s in S, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[1] == s) <= sup[s,p]);       # sum of everything that leaves supplier s
@constraint(model, dem_con[d in D, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[2] == d) >= dem[d,p]);       # sum of everything that arrives at demand d
@constraint(model, cap_con[a in A; cap[a[1],a[2]] < Inf], sum(x[a,p] for p in P if costs[p,a[1],a[2]] < Inf) <= cap[a[1],a[2]]);# arc capacity constraints  
# -

set_silent(model) # Actually works with GLPK
optimize!(model)
status = termination_status(model)
println(status)

## Saving the optimal value
obj = objective_value(model)

# Function ``dual()`` in the ``JuMP`` library gives the value of the dual variable associated to the constraint at the optimal solution, in other words the *marginal costs*. Here we need using the elemnt-wise operator ``.`` as we have multiple constraints (check the domains). The *marginal costs* value stands for how much adding one unit more to the constraint's RHS (in the case is a $\leq$ constraint) impacts the final result.<br/>
# <br/>
# One interpretation for the *marginal costs* in this problem is how much the company is willing to pay for expanding the supplies' or the arcs' capacity (depending on if we're analysing the dual of the supplies or the arcs constraints).

## Computing the duals to infer the marginal costs
mc_supply = dual.(sup_con);
mc_arcs = dual.(cap_con);
for s in S, p in P
    println("The marginal costs for the supply $s for the product $p is: $(mc_supply[s,p])")
end
for a in A
    if cap[a[1],a[2]] < Inf
        println("The marginal costs for the arc $a is: $(mc_arcs[a])")
    end
end
for a in A, p in P
    if costs[p,a[1],a[2]] < Inf
        println("The value of $(x[a,p]) is $(value(x[a,p]))")
    end
end

## Raising the supply availability of product 1 at the first supply node by 1
sup[1,1] = sup[1,1] + 1;

# +
## Define model using GLPK
model = Model(GLPK.Optimizer);

## Variables
@variable(model, x[a in A, p in P; costs[p,a[1],a[2]] < Inf] >= 0);

## OF
@objective(model, Min, sum(costs[p,a[1],a[2]]*x[a,p] for p in P, a in A if costs[p,a[1],a[2]] < Inf));

## Constraints
@constraint(model, sup_con[s in S, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[1] == s) <= sup[s,p]);       # sum of everything that leaves supplier s
@constraint(model, dem_con[d in D, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[2] == d) >= dem[d,p]);       # sum of everything that arrives at demand d
@constraint(model, cap_con[a in A; cap[a[1],a[2]] < Inf], sum(x[a,p] for p in P if costs[p,a[1],a[2]] < Inf) <= cap[a[1],a[2]]);# arc capacity constraints             
# -

set_silent(model)
optimize!(model)
status = termination_status(model)
println(status)

## New optimal value
new_obj = objective_value(model)

## Decrease in the optimal value
new_obj - obj

mc_supply[1,1]

# The marginal cost predicted the change in objective value correctly. Let's now try changing another bound.

## Back to the original supply availability
sup[1,1] = sup[1,1] - 1;

sup[1,1]

## Increasing the arc capacity by 1 for the arc from supplier 2 to demand node 2
cap[2, 2] = cap[2, 2] + 1;

# +
## Define model using GLPK
model = Model(GLPK.Optimizer);

## Variables
@variable(model, x[a in A, p in P; costs[p,a[1],a[2]] < Inf] >= 0);

## OF
@objective(model, Min, sum(costs[p,a[1],a[2]]*x[a,p] for p in P, a in A if costs[p,a[1],a[2]] < Inf));

## Constraints
@constraint(model, sup_con[s in S, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[1] == s) <= sup[s,p]);       # sum of everything that leaves supplier s
@constraint(model, dem_con[d in D, p in P], sum(x[a,p] for a in A if costs[p,a[1],a[2]] < Inf && a[2] == d) >= dem[d,p]);       # sum of everything that arrives at demand d
@constraint(model, cap_con[a in A; cap[a[1],a[2]] < Inf], sum(x[a,p] for p in P if costs[p,a[1],a[2]] < Inf) <= cap[a[1],a[2]]);# arc capacity constraints                            
# -

set_silent(model)
optimize!(model)
status = termination_status(model)
println(status)

## New optimal value
new_obj = objective_value(model)

new_obj-obj

mc_arcs[(2, 2)]

# Turns out the marginal cost that we calculated did not predict this change correctly. Comparing the solution below, we notice that one unit of product 1 to demand point 2 is now transported from supplier 2 instead of supplier 3 due to the increased capacity of arc 2->2. This shows that you should always check whether you can apply marginal costs. If the optimal basis changes because of a change in $b$, you can't apply marginal costs.

## Computing the duals to infer the marginal costs
mc_supply = dual.(sup_con);
mc_arcs = dual.(cap_con);
for s in S, p in P
    println("The marginal costs for the supply $s for the product $p is: $(mc_supply[s,p])")
end
for a in A
    if cap[a[1],a[2]] < Inf
        println("The marginal costs for the arc $a is: $(mc_arcs[a])")
    end
end
for a in A, p in P
    if costs[p,a[1],a[2]] < Inf
        println("The value of $(x[a,p]) is $(value(x[a,p]))")
    end
end

# ### Problem 5.5 - Complementary slackness
# Recall the paint factory problem introduced in Section 1.2.1. of the lecture notes

A = [6 4; 1 2; -1 1; 0 1]
b = [24; 6; 1; 2]
c = [5; 4];

using Cbc, JuMP

# #### Solving the primal and dual problems
# Formulate and solve the two problems to obtain their optimal solutions

#%
primalmodel = Model(Cbc.Optimizer)
@variable(primalmodel, x[1:2] >= 0)
@objective(primalmodel, Max, sum(c.*x))
@constraint(primalmodel, [m in 1:4], sum(A[m,:].*x) <= b[m])
optimize!(primalmodel)
value.(primalmodel[:x])
#%

#%
dualmodel = Model(Cbc.Optimizer)
@variable(dualmodel, p[1:4] >= 0)
@objective(dualmodel, Min, sum(b.*p))
@constraint(dualmodel, [n in 1:2], sum(A[:,n].*p) >= c[n])
optimize!(dualmodel)
value.(dualmodel[:p])
#%

# #### Complementary slackness
# Verify the optimality of your solutions using complementary slackness

# +
#%
x_opt = value.(primalmodel[:x])
p_opt = value.(dualmodel[:p]) 
lhs_primal = A*x_opt
lhs_dual = A'*p_opt

println("Primal feasibility:")
for m in 1:4
    println("$(round(lhs_primal[m],digits=2)) <= $(b[m])")
end

println("\nDual feasibility:")
for n in 1:2
    println("$(round(lhs_dual[n],digits=2)) >= $(c[n])")
end

println("\nComplementary slackness:")
for m in 1:4
    println("$(round((lhs_primal[m] - b[m])*p_opt[m],digits=2)) = 0")
end
for n in 1:2
    println("$(round((lhs_dual[n] - c[n])*x_opt[n],digits=2)) = 0")
end
#%
# -



StarcoderdataJupyterScriptsDedupFiltered