Paper

Annals of Operations Research (2018) 271:1023–1044 https://doi.org/10.1007/s10479-018-3045-2

ORIGINAL RESEARCH

A stochastic reverse logistics production routing model with environmental considerations

Yiqiang Zhang1 · Hussam Alshraideh2,3 · Ali Diabat2,4

Published online: 12 September 2018 © Springer Science+Business Media, LLC, part of Springer Nature 2018

Abstract Growing global concerns of environmental problems have led to the emergence of policies and regulations to control carbon emissions in the industrial sector. These regulations must be taken into consideration to obtain optimal operational decisions on production, inventory and routing in supply chain network models. In this study, we consider the reverse logistics supply chain model with a remanufacturing option to reduce carbon emissions. We aim at providing optimal production, inventory and delivery quantities along with delivery and pickup routes under a carbon cap-and-trade emissions policy. We provide a mathematical formulation of the problem that considers heterogeneous transportation fleets and allows for lost sales under the cap-and-trade carbon emissions policy. The proposed mathematical model is provided in a deterministic and a two-stage stochastic versions to account for demand uncertainty. Proposed formulations are demonstrated through a simulated reverse logistics supply chain with added sensitivity analysis to test for the effect of modeling parameters on the optimal problem solution. Simulation results indicate that carbon policies have significant effect on the supply chain performance with carbon price as the most significant parameter.

Keywords Carbon policy · Carbon emissions · Reverse logistics · Stochastic demand · Production routing

1 Introduction

Industry is the second largest Greenhouse Gases (GHG) emissions contributor worldwide, withelectricityandfossilfuelbasedfluidsastheasthemostcommonlyusedsourcesofenergy

B Ali Diabat Diabat@nyu.edu

1 Pan Shulun Honors College, Shanghai Lixin University of Accounting and Finance, 995 Shangchuan Road, Shanghai 201209, China

2 Division of Engineering, New York University Abu Dhabi, 129188 Saadiyat Island, Abu Dhabi, United Arab Emirates

3 Department of Industrial Engineering, Jordan University of Science and Technology, Irbid, Jordan

4 Department of Civil and Urban Engineering, Tandon School of Engineering, New York University, Brooklyn, NY 11201, USA

123

1024 Annals of Operations Research (2018) 271:1023–1044

(EPA 2017). Emissions from the industrial sector were estimated to be 21% of the global GHG emissions in 2010 and 22% of the US GHG emissions in 2016 (EPA 2017). The industrial sector consists of manufacturing, construction, mining, and agriculture. Manufacturing is the largest of the four industrial sectors and can be broken down into five main categories: paper, food, petroleum refineries, chemicals, and metal/mineral products. These categories account for the vast majority of the fossil fuel use and GHG emissions by this sector (IEA 2017). Manufacturing and industrial processes all combine to produce large amounts of each type of GHGs but specifically large amounts of carbon dioxide (CO2). This is because many manufacturing facilities directly use fossil fuels to create heat and steam needed at various stages of production.

TremendousamountsofharmfulGHGsareemittedasaconsequenceoffossilfuelburning. Such gases have extreme effects on human health and contribute significantly to global worming. GHG are defined as those types of gases that can absorb and emit energy within the thermal infrared range (EPA 2014). Examples of GHG include water vapor and CO2, which are the major constituents of fossil fuel burning outcomes. To reduce the effect of GHG emissions on human health, its emission levels must be controlled. Developed countries such as the US, China and EU countries have developed regulations and policies to control GHG emissions among which carbon cap-and-trade and carbon tax policies are the two basic ones (Jaber et al. 2013; Drake et al. 2016; Song et al. 2017). Under the cap-and-trade policy, the firm is subjected to an emissions cap on the total GHG emissions; extra carbon credits can be bought from the carbon trading market. Salvaged carbon credits can also be sold at the carbon trading market if emissions are lower than the cap. Under the carbon tax policy, a tax is imposed on every unit of the firm’s emissions (Drake et al. 2016; Song et al. 2017). It has been estimated that these regulations will result in about 50 to 200 million tons annual reduction of CO2 emissions (Ellerman and Buchner 2008).

Toavoidextracostofcarbonemissions,firmsmustcontroltheircarbonemissions.Accord- ing to the Environmental Protection Agency (EPA), carbon emissions can be reduced through four basic activities. The activities include energy efficiency through upgrading to more energy efficient industrial technologies, fuel switching through switching to fuels that result in less CO2 emissions but the same amount of energy when combusted, recycling through producing industrial products from materials that are recycled or renewable, rather than pro- ducing new products from raw materials, and training and awareness by making companies and workers aware of the steps to reduce or prevent emissions leaks from equipment (EPA 2017).

In the supply chain context, smart managerial decisions on production quantities, rout- ing, and inventory are needed to cut down carbon emissions while remaining operationally effective. To provide help in this matter, supply chain optimization problems have been con- sidered by many researchers. Examples of such optimization problems include the production inventory routing problem (PIRP) and the closed-loop supply chain problems. In the PIRP, joint optimization of production, inventory, distribution and routing is considered to allow for the best use of available resources while collection, reverse logistics and remanufacturing decisions are considered in closed-loop supply chain problems. Training and awareness to reduce emissions can be practiced through the consideration of production, inventory and routing related supply chain optimization problems. Recycling and remanufacturing as ways to reduce emissions and and to be environmentally friendly are best demonstrated through closed-loop supply chains. Recent advances in the supply chain optimization literature have also considered the closed-loop production routing problem with remanufacturing, simul- taneous pickups and deliveries known as the PRPRPD problem (Qiu et al. 2018). In this problem, a set of manufacturing and remanufacturing facilities are considered with the aim

123

Annals of Operations Research (2018) 271:1023–1044 1025

to reduce the supply chain environmental effects through remanufacturing of worn out prod- ucts. Modeling of the PRPRPD problem provides simultaneous decisions on production, inventory and delivery quantities along with optimal routes for deliveries of finished goods and pickups of worn out items.

Optimization problems of supply chain performance with environmental concerns have also been addressed extensively in the literature. Recent advances in this field include the work of Alkawaleet et al. (2014) where the inventory routing problem is considered along with carbon emissions, the work of Diabat et al. (2016) where the pollution inventory routing problem of perishable goods is considered, and the work of Qiu et al. (2017) where the PIRP has been considered under a carbon cap-and-trade emissions policy. To the best of our knowledge, the PRPRPD with carbon emission concerns has not been addressed in the literature and it is our aim to tackle this problem.

