An Integrated Approach for Sustainable Supply Chain Management with Replenishment , Transportation , and Production Decisions

Sustainable supply chain management is important for most firms in today’s competitive environment. This study considers a supply chain environment under which the firm needs to make decisions regarding from which supplier and what quantity of parts should be purchased, which vehicle with a certain emissions amount and transportation capacity should be assigned, and what kind of production mode should be used. The integrated replenishment, transportation, and production problem is concerned with coordinating replenishment, transportation, and production operations to meet customer demand with the objective of minimizing the cost. The problem considered in this study involves heterogeneous vehicles with different emission costs, various materials with dissimilar emission costs, and distinct production modes, each with their own emission costs. In addition, multiple suppliers with different quantity discount schemes are considered, different kinds of vehicles with different loading capacities and traveling distance limits are present, and different production modes with different production capacities and production costs are included. A mixed integer programming model is proposed first to minimize the total cost, which includes the ordering cost, purchase cost, transportation cost, emission cost, production cost, inventory-holding cost, and backlogging cost, while satisfying various constraints in replenishment, transportation, and production. A particle swarm optimization model is constructed next to deal with large-scale problems that are too complicated to solve by the mixed integer programming. The main advantage of the proposed models lies in their ability to simultaneously coordinate the replenishment, transportation, and production operations in a planning horizon. The proposed particle swarm optimization model could further identify a near-optimal solution to the complex problem in a very short computational time. To the best of the authors’ knowledge, this is the first paper that considers the sustainable supply chain management problem with multiple suppliers, multiple vehicles, and multiple production modes simultaneously. Case studies are presented to examine the practicality of the mixed integer programming and the particle swarm optimization models. The proposed models can be adopted by the management to make relevant supply chain management decisions.


Introduction
Good supply chain management (SCM) is important for firms to provide low-cost and high-quality products with greater flexibility in today's competitive market, and as a result, to survive and attain a reasonable profit.SCM should cover the management of business activities and relationships within a firm, with immediate suppliers, with the first and second-tier suppliers and customers along the supply chain, and with the entire supply chain [1].Coordination among location, inventory, transportation, and production in a firm and with other partners in a supply chain is necessary [2].Studies of individual topics in SCM have been done abundantly.Among them, inventory management has caught the most attention, and various inventory models and methodologies have been proposed.Transportation problems, such as the vehicle-routing problem, have also been studied, and problems that consider both the production and the transportation aspects have also been found.
In order to confront competition in most industries, firms today often need to satisfy dynamic customer demand spontaneously.As people are becoming more aware of environmental protection, firms also need to consider the greenhouse gas emission of the products.Thus, devising an appropriate integrated plan for replenishment, transportation, and production is essential.To the best of our knowledge, there is no paper that considers the sustainable supply chain management problem with multiple suppliers, multiple vehicles, and multiple production modes concurrently.In this study, an integrated decision-making model is proposed to study the replenishment, transportation, and production problems simultaneously.The materials that are needed by a factory are purchased from multiple suppliers, each with its own ordering cost, unit-purchase cost under different quantity discounts, and carbon emission cost.The plant needs to produce multiple products under different production modes, each with its own production capacity, unit-production cost, and carbon-emission cost.The inventory-holding cost of parts and finished goods, as well as the backlogging cost of parts and finished goods, are considered.Associated with the plant is a heterogeneous fleet of vehicles, each with its own fixed transportation cost, variable transportation cost, and carbon-emission cost.The problem is to determine the operation schedules to coordinate the replenishment, transportation, and production operations so that the customer demand, vehicle travel length, and loading constraints, plant production, and inventory and backlogging constraints are all satisfied, while the total cost (i.e., the sum of the replenishment, transportation, production, and emission costs) over a given planning horizon is minimized.
In this study, a mixed integer programming (MIP) model is first constructed to solve the sustainable supply chain management problem with multiple suppliers, multiple vehicles, and multiple production modes, and optimal solutions can be obtained.Particle swarm optimization (PSO) is proposed next to solve large problems by generating near-optimal solutions.The results show that the MIP model can obtain the optimal solution in a short computational time when the problem is small.However, when the problem becomes relatively large, the MIP may no longer solve the problem in a limited time frame.On the other hand, particle swarm optimization can obtain a near-optimal solution in a short computational time.By applying the proposed models, managers can make relevant supply chain decisions in replenishment, transportation, and production efficiently.
The rest of this paper is organized as follows.Some related works are reviewed in Section 2. The notations and assumptions are presented, and a mixed integer programming (MIP) model and a PSO model are constructed to solve the problem in Section 3. Some case studies are demonstrated in Section 4, and the conclusions are presented in the last section.

