19  Serial Systems

19.1 Local model

• I'_j: Stage j’s local on-hand inventory
• S'_j: Stage j’s local base-stock level
• IL'_j: Stage j’s local inventory level
• B'_j: Stage j’s local backorders

IL'_j = I'_j - B'_j

The expected cost function for a serial system can be expressed in terms of echelon base-stock levels as follows:

g'(\mathbf{S'}) = \mathbb{E} \left[ \sum_{j=1}^N h'_j (I'_j + IT_{j-1}) + p B'_1 \right]

19.2 Echelon model

• I_j: Stage j’s echelon on-hand inventory
• IT_j: Inventory in transit from stage j+1 to stage j
• S_j: Stage j’s echelon base-stock level

Stage j’s echelon on-hand inventory is the sum of the local on-hand inventory at stage j and the inventory in transit among all stages downstream of stage j:

I_j = \sum_{i=1}^j (I'_i + IT_{i-1})

where IT_{i-1} is the inventory in transit from stage i to stage i-1 and IT_0 = 0.

Stage j’s echelon inventory level is given by

IL_j = I_j - B'_1

where B'_1 is the backorders at the retailer (stage 1).

Stage j’s echelon holding cost is given by

h_j = h'_j - h'_{j+1}.

Here, we define h'_{N+1} = 0.

Also, we know that

h'_j = \sum_{k=j}^N h_k.

命題 19.1 The total holding cost across all stages is given by

\sum_{j=1}^N h_j I_j = \sum_{j=1}^N h'_j (I'_j + IT_{j-1})

証明. We know that h_j = h'_j - h'_{j+1}, h'_{N+1} = 0, so we have

\begin{aligned} \sum_{j=1}^N h_j I_j &= \sum_{j=1}^N (h'_j - h'_{j+1}) I_j \\ &= \sum_{j=1}^N h'_j I_j - \sum_{j=1}^N h'_{j+1} I_j \\ &= \sum_{j=1}^N h'_j I_j - \sum_{j=2}^{N+1} h'_j I_{j-1} \\ &= \sum_{j=1}^N h'_j I_j - \sum_{j=2}^{N} h'_j I_{j-1} \\ &= h'_1 I_1 + \sum_{j=2}^N h'_j I_j - \sum_{j=2}^{N} h'_j I_{j-1} \\ &= h'_1 I_1 + \sum_{j=2}^N h'_j (I_j - I_{j-1}) \end{aligned}

For j = 1, I_1 = I'_1 + IT_0. For j \geq 2, we have I_j = I'_j + IT_{j-1} + I_{j-1}. Therefore, we have I_j - I_{j-1} = I'_j + IT_{j-1} for j \geq 2. Thus, we have

\begin{aligned} \sum_{j=1}^N h_j I_j &= h'_1 (I'_1 + IT_0) + \sum_{j=2}^N h'_j (I'_j + IT_{j-1}) \\ &= \sum_{j=1}^N h'_j (I'_j + IT_{j-1}) \end{aligned}

Local base-stock levels can be derived from echelon base-stock levels as follows:

S'_j = S_j - S_{j-1}

Echelon base-stock levels can be expressed as the sum of local base-stock levels:

S_j = \sum_{i=1}^j S'_i

Here, we define S_0 = 0.

19.3 Cost Function

The expected cost function can also be expressed in terms of local base-stock levels as follows:

g(\mathbf{S}) = \mathbb{E} \left[ \sum_{j=1}^N h_j IL_j + (p + h'_1) IL^-_1 \right]

証明. We have g'(\mathbf{S'}) = \mathbb{E} \left[ \sum_{j=1}^N h'_j (I'_j + IT_{j-1}) + p B'_1 \right]. We want to show that g'(\mathbf{S'}) = g(\mathbf{S}).

We have g(\mathbf{S}) = \mathbb{E} \left[ \sum_{j=1}^N h_j IL_j + (p + h'_1) IL^-_1 \right]. We can express h_j and IL_j in terms of h'_j, I'_j, and B'_1 as follows:

\begin{aligned} h_j &= h'_j - h'_{j+1} \\ IL_j &= I_j - B'_1 = \sum_{i=1}^j (I'_i + IT_{i-1}) - B'_1 \end{aligned}

Substituting these expressions into g(\mathbf{S}), we have

\begin{aligned} g(\mathbf{S}) &= \mathbb{E} \left[ \sum_{j=1}^N (h'_j - h'_{j+1}) \left( \sum_{i=1}^j (I'_i + IT_{i-1}) - B'_1 \right) + (p + h'_1) IL^-_1 \right] \\ &= \mathbb{E} \left[ \sum_{j=1}^N (h'_j - h'_{j+1}) \sum_{i=1}^j (I'_i + IT_{i-1}) - \sum_{j=1}^N (h'_j - h'_{j+1}) B'_1 + (p + h'_1) IL^-_1 \right] \\ &= \mathbb{E} \left[ \sum_{j=1}^N (h'_j - h'_{j+1}) \sum_{i=1}^j (I'_i + IT_{i-1}) + p B'_1 + h'_1 IL^-_1 - \sum_{j=1}^N (h'_j - h'_{j+1}) B'_1 \right] \\ &= \mathbb{E} \left[ \sum_{j=1}^N (h'_j - h'_{j+1}) \sum_{i=1}^j (I'_i + IT_{i-1}) + p B'_1 \right] \\ &= g'(\mathbf{S'}) \end{aligned}

Here, we have used the fact that \sum_{j=1}^N (h'_j - h'_{j+1}) B'_1 = h'_1 B'_1 and IL^-_1 = \max\{-IL_1, 0\} = \max\{B'_1 - I_1, 0\} = B'_1.

Let ITP_j be the inventory position at stage j, which includes the echelon inventory level at stage j and the inventory in transit from stage j+1 to stage j. We have

ITP_j = IL_j + IT_j = IP_j - (IL'_{j+1})^-

D_j denotes the lead time demand at stage j. After L_j periods, all the inventory in transit from stage j+1 to stage j will arrive at stage j. Therefore, we have

IL_j(t + L_j) = ITP_j(t) - D_j

If IL_{j+1}(t) \geq S_j, then ITP_j(t) = S_j. If IL_{j+1}(t) < S_j, then ITP_j(t) = IL_{j+1}(t). Therefore, we have

ITP_j(t) = \min\{S_j, IL_{j+1}(t)\}

証明. omitted

At stage N, we have

IP_N(t) = ITP_N(t) = \min\{S_N, IL_{N+1}(t)\} = S_N,

since stage N does not stock out.

In steady state,

\begin{aligned} ITP_N &= S_N \\ IL_j &= ITP_j - D_j \\ ITP_j &= \min\{S_j, IL_{j+1}\} \end{aligned}

We introduce the following auxiliay functions

\begin{aligned} \hat{g}_j(x \mid \mathbf{S}) &= \mathbb{E} \left[ \sum_{i=1}^j h_i IL_i + (p + h'_1) IL^-_1 \mid IL_j = x \right] \\ g_{j}(y \mid \mathbf{S}) &= \mathbb{E} \left[ \sum_{i=1}^j h_i IL_i + (p + h'_1) IL^-_1 \mid ITP_j = y \right] \\ \underline{g}_j(x \mid \mathbf{S}) &= \mathbb{E} \left[ \sum_{i=1}^j h_i IL_i + (p + h'_1) IL^-_1 \mid IL_{j + 1} = x \right] \end{aligned}

\begin{aligned} \underline{g}_0 (x \mid \mathbf{S}) &= \mathbb{E} \left[ (p + h'_1) IL^-_1 \mid IL_1 = x \right] \\ &= (p + h'_1) x^- \end{aligned}

\begin{aligned} \hat{g}_1 (x \mid \mathbf{S}) &= \mathbb{E} \left[ h_1 IL_1 + (p + h'_1) IL^-_1 \mid IL_1 = x \right] \\ &= h_1 x + \underline{g}_0 (x \mid \mathbf{S}) \\ \end{aligned}

\begin{aligned} g_1 (y \mid \mathbf{S}) &= \mathbb{E} \left[ h_1 IL_1 + (p + h'_1) IL^-_1 \mid ITP_1 = y \right] \\ &= \mathbb{E}_{D_1} \left[ \mathbb{E} \left[ h_1 IL_1 + (p + h'_1) IL^-_1 \mid IL_1 = y - D_1 \right] \right] \\ &= \mathbb{E}_{D_1} \left[ \hat{g}_1 (y - D_1 \mid \mathbf{S}) \right] \\ \end{aligned}

\begin{aligned} \underline{g}_1 (x \mid \mathbf{S}) &= \mathbb{E} \left[ h_1 IL_1 + (p + h'_1) IL^-_1 \mid IL_2 = x \right] \\ &= \mathbb{E} \left[ h_1 IL_1 + (p + h'_1) IL^-_1 \mid ITP_1 = \min\{S_1, x\} \right] \\ &= g_1 (\min\{S_1, x\} \mid \mathbf{S}) \end{aligned}

\begin{aligned} \hat{g}_2 (x \mid \mathbf{S}) &= \mathbb{E} \left[ h_1 IL_1 + h_2 IL_2 + (p + h'_1) IL^-_1 \mid IL_2 = x \right] \\ &= h_2 x + \underline{g}_1 (x \mid \mathbf{S}) \\ \end{aligned}

In general, for j = 1, \ldots, N, given g_{j-1}(\cdot \mid \mathbf{S}), we have

\begin{aligned} \hat{g}_j (x \mid \mathbf{S}) &= h_j x + \underline{g}_{j-1} (x \mid \mathbf{S}) \\ g_j (y \mid \mathbf{S}) &= \mathbb{E} \left[ \hat{g}_j (y - D_j \mid \mathbf{S}) \right] \\ \underline{g}_j (x \mid \mathbf{S}) &= g_j (\min\{S_j, x\} \mid \mathbf{S}) \end{aligned}

The expected cost of the system can be expressed as

g(\mathbf{S}) = g_N (S_N \mid \mathbf{S}) = \mathbb{E} \left[ \sum_{j=1}^N h_j IL_j + (p + h'_1) IL^-_1 \mid ITP_N = S_N \right]

定理 19.1 Let \underline{g}_0 (x) = (p + h'_1) x^-, and for j = 1, \ldots, N, let

\begin{aligned} \hat{g}_j (x) &= h_j x + \underline{g}_{j-1} (x) \\ g_j (y) &= \mathbb{E} \left[ \hat{g}_j (y - D_j) \right] \\ S^*_j &= \arg\min g_j (y) \\ \underline{g}_j (x) &= g_j (\min\{S^*_j, x\}) \end{aligned}

Then, the optimal echelon base-stock levels are given by \mathbf{S}^* = (S^*_1, \ldots, S^*_N) and the optimal cost is given by g(\mathbf{S}^*).

19.4 Heuristic Approach

g_j(y) = \mathbb{E} \left[ \sum_{i=1}^j (\sum_{k=i}^j h_k) (I'_i + IT_{i-1}) + (p + h'_{j+1}) B'_1(y) \right]

The echelon holding cost for stage 1,\ldots, j can be represented as a vector \mathbf{h}^j = (h_1, \ldots, h_j). We consider two different ways to set h^u and h^l. h^u = \sum_{k=i}^j h_k and h^l = h_j.

The approximate \tilde{S_j} can be obtained by

\tilde{S_j} = \frac{1}{2} \left( \tilde{F}_j^{-1} \left( \frac{p + \sum_{i=j+1}^N h_i}{p + \sum_{i=j}^N h_i} \right) + \tilde{F}_j^{-1} \left( \frac{p + \sum_{i=j+1}^N h_i}{p + \sum_{i=1}^N h_i} \right) \right)

19.5 Python Implementation

例 19.1 D_1 \sim N(5, 1^2), L_1 = L_2 = 1, L_3 = 2, (h'_1, h'_2, h'_3) = (7, 4, 2), (h_1, h_2, h_3) = (3, 2, 2), p = 37.12.

Solve: \begin{aligned} \tilde{D}_1 &\sim N(5, 1^2) = N(5, 1) \\ \tilde{D}_2 &\sim N(5 \cdot 2, 1^2 \cdot 2) = N(10, 2) \\ \tilde{D}_3 &\sim N(5 \cdot 4, 1^2 \cdot 4) = N(20, 4) \end{aligned}

\begin{aligned} \tilde{S}_1 &= \frac{1}{2} \left( \tilde{F}_1^{-1} \left( \frac{37.12 + 4}{37.12 + 7} \right) + \tilde{F}_1^{-1} \left( \frac{37.12 + 4}{37.12 + 7} \right) \right) = 6.49 \\ \tilde{S}_2 &= \frac{1}{2} \left( \tilde{F}_2^{-1} \left( \frac{37.12 + 2}{37.12 + 4} \right) + \tilde{F}_2^{-1} \left( \frac{37.12 + 2}{37.12 + 7} \right) \right) = 12.03 \\ \tilde{S}_3 &= \frac{1}{2} \left( \tilde{F}_3^{-1} \left( \frac{37.12 + 0}{37.12 + 2} \right) + \tilde{F}_3^{-1} \left( \frac{37.12 + 0}{37.12 + 7} \right) \right) = 22.63 \end{aligned}

from scipy.stats import norm

def approximate_base_stock_levels(demand_means, demand_stds, echelon_holding_costs, backorder_cost):
    if not (len(demand_means) == len(demand_stds) == len(echelon_holding_costs)):
        raise ValueError("All input lists must have the same length.")

    total_holding_cost = sum(echelon_holding_costs)
    base_stock_levels = []

    for j, (mean_j, std_j) in enumerate(zip(demand_means, demand_stds)):
        stockout_cost = backorder_cost + sum(echelon_holding_costs[j + 1 :])
        holding_cost_upper = sum(echelon_holding_costs[j:])

        service_level_1 = stockout_cost / (backorder_cost + holding_cost_upper)
        service_level_2 = stockout_cost / (backorder_cost + total_holding_cost)

        s1 = norm.ppf(service_level_1, loc=mean_j, scale=std_j)
        s2 = norm.ppf(service_level_2, loc=mean_j, scale=std_j)
        base_stock_levels.append(0.5 * (s1 + s2))

    return base_stock_levels


# Parameters
demand_means = [5, 10, 20]
demand_stds = [1, np.sqrt(2), np.sqrt(4)]
echelon_holding_costs = [3, 2, 2]
backorder_cost = 37.12

base_stock_levels = approximate_base_stock_levels(
    demand_means,
    demand_stds,
    echelon_holding_costs,
    backorder_cost,
)

print("Approximate base-stock levels:", [round(s, 2) for s in base_stock_levels])

# case 2
# customer demand N(5, 1^2), L1 = 2, L2 = 2, L3 = 3, (h'_1, h'_2, h'_3) = (7, 4, 2), (h_1, h_2, h_3) = (3, 2, 2), p = 37.12
demand_means = [5 * 2, 5 * (2 + 2), 5 * (2 + 2 + 3)]
demand_stds = [1 * np.sqrt(2), 1 * np.sqrt(2 + 2), 1 * np.sqrt(2 + 2 + 3)]
echelon_holding_costs = [3, 2, 2]
backorder_cost = 37.12
base_stock_levels = approximate_base_stock_levels(
    demand_means,
    demand_stds,
    echelon_holding_costs,
    backorder_cost,
)   
print("Approximate base-stock levels for case 2:", [round(s, 2) for s in base_stock_levels])