# Research into Inventory Optimization of a Two-Echelon Distribution System Based on Particle Swarm Optimization (PSO)

Jingti Han, Mingliang Yu, Lin Liu

► To cite this version:

Jingti Han, Mingliang Yu, Lin Liu. Research into Inventory Optimization of a Two-Echelon Distribution System Based on Particle Swarm Optimization (PSO). 15th International Conference on Information and Semiotics in Organisations (ICISO), May 2014, Shanghai, China. pp.430-438, 10.1007/978-3-642-55355-4_45. hal-01350953

HAL Id: hal-01350953

https://inria.hal.science/hal-01350953

Submitted on 2 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research into Inventory Optimization of a Two-echelon
Distribution System based on Particle Swarm
Optimization (PSO)

HAN Jingti¹, YU Mingliang¹,², LIU Lin¹

¹School of Information Management and Engineering, Shanghai University of Finance and Economics, Shanghai 200433, China

²Labs of Economics and Management, University of Shanghai for Science and Technology, Shanghai 200093, China
{YU Mingliang, ymllw}@163.com

**Abstract.** Based on the assumption that distributor and retailers adopt different replenishment policies, this article gives consideration to the relationship between the inventory strategies of the two echelon. The authors construct an inventory model of a two-echelon distribution system under stochastic demand, and solve the problem based on algorithm of Particle Swarm Optimization (PSO) to quantitatively get the optimal ordering strategies of distribution center and retailers.

**Keywords:** Different replenishment policies, Two-echelon distribution system, Particle Swarm Optimization (PSO), Inventory model

# 1 Introduction

Inventory optimization of two-echelon distribution system is exactly the typical problem in multi-echelon inventory control of supply chain, as for the demand uncertainty of actual management, many experts make research on two-echelon inventory system based on random demand. [1] established a two-echelon inventory system jointly operation model including a single supplier (or distributor) and many retailers. [2] proposed an inventory control model of a two-echelon distribution system based on decentralized control strategy. [3] put forward that queuing theory can be used to simulate the process retailers waiting for the delivery under the assumption of Poisson distribution, thus a two-echelon inventory model based on (R,Q) inventory replenishment policy under stochastic demand was established. [4] put great emphasis on quantitative impact of retailers inventory strategy on demands faced by suppliers.

In this paper, we study a typical two-echelon distribution system, the distribution system consists a distributor and many retailers. Based on the assumption that distributor and retailers adopt different replenishment policies, we assume retailers with (R,Q) policy and distributor with (s,S) policy, the order quantity Q remains unchanged for each order. This article gives consideration to the relationship between the inventory strategies of the two echelon, aim to minimum the average total cost, establish an inventory model of a two-echelon distribution system and solve the
problem based on algorithm of Particle Swarm Optimization (PSO) with a numerical example.

# 2 Assumptions and Symbols

## 2.1 Assumptions

a. In the two-echelon distribution system, distributor sends ordering information to suppliers and provides single species products. b. The final demands of customers faced by various retailers are mutually independent random variables, and demands per unit time follow a Poisson distribution with the same parameters. c. The retailers adopt a continuous review (R,Q) strategy. Assuming $Q \ge R+1$, so that order crossover phenomenon will not occur within two successive cycles. d. Distributor adopt a continuous stock-taking (s,S) strategy and $s < S$, Where s and S are integer multiple of order quantity Q. e. Backorders are allowed in each echelon of the two-echelon system and are delay in delivery, backorders are sale without loss, namely Backorder mode, delay in delivery has a corresponding cost f. Lead time a random variable, the lead time is independent from unit time demand of retailers and it is composed of transport time and random delay of distributor. g. Assuming transportation time and cost are the same from distributor to different retailers.

## 2.2 Symbols

$N$ - The total number of retailers,

$L_0$ - Fixed lead time of distributor's purchases from suppliers,

$A_0$ - Fixed ordering cost per distributor order

$h_0$ - Holding cost of per unit goods within per unit time of distributor,

$b_0$ - Backorder cost of per unit goods within per unit time of distributor,

$L_1$ - Fixed shipping time from distributor to different suppliers,

