A Simulation Based Solution Approach for the Capacitated Lot-Sizing Problem With Remanufacturing

1. Introduction

The production planning with consideration of recycling options has received increasing attention in the last few years [1]. The literature describes two different categories of recycling production planning problems. In the first category (Cat. 1), the given external demand of the products has to be satisfied by remanufactured or newly produced goods Figure 1. Only two types of inventories are regulated: the inventory of new products and the inventory of remanufactured products, and this is the category that has been thoroughly investigated by numerous authors. It is assumed that new and remanufactured products are identical, and hence, both are regarded as serviceables, i.e., entities that can be used to satisfy demands. In the second category (Cat. 2) the additional feature of this problem is that in each period a deterministic amount of returned items enters the system. These returns can be remanufactured and used to satisfy the demand for remanufactured products in addition to the regular production of new items. Thus, there are three types of inventory: new product inventory, remanufactured product inventory and returned product inventory (Figure 1). This situation is typical in paper manufacturing. Demand for recovered paper has increased enormously worldwide, making separate collection more attractive again. Peter Probst, CEO of LEIPA Group, Germany (Viadrina University, 05.2019) explains the advantage: "paper recycling proves to be economically and ecologically efficient, because the production of recycled paper requires considerably less water and considerably less energy than the production of virgin fiber papers. In addition, the transport distances are generally shorter as well. Recycled paper comes from Germany, while the fresh fibers and fresh fiber papers are also imported to Germany in larger proportions". Paper production is also energy intensive. It takes as much energy to produce one ton of paper from fresh wood fibers as it does to produce one ton of steel ([2]).

2. Literature Review

Publications regarding Category 1: [3] analyzes several CLSP-RM formulations with the column generation method, where he uses the idea of Garfinkel and Nemhauser (1969), namely the Set Partitioning Problem (SPP). [4] compared MIP approaches to the economic lot-sizing with remanufacturing (ELSR) and proposes a shortest path formulation. [5] studied another multi-item variant of the problem and propose a variable neighborhood search heuristic. The polynomial-time heuristics of [6] also belongs to this category, where there is also an excellent compilation of the methods developed by [7,8,9,10,11].

Publications regarding category 2: The only publication known to us is that of [12] where he presents a Lagrange relaxation to solve the problem. Based on Zhang’s model we formulate the same problem but taking more general capacity constraints and doing more extensive numerical experiments.

The capacitated lot sizing problem model is NP-hard. The proof of this statement is to be found in [13,14]. [15] has shown that even the two-item problem with constant capacity is NP-hard. Many methods have been proposed to solve the CLSP  [16]. In general, these heuristics do not use random numbers. However, many heuristics cannot guarantee the generation of a feasible solucion  [15]. These heuristic methods used to solve a CLSP class problem are generally based on the definition of a multitude of rules derived from demand, capacity, setup costs, inventing costs and productions costs e.g. [9,10,17].The solution margin varies between 1% and 10% [12,17].

The contribution of this work is: independently of the input parameters of the problem, we organize the production in the framework of feasibility and randomly. A simple simulation with quasi-random numbers ( $\mathbb{QRN}$ ) will be presented. The heuristic will always generate a feasible solution, if the problem is solvable. After simulating N production plans we choose the cheapest one. The scheduling is very simple, the results are quite acceptable considering the size of the presented problems.

It is well known that the $\mathbb{QRN}$ of [18,19,20] have significant properties,[21]. The most important of these are: they are uniformly distributed. These are some allow us a systematic search of the solution in the set of admissible points with a few $\mathbb{QRN}$. We implement a search algorithm with Halton’s $\mathbb{QRN}$ and compare it with numerical solutions by the Gurobi solver.

In Section 3 we propose a new mathematical programming formulation for the problem. In Section 4 we present a relaxation of the problem without capacity constraints (Model A) and we will see that the solution obtained for the chosen class of tasks is very close to the optimal solution of the original problem. In Section 5, we present a heuristic solution method (Model B) for another class of tasks. In Section 6, we report the results of computational experiments. In Section 7 we present some conclusions and suggested directions for further research.

3. Problem description