The contributions of our work can be summarized as follows. First, we introduce a real-world variant of the PRPRPD problem that accounts for carbon emission regulations assuming a carbon cap-and-trade policy. Second, we consider the stochastic nature of cus- tomers’ demand and then we introduce a two-stage stochastic formulation that incorporates demand uncertainty. To simulate real-world supply chain systems, we allow for lost sales in our proposed model with penalty cost. Finally, we consider a heterogeneous fleet of trans- portation vehicles in terms of both vehicle capacity and emissions based on type of fuel used.

The rest of the paper is organized as follows. In Sect. 2, we provide a brief review of the literature related to our work. In Sect. 3, we introduce the PRPRPD problem considered in this work and provide the mathematical formulation of the problem in two versions, a deterministic and a two-stage stochastic formulation. Demonstration of our proposed work is shown through a simulated supply chain network in Sect. 4 with added sensitivity analysis. Finally, conclusions and future work remarks are provided in Sect. 5.

2 Literature review

Decisions of production, inventory and routing known as the PIRP problem was first con- sidered by Chandra (1993). He showed that decisions of production, inventory and routing must be made simultaneously to reduce the overall cost of supply chain. After the work of Chandra (1993), the PIRP problem has been extended by several researchers to cover capac- itated vehicles (Fukasawa et al. 2006), heterogeneous fleet (Taillard 1999), multiple plants (Adulyasak et al. 2013), perishable goods (Le et al. 2013), periodic distribution (Diabat et al. 2016) and others.

Environmental impacts of the PIRP problem have been recently considered by several researchers. In their work, Bektaş and Laporte (2011) consider the pollution-routing prob- lem and model emissions based on vehicle load and speed. Their work has been extended to a time version that considers departure time and vehicle speed optimization by Franceschetti et al. (2013). In a survey based review article, Palak et al. (2014) provide an excellent review of work that consider production and inventory management with environmental concerns and analyze the impact of carbon regulations on replenishment decisions in a biofuel supply chain model. In a similar work, Hammami et al. (2015) consider a multi-echelon production inventory model with lead time constraints that incorporate carbon emissions. Inventory rout- ing problem with carbon emissions has been introduced recently in the work of Alkawaleet et al. (2014) and Al Shamsi et al. (2014).

123

1026 Annals of Operations Research (2018) 271:1023–1044

A mathematical formulation of the PIRP has been introduced by Bard and Nananukul (2010) with a branch-and-price heuristic based solution approach. Their modeling approach has shown poor performance and tabu search algorithm was suggested instead (Bard and Nananukul 2010). The PIRP with emission concerns has been considered by Qiu et al. (2017). The authors consider the PIRP problem with carbon cap-and-trade regulation. Their solution approach is based on the use of shortest path dynamic programming algorithms combined with column generation and branch-and-price heuristics.

The PRPRPD problem has received the attention of several researchers as a common prob- lem of interest in electric appliances industry, returnable/reusable transport items industry and beverage industry (Qiu et al. 2018). Examples include the work of Battarra et al. (2014), Subramanian et al. (2013) and Qu and Bard (2014). Integration of inventory decisions with the vehicle routing problem with simultaneous pickups have been explored by Soysal (2016), Van Anholt et al. (2016) and Iassinovskaia et al. (2017).

The importance of remanufacturing in closed-loop supply chains was stressed by Savaskan et al. (2004). Extensions on the importance of remanufacturing have been introduced, such as inventory systems with remanufacturing in the work of DeCroix (2006), economic aspects of remanufacturing in the work of Geyer et al. (2007) and marketing in the work of Atasu et al. (2008). A novel approach to address the problem of closed-loop supply chain with inventory and routing decisions along with remanufacturing options and simultaneous pickups and deliveries was recently introduced by Qiu et al. (2018). Authors formulated the PRPRPD problem as an mixed integer linear program (MILP) and provided a hybrid algorithm that combines branch-and-cut and heuristic search to solve for the optimal solution.

To the best of our knowledge, the PRPRPD problem with carbon emissions has not been addressed yet in the literature. Therefore, we aim at bridging the literature gap by providing a mathematical formulation of the PRPRPD problem with carbon emission constraints under the stochastic demand nature.

3 Problem description and mathematical formulation

3.1 Notation and problem description

In this work, we consider the PIRP with remanufacturing, simultaneous pickups and delivery at the carbon cap-and-trade policy. We assume a closed-loop supply chain network with M manufacturing facilities, R remanufacturing facilities and C retailers/customers. Customers are responsible for selling manufactured and remanufactured products to end users while buying worn out items to be remanufactured at the remanufacturing facilities. Manufacturing facilities are responsible for producing brand new items from raw material to be sold at the customer facilities. A heterogeneous fleet of vehicles is available to transport products from manufacturing and remanufacturing facilities to customer facilities and to pick up worn out items from customer facilities to be delivered to remanufacturing facilities. A carbon emission cap is set for each of the manufacturing and remanufacturing facilities and the goal is to find the optimal production, inventory, pickup and delivery quantities and vehicle routes that maximize the total supply chain network profit for the considered planning horizon.

Three types of items are considered in this work, manufactured items, remanufactured items and worn out items. Manufactured items refer to products produced at the manufac- turing facilities directly from raw material and sold to end customers as brand new items. Remanufactured items are produced at the remanufacturing facilities through re-fabrication

123

Annals of Operations Research (2018) 271:1023–1044 1027

of worn out items and sold to end customers as refurbished products. Worn out items refer to used parts that need to be repaired/renovated at the remanufacturing facilities to be usable again. The following general assumptions are used throughout this work.

• A discrete planning horizon of length |T |. • An available estimate of the demand for both manufactured and remanufactured items.

Demand can be deterministic or stochastic. When stochastic, demand scenarios are assumed available.

• A penalty cost for unmet demand is present. • Each manufacturing and remanufacturing facility has a production capacity with associ-

ated fixed and variable production costs. • In both vehicle capacity and emissions per distance-weight traveled, the heterogeneous

transportation fleet allows for the consideration of various transportation modes • Manufactured items can be held at manufacturing and customer facilities while remanu-

factured and worn out items can be held at remanufacturing and customer facilities with storage capacity and associated holding cost.

