1. Introduction

2. Literature Review

3. Problem description

• 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, 1 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.

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 {\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)

Table 1. Data and parameters.

Table 1. Data and parameters.

Name	Paramenter
T	Number of time periods
	$t\in \left\{1,\cdots ,T\right\}$
${\displaystyle i_0^n},{\displaystyle i_0^r},{\displaystyle i_0^{re}}$	Initial inventory stocks
	Cost per unit and period t
${\displaystyle s_t^n}$	Setup cost for manufacturing new product
${\displaystyle s_t^r}$	Setup cost for remanufacturing product
${\displaystyle p_t^n}$	production cost of new product
${\displaystyle p_t^r}$	production cost of remanufactured product
${\displaystyle h_t^n}$	holding cost of new product
${\displaystyle h_t^r}$	holding cost of remanufactured product
${\displaystyle h_t^{re}}$	holding cost of returned product
	Demand and return in period t
${\displaystyle d_t^n}$	Demand of new product
${\displaystyle d_t^r}$	Demand of remanufactured product
${\displaystyle r_t^{re}}$	Quantity of returned product
	Available Capacities in period t
$c_t$	capacities for manufacturing and
	remanufacturing
	(capacity requirement for new product and
	recovery product is set to one).
	$d_t={\displaystyle d_t^n}+{\displaystyle d_t^r},{\displaystyle \forall t=1,\cdots ,T}$
	${\displaystyle d_{\left[T\right]}^n}={\sum }_{t=1}^T{\displaystyle d_t^n}$
	${\displaystyle d_{\left[T\right]}^r}={\sum }_{t=1}^T{\displaystyle d_t^r}$

Table 2. Decision variables.

Table 2. Decision variables.

Name	Paramenter
${\alpha }_t^n$	1, if new products are manufactured
	in period t; 0, otherwise
$x_t^n$	quantity of new products manufactured
	in period t
${\displaystyle i_t^n}$	inventory stock of new products
	at the end of period t
${\alpha }_t^r$	1, if returned products are remanufactured
	in period t; 0, otherwise
$x_t^r$	quantity of returned products remanufactured
	in period t
${\displaystyle i_t^r}$	inventory stock of remanufactured
	products at the end of period t
${\displaystyle i_t^{re}}$	inventory stock of returned products
	at the end of period t

Table 3. Original Data.

Table 3. Original Data.

t	${\mathit{r}}_{\mathit{t}}^{\mathit{r}\mathit{e}}$	${\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

4. Model A (Relaxation)

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 2 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]. 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.

Algorithm 1:Construction of Halton sequences

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 \text{is a Halton number}\end{array}$$

(35)

We generate a pseudorandom number $0\le g\le 1$ to decide, that values take z, see algorithm 2. 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).

Algorithm 2:Simulation of z with distribution (37)

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

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

Algorithm 3:Calculation of the objective function

6. Numerical Experiments

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, 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

References