We assume that a factory produces two types of products, one manufactured from raw materials and the other remanufactured from collected used products. The demands for these two products are separate, deterministic, and time varying during a finite planing horizon, and should be satisfied without backlogging. The costs consist of fixed setup cost, linear production cost proportional to the production quantity, and linear inventory holding costs. All cost components are considered for both manufacturing and remanufacturing activities per unit and periode  [12]. The following Table 1 summarize the notation used in this paper.

We make the following assumptions, which are necessary for the feasibility of the problem

• The demand for new products and remanufactured products are separate and backlog is not allowed.
• The manufacturing capacity is sufficient to meet the demands in each period, in particularly we have: the capacity can satisfy the demands for new products and remanufactured products simultaneously, i.e.

  

  (1)
• Initial and end inventory stocks are zero, 1 i.e.

  

  (2)
  The demand will be fully satisfied if the final inventories are zero.
• The quantity of returned products can satisfy the demand for remanufactured products i.e.

  

  (3)
• In economic terms, inventory holding cost of returned products is less than that of remanufactured products.

  

  (4)

  This hypothesis can be found in [22] too.

Hence, the problem can be formulated as

(5)

subject to

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

The objective function (5) minimizes the sum of setup cost, production cost, and inventory cost for new products and remanufactured products in all periods. Constraints (6),(7) and (8) are the inventory balance constraints for new products, remanufactured products and returned products. Constraints (9) represent capacity constraints for manufacturing and remanufacturing activities. Constraints (10), (11) allow production only with the according setups (i.e ${\alpha }_t^n=1,{\alpha }_t^r=1.$) (12) and (13) are the standard integrality and non-negative constraints.

3.1. Rewriting the optimization problem

$i_t^n,i_t^r,i_t^{re}$ are replaced in the objective function by

(14)

(15)

And with a bit of algebraic transformations, we obtain

$$\begin{array}{c}\hfill \sum _{t=1}^T(tn{x}_t^n+tntn)=(p_t^n+h_j^n)x_t^n-\sum _{t=1}^T\left(h_t^n\sum _{k=1}^t{d}_k^n\right)\\ \hfill \sum _{t=1}^T(tr{x}_t^r+trtr+h_t^{re}tre)=(tr+(h_j^r-h_j^{re}))x_t^r\\ \hfill +(-h_t^r{d}_k^r+h_t^{re}{r}_k^{re})\end{array}$$

With the notation

$$\begin{array}{cc}\hfill & v_t^n=p_t^n+h_j^n,\hfill \\ \hfill & v_t^r=p_t^r+(h_j^r-h_j^{re})\hfill \\ \hfill & K=\sum _{t=1}^T\left\{-h_t^n\sum _{k=1}^t{d}_k^n-h_t^r\sum _{k=1}^t{d}_k^r+h_t^{re}\sum _{k=1}^t{r}_k^{re}\right\}.\hfill \end{array}$$

our model is finally

(16)

Subject to

(17)

(18)

(19)

$$\begin{array}{cc}\hfill x_t^n+x_t^r& \le c_t,\forall t=1,\dots ,T\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill x_t^n& \le d_{[T]}^n{\alpha }_t^n,\forall t=1,\dots ,T\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill x_t^r& \le d_{[T]}^r{\alpha }_t^r,\forall t=1,\dots ,T\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill x_t^n,x_t^r& \ge 0,\forall t=1,\dots ,T\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\alpha }_t^n,{\alpha }_t^r& \in \{0,1\}\hfill \end{array}$$

(24)

4. Model A (Relaxation)

According to our conditions (1) and (3) the general problem has a solution, if the total pro-period demand is less than the pro-period capacity (i.e.$d_t<=c_t$). However there are situations, where conditions (1) and (3) are satisfied but $d_t<=c_t$ is not satisfied. We need a feasibility routine which ensures that all demand is satisfied without backlogging. Indeed there are periods (or could be) in which total demand exceeds total capacity. In this case some inventory will have to be build up in earlier periods which slack capacity. We explain how to shift excess demand to earlier periods in which slack capacity is available. We use and complement the idea of [23].

$$\begin{array}{cc}\hfill {\tilde{w}}_T& =0\hfill \\ \hfill {\tilde{w}}_t& =max\{d_{t+1}-c_{t+1}+{\tilde{w}}_{t+1};0\},t=T-1,\cdots ,1.\hfill \end{array}$$