• Carbon emissions are generated by manufacturing, remanufacturing, holding and trans- portation operations.

• A cap-and-trade carbon policy is followed such that excess carbon can be sold and extra credits can be bought.

• Not all worn out items can be restored/remanufactured and hence, a remanufacturing rate 0 ≤ ρ ≤ 1 is introduced. Notations for sets, parameters and decision variables of production, inventory and routing

are as follows: Sets:

M the set of manufacturing facilities indexed by m ∈ {1, 2, . . . , |M|}. R the set of remanufacturing facilities indexed by r ∈ {1, 2, . . . , |R|}. C the set of retailers/customers indexed by c ∈ {1, 2, . . . , |C|}. N the set of all network nodes N = M ∪ R ∪ C indexed by n ∈ {1, 2, . . . , |N|}. V the set of vehicles indexed by v ∈ {1, 2, . . . , |V |}. T the set of planning periods indexed by t ∈ {1, 2, . . . , |T |}. Parameters:

C Mm production capacity of manufacturing facility m. C Rr production capacity of remanufacturing facility r. CVv capacity of transportation vehicle v. C F M fixed setup cost for manufacturing facilities. C F R fixed setup cost for remanufacturing facilities. CV M unit manufacturing cost at manufacturing facilities. CV R unit remanufacturing cost at remanufacturing facilities. RM return per unit sold of manufactured products. RR return per unit sold of remanufactured products. P P penalty cost per unit of unmet demand. SMm storage capacity of manufactured items at manufacturing facility m. SRr storage capacity of remanufactured and worn out items at remanufactur-

ing facility r. SCc storage capacity of manufactured and remanufactured items at customer

c. SC Pc storage capacity of worn out items at customer c.

123

1028 Annals of Operations Research (2018) 271:1023–1044

H M unit holding cost of manufactured items at manufacturing facilities. H R unit holding cost of remanufactured items at remanufacturing facilities. HC unit holding cost of manufactured/remanufactured items at customer

facilities. H RP unit holding cost of worn out items at remanufacturing facilities. HC P unit holding cost of worn out items at customer facilities. CC Mm carbon-cap for manufacturing facility m for each time period t. CC Rr carbon-cap for remanufacturing facility r for each time period t. C E M carbon emission per unit produced at manufacturing facilities. C E R carbon emission per unit produced at remanufacturing facilities. C H F carbon emission per unit of holding manufactured/remanufactured items

for one time period. C H P carbon emission per unit of holding worn out products for one time

period. C E P carbon unit price. EVv emissions per unit distance-weight traveled by vehicle v. TCv cost of travel per unit distance-weight traveled by vehicle v. DMct demand of manufactured items at customer c in period t. DRct demand of remanufactured products at customer c in period t. P Rct amount of worn out products (pickups) at customer c in period t. Dnn′ travel distance from node n to node n

′ in the network {n, n′ ∈ N}. W unit weight of manufactured/remanufactured items. ρ remanufacturing rate from worn out items.

Decision variables:

QMmt production quantity at manufacturer m in time period t. QRrt production quantity at remanufacturer r in time t. I Mmt inventory of manufactured items at the end of period t at manufacturing

facility m. I Rrt inventory of remanufactured items at the end of period t at remanufac-

turing facility r. IC Mct inventory of manufactured items at the end of period t at customer c. IC Rct inventory of remanufactured items at the end of period t at customer c. I P Rrt inventory of worn out items at the end of period t at remanufacturing

facility r. I PCct inventory of worn out items at the end of period t at customer c. DLMmcvt quantity of manufactured items delivered to customer c from manufac-

turing facility m by vehicle v in period t. DL Rrcvt quantity of remanufactured items delivered to customer c from remanu-

facturing facility r by vehicle v in period t. DLCcrvt quantity of worn out items delivered from customer c to remanufacturing

facility r by vehicle v in period t. SQMct lost sales (unmet demand) of manufactured items at customer c at time

period t. SQRct lost sales (unmet demand) of remanufactured items at customer c at time

period t. xnn′vt a binary variable that equals one if vehicle v has traveled from node n to

node n′ at time period t and zero otherwise. Where n and n′ are aliases of network nodes m, r and c.

123

Annals of Operations Research (2018) 271:1023–1044 1029

ymt a binary variable that equals one if manufacturing facility m is operating in period t and zero otherwise.

zrt a binary variable that equals one if remanufacturing facility r is operating in period t and zero otherwise.

3.2 Deterministic mathematical formulation

In this section we formulate the considered PRPRPD problem with carbon cap-and-trade as a deterministic MILP model with an overall goal of maximizing total network profit. We refer to this problem as the PRPRPD-cap.

Maximize P = ∑

m

∑

t

(DMct − SQMct ) · RM + ∑

r

∑

t

(DRct − SQRct ) · RR (1a)

− ∑

m

∑

t

(QMmt · CV M + C F M · ymt ) − ∑

r

∑

t

(QRrt · CV R + C F R · zrt ) (1b)

− H M · ∑

m

∑

t

I Mmt − H R · ∑

r

∑

t

I Rrt − HC · ∑

c

∑

t

(IC Mct + IC Rct ) (1c)

− H RP · ∑

r

∑

t

I P Rrt − HC P · ∑

c

∑

t

I PCct (1d)

− ∑

v

∑

t

(∑

m

∑

c

Xmcvt · Dmc · TCv

+ ∑

r

∑

c

DL Rrcvt · Drc · TCv + ∑

c

∑

r

DLCcrvt · Dcr · TCv )

(1e)

− P P · ∑

c

∑

t

(SQMct + SQRct ) (1f)

− (∑

m

∑

t

(C E M · QMmt − CC Mm ) + ∑

r

∑

t

(C E R · QRrt − CC Rr ) )

· C E P (1g)

− C E P · W · ∑

v

∑

t

(∑

m

∑

c

DLMmcvt · Dmc · EVv

+ ∑

r

∑

c

DL Rrcvt · Drc · EVv + ∑

c

∑

r

DLCcrvt · Dcr · EVv )

(1h)

− C H F · C E P · (

∑

m

∑

t

I Mmt + ∑

r

∑

t

I Rrt + ∑

c

∑

t

(IC Mct + IC Rct ) )

(1i)

− C H P · C E P · (

∑

r

∑

t

I P Rrt + ∑

c

∑

t

I PCct

) (1j)

