Packages Used: No new packages are used in this notebook.
Minimum cost network flow (MCNF) is the most general network model. It can be used to represent a variety of network models that can be solved as an LP. Special procedures more efficient than LP were developed to solve MCNF and Transportation problems, e.g.:
It is usually easier to transform into an LP since current LP solvers are so good, with MCNF just aiding in the formulation of the problem. Special, very efficient procedures are only used for the shortest-path problem (Dijkstra).
Each row i of the arc-incidence matrix and the corresponding net supply represent the coefficients of the flow-balance equation for node i. The LP formulation of the MCNF problem is a follows:
$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{j \in E} c_j x_j \\ \quad \mbox{subject to} \quad \sum_{j \in E} a_{i,j} x_j &\ge& s_i\mbox{,} &\quad i \in V : s_i < 0 \\ \sum_{j \in E} a_{i,j} x_j &=& s_i\mbox{,} &\quad i \in V : s_i \ge 0 \\ x_j &\ge& 0 \mbox{,} \\ \end{eqnarray*} $
where,
$ \begin{eqnarray*} \quad G(V,E) &=& \mbox{network of nodes } V \mbox{ and arcs } E \ c_j &=& \mbox{arc cost per unit flow} \ xj &=& \mbox{arc flow} \ a{i,j} &=& \begin{cases} -1, & \mbox{node } i \mbox{ is source of arc } j \\ \:\:\, 1, & \mbox{node } i \mbox{ is destination of arc } j \end{cases} \ s_i &=& \mbox{net supply of node } i \ &=& \begin{cases} < 0, & \mbox{supply node } i \ = 0, & \mbox{transshipment node } i \
0, & \mbox{demand node } i \mbox{.} \end{cases} \ \end{eqnarray*} $
Note: Although the net flow $\sum a_{i,j} x_j $ cannot exceed supply $s_i$ for a supply node $i$, the flow-balance constraint is $\sum a_{i,j} x_j \ge s_i $ in the formulation because $a_{i,j} = -1$ and $s_i < 0$.
In Ex 4 from Net Mod 1, the LP model was created directly from the interlevel matrix C
. Since the transportation problem is a special case of the MCNF problem, C
can be converted into an incidence matrix, A
, and then solved as an MCNF instance.
using LightGraphs, SimpleWeightedGraphs
C = [8 6 10 9; 9 12 13 7; 14 9 16 5]
m,n = size(C)
W = [zeros(m,m) replace!(C, 0. => sqrt(eps())); zeros(n,m+n)]
s = [-55, -50, -40, 45, 20, 30, 30]
G = SimpleWeightedDiGraph(W)
A = incidence_matrix(G)
c = [weights(G)[src(e),dst(e)] for e in edges(G)]
Matrix([c' 0; A s])
8×13 Matrix{Float64}: 8.0 6.0 10.0 9.0 9.0 12.0 … 7.0 14.0 9.0 16.0 5.0 0.0 -1.0 -1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -55.0 0.0 0.0 0.0 0.0 -1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 -50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 -1.0 -1.0 -40.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 45.0 0.0 1.0 0.0 0.0 0.0 1.0 … 0.0 0.0 1.0 0.0 0.0 20.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 30.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 30.0
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[1:length(E)] )
@objective(model, Min, sum(c[j]*x[j] for j ∈ E) )
for i ∈ V
if s[i] < 0
@constraint(model, sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
else
@constraint(model, sum(A[i,j]*x[j] for j ∈ E) == s[i] )
end
end
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 7 rows, 12 columns, 24 non-zeros 0: obj = 0.000000000e+000 inf = 1.250e+002 (4) 6: obj = 1.100000000e+003 inf = 0.000e+000 (0) * 9: obj = 9.700000000e+002 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(970.0, [5.0, 20.0, 30.0, 0.0, 40.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 30.0])
An more concise means of specifying the MCNF model in JuMP
is using unnamed constaint containers with conditions, where
for i ∈ V
if s[i] < 0
@constraint(model, sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
else
@constraint(model, sum(A[i,j]*x[j] for j ∈ E) == s[i] )
end
end
is replaced with
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[1:length(E)] )
@objective(model, Min, sum(c[j]*x[j] for j ∈ E) )
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 7 rows, 12 columns, 24 non-zeros 0: obj = 0.000000000e+000 inf = 1.250e+002 (4) 6: obj = 1.100000000e+003 inf = 0.000e+000 (0) * 9: obj = 9.700000000e+002 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(970.0, [5.0, 20.0, 30.0, 0.0, 40.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 30.0])
Arc list information is used to create the incidence matrix for the network that can used as the node-balance constraints in the MCNF LP formulation:
i = [1, 1, 2, 2, 3, 4]
j = [2, 3, 4, 5, 5, 5]
w = [2, 3, 4, 5, 1, 3]
s = [-6, -4, 0, 0, 8]
G = SimpleWeightedDiGraph(i,j,w)
A = incidence_matrix(G)
c = [weights(G)[src(e),dst(e)] for e in edges(G)]
Matrix([c' 0; A s])
6×7 Matrix{Int64}: 2 3 4 5 1 3 0 -1 -1 0 0 0 0 -6 1 0 -1 -1 0 0 -4 0 1 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0 1 1 1 8
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[1:length(E)] )
@objective(model, Min, sum(c[j]*x[j] for j ∈ E) )
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 5 rows, 6 columns, 12 non-zeros 0: obj = 0.000000000e+000 inf = 8.000e+000 (1) 4: obj = 4.800000000e+001 inf = 0.000e+000 (0) * 6: obj = 3.400000000e+001 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(34.0, [0.0, 6.0, 0.0, 2.0, 6.0, 0.0])
Upper and lower arc bounds that differ from the default values of $\infty$ and $0$ can easily be added to the MCNF LP formulation.
Arc (4,5) is the only arc with a lower-bound $>0$, and arcs (2,4), (2,5), and (3,5) are the only arcs with upper-bounds $<\infty$:
i = [1, 1, 2, 2, 3, 4]
j = [2, 3, 5, 4, 5, 5]
w = [2, 3, 5, 4, 1, 3]
u = [Inf, Inf, 9, 3, 3, Inf]
l = zeros(length(i))
l[6] = 1
[["i","j","w","u","l"] [i'; j'; w'; u'; l']]
5×7 Matrix{Any}: "i" 1.0 1.0 2.0 2.0 3.0 4.0 "j" 2.0 3.0 5.0 4.0 5.0 5.0 "w" 2.0 3.0 5.0 4.0 1.0 3.0 "u" Inf Inf 9.0 3.0 3.0 Inf "l" 0.0 0.0 0.0 0.0 0.0 1.0
G = SimpleWeightedDiGraph(i,j,w)
@show A = incidence_matrix(G)
c = [weights(G)[src(e),dst(e)] for e in edges(G)]
[["i","j","w","c"] [i'; j'; w'; c']]
A = incidence_matrix(G) = -1 -1 ⋅ ⋅ ⋅ ⋅ 1 ⋅ -1 -1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ -1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ -1 ⋅ ⋅ ⋅ 1 1 1
4×7 Matrix{Any}: "i" 1 1 2 2 3 4 "j" 2 3 5 4 5 5 "w" 2 3 5 4 1 3 "c" 2 3 4 5 1 3
Note: In the above arc lists, edge 2 => 5 is defined before edge 2 => 4, but the column for edge 2 => 4 is first in A
. As a result, the weight vector w
used to define the network g
does not correspond to the order that the arcs appear in A
and, as a result, the arc cost vector c
needs to be explicitly determined by iterating over each edge (e.g., [weights(g)[src(e),dst(e)] for e in edges(g)]
). Since this cannot be done for the upper and lower arc bounds, instead, an index array idx
is created by iterating over each edge and then used to reorder the arc bounds:
s = [-6, -4, 0, 0, 8]
G = SimpleWeightedDiGraph(i, j, 1:length(i))
idx = [weights(G)[src(e),dst(e)] for e in edges(G)] # Edge order index
A = incidence_matrix(G)
c, u, l = w[idx], u[idx], l[idx]
Matrix([c'; u'; l'; A])
8×6 Matrix{Float64}: 2.0 3.0 4.0 5.0 1.0 3.0 Inf Inf 3.0 9.0 3.0 Inf 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 1.0 1.0 1.0
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, l[j] <= x[j = 1:length(E)] <= u[j] )
@objective(model, Min, sum(c[i]*x[i] for i ∈ E) )
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 5 rows, 6 columns, 12 non-zeros 0: obj = 3.000000000e+000 inf = 8.000e+000 (2) 4: obj = 5.000000000e+001 inf = 0.000e+000 (0) * 5: obj = 4.100000000e+001 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(41.0, [1.0, 3.0, 1.0, 4.0, 3.0, 1.0])
Matrix([xᵒ'; A])
6×6 Matrix{Float64}: 1.0 3.0 1.0 4.0 3.0 1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 1.0 1.0 1.0
Matrix(A)
5×6 Matrix{Int64}: -1 -1 0 0 0 0 1 0 -1 -1 0 0 0 1 0 0 -1 0 0 0 1 0 0 -1 0 0 0 1 1 1
The shortest path problem can be formulated as an MCNF problem by adding a supply/demand of one to starting/terminal nodes and a net supply of $0$ for the other (transshipment) nodes in the network:
i = [1, 1, 2, 2, 3, 3, 4, 4, 5]
j = [2, 3, 3, 4, 4, 5, 5, 6, 6]
w = [4, 2, 1, 5, 8, 10, 2, 6, 3]
G = SimpleWeightedDiGraph(i,j,w)
s = [-1, 0, 0, 0, 0, 1]
A = incidence_matrix(G)
c = [weights(G)[src(e),dst(e)] for e in edges(G)]
Matrix([c'; A])
7×9 Matrix{Int64}: 4 2 1 5 8 10 2 6 3 -1 -1 0 0 0 0 0 0 0 1 0 -1 -1 0 0 0 0 0 0 1 1 0 -1 -1 0 0 0 0 0 0 1 1 0 -1 -1 0 0 0 0 0 0 1 1 0 -1 0 0 0 0 0 0 0 1 1
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[1:length(E)] )
@objective(model, Min, sum(c[j]*x[j] for j ∈ E) )
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 6 rows, 9 columns, 18 non-zeros 0: obj = 0.000000000e+000 inf = 1.000e+000 (1) 5: obj = 1.500000000e+001 inf = 0.000e+000 (0) * 6: obj = 1.400000000e+001 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(14.0, [1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0])
Matrix([xᵒ'; A])
7×9 Matrix{Float64}: 1.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0
Additional constraints: The only reason to use the MCNF formulation for the shortest path problem would be if there needed to be additional constraints like arc bounds that need to be included in the model of the problem. For example, to find the shortest path in the directed network that always contains the arc from node 2 to node 3, solve as MCNF with LB of arc (2,3) = 1 to force its use:
l = zeros(length(i))
l[3] = 1
G = SimpleWeightedDiGraph(i, j, 1:length(i))
idx = [weights(G)[src(e),dst(e)] for e in edges(G)]
A = incidence_matrix(G)
c, l = w[idx], l[idx]
Matrix([c'; l'; A])
8×9 Matrix{Float64}: 4.0 2.0 1.0 5.0 8.0 10.0 2.0 6.0 3.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, x[i = 1:length(E)] >= l[i] )
@objective(model, Min, sum(c[j]*x[j] for j ∈ E) )
@constraint(model, [i ∈ V; s[i] < 0], sum(A[i,j]*x[j] for j ∈ E) >= s[i] )
@constraint(model, [i ∈ V; s[i] >= 0], sum(A[i,j]*x[j] for j ∈ E) == s[i] )
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 6 rows, 9 columns, 18 non-zeros 0: obj = 1.000000000e+000 inf = 3.000e+000 (3) 3: obj = 1.900000000e+001 inf = 0.000e+000 (0) * 5: obj = 1.800000000e+001 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(18.0, [1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
Matrix([xᵒ'; A])
7×9 Matrix{Float64}: 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0
If all of the arcs in a network have upper bounds $<\infty$, then determining the maximum amount of flow that can occur between source nodes and sink (or demand) nodes can be used to model many different problems.
In the case that there is a single source and a single sink node, an arc can be added from the sink to the source with a cost of $1$ and an upper bound of $\infty$. All of the other nodes have $0$ net supply and $0$ cost. The resulting flow on the added arc (18) represents the maximum flow in the network.
using LightGraphs, SimpleWeightedGraphs
i = [1, 1, 2, 2, 3, 3, 4, 4, 5, 6]
j = [2, 3, 3, 4, 4, 5, 5, 6, 6, 1]
w = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
u = [14, 12, 1, 5, 8, 10, 2, 16, 13, Inf]
G = SimpleWeightedDiGraph(i,j,1:length(i))
idx = [weights(G)[src(e),dst(e)] for e in edges(G)]
A = incidence_matrix(G)
c, u = w[idx], u[idx]
Matrix([c'; u'; A])
8×10 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 14.0 12.0 1.0 5.0 8.0 10.0 2.0 16.0 13.0 Inf -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 -1.0
using JuMP, GLPK
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[j = 1:length(E)] <= u[j] )
@objective(model, Max, sum(c[i]*x[i] for i ∈ E) )
@constraint(model, [i ∈ V], sum(A[i,j]*x[j] for j ∈ E) == 0 )
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 6 rows, 10 columns, 20 non-zeros * 0: obj = -0.000000000e+000 inf = 0.000e+000 (1) * 9: obj = 1.800000000e+001 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(18.0, [6.0, 12.0, 1.0, 5.0, 8.0, 5.0, 0.0, 13.0, 5.0, 18.0])
Matrix([xᵒ'; u'; A])
8×10 Matrix{Float64}: 6.0 12.0 1.0 5.0 8.0 5.0 0.0 13.0 5.0 18.0 14.0 12.0 1.0 5.0 8.0 10.0 2.0 16.0 13.0 Inf -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 -1.0
In this example of oil pipeline distribution network from Taha, when there are multiple source nodes, an additional node can be added as a new single-source node that is connected to the original source nodes with arcs of $0$ costs and upper bounds of $\infty$; similarly, a new single sink node can be added if there are multiple sink nodes:
It can be convienent to enter each arc as a row in a matric when entering network data manually:
IJU = [
1 4 20
2 4 10
2 5 20
2 6 50
3 5 15
4 5 20
4 6 10
4 7 10
5 6 30
5 8 30
6 7 50
6 8 20
7 9 Inf
8 9 Inf
9 10 Inf
10 1 Inf
10 2 Inf
10 3 Inf]
i = Int.(IJU[:,1])
j = Int.(IJU[:,2])
u = IJU[:,3]
w = zeros(length(i))
w[end-3] = 1
[i'; j'; u'; w']
4×18 Matrix{Float64}: 1.0 2.0 2.0 2.0 3.0 4.0 … 7.0 8.0 9.0 10.0 10.0 10.0 4.0 4.0 5.0 6.0 5.0 5.0 9.0 9.0 10.0 1.0 2.0 3.0 20.0 10.0 20.0 50.0 15.0 20.0 Inf Inf Inf Inf Inf Inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
G = SimpleWeightedDiGraph(i,j,1:length(i))
idx = [weights(G)[src(e),dst(e)] for e in edges(G)]
A = incidence_matrix(G)
c, u = w[idx], u[idx]
Matrix([c'; u'; A])
12×18 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 1.0 0.0 0.0 0.0 20.0 10.0 20.0 50.0 15.0 20.0 Inf Inf Inf Inf Inf Inf -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 -1.0 0.0 -1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 -1.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 1.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0
V = 1:size(A,1)
E = 1:size(A,2)
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x[j = 1:length(E)] <= u[j] )
@objective(model, Max, sum(c[i]*x[i] for i ∈ E) )
@constraint(model, [i ∈ V], sum(A[i,j]*x[j] for j ∈ E) == 0 )
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
xᵒ = value.(x)
TCᵒ, xᵒ
GLPK Simplex Optimizer 5.0 10 rows, 18 columns, 36 non-zeros * 0: obj = -0.000000000e+000 inf = 0.000e+000 (1) * 18: obj = 1.100000000e+002 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(110.0, [20.0, 10.0, 15.0, 50.0, 15.0, 10.0, 10.0, 10.0, 10.0, 30.0, 50.0, 20.0, 60.0, 50.0, 110.0, 110.0, 75.0, 15.0])
Matrix([xᵒ'; u'; A])
12×18 Matrix{Float64}: 20.0 10.0 15.0 50.0 15.0 10.0 10.0 … 50.0 110.0 110.0 75.0 15.0 20.0 10.0 20.0 50.0 15.0 20.0 10.0 Inf Inf Inf Inf Inf -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 -1.0 0.0 -1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 -1.0 -1.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0
Production-inventory planning models are the main use of MP in industry. They provide a means to make complex decisions over a rolling planning horizon. The following production-inventory planning model is for a single product, with a rolling horizon of t periods (e.g., months):
$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i=1}^{m} \sum_{j=1}^{t} c_i^p x_{ij} + \sum_{i=1}^{m} \sum_{j=2}^{t+1} c_i^i y_{ij} \\ \quad \mbox{subject to} \quad -x_{ij} + x_{(i+1)j} - y_{ij} + y_{i(j+1)} &=& 0 \mbox{,} &\quad i = 1, \dots, m-1; j = 1, \dots, t &\quad (a) \\ -x_{m,j} - y_{m,j} + y_{m(j+1)} &=& d_j \mbox{,} &\quad j = 1, \dots, t &\quad (b) \\ x_{ij} &\le& u_{ij} \mbox{,} &\quad i = 1, \dots, m; j = 1, \dots, t \\ y_{i,1} &=& y_i^0 \mbox{,} &\quad i = 1, \dots, m \\ y_{i(t+1)} &=& y_i^t \mbox{,} &\quad i = 1, \dots, m \\ x, y &\ge& 0 \mbox{,} \\ \end{eqnarray*} $
where,
$ \begin{eqnarray*} \quad m &=& \mbox{number of production stages} \\ t &=& \mbox{number of periods of production} \\ c_i^p &=& \mbox{production cost (dollar/ton) in stage } i \\ x_{ij} &=& \mbox{production (ton) at stage } i \mbox{ in period } j \\ c_i^i &=& \mbox{inventory cost (dollar/ton) in stage } i \\ y_{ij} &=& \mbox{stage-} i \mbox{ inventory (ton) from period } j-1 \mbox{ to } j \\ d_j &=& \mbox{demand (ton) in period } j \\ u_{ij} &=& \mbox{production capacity (ton) of stage } i \mbox{ in period } j \\ y_i^0 &=& \mbox{initial inventory (ton) of stage } i \\ y_i^t &=& \mbox{final inventory (ton) of stage } i \mbox{.} \\ \end{eqnarray*} $
In the objective function, inventory costs are summed from 2 to t + 1 because the period 1 initial inventory cost was accounted for in period 0 prior to period 1 (rolling horizon). Constraints (a) and (b) in the model are the flow-balance equations for each node, while the remaining constraints represent lower- and upper-bounds on the decision variables x and y.
A single product is produced in two stages. The planning period is three months, and the model is run each month using current demand forecasts and scheduled production capacity. The following shows the production-inventory planning network:
In the following MCNF formulation, the flow-balance constraints are directly implemented in the model instead of first determining the arc-incidence matrix:
U = [50 0 50; # capacity of stage m in period t (ton)
60 0 0];
d = [20, 10, 15] # demand in period t (ton)
cp = [200, 800] # production cost in stage m ($/ton)
ci = [5, 25] # inventory cost for stage m ($/ton)
y0 = [0, 0] # initial inventory of stage m (ton)
yt = [35, 0] # final inventory of stage m (ton)
m = length(cp) # number of production stages = 2
t = length(d) # number of periods of production = 3
Ylb = [y0 zeros(m,t-1) yt] # same LB and UB to fix inventory level
Yub = [y0 repeat([Inf],m,t-1) yt]
2×4 Matrix{Float64}: 0.0 Inf Inf 35.0 0.0 Inf Inf 0.0
using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, 0 <= X[i=1:m, j=1:t] <= U[i,j] )
@variable(model, Ylb[i,j] <= Y[i=1:m, j=1:t+1] <= Yub[i,j] )
@objective(model, Min, sum(cp[i]*X[i,j] for i=1:m, j=1:t ) +
sum(ci[i]*Y[i,j] for i=1:m, j=2:t+1 ) )
for j = 1:t
for i = 1:m-1
@constraint(model, -X[i,j] + X[i+1,j] - Y[i,j] + Y[i,j+1] == 0 ) # (a)
end
@constraint(model, -X[m,j] - Y[m,j] + Y[m,j+1] == -d[j] ) # (b)
end
print(model)
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
Xᵒ = value.(X)
Yᵒ = value.(Y)
TCᵒ, Xᵒ, Yᵒ
GLPK Simplex Optimizer 5.0 6 rows, 14 columns, 21 non-zeros 0: obj = 1.750000000e+002 inf = 8.000e+001 (4) 5: obj = 5.317500000e+004 inf = 0.000e+000 (0) * 6: obj = 5.317500000e+004 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(53175.0, [45.0 0.0 35.0; 45.0 0.0 0.0], [0.0 0.0 0.0 35.0; 0.0 25.0 15.0 0.0])
The following shows the optimal production and inventory flows (in blue) in the production-inventory planning network. Note that the flow is balanced at each node; that is, the flow into a node equals the combined production, inventory, and demand flow out of the node.
Optimal production flows (and required demands):
Matrix([Xᵒ; d'])
3×3 Matrix{Float64}: 45.0 0.0 35.0 45.0 0.0 0.0 20.0 10.0 15.0
Optimal inventory flows:
Matrix(Yᵒ)
2×4 Matrix{Float64}: 0.0 0.0 0.0 35.0 0.0 25.0 15.0 0.0
A single product is produced in a two-stage production process. Stage one has a capacity to produce 40 tons of the product per month, and there is a cost of \$200 for each ton of product produced. For stage two, capacity is 80 tons per month, and there is a cost of \\$800 per ton. Due to scheduled maintenance, no capacity will be available during month five for stage 1 and during months three and five for stage 2. Demand for the next six months is estimated to be 25, 15, 10, 50, 25, and 15 tons per month. Currently, five tons of material is in storage at stage 1 and 35 tons at stage 2; at the end of the planning period, ten tons of material should be in storage at stage 1 and 15 at stage 2. Up to 30 tons of material per month can be stored at each stage. Assuming that the inventory costs are \$12 and \\$46 per ton-month, determine the production plan that minimizes total costs over the planning horizon and what the total production costs will be, and what the total inventory costs will be.
Note: The maximum inventory storage constraint can be modeled as an upper bound on y.
U = [40 40 40 40 0 40; # capacity of each stage (ton)
80 80 0 80 0 80]
d = [25, 15, 10, 50, 25, 15] # demand in each period (ton)
cp = [200, 800] # variable production cost for each stage ($/ton)
ci = [12, 46] # inventory cost for each stage ($/ton)
y0 = [5, 25] # initial inventory of each stage (ton)
yt = [10, 15] # final inventory of each stage (ton)
m = length(cp) # number of production stages
t = length(d) # number of periods of production
Ylb = [y0 zeros(m,t-1) yt] # same LB and UB to fix inital and final inventory levels
Yub = [y0 repeat([30],m,t-1) yt] # max inventory storage level
2×7 Matrix{Int64}: 5 30 30 30 30 30 10 25 30 30 30 30 30 15
using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, 0 <= X[i=1:m, j=1:t] <= U[i,j] )
@variable(model, Ylb[i,j] <= Y[i=1:m, j=1:t+1] <= Yub[i,j] )
@objective(model, Min, sum(cp[i]*X[i,j] for i=1:m, j=1:t ) +
sum(ci[i]*Y[i,j] for i=1:m, j=2:t+1 ) )
for j = 1:t
for i = 1:m-1
@constraint(model, -X[i,j] + X[i+1,j] - Y[i,j] + Y[i,j+1] == 0 ) # (a)
end
@constraint(model, -X[m,j] - Y[m,j] + Y[m,j+1] == -d[j] ) # (b)
end
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
optimize!(model)
TCᵒ = objective_value(model)
Xᵒ = value.(X)
Yᵒ = value.(Y)
TCᵒ, Xᵒ, Yᵒ
GLPK Simplex Optimizer 5.0 12 rows, 26 columns, 42 non-zeros 0: obj = 8.100000000e+002 inf = 1.450e+002 (7) 13: obj = 1.344700000e+005 inf = 0.000e+000 (0) * 16: obj = 1.343000000e+005 inf = 0.000e+000 (0) OPTIMAL LP SOLUTION FOUND
(134300.0, [0.0 25.0 … 0.0 40.0; 0.0 30.0 … 0.0 30.0], [5.0 5.0 … 0.0 10.0; 25.0 0.0 … 0.0 15.0])
Optimal production flows (and required demands):
Matrix([Xᵒ; d'])
3×6 Matrix{Float64}: 0.0 25.0 30.0 40.0 0.0 40.0 0.0 30.0 0.0 70.0 0.0 30.0 25.0 15.0 10.0 50.0 25.0 15.0
Optimal inventory flows:
Matrix(Yᵒ)
2×7 Matrix{Float64}: 5.0 5.0 0.0 30.0 0.0 0.0 10.0 25.0 0.0 15.0 5.0 25.0 0.0 15.0
Calculate total production cost and total inventory cost:
TCp, TCi = sum(sum(cp.*Xᵒ)), sum(sum(ci.*Yᵒ))
(131000.0, 4510.0)
Sum of component costs exceeds total cost from solver because it includes first-period (time 0) inventory costs:
TCp+TCi, TCᵒ
(135510.0, 134300.0)
Exclude first period from inventory costs:
TCi = sum(sum(ci.*Yᵒ[:,2:end]))
3300.0
Total component costs now match total cost from solver:
TCp+TCi, TCᵒ
(134300.0, 134300.0)