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.