subject to:

QMmt ≤ C Mm · ymt ∀ t ∈ T (2) QRrt ≤ C RR · zrt ∀ t ∈ T (3) I Mmt = I Mmt−1 + QMmt −

∑

c

∑

v

DLMmcvt ∀ m ∈ M, t ∈ T (4)

I Rrt = I Rrt−1 + QRrt − ∑

c

∑

v

DL Rrcvt ∀ r ∈ R, t ∈ T (5)

IC Mct = IC Mct−1 + ∑

m

∑

v

DLMmcvt − DMct + SQMct ∀ c ∈ C, t ∈ T (6)

123

1030 Annals of Operations Research (2018) 271:1023–1044

IC Rct = IC Rct−1 + ∑

r

∑

v

DL Rrcvt − DRct + SQRct ∀ c ∈ C, t ∈ T (7)

I P Rrt = I P Rrt−1 + ∑

c

∑

v

DLCcrvt − QRrt ∀ r ∈ R, t ∈ T (8)

I PCct = I PCct−1 − ∑

r

∑

v

DLCcrvt + P Rct ∀ c ∈ C, t ∈ T (9)

I Mmt ≤ SMm ∀ m ∈ M, t ∈ T (10) I Rrt + I P Rrt ≤ SRr ∀ r ∈ R, t ∈ T (11)

IC Mct + IC Rct ≤ SCc ∀ c ∈ C, t ∈ T (12) I PCct ≤ SC Pc ∀ c ∈ C, t ∈ T (13)

DLMmcvt ≤ xmcvt · CVv ∀ m ∈ M, c ∈ C, v ∈ V , t ∈ T (14)

DL Rrcvt ≤ xrcvt · CVv ∀ r ∈ R, c ∈ C, v ∈ V , t ∈ T (15)

DLCcrvt ≤ xcrvt · CVv ∀ r ∈ R, c ∈ C, v ∈ V , t ∈ T (16)

∑

nn′ xnn′vt ≤ 1 ∀ v ∈ V , t ∈ T (17)

∑

n

xn′nvt ≤ ∑

n

xnn′vt−1 ∀ n′ ∈ N, v ∈ V , t ∈ T (18)

QRrt ≤ ρ · I P Rrt−1 ∀ r ∈ R, t ∈ T (19)

The objective function in (1a–1j) is composed of three main components, revenue from sold items, production, inventory and routing related costs, and emissions related costs. Equation (1a) represents revenue from sold items where the quantity of sold items in period t equals the demand minus the lost sales in that period. Equations (1b) through (1f) represent the production, inventory and routing related cost items and is composed of four components, production cost, inventory holding cost, transportation cost and lost sales penalty cost. Equa- tion (1b) shows the cost of production at both manufacturing and remanufacturing facilities and incorporates both fixed and variable production costs at each facility type. Equations (1c) and (1d) present the inventory holding cost where Eq. (1c) evaluates the cost of holding finished goods at the manufacturing, remanufacturing and customer facilities while Eq. (1d) evaluates the holding cost of worn out items at the remanufacturing and customer facilities. Transportation cost is included in Eq. (1e) and assumes a variable transportation cost as a function of distance traveled and weight of items transported. The transportation cost per vehicle parameter CVv allows for different types of vehicles to be included in the model as new transportation technologies, including full electric and hybrid vehicles, are emerging. Unmet demand penalty cost is included in Eq. (1f) to allow for lost sales as this is a good alternative at high carbon prices.

123

Annals of Operations Research (2018) 271:1023–1044 1031

Carbon emission related costs are incorporated by Eqs. (1g) through (1j) and consist of three main parts, emissions cost due to production, emissions cost due to transportation and emissions cost due to inventory holding. Production related emissions cost is presented in Eq. (1g) and shows the cost of buying extra carbon credits once carbon capacity is exceeded. The terms (C E M · QMmt −CC Mm) and (C E R · QRrt −CC Rr ) in Eq. (1g) can be negative when used emissions are lower than the carbon capacity. In this case, excess carbon will be sold and the return will be added to total revenue as a negative sign is generated. Cost of emissions due to transporting goods between network nodes is presented by Eq. (1h). In this equation, emissions are modeled as a function of distance traveled and weight of transported goods and allows for a heterogeneous fleet of transportation vehicles through the parameter EVv. Measuring units for this parameter are in weight-distance units. Finally, emissions cost due to inventory holding is presented in Eqs. (1i) and (1j). It is assumed that emissions due to holding finished goods are usually higher than those for holding worn out items and hence the two parameters C H F and C H P were introduced such that C H F ≥ C H P. Equation (1i) incorporates emissions cost due to holding finished goods at manufacturing, remanufacturing and customer facilities while Eq. (1j) incorporates emissions cost due to holding worn out items at remanufacturing and customer facilities.

Modeling constraints are shown in Eqs. (2) through (19). Equations (2) and (3) present the production capacity set of constraints for manufacturing and remanufacturing facilities, respectively.Equations(4)through(9)presentinventorybalancesetofconstraintsforfinished goods and worn out items at the manufacturing, remanufacturing and customer nodes of the network. Equations (10) through (13) present storage capacity set of constraints. In this set of constraints, it is assumed that finished goods and worn out items share the same storage place at the remanufacturing facilities while they use different storage places at the customer facilities. Usually, finished goods are displayed separately from storage areas at retail facilities. Vehicle capacity set of constraints are shown in Eqs. (14) through (16). Constraint set in Eq. (17) restricts the assignment of vehicles to one path at each time step t while constraint set (18) ensures that vehicles are available at the path start node before transport assignment. Finally, constraint set in Eq. (19) restricts production quantity at the remanufacturing facilities to available worn out items at a remanufacturing rate ρ.

3.3 Two-stage stochastic mathematical formulation

As demand is of a stochastic nature, its uncertainty must be included when operational deci- sions are to be made. To incorporate demand uncertainty, a two-stage stochastic formulation is proposed next. We assume that the demand of finished goods and worn out items changes over a set of scenarios S. For each scenario s ∈ S, estimates of the demand for manufactured items, remanufactured items and worn out items are available. The following parameters and decision variables are redefined to incorporate available demand scenarios.

Sets:

S the set of demand scenarios indexed by s ∈ {1, 2, . . . , |S|}. Parameters:

DMsct the demand of manufactured items at customer c in time period t under scenario s.