We define

$$\begin{array}{cc}\hfill w_1& ={\tilde{w}}_1\hfill \\ \hfill w_t& ={\tilde{w}}_t-{\tilde{w}}_{t-1},t=2,\cdots ,T.\hfill \end{array}$$

It is easy to see that the sum $w_{\left[T\right]}$ is null.

$$\begin{array}{cc}\hfill w_{\left[T\right]}& =w_1+w_2+\cdots +w_T\hfill \\ \hfill & ={\tilde{w}}_1+({\tilde{w}}_2-{\tilde{w}}_1)+\cdots +({\tilde{w}}_T-{\tilde{w}}_{T-1})\hfill \\ \hfill & ={\tilde{w}}_T=0.\hfill \end{array}$$

Remark 1.The vector $w^T=(w_1,\cdots ,w_T)$ is very useful. Because $w_t$ gives the amount of stock to accumulate ($w_t>0$) or reduce ($w_t<0$) in each period so that production does not exceed available capacity pro period. And allows us to determine a good permissible solution.

Table 3. Original Data.

Table 3. Original Data.

We see that the demand transformation done in Table 4 is a valid (permissible) solution to the problem (16)-(24). There are many ways to make this transformation, for this reason to transform the demand in an optimal way we formulate the following problem.

Table 4. Demand Transformation.

Table 4. Demand Transformation.

$$\begin{array}{cc}\hfill \phi (x,y)& \sum _{t=1}^T=[s_t^n{\alpha }_t^n+v_t^n{x}_t^n]\hfill \\ & +[s_t^r{\alpha }_t^r+v_t^r{x}_t^r]+K\to min\hfill \end{array}$$

(25)

Subject to

(26)

(27)

(28)

$$\begin{array}{cc}\hfill x_t^n+x_t^r-d_t^n-d_t^r& =w_t,\forall t=1,\dots ,T\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill x_t^n& \le d_{[T]}^n{\alpha }_t^n,\forall t=1,\dots ,T\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill x_t^r& \le d_{[T]}^n{\alpha }_t^n,\forall t=1,\dots ,T\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill x_t^n,x_t^r& \ge 0,\forall t=1,\dots ,T\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\alpha }_t^n,{\alpha }_t^r& \in \{0,1\}\hfill \end{array}$$

(33)

If the model (25)-(33) has no solution, it means the model (16)-(24) has no solution too.

We have a problem with no capacity restrictions. We solve this problem and compare it with the optimal solution of the initial problem (16)-(24).

Remark 2. 

It is important to clarify that this task only makes sense if the vector w is not equal to the null vector. If the vector w is equal to the null vector it means that in each period $d_t<=c_t.$ In this case a relaxation is not possible and the problem (16)-(24) will be solved with a heuristic method (Model B).

 5. Model B (Simulation) 

 5.1. Low discrepancy sequences 

The generation of random numbers with a computer is not possible Knuth [24]. As John von Neumann said: Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin, [25]. An excellent overview of the methods of generating pseudo-random numbers are available eg. [24,26,27].

In this section we explain the number-theoretical concept of discrepancy. Then, we introduce the Halton sequence which is probably the easiest low-discrepancy 2 number generation method to describe.

Definition 1. 

Let$z_1,\cdots ,z_N$a sequence of real numbers with$0<z_i<1,i=1,\cdots ,N$. The discrepancy  $D_N$ for the sequence is defined as

$$D_N=\underset{l}{sup}|S_N(l)-N|l||$$

(34)

where l is any subinterval $[a,b)\subseteq [0,1],$$\left|l\right|=b-a,$ and $S_N\left(l\right)$ denotes the number of elements of the sequence, that belongs to the interval $l.$

A measure for how a sequence of real numbers $\{z_1,\cdots ,z_N\}$ , $a<z_i<b,i=1,\cdots ,N$ is equidistributed on an interval $[a,b]$ is the discrepancy $D_N$. Low-discrepancy sequences, also known as quasirandom sequences, are numbers that are better equidistributed in a given volume than pseudo-random numbers.

Remark 3. 