Literature Review
Many production management and logistics problems have considered the greenhouse effects and costs related to the environmental impact of operations and transportation activities [3].Reduction in transportation implies that less fossil fuel is burned, and hence a smaller carbon footprint is incurred.Some recent works are as follows.Kannan et al. [4] studied a location/transportation decision problem for reverse logistics activities and constructed a mixed integer linear model for a carbon footprint-based reverse logistics network.Pan et al. [5] studied a freight consolidation to reduce greenhouse gas emissions.The problem considered the pooling of suppliers or retailers with similar flows in similar geographical areas using different transportation modes, such as road and rail, to reduce CO 2 emissions from freight transport.The authors studied the dependency of the emissions produced by the mode of transport on their loads, and defined the emissions function as a piecewise linear and discontinuous function.An optimization model with a piecewise linear objective function was used to calculate the reduction of emissions.Sarkar et al. [6] studied a sustainable integrated inventory problem that considered the unequal power between a vendor and a buyer, a discrete investment that could reduce the setup cost, fixed and variable transportation, and carbon emission costs.A Stackelberg game approach was applied to obtain the global optimum solution over a finite planning horizon.Sarkar et al. [7] studied a three-echelon supply chain problem that comprised a supplier, a manufacturer, and multiple retailers with variable transportation and carbon emission costs.An algebraic approach was applied to solve the problem of minimizing the joint total cost of the supply chain.Toro et al. [3] studied a green capacitated location-routing problem that considered greenhouse gas emissions.A mixed integer linear model was developed to solve the bi-objective problem so that the operational costs and the environmental effects could be minimized.Yuan et al. [8] compared the performances between the supply chain members and designed a contract to make the manufacturer disclose the carbon information.The effects of carbon price, carbon emissions, and carbon quota are quantified in the supply chain model for improving the green supply chain performance.Salehi et al. [9] studied a problem that considered green truck transportation scheduling and driver assignment.A bi-objective mixed integer nonlinear programming model was proposed to minimize the total transportation-related costs and total carbon emissions.A linearization technique was adopted, and a constructive heuristic approach was proposed to solve the problem efficiently.Soysal et al. [10] improved the traditional models for the one-to-one pickup and delivery problem with the consideration of sustainable logistics.Factors, including fuel consumption, variable vehicle speed, and road categorization were considered, and a mixed integer programming model was constructed.Table 1 compares and contrasts the reviewed works with this research.
The inventory-routing problem (IRP) integrates transportation activities with inventory management, and it basically considers the decisions regarding the replenishment quantities and the vehicle routes to visit all of the customers concurrently [11].In the past, such a problem was usually solved by breaking it into an inventory sub-problem and a vehicle routing sub-problem, and then solving each sub-problem independently [3].For example, Juan et al. [12] proposed a simheuristic algorithm for a single-period stochastic inventory-routing problem with stock-outs.The algorithm combined simulation with a randomized heuristics to solve the inventory-routing problem with several stochastic demand inventory problems.The goal was to determine the refill policies and routing plan that minimized the total costs of the system.Ghaniabadi and Mazinani [13] studied a dynamic lot sizing problem (DLSP) with multiple suppliers, backlogging, and quantity discounts by proposing a mixed integer linear programming (MILP) model and a forward dynamic programming (FDP) model.The execution times of the MILP models and FDP models obtained from the recursive formulations are presented and compared.The results demonstrate the efficiency of the FDP models, as they can solve even large-sized instances quite timely.
Particle swarm optimization (PSO) [14] has become one of the popular heuristics for solving production management problems in recent years.Modified versions of PSO have been introduced to attain different objectives.Jordehi and Jasni [15] analyzed existing strategies for setting PSO parameters and provided guidelines for setting parameters in research.Adewumi and Arasomwan [16] proposed two inertia weight strategies, called the swarm success rate descending inertia weight and the swarm success rate random inertia weight, to improve the convergence speed, global search ability, and solution accuracy of the algorithm.Aminbakhsh and Sonmez [17] studied the discrete time-cost trade-off problem by developing a discrete PSO based on the principles for representation, initialization, and position updating of the particles.Liu et al. [18] proposed a hybrid non-parametric PSO algorithm for selecting suitable parameters.Operations, including a multi-crossover operation, a vertical crossover, and an exemplar-based learning strategy, were combined to improve the global and the local exploration capabilities.

Model Development
In this section, assumptions and notations are introduced first, various costs for determining the total cost in a system are presented next, an MIP model is constructed, and the PSO procedure is described last.

Assumptions
The research proposes a decision model for the replenishment of multiple parts from multiple suppliers using multiple vehicles and the production of multiple products using multiple production modes in multiple periods.The assumptions are as follows [19][20][21][22][23]: 1.
Demand of each part in each period is known.The production is based on make-to-order.

2.
There is no beginning inventory in the first period.

3.
Ordering lead time, purchase lead time, transportation lead time, and production lead time are known and set to zero. 4.
Transportation distance and the transportation loading size of each vehicle are fixed and known.

5.
A larger vehicle produces a higher amount of emissions.6.
Each kind of part can be purchased from at least two suppliers and can be purchased from only one supplier in a period.7.
Quantity discount is available.The unit-purchase cost of each kind of part is determined by the quantity of the part purchased in that period.8.
The purchased amount of each kind of part must be delivered in a single batch in a period.9.
The transportation of the ordered parts in a period must be complete in that period.10.At most, one vehicle can travel to and out of a shipment point (supplier) in each period.11.Products can be produced in advance, and backlogging is allowed.12. Different materials incur different amounts of emissions depending on when they were made.13.Different production modes incur different amounts of emissions.