DRsrt the demand of remanufactured items at customer c in time period t under scenario s.

123

1032 Annals of Operations Research (2018) 271:1023–1044

P Rsct the amount of worn out items available for pick up at customer c in time period t under scenario s.

ps occurrence probability of scenario s.

Decision variables:

IC Msct inventory of manufactured items at customer c in time period t under sce- nario s.

IC Rsct inventory of remanufactured items at customer c in time period t under scenario s.

I P Rsrt inventory of worn out items at remanufacturing facility r in time period t under scenario s.

I PCsct inventory of worn out items at customer c in time period t under scenario s. SQMsct demand shortage of manufactured items at customer c in time period t under

scenario s. SQRsct demand shortage of remanufactured items at customer c in time period t

under scenario s. DLCscrvt quantity of worn out items shipped from customer c to remanufacturing

facility r in time period t under scenario s.

Based on these redefined parameters and decision variables along with the previously defined ones, the following two-stage stochastic formulation of the PRPDP-cap is proposed.

Maximize P = ∑

s

ps · (

∑

m

∑

t

( DMct − SQMsct

) · RM + ∑

r

∑

t

( DRct − SQRsct

) · RR )

− ∑

m

∑

t

(QMmt · CV M + C F M · ymt ) − ∑

r

∑

t

(QRrt · CV R + C F R · zrt )

− H M · ∑

m

∑

t

I Mmt − H R · ∑

r

∑

t

I Rrt − HC · ∑

s

ps · ∑

c

∑

t

( IC Msct + IC Rsct

)

− H RP · ∑

s

ps · ∑

r

∑

t

I P Rsrt − HC P · ∑

s

ps · ∑

c

∑

t

I PCsct

− ∑

v

∑

t

(∑

m

∑

c

Xmcvt · Dmc · TCv

+ ∑

r

∑

c

DL Rrcvt · Drc · TCv + ∑

s

ps · ∑

c

∑

r

DLCscrvt · Dcr · TCv )

− P P · ∑

s

ps · ∑

c

∑

t

(SQMct + SQRct )

− (∑

m

∑

t

(C E M · QMmt − CC Mm ) + ∑

r

∑

t

(C E R · QRrt − CC Rr ) )

· C E P

− C E P · W · ∑

v

∑

t

(∑

m

∑

c

DLMmcvt · Dmc · EVv

+ ∑

r

∑

c

DL Rrcvt · Drc · EVv + ∑

s

ps · ∑

c

∑

r

DLCscrvt · Dcr · EVv )

− C H F · C E P · (

∑

m

∑

t

I Mmt + ∑

r

∑

t

I Rrt + ∑

s

ps · ∑

c

∑

t

( IC Msct + IC Rsct

) )

− C H P · C E P · ∑

s

ps · (

∑

r

∑

t

I P Rsrt + ∑

c

∑

t

I PCsct

) (20)

123

Annals of Operations Research (2018) 271:1023–1044 1033

subject to:

QMmt ≤ C Mm · ymt ∀ t ∈ T (21) QRrt ≤ C RR · zrt ∀ t ∈ T (22) I Mmt = I Mmt−1 + QMmt −

∑

c

∑

v

DLMmcvt ∀ m ∈ M, t ∈ T (23)

I Rrt = I Rrt−1 + QRrt − ∑

c

∑

v

DL Rrcvt ∀ r ∈ R, t ∈ T (24)

IC Msct = IC Msct−1 + ∑

m

∑

v

DLMmcvt − DMsct + SQMsct ∀ c ∈ C, t ∈ T (25)

IC Rsct = IC Rsct−1 + ∑

r

∑

v

DL Rrcvt − DRsct + SQRsct ∀ c ∈ C, t ∈ T (26)

I P Rsrt = I P Rsrt−1 + ∑

c

∑

v

DLCscrvt − QRrt ∀ r ∈ R, t ∈ T (27)

I PCsct = I PCsct−1 − ∑

r

∑

v

DLCscrvt + P Rsct ∀ c ∈ C, t ∈ T (28)

I Mmt ≤ SMm ∀ m ∈ M, t ∈ T (29) I Rrt +

∑

s

ps · I P Rsrt ≤ SRr ∀ r ∈ R, t ∈ T (30)

IC Msct + IC Rsct ≤ SCc ∀ c ∈ C, t ∈ T (31) I PCsct ≤ SC Pc ∀ c ∈ C, t ∈ T (32)

DLMmcvt ≤ xmcvt · CVv ∀ m ∈ M, c ∈ C, v ∈ V , t ∈ T (33)

DL Rrcvt ≤ xrcvt · CVv ∀ r ∈ R, c ∈ C, v ∈ V , t ∈ T (34)

DLCscrvt ≤ xcrvt · CVv ∀ r ∈ R, c ∈ C, v ∈ V , t ∈ T (35)

∑

nn′ xnn′vt ≤ 1 ∀ v ∈ V , t ∈ T (36)

∑

n

xn′nvt ≤ ∑

n

xnn′vt−1 ∀ n′ ∈ N, v ∈ V , t ∈ T (37)

QRrt ≤ ρ · I P Rrt−1 ∀ r ∈ R, t ∈ T (38)

4 Simulation study

To demonstrate how our proposed modeling approach works, we consider a reverse logis- tics supply chain with three manufacturing facilities, three remanufacturing facilities and three customer facilities as shown in Fig. 1. As indicated in Fig. 1, production capacities for manufacturing facilities are 1000, 1100 and 700 units while production capacities of reman- ufacturing facilities are 900, 1000, and 800 units. Storage capacity of all considered facilities are shown also in Fig. 1. We assume that manufacturing facilities have carbon caps of 150, 160, 100 metric tons CO2 while remanufacturing facilities have carbon caps of 100, 70 and 80 metric tons CO2 for each time period t ∈ T . There are ten transportation vehicles that can transport goods between network nodes with capacities and initial locations as shown in Fig. 1. Used network modeling parameters are listed in Tables 7, 8, 9 and 10 in “Appendix A”. For simplicity, it was assumed that the distance between any two nodes in the network is 20 miles.

123

1034 Annals of Operations Research (2018) 271:1023–1044

Fig. 1 Simulated reverse logistics supply chain network

Table 1 Values of the revenue and cost items at the optimal solution

Item Optimal value ($)

Total revenue 2,136,000.00

