# A Multiobjective Optimization Model in Automotive Supply Chain Networks

Authors: Abdolhossein Sadrnia; Napsiah Ismail; Norzima Zulkifli; M. K. A. Ariffin; Hossein Nezamabadi-pour; Hamed Mirabi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

In the new decade, green investment decisions are attracting more interest in design supply chains due to the hidden economic benefits and environmental legislative barriers. In this paper, a supply chain network design problem with both economic and environmental concerns is presented. Therefore, a multiobjective optimization model that captures the trade-off between the total logistics cost and CO2 emissions is proposed. With regard to the complexity of logistic networks, a new multiobjective swarm intelligence algorithm known as a multiobjective Gravitational search algorithm (MOGSA) has been implemented for solving the proposed mathematical model. To evaluate the effectiveness of the model, a comprehensive set of numerical experiments is explained. The results obtained show that the proposed model can be applied as an effective tool in strategic planning for optimizing cost and CO2 emissions in an environmentally friendly automotive supply chain.


---

## Body

## 1. Introduction

In conventional supply chain management, the emphasis of the supply chain network is usually a single objective, such as profit maximization or cost minimization. Despite waterway pollution and the massive loss of biodiversity, global warming and carbon dioxide (CO2) emissions are considered to be the types of environmental impacts that have to be considered with regard to the coordinated activity of organizations in a supply chain. Additionally, attention must be given to the greenhouse gases that are increasingly facilitated by global warming, and which are the result of human activities since the industrial revolution. The most significant anthropogenic greenhouse gas is CO2. The situation is that global CO2 emissions have increased by about 80% during the period of 1970–2004, mainly as a result of using fossil fuels [1].The distances between facilities in a distribution network have grown dramatically due to the globalisation of supply chains, which leads to increased vehicle emissions as part of transportation. In the US, approximately 33% of CO2 emissions and 73% of carbon monoxide (CO) are based on the distribution of products and transportation [2]. Consequently, there is the need to effectively and efficiently design environmentally friendly supply chain networks, to both improve environmental conditions as well as the economic dimension of an organisation. Elhedhli and Merrick [3] expressed network design as a logical place to start when looking towards greening a supply chain design. Usually, the planning of a logistics network requires making decisions regarding the followings [4]:(i)
the number, location, capacity, and technology of manufacturing plants, warehouses, and retailers;(ii)
selection of suppliers;(iii)
assignment of product ranges to manufacturing plants and warehouses;(iv)
selection of distribution channels and transportation modes;(v)
flow of raw materials through the network.Over the last decade, besides the economic concerns, environmental issues have interested many scholars when they design supply chain networks. The network design problem is one of the most comprehensive strategic decision problems that should be optimised for long-term efficient operations throughout the supply chain. The problem includes a wide range of mathematical modelling ranging from linear deterministic models [5–7] to complex nonlinear stochastic models [8]. In the literature, several attempts have been made to solve the network design problems of supply networks, and these studies have been surveyed by Beamon and Fernandes [9], Chopra [10], and  Peidro et al. [11]. Scholars have attempted to develop the classical quantitative model by incorporating various factors such as transportation modes [12, 13], vehicle speed [14] and greenhouse gasses [15, 16].Despite the volume of work undertaken in respect to the supply chain network design problem, the modelling is usually done as a single objective problem [17, 18] while the “design” in real cases typically involves compromises and trade-offs between different incompatible objectives. For example, the optimum design for a supply chain network should not only concentrate on total cost minimisation, but it also has to satisfy other objectives such as delivery flexibility [19], minimising waste [20], and minimising CO2 emissions [21].Satisfying environmental issues and minimizing the total variable and fixed costs causes the logistics network design problem to be an NP-hard problem [22, 23]. Solving this problem on a large scale by exact algorithms is very time consuming and sometimes in the practical case is impossible. Therefore, many heuristics and metaheuristics have been developed to obtain optimal solutions for these kinds of problems [24]. In this regard, during the last decade, the implementation of meta-heuristics instead of exact and exhaustive methods has been growing. Approaches such as bioinspired algorithms such as ant colony optimization (ACO) [25–27], particle swarm optimization (PSO) [28], artificial weed colonies [29], and also evolutionary algorithms such as genetic algorithms (GA) have been used to solve a variety of multi-objective problems in production and supply chain management specially in complex environment [29, 30].Recently, the consideration of supply network design with multi-objective optimisation (generally incompatible objectives) is a new trend worthy of study. The solution to the above problem is not unique, but a set of Pareto optimal points. A Pareto optimal solution is the best solution vector out of several numbers of solution vectors that could be achieved without disadvantaging other objectives. Each solution set is known as nondominated solution set. An array containing all the Pareto optimal solutions of a multi-objective problem is called the Pareto optimal set. Thus, instead of being a unique solution to the problem, the solution of a multi-objective problem is a possibly infinite set of Pareto points [31]. Compared to studies working with a single objective, a multi-objective approach is more reasonable and more practical in terms of actual applications such that it has been considered by different researchers in the literature. For example, a multi-objective programming model has been proposed [32] to minimise the total cost of plant and its distribution centres and also to maximise customer delivery time. In addition, Pishvaee et al. [33] proposed a biobjective model that involved a trade-off between the total costs and the responsiveness of the supply chain network. Alçada-Almeida et al. [34] proposed a multi-objective programming method to determine the most appropriate location allocation of hazardous material incineration facilities and to make a mutual balance between the societal, economic, and environmental impacts. Paksoya et al. [35] conducted another related study in which the green impact on a close-looped logistics supply chain network was investigated. The authors attempted to decrease and prevent greater CO2 gas emissions and also to encourage customers to use recyclable goods. They have presented different transportation choices between echelons according to CO2 emissions. Another related study was conducted by Wang et al. [36], who presented a mathematical model to minimize CO2 emissions and cost to design a green supply chain. They defined environmental protection level variable and CO2 emissions between supply chain facilities.Among recently published studies, gaps can be found in adaptation research, and there is a deficit of studies in most developing regions. This is especially true for determining alternative solutions for decision makers that can make the trade-off between environmental impacts and economic activities based on Pareto optimality that is known as quite a new problem in supply chain network design [37]. Among the related researches that have been presented on the mathematical model in green supply chain, Wang et al. [36] used the Pareto set concept to make a solution set, but they neglected the truck type for transporting material and products in modeling. Another multi-objective green supply chain modeling has been done by Paksoya et al. [35] which focuses on optimizing CO2 and cost, but it neglects to present the Pareto optimality set. In this research, a multi-objective mixed-integer formulation is provided for the supply chain network design problem. The multi-objective model explicitly considers environmental concerns by introducing the three groups of tier-1 suppliers as the new concept of supply chain thinking, which is compatible with the automotive supply chain logistics. Furthermore, in this research, choosing truck type and facility location among the potential locations have been considered to find Pareto optimal set solutions and trade-off cost and CO2 emission in the logistics network. Another contribution of the research is implementation of the multi-objective gravitational search algorithm (MOGSA) as a new swarm algorithm for solving complex logistical mathematical models to find a set of Pareto optimal solutions. The result obtained can be easily applied to the decision support systems, which the industry needs. Furthermore, the last contribution of this research is the evaluation of the proposed model by using a real case study for one of the largest car manufacturers in the Middle East.The rest of this research is organised as follows. In the next section, the problem is presented in detail including mathematical modelling with two objectives: cost and CO2 minimisation. In Section 3, GSA and MOGSA theories are introduced. For a real case the objective functions are solved by using MOGSA in MATLAB software in Section 4. Finally, the paper is concluded in Section 5, and some interesting possible future research is also introduced.