Various Costs
Equation (1) shows the ordering cost, where o vr is the ordering cost of part r from supplier v for each purchase, and α tvr indicates whether an order of part r from supplier v in period t is placed [19]: Equation ( 2) is the purchase cost.Based on the all-units discount brackets from suppliers and the purchase quantity in each period, the total purchase cost of the parts over the horizon can be calculated [19]: where P(Q tvr ) is the unit purchase cost of part r from supplier v in period t, Q tvr is the quantity of part r purchased from supplier v in period t, and α tvr indicates whether an order of part r from supplier v in period t is placed.Equation (3) calculates the transportation cost, which comprises the fixed cost and variable cost.The fixed cost incurs whenever a vehicle is dispatched, and the variable cost depends on the distance the vehicle travels [21]: where FC e is the fixed cost of vehicle e, ϕ e t0j indicates whether vehicle e is dispatched from shipment point 0 (the factory) to shipment point j in period t, c ij is the transportation cost from shipment point i to shipment point j, and ϕ e tij indicates whether vehicle e travels from shipment point i to shipment point j in period t.
Equation (4) calculates the production cost.Based on the production quantity of the finished good g in a period, the production cost in a different production mode s can be calculated [22]: where P ϑ tgs is the unit-production cost of finished good g under production mode s in period t, ϑ tgs is the production quantity of finished good g under production mode s in period t, and ω tgs indicates whether finished good g is manufactured under production mode s in period t.Equation ( 5) calculates the carbon emission cost, which comprises the carbon-emission cost of the vehicle during transportation, the carbon-emission cost of the material, and the carbon-emission cost during production [23].
where u ij is the distance from shipment point i to shipment point j, ϕ e tij indicates whether vehicle e travels from shipment point i to shipment point j in period t, θ e 1 is the carbon-emission cost of vehicle e per distance, Q tvr is the quantity of part r purchased from supplier v in period t, α tvr indicates whether an order of part r from supplier v in period t is placed, θ r 2 is the carbon-emission cost per unit of material r, ϑ tgs is the production quantity of finished good g under production mode s in period t, ω tgs indicates whether finished good g is manufactured under production mode s in period t, and θ s 3 is the carbon-emission cost per unit of product under production mode s.
Equation ( 6) calculates the inventory-holding cost of parts and finished goods.The ending inventory of part r in a period is the sum of the beginning inventory of part r in the period and the purchase quantity of part r in the period, minus the quantity of part r that was used in production in the period.The ending inventory of finished good g in a period is the sum of the beginning inventory of good g in the period and the production quantity of good g in the period, minus the quantity of good g demanded in the period.The inventory-holding cost is as follows [22]: where F + tr is the ending inventory of part r in period t, h r is the unit-holding cost of part r per period, FG + tg is the ending inventory of finished good g in period t, and h g is the unit-holding cost of finished good g per period.
Equation (7) calculates the backlogging cost of parts and finished goods [20]: where F − tr is the backlogging of part r in period t, ε r is the unit-backlogging cost of part r per period, FG − tg is the backlogging of finished good g in period t, and ε g is the unit-backlogging cost of finished good g per period.