$x_r$ - Retailers' demand per unit time, $x_r \sim P(\lambda_r)$,

$h_r$ - Holding cost of per unit goods within per unit time of retailers,

$A_r$ - Fixed ordering cost per retailer order,

$b_r$ - Backorder cost of per unit goods within per unit time of retailers,

$\lambda_0$ - The average arrival rate of distributor's demand,

$x_0$ - The number of orders reaches distributor per unit time,

$D(t, t+\tau) = D(\tau)$ - Random demand of distributor from time period $t$ to $t+\tau$,

$\delta_j(t)$ - The probability of the demand arrived at distribution center at moment $t$,

$\omega_j(\tau)$ - The probability that the number of orders arriving at distributor is $j$ in time period $\tau$,

$L_d$ - Random delay of retailers' orders at distribution center,
$L_r$ -The lead time of retailers' orders.

# 3 An inventory model of a two-echelon distribution system

## 3.1 An inventory model of distribution center

### A. Demand analysis on distribution center

The average arrival rate of distributor's demand is

$$ \lambda_0 = \sum_{i=1}^{N} \frac{\lambda_i}{Q} = \frac{N\lambda_r}{Q_r} \quad (3.1) $$

The orders received by distribution center $x_0$ obey Poisson distribution with parameter $\lambda_0$, that is $x_0 \sim P(\lambda_0)$. As the batch of each order is $Q_r$, the probability of the demand arrived at distribution center at moment $t$ is

$$ \delta_j = \frac{\lambda_0^j}{Q_r} e^{-\lambda_0 \tau}, \quad j \mid Q_r \in \mathbb{N}^+ \quad (3.2) $$

The orders arrive as Poisson flow, thus the probability that the number of orders arriving at distributor equal to $j$ within a time interval $t$ $\tau$ is

$$ \omega_j(\tau) = \frac{(\lambda_0 \tau)^j}{j!} e^{-\lambda_0 \tau}, \quad j = 0,1,2,\dots \quad (3.3) $$

As the batch of each order is $Q_r$, the probability of the demand equal to $d$ within a time interval $t$ $\tau$ is

$$ P\{D(\tau) = d\} = \omega_{d/\tau}(\tau) = \frac{(\lambda_0 \tau)^{d/\tau}}{d^{d/\tau}} e^{-\lambda_0 \tau} \quad (3.4) $$

$\mu$ presents for average demand per unit time $d/\tau$ of distribution center, then

$$ \mu = \lambda_0 Q_r = N \lambda_r \quad (3.5) $$

### B. Analysis of inventory level

Suppose the holding inventory at moment $t$ is $\{IP(t), t \ge 0\}$ and inventory level is $\{IL(t), t \ge 0\}$. So we can get the following properties,

$$ IL(t + L_0) = IP(t) - D(t, t + L_0) = IP(t) - D(L_0) \quad (3.6) $$

$\{IP(t), t \ge 0\}$ is a random process with continuous parameters and discrete values, and the scope of $IP(t)$'s value is $\{s_0 + Q_r, s_0 + 2Q_r, \dots, S_0\}$. We can see from the assumption (5) that, $s_0$ and $S_0$ are integral multiples of $Q_r$, assume $s_0 = k_0Q_r, S_0 = K_0Q_r$, and $k_0, K_0 \in \mathbb{N}^+$ here.

Suppose $\rho_j$ to be the probability that the holding inventory $IP(t)$ reach $j$, $j = mQ_r, k_0+1 \le m \le K_0$. At the beginning of each order cycle, the orders of distribution centers ensure the holding inventory reach $S_0$, so . In an order cycle, each time the arrival of the demand will cause transitions of holding inventory, and each holding inventory point can be reached at most once. Obviously, when $IP(t) = f > j$, if the next demand arriving at distribution center is $f - j$, then
the holding inventory $IP(t) = j$, so, there is

$$ \rho_j = \sum_{f>j} \rho_f \delta_{f-j} = \sum_{i=m+1}^{K_0} \rho_{iQ_r} \delta_{(i-m)Q_r}, \quad \text{that is} \quad \rho_j = \sum_{i=\frac{j}{Q_r}+1}^{K_0} \rho_{iQ_r} \delta_{iQ_r-j} \qquad (3.7) $$

