Network Flow Problems¶
The network flow problem is defined as
Rewriting into standard form gives
Affixing (dual) slack variables to the foregoing formulation yields
which differs from the textbook’s notation by a negative sign.
Given a network flow problem, any selection of primal flow values that satisfies the balance equations at every node will be called a balanced flow. If all the flows are nonnegative, then a balanced flow is called a feasible flow.
Given a spanning tree, a balanced flow that assigns zero flow to every arc not on the spanning tree will be called a tree solution. From duality theory, the tree solution is optimal if all the primal flows are nonnegative and if all the dual slacks are nonnegative.
[3]:
from operator import itemgetter
from pprint import pprint
from pymprog import model
import itertools
import numpy
class NetworkFlowSolver:
def __init__(self, nodes, edges):
_ = sorted(edges, key=itemgetter(0, 2))
self.arcs = [edge[0] + edge[2] for edge in _]
self.c = [edge[1] for edge in _]
self.n = len(self.arcs)
self.m = len(nodes)
self.nodes = nodes
self.edges = edges
self.A = numpy.zeros((self.m, self.n))
for i, (name, _) in enumerate(nodes):
for j, arc in enumerate(self.arcs):
if name == arc[0]:
self.A[i, j] = -1
elif name == arc[1]:
self.A[i, j] = 1
self.b = numpy.asarray([_[1] for _ in nodes])
def solve_primal(self):
lp = model('Original Primal LP')
#lp.verbose(True)
x = lp.var('x', self.n, bounds=(0, None))
lp.minimize(sum(x_j * c_j for x_j, c_j in zip(x, self.c)))
for i, b_i in enumerate(self.b):
sum(a_ij * x_j for a_ij, x_j in zip(self.A[i, :], x)) == -b_i
lp.solve()
lp.sensitivity()
lp.end()
self.lp = lp
self.x = x
@property
def primal_flows(self):
_ = []
for x_i, arc_name in zip(self.x, self.arcs):
if not numpy.isclose(x_i.primal, 0):
_.append((arc_name, x_i.primal))
groups = {}
for k, v in itertools.groupby(_, key=lambda x: x[0][0]):
groups[k] = list(v)
return groups
@property
def dual_variables(self):
_ = []
for con in range(1, self.lp.get_num_rows() + 1):
_.append((self.nodes[con - 1][0], self.lp.get_row_dual(con)))
return _
@property
def dual_slacks(self):
_ = []
for x_i, arc_name in zip(self.x, self.arcs):
if numpy.isclose(x_i.primal, 0):
_.append((arc_name, x_i.dual))
groups = {}
for k, v in itertools.groupby(_, key=lambda x: x[0][0]):
groups[k] = list(v)
return groups
Exercise 14.1¶
The proposed tree solution is not optimal. It is neither primal feasible nor dual feasible.
[4]:
# supplies/demands
nodes = [('a', -9),
('b', 0),
('c', 3),
('d', 17),
('e', -8),
('f', 0),
('g', -8),
('h', -3),
('i', 8),]
# costs
edges = [('a', -3, 'b'), ('a', 0, 'f'),
('b', 1, 'g'), ('b', 5, 'h'),
('c', 10, 'b'), ('c', 12, 'a'),
('d', 1, 'h'), ('d', 3, 'e'), ('d', 1, 'a'),
('e', 6, 'a'),
('f', 5, 'c'), ('f', 11, 'e'),
('g', 3, 'h'), ('g', 5, 'i'),
('h', 6, 'a'),
('i', 7, 'b'),
('i', 3, 'c')]
nfs = NetworkFlowSolver(nodes, edges)
nfs.solve_primal()
PyMathProg 1.0 Sensitivity Report Created: 2018/03/04 Sun 15:49PM
================================================================================
Variable Activity Dual.Value Obj.Coef Range.From Range.Till
--------------------------------------------------------------------------------
*x[0] 0 0 -3 -4 -2
x[1] 0 9 0 -9 inf
*x[2] 8 0 1 0 1.79769e+308
x[3] 0 2 5 3 inf
*x[4] 3 0 12 7 13
x[5] 0 1 10 9 inf
*x[6] 6 0 1 0 7
*x[7] 8 0 3 -5 12
*x[8] 3 0 1 -5 2
x[9] 0 8 6 -2 inf
x[10] 0 8 5 -3 inf
*x[11] 0 0 11 2 19
x[12] 0 1 3 2 inf
x[13] 0 13 5 -8 inf
x[14] 0 6 6 0 inf
*x[15] 8 0 7 -6 12
x[16] 0 5 3 -2 inf
================================================================================
================================================================================
Constraint Activity Dual.Value Lower.Bnd Upper.Bnd RangeLower RangeUpper
--------------------------------------------------------------------------------
R1 9 10 9 9 9 9
R2 0 7 0 0 0 0
R3 -3 -2 -3 -3 -3 -3
R4 -17 9 -17 -17 -17 -17
R5 8 12 8 8 8 8
R6 0 1 0 0 0 0
R7 8 8 8 8 8 8
R8 3 10 3 3 3 3
*R9 -8 0 -8 -8 -8 -8
================================================================================
(a)¶
Recall that each node \(k \in \mathcal{N}\) must satisfy the flow balance constraint
The primal flows for the proposed spanning tree can be obtained by starting at the leaves and working inward:
The remaining non-tree arcs have zero flow.
The optimal primal flows for each tree arc correspond to the LP’s basic variables.
[5]:
pprint(nfs.primal_flows)
{'b': [('bg', 8.0)],
'c': [('ca', 3.0)],
'd': [('da', 6.0), ('de', 8.0), ('dh', 3.0)],
'i': [('ib', 8.0)]}
(b)¶
The dual feasibility conditions are
By complementarity, \(z_{ij} = 0\) for each \((i, j)\) in the spanning tree \(\mathcal{T}\).
Let node \(d\) be the root node in the proposed spanning tree. The dual variables can be obtained by starting at the root and working outward:
The optimal dual variables for each node correspond to the duals of the primal slacks.
[6]:
pprint(nfs.dual_variables)
[('a', 10.0),
('b', 7.0),
('c', -2.0),
('d', 9.0),
('e', 12.0),
('f', 1.0),
('g', 8.0),
('h', 10.0),
('i', 0.0)]
(c)¶
Given the dual variables, the dual slacks must satisfy
The dual slacks for the proposed spanning tree are
The optimal dual slacks for each non-tree arc correspond to the duals of the nonbasic variables.
[7]:
pprint(nfs.dual_slacks)
{'a': [('ab', 0.0), ('af', 9.0)],
'b': [('bh', 2.0)],
'c': [('cb', 1.0)],
'e': [('ea', 8.0)],
'f': [('fc', 8.0), ('fe', 0.0)],
'g': [('gh', 1.0), ('gi', 13.0)],
'h': [('ha', 6.0)],
'i': [('ic', 5.0)]}
Exercise 14.2¶
(a)¶
The leaving arc with the most infeasible primal flow is \((c, a)\).
(b)¶
Removing the leaving arc disconnects the spanning tree into two disjoint subtrees \(\mathcal{T}_c = \{b, c, d, g, h, i\}\) and \(\mathcal{T}_a = \{a, e, f\}\).
The bridging arcs are \(\left\{ (f, c), (a, b), (h, a), (d, a), (d, e) \right\}\).
The arcs in the opposite direction from the leaving arc \((c, a)\) are \(\left\{ (f, c), (a, b) \right\}\). Here opposite direction means the arc must start from \(\mathcal{T}_a\) and enter \(\mathcal{T}_c\).
The entering arc with the smallest dual slack is \((a, b)\).
(c)¶
The new tree solution after the pivot is given by steps (9), (10), and (11) of the self-dual network simplex method.
Update Primal Flows¶
The entering arc creates a cycle containing nodes \(c, a, b\). Let \(t\) denote the flow of the leaving arc \((c, a)\).
\((a, b)\) is oriented in the same direction as the leaving arc \((c, a)\), so its flow is decreased by \(t\). \((c, b)\) is in the opposite direction, so its flow increased by \(t\).
Update Dual Slacks¶
All dual slacks remain unchanged except for those associated with nontree arcs that bridge the two subtrees \(\mathcal{T}_a\) and \(\mathcal{T}_c\).
The old dual slack of the entering arc is \(z_{ab} = 8\). The dual slacks corresponding to those arcs that bridge in the same direction as the entering arc get decremented by \(z_{ab}\) whereas those that correspond to arcs bridging in the opposite direction get incremented by this amount.
Update Dual Variables¶
The dual variables for nodes in \(\mathcal{T}_a\) remains unchanged. The dual variables for nodes in \(\mathcal{T}_c\) get incremented by the old dual slack of the entering arc \(z_{ab}\).
Exercise 14.3¶
The proposed tree solution is not optimal. It is primal feasible, but not dual feasible.
(a)¶
The primal flows for the proposed spanning tree are
The remaining non-tree arcs have zero flow.
(b)¶
Let node \(c\) be the root node in the proposed spanning tree. The dual variables are
(c)¶
The dual slacks are
Exercise 14.4¶
Exercise 14.5¶
Exercise 14.6¶
[9]:
nodes = [('a', 2),
('b', 1),
('c', -2),
('d', 6),
('e', -1),
('f', 0),
('g', -3),
('h', -3),]
edges = [('a', 3, 'b'), ('a', 1, 'd'),
('b', 0, 'c'),
('c', 0, 'a'), ('c', 1, 'f'),
('d', 2, 'e'), ('d', 4, 'g'),
('e', 2, 'b'), ('e', 3, 'h'),
('f', 1, 'g'),
('g', 0, 'h'),
('h', 1, 'f')]
nfs = NetworkFlowSolver(nodes, edges)
nfs.solve_primal()
PyMathProg 1.0 Sensitivity Report Created: 2018/03/04 Sun 15:56PM
================================================================================
Variable Activity Dual.Value Obj.Coef Range.From Range.Till
--------------------------------------------------------------------------------
*x[0] 1 0 3 3 5
*x[1] 1 0 1 -1 1
*x[2] 2 0 0 0 2
x[3] 0 3 0 -3 inf
*x[4] 0 0 1 1 3
*x[5] 1 0 2 1 1.79769e+308
*x[6] 6 0 4 2 4
x[7] 0 2 2 0 inf
x[8] 0 1 3 2 inf
x[9] 0 0 1 1 inf
*x[10] 3 0 0 -2 1
x[11] 0 2 1 -1 inf
================================================================================
================================================================================
Constraint Activity Dual.Value Lower.Bnd Upper.Bnd RangeLower RangeUpper
--------------------------------------------------------------------------------
R1 -2 -5 -2 -2 -2 -2
R2 -1 -2 -1 -1 -1 -1
R3 2 -2 2 2 2 2
R4 -6 -4 -6 -6 -6 -6
R5 1 -2 1 1 1 1
R6 0 -1 0 0 0 0
R7 3 0 3 3 3 3
*R8 3 0 3 3 3 3
================================================================================
[11]:
pprint(nfs.primal_flows)
{'a': [('ab', 1.0), ('ad', 1.0)],
'b': [('bc', 2.0)],
'd': [('de', 1.0), ('dg', 6.0)],
'g': [('gh', 3.0)]}
[12]:
pprint(nfs.dual_variables)
[('a', -5.0),
('b', -2.0),
('c', -2.0),
('d', -4.0),
('e', -2.0),
('f', -1.0),
('g', 0.0),
('h', 0.0)]
[10]:
pprint(nfs.dual_slacks)
{'c': [('ca', 3.0), ('cf', 0.0)],
'e': [('eb', 2.0), ('eh', 1.0)],
'f': [('fg', 0.0)],
'h': [('hf', 2.0)]}
[13]:
#http://orfe.princeton.edu/~rvdb/307/lectures/lec14_show.pdf
#http://orfe.princeton.edu/~rvdb/307/lectures/lec13_show.pdf