A sequence $\{z_1,\cdots ,z_N\},0<z_i<1,i=1,\cdots ,N$ of real numbers is said equidistributed on the interval $[0,1]$ if $D_N=o\left(N\right),N\to \infty ,$ [18].

The ( $\mathbb{QRN}$ ) of Halton, Sobol and Niederreiter have a low discrepancy $D_N=O(ln\left(N\right)/N)$. While pseudorandon sequences have a discrepancy $D_N=O(1/\sqrt{N}),$ [26]. Figure uses two-dimensional projection of a pseudorandom sequence and of a low-discrepancy (Halton) sequence to demonstrate the fundamental difference between the two classes of sequences.

The desirable properties of a sequence of this ( $\mathbb{QRN}$ ) may be summarized as follows ( [20]):

• the least period length should be sufficiently large,
• it should have littie intrinsic structure (such as lattice structure),
• it should have good statistical properties,
• the algorithm generating the sequence should be reasonably efficient.

It’s easy to generate sequences of Halton with the following algorithm 1, [20].

Algorithm 1: Construction of Halton sequences

The following Halton sequences of Table 5 are constructed according to algorithm 1 that uses a prime number as its base.

Table 5. Halton Numbers.

Table 5. Halton Numbers.

Remark 4. 

To generate the n-th Halton point in a sequence consider the base $b-$ary expansion of a $n={\sum }_{i=0}^{\infty }{a}_i{b}^i$ where the $b-$ary coefficients $a_i\in \{0,\cdots ,b-1\}.$ Then the n-th Halton point is $H\left(n\right)={\sum }_{i=0}^{\infty }{a}_i{b}^{-i-1}.$ It’s easy to build a Halton-sequence with the following observation: If $a_0<b-1$ then $H(n+1)=H\left(n\right)+1/b$ else if $a_0=b-1$ then $H(n+1)=H\left(n\right)-(1-b^k-b^{k+1})$ where $k=mini\ge 0:a_i\ne b-1$ (details see  [21]). This method is very efficient and will be used in this paper.

 5.2. Notation 

We use the following notation $\forall t=1,\cdots ,T$

The notation $x\in {[a,b]}_p$ means $x=(b-a)z+a,a\le b,0\le z\le 1$ and the number z is simulated with the following distribution

$$\begin{array}{cc}\hfill z=& \left(\begin{array}{ccc}{z}_h& 0& 1\\ p& q& q\end{array}\right)\hfill \\ \hfill 0& \le p\le 1,q=\frac{1}{2}(1-p),z_h\hfill & \hfill isaHaltonnumber\end{array}$$

(35)

We generate a pseudorandom number $0\le g\le 1$ to decide, that values take z, see algorithm 2.

Algorithm 2: Simulation of z with distribution (35)

In Model B (Section 5) we will investigate the class of problems with the condition

$$\begin{array}{cc}\hfill d_t^n+d_t^r=d_t& \le c_t,\forall t\hfill \end{array}$$

(36)

If this condition is not satisfied, the problem is easily solved with Model A (Section 4).

5.3 Simulation 

The simulation is based on the following lemma.

Lemma 1. 

Let $x_t=x_t^n+x_t^r$ and $d_t\le c_t,$ $\forall t=1,\cdots ,T$ Then

(37)

Proof. The proof proceeds by induction on t. In fact, if $t=1$ because x_[0] = 0 and $d_1\le c_1$, we can choose $x_1$ such that $d_1\le x_1\le c_1$.

$$\begin{array}{cc}\hfill d_{[t+1]}-x_{\left[t\right]}& =d_{t+1}-x_t+d_{\left[t\right]}-x_{[t-1]}\hfill \end{array}$$

By induction hypothesis

$$\begin{array}{c}\le d_{t+1}-x_t+x_t\hfill \end{array}$$

$$\begin{array}{cc}& \le c_{t+1}\hfill \end{array}$$

Then we can $x_{t+1}$ choose such that $d_{[t+1]}-x_{\left[t\right]}\le x_{t+1}\le c_{t+1}.$    □

Remark 5. 

If in lemma 1 $d_{\left[t\right]}-x_{[t-1]}<0$ then production in period t is $x_t=0$ because production up to period $t-1$ satisfies demand up to period t.

