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.

  $${\displaystyle \sum _{j=1}^t}({\displaystyle d_j^n}+{\displaystyle d_j^r})\le {\displaystyle \sum _{j=1}^t}{c}_j,{\displaystyle \forall t=1,\cdots ,T}$$

  (1)
• Initial and end inventory stocks are zero, (We can always transform a problem with non zero initial or final stock by adapting the demand), i.e.

  $$\begin{array}{ccc}\hfill {\displaystyle i_0^n}=0,& {\displaystyle i_0^r}=0,\hfill & \hfill {\displaystyle i_0^{re}}=0\\ \hfill {\displaystyle i_T^n}=0,& {\displaystyle i_T^r}=0\hfill & \hfill \end{array}$$

  (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.

  $$\sum _{j=1}^t{\displaystyle d_j^r}\le \sum _{j=1}^t{\displaystyle r_j^{re}}{\displaystyle \forall t=1,\cdots ,T}.$$

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

  $$\sum _{j=t}^T{\displaystyle h_j^{re}}\le \sum _{j=t}^T{\displaystyle h_j^r},{\displaystyle \forall t=1,\cdots ,T}$$

  (4)
  This hypothesis can be found in [22] too.

Hence, the problem can be formulated as

$$\begin{array}{cc}\hfill f(x^n,x^r)& ={\displaystyle \sum _{t=1}^T}({\displaystyle s_t^n}{\alpha }_t+{\displaystyle p_t^n}{x}_t^n+{\displaystyle h_t^n}{\displaystyle i_t^n})\hfill \\ & +{\displaystyle \sum _{t=1}^T}({\displaystyle s_t^r}{\alpha }_t^r+{\displaystyle p_t^r}{x}_t^r+{\displaystyle h_t^r}{\displaystyle i_t^r})\hfill \\ & +{\displaystyle \sum _{t=1}^T}\left({\displaystyle h_t^{re}}{\displaystyle i_t^{re}}\right)⟶min\hfill \end{array}$$

(5)

subject to

$$\begin{array}{cc}\hfill 1-1{\displaystyle i_t^n}& ={\displaystyle i_{t-1}^n}+x_t^n-{\displaystyle d_t^n},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\displaystyle i_t^r}& ={\displaystyle i_{t-1}^r}+x_t^r-{\displaystyle d_t^r},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\displaystyle i_t^{re}}& ={\displaystyle i_{t-1}^{re}}+{\displaystyle r_t^{re}}-x_t^r,{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(8)

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

(9)

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

(10)

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

(11)

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

(12)

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

(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)–(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

$$\begin{array}{ccccc}\hfill {\displaystyle i_t^n}& =\sum _{j=1}^t\left(x_j^n-{\displaystyle d_j^n}\right),\hfill & \hfill {\displaystyle i_t^r}& =\sum _{j=1}^t\left(x_j^r-{\displaystyle d_j^r}\right),\hfill & \hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\displaystyle i_t^{re}}& =\sum _{j=1}^t\left({\displaystyle r_j^{re}}-x_j^r\right)\hfill \end{array}$$

(15)

And with a bit of algebraic transformations, we obtain

$$\begin{array}{c}\hfill \sum _{t=1}^T({\displaystyle p_t^n}{x}_t^n+{\displaystyle h_t^n}{\displaystyle i_t^n})={\displaystyle \sum _{t=1}^T}(p_t^n+{\displaystyle \sum _{j=t}^T}{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({\displaystyle p_t^r}{x}_t^r+{\displaystyle h_t^r}{\displaystyle i_t^r}+h_t^{re}{\displaystyle i_t^{re}})={\displaystyle \sum _{t=1}^T}({\displaystyle p_t^r}+{\displaystyle \sum _{j=t}^T}(h_j^r-h_j^{re}))x_t^r\\ \hfill +{\displaystyle \sum _{t=1}^T}(-h_t^r{\displaystyle \sum _{k=1}^t}{d}_k^r+h_t^{re}{\displaystyle \sum _{k=1}^t}{r}_k^{re})\end{array}$$

With the notation

$$\begin{array}{cc}\hfill & v_t^n=p_t^n+{\displaystyle \sum _{j=t}^T}{h}_j^n,\hfill \\ \hfill & v_t^r=p_t^r+{\displaystyle \sum _{j=t}^T}(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

$$\begin{array}{cc}\hfill \phi (x^n,x^r)& ={\displaystyle \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⟶min\hfill \end{array}$$

(16)

Subject to

$$\begin{array}{cc}\hfill 1-1{\displaystyle \sum _{j=1}^t}{x}_j^n& \ge {\displaystyle \sum _{j=1}^t}{\displaystyle d_t^n},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{j=1}^t}{x}_j^r& \ge {\displaystyle \sum _{j=1}^t}{\displaystyle d_t^r},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{j=1}^t}{\displaystyle r_j^{re}}& \ge {\displaystyle \sum _{j=1}^t}{x}_j^r,{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(19)

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

(20)

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

(21)

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

(22)

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

(23)

$$\begin{array}{cc}\hfill {\alpha }_t^n,{\alpha }_t^r& \in \left\{0,1\right\}\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\left\{d_{t+1}-c_{t+1}+{\tilde{w}}_{t+1};0\right\},t=T-1,\cdots ,1.\hfill \end{array}$$

We define

$$\begin{array}{cc}\hfill w_1& ={\tilde{w}}_1\hfill \end{array}$$

$$\begin{array}{cc}\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.

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.

$$\begin{array}{cc}\hfill \phi (x,y)& ={\displaystyle \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⟶min\hfill \end{array}$$

(25)

Subject to

$$\begin{array}{cc}\hfill 1-1{\displaystyle \sum _{j=1}^t}{x}_j^n& \ge {\displaystyle \sum _{j=1}^t}{\displaystyle d_t^n},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{j=1}^t}{x}_j^r& \ge {\displaystyle \sum _{j=1}^t}{\displaystyle d_t^r},{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{j=1}^t}{\displaystyle r_j^{re}}& \ge {\displaystyle \sum _{j=1}^t}{x}_j^r,{\displaystyle \forall t=1,\cdots ,T}\hfill \end{array}$$

(28)

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

(29)

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

(30)

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

(31)

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

(32)

$$\begin{array}{cc}\hfill {\alpha }_t^n,{\alpha }_t^r& \in \left\{0,1\right\}\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).

Table 3. Original Data.

t	${\mathit{r}}_{\mathit{t}}^{\mathbf{re}}$	${\mathit{d}}_{\mathit{t}}^{\mathit{n}}$	${\mathit{d}}_{\mathit{t}}^{\mathit{r}}$	${\mathit{d}}_{\mathit{t}}$	${\mathit{c}}_{\mathit{t}}$	w
1	198	153	183	336	609	57
2	806	84	302	386	632	246
3	223	100	146	246	101	-145
4	283	100	127	227	295	68
5	500	248	598	846	620	-226
6	500	0	0	0	561	0
Sum	2510	685	1356	2041	2818	0

Table 3. Original Data.

t	${\mathit{r}}_{\mathit{t}}^{\mathbf{re}}$	${\mathit{d}}_{\mathit{t}}^{\mathit{n}}$	${\mathit{d}}_{\mathit{t}}^{\mathit{r}}$	${\mathit{d}}_{\mathit{t}}$	${\mathit{c}}_{\mathit{t}}$	w
1	198	153	183	336	609	57
2	806	84	302	386	632	246
3	223	100	146	246	101	-145
4	283	100	127	227	295	68
5	500	248	598	846	620	-226
6	500	0	0	0	561	0
Sum	2510	685	1356	2041	2818	0

Table 4. Demand Transformation.

t	$\widetilde{{\mathit{d}}_{\mathit{t}}^{\mathit{n}}}$	$\widetilde{{\mathit{d}}_{\mathit{t}}^{\mathit{r}}}$	${\tilde{\mathit{d}}}_{\mathit{t}}$	${\mathit{c}}_{\mathit{t}}$	${\tilde{\mathit{d}}}_{\left[\mathit{t}\right]}$	${\mathit{c}}_{\left[\mathit{t}\right]}$
1	195	198	393	609	393	609
2	42	590	632	632	1025	1241
3	101	0	101	101	1126	1342
4	295	0	295	295	1421	1637
5	52	568	620	620	2041	2257
6	0	0	0	561	2041	2818
Sum	685	1356	2041	2818

Table 4. Demand Transformation.

t	$\widetilde{{\mathit{d}}_{\mathit{t}}^{\mathit{n}}}$	$\widetilde{{\mathit{d}}_{\mathit{t}}^{\mathit{r}}}$	${\tilde{\mathit{d}}}_{\mathit{t}}$	${\mathit{c}}_{\mathit{t}}$	${\tilde{\mathit{d}}}_{\left[\mathit{t}\right]}$	${\mathit{c}}_{\left[\mathit{t}\right]}$
1	195	198	393	609	393	609
2	42	590	632	632	1025	1241
3	101	0	101	101	1126	1342
4	295	0	295	295	1421	1637
5	52	568	620	620	2041	2257
6	0	0	0	561	2041	2818
Sum	685	1356	2041	2818

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] or [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 (low discrepancy sequences are called quasi-random sequences) number generation method to describe.

Definition 1. 

Let $\left\{z_1,\cdots ,z_N\right\}$ 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}\left|S_N\left(l\right)-N\left|l\right|\right|$$

(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 $\left\{z_1,\cdots ,z_N\right\}$ , $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 $\left\{z_1,\cdots ,z_N\right\},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 2 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.

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 \left\{0,\cdots ,b-1\right\}.$ 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=min\left\{i\ge 0:a_i\ne b-1\right\}$ (details see [21]). This method is very efficient and will be used in this paper.

5.2. Notation

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

Name	Meaning
${\displaystyle d_{\left[t\right]}^n}$ =	${\displaystyle \sum _{j=1}^t}{d}_j^n$	${\displaystyle d_{\left[0\right]}^n}$ =	0
${\displaystyle d_{\left[t\right]}^r}$ =	${\displaystyle \sum _{j=1}^t}{d}_j^r$	${\displaystyle d_{\left[0\right]}^r}$ =	0
${\displaystyle r_{\left[t\right]}^{re}}$=	${\displaystyle \sum _{j=1}^t}{r}_j^{re}$	${\displaystyle r_{\left[0\right]}^{re}}$=	0
${\displaystyle x_{\left[t\right]}^n}$ =	${\displaystyle \sum _{j=1}^t}{x}_j^n$	${\displaystyle x_{\left[0\right]}^n}$ =	0
${\displaystyle x_{\left[t\right]}^r}$ =	${\displaystyle \sum _{j=1}^t}{x}_j^r$	${\displaystyle x_{\left[0\right]}^r}$ =	0

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,{\displaystyle \forall t=1,\cdots ,T}.$ Then

$${\displaystyle d_{\left[t\right]}}-{\displaystyle x_{\left[t-1\right]}}\le x_t\le c_t,{\displaystyle \forall t=1,\cdots ,T}$$

(37)

Proof. 

The proof proceeds by induction on t. In fact, if $t=1$ because ${\displaystyle x_{\left[0\right]}}=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 {\displaystyle d_{\left[t+1\right]}}-{\displaystyle x_{\left[t\right]}}& =d_{t+1}-x_t+{\displaystyle d_{\left[t\right]}}-{\displaystyle x_{\left[t-1\right]}}\hfill \end{array}$$

By induction hypothesis

$$\begin{array}{cc}\hfill 1-1& \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 ${\displaystyle d_{\left[t+1\right]}}-{\displaystyle x_{\left[t\right]}}\le x_{t+1}\le c_{t+1}.$    □

Remark 5. 

If in Lemma 1 ${\displaystyle d_{\left[t\right]}}-{\displaystyle x_{\left[t-1\right]}}<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

$$\begin{array}{cc}\hfill x_t^n\ge u_t^n=& max\left\{{\displaystyle d_{\left[t\right]}^n}-{\displaystyle x_{\left[t-1\right]}^n};0\right\}\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill x_t^r\ge u_t^r=& max\left\{{\displaystyle d_{\left[t\right]}^r}-{\displaystyle x_{\left[t-1\right]}^r};0\right\}\hfill \end{array}$$

(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

$${\displaystyle x_{\left[t\right]}^n}-{\displaystyle d_{\left[t\right]}^n}=x_t+({\displaystyle x_{\left[t-1\right]}^n}-{\displaystyle d_{\left[t\right]}^n})\ge 0$$

(40)

If $({\displaystyle x_{\left[t-1\right]}^n}-{\displaystyle d_{\left[t\right]}^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 $({\displaystyle x_{\left[t-1\right]}^n}-{\displaystyle d_{\left[t\right]}^n})<0,$ then the production $x_t^n$ has a lower bound. From (40) we obtain

$$x_t^n\ge {\displaystyle d_{\left[t\right]}^n}-{\displaystyle x_{\left[t-1\right]}^n}.$$

(41)

   □

Proposition 2. 

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

$${\displaystyle r_{\left[t\right]}^{re}}-{\displaystyle x_{\left[t-1\right]}^r}\ge x_t^r\ge u_t^r.$$

(42)

Proof. 

From (18) and (19)

$${\displaystyle r_{\left[t\right]}^{re}}\ge x_t^r+{\displaystyle x_{\left[t-1\right]}^r}\ge {\displaystyle d_{\left[t\right]}^r}.$$

(43)

If $({\displaystyle x_{\left[t-1\right]}^r}-{\displaystyle d_{\left[t\right]}^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 $({\displaystyle x_{\left[t-1\right]}^r}-{\displaystyle d_{\left[t\right]}^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,{\displaystyle d_{\left[T\right]}^n}-{\displaystyle x_{\left[t-1\right]}^n}\right\}\hfill \\ \hfill o_t^r& =min\left\{c_t-x_t^n,{\displaystyle d_{\left[T\right]}^r}-{\displaystyle x_{\left[t-1\right]}^r},{\displaystyle r_{\left[t\right]}^{re}}-{\displaystyle x_{\left[t-1\right]}^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)

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

(47)

From From (46) and (47)

$$x_t^n\le o_t^n$$

(48)

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

$$x_t^n\in {[u_t^n,o_t^n]}_p,{\alpha }_t^n=1.$$

(49)

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

$$c_t-x_t^n\ge x_t^r$$

(50)

and then from (2) we obtain

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

(51)

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

$$x_t^r\le o_t^r$$

(52)

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

$$x_t^r\in {[u_t^r,o_t^r]}_p,{\alpha }_t^r=1.$$

(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. 

• The generation of the matrix $R[N\times T]$ requires $O(N\times T)$ operations.
• 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$.

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.

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.

T	$\in \left\{15,30,\cdots ,300\right\}$
m	$\in Count=\left\{172,155,\cdots ,170\right\}$
${\phi }_{Ti}(x^n,x^r)$	$i=1,\cdots ,m$ solution of problem (25)-(33)
${\phi }_{Ti}^{*}(x^n,x^r)$	$i=1,\cdots ,m$ solution of problem (16)-(24)
${\mu }_T$	$=\frac{{\sum }_{i=1}^m{\phi }_{Ti}(x^n,x^r)}{m}$
${\mu }_T^{*}$	$=\frac{{\sum }_{i=1}^m{\phi }_{Ti}^{*}(x^n,x^r)}{m}$
$CP{U}_T$	$=\frac{{\sum }_{i=1}^{m}CP{U}_{Ti}}{m}$
$CP{U}_T^{*}$	$=\frac{{\sum }_{i=1}^{m}CP{U}_{Ti}^{*}}{m}$

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.

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).

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

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 and Figure A5).

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

T	$\in \left\{15,30,\cdots ,300\right\}$
m	$=200$
$S_{Ti}(x^n,x^r)$	$i=1,\cdots ,m$
$G_{Ti}(x^n,x^r)$	$i=1,\cdots ,m$
${\mu }_T$	$=\frac{{\sum }_{i=1}^m{S}_{Ti}(x^n,x^r)}{m}$
${\mu }_T^{*}$	$=\frac{{\sum }_{i=1}^m{G}_{Ti}(x^n,x^r)}{m}$
$CPU\_S_T$	$=\frac{{\sum }_{i=1}^{m}CP{U}_{Ti}}{m}$ with Simulation
$CP{U}_T^{*}$	$=\frac{{\sum }_{i=1}^{m}CP{U}_{Ti}^{*}}{m}$ with Gurobi

In Table A3 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.

Appendix A.1. Model A

Figure A1. Model A: Average Totalcost and CPU Time.

Figure A1. Model A: Average Totalcost and CPU Time.

Figure A2. Model A: Average Inventory cost.

Figure A2. Model A: Average Inventory cost.

Appendix A.2. Model B

Visualization of the results and average costs.

Figure A3. Model B: Average Total cost and CPU time.

Figure A3. Model B: Average Total cost and CPU time.

Figure A4. Model B: Total costs.

Figure A4. Model B: Total costs.

Figure A5. Model B: Total costs.

Figure A5. Model B: Total costs.

Appendix B.1. Model A

Table A1. Model A: Average Total Costs.

Total Costs	T	15	30	60	90	120	150
	Count	172	155	167	162	166	164
TD	Mean	274061	665537	1798127	3348945	5552436	8062797
	Std.	87161	151370	389288	660684	906092	1349699
Optimal	Mean	268593	659453	1792494	3344788	5546715	8058311
	Std.	83431	148783	385859	659403	904686	1348457
Relative error	Tc	2,04%	0,92%	0,31%	0,12%	0,10%	0,06%
Total costs	T	150	180	210	240	270	300
	Count	164	164	168	168	170	168
TD	Mean	8062797	11087631	14620895	19015245	23243458	28160753
	Std.	1349699	1811439	2199306	2911556	3384673	3870830
Optimal	Mean	8058311	11081145	14614610	19008049	23237758	28154353
	Std.	1348457	1808995	2198938	2910200	3382841	3868995
Relative error	Tc	0,06%	0,06%	0,04%	0,04%	0,02%	0,02%

Table A1. Model A: Average Total Costs.

Total Costs	T	15	30	60	90	120	150
	Count	172	155	167	162	166	164
TD	Mean	274061	665537	1798127	3348945	5552436	8062797
	Std.	87161	151370	389288	660684	906092	1349699
Optimal	Mean	268593	659453	1792494	3344788	5546715	8058311
	Std.	83431	148783	385859	659403	904686	1348457
Relative error	Tc	2,04%	0,92%	0,31%	0,12%	0,10%	0,06%
Total costs	T	150	180	210	240	270	300
	Count	164	164	168	168	170	168
TD	Mean	8062797	11087631	14620895	19015245	23243458	28160753
	Std.	1349699	1811439	2199306	2911556	3384673	3870830
Optimal	Mean	8058311	11081145	14614610	19008049	23237758	28154353
	Std.	1348457	1808995	2198938	2910200	3382841	3868995
Relative error	Tc	0,06%	0,06%	0,04%	0,04%	0,02%	0,02%

Table A2. Model A: Average CPU time.

CPU Time	T	15	30	60	90	120	150
	Count	172	155	167	162	166	164
TD	Mean	0,10	0,11	0,21	0,31	0,50	0,78
	Std.	0,06	0,05	0,10	0,16	0,26	0,39
Optimal	Mean	0,14	0,28	0,41	0,59	1,12	1,79
	Std.	0,06	0,11	0,17	0,30	0,51	1,42
Relative error	CPU	-25,27%	-59,97%	-49,61%	-47,53%	-55,39%	-56,26%
CPU Time	T	150	180	210	240	270	300
	Count	164	164	168	168	170	168
TD	Mean	0,78	1,00	1,25	1,59	1,77	2,36
	Std.	0,39	0,47	0,69	0,82	0,90	1,39
Optimal	Mean	1,79	2,44	3,52	5,03	5,33	7,60
	Std.	1,42	1,73	3,00	3,50	3,98	5,97
Relative error	CPU	-56,26%	-59,00%	-64,53%	-68,46%	-66,87%	-68,94%

Table A2. Model A: Average CPU time.

CPU Time	T	15	30	60	90	120	150
	Count	172	155	167	162	166	164
TD	Mean	0,10	0,11	0,21	0,31	0,50	0,78
	Std.	0,06	0,05	0,10	0,16	0,26	0,39
Optimal	Mean	0,14	0,28	0,41	0,59	1,12	1,79
	Std.	0,06	0,11	0,17	0,30	0,51	1,42
Relative error	CPU	-25,27%	-59,97%	-49,61%	-47,53%	-55,39%	-56,26%
CPU Time	T	150	180	210	240	270	300
	Count	164	164	168	168	170	168
TD	Mean	0,78	1,00	1,25	1,59	1,77	2,36
	Std.	0,39	0,47	0,69	0,82	0,90	1,39
Optimal	Mean	1,79	2,44	3,52	5,03	5,33	7,60
	Std.	1,42	1,73	3,00	3,50	3,98	5,97
Relative error	CPU	-56,26%	-59,00%	-64,53%	-68,46%	-66,87%	-68,94%

Appendix B.2. Model B

Table A3. Model B: Average Total Costs.

Total Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	185971	342933	815220	4576352	7829443	3979678
	Std.	21273	45909	155884	858442	1559305	397507
Gurobi	Mean	182040	332163	784654	4498331	7718552	3845522
	Std.	21125	45150	155097	882996	1594738	393440
Relative error	Tc	2,16%	3,24%	3,90%	1,73%	1,44%	3,49%
Total costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	3979678	5477237	8177244	19065355	18786745	20708665
	Std.	397507	593444	1275892	4030464	3201422	3672141
Gurobi	Mean	3845522	5318346	7905438	18570461	18146172	19992493
	Std.	393440	590176	1265172	4055425	3179322	3670651
Relative error	Tc	3,49%	2,99%	3,44%	2,66%	3,53%	3,58%

Table A3. Model B: Average Total Costs.

Total Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	185971	342933	815220	4576352	7829443	3979678
	Std.	21273	45909	155884	858442	1559305	397507
Gurobi	Mean	182040	332163	784654	4498331	7718552	3845522
	Std.	21125	45150	155097	882996	1594738	393440
Relative error	Tc	2,16%	3,24%	3,90%	1,73%	1,44%	3,49%
Total costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	3979678	5477237	8177244	19065355	18786745	20708665
	Std.	397507	593444	1275892	4030464	3201422	3672141
Gurobi	Mean	3845522	5318346	7905438	18570461	18146172	19992493
	Std.	393440	590176	1265172	4055425	3179322	3670651
Relative error	Tc	3,49%	2,99%	3,44%	2,66%	3,53%	3,58%

Table A4. Model B: Average CPU time.

CPU Time	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	0,79	2,28	5,94	11,68	19,87	34,12
	Std.	0,13	0,13	0,16	0,18	0,23	0,54
Gurobi	Mean	0,12	1,57	90,38	111,93	204,94	597,03
	Std.	0,08	1,53	156,10	186,04	245,47	36,02
Relative error	CPU	552,92%	44,87%	-93,43%	-89,56%	-90,31%	-94,28%
CPU Time	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	34,12	46,26	62,84	77,96	101,69	130,17
	Std.	0,54	0,58	1,22	2,63	2,74	3,78
Gurobi	Mean	597,03	596,88	597,47	591,51	596,54	593,22
	Std.	36,02	4,06	3,92	33,19	3,80	7,14
Relative error	CPU	-94,28%	-92,25%	-89,48%	-86,82%	-82,95%	-78,06%

Table A4. Model B: Average CPU time.

CPU Time	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	0,79	2,28	5,94	11,68	19,87	34,12
	Std.	0,13	0,13	0,16	0,18	0,23	0,54
Gurobi	Mean	0,12	1,57	90,38	111,93	204,94	597,03
	Std.	0,08	1,53	156,10	186,04	245,47	36,02
Relative error	CPU	552,92%	44,87%	-93,43%	-89,56%	-90,31%	-94,28%
CPU Time	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	34,12	46,26	62,84	77,96	101,69	130,17
	Std.	0,54	0,58	1,22	2,63	2,74	3,78
Gurobi	Mean	597,03	596,88	597,47	591,51	596,54	593,22
	Std.	36,02	4,06	3,92	33,19	3,80	7,14
Relative error	CPU	-94,28%	-92,25%	-89,48%	-86,82%	-82,95%	-78,06%

Table A5. Model B: Inventory costs.

Inventory Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	28307	52842	112536	253095	340599	260657
	Std.	5341	8170	12391	86118	107587	28744
Gurobi	Mean	27801	57117	132943	509248	790970	318360
	Std.	6907	12243	29889	194930	324335	43350
Relative error	Inv.	1,82%	-7,49%	-15,35%	-50,30%	-56,94%	-18,13%
Inventory costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	260657	345666	542344	990715	1147482	1282444
	Std.	28744	38613	75292	175805	161861	196713
Gurobi	Mean	318360	507104	835915	1859977	1730800	1959659
	Std.	43350	89630	191367	486708	372666	455506
Relative error	Inv.	-18,13%	-31,84%	-35,12%	-46,74%	-33,70%	-34,56%

Table A5. Model B: Inventory costs.

Inventory Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	28307	52842	112536	253095	340599	260657
	Std.	5341	8170	12391	86118	107587	28744
Gurobi	Mean	27801	57117	132943	509248	790970	318360
	Std.	6907	12243	29889	194930	324335	43350
Relative error	Inv.	1,82%	-7,49%	-15,35%	-50,30%	-56,94%	-18,13%
Inventory costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	260657	345666	542344	990715	1147482	1282444
	Std.	28744	38613	75292	175805	161861	196713
Gurobi	Mean	318360	507104	835915	1859977	1730800	1959659
	Std.	43350	89630	191367	486708	372666	455506
Relative error	Inv.	-18,13%	-31,84%	-35,12%	-46,74%	-33,70%	-34,56%

Table A6. Model B: Return stock cost.

Return Inventory Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	45867	92667	249888	3690209	6646915	821946
	Std.	18040	35850	137395	1028216	1793396	208315
Gurobi	Mean	41356	74641	195294	3346068	6070949	700634
	Std.	17493	33456	133042	1109603	1957840	207573
Relative error	Ret. Inv.	10,91%	24,15%	27,96%	10,28%	9,49%	17,31%
Return Inventory costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	821946	1548362	2628642	9347966	7251886	8608929
	Std.	208315	482657	1041268	4406311	2744151	3706917
Gurobi	Mean	700634	1293559	2169556	8142232	6309891	7524531
	Std.	207573	477084	1026872	4293352	2685984	3610684
Relative error	Ret. Inv.	17,31%	19,70%	21,16%	14,81%	14,93%	14,41%

Table A6. Model B: Return stock cost.

Return Inventory Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	45867	92667	249888	3690209	6646915	821946
	Std.	18040	35850	137395	1028216	1793396	208315
Gurobi	Mean	41356	74641	195294	3346068	6070949	700634
	Std.	17493	33456	133042	1109603	1957840	207573
Relative error	Ret. Inv.	10,91%	24,15%	27,96%	10,28%	9,49%	17,31%
Return Inventory costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	821946	1548362	2628642	9347966	7251886	8608929
	Std.	208315	482657	1041268	4406311	2744151	3706917
Gurobi	Mean	700634	1293559	2169556	8142232	6309891	7524531
	Std.	207573	477084	1026872	4293352	2685984	3610684
Relative error	Ret. Inv.	17,31%	19,70%	21,16%	14,81%	14,93%	14,41%

Table A7. Model B: Average Setup costs.

Setup Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	50033	106876	238570	293071	388685	1521814
	Std.	10663	17638	31424	27563	34664	212705
Gurobi	Mean	55707	114588	248504	307481	408104	1458747
	Std.	10343	17268	30443	27399	37230	201433
Relative error	Setup	-10,19%	-6,73%	-4,00%	-4,69%	-4,76%	4,32%
Setup Costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	1521814	1699197	2214098	2828660	3716988	4920930
	Std.	212705	233148	330991	404730	467555	902775
Gurobi	Mean	1458747	1641182	2118141	2669355	3397775	4526298
	Std.	201433	223689	310292	367954	412785	200
Relative error	Setup	4,32%	3,53%	4,53%	5,97%	9,39%	8,72%

Table A7. Model B: Average Setup costs.

Setup Costs	T	15	30	60	90	120	150
	Count	200	200	200	200	200	200
Simulation	Mean	50033	106876	238570	293071	388685	1521814
	Std.	10663	17638	31424	27563	34664	212705
Gurobi	Mean	55707	114588	248504	307481	408104	1458747
	Std.	10343	17268	30443	27399	37230	201433
Relative error	Setup	-10,19%	-6,73%	-4,00%	-4,69%	-4,76%	4,32%
Setup Costs	T	150	180	210	240	270	300
	Count	200	200	200	200	200	200
Simulation	Mean	1521814	1699197	2214098	2828660	3716988	4920930
	Std.	212705	233148	330991	404730	467555	902775
Gurobi	Mean	1458747	1641182	2118141	2669355	3397775	4526298
	Std.	201433	223689	310292	367954	412785	200
Relative error	Setup	4,32%	3,53%	4,53%	5,97%	9,39%	8,72%