Mixed Integer Programming (MIP)
In this sub-section, a mixed integer programming (MIP) model is developed to solve the multiple-supplier, multiple-product replenishment problem, and devise the replenishment plan and production mode in each period in the planning horizon.The MIP model is as follows: Subject to ∑ I i=0 ϕ e tij − ∑ I i=0 ϕ e tji = 0, for all t, e, j = 1, 2 . . ., J, j = i (18) ∑ I i=1 δ e ti0 ≤ w e , for all t, e S tg = ∑ S s=1 ϑ tgs × ω tgs , for all t, g (28) α tvr ∈ {0, 1}, for all t, v, r (31) ϕ e tij ∈ {0, 1}, for all t, i, j, e φ e ti ∈ {0, 1}, for all t, i, e (34) Objective function ( 8) is to minimize the total cost (TC) in the planning horizon.The costs include nine kinds of costs: ordering cost, purchase cost, transportation cost, production cost, three types of carbon-emission costs, holding cost, and backlogging cost.Equation ( 9) calculates the total quantity of part r purchased in period t, S tr , where Q tvr is the quantity of part r purchased from supplier v in period t, and α tvr indicates whether an order of part r from supplier v in period t is placed.Equation ( 10) makes sure that an order of part r from supplier v in period t must be placed so that the parts can be purchased from that supplier in that period.Equation ( 11) determines the unit-purchase cost of part r from supplier v in period t, P(Q tvr ), which is based on the quantity discount brackets from that supplier and whether an order is placed from that supplier.Equation ( 12) ensures that only one quantity discount bracket can be applied for an order of part r from supplier v in period t.Equation ( 13) makes sure that the quantity of part r purchased from supplier v in period t is in the correct quantity discount bracket x.Equation ( 14) ensures that the quantity of all of the parts purchased from a supplier in a period must be the purchase size of the shipment from that supplier in that period.Equation (15) shows that only one vehicle can travel to shipment point j (excluding the factory, i.e., j = 0) in period t.Equation ( 16) makes sure that at most, one vehicle can travel from shipment point i (excluding the factory, i.e., i = 0) to one single shipment point in period t.Equation (17) ensures that each vehicle, if dispatched, can only travel starting from the factory (i = 0) in each period.Equation (18) shows that all of the vehicles that travel from the factory (i = 0) will go back to the factory in each period.Equation (19) ensures that if vehicle e travels from shipment point i (excluding the factory, i.e., i = 0) it will travel to only one shipment point j directly in each period.Equation (20) ensures that in each period, at most, one vehicle can travel from shipment point i (excluding the factory, i.e., i = 0).Equation ( 21) calculates the loading size from shipment point i in period t, Y ti , which is the added loading of vehicle e in that shipment point.That is, it measures the difference in the loading sizes of the vehicle between the departure and arrival.Equation (22) ensures that for each vehicle e in each period, the sum of the purchase sizes from all of the shipment points must be equal to the loading size of vehicle e from all of the shipment points back to the factory.Equation (23) ensures that for each period, the loading size of all of the vehicles from all of the shipment points back to the factory must be equal to the sum of the purchase sizes from all of the shipment points.Equation (24) states that the total loading of vehicle e back to the factory in period t must be less than or equal to the maximum loading size of that vehicle.Equation (25) states that the total traveling length of vehicle e in period t must be less than or equal to the maximum traveling length of that vehicle.Equation (26) ensures that the production quantity of finished good g under production mode s in period t must be less than or equal to the quantity that can be produced under production mode s.Equation ( 27) calculates the ending inventory (FG + tg ) or the backlogging (FG − tg ) of finished good g in period t by deducting the accumulated demand of that good from the first period to period t from the accumulated production quantity of that good from the first period to period t.Equation (28) calculates the total quantity of finished good g that was produced in period t (S tg ) by summing up the production quantity of that good under all of the production modes in that period.Equation (29) calculates the ending inventory (F + tr ) or the backlogging (F − tr ) of part r in period t by deducting the accumulated demand of that part from the first period to period t from the accumulated purchase quantity of that part from the first period to period t.Equation (30) calculates the demand of part r in period t (d tr ) by multiplying the total quantity of finished good g produced in period t (S tg ) by the units of material r that are required to produce product g (ρ gr ).Equations (31) to (35) define some variables as binary variables.Equation ( 31) is to set α tvr , which indicates whether an order of part r from supplier v in period t is placed, as a binary variable.Equation ( 32) is to set β tvrx , which indicates whether an order of part r under a quantity discount bracket x from supplier v in period t is placed, as a binary variable.Equation ( 33) is to set ϕ e tij , which indicates whether vehicle e travels from shipment point i to shipment point j in period t, as a binary variable.Equation ( 34) is to set φ e ti , which indicates whether vehicle e travels from shipment point i in period t, as a binary variable.Equation ( 35) is to set ω tgs , which indicates whether finished good g is manufactured under production mode s in period t, as a binary variable.

Particle Swarm Optimization (PSO)
In this research, the PSO procedure is developed based on the constriction factor proposed by Kennedy and Eberhart [14] and linear decreasing inertia weight proposed by Shi and Eberhart [24].Based on the procedure, the speed of convergence can be increased and the local and global search capabilities can be improved.The steps are as follows [25,26]: Step 1. Initialize particles with random positions and velocities.With a search space of d-dimensions, a set of random particles (solutions) is first initialized.Let the lower and the upper bounds on the variables' values be λ min and λ max .We can randomly generate the positions, λ τ n (the superscript denotes the τ th particle, and the subscript denotes the n th iteration), and the exploration velocities, µ τ n , of the initial swarm of particles: where the positions and exploration velocities are in a vector format, rand is a random number between 0 and 1, and ∆ is the constant time increment, and is assumed to be 1.

•
Step 2. Evaluate the fitness of all of the particles.The performance of each solution is evaluated with the fitness function.

•
Step 4. Keep track of the locations where each individual has its highest fitness.

•
Step 5. Keep track of the position with the global best fitness.

•
Step 6. Update the velocity of each particle: where ω n is the inertia factor, µ τ n is the velocity of the τ th particle at the n th iteration, ϕ 1 and ϕ 2 are the acceleration constants toward pbest and gbest, rand 1 and rand 2 are random numbers between 0 and 1, pbest e n is the best searching experience of the τ th particle so far at the n th iteration, gbest n is the best result obtained among all of the particles at the n th iteration, λ τ n is the current position of the τ th particle, and ω n can be set as a constant value or a variable changing in all of the iterations.

•
Step 7. Update the position of each particle:

•
Step 9. Terminate the process if a maximum number of iterations is attained.Otherwise, go to Step 2.

Case Studies
The proposed MIP and PSO models for joint replenishment, transportation, and production planning are applied in case studies here.The MIP model is implemented using the software LINGO 10, and the PSO is implemented using the software MATLAB.The proposed models aim to determine the most appropriate replenishment amounts of the parts from various suppliers in different periods, the routing and loading of vehicles in different periods, and the amount of each kind of product manufactured under different production modes in different periods.