Here, $j = s_0 + Q_r, s_0 + 2Q_r, \dots, S_0$. Only one solution can be get by solving equation (3.7), denoted by $\rho_j^*(j = s_0 + Q_r, s_0 + 2Q_r, \dots, S_0)$. In fact, $\rho_j^*$ also can be seen as the expected frequencies that holding inventory reaches $j$ in one order cycle. Therefore, the frequency of average total demand in one order cycle is $\sum_j \rho_j^* = \rho_{s_0+Q_r}^* + \rho_{s_0+2Q_r}^* + \dots + \rho_{S_0}^* = \sum_{i=k_0+1}^{K_0} \rho_{iQ_r}^*$. As orders to distribution center obey

compound Poisson process, and the average time intervals of demands to systems are always the same, therefore, under the stable state, the probability distribution of inventory location can be written as follow,

$$ P\{IP = y\} = \frac{\rho_y^*}{\sum_j \rho_j^*} = \frac{\rho_y^*}{\sum_{i=k_0+1}^{K_0} \rho_{iQ_r}^*}, \quad y = s_{0+}Q_r, s_0 + 2Q_r, \dots, S_0 \qquad (3.8) $$

### C. Holding cost and backorder cost of anticipation inventory

For the given inventory level $IL$, inventory holding cost can be expressed as $h_0 \cdot E(\max\{IL,0\})$, and $\max\{IL,0\}$ as $(IL)^+$, then inventory cost can be expressed as $h_0 E(IL)^+$. In addition, the delayed delivery cost can be expressed as $b_0 \cdot E(\max(-IL,0))$, and $\max(-IL,0)$ as $(IL)^-$, then delayed delivery cost can be expressed as $b_0 E(IL)^-$. Assume inventory position of the distribution center is always $IP = y$, the following can be get from property(3.6).

$$ P\{IL = x\} = P\{D(L_0) = y - x\}, \quad x \le y \qquad (3.9) $$

When $IP = y$, $G(y)$ represents the sum of anticipation inventory cost and anticipation delayed delivery cost, so

$$ G(y) = b_0 E(IL)^- + h_0 E(IL)^+ \qquad (3.10) $$

Because $(IL)^- = (IL)^+ - IL$,

$$ E(IL)^- = E((IL)^+ - IL) = E(IL)^+ - E(IL) \qquad (3.11) $$

$$ \text{So, } G(y) = -b_0(y - \mu L_0) + (h_0 + b_0) \sum_{i=1}^{\frac{y}{Q_r}} i Q_r \cdot P\{D(L_0) = y - i Q_r\} \qquad (3.12) $$

### D. Anticipation total cost per unit time of the distribution center

We can know from the analysis of (B) that the average number of holding inventories' changes in one order cycle is $\sum_j \rho_j^* = \sum_{i=k_0+1}^{K_0} \rho_{iQ_r}^*$, and every demand arrival can lead to changes of inventory levels, thus the average length of one order cycle is

$$ T_0 = \sum_{i=k_0+1}^{K_0} \rho_{iQ_r}^* = \sum_{i=\frac{s_0}{Q_r}+1}^{\frac{s_0}{Q_r}} \rho_{iQ_r}^* \qquad (3.13) $$
From the above analysis, we can know that anticipation total cost per unit time of the distribution center is

$$ATC_0(s_0, S_0) = \frac{A_0}{T_0} + \sum_{k=s_0/Q_r+1}^{s_0/Q_r} P\{IP = kQ_r|G(kQ_r), in \frac{s_0}{Q_r}, \frac{S_0}{Q_r} \in N^+\} \quad (3.14)$$

## 3.2 An inventory model of retailers

### A. The analysis on lead time of retailers’ order

We can know from assumption(6) that the lead time of retailers’ order consists of two parts, of which the fixed shipping time from distributor to different suppliers $L_r$ is constant while random delay of retailers’ orders at distribution center $L_d$ is a random variable. Namely,

$$L_r = L_t + L_d \quad (3.15)$$