Remark 6. 

Lemma 1 gives the following lower bound for the production of the products

(38)

(39)

 5.4. Basis of the simulation 

The production plan is created step by step starting from period $t=1$. To determine the production in period t, we know the selected production until period $t-1$. In each period, using the constraints of the task (16)-(24) for the production of the products, we determine a lower and an upper bound. Then the production quantity is chosen randomly between the lower and upper limits. These production quantities affect the lower and upper bounds of the future period. We continue in this way until the period $t=T$. Then we calculate the value of the cost function( i.e. objective function). We repeat this procedure (N times) and choose the production plan with the lowest cost. The advantages of this method, we do not have to worry about the inventory, production or setup costs. The method is now presented in more detail.

Proposition 1. 

If $u_t^n=0$ then $x_t^n=0,{\alpha }_t^n=0,$ else $x_t^n\ge u_t^n,{\alpha }_t^n=1.$

Proof. From (17) we get

(40)

If $(x_{[t-1]}^n-d_{[t]}^n)\ge 0$ then total production of new products up to period $(t-1)$ satisfies total demand up to period t and therefore nothing is produced in period $t,$ i.e. $x_t^n=0,{\alpha }_t^n=0$.

If $(x_{[t-1]}^n-d_{[t]}^n)<0$ then the production $x_t^n$ has a lower bound. From (40) we obtain

(41)

□

Proposition 2. 

If $u_t^r=0$ then $x_t^r=0,{\alpha }_t^r=0,$ else ${\alpha }_t^r=1$ and

(42)

.

Proof. From (18) and (19)

(43)

. If $(x_{[t-1]}^r-d_{[t]}^r)\ge 0$, then total production up to period $(t-1)$ satisfies total demand up to period t and therefore nothing is produced in period $t,$ i.e. $x_t^r=0,{\alpha }_t^r=0$.

If $(x_{[t-1]}^r-d_{[t]}^r)<0$ then the production $x_t^r$ has an upper and lower bound. From (43) results the assertion.    □

Proposition 3. 

If $u_t^n>0$ and $u_t^r>0,$ then

$$\begin{array}{cc}\hfill o_t^n& =min\left\{c_t-u_t^r,d_{\left[T\right]}^n-x_{[t-1]}^n\right\}\hfill \\ \hfill o_t^r& =min\left\{c_t-x_t^n,d_{\left[T\right]}^r-x_{[t-1]}^r{r}_{\left[t\right]}^{re}-x_{[t-1]}^r\right\}\hfill \end{array}$$

And

$$\begin{array}{cc}\hfill u_t^n& \le x_t^n\le o_t^n\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill u_t^r& \le x_t^r\le o_t^r\hfill \end{array}$$

(45)

Proof. From (20) and remark 6

$$\begin{array}{cc}\hfill c_t& \ge x_t^n+x_t^r\ge x_t^n+u_t^r\hfill \\ \hfill c_t-u_t^r& \ge x_t^n\hfill \end{array}$$

(46)

From (2)

(47)

From (46) and (47)

$$\begin{array}{c}{x}_t^n\le o_t^n\hfill \end{array}$$

(48)

xntont Therefore we simulate the production $x_t^n$ according to (35)

(49)

From (20) using $x_t^n$ of (49) we get

(50)

and then from (2) we obtain

(51)

From (50), (51) and (42) results

(52)

Therefore we simulate the production $x_t^r$ according to (35)

(53)

   □

 5.5. Simulation of the objective function 

Let R be a matrix with Halton’s $\mathbb{QRN}$

$$R=\left(\begin{array}{ccc}R(1,1)& \cdots & R(1,T)\\ R(2,1)& \cdots & R(2,T)\end{array}\right)=\left(\begin{array}{c}{R}_1\\ R_2\end{array}\right)$$

We simulate the production starting in period t=1. The calculation of the objective function is carried out with algorithm 3 using proposition (1), (2) and (3).

Algorithm 3: Calculation of the objective function

Further information on the complexity-theoretic approach to randomness can be found in  [28,29,30].

Remark 7. 

1. The generation of the matrix R[N × T] requires O(N × T) operations.