Data
In the case studies, a machine tool manufacturer in Taichung, Taiwan, is used as an example.The manufacturer is selected because its production environment is suitable to the proposed model, and the data is collected and revised based on the assumptions of the problem setting.The manufacturer needs to devise its integrated replenishment, transportation, and production plan for spindles.Three major parts need to be purchased: spindle shaft, shaft sleeve, and bearing.Each part can be purchased from two suppliers.The information of the cases is as follows.Table 2 shows the ordering cost from each supplier.Table 3 shows the unit cost under various quantity discounts from each supplier.Table 4 shows the inventory-holding cost of each material and finished good.Table 5 shows the fixed cost, maximum loading size, and maximum traveling length of each vehicle.Table 6 shows the distance and transportation cost matrix among the factory and the suppliers.Table 7 shows the production cost of spindles under different production modes.Table 8 shows the part requirements for the finished goods.In addition, the carbon emission costs of a small vehicle (e = 1) and large vehicle (e = 2) are $90 and $100 per kilometer, respectively.The carbon emission costs per unit of spindle shaft (r = 1), shaft sleeve (r = 2), and bearing (r = 3) are $10, $11, and $12, respectively.The carbon emission costs per unit of product under normal production (s = 1), overtime production (s = 2), and outsourcing (s = 3) are $20, $20, and $30, respectively.The unit backlogging cost of a part was $50 per period, and for a finished good was $400 per period.

Case I
A simplified case is presented here.The planning horizon contains three periods, and the demand of the finished good in each period is shown in Table 9 Case I. Assume that the firm only produces one kind of product, a basic spindle (g = 1), which require two kinds of parts: the spindle shaft (r = 1) and the shaft sleeve (r = 2).The spindle shaft (r = 1) can be purchased from supplier 1 (v = 1) or supplier 2 (v = 2), and the shaft sleeve (r = 2) can be purchased from supplier 3 (v = 3) or supplier 4 (v = 4).For the transportation, a small vehicle and a large vehicle can be assigned.Three production modes are available: normal (s = 1), overtime (s = 2), and outsourcing (s = 3).For the PSO, the number of particles is set to be 150, and the number of iterations is set to be 1000.The results from the MIP and the PSO are the same, as shown in Table 10.Table 10 shows that an order of spindle shaft (r = 1) from supplier 1 (v = 1) in period 1 (t = 1) is placed, i.e., α 111 = 1, and 360 units of spindle shaft (r = 1) are purchased from supplier 1 in period 1, i.e., Q 111 = 360.A quantity discount can be obtained, and the ordering cost and transportation cost can be reduced as a result.In addition, an order of shaft sleeve (r = 2) from supplier 3 (v = 1) in period 1 (t = 1) is placed, and 360 units of shaft sleeve (r = 2) are purchased from supplier 3 in period 1.In period 1, 100 units of basic spindle (g = 1) are produced under normal production (ϑ 111 = 100), and 30 units of basic spindle (g = 1) are produced under overtime production (ϑ 112 = 30).The same applies to period 2; that is, ϑ 211 = 100, ϑ 212 = 30.In period 3, 100 units of basic spindle (g = 1) are produced under normal production, i.e., ϑ 311 = 100.A large vehicle (v = 2) travels from the factory to supplier 3, then to supplier 1, and back to the factory in period 1, i.e., ϕ 2 103 = 1, ϕ 2 131 = 1, ϕ 2 110 = 1.In addition, the loading size of the large vehicle from the factory to supplier 3 is 0, from supplier 3 to supplier 1 is 360 units, and from supplier 1 to the factory is 720, i.e., δ 2 103 = 0, δ 2 131 = 360, δ 2 110 = 720.The ending inventory of the spindle shaft (r = 1) at the end of period 1 is 230 units, i.e., F + 11 = 230, and that of shaft sleeve (r = 2) is also 230, i.e., F + 12 = 230.The ending inventory of the spindle shaft and shaft sleeve at the end of period 2 are both 100 units.Finally, the ending inventory (FG + 11 ) of finished good 1 in period 1 is 18 units, and the backlogging (FG − 21 ) of the finished good in period 2 is 13 units.Table 10 also shows each kind of cost and the total cost of the firm in the horizon.Since both the MIP and the PSO models lead to the same results, the total cost is the same, i.e., $7,954,180.Figure 1 shows the vehicle routing under Case I.The computational time for the MIP is 21 seconds.
The PSO model is implemented using the software MATLAB (2015).The variant is distinguished in the literature due to its theoretical properties that imply the following explicit selection of the parameters by Clerc and Kennedy [27].The criteria for a performance evaluation of the PSO algorithm include: minimize the computational time; minimize the error; and reduce the use of variables, such as the decision variable, auxiliary variables, and 0-1 variables.A penalty term is added to the objective function, and infeasible solutions will be penalized.In this case, the settings ω n = 0.729 and ϕ 1 = ϕ 2 = 2.05 are currently considered as the default parameter set of the constriction coefficients of the variant.The initialized particle swarm size is selected as 50, and the maximum evolution number of the particle swarm is 1000.The PSO converges after the 19 th generation.The computational time for the PSO is 227 s. constriction coefficients of the variant.The initialized particle swarm size is selected as 50, and the maximum evolution number of the particle swarm is 1000.The PSO converges after the 19 th generation.The computational time for the PSO is 227 s.