1. M. Thierry, M. Salomon, J. Van Nunen, and L. Van Wassenhove. Strategic issues in product recovery management. California Management Review, 37(2):114–135, 1995. [CrossRef]
2. Umweltbundesamt. Papier und druckerzeugnisse, 2020. URL https://www.umweltbundesamt.de/papier-druckerzeugnisserechnungen-amp-prognosen-der-studie-des-ifeu-instituts (Access date: 2020-09-14).
3. F. Sahling. A column-generation approach for a short-term production planning problem in closed-loop supply chains. BuR- Business Research, 6(1):55–75, 2016. [CrossRef]
4. R. Helmrich, M.J. Jans, W. van den Heuvel, and A.P.M. Wagelmans. Economic lot-sizing with remanufacturing: Complexity and efficient formulations. IISE Transactions, 46(1):67–86, 2014. [CrossRef]
5. A. Sifaleras and I. Konstantaras. Variable neighborhood descent heuristic for solving reverse logistics multi-item dynamic lot-sizing problems. Electronic Notes in Discrete Mathematics, 47:69–76, 2015. [CrossRef]
6. O.A. Kilic and W. van den Heuvel. Economic lot sizing with remanufacturing: Structural properties and polynomial-time heuristics. IISE Transactions, 51:12:1318–1331, 2019. [CrossRef]
7. K. Richter and M. Sombrutzki. Remanufacturing planing for the reverse wagner/whitin models. Journal of Operational Research, 121:304–315, 2000.
8. K. Richter and J. Weber. Thre reverse wagner/whitin model with variable manufacturing and remanufacturing cost. International Journal of Production Economics, 71:447–456, 2001. [CrossRef]
9. R.H. Teunter, Z.P. Bayindir, and W. van den Heuvel. Dynamic lot sizing with product returns and remanufacturing. Int. Journal of Production Research, pages 4377–4400, 2006. [CrossRef]
10. T. Schulz. A new silver-meal basic heuristic for the single-item dynamic lot sizing problem with returns and remanufacturing. International Journal of Production Research, page 2519–2533, 2011. [CrossRef]
11. J.O. Cunha and R.A. Melo. A computational comparison of formulations for the economic lot-sizing with remanufacturing. Computers & Industrial Engineering, 92:72–81, 2016. [CrossRef]
12. Z.H. Zhang, H. Jiang, and X. Pan. A lagrangian relaxation based approach for the capacitated lot sizing problem in closed-loop supply chain. Int. J. Productions Economics, 140:249–255, 2012. [CrossRef]
13. M. Florian, J.K. Lenstra, and K. Rinnooy. Deterministic production planning algorithms and complexity. Management Science, 26(7):669–679, 1980. [CrossRef]
14. G.R. Bitran and H.H. Yanasse. Computational complexity of the capacitated lot size problem. Management Science, 28(10):1174–1186, 1982. [CrossRef]
15. P.S. Dixon. Multi-Item Lot-Sizing with Limited Capacity. dissertation, University of Waterloo, Ontario, 1979.
16. B. Karimi, S. Fatemi Ghomi, and J. M. Wilson. The capacitated lot sizing problem review of models and algorithms. OMEGA, 31:365–378, 2003. [CrossRef]
17. J. Maes, JO. McClain, and LN. Van Wassenhove. Multilevel capacitated lotsizing complexity and lp-based heuristics. Journal of Operational Research, 53:131–148, 1991. [CrossRef]
18. I.M. Sobol. Calculation of improper integrals using uniformly distributed sequences. Soviet Math. Dokl., 14:734–738, 1973. [CrossRef]
19. J. H. Halton. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numer. Math., 2:84–90, 1960. [CrossRef]
20. H. Niederreiter. Quasi-monte carlo methods and pseudo-random numbers. Bull. Am. Math. Soc., 84:957–1041, 1978. [CrossRef]
21. B. Klinger. Numerical integration of singular integrands using low-discrepancy sequences. Computing, 59(3):223–236, 1997. [CrossRef]
22. W.V.D. Heuvel and A.P.M. Wagelmans. Four equivalent lot-sizing models. Operations Research Letters, 36:465–470, 2008. [CrossRef]
23. J. Maes and L.N. Van Wassenhove. A simple heuristic for the multi item single level capacitated lotsizing problem. Operations research letters, 4:265–273, 1986. [CrossRef]
24. Donald E. Knuth. The Art of Compute Programming: Seminumerical Algorithms. Volume 2. Addison-Wesley, 1981.
25. J. von Neumann. Various techniques used in connection with random digits. Monte Carlo Method, Appl. Math. Series, 12:36–38, 1951.
26. H. Niederreiter. Random Number Generation and Quasi-Monte Carlo Methods. SIAM, 1992. ISBN 10:0-89871-295-5.
27. S. Tezuka. Uniform Random Numbers: Theory and Practice. Kluwer Academic, 1995.
28. G.J. Chaitin. Information, randomness and incompleteness. Papers on Algorithmic Information Theory., page 236, 1987. [CrossRef]
29. A.N. Kolmogorov and V.A. Uspenskii. Algorithms and randomness. Teor. Veroyatnost i Primenen., 32:425–455, 1987.
30. C.P. Schnorr. Zufälligkeit und wahrscheinlichkeit. Lecture Notes in Math., 218:109, 1971.