Manufacturing cost 900,000.00

Remanufacturing cost 461,000.00

Inventory holding cost 4812.50

Transportation cost 0.03

Shortage penalty cost 550.00

Production carbon emission cost 16,780.40

Transportation carbon emission cost 3.36

Inventory holding carbon emission cost 5605.76

The problem was coded in GAMS software (GAMS Development Corporation 2013) and solved using CPLEX solver on a 2.8 GHz Windows 10 quad core machine with 16 GB of ram. Optimal solution was achieved in 1390.09 s with a maximum profit of $747,247.96. Optimal cost values obtained are shown in Table 1 while optimal production and delivery quantities are shown in Tables 2 and 3 . As emissions are important here, extra emissions produced at the manufacturing and remanufacturing facilities due to production are shown in Table 4 with a negative sign indicating lower than the carbon cap emissions. These results show that both manufacturing and remanufacturing facilities are working at their almost full capacity to fulfill customer demand. This is because on average the profit from selling finished goods is higher than the price paid for extra emissions used. A total of 585.00 tons CO2 extra emissions was produced by the manufacturing facilities while a total of 350.00 tons CO2 extra emissions was produced by the remanufacturing facilities all at an additional cost of $16,780.40. A graphical representation of the optimal solution obtained is shown in Fig. 2 for time period t = 1.

123

Annals of Operations Research (2018) 271:1023–1044 1035

Table 2 Optimal production quantities at the manufacturing and remanufacturing facilities for the planning period

Facility Time

t1 t2 t3 t4

m1 1000 1000 900 0

m2 800 1100 1100 900

m3 600 700 700 0

r1 500 900 900 900

r2 700 1000 1000 1000

r3 500 0 800 800

Table 3 Optimal delivery and pickup quantities for the planning period

From To Time

t1 t2 t3 t4

m1 c1 1000

c2 900

c3 1500

m2 c1 1000

c2 1000 900

c3 1500

m3 c1 1500

c2 c3 1000

r1 c1 1000 800 900

c2 1000

c3 r2 c1

c2 1000 1000

c3 1200 1000

r3 c1 900

c2 c3 900 800

c1 r1 r2 1100 900

r3 900

c2 r1 950

r2 1000

r3 c3 r1 850 900

r2 r3 700

Three main emission related parameters play a key role in the optimal solution of the PRPDP-cap considered problem. These parameters are the carbon price, carbon capacity and emissions per unit produced. To assess the effect of each of these three parameters on

123

1036 Annals of Operations Research (2018) 271:1023–1044

Table 4 Extra emissions due to production at the manufacturing and remanufacturing facilities for the planning period in metric tons of CO2

Facility Time

t1 t2 t3 t4

m1 100 100 75 − 150 m2 40 115 115 65

m3 50 75 75 − 100 r1 − 25 35 35 35 r2 35 80 80 80

r3 − 5 − 80 40 40

Fig. 2 Optimal solution of the considered supply chain network for time period t = 1

the optimal solution, sensitivity analysis was performed and results are shown in Figs. 3, 4 and 5. Figure 3 shows the profit and total network emissions as a function of carbon price. As it can be seen, a nonlinear relationship between carbon price and profit is present. The profit decreases as the carbon price increases until the average profit from manufacturing facilities equals the price that will be paid for carbon emitted. At this point further, the better decision is to sell the available carbon as it generates more revenue and the profit starts to increase again. Similarly, the total emissions curve decreases incrementally until the carbon price is high enough to be more profitable to sell extra carbon than fulfilling customer demand.

The effect of carbon capacity of production facilities on the total network profit is shown in Fig. 4. The available carbon capacity for each of the production facilities was assumed to be a multiplier (λ) of the current carbon capacity indicated in Fig. 1. Results in Fig. 4 indicate that the profit increases as the carbon capacity increases due to the revenue from selling extra carbon available. The total emissions curve shows no effect of carbon capacity on the amount of production as in all cases at the current carbon price, fulfilling customer demand generates more profit. Figure 5 shows the effect of emissions per unit produced on

123

Annals of Operations Research (2018) 271:1023–1044 1037

Fig. 3 Optimal supply chain network profit and emissions as a function of carbon price

Fig. 4 Optimal supply chain network profit and emissions as a function of carbon capacity

the optimal solution in terms of total profit and total emissions of the supply chain network. Similar to the carbon capacity scenario, emissions per unit produced was set at a multiplier (λ) of the current value of the emissions per unit as indicated in Fig. 1. The network profit curve indicates a nonlinear relationship between emissions per unit and optimal network profit. Profit decreases as emissions per unit increases due to the increased carbon cost of extra emissions generated until the emissions per unit are high enough to generate no profit from selling finished goods. Beyond this point, emissions will be restricted to the carbon caps as no profit will be made by exceeding the carbon caps. Similarly, total emissions increase as emissions per unit increase until emissions per unit are high enough to generate no profit from production. Emissions drop immediately after this point and will be restricted to carbon caps.

123

1038 Annals of Operations Research (2018) 271:1023–1044

Fig. 5 Optimal supply chain network profit and emissions as a function of emissions per unit of manufac- tured/remanufactured product

Table 5 Cost values at the optimal solution with stochastic demand

Item Optimal value

Total revenue 882,000.00

Manufacturing cost 179,000.00

Remanufacturing cost 113,000.00

Inventory holding cost 5095.50

Transportation cost 0.01

Shortage penalty cost 9900.00

Production carbon emission cost − 34,528.90 Transportation carbon emission cost 1.30

Inventory holding carbon emission cost 4995.40

To examine the effect of demand uncertainty on the optimal solution, three demand sce- narios were simulated as shown in Table 11. Scenario 1 shows an optimistic demand estimate while scenario 3 shows a pessimistic demand estimate. The proposed two-stage stochastic formulation in Eqs. (20) through (38) was implemented in GAMS (GAMS Development Corporation 2013) and solved using CPLEX solver. The problem was solved in 4.47 seconds with an optimal profit of $604,536.70. Cost values at the optimal solution are shown in Table 5 while optimal production quantities are shown in Table 6. Table 5 shows a negative carbon emissions cost indicating a less than allowed emissions and hence, a revenue is generated by selling the unused carbon credits.

5 Conclusions and future work

As industry is the second largest contributor to carbon emissions worldwide, developed countries have set strict polices and regulations that control carbon emissions in this sector. Firms must take those regulations into consideration at their operational decisions level to