2. Die evaluation of the function $\phi (x,y)$ with algorithm 3 requires $O\left(T\right)$ operations Then, the computational complexity of the simulation with N points and T periods is $O(N\times T).$

 6. Numerical Experiments 

 6.1. Test Design 

We analyse the quality of our solution approach of model A and model B by defining 11 problem classes (PC) by varying the number of periods, see Table 6. The planning horizon T is made very large because in the paper industry planning is done daily. Each PC consists of 200 test instances (TI). In model A, 1824 of the 2200 TI were solvable, in model B 2200. In total, we examined 4024 TI. Model A and Model B are implemented on a computer with $Intel\left(R\right)$$Core\left(TM\right)iT-9700K$, $CPU@3.60GHz$,$3600MHz$.

Table 6. Problem classes.

Table 6. Problem classes.

We vary different parameters to define the TI, e.g., the time between orders (TBO) to determine setup costs. The specifications of the parameters are designed in an exaggerated form, which may not occur in practice, to make the TI as difficult as possible.

The parameters for generating data sets (see Table 7) use the following notation $x\in [a;b]\iff x=(b-a)\theta +a,a\le b,$ where $0\le \theta \le 1$ is a random number, that means the values are uniformly distributed on the interval $[a,b].$

 6.2. Model A 

We randomly generate capacities $c_t$, then with condition (1) randomly generate aggregate demand $d_t$. This latter is then randomly cut into $d_t^n$ and $d_t^r$. The remaining parameters according to Table 7.

Table 7. Original Data.

Table 7. Original Data.

6.2.1. Results of Model A

The solution of problem (25)-(33) (TD) and the optimal solution of problem (16)-(24) (OP) were found with the Gurobi solver version 9.0.3.

We randomly generated 200 instances per PC and usually between 10 and approx. 20% of the instances have no solution (see line Count). This justifies the fact that the problem is not standard. For example, in Table A1 we see that of the 200 instances for problem class $T=60$, only 167 had a solution.

With the following notation, we can better understand the results of Table A1.

In Table A1 the relative error for Total costs (Tc) and CPU-time for every problem class was calculated as

$$\begin{array}{cc}\hfill RelativeAvg.Error\left(Tc\right)& =\frac{{\mu }_T-{\mu }_T^{*}}{{\mu }_T^{*}}\hfill \\ \hfill RelativeAvg.Error\left(CPU\right)& =\frac{CP{U}_T-CP{U}_T^{*}}{CP{U}_T^{*}}\hfill \end{array}$$

This is exactly how we calculated the relative errors in inventory cost and setup cost.

Attached in Appendix A are the results of the average total cost, average CPU time, average Inventory costs

There is hardly any difference between the cost of relaxation (TD) and the original task (OP). Only the computation time for relaxation is faster. The longer the planning horizon, the smaller the difference between the optimal solutions of the problems TD and original optimization task OP. This feature applies to the stocks of the return and setup costs too. On the other hand, the inventory costs for problem TD are always smaller than the inventory costs of problems OP. For more details, please see Figure A1 and Appendix A.1 for the exact calculations. However, we see beyond doubt that this class of problems can be solved very well either with the relaxation TD or directly with a standard solver (here Gurobi).

 6.3. Model B 

[13] have shown that several families of CLSP are Np-hard. For the construction of the Np-hard instances (2200 instances) we follow the findings of  [14]. They use the following notation $Nr/\alpha /\beta /\gamma /\sigma$, where $Nr,\alpha ,\beta ,\gamma$ and $\sigma$ specify respectively the number of items, a special structure for the setup costs, the holding costs, production costs, and capacities. In this paper  [14] they show that the following class 2/C/G/A/C is NP-hard. For this reason we have created 2200 instances, where the set-up costs per product and capacities per instance are constant. The holding cost do not necessarily follow a specified pattern, the production costs can be chosen arbitrarily. The maximum calculation time for Gurobi is 600 seconds. The heuristic operates according to the number of simulations, which gradually increases with the number of periods. The largest group T300 uses 130 seconds. The parameters for generating data sets (see Table 8) use the following notation $x\in [a;b]\iff x=(b-a)\theta +a,a\le b,$ where $0\le \theta \le 1$ is a random number, that means the values are uniformly distributed on the interval $[a,b].$

