A Multiobjective Mathematical Model of Reverse Logistics for Inventory Management with Environmental Impacts: An Application in Industry
aa r X i v : . [ m a t h . O C ] J a n A Multiobjective Mathematical Model of Reverse Logisticsfor Inventory Management with Environmental Impacts:An Application in Industry
M. Forkan ∗ M. M. Rizvi (cid:0) † M. A. M. Chowdhury ‡ January 26, 2021
Abstract
We propose new mathematical models of inventory management in a reverselogistics system. The proposed models extend the model introduced by Nahmias andRivera with the assumption that the demand for newly produced and repaired (reman-ufacturing) items are not the same. We derive two mathematical models and formulateunconstrained and constrained optimization problems to optimize the holding cost. Wealso introduce the solution procedures of the proposed problems. The exactness of theproposed solutions has been tested by numerical experiments. Nowadays, it is an essentialcommitment for industries to reduce greenhouse gas (GHG) emissions as well as energyconsumption during the production and remanufacturing processes. This paper also ex-tends along this line of research, and therewith develops a three-objective mathematicalmodel and provides an algorithm to obtain the Pareto solution.
Keywords
Multiobjective programming problems, Inventory management, Reverse lo-gistics system, Scalarization method, Pareto front.
AMS subject classifications. 90B05, 90C25, 90C29, 90C30.
Reverse logistics (RL) has been defined as a term that refer to the role of logistics inproduct returns, a recycling, the materials substitution, a reuse of materials, a waste dis- ∗ Department of Mathematics, University of Chittagong, Chittagong-4331, Bangladesh, email:[email protected] , † Department of Mathematics, University of Chittagong, Chittagong-4331, Bangladesh, and Universityof South Australia-STEM, Adelaide email: [email protected] (corresponding author), ‡ JNIRCMPS, University of Chittagong, Chittagong-4331, Bangladesh, email: [email protected] osal, and the refurbishing [1]. Moreover, inventory management in reverse logistics, whichincorporates joint manufacturing and remanufacturing options, has received increasing at-tention in recent years. However, fast developments in technology and mass appearance ofnew industrial products, which are coming to the market, have resulted in an increasingnumber of idle products and caused growing environmental problems worldwide. There-fore, increasing ecological concerns, end user awareness, economic considerations, andlegislation, related to waste disposal, encourage manufacturers to take back products af-ter customer have used them. Recently, growing interest and realizations in the reverselogistics processes, such as the recovery of the returned products, have become one of theways, in which businesses endeavor to retain and increase competitiveness in the globalmarket.A number of publications with various models have appeared in the literature aimed tooptimize holding cost in the process of reverse logistics system. McNall [2] and Schrady [3]were the first to address the inventory problem for repairable (recoverable) items. Schrady[3] first explored a deterministic reverse logistic Economic Order Quantity (EOQ) modelfor repairable items with multiple repair cycles and one production cycle. The model ofSchrady [3] was extended by Nahmias and Rivera [4] with inclusion of the case of finiterepair rate. Richter ([5]-[7]) proposed an EOQ model with waste disposal and lookedover the optimal figure of production and remanufacturing batches, depending on the rateof return. Richter [8], Richter and Dobos [9, 10] investigated whether a policy of eithertotal waste disposal or no waste disposal is optimal. Teunter [11] considered multipleproductions and remanufacturing cycles and generalized the results from Schrady [3].Dobos and Richter [12] developed a production/recycling setup with constant demandthat is satisfied by non-instantaneous production and recycling with a single repair and asingle production batch in an interval of time. Later on, Dobos and Richter [13] generalizedtheir earlier model [12] by considering multiple refurbish/repair and production batches ina time interval. Along the same line of study, Dobos and Richter [14] further extended themodel and assumed that the quality of collected used/returned items is not always suitablefor further recycling. Later on, Jaber and El Saadany [15, 16] extended the work of Richter[6, 7] by assuming that the remanufactured items are considered by the customers to beof lower quality than the new ones. Alamri [17] put forward a general reverse logisticsinventory model for deteriorating items by considering the acceptable returned quantityas a decision variable. Singh and Saxena [18] proposed a reverse logistics inventory modelallowing for back-orders. Hasanov et al. [19] extended the work of Jaber and El Saadany[15] by assuming that unfulfilled demand of remanufactured and produced items is eitherfully or partially backordered. El Saadany et al. [20] discussed an inventory model with thequestion as to how many times a product can be remanufactured. Singh and Sharma [21]explored an integrated model with variable production and demand rates under inflation.Later, Singh and Sharma [22] established a production reliability model for deteriorating roducts with random demand and inflation. Recently, Bazan et al. [23, 24] presented amathematical inventory models for reverse logistics with environmental perspectives.In the existing literature, most of the research articles are developed with the assumptionthat the produced and recovered items are not of different quality. In many practical busi-ness situations this hypothesis is not adequate, as the repaired (remanufactured) items areconsidered of secondary quality by the customers, for instance, wheel tyre and computer,etc. Therefore, in this study, it is assumed that newly produced and remanufactured itemsare different in quality. This paper is an extension of the work of Nahmias and Rivera [4]for the case of finite repair rate, production and remanufacturing. Moreover, the solutionapproach of our proposed model enriched and propagated to the case where constraintsof confined storage space in the repair and supply depots are imposed. Note that theportions of products to be procured, repaired, and disposed of in per time unit are fixed,so that, the unit cost of these does not have any impact on optimizing the reverse logisticsmodel. Therefore, our proposed model is associated only with setup and holding costs ofitems of supply and repair depots. Numerical experiments are provided to illustrate theproposed models. The behaviors of the total average holding cost functions for differentstock-out cases are presented with respective tables and graphs.In addition, we developed a multiobjective mathematical model of the reverse logisticssystem that satisfies not only holding cost objective but also environmental requirements.The model considers two environment perspectives; one is to minimize the greenhouse gasemission during the production process, and the other is to optimize energy used in theproduction and remanufacturing processes. Since the problem has complex nature; theeffective methodology and efficient algorithms are needed to approximate the solutionson the front. So far we know that very few researchers have been devoted to solvingmultiobjective reverse logistics model [23, 24]. One of the reasons for not having enoughliterature about the solution approach is that the feasible set is not convex and mighteven be disconnected due to the conflicting nature of the multiple objectives. This posesdifficulties for any technique to approximate the Pareto solutions. We initiate a well-known scalarization approach and algorithms [25] to solve the proposed three-objectivereverse logistics model. Extensive numerical experiments are conducted to approximatethe non-dominated solution of the multiobjective model and the Pareto solutions havebeen demonstrated.In this paper, we have also integrated the applications of our models in the tyre indus-try. In the computational experiments, we consider three numerical examples that arepresented in such a way so that they can align with the industries where the demandsof new and repaired items are different (for example, tyre and IT industries, etc.). Atthe same time, we would like to test our proposed models’ capabilities and analyze theobtained solutions under different parameter settings. We are aware that in the tyre in- ustry the requirements of new and repaired tyres are always different. Customers havedifferent choices to use new and repaired tyres. Some customers prefer the tyre’s quality,and they do not compromise with the safest driving conditions, so they can make a choicefor a new tyre, and others can choose the repaired tyre as the repaired tyre is cheaper andfulfill buyers’ requirements.The rest of this paper comprises the following sections. In Section 2, the new mathe-matical reverse logistics models for unconstrained and constrained optimization problemsand their solutions are illustrated. The extension of these models in the context of multi-objective case is described in Section 3 and in Section 3.1, a solution approach for threeobjectives problem is proposed. In Section 4, numerical experiments are conducted andwe provide the results and discussions based on numerical experiments. The last sectionpresents the conclusion of the paper. We consider two type of depots in the proposed model. The first depot is the supply depotwhich stores the new procurement and repaired items. The user’s demand of primary andsecondary markets can be satisfied from this supply depot. When the items are repaired,first they are stored at the repair depot (second depot) and subsequently shipped in batchesto the supply depot. We assume that shortage is not allowed in this stock point and thelead time is zero.The demand of the new and repaired items are fixed in time but maynot be equal, as the repaired items sell at lower price in the secondary market. Let D p and D r are the demands for new and repaired items respectively in per unit time. Notethat the used items are sent back from user to the overhaul and then repair depot with aconstant rate. The material flow of the model is depicted in Figure 1.It is assumed that the procurement and repair batch sizes are Q p and Q r , respectivelyover the time cycle T . The following parameters have been used in the model. Thecollection percentage of available returns of used items is p (0 < p < r , waste disposal rate 1 − r , λ be the repair rate per unit time with λ > D r , fixedprocurement cost A p , fixed repair batch induction cost A r , holding costs at supply andrepair depots are h and h per unit per unit of time, respectively.In the proposed models, we consider n cycles repair and single-cycle production overthe time T . Therefore, the relations between in and outflows of the stocking points in aprocurement and n repair cycles are Q p + n.Q r = ( D p + D r ) .T and n.Q r = r.p.D p .T he behavior of inventory for produced, collected used and repaired items over the timeinterval T is illustrated in Figure 2. According to the Figure 2, the inventory of the supplydepot always decrease at rates D p and D r for new and repaired items, respectively. Theinventory of the repair depot increase at rate rpDp, and pDp is the portion of the demandwhich is returned to the system for either repair or dumping. The portion of demand(1 − p ) D p never return from primary market to the system, so we assume that it doesnot have any influence in the system. The returned items are usually of varying quality.It is considered that a returned item with a quality less than the acceptance quality isrejected. In the process, the amount of the failed/waste item is (1 − r ) pD p . The decreaseof repaired items during repair is L = (cid:18) Q r λ (cid:19) ( λ − rpD p ). It is mentioned here that, we r Market MarketSecondary PrimaryDepotDepotSupplyRepairProduction r D D p p p pD Dp p D Waste(1-r)p D p D p rp D p RawMaterials
Figure 1: Material flow for a new procurement and repair system.offen introduce various notations in the text to derive the simplest expression for holdingcost function. Let, C = (cid:18) − rpD p λ (cid:19) , where λ > rpD p , this follows that L = C Q r .The repaired item falls below the point S , the repair is suspended, and the inventoryin the repair depot grows at a rate of rpD p for a time T (see, Figure 2). We obtain L − T D r = 0, thus, T = C Q r D r . At the time T a new procurement Q p is received in the supply depot and is used to meetdemand whereas return items accumulated at the repair depot. At the time T when p D p L Q p p r TT DD pD P Q T r p D p B CDE E T L R P P supply depotrepair depotS L Figure 2: The behaviour of inventory for produced, collected used and repaired items overinterval T .the new procurement is all distributed, repaired items are initiated at the supply depot.Therefore, we have T = P Q + QR = T + Q p D p , it follows that, T = C Q r D r + Q p D p . The cycle length T be the time between two new procurements which is also the sameas the time between two successive suspensions of repair. The items fail to repair foreach successive induction at a constant rate (1 − r ) pD p for a period of time P P , that is, T + Q r λ . In this period the net loss is L − L (see, Figure-2) which is C Q r − T rpD p .The total loss of returned item from the system due to item failed to qualify the qualitytest is ( n − C Q r − T rpD p ). It can be seen from Figure 2 that the following relationholds as ( n − C Q r − T rpD p ) = rpD p T − Q r C , his gives n = C Q p Q r , (1)where C = rpC (cid:16) − rpD p D r (cid:17) . The total number of units giving up the supply depot in acycle T is exactly T ( D p + D r ) which must be equal to the total number of units ingoinginto the supply depot in the same time, which is Q p + nQ r , therefore, T = Q p + nQ r D p + D r . (2)From (1) and (2) we have, T = C Q p , where C = 1 + C D p + D r .Now we evaluate the area of the inventory curve of the supply depot during a cycle T .Let A be the area bounded by the inventory curve of the supply depot over the time T.Thus, from Figure 2 we get A = Q p D p + n L (cid:20) T + Q r λ (cid:21) , (3)after substituting L , T and n in (3), gives the area of the supply depot as A = Q p D p + C C Q p Q r (cid:20) C D r + 1 λ (cid:21) . We assume the inventory level in the repair depot is zero when the repair process issuspended. Now we compute the area formed by the inventory curve in the repair depotto find the associate holding costs. Since we are adopting that starting and ending valueson this curve are zero, it ensues that the bounded area in the repair depot can be partedinto triangles and rectangles as displayed in Figure 2.The area of triangle B is 12 rpD p T , by substituting T , we have B = 12 rpD p (cid:18) C Q r D r + Q p D p (cid:19) . Since n cycles of repair inducted in a interval of time T , therefore, there are n triangles C with the total area n Q r λ L , and denoted by, C ′ = 12 λ C C Q p Q r . e also have n − D with the total area n − rpD p T , and denoted by, D ′ = 12 rpD p (cid:18) C Q p Q r − (cid:19) (cid:18) C Q r D r (cid:19) . Rectangle E has an area ( rpD p T − L ) (cid:18) Q r λ + T (cid:19) , this follows, E = Q r (cid:18) λ + C D r (cid:19) (cid:18) rpD p C Q r D r + rpQ p − C Q r (cid:19) . Rectangle E has an area ( L − rpD p T ) (cid:18) Q r λ + T (cid:19) , thus, E = Q r (cid:18) λ + C D r (cid:19) (cid:18) C − rpD p C D r (cid:19) . Total area of the repair depot is A = B + C ′ + D ′ + E + E . It is now pursues that the total average cost for set-up and holding incurred in one cycle,say f ( Q p , Q r ), is given by f ( Q p , Q r ) = A p + nA r + h A + h A . Now we derive an expression of the average holding cost per unit per unit time is deter-mined by first finding the total cost incurred in a cycle and then dividing by the cyclelength. Hence, the average holding cost per unit per unit time, say f ( Q p , Q r ), is createdby taking f ( Q p , Q r ) /T . Therefore, f ( Q p , Q r ) = 1 C Q p ( A p + nA r + h A + h A ) . (4) Theorem 2.1 ( Q ∗ p , Q ∗ r ) is the global minimum of (4) occurs at Q ∗ p = s A p D p h + h pr , (5) and Q ∗ r = vuut λC A r D r C C D r ( h + h ) + 2 D r h pr + λC (cid:16) C C h + 4 h pr + C C D p h prD r (cid:17) . (6) roof. Let us now take the partial derivatives of (4) with respect to Q p and Q r , andthen one can find the optimal procurement and repair batches by making ∂f ∂Q p = 0 and ∂f ∂Q r = 0, these give (5) and (6). ✷ Let us now formulate a mathematical model that minimize the average cost per unit perunit time subject to the constraints. These constraints include available floor space for therepair and supply depots. The problem needs to fulfil these requirements. Suppose p and p be the amount of square feet required to each of the item in supply and repair depots,respectively. We also assume that the availability of maximum number of floor space forsupply and repair depots are k and k , respectively. The highest level of supply inventoryis Q p , thus, the useable floor space is p Q p which is the less than equal of k . Similarly,the maximum level of inventory in the repair depot is the (cid:18) C Q r D r + Q p D p (cid:19) rpD p is less thanequal k . Thus, the proposed mathematical model of the problem with constraints is asfollows. min f ( Q p , Q r )subject to the constraints p Q p ≤ k ,p (cid:18) C Q r D r + Q p D p (cid:19) rpD p ≤ k . (7)The problem is convex, therefore, sufficiency of Karush Kuhn-Tucker (KKT) conditionsholds and we do not need any further conditions [26]. Note that the objective function f ( Q p , Q r ) and constraints are continuously differentiable. Therefore, the KKT conditionsof the Problem (7) is presented as below. Suppose that there exist multipliers λ ≥ λ ≥ ∇ f ( Q p , Q r ) + λ ∇ ( p Q p − k ) + λ ∇ (cid:18) p (cid:18) C Q r D r + Q p D p (cid:19) rpD p − k (cid:19) = 0 , (8) λ ( p Q p − k ) = 0 , (9) λ (cid:18) p (cid:18) C Q r D r + Q p D p (cid:19) rpD p − k (cid:19) = 0 . (10)Taking partial derivatives of the equation (8) with regard to Q p and Q r and then equatingwith zero, thus the solutions form as Q ∗∗ p = s A p D p h + 2 λ D p C p + h pr + 2 λ D p C pp r , (11) nd Q ∗∗ r = s λC A r D r C C D r ( h + h ) + 2 D r h pr + λC C , (12)where C = C C h + 4 h pr + 2 λ D p C pp r + C C D p h prD r . Now we can look on the following four casesCase I: Putting λ = λ = 0 in (11) and (12), we have (5) and (6). We can claim that (5)and (6) are optimal solutions of (7) if the solutions satisfy the constraints of the problem.Case II: If λ = 0 and λ = 0, then we obtain λ = 2 A p D p − Q p ( h + h pr )2 Q p D p C p ,Q ∗∗ p = k p , and Q ∗∗ r is as same as (6). If λ >
0, then we can claim that these are global optimalsolutions of (7) if the solutions satisfy constraints of (7).Case III: If λ = 0 and λ = 0, then we make Q ∗∗ p = s A p D p h + h pr + 2 λ D p C pp r , (13)and Q ∗∗ r = D r ( k − rpQ ∗∗ p p ) rpD p C p . (14)Plugging the solutions (13) and (14) into (8) evaluates λ . If λ > λ = 0 and λ = 0, thereafter by complementary slackness conditions (9) and(10), we reach Q ∗∗ p = k p , (15)and Q ∗∗ r = D r ( k − rpQ ∗∗ p p ) rpD p C p . (16)Combining the solutions (15) and (16) into the system (11) and (12) to assess λ and λ .If λ , λ > Multiobjective Model and Solution Approach
Many environmental factors can arise in the reverse logistics inventory system. For ex-ample, greenhouse gas (GHG) emissions, substantial waste disposal pollution, and energyconsumption can occur during the production of the products. It is required to controlthese environmental effects when one needs to optimize the average holding cost of the re-verse logistics system. In our analysis, we initiate two more objectives such as minimizationof greenhouse gas emission and energy used during the production and remanufacturingprocesses along with the objective of the total average holding cost described in Section 2.Two variables are introduced to construct the mathematical model of the multiobjectivereverse logistics problem, and these are Q p : Procurement batch size , Q r : Repair batch size .Three objective functions for the above model are considered, one is introduced in (4) as: f ( Q p , Q r )= 1 C Q p ( A p + nA r + h A + h A ),and other two functions listed in [24] that illustrated as follows:The greenhouse gas (GHG) emissions (ton per unit) are produced in the production processat a rate P = D p M , where M = 1 − A p D p h Q p > f ( Q p )= a p (cid:18) D p M (cid:19) − b p (cid:18) D p M (cid:19) + c p ! ,where a p an emissions function parameter (ton year /unit ), b p an emissions functionparameter (ton year/unit ) and c p an emissions function parameter (ton year/unit) formanufacturing/production. It is here mentioned that, since repair/remanufacturing rate λ is constant so the second objective function f does not have any effect by Q r , therefore,it is only a function of Q p .And, the third objective function for energy consumed (KWh/year) in the production andremanufacturing process studied in [23] and [28]–[31] is defined as follows: f ( Q p , Q r )= 1 T (cid:20)(cid:18) M W p D p + K p (cid:19) Q p + (cid:18) W r λ + K r (cid:19) nQ r (cid:21) , here W p and K p are respectively, idle power ( KW /year) and energy used (KWh/unit) forthe manufacturing/production and also W r and K r are respectively, idle power ( KW /year)and energy used (KWh/unit) for the remanufacturing process.Thus, we construct the reverse logistics problem as a nonlinear multiobjective optimiza-tion problem as follows: min [ f , f , f ]subject to the constraints p Q p ≤ k ,p (cid:18) C Q r D r + Q p D p (cid:19) rpD p ≤ k , − A p D p h Q p > . (17)The solutions of (17) are called efficient points [32], or Pareto points [26]. A less restrictiveconcept of solution of (17) is the one of a weak Pareto point. We use the followingdefinitions of Pareto point and weak Pareto point. [see details in [33] and [26]].We intorduce the following standard notation to derive the definitions of efficient pointand weak efficient point [34]. Let u, v ∈ R ℓ we write that u ≦ v if and only if u i ≤ v i ∀ i = 1 , . . . , ℓ and ∃ j such that u j < v j ; u < v if and only if u i < v i ∀ i = 1 , . . . , ℓ . Let X be a feasible set of Problem (17). A point ¯ x ∈ X is said to be Pareto for Problem(17) iff there is no x ∈ X , x = ¯ x such that f ( x ) ≦ f (¯ x ). Moreover, a point ¯ x ∈ X is saidto be weak Pareto for Problem (17) iff there is no x ∈ X such that f ( x ) < f (¯ x ). In this section, we recall the scalarization approach introduced in [25] named the Objective-Constraint Approach to solve the proposed Problem (17). The method seems to performefficiently to construct boundary and interior of the Pareto front when the front is discon-nected. For more details on these techniques, see [25] and references therein.
The Objective-Constraint Approach : For ˆ x ∈ X , define the weight w i := 1 /f i (ˆ x ) ℓ X j =1 /f j (ˆ x ) . (18) hen w ∈ W and w k f k (ˆ x ) = w i f i (ˆ x ) , ∀ i = 1 , , ..., ℓ, i = k. The associated scalar problemis defined as( P k ˆ x ) min x ∈ X w k f k ( x ) , subject to w i f i ( x ) ≤ w k f k (ˆ x ) , i=1,2,..., ℓ, i = k. Note that if ¯ x is an efficient solution of Problem (17) and w is as in (18), then ¯ x is thesolution of ( P k ¯ x ) for all k . However, the converse is not true. For weak efficient point theabove property holds for both necessary and sufficient conditions. In this section, we test our proposed models that are stated in (4), (7) and (17), and thesolution formulas introduced in (5)–(6) and (11)–(12) for the range of input parametersettings in the tyre industry. We first introduce Example-4.1 to test the average holdingcost function (4) which depends on a set of parameters. In this problem, we consider theunlimited floor spaces in supply and repair depots, and calculate average holding costs ofprocurement and repair items in per unit time of cycle T . Example-4.2 is commenced forcertain situations where the restriction of floor spaces in supply and repair depots needto be considered. Therefore, we test constraint-model (7) for the same input parametersthat are used in Example-4.1. We also introduce a more challenging problem for themodel (17) where multiple objective functions are considered. As far we know, not enoughmodels and algorithms are presented in the literature of the reverse logistics systems wheremore objective functions are taken, and solutions are reported. By considering this fact,we presented here Problem-4.3 to analyze the proposed multiobjective model (17), andapproximate the solution set of the Pareto front in Figure 7. Example 4.1
Suppose that in the tyre industry, tyre procurement and repair setup costsare A p = $ 10 and A r = $ 30 for batch sizes Q p and Q r , respectively. The demand ofthe new and repaired tyres are D p = 100 and D r = 43 , respectively, over cycle T . Wealso assume that the collection percentage of available returns of used items is p = 0 . ,recovery rate is r = 0 . and repair rate is λ . It is here mention that D r is smaller than λ ,and they both are always larger than rpD p according to our proposed model and as above rpD p = 42 . Let us also assume that holding costs per unit per unit of time for supplyand repair depots are h = $ 1 . and h = $ 1 . , respectively. We intend to find Q p and Q r so that the average holding cost would be the minimum over per unit time of cycle T .Therefore, we aim to minimize the model (4) under the above parameter settings. Now we employ the obtained solution formulas (5) and (6) for the Example-4.1. Accordingto the assumptions of the model (4), in our experiment, we first set λ = 45 which is always able 1: Example-4.1 – Numerical performance of model (4) , and formulas (5) and (6).
Repair Repair Demand for Procurement Repaired Holding Repair Time CPUportions rates repaired items batch sizes batch sizes costs cycles cycles time rpD p λ D r Q ∗ p Q ∗ r f ( Q ∗ p , Q ∗ r ) n T [sec]42 45 43 30.83 115.10 74.61 72.56 58.62 1.1142 60 43 30.83 54.53 156.81 34.15 13.20 0.7342 75 43 30.83 44.92 188.68 28.17 9.07 0.8842 90 43 30.83 40.83 206.80 25.75 7.52 0.6342 105 43 30.83 38.51 218.63 24.10 6.70 0.61 larger than rpD p and D r . As a result, when D r = 43 and λ = 45 we attain the optimalprocurement and repaired batches are Q ∗ p = 30 .
83 and Q ∗ r = 115 .
10 with the optimalaverage holding cost 74 .
61, which are dipicted in Table 1 . In our analysis, we observethat, the optimal cost increases when λ increases ( λ > D r = 43, when λ increases gradually to 105, the optimal solutionand ’average holding ¸costs’ are changes to ( Q ∗ p , Q ∗ r , holdingcosts ) = (30 . , . , . λ = 45 and D r = 43, the repaircycles is 72 (consider the integer number) and the time cycle is 58 .
62 (see, in Table 1),however, when λ increases these two cycles are decreases (see, Figure 3(a)). The aboveanalysis indicates that, when repair rate increases for fixed D r then the inventory ofrepaired items in supply depot increases, and therefore the average holding costs of theitems are also increases. Here, decision makers have the option of choosing an appropriate λ to obtain average holding costs and cycle loops, which may require a trade-off with theoptimum average holding costs. Similar results obtain when D r = 60, this can be seen inFigure 3(c)(d).In our experiments, we employ other kinds of techniques to approximate the solutionof the Example-4.1, for instance, Example-4.1 is solved using a range of nonlinear solverssuch as fmincon with sequential quadratic programming, Ipopt [35], and SCIP [36]. Weobserve that the approximated results vary for solvers and initial guesses. One can beinterested to see the whole set of solutions of Example-4.1. Therefore, we utilize a searchbase algorithm, like Brute Force algorithm, to approximate the solution set. In Figure 4,the straight line formed by the small circles (red) is the approximation of the solutions ofExample-4.1 obtained by search based algorithm.Now we introduce constrained Example-4.2 and test the model (7) and its proposedsolution (11) and (12). Example 4.2
Consider the input parameters that are used in Example-4.1. Moreover, we (a) Holding costs, Repair cycle n , andTime cycle T at D r = 43. (b) Optimal solution ( Q ∗ p , Q ∗ r , costs )at D r = 43 and for 44 < λ <
55 60 65 70 75 80 85 90 95 100 105050100150200250 (c) Holding costs, Repair cycle n , andTime cycle T at D r = 60. (d) Optimal solution ( Q ∗ p , Q ∗ r , costs )at D r = 60 and for 61 < λ < Figure 3: Optimal solutions with average holding costs, repair and time cycles approxi-mations over the selection of repair rates for Example–4.1. Figures (a)(c) are presentedrelations among costs, repair cycle n and time cycle T when D r = 43 and D r = 60, respec-tively, whereas optimum solutions and associated average holding costs ( Q ∗ p , Q ∗ r , costs ) aredepicted in Figures (b) (d) that are approximated for D r = 43 and D r = 60, over thearbitrary interval of λ . a) Solution position on the surfacewhen λ = 45 and D r = 43 .PSfrag replacements f f (b) Solution position on the surfacewhen λ = 60 and D r = 43 . (c) Solution position on the surfacewhen λ = 75 and D r = 43 .Figure 4: (Green) surface produce by (4). (Red) circles represent the solution points of Example-4.1 were obtained by solvers and (Blue) big circle indicates a solution obtained by (5) and (6). assume that the amount of square feets required to each of the item in supply and repairdepots are p = 0 . and p = 0 . , respectively. We also assume that the availability ofmaximum number of floor spaces for supply and repair depots are, k = 20 , k = 10 , infeets, respectively. We intend to find Q p and Q r so that the total average holding cost Sfrag replacements f f (a) Solution position on the surfacewhen λ = 45 and D r = 43. PSfrag replacements f f (b) Solution position on the surfacewhen λ = 75 and D r = 43. Figure 5: (Green) surface is the feasible space of Example-4.2 whereas black surface obtainedfrom Example-4.1 is discarded here because of the floor restrictions. (Red) circles represent thesolution points were obtained by solvers and (Blue) big circle indicates a solution obtained by (11)and (12). would be the minimum per unit time of cycle T . Therefore, we are interested to minimizethe model (7) under the above parameter settings. In the first instance, if we consider D r = 43, λ = 45 and rpD p = 42, then formulas(11)–(12) give the solution of Example-4.2, and thus Case I of the model (7) provides theoptimal procurement and repaired batches which are Q ∗ p = 29 .
77 and Q ∗ r = 115 .
09, andthe optimal average holding cost is 74 .
61 (see Table 2). Whereas we obtain an alternativesolution (30 . , .
5) with the approximate optimal average holding cost 74 .
78 from CaseII. Lagrange multiplier conditions are not satisfied for Cases III and IV therefore, thesolutions are discarded. Furthermore, for fixed D r = 43 and the optimal batch sizes Q ∗ p and Q ∗ r , we also calculate repair and time cycles which are n = 70 (taken only integerpart) and T = 56 .
61. When λ increases within an arbitrary interval 44 ≤ λ ≤
105 thenoptimal average holding costs increases, however two cycles n and T are decreases whichare depicted in Figures 6(a)(b). Therefore, one can cautiously choose the repair ratecompare to the demand rate D r of repaired items, this might be required to trade-offamong average holding costs, repair and time cycles. On the other hand, for D r = 43 and λ that are chosen randomly from { , , , } , the optimal solutions and associatedaverage holding costs are changed that are demonstrated in Table 2. The above analysisindicates that when the repair rate is gradually increased that are far larger than D r ,then the optimal average holding costs would be increased, whereas both repair and time (a) Holding costs, Repair cycly n , andTime cycle T . (b) Optimal solution ( Q ∗ p , Q ∗ r , costs ). Figure 6: Figure (a) is presented optimal average holding costs, repair cycle n and timecycle T over the interval for repair rate 44 < λ <
105 when D r = 43 for Example–4.2,whereas optimum solutions ( Q ∗ p , Q ∗ r , costs ) are depicted in Figure (b) for the same λ ’s and D r .cycles significantly are decreased (see Figure 6(a)). These experiments have also beenconducted by using solvers such as fmincon with sequential quadratic programming andIpopt [35], and we obtain the same set of approximations of the solution. Brute Forcealgorithm is also used to obtain the whole set of solutions. The solution set obtained bythe experiments with Table 2 is depicted in Figure 5 as (red) small circles. Note that,in these input parameter settings, the obtained feasible space and optimal solutions ofconstrained model (7) are not the same as Example-4.1, because floor restrictions areapplied in model (7).In addition, there are requirements to minimize the environmental pollution that arisesduring the production and remanufacturing process in tyre industry. The environmentalimpacts from tyre production include greenhouse gas emissions, energy consumption, dustemission and solvent emission etc. We consider two additional objective functions withExample-4.2 to minimize greenhouse gas emission and energy consumption during theproduction and remanufacturing process, which is presented here as Example-4.3.In Example-4.3, three objective functions are considered which are average holding costfunction f ( Q p ; Q r ), greenhouse gas emissions function f ( Q p ) and energy consumptionfunction f ( Q p ; Q r ). Three objective functions with the constraints set make the Example able 2: Example-4.2 – Numerical performance of model (7) , and formulas (11) and (12).
Repair Repair Demand for Procurement Repaired Holding Repair Time CPUportions rates repair items batch sizes batch sizes costs cycles cycles time rpD p λ D r Q ∗ p Q ∗ r f ( Q ∗ p , Q ∗ r ) n T [sec]42 45 43 29.77 115.09 74.61 70.08 56.61 1.5942 60 43 11.13 52.29 157.78 12.82 4.77 1.0142 75 43 7.28 39.42 193 7.58 2.14 0.6542 90 43 6.26 33.35 215.15 6.35 1.53 1.142 105 43 5.82 30 230.7 5.85 1.27 0.914 (4.3) of model (17) hard to solve. Moreover, the problem has disconnected Pareto front,therefore, we need to choose a scalarization method carefully to solve the problem. Wechoose ( P k ˆ x ) to solve the problem as the method is efficient to approximate the Paretofront when the Pareto front is disconnected. Example 4.3
We make changes in the parameter settings that are used in Examples 4.1and 4.2 as the available differtial solvers have difficulties to approximate solutions forcertain settings. In this Example we assume that r = 0 . , p = 0 . , D p = 1000 , D r = 422 ; λ = 450 , A p = 50 , A r = 100 , h = 20 , h = 10 , p = 1 , p = 1 , k = 2000 , k = 2000 .We set emission function parameters as a p = 0 . (ton year /unit ), b p = 0 . (ton year/unit ) and c p = 1 . (ton year/unit) for production process. Let, the idle powers W p = 120 ( KW/year ) and W r = 80 ( KW/year ) for production and remanufacturingprocess, respectively. Moreover, we set energy K p = 5 . KW h/unit ) , and K r = 2 . KW h/unit ) are used in the production and remanufacturing process, respectively. Weaim to approximate the Pareto points of the multiobjective optimization model (17) underthe above parameter settings. We now introduce the following Algorithm to solve Example-4.3.
Algorithm
We use the objective-constraint approach ( P k ˆ x ) in the algorithm. In Step 2 of Algorithm,each objective function is minimized subject to the original constraints of the Example-(4.3). These individual optimum point are used in Step 3 to form weighted grids. Eachgrid point corresponds to a weight vector in R . In Step 4, three sub-problems are solvedat each grid point to generate Pareto points.In Step 4(b), we calculate the efficient and weak efficient points using the fact that Sfrag replacements f f Figure 7:
The Pareto points are obtained for Example 4.3. if ¯ x = ¯ x = ¯ x =: ¯ x (say) holds, then the solution ¯ x is an efficient. Here, ¯ x k are thesolutions of ( P k ˆ x ) for k = 1 , ,
3. On the other hand, if ¯ x = ¯ x = ¯ x does not hold, thenany dominated point is removed from the set { ¯ x , ¯ x and ¯ x } (see Step 4(b)), and theseare all weak efficient points. The latter case is typically encountered when the Pareto frontand/or the domain is disconnected. Therefore, this algorithm is efficient in finding Paretopoints even when the feasible set is discrete or disconnected. Step 1 (Input)
Set all parameters that stated in Example-4.3.
Step 2 (Determine the individual minima)
Solve Problem min f i , i = 1 , ,
3, subject to the constraints of Example-4.3 thatgive the solutions (cid:0) ¯ Q p i , ¯ Q r i (cid:1) , for i = 1 , ,
3, respectively.
Step 3 (Generate weighted parameters)
Weighted parameters are generated in this step. We generate weighted grids asintroduced in [37, Step 3 of Algorithm 3]
Choose w = ( w , w , − w − w ), which generated from Step 3. (a) Find ˆ x k = ( Q p k , Q r k ) that solves Problems ( P k ˆ x ), k = 1 , , . (b) Determine weak efficient points : (i)
If ¯ x = ¯ x = ¯ x , then set ¯ x = ¯ x (an efficient point)and Record the points. (ii) If ¯ x = ¯ x = ¯ x does not hold, then, any dominated point is discarded bycomparing these 3 solutions.Record non dominated points. Step 5 (Output)All recorded points are Pareto point of Problem-4.3.We write code in MATLAB. We test the range of solvers in Steps 3 and 4(a), theseinclude differential and non-differential solvers such as fmincon with sequential quadraticprogramming algorithm, SCIP [36] and SolvOpt. We take the input parameter valuesthat are listed in the problem statement and implement the proposed algorithm. Theparameter setting plays an important role in Example- 4.3. We choose input parametersin our problem randomly, and for the given parameter setting, we provide 1325 weight-grids ( w ) into the algorithm. As a result, algorithm approximates 2527 Pareto pointswhich are depicted in Figure 7. The elapsed CPU time was about 5 minutes. Note thatthe computations have been performed on a DELL Inspiron 15 7000 laptop with 8 GBRAM and core i7 at 4.6GHz. We developed a new material flow model in the reverse logistics system to measure holdingcost under assumptions that the demands of new and repaired items are different anddeterministic. We established mathematical expressions to compute the holding cost andprovided formulae to optimize the holding cost. We extended the proposed model tothe case where constraints were imposed, and the solution approaches were established.Extensive computational experiments were conducted to demonstrate the efficiency of theproposed models. MATLAB was used to perform the tests, and a wide range of solvers havebeen utilized to solve the complex nature problem. Moreover, we extended the proposedsingle objective problem into a three-objective problem that simultaneously optimizesholding cost, greenhouse gas emissions, and energy consumptions during the productionand remanufacturing processes. Well-known scalarization method employed to solve this hree-objective optimization problem. We successfully approximated the Pareto front ofthis problem, and the computational time has been reported. Acknowledgments:
The first author is supported by NST (National Science andTechnology) fellowship, reference no. 120005100-3821117 in the session: 2019-2020, throughthe Ministry of National Science and Technology, Bangladesh. This financial support ofNST fellowship is gratefully acknowledged.
References
11] Teunter, R. H. (2001) Economic ordering quantities for recoverable item inventorysystems. Naval Research Logistics. 48(6) 484–495.[12] Dobos I, Richter K (2003) A production/recycling model with stationary demand andreturn rates, Central European Journal of Operations Research. 11(1) 35–46.[13] Dobos I, Richter K. (2004) An extended production/recycling model with station-ary demand and return rates, International Journal of Production Economics. 90(3)311–323.[14] Dobos I, Richter K (2006) A production/recycling model with quality consideration,International Journal of Production Economics. 104(2) 571–579.[15] Jaber MY, El Saadney AMA (2009) The production, remanufacture and waste dis-posal model with lost sales, Int. J. Production Economics. 120 115-124.[16] El Saadney AMA, Jaber MY (2010) A production/remanufacturing inventory modelwith price and quality dependant return rate, Computers & Industrial Engineering.58 352-362.[17] Alamri, A. A. (2011) Theory and methodology on the global optimal solution toa general reverse logistics inventory model for deteriorating items and time-varyingrates. Computers & Industrial Engineering, 60(2) 236–247.[18] Singh, S.R. , Saxena, N. (2012) An optimal returned policy for a reverse logisticsinventory model with backorders. Advances in Decision Sciences. Article ID 386598,21 pages, DOI:10.1155/2012/386598.[19] Hasanov, P., Jaber, M.Y., Zolfaghari, S. (2012) Production remanufacturing andwaste disposal model for the cases of pure and partial backordering. Applied Mathe-matical Modelling. 36 (11) 5249–5261.[20] El Saadney AMA, Jaber MY (2013) How many times to remanufacture?, InternationalJournal of Production Economics. 143(2) 598-604.[21] Singh, S.R., Sharma, S. (2013b) An integrated model with variable production anddemand rate under inflation. International Conference on Computational Intelli-gence: Modelling, Techniques and Applications (CIMTA-2013). Procedia Technol-ogy.10 381–391.[22] Singh, S.R., Sharma, S. (2016) A production reliable model for Deteriorating prod-ucts with random demand and inflation. International Journal of Systems Science:Operations & Logistics. Pages 330-338, DOI:10.1080/23302674.2016. 1181221.
23] Ehab Bazan , Mohammad Y. Jaber , Ahmed M. A. El Saadany (2015) Carbon emis-sions and energy effects on manufacturing-remanufacturing inventory models. Com-puters & Industrial Engineering. 88 307-316.[24] Bazan E, Jaber MY, Zanoni S (2016) A review of mathematical inventory modelsfor reverse logistics and the future of its modeling: An environmental perspective,Applied Mathematical Modelling. 40 4151-4178.[25] Burachik RS, Kaya CY, Rizvi MM (2017) A new scalarization technique and newalgorithms to generate Pareto fronts, SIAM Journal on Optimization. 27(2) 1010–1034.[26] Miettinen KM (1999) Nonlinear multiobjective optimization, Boston: Kluwer Aca-demic.[27] Jaber MY, Glock CH, El Saadney AMA (2013a)Supply chain coordination with emis-sions reduction incentives, International Journal of Production Research. 51(1) 69-82.[28] Bogaschewsky R (1995) Nat¨urliche Umwelt und Produktion, Wiesbaden: Gabler-Verlag.[29] Gutowski T, Dahmus J, Thiriez A (2006) Electrical energy requirements for manufac-turing processes, Proceedings of 13th CIRP International Conference on Life CycleEngineering, Leuven, Belgium. 560-564.[30] Mouzon G, Yildirim MB (2008) A framework to minimise total energy consumptionand total tardiness on a single machine, International Journal of Sustainable Engi-neering. 1 105-116.[31] Nolde K, Morari M (2010) Electrical load tracking scheduling of a steel plant, Comput.Ind. Eng. 34(11) 1899–1903.[32] Yu PL (1985) Multicriteria decision making: concepts, techniques, and extensions,New York: Plenum Press.[33] Chankong V, Haimes YY (1983) Multiobjective decision making: Theory andMethodology, Amsterdam: North-Holland.[34] Burachik RS, Kaya CY, Rizvi MM (2014) A new scalarization technique to approxi-mate Pareto fronts of problems with disconnected feasible sets, Journal of Optimiza-tion Theory and Applications. 162(2) 428–446.[35] W¨ a chter A, Biegler LT (2006) On the implementation of a primal-dual interior pointfilter line search algorithm for large-scale nonlinear programming, Mathematical Pro-gramming. 106(1) 25–57.
36] Achterberg T (2009) SCIP: Solving constraint integer programs, Mathematical Pro-gramming Computation. 1(1) 1-41.[37] Burachik RS, Kaya CY, Rizvi MM (2019) Algorithms for Generating ParetoFronts of Multi-objective Integer and Mixed-Integer Programming Problems,arXiv:1903.07041v1.36] Achterberg T (2009) SCIP: Solving constraint integer programs, Mathematical Pro-gramming Computation. 1(1) 1-41.[37] Burachik RS, Kaya CY, Rizvi MM (2019) Algorithms for Generating ParetoFronts of Multi-objective Integer and Mixed-Integer Programming Problems,arXiv:1903.07041v1.