123

Annals of Operations Research (2018) 271:1023–1044 1039

Table 6 Optimal production quantities for the stochastic demand scenario

Facility Time

t1 t2 t3 t4

m1 750 1000

m2 m3 r1 500

r2 1000

r3 700

cut down cost incurred due to extra emissions generated. In this work, we aim at providing a mathematical model that optimizes available supply chain resources to obtain maximum net- work profit. We consider the PRPRPD supply chain model with added carbon cap-and-trade emissions policy at an overall goal of maximizing total network profit. Optimal simultane- ous decisions on production, inventory and delivery quantities are provided by the model along with optimal delivery and pickup routes under the cap-and-trade carbon emissions policy.

Contributions of this work are summarized as follows. A real-world variant of the PRPRPD supply chain model that accounts for carbon emissions control policy is introduced. The pro- posed model considers also the stochastic nature of demand and, hence, provides both a deterministic and a stochastic versions to help managers in optimizing their supply chain network based on available demand scenarios. To mimic reality, lost sales are allowed in the model with added penalty cost. Finally, the model considers heterogeneous fleet of trans- portation vehicles in both size and emissions based on fuel type.

The proposed formulation has been illustrated through a simulated supply chain network with three manufacturing facilities, three remanufacturing facilities, three customer facilities and ten transportation vehicles. Optimization results have shown that carbon policy is an essential item to be included in the mathematical model with carbon price as the parameter with the highest effect on the total supply chain profit. Emissions per unit parameter had also a significant negative effect on the network profit.

As for future research, several extensions are possible for this work. In one direction, other carbon emission control policies could be investigated to choose the best policy to follow when firms can choose between different emission policies. In another direction, solution methods to the provided mathematical formulation need to be investigated as the model is hard to solve for large supply chain networks.

Appendix A: Simulated network modeling parameters

See Tables 7, 8, 9, 10 and 11.

123

1040 Annals of Operations Research (2018) 271:1023–1044

Table 7 Used values for the simulated network modeling scalar parameters

Parameter Value Units

CFM 2000 $

CFR 1000 $

CVM 100 $

CVR 50 $

RM 120 $

RR 60 $

PP 1 $

HM 0.25 $

HR 0.25 $

HC 0.25 $

HRP 0.2 $

HCP 0.2 $

CEM 0.25 Tons CO2 per unit

CER 0.15 Tons CO2 per unit

CEP 18.44 $ per ton CO2 CHF 0.02 Tons CO2 per unit per time period

CHP 0.01 Tons CO2 per unit per time period

W 20 Kilograms

ρ 1

Table 8 Transportation vehicles emission and travel cost parameters

Vehicle EVv TCv

v1 161.8 0.000024

v2 161.8 0.000024

v3 161.8 0.000024

v4 161.8 0.000024

v5 161.8 0.000024

v6 161.8 0.000024

v7 161.8 0.000024

v8 161.8 0.000024

v9 161.8 0.000024

v10 161.8 0.000024

Average transportation trucks emissions per unit distance weight (EVv) were used and are measured in gCO2 per ton-mile of travel (Mathers et al. 2014). Transportation cost (TCv) values are measured in $ per kilogram- mile of travel assuming an average of 6.64 mpg for transportation trucks, an average gas price of $3.5 per gallon and 20 tons average truck weight (US Department of Energy 2018; Gas Prices 2018)

123

Annals of Operations Research (2018) 271:1023–1044 1041

Table 9 Demand for the four time periods planning horizon

Item type Customer Time

t1 t2 t3 t4

Brand new c1 1200 950 1000 900

c2 700 800 900 900

c3 1300 1100 1200 1000

Remanufactured c1 900 1200 1100 900

c2 850 900 1000 1000

c3 1200 1250 1100 1000

Worn out (pickups) c1 800 900 900 950

c2 900 850 600 700

c3 600 850 700 650

Table 10 Initial inventory at the manufacturing, remanufacturing and customer facilities

Parameter Facility Initial value (units)

IM m1 500

m2 500

m3 500

IR r1 500

r2 500

r3 500

IPR r1 500

r2 700

r3 500

ICM c1 500

c2 500

c3 500

ICR c1 500

c2 500

c3 500

IPC c1 300

c2 300

c3 300

123

1042 Annals of Operations Research (2018) 271:1023–1044

Table 11 Demand under the three scenarios used for the four time periods planning horizon

Scenario Item type Customer Time

t1 t2 t3 t4

1 Brand new c1 1200 950 1000 900

c2 700 800 900 900

c3 1300 1100 1200 1000

Remanufactured c1 900 1200 1100 900

c2 850 900 1000 1000

c3 1200 1250 1100 1000

Worn out (pickups) c1 800 900 900 950

c2 900 850 600 700

c3 600 850 700 650

2 Brand new c1 1000 750 800 700

c2 500 600 700 700

c3 1100 900 1000 800

Remanufactured c1 650 950 850 650

c2 600 650 750 750

c3 950 1000 850 750

Worn out (pickups) c1 500 600 600 650

c2 600 550 300 400

c3 300 850 400 350

3 Brand new c1 600 350 400 300

c2 100 200 300 300

c3 700 500 600 400

Remanufactured c1 300 600 500 300

c2 250 300 400 400

c3 600 650 500 400

Worn out (pickups) c1 400 500 500 550

c2 500 550 200 300

c3 200 650 200 150

References

Adulyasak, Y., Cordeau, J.-F., & Jans, R. (2013). Formulations and branch-and-cut algorithms for multivehicle production and inventory routing problems. INFORMS Journal on Computing, 26(1), 103–120.

Al Shamsi, A., Al Raisi, A., & Aftab, M. (2014). Pollution-inventory routing problem with perishable goods. In: Logistics operations, supply chain management and sustainability (pp. 585–596). Springer.

Alkawaleet, N., Hsieh, Y.-F., & Wang, Y. (2014). Inventory routing problem with co2 emissions consideration. In: Logistics operations, supply chain management and sustainability (pp. 611–619). Springer.

Atasu, A., Sarvary, M., & Van Wassenhove, L. N. (2008). Remanufacturing as a marketing strategy. Manage- ment Science, 54(10), 1731–1746.

Bard, J. F., & Nananukul, N. (2010). A branch-and-price algorithm for an integrated production and inventory routing problem. Computers & Operations Research, 37(12), 2202–2217.