## 2. Problem Modeling

In this section, the mathematical formulation of the design model for a supply chain network will be presented in detail. In automotive industries, the general structure of the proposed supply chain logistic network is illustrated in Figure1. In the forward direction, the suppliers are responsible for providing all the parts needed by the car manufacturing facilities. Suppliers have been grouped into three categories as tier-1 suppliers. The first group of tier-1 suppliers is responsible for supplying metal-based parts such as wheels, the main body, and springs. The second group of suppliers is responsible for supplying electric parts such as wiring harnesses, instrument panels, and engine control units. The last tier-1 group supplies plastics- and latex-based parts such as tires, interior parts, and dashboard panels. The new products are conveyed from plants to retailers via distribution centres to meet customer demand. All assumptions and limitations, sets, parameters, and variables are defined and follow a model-detailed description.Figure 1
General structure of automotive supply chain logistic network.


### 2.1. Assumptions and Limitations

(i)
The model is deterministic so that demand is fixed and known and must be satisfied.
(ii)
The model selects plants, warehouses, and distributions from a given set of potential locations. Moreover, the model also selects the supplier and retailer from a set of fixed locations.
(iii)
Suppliers, plants, warehouse, and distributors have a capacity limit for the product.
(iv)
Cost is independent of product load but is varied for any kind of transportation vehicles such as the truck road classes.
(v)
Carbon emission is dependent on product load (capacity of truck type) and also is varied for any kind of transportation vehicles.
(vi)
All truck types in transportation are turned off when unloading product.
(vii)
The flow is only allowed to be transferred between two sequential echelons.(1)Indices and SetsI={1,…,Nv}: set of fixed location of metal-parts suppliers’ tier-1.J={1,…,Ng}: set of fixed location of electric-parts suppliers’ tier-1.N={1,…,Nd}: set of fixed location of plastic-parts suppliers’ tier-1.K={1,…,Nf}: set of potential available for plants.L={1,…,Nw}: set of potential available for warehouses and distributors.M={1,…,Nr}: set of fixed location of retailers.(2)Parametersεt: maximum capacity in truck type t (tone), t∈T.τt: transportation cost mode t per km ($/km), t∈T.φt: carbon dioxide emission by truck type t (gram/km), t∈T.∂ikvf: distance between metal-parts suppliers’ tier-1 i to plant k (km) that i∈I, k∈K.∂jkgf: distance between electric-parts suppliers’ tier-1 j to plant k (km) that i∈I, k∈K.∂nkdf: distance between plastic-parts suppliers’ tier-1 n to plant k (km) that n∈N, k∈K.∂klfw: distance between plant k to distributor l (km) that k∈K, l∈L.∂lmwr: distance between warehouse and distributor l to retailer m (km) that l∈L, k∈K.ωp: weight of product (tone).ρm: rate of metal parts contribution in the product (0-1).ρe: rate of electrical parts contribution in the product (0-1).ρp: rate of plastic parts contribution in the product (0-1).ηm: product demand retailer m.μiv: maximum capacity metal-parts suppliers’ tier-1 i (tone) i∈I.μjg: maximum capacity electric-parts suppliers’ tier-1 j∈J.μnd: maximum capacity plastic-parts suppliers’ tier-1 n (tone) n∈N.μkf: maximum product capacity for plant k that k∈K.μlw: maximum product capacity for warehouse and distributor l that l∈L.γkf: unit Production and material handling cost of product in plant k ($) that k∈K.γlw: unit operation cost of product in distributor l ($) that l∈L.γmr: unit operation cost of product in retailer m ($) that m∈M.αkf: fixed cost of opening plant k ($) that k∈K.αlw: fixed cost of opening distribution center l ($) that l∈L.(3)Decision VariablesXikt1:quantity of product produced in plant k using metal parts of metal-parts suppliers’ tier-1 i by truck type t. k∈K, i∈I, t∈T.Xjkt2:quantity of product produced in plant k using electrical parts of electric-parts suppliers’ tier-1 j by truck type t. k∈K, j∈J, t∈T.Xnkt3:quantity of product produced in plant k using plastic parts of plastic-parts suppliers’ tier-1 n by truck type t. n∈N, k∈K, t∈T.Xklt4:quantity of product shipped from plant k todistribution center l by truck type t. k∈K, l∈L, t∈T.Xlmt5:quantity of product shipped fromdistribution center l to retailer m by truck type t. l∈L, m∈M, t∈T.
(1)Yk1={1iffactorykisopened,k∈K,0otherwise,Yl2={1ifwarehouselisopened,l∈L,0otherwise.Objective Functions. With regards to the above description, the first objective function that minimizes the total cost of supply chain network is as follows:
(2)TotalCost=TransportationCost+OperationCost+InitialFacilityCost.In objective function (3), the first five terms define the total transportation cost of the products between each pair suppliers, plants, and distributors. Next two terms are referred to operation cost in plants and distributors. Fixed setup cost for opening the facilities (factories, distributor centers, resp.) can be defined by the last two terms:
(3)TotalCost=∑t∑k∑i(ωp·ρm·xikt1εt)·∂ikvf·τt+∑t∑k∑j(ωp·ρe·xjkt2εt)·∂jkgf·τt+∑t∑k∑n(ωp·ρp·xnkt3εt)·∂nkdf·τt+∑t∑l∑k(ωp·xklt4εt)·∂klfw·τt+∑t∑m∑l(ωp·xlmt5εt)·∂lmwr·τt+∑t∑l∑kxklt4·γkf+∑t∑m∑lxlmt5·γlw+∑k∈Kαkf·Yk1+∑l∈Lαlw·Yl2.The second objective function that should be minimized is total CO2 emission in supply chain network due to transportation between all facilities:
(4)CarbonEmission=∑t∑k∑i(ωp·ρm·xikt1εt)·∂ikvf·φt+∑t∑k∑j(ωp·ρn·xjkt2εt)·∂jkgf·φt+∑t∑k∑n(ωp·ρp·xnkt3εt)·∂nkdf·φt+∑t∑l∑k(ωp·xklt4εt)·∂klfw·φt+∑t∑m∑l(ωp·xlmt5εt)·∂lmwr·φt.Constrains. This section is a representation of the constraints of the proposed model.
(5)∑t∑lxlmt5=ηm∀m∈M,(6)ωp·ρm·∑u∑t∑ixikt1+ωp·ρe·∑t∑jxjkt2+ωp·ρp∑t∑nxnkt3=∑t∑lωp·xklt4∀k∈K,(7)∑t∑kxlktu4=∑t∑mxlmt5∀l∈L,(8)ρm·∑t∑kxikt1≤μiv∀i∈I,(9)ρe·∑t∑kxjkt2≤μjd∀j∈J,(10)ρp·∑t∑kxnkt3≤μjd∀n∈N,(11)∑t∑lxkltu4≤μkf·Yk1∀k∈K,(12)∑t∑mxlmtu5≤μlw·Yl2∀l∈L,(13)Xikt1,Xjkt2,Xnkt3,Xklt4,Xlmt5≥0∀i,j,k,l,m,t,(14)Yk1,Yl2∈{0,1}∀k,l.Equation (5) states that, for each product, the flow exiting the distribution centers must satisfy the demand of all customers. Equation (6) shows that, for each product, the sum of the flow entering each plant from all suppliers is equal to the flow leaving the facility to the distribution centers. Equation (7) insures that the sum of the product flow entering each distribution centre from all plants is equal to the flow leaving this distribution centre to retailers. Equation (8) ensures that, for each product, the sum of the flow exiting each metal-parts supplier to all plants does not exceed the supplier’s capacity. Equation (9) ensures that, for each product, the sum of the flow exiting each electric-parts supplier to all plants does not exceed the supplier’s capacity. Equation (10) ensures that, for each product, the sum of the flow exiting each electric-parts supplier to all plants does not exceed the supplier’s capacity. Equation (11) states that the sum of the flow exiting each plant to all distribution counters does not exceed the plant’s capacity. And finally, (12) shows that the sum of the flow exiting from each distribution centre to all customers does not exceed the distributor’s capacity. Constraints (13) and (14) impose binary and nonnegativity restrictions on the corresponding decision variables, respectively.


## 2.1. Assumptions and Limitations

(i)
The model is deterministic so that demand is fixed and known and must be satisfied.
(ii)
The model selects plants, warehouses, and distributions from a given set of potential locations. Moreover, the model also selects the supplier and retailer from a set of fixed locations.
(iii)
Suppliers, plants, warehouse, and distributors have a capacity limit for the product.
(iv)
Cost is independent of product load but is varied for any kind of transportation vehicles such as the truck road classes.
(v)
Carbon emission is dependent on product load (capacity of truck type) and also is varied for any kind of transportation vehicles.
(vi)
All truck types in transportation are turned off when unloading product.
(vii)
The flow is only allowed to be transferred between two sequential echelons.(1)Indices and SetsI={1,…,Nv}: set of fixed location of metal-parts suppliers’ tier-1.J={1,…,Ng}: set of fixed location of electric-parts suppliers’ tier-1.N={1,…,Nd}: set of fixed location of plastic-parts suppliers’ tier-1.K={1,…,Nf}: set of potential available for plants.L={1,…,Nw}: set of potential available for warehouses and distributors.M={1,…,Nr}: set of fixed location of retailers.(2)Parametersεt: maximum capacity in truck type t (tone), t∈T.τt: transportation cost mode t per km ($/km), t∈T.φt: carbon dioxide emission by truck type t (gram/km), t∈T.∂ikvf: distance between metal-parts suppliers’ tier-1 i to plant k (km) that i∈I, k∈K.∂jkgf: distance between electric-parts suppliers’ tier-1 j to plant k (km) that i∈I, k∈K.∂nkdf: distance between plastic-parts suppliers’ tier-1 n to plant k (km) that n∈N, k∈K.∂klfw: distance between plant k to distributor l (km) that k∈K, l∈L.∂lmwr: distance between warehouse and distributor l to retailer m (km) that l∈L, k∈K.ωp: weight of product (tone).ρm: rate of metal parts contribution in the product (0-1).ρe: rate of electrical parts contribution in the product (0-1).ρp: rate of plastic parts contribution in the product (0-1).ηm: product demand retailer m.μiv: maximum capacity metal-parts suppliers’ tier-1 i (tone) i∈I.μjg: maximum capacity electric-parts suppliers’ tier-1 j∈J.μnd: maximum capacity plastic-parts suppliers’ tier-1 n (tone) n∈N.μkf: maximum product capacity for plant k that k∈K.μlw: maximum product capacity for warehouse and distributor l that l∈L.γkf: unit Production and material handling cost of product in plant k ($) that k∈K.γlw: unit operation cost of product in distributor l ($) that l∈L.γmr: unit operation cost of product in retailer m ($) that m∈M.αkf: fixed cost of opening plant k ($) that k∈K.αlw: fixed cost of opening distribution center l ($) that l∈L.(3)Decision VariablesXikt1:quantity of product produced in plant k using metal parts of metal-parts suppliers’ tier-1 i by truck type t. k∈K, i∈I, t∈T.Xjkt2:quantity of product produced in plant k using electrical parts of electric-parts suppliers’ tier-1 j by truck type t. k∈K, j∈J, t∈T.Xnkt3:quantity of product produced in plant k using plastic parts of plastic-parts suppliers’ tier-1 n by truck type t. n∈N, k∈K, t∈T.Xklt4:quantity of product shipped from plant k todistribution center l by truck type t. k∈K, l∈L, t∈T.Xlmt5:quantity of product shipped fromdistribution center l to retailer m by truck type t. l∈L, m∈M, t∈T.
(1)Yk1={1iffactorykisopened,k∈K,0otherwise,Yl2={1ifwarehouselisopened,l∈L,0otherwise.Objective Functions. With regards to the above description, the first objective function that minimizes the total cost of supply chain network is as follows:
(2)TotalCost=TransportationCost+OperationCost+InitialFacilityCost.In objective function (3), the first five terms define the total transportation cost of the products between each pair suppliers, plants, and distributors. Next two terms are referred to operation cost in plants and distributors. Fixed setup cost for opening the facilities (factories, distributor centers, resp.) can be defined by the last two terms:
(3)TotalCost=∑t∑k∑i(ωp·ρm·xikt1εt)·∂ikvf·τt+∑t∑k∑j(ωp·ρe·xjkt2εt)·∂jkgf·τt+∑t∑k∑n(ωp·ρp·xnkt3εt)·∂nkdf·τt+∑t∑l∑k(ωp·xklt4εt)·∂klfw·τt+∑t∑m∑l(ωp·xlmt5εt)·∂lmwr·τt+∑t∑l∑kxklt4·γkf+∑t∑m∑lxlmt5·γlw+∑k∈Kαkf·Yk1+∑l∈Lαlw·Yl2.The second objective function that should be minimized is total CO2 emission in supply chain network due to transportation between all facilities:
(4)CarbonEmission=∑t∑k∑i(ωp·ρm·xikt1εt)·∂ikvf·φt+∑t∑k∑j(ωp·ρn·xjkt2εt)·∂jkgf·φt+∑t∑k∑n(ωp·ρp·xnkt3εt)·∂nkdf·φt+∑t∑l∑k(ωp·xklt4εt)·∂klfw·φt+∑t∑m∑l(ωp·xlmt5εt)·∂lmwr·φt.Constrains. This section is a representation of the constraints of the proposed model.
(5)∑t∑lxlmt5=ηm∀m∈M,(6)ωp·ρm·∑u∑t∑ixikt1+ωp·ρe·∑t∑jxjkt2+ωp·ρp∑t∑nxnkt3=∑t∑lωp·xklt4∀k∈K,(7)∑t∑kxlktu4=∑t∑mxlmt5∀l∈L,(8)ρm·∑t∑kxikt1≤μiv∀i∈I,(9)ρe·∑t∑kxjkt2≤μjd∀j∈J,(10)ρp·∑t∑kxnkt3≤μjd∀n∈N,(11)∑t∑lxkltu4≤μkf·Yk1∀k∈K,(12)∑t∑mxlmtu5≤μlw·Yl2∀l∈L,(13)Xikt1,Xjkt2,Xnkt3,Xklt4,Xlmt5≥0∀i,j,k,l,m,t,(14)Yk1,Yl2∈{0,1}∀k,l.Equation (5) states that, for each product, the flow exiting the distribution centers must satisfy the demand of all customers. Equation (6) shows that, for each product, the sum of the flow entering each plant from all suppliers is equal to the flow leaving the facility to the distribution centers. Equation (7) insures that the sum of the product flow entering each distribution centre from all plants is equal to the flow leaving this distribution centre to retailers. Equation (8) ensures that, for each product, the sum of the flow exiting each metal-parts supplier to all plants does not exceed the supplier’s capacity. Equation (9) ensures that, for each product, the sum of the flow exiting each electric-parts supplier to all plants does not exceed the supplier’s capacity. Equation (10) ensures that, for each product, the sum of the flow exiting each electric-parts supplier to all plants does not exceed the supplier’s capacity. Equation (11) states that the sum of the flow exiting each plant to all distribution counters does not exceed the plant’s capacity. And finally, (12) shows that the sum of the flow exiting from each distribution centre to all customers does not exceed the distributor’s capacity. Constraints (13) and (14) impose binary and nonnegativity restrictions on the corresponding decision variables, respectively.


## 3. Solving Method

In this section, we first introduce the Gravitational Search Algorithm (GSA) and then explain the extended GSA, which is used to solve multi-objective optimisation, known as the multi-objective gravitational search algorithm (MOGSA).


### 3.1. GSA

One of the novel meta-heuristic stochastic optimisation algorithms inspired by the law of gravity and mass interactions is called the Gravitational Search Algorithm (GSA) [38]. In GSA, the individuals are called agents and are a collection of objects that interact with each other based on Newtonian gravity and the laws of motion. The agents share information using gravitational force to guide the search towards the best location in the search space. The high performance and the global search ability of GSA in solving various nonlinear functions infers from the results of experiments undertaken previously [38–42]. At GSA, searcher agents are a collection of objects. All these objects attract one another by gravitational force, and this force causes a global movement of all objects toward the objects with heavier masses, as Figure 2 shows. Hence, objects cooperate using a direct form of communication, through gravitational force. The heavy objects that correspond to good solutions move more slowly than the lighter objects, which guarantee the exploitation ability of the algorithm.Figure 2
Each object accelerates toward the result force that attracts it from the other objects.The position of the objects corresponds to a solution of the problem, and its mass is determined using a map of fitness function. By changing the velocities over time, the agents are likely to move towards the global optima. In this algorithm, there are objects. The position of theith agent is defined by
(15)Xi=(xi1,…,xid,…,xin),i=1,2,…,N,
where xid presents the position of the ith agent in the dth dimension. Based on the Newton gravitation theory, the force acting on the ith object from the jth object is
(16)Fijd(t)=G(t)Mi(t)×Mj(t)Rij(t)+ε(xjd(t)-xid(t)),(17)G(t)=G0*exp(-atT),
where Mi is the mass of object i, Mj is the mass of object j, and  G(t) is the gravitational constant at time t. Notably, the gravitational constant G(t) is important in determining the performance of GSA; G(t) is defined as a decreasing function (17) of time, which is set to G0 at the beginning and decreased exponentially to zero as time progresses.To compute the acceleration of an agent, total forces from a set of heavier masses that apply on it should be considered based on law of gravity (18) which is followed by calculation of agent acceleration using law of motion (19). Afterward, the next velocity of an agent is calculated as a fraction of its current velocity added to its acceleration (20). Then, its position could be calculated using
(18)Fid(t)=∑j∈kbest,j≠irandjG(t)Mj(t)Mi(t)Rij(t)+ε(xjd(t)-xid(t)),(19)aid(t)=Fid(t)Mi(t)=∑j∈kbest,j≠irandjG(t)Mj(t)Rij(t)+εrandj××(xjd(t)-xid(t)),(20)vid(t+1)=randi×vid(t)+aid(t),(21)xid(t+1)=xid(t)+vid(t+1),
where randi and randj are two uniform randoms in the interval [0,1], ε is a small value, and Rij(t) is the Euclidian distance between two agents i and j defined as Rij(t)=∥Xi(t),Xj(t)∥2. kbest is the set of first K agents with the best fitness value and biggest mass. kbest is a function of time, initialized to K0 at the beginning and decreasing with time. Here, K0 is set to N (total number of agents) and is decreased linearly to 1.Based on [39], the mass of each agent is calculated after computing current population’s fitness as follows:
(22)qi(t)=fiti(t)-worst(t)best(t)-worst(t),Mi(t)=qi(t)∑j=1sqj(t),
where Mi(t) and fiti(t) represent the mass and the fitness value of the agent i at t, and worst(t) and best(t) are defined as follows (for a minimization problem):
(23)best(t)=minj∈{1,..,s}fitj(t),worst(t)=maxj∈{1,..,s}fitj(t).
Figure 3 illustrates the flow chart of the GSA.Figure 3
General flowchart of GSA.


### 3.2. MOGSA

Suppose thatX=(x1,x2,…,xn) is the vector of decision variable, fi:Rn→R, i=1,…,k are the objective functions and gi,hi:Rn→R, i=1,…,m, j=1,…,p are the constraints of the problem. Pareto domain-based algorithm, MOGSA, that is used in this paper, presents this problem and definition as follows:(i)
Minimizef(X)=[f1(X),f2(X),…,fk(X)](ii)
Subject togi(X)≤0, i=1,2,…,m(iii)
hj(X)=0, j=1,2,…,p.Definition 1.
Given two vectorsX,Y∈Rk, we say that X≤Y if xi≤yi, i=1,…,k and that X dominates Y (or presented  by X≺Y) if X≤Y and X≠Y.Definition 2.
We sayX∈χ⊂Rk is nondominated with respect to χ if another X′∈χ  does not exist such that f(X′)≺f(X).Selected non-dominated solutions are stored in an external archive. Archive updates in each iteration of the algorithm. When the running algorithm is terminated, external archive presents final results.Recently, the multi-objective gravitational search algorithm (MOGSA) [43–45] has been developed to solve multi-objective problems by using the Uniform Mutation Operator [46] and an Elitist Policy [47] as well-known operators in genetic algorithm (GA) into GSA. MOGSA also uses an external archive to reserve the non-dominated solutions and update it, the same as Simple Multi-Objective PSO (SMOPSO) [47]. The first operator selects one of the object’s dimensions randomly and adds a signed random term to it. The total value should be in the permitted range otherwise it will be truncated to its nearest feasible value. The non-dominated solutions found are stored by the Elitist Policy in an archive that has a grid structure. The grid structure is created as follows: each dimension in the objective space is divided into 2ni equal divisions where i denotes dimension index and hence for a kth objective optimisation problem there will be 2(∑i=1kni) different segments. As long as the archive is not full, new non-dominated solutions are added to it. As soon as the number of elements in the archive meets its limit, one of the elements should be omitted from it. The strategy to remove an element from the archive is to find the most crowded hypercube and randomly select one of the elements in it to extract. For this purpose, objective functions and solution space are divided into hyper-rectangles. When archive overflow is reached, one particle from the most crowded hyper-rectangle is randomly selected and removed. In MOGSA, the distance between each particle to its nearest neighbor is calculated, and its mass is updated. Distribution of archived particles is done by using the niching technique. In this algorithm, the archived particles apply gravitational force. The list below shows the pseudocode of MOGSA that is similar to Hassanzadeh and Rouhani [43]. (1)
Set the population size,G0, a and ε.(2)
Initialize population of solutions.(3)
Evaluate population (with regards to objective functions: CO2 and Cost).(4)
Update gravitational constant.(5)
Calculate the masses.(6)
Calculate force and acceleration between objects.(7)
Update velocity and position per each object.(8)
Select set of best solutions (Pareto front)  //Update external archive//.(9)
Generate new population.(10)
If the fitness function is satisfied ort=Tmax go to 4 otherwise continued.(11)
Figure set of Pareto solutions and end.


## 3.1. GSA

One of the novel meta-heuristic stochastic optimisation algorithms inspired by the law of gravity and mass interactions is called the Gravitational Search Algorithm (GSA) [38]. In GSA, the individuals are called agents and are a collection of objects that interact with each other based on Newtonian gravity and the laws of motion. The agents share information using gravitational force to guide the search towards the best location in the search space. The high performance and the global search ability of GSA in solving various nonlinear functions infers from the results of experiments undertaken previously [38–42]. At GSA, searcher agents are a collection of objects. All these objects attract one another by gravitational force, and this force causes a global movement of all objects toward the objects with heavier masses, as Figure 2 shows. Hence, objects cooperate using a direct form of communication, through gravitational force. The heavy objects that correspond to good solutions move more slowly than the lighter objects, which guarantee the exploitation ability of the algorithm.Figure 2
Each object accelerates toward the result force that attracts it from the other objects.The position of the objects corresponds to a solution of the problem, and its mass is determined using a map of fitness function. By changing the velocities over time, the agents are likely to move towards the global optima. In this algorithm, there are objects. The position of theith agent is defined by
(15)Xi=(xi1,…,xid,…,xin),i=1,2,…,N,
where xid presents the position of the ith agent in the dth dimension. Based on the Newton gravitation theory, the force acting on the ith object from the jth object is
(16)Fijd(t)=G(t)Mi(t)×Mj(t)Rij(t)+ε(xjd(t)-xid(t)),(17)G(t)=G0*exp(-atT),
where Mi is the mass of object i, Mj is the mass of object j, and  G(t) is the gravitational constant at time t. Notably, the gravitational constant G(t) is important in determining the performance of GSA; G(t) is defined as a decreasing function (17) of time, which is set to G0 at the beginning and decreased exponentially to zero as time progresses.To compute the acceleration of an agent, total forces from a set of heavier masses that apply on it should be considered based on law of gravity (18) which is followed by calculation of agent acceleration using law of motion (19). Afterward, the next velocity of an agent is calculated as a fraction of its current velocity added to its acceleration (20). Then, its position could be calculated using
(18)Fid(t)=∑j∈kbest,j≠irandjG(t)Mj(t)Mi(t)Rij(t)+ε(xjd(t)-xid(t)),(19)aid(t)=Fid(t)Mi(t)=∑j∈kbest,j≠irandjG(t)Mj(t)Rij(t)+εrandj××(xjd(t)-xid(t)),(20)vid(t+1)=randi×vid(t)+aid(t),(21)xid(t+1)=xid(t)+vid(t+1),
where randi and randj are two uniform randoms in the interval [0,1], ε is a small value, and Rij(t) is the Euclidian distance between two agents i and j defined as Rij(t)=∥Xi(t),Xj(t)∥2. kbest is the set of first K agents with the best fitness value and biggest mass. kbest is a function of time, initialized to K0 at the beginning and decreasing with time. Here, K0 is set to N (total number of agents) and is decreased linearly to 1.Based on [39], the mass of each agent is calculated after computing current population’s fitness as follows:
(22)qi(t)=fiti(t)-worst(t)best(t)-worst(t),Mi(t)=qi(t)∑j=1sqj(t),
where Mi(t) and fiti(t) represent the mass and the fitness value of the agent i at t, and worst(t) and best(t) are defined as follows (for a minimization problem):
(23)best(t)=minj∈{1,..,s}fitj(t),worst(t)=maxj∈{1,..,s}fitj(t).
Figure 3 illustrates the flow chart of the GSA.Figure 3
General flowchart of GSA.


## 3.2. MOGSA

Suppose thatX=(x1,x2,…,xn) is the vector of decision variable, fi:Rn→R, i=1,…,k are the objective functions and gi,hi:Rn→R, i=1,…,m, j=1,…,p are the constraints of the problem. Pareto domain-based algorithm, MOGSA, that is used in this paper, presents this problem and definition as follows:(i)
Minimizef(X)=[f1(X),f2(X),…,fk(X)](ii)
Subject togi(X)≤0, i=1,2,…,m(iii)
hj(X)=0, j=1,2,…,p.Definition 1.
Given two vectorsX,Y∈Rk, we say that X≤Y if xi≤yi, i=1,…,k and that X dominates Y (or presented  by X≺Y) if X≤Y and X≠Y.Definition 2.
We sayX∈χ⊂Rk is nondominated with respect to χ if another X′∈χ  does not exist such that f(X′)≺f(X).Selected non-dominated solutions are stored in an external archive. Archive updates in each iteration of the algorithm. When the running algorithm is terminated, external archive presents final results.Recently, the multi-objective gravitational search algorithm (MOGSA) [43–45] has been developed to solve multi-objective problems by using the Uniform Mutation Operator [46] and an Elitist Policy [47] as well-known operators in genetic algorithm (GA) into GSA. MOGSA also uses an external archive to reserve the non-dominated solutions and update it, the same as Simple Multi-Objective PSO (SMOPSO) [47]. The first operator selects one of the object’s dimensions randomly and adds a signed random term to it. The total value should be in the permitted range otherwise it will be truncated to its nearest feasible value. The non-dominated solutions found are stored by the Elitist Policy in an archive that has a grid structure. The grid structure is created as follows: each dimension in the objective space is divided into 2ni equal divisions where i denotes dimension index and hence for a kth objective optimisation problem there will be 2(∑i=1kni) different segments. As long as the archive is not full, new non-dominated solutions are added to it. As soon as the number of elements in the archive meets its limit, one of the elements should be omitted from it. The strategy to remove an element from the archive is to find the most crowded hypercube and randomly select one of the elements in it to extract. For this purpose, objective functions and solution space are divided into hyper-rectangles. When archive overflow is reached, one particle from the most crowded hyper-rectangle is randomly selected and removed. In MOGSA, the distance between each particle to its nearest neighbor is calculated, and its mass is updated. Distribution of archived particles is done by using the niching technique. In this algorithm, the archived particles apply gravitational force. The list below shows the pseudocode of MOGSA that is similar to Hassanzadeh and Rouhani [43]. (1)
Set the population size,G0, a and ε.(2)
Initialize population of solutions.(3)
Evaluate population (with regards to objective functions: CO2 and Cost).(4)
Update gravitational constant.(5)
Calculate the masses.(6)
Calculate force and acceleration between objects.(7)
Update velocity and position per each object.(8)
Select set of best solutions (Pareto front)  //Update external archive//.(9)
Generate new population.(10)
If the fitness function is satisfied ort=Tmax go to 4 otherwise continued.(11)
Figure set of Pareto solutions and end.


## 4. Experimental Results and Discussion

In this section, a multi-objective mixed integer linear programming supply chain logistics network model is applied to a real world company. The case study is conducted to illustrate the efficiency and consistency of the proposed method to solve multi-objective optimization in supply chain networks. This case is motivated by the Pars Khodro Company, which is located in Tehran, Iran and is one the largest car manufacturers in the Middle East. The company has been under pressure to satisfy recently enacted environmental legislation. Our purpose is to provide some strategic viewpoints and recommendations for the supply chain network design where the CO2 emission issue is being considered. In this regard, we set up a bi-objective optimization model and conducted a case analysis based on a network consisting of the following: the three groups of tier-1 suppliers (metal-based parts, electrical parts and plastic parts) which include 17, 22, and 17 tier-1 suppliers in each group, respectively; the 9 potential locations for establishing the manufacturer; the 23 potential locations for establishing the distribution centres; the 47 retailers; the three trucks type, including a new 12-ton truck, an 8-ton truck and an old 12-ton truck. Figure 4 illustrates some of the potential networks in the supply chain logistics network. The formulation of this case contains 31 binary variables, 1792 continuous variables, and 302 constraints.Figure 4
Potential logistics network in the case study supply chain.In Table1, the parameters of the problem are given. The transportation costs and CO2 emissions for all truck types are defined as corresponding to every kilometer between facilities on each layer of the supply chain network. Moreover, the fixed costs and capacity of the same facilities and the demand of each retailer are discriminated.Table 1
Value of parameters in the numerical experiments.
Parameter
Value
Parameter
Value
Parameter
Value
ε
1, ε3
12
ρ
e
0.3
φ
1
600
ε
2
6
ρ
m
0.5
φ
2
400
τ
1
0.2
η
1
-
η
47
12000
φ
3
800
τ
2, τ3
0.15
γ
1
w
-
γ
23
w
80
α
1
f
-
α
9
f
50,000,000
γ
1
f
1000
γ
1
r
-
γ
47
r
65
α
1
w
-
α
23
w
5,000,000
γ
2
f
1200
γ
3
f
-
γ
9
f
900
ω
p
2For MOGSA implementation as a solving method, material and product flows between facilities in automotive supply chains such as suppliers, manufacturers, distributors, and retailers as objects with some respective masses should be considered. And each particle is as a vector with 7 dimensions,M=[Xikt1,Xjkt2,Xnkt3,Xklt4,Xlmt5,Yk1,Yl2]. In this case, we set MOGSA regarding to the parameter setting that was mentioned in Table 2. The problem is solved by the MOGSA method, and it is implemented using MATLAB software. All the experiments are conducted using a laptop with an Intel Core 3 Duo 2.19 GHz and 1 GB RAM.Table 2
Parameters values of MOGSA algorithm.
G
constant
8
Mutation prop
0.5
Number  of  divisions  in each dimension
[3 2 2]Since the Pareto frontier can provide the decision maker with a portfolio of alternative “optimal solutions,” the solutions to the assignment problem estimated previously were the marginal points, as shown in Figure5. The minimization of total logistics cost, including freight transportation system costs, operation costs, and the initial cost of opening and running all the facilities is shown on the upper left side, and the minimisation of CO2 emissions is shown on the lower right side. It clearly demonstrates the trade-off between the total cost and the total CO2 emissions. It coincides with our intuition that a lower CO2 emission can only be reached by applying more investment. Figure 5 shows 20 solutions, but these are only a subset of the Pareto optimal solutions. However, the algorithm can estimate less than 20 solutions because there might be unfeasible solutions in the iterations. The relationship between cost and CO2 emissions in the entire network is not exactly linear. The linearity of the Pareto optimal solution is not necessary, even if all the objective functions are linear. Specifically, the changed cost amounts are not necessarily proportional as the CO2 emissions allowed are decreased. Also, we can find that the absolute value of the slope of the curve decreases, and the curves became flatter as the total cost increases.Figure 5
Pareto optimal solution to the case study.Solution selection is depended on how much total cost the company would like to invest for environmental issue. For example, if the company prefers to produce minimum CO2 emission (180 tones) should pay about 990,560 thousand US dollar for total logistic cost. If the company is not under government and customer pressure to decrease CO2 emission, it can select the cheaper solution such as the first solution. In this case total logistic cost would be 290,260 thousand US dollar but more than 680 tones CO2 emission is emitted in environment. Decision maker compare the Pareto solution set with the real situation for the company regarding the total cost and environmental policy and finally choose one of the solutions to trade-off these limitations.


## 5. Conclusion

The quantitative relationship between CO2 emissions and logistics costs has been gaining   importance in the logistics field because of global warming, as well as rapidly increasing fuel costs. This study is an attempt to estimate the relationship by using multi-objective linear programming mathematical modelling in the automotive supply chain network. This model will have an important application in the regional or global automotive supply chain network design with green considerations. The model is a multi-objective model that consists of minimising the total cost and the CO2 emission—the type of pollution human kind is most responsible for. The MOGSA method is used as a new swarm intelligence algorithm to solve the model and obtain a Pareto optimal set. After that, the model was tested by a real case study to illustrate the model and MOGSA. The Pareto optimal solutions given by the model provide a portfolio of configuration for decision makers, and the computation experiments show that the model can serve as an effective tool in designing a green supply chain network. The results show that the trade-off curves have almost a linear relationship, in that logistics costs should be higher, as a greater reduction of CO2 emissions is required. The quantity of the relationship varies, ranging from 290 to 990 billion US dollars in terms of the CO2 emission levels. Also the absolute value of the slope of the curve decreases, and the curves became flatter as the total cost increases.The research can be extended to minimise the lead time or ensuring just-in-time arrival that is not taken into account in this study. Furthermore, since considering reverse logistics in the supply chain network is attracting great interest among researchers and it can be added to the network modelling, authors are starting to extend the model as a closed loop. Accordingly, nonlinear programming would be a better option as a modelling method to find a more accurate solution, which should be investigated more in the future. Also comparing the MOGSA result with the other multi-objective optimization algorithm such as Nondominated Sorting Genetic Algorithm (NSGA) or Multi-objective Particle Swarm Optimization (MOPSO) can be one of the future research.



---

*Source: 823876-2013-09-18.xml*