The important assumption in model B is: the capacities for each TI is constant, the set-up costs in each TI and for each product are constant. These parameters vary between a minimum and a maximum depending on the T parameter.

The capacities for each TI vary according to the parameter T. For example if $T=15$ the capacities are between 600 and 800. If $T=300$ the capacities vary between 3000 and 5500. Analogously the other parameters. The remaining parameters according to Table 8.

Table 8. Model B: Parameters.

Table 8. Model B: Parameters.

For the simulation we used following parameters:

$$N_T=2^{15}T$$

$N_T$ is the Halton’s numbers used for T periods (see Table 9).

Table 9. Periods, Halton’s numbers.

Table 9. Periods, Halton’s numbers.

All instances for $T=15,30\cdots ,300$ have the same schema Step 1 until Step 5. We used $\tau =4$ (see Table 10).

Table 10. Schematic of the simulation.

Table 10. Schematic of the simulation.

Algorithm 4: Heuristic: blind search

6.3.1. Results of Model B

We generated 200 random instances for each problem class (PC) and with a Box-plot we compared the feasible solutions $G_{Ti}(x^n,x^r)$ of the problem (16)-(24) found by Gurobi 9.0.3 with the solution $S_{Ti}(x^n,x^r)$ of the simulation presented in this paper and clearly see the similarity of the results found (see Figure A4, A5).

With the following notation we present the results (see Table A6).

In Table A6 the relative error for Total costs (Tc) and CPU-time for every problem class was calculated as

$$\begin{array}{cc}\hfill RelativeAvg.Error\left(Tc\right)& =\frac{{\mu }_T-{\mu }_T^{*}}{{\mu }_T^{*}}\hfill \\ \hfill RelativeAvg.Error\left(CPU\right)& =\frac{CPU\_S_T-CP{U}_T^{*}}{CP{U}_T^{*}}\hfill \end{array}$$

The graphical comparison is shown in Figure A3.

This is exactly how we calculated the relative errors in inventory cost and setup cost. Attached in Appendix B are the results of the average total cost, average CPU time, average Inventory costs. Figure A3 (average CPU time) clearly shows that the Gurobi admissible solution for the PC from T150 to T300 are not optimal and that the CPU of the simulation is much faster.

The simulation in average determines the setup costs and return inventory costs always higher than Gurobi. On the other hand, the inventory costs of Gurobi are higher than the simulation. What can we say about the quality of the solution of the problems? the simulation could not give a better solution than Gurobi’s solution. Gurobi solved the (PC) problems up to T120 in an optimal way. The problems from T150 to T300 were not solved optimally by Gurobi, since it would take too much time due to the Np-hard category of the problem. The simulation found feasible solutions much faster than Gurobi. Here is the advantage of the simulation, the simplicity of its implementation and the speed in finding an acceptable solution.

 7. Conclusions and Outlook 

We have analyzed a problem that belongs to the NP-hard class. However, the choice of the parameters is very important to obtain a problem that is really Np-hard. With the choice of parameters made in model A, we see that this class of problems is easily solved with a standard solver. By doing a relaxation of the problem, the solution is found more quickly. The error rate is between 0.02% and 2% (taking into account more than 1800 instances).

If we choose the parameters according to model B, Gurobi needs a lot of time to find the optimal solution. In this kind of problem the presented simulation can help a lot in finding a good solution. We have seen that the error rate is between 1.7% and 3.5% (taking into account more than 2000 instances). The simulation is sometimes better than Gurobi, but on average Gurobi solves the problem with a maximum time of 600 seconds better than the simulation. The great advantage of the simulation is that the calculation is extremely fast and easy.Thanks to the Halton numbers, few simulations are needed to obtain a very good approximation of the solution.

The simulation has a general character and can be adapted to investigate even more complex problems of the Np-hard category ( for example, researching CLSP problems with n products).

The quality of the solutions can be improved by increasing the number of simulations but it is necessary to have a fairly fast computer. In this work we use at most 10 million simulations. Another parameter that influences the quality of the solutions is the correct choice of the probability p (see (35)). Is there an optimal probability? This is a question for further research.