Just because the uncertainty distribution of $L_d$, its mean value is often used to represent the random variable in practice, we can get equation (3.16) from Little’s formula in Queuing Theory,

$$\bar{L}_d = E(\bar{I}_{L_0})^{-} / \mu_0 \quad (3.16)$$

Here, the numerator represents the average backorder level of distribution center, the denominator represents the average demand per unit time of distribution center. Substitute equation (3.5) and equation (3.11) into equation (3.16),

$$\bar{L}_d = \frac{1}{N\lambda_r} \left\{ \sum_{k=s_0/Q_r+1}^{s_0/Q_r} P\{IP = kQ_r\} \cdot \left[ \sum_{i=1}^{k} iQ_r \cdot P\{D(L_0) = (k-i)Q_r\} - (kQ_r - N\lambda_r L_0) \right] \right\} \quad (3.17)$$

$$\text{And } \bar{L}_r \text{ satisfies the following conditions, } \bar{L}_r = L_t + \bar{L}_d \quad (3.18)$$

### B. Holding cost of inventory

From assumption (3), the demand during lead time $x_r$ is a discrete random variable taking positive integers, provided the cumulative distribution function is $\phi(x)$, the probability distribution is $P(X_r = x) = \phi(x)$, the expected value (average) is $l$. The average inventory in one cycle is $\bar{I}_r = \frac{1}{2}[R_r - l + Q_r + R_r - l] = \frac{Q_r}{2} + R_r - l$.

As $x_r \sim P(\lambda_r)$, according to the characteristics of Poisson distribution, the average demand per unit time is, $\lambda_r$ and the order quantity is $Q_r$, then the average order times per unit time is $\frac{\lambda_r}{Q_r}$, thus the average length of one order cycle is $\frac{Q_r}{\lambda_r}$. Therefore, the average inventory holding cost of the retailers in one cycle is,

$$C_r^l = h_r \frac{Q_r}{\lambda_r} \left( \frac{Q_r}{2} + R_r - l \right) \quad (3.19)$$

### C. Backorder cost

Backorders occur when demands in the lead time exceed a reorder point, so the average value of backorders in one cycle is,

Here, $\phi(x)$ is the probability density function of demand $x_r$ in the lead time, and $\phi(x)$ is a cumulative distribution function of $X_r$.
Thus the average total backorder cost in a cycle is

$$C_r^S = b_r \left\{ \sum_{x=R_r}^{+\infty} x\phi(x) - R_r(1 - \phi(R_r)) \right\} \quad (3.20)$$

### D. The average total cost per unit time

From equation (3.19) and (3.20), the average total cost of the retailers and since the average time of one ordering cycle is $\frac{Q_r}{\lambda_r}$, the average total cost of retailers per unit time is

$$ATC_r(R_r, Q_r) = \frac{\lambda_r}{Q_r} A_r + h_r \left( \frac{Q_r}{2} + R_r - l \right) + b_r \frac{\lambda_r}{Q_r} \left\{ \sum_{x=R_r}^{+\infty} x\phi(x) - R_r[1 - \phi(R_r)] \right\} \quad (3.21)$$

Retailers' demand per unit of time $x_r \sim P(\lambda_r)$, as convolution formula of Poisson distribution, demands during lead time meet $X_r \sim P(\bar{L}_r \lambda_r)$, so we have the followings in equation (3.21),

$$\phi(x) = \frac{(\bar{L}_r \lambda_r)^x}{x!} e^{-\bar{L}_r \lambda_r} \quad (3.22)$$

$$\phi(R_r) = \sum_{x=0}^{R_r} \phi(x) \quad (3.23)$$

$$l = \bar{L}_r \lambda_r \quad (3.24)$$

$\bar{L}_r$ can be get from (3.17) and (3.18).

Refer to equation (3.14) and (3.21), the function of anticipation total cost of the overall two-echelon distribution system per unit time can be expressed as follow,

$$ATC(s_0, S_0, R_r, Q_r) = ATC_0(s_0, S_0, Q_r) + N \cdot ATC_r(R_r, Q_r, \bar{L}_d(s_0, S_0)) \quad (3.25)$$