Battarra, M., Cordeau, J.-F., & Iori, M. (2014). Chapter 6: Pickup-and-delivery problems for goods transporta- tion. In: Vehicle routing: Problems, methods, and applications(2nd ed., pp. 161–191). SIAM.

Bektaş, T., & Laporte, G. (2011). The pollution-routing problem. Transportation Research Part B: Method- ological, 45(8), 1232–1250.

123

Annals of Operations Research (2018) 271:1023–1044 1043

Chandra, P. (1993). A dynamic distribution model with warehouse and customer replenishment requirements. Journal of the Operational Research Society, 44(7), 681–692.

DeCroix, G. A. (2006). Optimal policy for a multiechelon inventory system with remanufacturing. Operations Research, 54(3), 532–543.

Diabat, A., Abdallah, T., & Le, T. (2016). A hybrid tabu search based heuristic for the periodic distribution inventory problem with perishable goods. Annals of Operations Research, 242(2), 373–398.

Drake, D. F., Kleindorfer, P. R., & Van Wassenhove, L. N. (2016). Technology choice and capacity portfolios under emissions regulation. Production and Operations Management, 25(6), 1006–1025.

Ellerman, A. D., & Buchner, B. K. (2008). Over-allocation or abatement? A preliminary analysis of the eu ets based on the 2005–06 emissions data. Environmental and Resource Economics, 41(2), 267–287.

Environmental Protection Agency (EPA). (2014). Light-duty automotive technology, carbon dioxide emissions, and fuel economy trends: 1975–2014. Washington, DC. Accessed July 23, 2018, from www.epa.gov/ otaq/fetrends.htm.

Environmental Protection Agency (EPA). (2017). Global greenhouse gas emissions data. Accessed July 23, 2018, from https://www.epa.gov/ghgemissions/global-greenhouse-gas-emissions-data.

Franceschetti, A., Honhon, D., Van Woensel, T., Bektaş, T., & Laporte, G. (2013). The time-dependent pollution-routing problem. Transportation Research Part B: Methodological, 56, 265–293.

Fukasawa, R., Longo, H., Lysgaard, J., de Aragão, M. P., Reis, M., Uchoa, E., et al. (2006). Robust branch-and- cut-and-price for the capacitated vehicle routing problem. Mathematical Programming, 106(3), 491–511.

GAMS Development Corporation. (2013). General algebraic modeling system (GAMS) release 24.2.1. Wash- ington, DC, USA. http://www.gams.com/.

Gas Prices. (2018). National average gas prices. Accessed July 18, 2018, from https://gasprices.aaa.com. Geyer, R., Van Wassenhove, L. N., & Atasu, A. (2007). The economics of remanufacturing under limited

component durability and finite product life cycles. Management Science, 53(1), 88–100. Hammami, R., Nouira, I., & Frein, Y. (2015). Carbon emissions in a multi-echelon production-inventory model

with lead time constraints. International Journal of Production Economics, 164, 292–307. Iassinovskaia, G., Limbourg, S., & Riane, F. (2017). The inventory-routing problem of returnable transport

items with time windows and simultaneous pickup and delivery in closed-loop supply chains. Interna- tional Journal of Production Economics, 183, 570–582.

International Energy Agency (IEA). (2017). CO2 emissions from fuel combustion 2017. Accessed July 23, 2018, from https://www.iea.org/publications/freepublications/publication/ CO2EmissionsfromFuelCombustionHighlights2017.pdf.

Jaber, M. Y., Glock, C. H., & El Saadany, A. M. (2013). Supply chain coordination with emissions reduction incentives. International Journal of Production Research, 51(1), 69–82.

Le, T., Diabat, A., Richard, J.-P., & Yih, Y. (2013). A column generation-based heuristic algorithm for an inventory routing problem with perishable goods. Optimization Letters, 7(7), 1481–1502.

Mathers, J., Wolfe, C., Norsworthy, M., & Craft, E. (2014). The green freight handbook. Environmental Defense Fund.

Palak, G., Ekşioğlu, S. D., & Geunes, J. (2014). Analyzing the impacts of carbon regulatory mechanisms on supplier and mode selection decisions: An application to a biofuel supply chain. International Journal of Production Economics, 154, 198–216.

Qiu, Y., Ni, M., Wang, L., Li, Q., Fang, X., & Pardalos, P. M. (2018). Production routing problems with reverse logistics and remanufacturing. Transportation Research Part E: Logistics and Transportation Review, 111, 87–100.

Qiu, Y., Qiao, J., & Pardalos, P. M. (2017). A branch-and-price algorithm for production routing problems with carbon cap-and-trade. Omega, 68, 49–61.

Qu, Y., & Bard, J. F. (2014). A branch-and-price-and-cut algorithm for heterogeneous pickup and delivery problems with configurable vehicle capacity. Transportation Science, 49(2), 254–270.

Savaskan, R. C., Bhattacharya, S., & Van Wassenhove, L. N. (2004). Closed-loop supply chain models with product remanufacturing. Management Science, 50(2), 239–252.

Song, S., Govindan, K., Xu, L., Du, P., & Qiao, X. (2017). Capacity and production planning with carbon emission constraints.Transportation Research Part E: Logisticsand Transportation Review,97,132–150.

Soysal, M. (2016). Closed-loop inventory routing problem for returnable transport items. Transportation Research Part D: Transport and Environment, 48, 31–45.

Subramanian, A., Uchoa, E., Pessoa, A. A., & Ochi, L. S. (2013). Branch-cut-and-price for the vehicle routing problem with simultaneous pickup and delivery. Optimization Letters, 7(7), 1569–1581.

Taillard, É. D. (1999). A heuristic column generation method for the heterogeneous fleet vrp. RAIRO- Operations Research, 33(1), 1–14.

US Department of Energy. (2018). Average fuel economy of major vehicle categories. Accessed July 18, 2018, from https://www.afdc.energy.gov/data/10310.

123

1044 Annals of Operations Research (2018) 271:1023–1044

Van Anholt, R. G., Coelho, L. C., Laporte, G., & Vis, I. F. (2016). An inventory-routing problem with pickups and deliveries arising in the replenishment of automated teller machines. Transportation Science, 50(3), 1077–1091.

123

Reproduced with permission of copyright owner. Further reproduction prohibited without permission.