(I) SIMULATION.Our first approach is to estimatepˇj using the transfer statistics from the central supplier to location j obtained in simulation runs of the multi-product model. Thereafter these estimated values are inserted into the analytical formulas of the simpler model. Clearly, this is not a practical recipe, but a way to get insights into possible similarities of both models. Great differences e.g. in inventory sizes or queue lengths would recommend not to use the simpler model. But fortunately enough, these differences are small.
For this purpose we construct a fictional system with deterministic routing with two locations, J ={1,2}, and parameters λ1 = 2, λ2 = 5, µ1 = 4, µ2 = 10, b1 = 4, b2 = 12 andν in a range from0.5 to 40with step size 0.25. We have chosen the service ratesµj
greater than the demand rates λj to keep the system ergodic, and we have chosen the demand rates λj and the base stock levels bj in such a way that the ratios λ1/λ2 and b1/b2 are different. A larger base stock level b2 is attributed to a system with a larger demand stream λ2. Furthermore, we have chosen a very large run time T = 100000to obtain the results close to the steady state solution in a single simulation run. We stored the numberdj of finished items sent from the central supplier to locationj. The results of the simulation are plotted in Figure 2.4.1.
From Figure 2.4.1(a) we see, ifν is small,d1/(d1+d2)seems to approachb1/(b1+b2), and ifν is large,d1/(d1+d2)approaches λ1/(λ1+λ2). The intuitive explanation of this behaviour with the help of Figures 2.4.1(b) - 2.4.1(d) is as follows: If the central supplier is much slower than the production systems at the locations, then the central supplier’s queue is almost always full and the inventories are almost always depleted. This means, in the queue of the central supplier there are approximatelyb1 orders from location1and b2 orders from location2 in random order.
If the central supplier is much faster than the locations, then the orders pass the central supplier almost immediately. The order streams from locations1and 2 behave similarly to the superposition of two stochastically independent Poisson streams with ratesλ1 and λ2. When two independent Poisson streams pass the central supplier with no delay, this is stochastically the same, as to input a Poisson stream with non-distinguishable orders to the central supplier with a rateλ1+λ2 and then randomly decide of which type (location) it is: of type 1 with probabilityλ1/(λ1+λ2)and of type 2 with probability λ2/(λ1+λ2).
Figures 2.4.1(e) - 2.4.1(f) show that the analytically obtained average queue size almost perfectly matches the simulated (true) values of the average queue sizes at the locations.
Furthermore, we can see that the average queue size is independent of the service rate ν at the central supplier (see equation (2.3.7)), which is predicted by the product form stationary distribution (2.3.5).
The comparison of the performance metrics of the simulated more complicated model with those of our analytically obtained results shows that the analytically obtained values can be used as an approximation for the multi-product and more complicated system’s metrics. As noted before, the much more extensive comparisons of [GKAR12] support such substitutions of complex systems by suitably chosen product form systems as well.
2.4. Comparison with a more complex model with location specific items
ν 0.24
0.25 0.26 0.27 0.28
0 10 20 30 40
b1(b1+b2)
λ1(λ1+λ2) d1(d1+d2)
(a)d1/(d1+d2) — relative departure portion of finished items for location 1 from the central supplier in the system with deterministic rout-ing
ν 0
5 10 15
0 10 20 30 40
simulation approximation
(b) Average queue size at the central supplier
ν 2
4 6 8 10
0 10 20 30 40
simulation approximation
(c) Average inventory at location 1
ν 2
4 6 8 10
0 10 20 30 40
simulation approximation
(d) Average inventory at location 2
ν 0.2
0.4 0.6 0.8 1.0
0 10 20 30 40
simulation approximation
(e) Average queue size at location 1
ν 0.2
0.4 0.6 0.8 1.0
0 10 20 30 40
simulation approximation
(f) Average queue size at location 2
Figure 2.4.1.: Comparison between the simulated results of the system with deterministic routing and the analytic results of the system with random routing
(II) ITERATIVE ALGORITHM.In our second approach we use a queueing model to set the values of the routing probabilitiespj,j∈J. In this procedure the routing prob-abilities are determined as being proportional to the effective arrival rates from location j at the central supplier. We consider the following model, where an admitted customer at locationj represents an outstanding replenishment order from locationj. The process time of such orders is complex (consisting of waiting and service times at the central sup-plier). Therefore, the admitted customer’s service time is modelled as a random variable with a general distribution. Consequently, each location j, j ∈ J, is considered as an M/G/bj/0-Erlang-loss system3 as depicted in Figure 2.4.2.
The arrival rate λj, j = 1, . . . , J, is diminished by the loss probability qj when all bj
units are on order at the central supplier. The value ofqj will be determined iteratively.
Furthermore, the total arrival rateλJ+1 at the central supplier equals PJ
j=1λj·(1−qj).
Lost customers
1 2
Figure 2.4.2.: Location j approximated as an M/G/bj/0-FCFS queue The central supplier is modelled as anM/M/1/(b−1)-FCFS queue withb=P
j∈Jbj, service rate ν, and arrival rate λJ+1 of orders generated by admitted customers at the localM/G/bj/0queues (Figure 2.4.3).
Soujourn time T Lost
customers
Figure 2.4.3.: Central supplier approximated as an M/M/1/(b−1)-FCFS queue To determine the arrival rate at the central supplier, the blocking probabilities (loss probabilities) qj of the M/G/bj/0 queues need to be known and to determine these blocking probabilities we need to know the sojourn time T of the replenishment orders at the central supplier. This can be solved iteratively where the algorithm will stop if the blocking probabilitiesqj remain unchanged under further iterations. As a result the routing probabilities pj, j ∈ J, can be calculated because they are proportional to the effective customer arrival rates at the central supplier. That is
pj = (1−qj)·λj P
k∈J(1−qk)·λk.
3 bjservice channels, no waiting room, Poisson arrivals, general service time distribution
2.4. Comparison with a more complex model with location specific items
AlgorithmCalculation ofpj,j∈J
Input: number of locations J
service rate of the central supplier ν arrival rates at the locationsj∈J λj
base stock levels at the locationsj∈J bj
stop criterion ε
Output: routing probabilitiespj,j∈J
Initialize: blocking probabilitiesqj= 0, ∀j ∈J
(1)Calculate the effective arrival rate of replenishment orders at the central supplier λJ+1=
J
X
j=0
(1−qj)·λj
and the average sojourn timeE[T]of a replenishment order at the central supplier
E[T] = 1
Pb
`=0
λ
J+1
ν
`·
b−1
X
n=0
λJ+1 ν
n
·n+ 1
ν withb=X
j∈J
bj.
(2)Determine the new blocking probabilitiesqj(new)at the locationsj∈J qj(new)= (λj·E[T])bj
Pbj
n=0(λj·E[T])n.
(3)If X
j∈J
|qj−qj(new)|> ε, then
qj ←qj(new) ∀j∈J and return to(1), else calculate the routing probabilitiespj,j∈J,
pj= (1−q(new)j )·λj
P
j∈J(1−q(new)j )·λj
.
The iterative algorithm can be modelled in R. The R code is presented in Appendix B.1 on page 267.
We can use the simulation results of the multi-product system again and compare the values ofpj with the results of the iterative procedure. We again consider two locations J = {1,2} with the same parameter values, except we do not need the service rates µj
for the algorithm. We have chosen a stop criterion ε= 0.001. Furthermore, we stopped the algorithm when no convergence was reached within500iterations.
The results are plotted in Figure 2.4.4. We see that the resulting routing probabilities pj of the algorithm show a good approximation if the service rate of the central supplierν is larger than or equal to some valueν∗. In Figure 2.4.4 the value ofν∗ is approximately 5.5. In the grey area of Figure 2.4.4, the iterative algorithm did not satisfy the stop criterionεin less than 500 iterations for about 85% of the instances. The resulting values ofpj are not good approximations for the actual values. A series of experiments withλ1
andλ2 in a range from 0.5 to 20 with step size 0.25 supports our conjecture that such a valueν∗ exists in general. The existence of ν∗ is an open problem.
ν p1
0.24 0.25 0.26 0.27 0.28
0 5 10 15 20 25 30 35 40
algorithm simulation
Figure 2.4.4.: Comparison between the simulated results and the results of the algorithm We found it necessary to investigate the relationship betweenν∗ and other parameters of the system,4 i.e. to determine the range of ν-values where the algorithm converges.
Forλ1 =λ2in a range from0.5to20with step size0.25(and all other parameters equal to the above) we see in Figure 2.4.5 that there is an approximately linear relationship between the arrival rates at the locations andν∗.
λ1
ν*
0 10 20 30
0 5 10 15 20
Figure 2.4.5.: Relationship betweenν∗ and λ1 =λ2
4The investigation of the relationship between ν∗and other parameters of the system is an improved version of [OKD16].
2.4. Comparison with a more complex model with location specific items
We mention that the optimal base stock levels bj, which are the final decision vari-ables, are determined in Theorem 2.5.2. The values pj are needed in γj = νpλj
j to cal-culate P(Yj = 0) and E(Yj) in the cost function gj(bj). The decision in the grey area seems to be relatively robust because ν is small in this area. However, it depends on the combination of λj and ν. 5We can use the simulation results of the multi-product system again and compare the values with the results of the iterative procedure with the routing probabilitiespj in the grey area where ν < ν∗. In our example the value of ν∗ is approximately5.5. The results are plotted in Figure 2.4.6.
ν γ1
0.2 0.4 0.6
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(a)γ1= νpλ1
1
ν γ2
0.5 1.0 1.5 2.0
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(b) γ2= νpλ2
2
ν P(Y1=0)
0.4 0.5 0.6 0.7 0.8 0.9
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(c)P(Y1= 0)
ν P(Y2=0)
0.0 0.2 0.4 0.6 0.8
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(d) P(Y2= 0)
ν E(Y1)
0.2 0.4 0.6 0.8 1.0 1.2
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(e)E(Y1)
ν E(Y2)
0 2 4 6 8 10
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
algorithm simulation
(f)E(Y2)
Figure 2.4.6.: Comparison between the simulated results and the results of the algorithm in the grey area (i.e. ν in a range from0.5 to 5.25 with step size 0.25)
5The following investigations are an improved version of [OKD16].
Furthermore, we have plotted the relative errors |value of simulation−value of algorithm|
value of simulation for γj,P(Yj = 0)andE(Yj) in Figure 2.4.7. We see the greatest relative error forP(Y2 = 0) in Figure 2.4.7(d) which was not clearly visible in Figure 2.4.6(d). However, it seems to be relatively robust since the value of P(Y2 = 0) is small in this area (see Figure 2.4.6(d)). Consequently, the resulting values forγj,P(Yj = 0) and E(Yj) show a good approximation for the actual values in our example.
However, the decision of the optimal base stock levels also depends on the specific cost values. Therefore, further studies are still needed to make a statement about the robustness of the optimal base stock levels in the grey area.
ν 0.05
0.10 0.15
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(a)γ1=νpλ1
1
ν 0.05
0.10 0.15
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(b) γ2=νpλ2
2
ν 0.1
0.2 0.3 0.4
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(c)P(Y1 = 0)
ν 0.1
0.2 0.3 0.4
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(d) P(Y2= 0)
ν 0.05
0.10 0.15
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(e)E(Y1)
ν 0.05
0.10 0.15
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(f) E(Y2)
Figure 2.4.7.: Relative errors |value of simulation−value of algorithm|
value of simulation