|
1 | 1 | """
|
2 |
| - mincost_flow(graph, node_demand, edge_capacity, edge_cost, solver,<keyword arguments>) |
| 2 | + function mincost_flow(g::AbstractGraph, |
| 3 | + node_demand::AbstractVector, |
| 4 | + edge_capacity::AbstractMatrix, |
| 5 | + edge_cost::AbstractMatrix, |
| 6 | + optimizer; |
| 7 | + edge_demand::Union{Nothing,AbstractMatrix} = nothing, |
| 8 | + source_nodes = (), |
| 9 | + sink_nodes = () |
| 10 | + ) |
3 | 11 |
|
4 |
| -Find a flow satisfying the `node_demand` at each node and `edge_capacity` constraints for each edge |
5 |
| -while minimizing the `sum(edge_cost.*flow)`. |
| 12 | +Find a flow over a directed graph `g` satisfying the `node_demand` at each |
| 13 | +node and `edge_capacity` constraints for each edge while minimizing `dot(edge_cost, flow)`. |
6 | 14 |
|
7 |
| --`node_demand` is a vector of nodal values, which should be positive for sources and |
8 |
| - negative at sink nodes, any node with zero demand is a transport node |
| 15 | +Returns a sparse flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc. |
9 | 16 |
|
10 |
| -- The problem can be seen as a linear programming problem and uses a LP |
11 |
| -solver under the hood. We use Clp in the examples and tests. |
12 |
| -
|
13 |
| -Returns a flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc. |
14 | 17 | # Arguments
|
15 |
| -- `edge_demand::AbstractMatrix`: demand a minimum flow for a given edge. |
16 |
| -- `edge_demand_exact=true`: changes the capacity of a non zero edge to the demanded value |
17 |
| -- `source_nodes::AbstractVector` Sources at which the nodal netflow is allowed to be greater than nodal demand ** |
18 |
| -- `sink_nodes::AbstractVector` Sinks at which the nodal netflow is allowed to be less than nodal demand ** |
19 |
| - ** source_nodes & sink_nodes are only needed when nodal flow are not explictly set in node_demand |
| 18 | +
|
| 19 | +- `g` is a directed `LightGraphs.AbstractGraph`. |
| 20 | +-`node_demand` is a vector of nodal demand values, which should be negative for sources, |
| 21 | +positive for sink nodes, and zero for all other nodes. |
| 22 | +- `edge_capacity::AbstractMatrix` sets an upper bound on the flow of each arc. |
| 23 | +- `edge_cost::AbstractMatrix` the cost per unit of flow on each arc. |
| 24 | +- `optimizer` is an optimizer factory as defined in JuMP, passed at the construction of the `JuMP.Model`. |
| 25 | +
|
| 26 | +# Keyword arguments |
| 27 | +
|
| 28 | +- `edge_demand::Union{Nothing,AbstractMatrix}`: require a minimum flow for edges, or nothing. |
| 29 | +- `source_nodes` Collection of sources at which the nodal netflow is allowed to be greater than nodal demand, defaults to an empty tuple. |
| 30 | +- `sink_nodes` Collection of sinks at which the nodal netflow is allowed to be less than nodal demand, defaults to an empty tuple. |
| 31 | +
|
| 32 | +`source_nodes` & `sink_nodes` are only needed when nodal flow are not explictly set in node_demand |
20 | 33 |
|
21 | 34 | ### Usage Example:
|
22 | 35 |
|
23 | 36 | ```julia
|
24 |
| -julia> using Clp: ClpSolver # use your favorite LP solver here |
25 |
| -julia> g = lg.DiGraph(6) # Create a flow-graph |
26 |
| -julia> lg.add_edge!(g, 5, 1) |
27 |
| -julia> lg.add_edge!(g, 5, 2) |
28 |
| -julia> lg.add_edge!(g, 3, 6) |
29 |
| -julia> lg.add_edge!(g, 4, 6) |
30 |
| -julia> lg.add_edge!(g, 1, 3) |
31 |
| -julia> lg.add_edge!(g, 1, 4) |
32 |
| -julia> lg.add_edge!(g, 2, 3) |
33 |
| -julia> lg.add_edge!(g, 2, 4) |
34 |
| -julia> w = zeros(6,6) |
35 |
| -julia> w[1,3] = 10. |
36 |
| -julia> w[1,4] = 5. |
37 |
| -julia> w[2,3] = 2. |
38 |
| -julia> w[2,4] = 2. |
| 37 | +julia> import LightGraphs |
| 38 | +julia> const LG = LightGraphs |
| 39 | +julia> using LightGraphsFlows: mincost_flow |
| 40 | +julia> import Clp # use your favorite LP solver here |
| 41 | +julia> import JuMP |
| 42 | +julia> using SparseArrays: spzeros |
| 43 | +julia> g = LG.DiGraph(6) # Create a flow-graph |
| 44 | +julia> LG.add_edge!(g, 5, 1) |
| 45 | +julia> LG.add_edge!(g, 5, 2) |
| 46 | +julia> LG.add_edge!(g, 3, 6) |
| 47 | +julia> LG.add_edge!(g, 4, 6) |
| 48 | +julia> LG.add_edge!(g, 1, 3) |
| 49 | +julia> LG.add_edge!(g, 1, 4) |
| 50 | +julia> LG.add_edge!(g, 2, 3) |
| 51 | +julia> LG.add_edge!(g, 2, 4) |
| 52 | +julia> cost = zeros(6,6) |
| 53 | +julia> cost[1,3] = 10 |
| 54 | +julia> cost[1,4] = 5 |
| 55 | +julia> cost[2,3] = 2 |
| 56 | +julia> cost[2,4] = 2 |
39 | 57 | julia> demand = spzeros(6)
|
40 |
| -julia> demand[5] = 2 |
41 |
| -julia> demand[6] = -2 |
| 58 | +julia> demand[5] = -2 |
| 59 | +julia> demand[6] = 2 |
42 | 60 | julia> capacity = ones(6,6)
|
43 |
| -julia> flow = mincost_flow(g, demand, capacity, w, ClpSolver()) |
| 61 | +julia> flow = mincost_flow(g, demand, capacity, cost, JuMP.with_optimizer(Clp.Optimizer)) |
44 | 62 | ```
|
45 | 63 | """
|
| 64 | +function mincost_flow end |
46 | 65 |
|
47 |
| - |
48 |
| -function mincost_flow(g::lg.DiGraph, |
| 66 | +@traitfn function mincost_flow(g::AG::lg.IsDirected, |
49 | 67 | node_demand::AbstractVector,
|
50 | 68 | edge_capacity::AbstractMatrix,
|
51 | 69 | edge_cost::AbstractMatrix,
|
52 |
| - solver::AbstractMathProgSolver; |
| 70 | + optimizer; |
53 | 71 | edge_demand::Union{Nothing,AbstractMatrix} = nothing,
|
54 |
| - edge_demand_exact::Bool = false, #If true changes capacity of that node to the value in edge demand |
55 |
| - source_nodes::AbstractVector{<:Integer} = Vector{Integer}(), #Source nodes at which to allow a netflow greater than nodal demand |
56 |
| - sink_nodes::AbstractVector{<:Integer} = Vector{Integer}() #Sink nodes at which to allow a netflow less than nodal demand |
57 |
| - ) |
58 |
| - T = eltype(g) |
59 |
| - VT = promote_type(eltype(edge_capacity),eltype(node_demand),eltype(edge_cost)) |
60 |
| - |
61 |
| - nv = lg.nv(g) |
62 |
| - ne = lg.ne(g) |
63 |
| - b = Vector{VT}(undef,nv) |
64 |
| - AI = Vector{T}(undef,2*ne) |
65 |
| - AJ = Vector{T}(undef,2*ne) |
66 |
| - AV = Vector{VT}(undef,2*ne) |
67 |
| - edge_costs = Vector{VT}(undef,ne) |
68 |
| - upperbounds = Vector{VT}(undef,ne) |
69 |
| - lowerbounds = Vector{VT}(undef,ne) |
70 |
| - sense = ['=' for _ in 1:nv] |
71 |
| - |
72 |
| - has_edge_demand= edge_demand !== nothing |
73 |
| - |
74 |
| - inds = [-1,0] |
75 |
| - for (n,e) in enumerate(lg.edges(g)) |
76 |
| - inds .+= 2 |
77 |
| - s = lg.src(e) |
78 |
| - t = lg.dst(e) |
79 |
| - AI[inds] = [s,t] #the node index |
80 |
| - AJ[inds] = [n,n] #the edge index |
81 |
| - AV[inds] = [1,-1] |
82 |
| - upperbounds[n] = edge_capacity[s,t] |
83 |
| - lowerbounds[n] = has_edge_demand ? edge_demand[s,t] : zero(VT) |
84 |
| - if edge_demand_exact && has_edge_demand && lowerbounds[n] != 0 |
85 |
| - upperbounds[n] = lowerbounds[n] |
86 |
| - end |
87 |
| - |
88 |
| - edge_costs[n] = edge_cost[s,t] |
89 |
| - end |
90 |
| - for n = 1:nv |
91 |
| - b[n] = node_demand[n] |
92 |
| - end |
93 |
| - |
94 |
| - for n in source_nodes |
95 |
| - sense[n] = '>' |
| 72 | + source_nodes = (), # Source nodes at which to allow a netflow greater than nodal demand |
| 73 | + sink_nodes = () # Sink nodes at which to allow a netflow less than nodal demand |
| 74 | + ) where {AG <: lg.AbstractGraph} |
| 75 | + |
| 76 | + m = JuMP.Model(optimizer) |
| 77 | + vtxs = vertices(g) |
| 78 | + |
| 79 | + source_nodes = [v for v in vtxs if v in source_nodes || node_demand[v] < 0] |
| 80 | + sink_nodes = [v for v in vtxs if v in sink_nodes || node_demand[v] > 0] |
| 81 | + |
| 82 | + @variable(m, 0 <= f[i=vtxs,j=vtxs; (i,j) in lg.edges(g)] <= edge_capacity[i, j]) |
| 83 | + @objective(m, Min, sum(f[src(e),dst(e)] * edge_cost[src(e), dst(e)] for e in lg.edges(g))) |
| 84 | + |
| 85 | + for v in lg.vertices(g) |
| 86 | + if v in source_nodes |
| 87 | + @constraint(m, |
| 88 | + sum(f[v, vout] for vout in outneighbors(g, v)) - sum(f[vin, v] for vin in lg.inneighbors(g, v)) >= -node_demand[v] |
| 89 | + ) |
| 90 | + elseif v in sink_nodes |
| 91 | + @constraint(m, |
| 92 | + sum(f[vin, v] for vin in lg.inneighbors(g, v)) - sum(f[v, vout] for vout in outneighbors(g, v)) >= node_demand[v] |
| 93 | + ) |
| 94 | + else |
| 95 | + @constraint(m, |
| 96 | + sum(f[vin, v] for vin in lg.inneighbors(g, v)) == sum(f[v, vout] for vout in outneighbors(g, v)) |
| 97 | + ) |
| 98 | + end |
96 | 99 | end
|
97 |
| - for n in sink_nodes |
98 |
| - sense[n] = '<' |
99 |
| - end |
100 |
| - |
101 | 100 |
|
102 |
| - A = sparse(AI,AJ,AV,nv,ne) |
103 |
| - # return (edge_costs, A, sense, b, lowerbounds, upperbounds, solver) |
104 |
| - sol = linprog(edge_costs, A, sense, b, lowerbounds, upperbounds, solver) |
105 |
| - sol_sparse = sparse(view(AI,1:2:2*ne),view(AI,2:2:2*ne),sol.sol,nv,nv) |
106 |
| - return sol_sparse |
| 101 | + if edge_demand isa AbstractMatrix |
| 102 | + for e in lg.edges(g) |
| 103 | + (i,j) = Tuple(e) |
| 104 | + JuMP.set_lower_bound(f[i,j], edge_demand[i,j]) |
| 105 | + end |
| 106 | + end |
| 107 | + optimize!(m) |
| 108 | + ts = termination_status(m) |
| 109 | + result_flow = spzeros(nv(g), nv(g)) |
| 110 | + if ts != MOI.OPTIMAL |
| 111 | + @warn "Problem does not have an optimal solution, status: $(ts)" |
| 112 | + return result_flow |
| 113 | + end |
| 114 | + for e in lg.edges(g) |
| 115 | + (i,j) = Tuple(e) |
| 116 | + result_flow[i,j] = JuMP.value(f[i,j]) |
| 117 | + end |
| 118 | + return result_flow |
107 | 119 | end
|
108 |
| - |
109 |
| -@deprecate mincost_flow(g::lg.DiGraph, |
110 |
| - capacity::AbstractMatrix, |
111 |
| - demand::AbstractMatrix, |
112 |
| - cost::AbstractMatrix, |
113 |
| - solver::AbstractMathProgSolver, |
114 |
| - source::Int, # if source and/or sink omitted or not in nodes, circulation problem |
115 |
| - sink::Int) mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand,source_nodes=[source],sink_nodes=[sink]) |
116 |
| - |
117 |
| -@deprecate mincost_flow(g::lg.DiGraph, |
118 |
| - capacity::AbstractMatrix, |
119 |
| - demand::AbstractMatrix, |
120 |
| - cost::AbstractMatrix, |
121 |
| - solver::AbstractMathProgSolver |
122 |
| - ) mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand) |
0 commit comments