Case II
In Case II, the planning horizon contains five periods, and the hybrid spindle (g = 2) is produced.The demand of the finished good in each period is shown in Table 9.The firm needs to purchase three kinds of material, i.e., spindle shaft (r = 1), shaft sleeve (r = 2), and bearing (r = 3).Each kind of material can be purchased from two suppliers, that is, supplier 1 (v = 1) and supplier 2 (v = 2) for spindle shaft (r = 1), supplier 3 (v = 3) and supplier 4 (v = 4) for shaft sleeve (r = 2), and supplier 5 (v = 5) and supplier 6 (v = 6) for bearing (r = 3).For the transportation, a small vehicle and a large vehicle can be assigned.Three production modes are available.For the PSO, the number of particles is set to be 200, and the number of iterations is set to be 1000.The results from the MIP and the PSO models are different.The results from the MIP are shown in Table 11, and those from the PSO are shown in Table 12.For example, an order of spindle shaft (r = 1) from supplier 1 (v = 1) in period 1 (t = 1) is placed with an amount of 121 units under the MIP model; however, the order is placed from supplier 2 (v = 1) with an amount of 100 units under the PSO model.Nevertheless, the production quantities of the finished good under different production modes in different periods, i.e., , are the same under the MIP model and the PSO model.The selection of the vehicle type, the routings of the vehicles, and the loading of the vehicles in various periods are different under the two models.The ending inventories of various parts in various periods are different under the

Case II
In Case II, the planning horizon contains five periods, and the hybrid spindle (g = 2) is produced.The demand of the finished good in each period is shown in Table 9.The firm needs to purchase three kinds of material, i.e., spindle shaft (r = 1), shaft sleeve (r = 2), and bearing (r = 3).Each kind of material can be purchased from two suppliers, that is, supplier 1 (v = 1) and supplier 2 (v = 2) for spindle shaft (r = 1), supplier 3 (v = 3) and supplier 4 (v = 4) for shaft sleeve (r = 2), and supplier 5 (v = 5) and supplier 6 (v = 6) for bearing (r = 3).For the transportation, a small vehicle and a large vehicle can be assigned.Three production modes are available.For the PSO, the number of particles is set to be 200, and the number of iterations is set to be 1000.The results from the MIP and the PSO models are different.The results from the MIP are shown in Table 11, and those from the PSO are shown in Table 12.For example, an order of spindle shaft (r = 1) from supplier 1 (v = 1) in period 1 (t = 1) is placed with an amount of 121 units under the MIP model; however, the order is placed from supplier 2 (v = 1) with an amount of 100 units under the PSO model.Nevertheless, the production quantities of the finished good under different production modes in different periods, i.e., Q tvr , are the same under the MIP model and the PSO model.The selection of the vehicle type, the routings of the vehicles, and the loading of the vehicles in various periods are different under the two models.The ending inventories of various parts in various periods are different under the two models, too.The ending inventories or backlogging of the finished good in various periods are the same under the two models.The total cost under the MIP model is $15,244,600.The total cost under the PSO model is $15,663,510, which is 2.75% higher than that under the MIP model.The vehicle-routing results from the MIP and the PSO models are also shown in Figures 2 and 3.The computational time for the MIP model is 417 seconds.In addition, Figure 4 shows the PSO result generated from MATLAB; the optimal solution is obtained at the 113 th generation, and the computational time is 325 s.

Case III
The planning horizon contains nine periods.The firm produces two products, a basic spindle (g = 1) and a hybrid spindle (g = 2), and the demand of the finished good in each period is shown in Table 9.The firm needs to purchase three kinds of material: spindle shaft (r = 1), shaft sleeve (r = 2), and bearing (r = 3).Each kind of material can be purchased from two suppliers, that is, supplier 1 (v = 1) and supplier 2 (v = 2) for the spindle shaft (r = 1), supplier 3 (v = 3) and supplier 4 (v = 4) for the shaft sleeve (r = 2), and supplier 5 (v = 5) and supplier 6 (v = 6) for the bearing (r = 3).For the transportation, one small vehicle and one large vehicles can be assigned.Three production modes are available.Due to the increasing number of variables and constraints, the problem becomes non-deterministic polynomial hard (NP-hard), and the MIP cannot obtain the optimal solutions.Therefore, only the PSO is applied in Case III.The number of particles is set to be 250, and the number of iterations is set to be 1000.The results from the PSO are shown in Table 13 and Figure 5.The execution result from the PSO shows that the best generation occurs at the 213 th generation.The computational time for the PSO is 471 s.

