11.5 Dynamic programming in software
Software 11.5.1 (The knapsack table). The table \(f_k(c)\) for the knapsack problem, filled row by row; its last row; the items chosen, read back from the table; and a brute-force check over all \(2^4\) subsets:
>>> import numpy as np >>> from itertools import combinations, permutations >>> v = [10, 13, 7, 8]; w = [4, 6, 3, 5]; W = 10 # the knapsack of Section 9 >>> n = len(v) >>> f = np.zeros((n + 1, W + 1), int) # f[k, c]: first k items, room c >>> for k in range(1, n + 1): ... for c in range(W + 1): ... f[k, c] = f[k - 1, c] # skip item k ... ... if w[k - 1] <= c: # ... or take it ... f[k, c] = max(f[k, c], v[k - 1] + f[k - 1, c - w[k - 1]]) ... >>> f[n] # best value for every capacity 0..10 array([ 0, 0, 0, 7, 10, 10, 13, 17, 17, 20, 23]) >>> take, c = [], W # trace back the decisions >>> for k in range(n, 0, -1): ... if f[k, c] != f[k - 1, c]: take.append(k); c -= w[k - 1] ... >>> sorted(take), f[n, W] ([1, 2], 23) >>> subsets = [S for r in range(n + 1) for S in combinations(range(n), r)] >>> best = max((sum(v[i] for i in S), S) for S in subsets ... if sum(w[i] for i in S) <= W) >>> best # brute force over all 16 subsets (23, (0, 1))
Reading the output. The last row gives the best value for every capacity from \(0\) to \(10\); at \(10\) it is \(23\), taking items \(1\) and \(2\), which agrees with the enumeration in the section on integer programming and with the brute force (items \(0\) and \(1\) in Python’s numbering).
MATLAB/Octave: the same double loop over a matrix f.
Software 11.5.2 (Allocating salesmen). Example 11.3.1, with the value function stored in a dictionary indexed by stage and state, computed backwards, and the decisions read forwards; then the value functions of stages \(3\) and \(2\):
>>> R = np.array([[0, 20, 33, 42, 48, 52, 54], # sales (Rs. 1000), region 1
... [0, 15, 28, 39, 48, 55, 60], # region 2
... [0, 25, 38, 46, 51, 54, 56]]) # region 3; 0..6 salesmen
>>> N = 6
>>> F = {(4, s): 0 for s in range(N + 1)} # stage 4: no regions left, value 0
>>> choice = {}
>>> for k in (3, 2, 1): # backward: regions k..3, s salesmen
... for s in range(N + 1):
... val = lambda x: R[k - 1, x] + F[(k + 1, s - x)]
... choice[(k, s)] = max(range(s + 1), key=val)
... F[(k, s)] = val(choice[(k, s)])
...
>>> plan, s = [], N # forward: read off the decisions
>>> for k in (1, 2, 3):
... plan.append(choice[(k, s)]); s -= choice[(k, s)]
...
>>> plan, F[(1, N)]
([2, 2, 2], 99)
>>> [F[(3, s)] for s in range(N + 1)], [F[(2, s)] for s in range(N + 1)]
([0, 25, 38, 46, 51, 54, 56], [0, 25, 40, 53, 66, 77, 86])
Reading the output. Two salesmen to each region, total sales \(99\). \(f_2(4) = 66\) is the best that regions \(2\) and \(3\) can do with four salesmen (\(28 + 38\)), and the full table \(f_2\) answers the question for every number left over by region \(1\).
MATLAB/Octave: a matrix F(k, s+1) and the same loops.
Software 11.5.3 (Johnson’s rule). The makespan of a two-machine flow shop for a given order, Johnson’s order for Example 11.4.3, and the best and worst makespans over all \(720\) orders, with that of the order \(J_1, \ldots , J_6\):
>>> p = {'J1': (5, 2), 'J2': (1, 6), 'J3': (9, 7), 'J4': (3, 8),
... 'J5': (10, 4), 'J6': (7, 3)} # times on machine 1, then machine 2
>>> def makespan(order):
... t1 = t2 = 0
... for j in order:
... t1 += p[j][0]; t2 = max(t2, t1) + p[j][1]
... return t2
...
>>> first = sorted([j for j in p if p[j][0] < p[j][1]], key=lambda j: p[j][0])
>>> last = sorted([j for j in p if p[j][0] >= p[j][1]],
... key=lambda j: -p[j][1])
>>> johnson = first + last # Johnson's rule
>>> johnson, makespan(johnson)
(['J2', 'J4', 'J3', 'J5', 'J6', 'J1'], 37)
>>> spans = [makespan(o) for o in permutations(p)] # all 720 orders
>>> min(spans), max(spans), makespan(list(p))
(37, 52, 38)
Reading the output. Johnson’s order has makespan \(37\), the minimum over all orders; the worst order takes \(52\), and the natural order \(J_1, \ldots , J_6\) happens to take \(38\). The rule needed one sort; the enumeration, \(720\) evaluations, which would be \(2.4 \times 10^{18}\) for \(20\) jobs.
MATLAB/Octave: perms(1:6) for the enumeration.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.