The goal of the model is to find the optimal $s_0^*, S_0^*, R_r^*, Q_r^*$, so that the anticipation total cost of the overall system per unit time to be minimum. So the problems to be solved by the model can be described as

$$\min \quad ATC = ATC(s_0, S_0, R_r, Q_r) \quad (3.26)$$

$$\begin{array}{l} \left\{ \begin{array}{l} R_r, Q_r \in N \\ 0 < R_r + 1 \le Q_r \\ s_0/Q_r, s_0/Q_r \in N^+ \\ s_0 < S_0 \end{array} \right. \\ \text{s.t.} \\ \text{Equation (3.1) } \sim (3.25) \end{array}$$

## 4 Solution method for model based on Particle Swarm Optimization (PSO)

### 4.1 Coding and fitness function

Set particle size M=20, the maximum iteration Kmax=30, adopt real number coding on the respective components of each vector's position Xᵢ, the position vector
$X_i$ contains four dimensions, corresponding $s_0$, $S_0$, $R_r$, $Q_r$ respectively, so the current position of each particle forms a solution vector, that is $X_i = (s_i, S_i, R_i, Q_i)$, $i = 1, 2, \dots, 20$.

The goal of the two-echelon system joint inventory optimization is to optimize the total cost of distribution centers and the retailers, therefore, fitness function of algorithm can be defined as the sum of the cost of both, that is

$$ \text{fitness} = \text{ATC}_0(s_i, S_i, Q_i) + N \cdot \text{ATC}_r(s_i, S_i, R_i, Q_i) $$

## 4.2 The solving process of Particle Swarm Optimization (PSO)

Step1: Its scale $M = 20$, Maximum iteration $k_{max} = 30$.

Step2: Determine the initial position of each particle $X_i^{(0)} = (s_i, S_i, R_i, Q_i)$, $i = 1, 2, \dots, M$, the initial velocity of each particle $v_i^{(0)} = (v_i^S, v_i^S, v_i^R, v_i^Q)$.

Step3: Measure the fitness of each particle based on (4.1) $f_i^{(0)}$, the optimal position of the each particle $p_i^{(0)} = f_i^{(0)}$.

Step4: Find out the global optimum $p_g^{(0)}$ from $p_g^{(0)} = \min\{p_1^{(0)}, p_2^{(0)}, \dots, p_M^{(0)}\}$.

Step5: Let iterations $k=0$ and $k< -k+1$, update particle velocity $v_i^{(k)}$ and particle position $X_i^{(k)}$ according to the following three equations. $X_i^k = X_i^{k-1} + v_i^{(k)}$

$$ v_i^{(k)} = w^k v_i^{(k-1)} + c_1 r_1 (p_i - X_i^{(k-1)}) + c_2 r_2 (p_g - X_i^{(k-1)}), \quad w^k = w_{\max} - \frac{w_{\max} - w_{\min}}{k_{\max}} \cdot k $$

Here, $r_1, r_2$ are random numbers among [0,1], accelerated factor $c_1, c_2 = 0.7$, the maximum value of inertia factor $w_{\max} = 0.9$, minimum value $w_{\min} = 0.7$.

Step6: Measure the fitness of each particle fitness $f_i^{(k)}$ according to fitness function (4.1).

Step7: For each particle, compare its fitness $f_i^{(k)}$ with its historical optimal position $p_i^{(k-1)}$, place it as the current optimal position $p_i^{(k)}$ if it's better than $p_i$.

Step8: For each particle, compare its optimal position $p_i^{(k)}$ with the global optimal location of the whole group $p_g^{(k-1)}$ place it as the optimal location of the whole group $p_g^{(k)}$ if it's better than $p_g^{(k)}$, then $p_g^{(k)} = \min\{p_1^{(k)}, p_2^{(k)}, \dots, p_M^{(k)}\}$.

Step9: Check the termination conditions, jump to step 6 if $k < k_{max}$; Otherwise, the maximum iteration, termination.

## 5 Examples

The distribution center face 40 retailers, that is $N = 40$, retailer's daily demand obey Poisson distribution, average demand $\lambda_r = 9$ units. Daily per holding cost of the distribution centers and retailers is $h_0 = h_r = 2$ yuan, unit backorder cost of the distribution centers is $b_0 = 4$ yuan, unit backorder cost of the retailers is $b_r = 8$ yuan,
ordering fee of retailers each time is $A_r = 96$ yuan, ordering fee of the distribution centers each time $A_0 = 1200$ yuan, the lead time of distribution center is $L_0 = 2$ days, transportation time from the distribution center to the retailers is $L_r = 1$ day.

Take various parameters in the above application examples into the two-echelon distribution system and solve the problem based on algorithm of Particle Swarm Optimization (PSO) with the software MATLAB (R2008a). We will record once the current optimal solution every three times of iterations, the results refer to Table 1.

From the table1 we can clearly see the convergence process of particle swarm, the difference of the optimal fitness between adjacent two iterations is gradually reduced, and the gap between the average fitness of each particle and the optimal value is gradually narrowing, indicating that each particle converges to the optimal position, its visible that the algorithm has preferable convergence.

<table><caption>Table 1 Example 5.1 Iterative process records algorithm of PSO</caption><thead><tr><th>Number</th><th>s<sub>0</sub></th><th>S<sub>0</sub></th><th>R<sub>r</sub></th><th>Q<sub>r</sub></th><th>Optimal fitness</th><th>Average fitness</th></tr></thead><tbody><tr><td>0</td><td>663</td><td>1938</td><td>9</td><td>51</td><td>2755.69</td><td>7947.14</td></tr><tr><td>3</td><td>480</td><td>1440</td><td>11</td><td>40</td><td>2672.36</td><td>6148.03</td></tr><tr><td>6</td><td>480</td><td>1440</td><td>11</td><td>40</td><td>2672.36</td><td>4200.38</td></tr><tr><td>9</td><td>396</td><td>1672</td><td>11</td><td>44</td><td>2671.7</td><td>3345.47</td></tr><tr><td>12</td><td>540</td><td>1575</td><td>8</td><td>45</td><td>2651.01</td><td>2957.4</td></tr><tr><td>15</td><td>540</td><td>1575</td><td>8</td><td>45</td><td>2651.01</td><td>2870.4</td></tr><tr><td>18</td><td>484</td><td>1672</td><td>8</td><td>44</td><td>2648.03</td><td>2869.76</td></tr><tr><td>21</td><td>495</td><td>1620</td><td>9</td><td>45</td><td>2645.65</td><td>2711.55</td></tr><tr><td>24</td><td>484</td><td>1628</td><td>9</td><td>44</td><td>2644.08</td><td>2762.3</td></tr><tr><td>27</td><td>484</td><td>1584</td><td>9</td><td>44</td><td>2643.91</td><td>2884.25</td></tr><tr><td>30</td><td>473</td><td>1591</td><td>9</td><td>43</td><td>2643.52</td><td>2674.47</td></tr></tbody></table>

# 6 Summary

This paper selects an important part of the multi-echelon inventory system-distribution system as the research object, the optimization of the two-echelon inventory system is studied under centralized control strategy. Distribution centers adopt (s,S) ordering strategy while retailers adopt (R,Q) ordering strategy, the lead time of distribution centers to retailers is a random variable while each retailer's demand is independent Poisson process, the authors construct an inventory model of a two-echelon distribution system under this assumption, and solve the problem based on algorithm of Particle Swarm Optimization (PSO) with inertia weights to achieve effective exploration on swarm intelligence technology on inventory optimization problems.
References

1. Federgruen A, Zipkin P H. Approximations of dynamic multi-location production and inventory problems [J]. Management Science, 1984(30):69~84

2. Axsäter S. A simple decision rule for decentralized two-echelon inventory control[J]. International Journal of Production Economics, 2005(93-94):53~59

3. Olsson R J. Hill RM. A two-echelon base stock inventory model with Poisson demand and the sequential processing of orders at the upper echelon[J]. European Journal of Operational Research, 2007(177): 310~324

4. Thangam A, Uthayakumar R. A two-level distribution inventory system with stochastic lead time at the lower echelon[J].The International Journal of Advanced Manufacturing Technology, 2009(41):1208~1220