Case III
The planning horizon contains nine periods.The firm produces two products, a basic spindle (g = 1) and a hybrid spindle (g = 2), and the demand of the finished good in each period is shown in Table 9.The firm needs to purchase three kinds of material: spindle shaft (r = 1), shaft sleeve (r = 2), and bearing (r = 3).Each kind of material can be purchased from two suppliers, that is, supplier 1 (v = 1) and supplier 2 (v = 2) for the spindle shaft (r = 1), supplier 3 (v = 3) and supplier 4 (v = 4) for the shaft sleeve (r = 2), and supplier 5 (v = 5) and supplier 6 (v = 6) for the bearing (r = 3).For the transportation, one small vehicle and one large vehicles can be assigned.Three production modes are available.Due to the increasing number of variables and constraints, the problem becomes non-deterministic polynomial hard (NP-hard), and the MIP cannot obtain the optimal solutions.Therefore, only the PSO is applied in Case III.The number of particles is set to be 250, and the number of iterations is set to be 1000.The results from the PSO are shown in Table 13 and Figure 5.The execution result from the PSO shows that the best generation occurs at the 213 th generation.The computational time for the PSO is 471 s.

Results and Discussions
The aim of this research is to develop an integrated operations plan that minimizes the replenishment, transportation, production, and emission costs, subject to all of the constraints involved.In the past, only several kinds of costs were considered in relevant problems [19][20][21][22][23].However, in this study, we take seven different kinds of costs into consideration, i.e., ordering cost, purchase cost, transportation cost, production cost, carbon-emission cost, holding cost, and backlogging cost.As shown in Table 1, this research is more comprehensive compared to other works.It considers the green supply chain, the transportation problem, and the emission issue, and constructs a mathematical model and an algorithm aiming to find the global optimum.Therefore, integrated models for a sustainable supply chain are developed in this study.For small to medium-scale problems, the MIP models are constructed to solve the problem using the software LINGO.The solutions are the global optimums.However, when the scale of the problem increases, the problem becomes NP-hard.Thus, for medium to large-scale problems, the PSO can be applied to solve the problems efficiently and find near-optimal solutions.Sensitivity analysis is applied to Case II to show the robustness of the proposed MIP model.Table 14 shows the effects of parameter changes on the total cost of the system.When the value of , the fixed cost of vehicle e per trip, changes, the assignment of different vehicles in each period changes.Thus, when increases by 50%, the total cost increases from $15,244,600 to $15,249,250.

Results and Discussions
The aim of this research is to develop an integrated operations plan that minimizes the replenishment, transportation, production, and emission costs, subject to all of the constraints involved.In the past, only several kinds of costs were considered in relevant problems [19][20][21][22][23].However, in this study, we take seven different kinds of costs into consideration, i.e., ordering cost, purchase cost, transportation cost, production cost, carbon-emission cost, holding cost, and backlogging cost.As shown in Table 1, this research is more comprehensive compared to other works.It considers the green supply chain, the transportation problem, and the emission issue, and constructs a mathematical model and an algorithm aiming to find the global optimum.Therefore, integrated models for a sustainable supply chain are developed in this study.For small to medium-scale problems, the MIP models are constructed to solve the problem using the software LINGO.The solutions are the global optimums.However, when the scale of the problem increases, the problem becomes NP-hard.Thus, for medium to large-scale problems, the PSO can be applied to solve the problems efficiently and find near-optimal solutions.Sensitivity analysis is applied to Case II to show the robustness of the proposed MIP model.Table 14 shows the effects of parameter changes on the total cost of the system.When the value of FC e , the fixed cost of vehicle e per trip, changes, the assignment of different vehicles in each period changes.Thus, when FC e increases by 50%, the total cost increases from $15,244,600 to $15,249,250.When the value of θ e 1 , θ r 2 , θ s 3 , c ij , o vr , h r or h g increases (decreases) by up to 50% individually, the values for the decision variables do not change, and the total cost changes due to the changes in the specific parameter.When the value of ε r , the unit-backlogging cost of part r per period, changes, the total cost does not change.This is because there is no backlogging of parts in all of the periods.When the value of ε g , the unit-backlogging cost of finished good g per period, changes, the production quantity, production mode, and backlogging of finished goods in different periods change, too.Based on the proposed models, managers can devise their supply chain plans.For the replenishment decisions, the managers can determine when the firm should purchase their parts, from which supplier(s), with what quantity, and at what unit-purchase cost.For the production decisions, the managers can determine when and how many units the firm should produce of their products, what kind of production mode to use, and at what unit-production cost.For the transportation decisions, the managers can determine when and which vehicle(s) should be assigned to certain shipment points.The solutions can enable managers to minimize the total costs of the firm.

