Software 15.5.1 (Chains). Example 15.1.1: two-step probabilities, the stationary distribution as the left eigenvector of \(P\) for eigenvalue \(1\), and \(P^{20}\); then the fundamental matrix of Example 15.2.2, the expected times to absorption, and the absorption probabilities:

>>> import numpy as np
>>> P = np.array([[0.80, 0.10, 0.10],    # brand switching: row = this month,
...               [0.20, 0.70, 0.10],      # columns = next month (A, B, C)
...               [0.10, 0.30, 0.60]])
>>> np.linalg.matrix_power(P, 2)[0]        # two months after buying A
array([0.67, 0.18, 0.15])
>>> w, V = np.linalg.eig(P.T)              # pi P = pi: left eigenvector for 1
>>> pi = np.real(V[:, np.argmin(abs(w - 1))])
>>> pi = pi / pi.sum(); np.round(pi, 4)
array([0.45, 0.35, 0.2 ])
>>> np.round(np.linalg.matrix_power(P, 20), 4)   # every row tends to pi
array([[0.45, 0.35, 0.2 ],
       [0.45, 0.35, 0.2 ],
       [0.45, 0.35, 0.2 ]])
>>> Q = np.array([[0.0, 0.6], [0.0, 0.0]])  # accounts: 0-30 days, 31-60 days
>>> R = np.array([[0.4, 0.0], [0.7, 0.3]])  # -> paid, bad debt
>>> N = np.linalg.inv(np.eye(2) - Q)        # fundamental matrix
>>> N, N @ np.ones(2), N @ R                # visits, time to absorption, B
(array([[1. , 0.6],
       [0. , 1. ]]), array([1.6, 1. ]), array([[0.82, 0.18],
       [0.7 , 0.3 ]]))

Reading the output. The stationary distribution is \((0.45, 0.35, 0.2)\) and \(P^{20}\) has three identical rows equal to it. For the accounts, a new account is resolved in \(1.6\) months on average, and is paid with probability \(0.82\); an account already \(31\)–\(60\) days old is paid with probability \(0.7\).

MATLAB/Octave: [V, D] = eig(P’); N = inv(eye(2) - Q).

Software 15.5.2 (When to replace a machine). Example 15.3.2: the average annual cost for replacement every \(1, \ldots , 8\) years and its minimiser; then the same with a \(10\%\) discount rate, as the equivalent annual cost:

>>> price = 6000.0                   # machine: price, running costs, resale
>>> run = np.array([1000, 1200, 1400, 1800, 2300, 2800, 3400, 4000.])
>>> resale = np.array([3000, 1500, 750, 375, 200, 200, 200, 200.])
>>> years = np.arange(1, 9)
>>> avg = (price - resale + np.cumsum(run)) / years   # average cost a year
>>> np.round(avg, 1), int(years[np.argmin(avg)])
(array([4000. , 3350. , 2950. , 2756.2, 2700. , 2716.7, 2814.3, 2962.5]), 5)
>>> r = 0.10; d = 1 / (1 + r)        # with discounting: equivalent annual
>>> pv = np.array([price - resale[n - 1] * d**n + sum(run[k] * d**(k + 1)
...                for k in range(n)) for n in years])
>>> annuity = pv * r / (1 - d**years)      # cost per year with the same value
>>> np.round(annuity, 1), int(years[np.argmin(annuity)])
(array([4600. , 3838.1, 3373.4, 3131.4, 3030. , 3002.8, 3046.7, 3131.9]), 6)

Reading the output. Without discounting the average cost falls to \(2700\) at \(5\) years and then rises; with discounting the equivalent annual cost is lowest, \(3003\), at \(6\) years, though \(5\) years costs only \(27\) a year more. Both curves are flat near the optimum, as the EOQ cost was, so the decision is robust to small errors in the data.

MATLAB/Octave: (price - resale + cumsum(run))./(1:8).

Software 15.5.3 (Value and policy iteration). Example 15.4.2: value iteration to a tolerance of \(10^{-9}\), with the number of iterations, the optimal costs and the optimal action in each condition (\(1\) = replace); then policy iteration from the policy of never replacing:

>>> from scipy.optimize import linprog
>>> T = {0: np.array([[0.7, 0.2, 0.1, 0.0], [0.0, 0.6, 0.3, 0.1],
...                   [0.0, 0.0, 0.5, 0.5], [0.0, 0.0, 0.0, 1.0]]),
...      1: np.tile([0.7, 0.2, 0.1, 0.0], (4, 1))}    # replace: as if new
>>> cost = {0: np.array([0., 1000, 3000, 6000]),   # running cost by condition
...         1: np.full(4, 4000.)}                     # replacement cost
>>> beta = 0.9
>>> v = np.zeros(4)
>>> for it in range(500):                  # value iteration
...     q = np.array([cost[a] + beta * T[a] @ v for a in (0, 1)])
...     v_new = q.min(axis=0)
...     if np.abs(v_new - v).max() < 1e-9: break
...     v = v_new
...
>>> it, np.round(v, 1), q.argmin(axis=0)   # 1 = replace in that condition
(264, array([ 9450., 12700., 13450., 13450.]), array([0, 0, 1, 1]))
>>> pol = np.zeros(4, int)           # policy iteration from "always keep"
>>> for k in range(10):
...     Tp = np.array([T[pol[s]][s] for s in range(4)])
...     cp = np.array([cost[pol[s]][s] for s in range(4)])
...     vp = np.linalg.solve(np.eye(4) - beta * Tp, cp)
...     qp = np.array([cost[a] + beta * T[a] @ vp for a in (0, 1)])
...     new = qp.argmin(axis=0)
...     if (new == pol).all(): break
...     pol = new
...
>>> k, pol, np.round(vp, 1)
(2, array([0, 0, 1, 1]), array([ 9450., 12700., 13450., 13450.]))

Reading the output. Value iteration takes \(264\) iterations, as the contraction factor \(0.9\) predicts for this tolerance; policy iteration needs \(2\) improvements. Both give costs \((9450, 12\,700, 13\,450, 13\,450)\) and the policy keep, keep, replace, replace. The cost in conditions \(2\) and \(3\) is the same because the machine is replaced at once in both.

MATLAB/Octave: the same loops; the MDP toolbox has mdp_value_iteration.