Conclusions
This research proposes two models to the replenishment, transportation, and production problem where the vehicle routing must be optimized together with the production modes and the replenishment policies.This research aims to reflect the production environment in real practice, including quantity discounts from different suppliers, different sizes of vehicles with different emissions, and production-capacity constraints.Through the proposed models, operation schedules are devised to coordinate the replenishment, transportation, and production activities and minimize the total cost (i.e., the sum of replenishment, transportation, production, and emission cost) over a given planning horizon, while the customer demand, vehicle travel length, and loading constraints, plant production, and inventory and backlogging constraints are all satisfied.Both a mixed integer programming (MIP) model and a particle swarm optimization (PSO) model are constructed to solve this sustainable supply chain management problem to minimize the total cost in the system during a planning horizon.When the problem scale is small, both the MIP and the PSO can lead to optimal solutions within a reasonable time.However, when the problem becomes complicated and reflects real application more, the MIP model may require a long computational time and may still not obtain the optimal solution.On the other hand, the PSO can be an efficient model for obtaining a near-optimal solution.Based on the results of the models, managers can devise their supply chain plans, including the purchase, production, and transportation decisions, efficiently.
In this research, a lot of aspects are presumed to be known and fixed.For example, the demand, lead time, quantity discount information, and transportation issues are known.However, in practice, a supply chain environment usually consists of uncertain demand, variable lead time, different ranks of orders, different quantity discounts under different situations, uncertain transportation conditions, etc.Additional cost and time issues can be considered to meet the practical requirements; these could include for example, the lead-time constraint for part procurement, crashing cost for shipment, limited storage capacity, and fixed setup cost under different production modes.Thus, some of the assumptions that were made in the study can be relaxed to consider these issues and better represent a real sustainable supply chain setting.The problem can be constructed as a multi-objective problem to consider multiple objectives, such as decreasing the monetary cost, increasing the sustainability, and increasing the social responsibility of the system.Both a multi-objective programming model and a PSO model may be constructed to solve the problem.In addition, other heuristics, such as a genetic algorithm, ant colony system, and artificial immune system may be applied, and a comparison of the methods can be performed to examine which method is more suitable for solving the problem.

= 3 Figure 1 .
Figure 1.Vehicle routing under Case I using the MIP and the PSO models.

Figure 1 .
Figure 1.Vehicle routing under Case I using the MIP and the PSO models.

Figure 3 .
Figure 3. Vehicle routing under Case II using the PSO model.Figure 3. Vehicle routing under Case II using the PSO model.

Figure 4 .
Figure 4. Execution result of Case II using the PSO model.

Figure 4 .
Figure 4. Execution result of Case II using the PSO model.

= 9 Figure 5 .
Figure 5. Vehicle routing under Case III using the PSO model.

Figure 5 .
Figure 5. Vehicle routing under Case III using the PSO model.

Table 1 .
Comparison of relevant works.

Table 2 .
Ordering cost (o vr ) from each supplier.

Table 3 .
Unit cost under quantity discount from each supplier.

Table 4 .
Inventory-holding cost for each part and finished product.

Table 5 .
Data for each vehicle.

Table 6 .
Distance (u ij ) and transportation cost (c ij ) for vehicles.

Table 7 .
Production cost under different production modes.

Table 8 .
Part requirements for finished good.

Table 9 .
Demand of finished good (d tg ).

Table 10 .
Relevant results in each period under Case I using the mixed integer programming (MIP) and the particle swarm optimization (PSO) models.

Table 11 .
Relevant results in each period under Case II using the MIP model.

Table 12 .
Relevant results in each period under Case II using the PSO.

Table 13 .
Relevant results in each period under Case III using the PSO model.

Table 14 .
Effects of parameter changes on the system under Case II.
Unit purchase cost of part r under quantity discount bracket x from supplier v in period t l vrx Maximum quantity of part r under quantity discount bracket x from supplier v z gs Maximum accumulated quantity of finished good g that can be produced from production mode 1 to s k e Maximum travelling length of vehicle e u ij Distance from shipment point i to shipment point j w e Maximum loading size of vehicle e ρ gr Units of material r required to produce product g.c ij Transportation cost from shipment point i to shipment point j FC e Fixed cost of vehicle e per trip P(Q tvr ) Unit purchase cost of part r from supplier v in period t Q tvr Quantity of part r purchased from supplier v in period t S tr Total quantity of part r purchased in period t P ϑ tgs Unit production cost of finished good g under production mode s in period t.Depending on the quantity manufactured, the unit production cost will be based on the production mode.P ϑ tgs Production quantity of finished good g under production mode s in period t S tg Total quantity of finished good g produced in period t Y ti Purchase size from shipment point i in period t Backlogging of finished good g in period t α tvr Binary variable, 1 indicates that an order of part r from supplier v in period t is placed, and 0 indicates that no order is placed β tvrx Binary variable, 1 indicates that an order of part r under quantity discount bracket x from supplier v in period t is placed, and 0 indicates that no order is placed ϕ e tij Binary variable, 1 indicates that vehicle e travels from shipment point i to shipment point j in period t, and 0 indicates that no travel is incurred φ e ti Binary variable, 1 indicates that vehicle e travels from shipment point i in period t, and 0 indicates that no travel is incurred ω tgs Binary variable, 1 indicates that finished good g is manufactured under production mode s in period t, and 0 indicates that no product is manufactured tij Loading size of vehicle e from shipment point i to shipment point j in period t F + tr Ending inventory of part r in period t FG + tg Ending inventory of finished good g in period t F − tr Backlogging of part r in period t FG − tg