# Designing an Agile Closed-Loop Supply Chain with Environmental Aspects Using a Novel Multiobjective Metaheuristic Algorithm

Authors: Liu Kai; Ramina Malekalipour Kordestanizadeh
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Success in supply chain implementation depends on the way of dealing with market changes and customer needs. Agility is a concept that has been introduced in recent years to improve the supply chain. On the other hand, paying attention to environmental problems is another issue, and chains are trying to increase their popularity by focusing on this issue. Considering the importance of this issue, designing a multiobjective closed-loop supply chain network has been discussed in this research. The main contribution of this research is the integration of green and agility concepts in supply chain design. In this regard, a mathematical model is presented with economic, environmental, and agility objectives. First, the mathematical model is solved using the Epsilon constraint method, and then, the multiobjective weed algorithm is proposed to solve the model. The results of comparisons between the two methods show that the multiobjective weed algorithm has performed well in terms of various metrics of NPS, SNS, and Max Spread. In terms of the solving time, the average solving time of this algorithm was about 0.1% of the solving time of the Epsilon constraint method. Moreover, all cases show the superiority of the multiobjective weed algorithm over the Epsilon constraint method in solving the proposed mathematical model.


---

## Body

## 1. Introduction

The inability to deal with challenges can lead to falling and loss in competition among supply chains. This may seem concerning to some people, but it more or less emphasizes the need to develop new ways of thinking, finding solutions, and new criteria to maintain its potential ability in solving changing and unpredictable market challenges. Many analyses have proposed different solutions to deal with these unstable changes more successfully: networking, reengineering, modular organizations, virtual cooperation, high-efficient organizations, staff development, and flexible and just-in-time manufacturing systems [1]. Among them, the concepts of the agile institution and lean manufacturing system are more influential.An agile supply chain means adaptability and flexibility and plans the chain in such a way that it uses a new way to deliver services and goods to respond to customer needs and reduce premium costs quickly. The customer who uses the goods or service at least wants to use innovative technologies and products which are different from previous products and services [2]. An agile chain changes its work environment using information technology so that in addition to responding to changes in customer needs, it can reduce prices effectively, forecast these changes to some extent, and deal with them with appropriate decisions. The organization’s human resources and innovational capabilities are two essential and constituent factors in creating an agile company [3]. An agile supply chain emphasizes the flexibility and improvement of the quality of manufacturing and responds quickly to market changes and customer needs. Each supply chain, which is considered an agile supply chain, makes its own operational changes. In the changes section, new ideas and solutions resulting from innovation in the organization are provided by the human resources. This creativity in human resources and innovation cause improvement in the technology of products and services that can compete in a competitive market. In this regard, the use of agility criteria by a supply chain not only maintains it but also distinguishes it from other supply chains. Due to constant changes in the market, an agile organization must also constantly change over time in order to adapt to market changes [4].Paying attention to the environmental problems of supply chains is one of the important concerns of communities. Pollution caused by manufacturing and transportation causes the most damage to the environment. Therefore, in addition to trying to achieve more profits and more agility to deal with market changes, supply chains must consider minimizing their environmental pollution [5]. Accordingly, this study deals with the integration of agile and green supply chains. In the second section, the newest relevant studies will be presented. In the third section, mathematical models and, in the fourth section, optimization methods will be introduced. The fifth section will present the numerical results of model optimization, and the research will be summarized in the seventh section.The trend toward a green and agile supply chain creates a vital opportunity to simultaneously address issues of sustainability and business and environmental performance. The chain has led companies to redesign their products to fit the environment, in addition to the obvious costs, minimizing intangible costs that are not included in the cost of the product and can have a negative impact on the environment. Reviewing the studies conducted in the field of supply chain network design, it can be seen that despite the importance of both agile and green aspects, these two aspects have not been considered together and only separate studies have been done in this regard. Therefore, in order to cover the research gap, this study will design a green and agile closed-loop supply chain network for the first time. Also, the metaheuristic multiobjective invasive weed algorithm will be used to optimize this problem, which is one of the newest multiobjective solution methods.


## 2. Literature Review

The green supply chain has received much attention in recent years. Liang and Questa [6] designed the green supply chain for fuel transportation. They introduced energy optimization and emissions reduction as the most important objectives of this chain. This study that was conducted in Japan offers a solution that balances energy costs and environmental pollution. Yadav et al. [7] designed a two-objective optimization supply chain network with a sustainability approach. In the mathematical model, the reduction of pollutants is considered along with the reduction of costs. GAMS software has been used to solve this mathematical model. Hasani et al. [8] presented a multiobjective optimization approach for green and resilient supply chain network design. In the proposed model of this research, cost reduction along with environmental goals are included as objective functions in the mathematical model, and a robust multiobjective optimization model is designed. Strength Pareto Evolutionary Algorithm 2 (SPEA2) has been used to solve the proposed model. The results indicate that the designed supply chain network will be able to meet the needs of customers in an agile and green form.Recycling is an industrial process in which disposable products are recycled with the aim of reusing and returning products to the consumption cycle. Paying attention to this concept along with the green supply chain will strengthen the supply chain as much as possible. Prakash et al. [9] designed a closed-loop supply chain using an integrated robust optimization approach. The proposed model is considered based on risk and demand uncertainty. Mixed-integer programming has been used in the design of this network. The goods are distributed directly from the factory to the customer through distribution channels. In addition, supply risk and transportation risk are included in the model, along with uncertain customer demand. The proposed model is implemented in an Indian e-commerce company. Boronoos et al. [10] designed a multiobjective closed-loop green supply chain network. Reducing environmental pollutants along and reducing costs are considered objective functions in the mathematical model. Flexibility and uncertainty are included as innovations in the mathematical model. In order to show the efficiency of the proposed mathematical model, it has been implemented in the printing industry, the results of which show the efficiency of this model.Yavari and Geraeli [11] used a metaheuristic method to robust optimize the green closed-loop supply chain network using an innovative method. The intended product is perishable products. The proposed model was multiperiod and included four echelons: suppliers, manufacturers, warehouses, retailers, and product collection centers. Paksoy et al. [12], in their proposed model, considered a green closed-loop supply chain along with the lean concept. This model simultaneously minimizes transportation costs, purchase, and operating costs and also minimizes fixed facility costs, environmental costs, shipping costs, and late delivery. The fuzzy analytic hierarchy process method is used to solve a numerical example to confirm the validity of the proposed model. Gholizadeh and Fazlollahtabar [13] investigated the robust optimization of a green closed-loop supply chain scheme and solved it using a genetically modified algorithm. The proposed model is investigated under an uncertain environment. The intended objective is to reduce costs and environmental pollutants. A case study was conducted in the ironworks industry to bring the problem closer to the real world. The results of the implementation of the model with a genetic algorithm indicated the effectiveness and application of the proposed model. Ghahremani-Nahr et al. [14] designed a green closed-loop supply chain network considering the discount. The proposed model is presented in terms of uncertainty.Nowadays, many organizations and companies face a lot of competition and an uncertain environment that has intensified due to technological innovations and changing customer needs [4]. As mentioned, one of the ways to deal with such challenges is agility. Moradi et al. [1] introduced a multiobjective model for agile supply chain design. In this mathematical model, the pricing of products was carried out in uncertain conditions. The objectives of this mathematical model were to reduce total costs and total delivery time. A fuzzy programming approach was used to solve this mathematical model. Numerical examples are used to show the validity of the proposed model. Mahmoodi [15] designed a multiobjective lean supply chain network with transportation constraints. In this mathematical model, supply chain agility and manufacturing and distribution risks were considered. The objectives of this mathematical model included minimizing total costs, minimizing total supply chain risk, and maximizing supply chain flexibility. Using the NSGA-II algorithm, the proposed model is solved. The results showed the effect of agility on supply chain design and better market policies. Pahlevan et al. [16] proposed a mathematical model to design the closed-loop supply chain of the aluminum industry. They used Red Deer Algorithm (RDA) to optimize this model. Du et al. [17] presented a systematic literature review on supply chain agility. In this study, the evolution of the concept of the agile supply chain has been investigated, and a systematic literature review approach has been used to understand the research gap. The results show that, in the current competitive environment, supply chain agility increases customer satisfaction and thus profit. On the other hand, achieving the concept of agility in the organization requires continuous planning, so the allocation of resources to improve aspects of agility is one of the requirements of organizations in global markets. Therefore, in this research, methods to achieve supply chain agility are mentioned.Table1 provides a summary of studies in this field.Table 1
Summary of the works related to the topic.
ResearcherYearLocationAllocationDistributionMultiproductMultiperiodEnvironmental pollutionAgilityLiang and Quesada2019✓✓✓Moradi et al.2019✓✓✓✓Mahmoodi2019✓✓✓✓Paksoy et al.2019✓✓✓Yadav et al.2019✓✓✓✓Yavari, and Geraeli2019✓✓✓Ghahremani-Nahr et al.2020✓✓✓✓Gholizadeh and Fazlollahtaba2020✓✓Hasani et al.2021✓✓✓✓✓Boronoos et al.2021✓✓✓Pahlevan et al.2021✓✓✓✓Current research2021✓✓✓✓✓✓✓As can be seen in Table1, this research gap includes the integration of green supply chain and agile supply chain. None of the reviewed studies have investigated the simultaneous optimization of agility and environmental pollution. Also, the metaheuristic multiobjective invasive weed algorithm will be used to optimize this problem, which is one of the newest multiobjective solving methods.


## 3. Mathematical Model of the Research

This study examines the supply chain of dairy products. In this supply chain, milk and other raw materials are received from various suppliers. Excess demand is also supplied by the milk collection station, which is an intermediary. After the products were produced, they will be sent to the distribution centers, and they maintain these products after inspecting them and distribute them among the customers at the proper time. Distributors have the ability to deliver stored products in the future. These conditions cause the corruption of some of the products in the warehouse, which need to be returned to the factory. The manufacturing factory processes a percentage of the returned products and takes them back to the distribution and sales cycle. The rest of these products, which were returned, are disposed of as waste.Other assumptions of the mathematical model are as follows:(i)
The supply chain proposed in this research consists of echelons. These four echelons are supplier, manufacturer, distributor, and customers(ii)
In order to supply the excess required milk, in addition to the suppliers, a milk collection station has been considered(iii)
There are different transportation systems among the chain members; each has the fixed costs and variable costs(iv)
Only access to one transportation system between each supplier and distributor as well as the between each distributor and customer is allowed(v)
The chain must meet the total forecasted demand of customers(vi)
There are a number of potential centers for distribution centers and several potential centers for factories, in which the network is required to establish one or more(vii)
The number of centers required for the establishment is not known, but the mathematical model must choose the one that is the best possible, according to the demand of the network and the costs of establishment and transportation(viii)
The objectives of the proposed mathematical model of this research are reducing total costs and reducing greenhouse gases and optimizing supply chain flexibility


### 3.1. Indices

S: index of fixed location of suppliers (s = 1, 2, …, S)i: index of fixed location of factories (i = 1, 2, …, I)j: index of potential locations for distribution centers (j = 1, 2, …, J)c: index of fixed locations of customers (c = 1, 2, …, C)p: index of products (p = 1, 2, …, P)r: index of raw materials (r = 1, 2, …, R)l: index of product transportation systems (l = 1, 2, …, L)t: index of time periods (t = 1, 2, …, T)


### 3.2. Parameters

Dcpt: demand of customer c for product p during period tSCsrt: the purchasing cost of one unit of raw material r from supplier s during period tMC˙ipt: the manufacturing cost of each unit of product p in factory i during period tICjpt: the inspection and recycling cost of each unit of product p in distribution center j during period tHCjPt: maintaining the cost of each unit of product p in distribution center j during period tFXit: the fixed cost of setting up the factory i during period tFYjt: the fixed cost of setting up the distribution center j during period tCSsrt: the capacity of supplier s for supplier r during period tCXit: the manufacturing capacity in factory i during period tCYjt: the capacity of the distribution center j during period tCSIsirt: the transportation cost of each unit of raw material r from supplier s to factory i during period tCIJijplt: the transportation cost of each unit of product p from factory i to distribution center j during period t with transportation system lCJCjcplt: the transportation cost of each unit of product p from distribution center j to customer c with transportation system l during period tCJIjiplt: the transportation cost of each unit of product p from distribution center j to factory i during period t with transportation system lCTRlt: the fixed cost of using transportation system l during period tESIsirt: the unit of CO2 pollution resulting from the transportation of raw material r from supplier s to factory i during period tEIJijplt: the unit of CO2 pollution resulting from the transportation of product p from factory i to distribution center j during period t with transportation system lEJIjiplt: the unit of CO2 pollution resulting from the transportation of product p from distribution center j to factory i during period t with transportation system lEJCjcplt: the unit of CO2 pollution resulting from the transportation of product p from distribution center j to customer c with transportation system l during period tEC˙ipt: the unit of CO2 pollution resulting from the manufacture of product p in factory i during period tnrp: the consumption coefficient of raw material r in product pmp: the rate of capacity utilization in manufacture of product pRp: the percentage of return of product p to the factoryRpp: the percentage of processing product p from the returned products to the factoryBM: a very large numberW1: the impact (weight) of agility of factories on the agility of the whole supply chainW2: the impact (weight) of agility of distribution centers on the agility of the whole supply chainW3: the impact (weight) of agility of suppliers on the agility of the whole supply chain


### 3.3. Decision Variables

QSIsirt: the amount of raw material r sent from supplier s to factory i during period tQS0irt: the amount of raw material r sent from the raw material supply station to factory i during period tQIipt: the amount of product p manufactured in factory i during period tQIJijplt: the amount of product p sent from factory i to distribution center j with transportation system l during period tINVjpt: the inventory of product p in distribution center j at the end of period tQJCjcplt: the amount of product p transported from the distribution center j to the customer c with the transport system l during period tQJIjiplt: the amount of product p sent from the distribution center j to factory i with the transportation system l during period tWst: a binary variable equal to 1 if it is selected from the supplier s during period tXit: a binary variable equal to 1 if a factory is established at point i during period tYJt: a binary variable equal to 1 if the distribution center is set up at point j during period tAijlt: a binary variable equal to 1 if the transportation system l connects factory i and distribution center j during period tBjclt: a binary variable equal to 1 if the transportation system l connects the distribution center j to the customer c during period t


### 3.4. Mathematical Model Relations

(1)MinZ1=∑t∑jFYjtYjt+∑t∑iFXitXit+∑t∑s∑i∑rQSIsirt+QS0irtSCsrt+∑i∑p∑tQIiptMCipt+∑t∑c∑j∑p∑lQJCjcpltICJpt+∑t∑j∑pINVjptHCpt+∑t∑s∑i∑rQSIsirtCSIsirt+∑t∑i∑j∑p∑lQIJijpltCIJijplt+∑t∑j∑c∑p∑lQJCjcpltCJCjcplt+∑t∑i∑j∑p∑lQJIjipltCJIjiplt+∑t∑j∑i∑lAijltCTRlt+∑t∑j∑c∑lBjcltCTRlt,(2)MinZ2=∑i∑p∑tQIiptECipt+∑t∑s∑i∑rQSIsirtESIsirt+∑t∑i∑j∑p∑lQIJijpltEIJijplt+∑t∑j∑c∑p∑lQJCjcpltEJCjcplt+∑t∑i∑j∑p∑lQJIjipltEJIjiplt,(3)MinZ3=∑tw1AIt+∑tw2AJt+∑tw3ASt,(4)AIt=∑iCXitXit−∑c∑pDcpt,(5)ASt=∑s∑rCSsrtWst−∑c∑pDcpt,(6)AJt=∑jCYjtYjt−∑c∑pDcpt,(7)∑j∑P∑lnrpQIJijplt=∑sQSIsirt,∀i,r,t,(8)INVjpt−1+∑i∑lQIJijplt=INVjpt+∑c∑lQJCjcplt+∑i∑lQJIjiplt,∀j,p,t,(9)∑j∑lQJCjcplt=Dcpt,∀c,p,t,(10)∑iQSIsirt≤CSsrtWst,∀s,r,t,(11)∑j∑p∑lmpQIJijplt≤CXitXit,∀i,t,(12)∑pINVjpt+∑c∑p∑lQJCjcplt≤CYjtYjt,∀j,t,(13)∑lAijlt≤1,∀i,j,t,(14)∑lBjclt≤1,∀j,c,t,(15)Aijlt≤∑pQIJijplt,∀i,j,l,t,(16)Bjclt≤∑pQJCjcplt,∀j,c,l,t,(17)∑pQIJijplt+QJIjiplt≤BM∗Aijlt,∀i,j,l,t,(18)∑pQJCjcpt≤BM∗Bjclt,∀j,c,l,t,(19)∑iQJIjiplt=Rp∗INVjpt,∀j,p,l,t,(20)∑jQIJijplt=QIipt+Rpp∗∑jQJIjiplt,∀i,p,l,t.Equation (1) shows the total costs in each period of the chain. These costs are fixed costs related to the establishment of distribution centers and purchases from suppliers, manufacturing costs, operating costs formed in distribution centers, inventory costs in distribution centers, and transmission costs for different transport systems in the chain. Equation (2) also minimizes the total amount of CO2 emitted from the manufacture of products, transportation from supplier to the factory, transportation from factory to distribution centers, and transportation from distribution centers to customers. Equation (3) determines the total amount of supply chain agility. Therefore, because the agility of supply, manufacture, and distribution may not be equally significant, the weighted sum of agility of supply, manufacture, and distribution is introduced as the agility of the whole supply chain. Equation (4) shows the agility of factories or the level of manufacture. Equation (5) shows the agility of suppliers or the level of supply. Equation (6) shows the agility of the distributors or the level of distribution. Equation (7) shows that, in each period, the amount of raw material entering each factory is equal to the amount of output from that factory in the same period. Equation (8) guarantees that, in each period and for each product, the amount of input to each of the distribution centers and the remaining inventory from the previous period is equal to the amount sent to customers and the rest of the inventory at the end of the period. This equation is known as the inventory balance equation. Equation (9) states that, in each period and for each product, the inventory available in each of the distribution centers must be able to meet the demand for the product in question. Equation (10) ensures that the amount of each of the raw materials sent from the suppliers does not exceed their capacity. Equation (11) expresses the capacity of materials in factories similar to suppliers. Equation (12) states that inventory in each distribution center should not exceed the capacity of the center. Equations (13) and (14) state that only one transportation system can be used in each member of the chain. Equations (15) and (16) express that members of the chain who work together are allowed to use the transportation system. Equations (17) and (18) express that unrelated members of the chain do not send any products. Equation (19) determines the amount of returned products. This amount is determined as a percentage of the distributor’s inventory. Equation (20) indicates that the amount of products sent from the factory to the distribution centers is equal to the amount manufactured in the factory, which is a percentage of returned products that is reprocessed.


## 3.1. Indices

S: index of fixed location of suppliers (s = 1, 2, …, S)i: index of fixed location of factories (i = 1, 2, …, I)j: index of potential locations for distribution centers (j = 1, 2, …, J)c: index of fixed locations of customers (c = 1, 2, …, C)p: index of products (p = 1, 2, …, P)r: index of raw materials (r = 1, 2, …, R)l: index of product transportation systems (l = 1, 2, …, L)t: index of time periods (t = 1, 2, …, T)


## 3.2. Parameters

Dcpt: demand of customer c for product p during period tSCsrt: the purchasing cost of one unit of raw material r from supplier s during period tMC˙ipt: the manufacturing cost of each unit of product p in factory i during period tICjpt: the inspection and recycling cost of each unit of product p in distribution center j during period tHCjPt: maintaining the cost of each unit of product p in distribution center j during period tFXit: the fixed cost of setting up the factory i during period tFYjt: the fixed cost of setting up the distribution center j during period tCSsrt: the capacity of supplier s for supplier r during period tCXit: the manufacturing capacity in factory i during period tCYjt: the capacity of the distribution center j during period tCSIsirt: the transportation cost of each unit of raw material r from supplier s to factory i during period tCIJijplt: the transportation cost of each unit of product p from factory i to distribution center j during period t with transportation system lCJCjcplt: the transportation cost of each unit of product p from distribution center j to customer c with transportation system l during period tCJIjiplt: the transportation cost of each unit of product p from distribution center j to factory i during period t with transportation system lCTRlt: the fixed cost of using transportation system l during period tESIsirt: the unit of CO2 pollution resulting from the transportation of raw material r from supplier s to factory i during period tEIJijplt: the unit of CO2 pollution resulting from the transportation of product p from factory i to distribution center j during period t with transportation system lEJIjiplt: the unit of CO2 pollution resulting from the transportation of product p from distribution center j to factory i during period t with transportation system lEJCjcplt: the unit of CO2 pollution resulting from the transportation of product p from distribution center j to customer c with transportation system l during period tEC˙ipt: the unit of CO2 pollution resulting from the manufacture of product p in factory i during period tnrp: the consumption coefficient of raw material r in product pmp: the rate of capacity utilization in manufacture of product pRp: the percentage of return of product p to the factoryRpp: the percentage of processing product p from the returned products to the factoryBM: a very large numberW1: the impact (weight) of agility of factories on the agility of the whole supply chainW2: the impact (weight) of agility of distribution centers on the agility of the whole supply chainW3: the impact (weight) of agility of suppliers on the agility of the whole supply chain


## 3.3. Decision Variables

QSIsirt: the amount of raw material r sent from supplier s to factory i during period tQS0irt: the amount of raw material r sent from the raw material supply station to factory i during period tQIipt: the amount of product p manufactured in factory i during period tQIJijplt: the amount of product p sent from factory i to distribution center j with transportation system l during period tINVjpt: the inventory of product p in distribution center j at the end of period tQJCjcplt: the amount of product p transported from the distribution center j to the customer c with the transport system l during period tQJIjiplt: the amount of product p sent from the distribution center j to factory i with the transportation system l during period tWst: a binary variable equal to 1 if it is selected from the supplier s during period tXit: a binary variable equal to 1 if a factory is established at point i during period tYJt: a binary variable equal to 1 if the distribution center is set up at point j during period tAijlt: a binary variable equal to 1 if the transportation system l connects factory i and distribution center j during period tBjclt: a binary variable equal to 1 if the transportation system l connects the distribution center j to the customer c during period t


## 3.4. Mathematical Model Relations

(1)MinZ1=∑t∑jFYjtYjt+∑t∑iFXitXit+∑t∑s∑i∑rQSIsirt+QS0irtSCsrt+∑i∑p∑tQIiptMCipt+∑t∑c∑j∑p∑lQJCjcpltICJpt+∑t∑j∑pINVjptHCpt+∑t∑s∑i∑rQSIsirtCSIsirt+∑t∑i∑j∑p∑lQIJijpltCIJijplt+∑t∑j∑c∑p∑lQJCjcpltCJCjcplt+∑t∑i∑j∑p∑lQJIjipltCJIjiplt+∑t∑j∑i∑lAijltCTRlt+∑t∑j∑c∑lBjcltCTRlt,(2)MinZ2=∑i∑p∑tQIiptECipt+∑t∑s∑i∑rQSIsirtESIsirt+∑t∑i∑j∑p∑lQIJijpltEIJijplt+∑t∑j∑c∑p∑lQJCjcpltEJCjcplt+∑t∑i∑j∑p∑lQJIjipltEJIjiplt,(3)MinZ3=∑tw1AIt+∑tw2AJt+∑tw3ASt,(4)AIt=∑iCXitXit−∑c∑pDcpt,(5)ASt=∑s∑rCSsrtWst−∑c∑pDcpt,(6)AJt=∑jCYjtYjt−∑c∑pDcpt,(7)∑j∑P∑lnrpQIJijplt=∑sQSIsirt,∀i,r,t,(8)INVjpt−1+∑i∑lQIJijplt=INVjpt+∑c∑lQJCjcplt+∑i∑lQJIjiplt,∀j,p,t,(9)∑j∑lQJCjcplt=Dcpt,∀c,p,t,(10)∑iQSIsirt≤CSsrtWst,∀s,r,t,(11)∑j∑p∑lmpQIJijplt≤CXitXit,∀i,t,(12)∑pINVjpt+∑c∑p∑lQJCjcplt≤CYjtYjt,∀j,t,(13)∑lAijlt≤1,∀i,j,t,(14)∑lBjclt≤1,∀j,c,t,(15)Aijlt≤∑pQIJijplt,∀i,j,l,t,(16)Bjclt≤∑pQJCjcplt,∀j,c,l,t,(17)∑pQIJijplt+QJIjiplt≤BM∗Aijlt,∀i,j,l,t,(18)∑pQJCjcpt≤BM∗Bjclt,∀j,c,l,t,(19)∑iQJIjiplt=Rp∗INVjpt,∀j,p,l,t,(20)∑jQIJijplt=QIipt+Rpp∗∑jQJIjiplt,∀i,p,l,t.Equation (1) shows the total costs in each period of the chain. These costs are fixed costs related to the establishment of distribution centers and purchases from suppliers, manufacturing costs, operating costs formed in distribution centers, inventory costs in distribution centers, and transmission costs for different transport systems in the chain. Equation (2) also minimizes the total amount of CO2 emitted from the manufacture of products, transportation from supplier to the factory, transportation from factory to distribution centers, and transportation from distribution centers to customers. Equation (3) determines the total amount of supply chain agility. Therefore, because the agility of supply, manufacture, and distribution may not be equally significant, the weighted sum of agility of supply, manufacture, and distribution is introduced as the agility of the whole supply chain. Equation (4) shows the agility of factories or the level of manufacture. Equation (5) shows the agility of suppliers or the level of supply. Equation (6) shows the agility of the distributors or the level of distribution. Equation (7) shows that, in each period, the amount of raw material entering each factory is equal to the amount of output from that factory in the same period. Equation (8) guarantees that, in each period and for each product, the amount of input to each of the distribution centers and the remaining inventory from the previous period is equal to the amount sent to customers and the rest of the inventory at the end of the period. This equation is known as the inventory balance equation. Equation (9) states that, in each period and for each product, the inventory available in each of the distribution centers must be able to meet the demand for the product in question. Equation (10) ensures that the amount of each of the raw materials sent from the suppliers does not exceed their capacity. Equation (11) expresses the capacity of materials in factories similar to suppliers. Equation (12) states that inventory in each distribution center should not exceed the capacity of the center. Equations (13) and (14) state that only one transportation system can be used in each member of the chain. Equations (15) and (16) express that members of the chain who work together are allowed to use the transportation system. Equations (17) and (18) express that unrelated members of the chain do not send any products. Equation (19) determines the amount of returned products. This amount is determined as a percentage of the distributor’s inventory. Equation (20) indicates that the amount of products sent from the factory to the distribution centers is equal to the amount manufactured in the factory, which is a percentage of returned products that is reprocessed.


## 4. Solution Methods

The proposed mathematical model in this research is multiobjective. Therefore multiobjective optimization methods are needed to solve it. Therefore, the Epsilon constraint method and multiobjective invasive weed metaheuristic algorithm have been used, which are described in the following.


### 4.1. Epsilon Constraint Method

The Epsilon constraint method is an exact solution method for solving mathematical models with multiobjective functions. This method leads to the optimization of multiobjective mathematical models by maintaining only one objective function and transferring other objective functions to constraints. This general model is as follows.(21)Minimizefμxsubject tofmm≤εm,m=1,2,…,Mandm≠μgjx≥0,j=1,2,…,J.In the above model, the parameterεm indicates the upper limit of the value fm. To show how this method works, consider the problem with the two objectives f1 and f2. Also, suppose that f2remains in the objective function and f1 is added to the set of constraintsf1x≤ε1.Assume that initiallyε1=ε1c. The problem resulting from considering this constraint divides the main feasible objective space into two parts of f1x≤ε1c and f1x>ε1c. The left part of the feasible space obtains the feasible solution to the problem. Now, the mission of the resulting problem is to find a solution that is the minimum value of this feasible space. It is obvious from the figure that c is the minimum solution. According to this figure, intermediate solutions can be obtained in problems with nonconvex objective space using the ε-constraint method.


### 4.2. Multiobjective Invasive Weed Optimization Algorithm

The invasive weed optimization algorithm (IWO) is a nature-inspired metaheuristic algorithm that is inspired by weed behavior and its growth. They grow widely and strikingly so that they cannot be removed and controlled by humans. One claim about weeds is that weeds always win. In general, the reasons for this claim can be stated as follows:(i)
Existence of weeds after thousands of years of agriculture(ii)
Existence of weeds even after using various pesticides(iii)
The emergence of new species of weeds widely on the ground(iv)
Adapting to the environmentThe above characteristics show that weeds are vigorous and interfere with agriculture. On the other hand, weeds adapt themselves to the environment and change their behavior to grow. The success of weeds depends on their ecology and biology. In the following, optimization based on invasive weed behavior is investigated. The stages of the invasive weed optimization algorithm are as follows:(i)
First stage: spreading the seeds in the desired area(ii)
Second stage: seed growth according to utility (sprouting) and spatial dispersal(iii)
Third stage: survival of weeds with more utility (competitive elimination)(iv)
Fourth stage: continuing the process until reaching the plants with the most utilityDetails of the steps of the invasive weeds optimization algorithm are as follows:First step: generating the random initial population and evaluating their objective functionAccording to this stage, a certain number of seeds are produced from each plant. The seeds produced depend on their fitness level as well as the maximum and minimum fitness of the colony of the weeds, and this number will increase linearly with increasing the amount of fitness. On the other hand, the number of seeds produced by each plant depends on the proportion of the plant to its environment; this relationship is linear, which means that the best solution from the current population gives the highest new solution and the worst solution leads to the generation of the lowest new solution. The number of solutions varies linearly between these two limits.Second stage: reproduction based on fitness and updating of standard deviationInvasive weeds may reproduce with or without using sex cells, depending on the type of the plant. In the method of sexual reproduction, seeds or spores are used. In this case, the plant is born and begins its life by being fertilized. Then, the seeds are dispersed by various factors such as wind, water, and animals, until the seed has a chance to grow. If the conditions are good, the seeds will start to germinate and grow. The process of growth will continue with other neighbor plants until they become mature plants. In the last stage of the life cycle, they become flowering plants and cause seed production.Third stage: competitive exclusionIf a plant does not reproduce, it will die. Therefore, competition between plants is needed to limit the maximum number of plants in a colony. After producing seeds around each weed, we can only transfer the predetermined maximum number of plants (Pmax) of the total weeds and seeds to the next generation. Plants that have a chance of surviving reproduce and repeat the above stages so that the solutions obtained in each of the iterations would be fitter. This mechanism gives plants with low proportions a chance to reproduce, and if the seeds produced by them have a good proportion in the colony, then they can survive. This algorithm stops when the number of iterations reaches the maximum number allowed. The maximum number of preset plants can be equal to the number of the initial population, in which case one of the parameters of the algorithm is removed.Fourth stage: checking the termination conditionsTypes of termination conditions in metaheuristic methods are as follows:(i)
Obtaining an acceptable minimum of solutionsIn this case, which is one of the termination conditions, assume that the total cost of a company is 1000 monetary units and the manager of the company tends to reduce this cost to 800 units, which is a desirable amount from the manager’s point of view. It is noteworthy that the value of 800 is not the optimal point, and the value of the objective function may become less than 800 units. For example, in Figure1, point A shows the desired point from the manager’s point of view, which is one of the termination conditions.(ii)
The passage of time or determined iterationIn this case, a determined iteration is introduced for the algorithm to end the algorithm. For example, in iteration 100, the algorithm ends, and there would be no need for iterations after 100, although the value of the objective function may improve after the determined iteration.(iii)
The passage of time or determined iteration without observing much improvement in the resultIf, after implementing the algorithm in successive iterations, the value of the objective function does not improve or improves insignificantly, the algorithm will stop. For example, if it does not improve from pointa to b with n iterations or improves insignificantly, the algorithm will stop. In Figure 2, the value of the objective function is constant in iteration a to iteration b that causes the algorithm to stop.Figure 1
Pseudocode of the fast nondominated sorting procedure.The general structure of MOIWO and IWO is similar. The only difference between them is that, in each iteration, the nondominated solutions are separated. This method is presented in Section4.2.1. The pseudocode of the MOIWO algorithm is shown in Figure 3.Figure 2
Comparison of solving methods based on Max Spread metric.Figure 3
Pseudocode of MOIW algorithm.


#### 4.2.1. Nondominated Sorting

First, two parametersnp and Sp must be calculated; np is actually the number of solutions that have dominated the solution p, and Sp is the set of solutions that have been dominated by the solution p. Then, a set of solutions can be ranked, and they can be placed on different fronts in terms of nondomination degrees. All solutions with np=0 take the first rank and are placed in the first front. Then, for each solution p with the first rank, each member (q) from the set of Sp is met, and its nq is reduced to one unit. Thus, if nq=0, then the solution q takes the second rank and is placed in the second front. In the same way, the next fronts are formed. In the following, the fast nondominated sorting is presented in the form of pseudocode and figure. Figure 1 shows the nondominated sorting of a population of solutions in the MOIWO metaheuristic algorithm.


## 4.1. Epsilon Constraint Method

The Epsilon constraint method is an exact solution method for solving mathematical models with multiobjective functions. This method leads to the optimization of multiobjective mathematical models by maintaining only one objective function and transferring other objective functions to constraints. This general model is as follows.(21)Minimizefμxsubject tofmm≤εm,m=1,2,…,Mandm≠μgjx≥0,j=1,2,…,J.In the above model, the parameterεm indicates the upper limit of the value fm. To show how this method works, consider the problem with the two objectives f1 and f2. Also, suppose that f2remains in the objective function and f1 is added to the set of constraintsf1x≤ε1.Assume that initiallyε1=ε1c. The problem resulting from considering this constraint divides the main feasible objective space into two parts of f1x≤ε1c and f1x>ε1c. The left part of the feasible space obtains the feasible solution to the problem. Now, the mission of the resulting problem is to find a solution that is the minimum value of this feasible space. It is obvious from the figure that c is the minimum solution. According to this figure, intermediate solutions can be obtained in problems with nonconvex objective space using the ε-constraint method.


## 4.2. Multiobjective Invasive Weed Optimization Algorithm

The invasive weed optimization algorithm (IWO) is a nature-inspired metaheuristic algorithm that is inspired by weed behavior and its growth. They grow widely and strikingly so that they cannot be removed and controlled by humans. One claim about weeds is that weeds always win. In general, the reasons for this claim can be stated as follows:(i)
Existence of weeds after thousands of years of agriculture(ii)
Existence of weeds even after using various pesticides(iii)
The emergence of new species of weeds widely on the ground(iv)
Adapting to the environmentThe above characteristics show that weeds are vigorous and interfere with agriculture. On the other hand, weeds adapt themselves to the environment and change their behavior to grow. The success of weeds depends on their ecology and biology. In the following, optimization based on invasive weed behavior is investigated. The stages of the invasive weed optimization algorithm are as follows:(i)
First stage: spreading the seeds in the desired area(ii)
Second stage: seed growth according to utility (sprouting) and spatial dispersal(iii)
Third stage: survival of weeds with more utility (competitive elimination)(iv)
Fourth stage: continuing the process until reaching the plants with the most utilityDetails of the steps of the invasive weeds optimization algorithm are as follows:First step: generating the random initial population and evaluating their objective functionAccording to this stage, a certain number of seeds are produced from each plant. The seeds produced depend on their fitness level as well as the maximum and minimum fitness of the colony of the weeds, and this number will increase linearly with increasing the amount of fitness. On the other hand, the number of seeds produced by each plant depends on the proportion of the plant to its environment; this relationship is linear, which means that the best solution from the current population gives the highest new solution and the worst solution leads to the generation of the lowest new solution. The number of solutions varies linearly between these two limits.Second stage: reproduction based on fitness and updating of standard deviationInvasive weeds may reproduce with or without using sex cells, depending on the type of the plant. In the method of sexual reproduction, seeds or spores are used. In this case, the plant is born and begins its life by being fertilized. Then, the seeds are dispersed by various factors such as wind, water, and animals, until the seed has a chance to grow. If the conditions are good, the seeds will start to germinate and grow. The process of growth will continue with other neighbor plants until they become mature plants. In the last stage of the life cycle, they become flowering plants and cause seed production.Third stage: competitive exclusionIf a plant does not reproduce, it will die. Therefore, competition between plants is needed to limit the maximum number of plants in a colony. After producing seeds around each weed, we can only transfer the predetermined maximum number of plants (Pmax) of the total weeds and seeds to the next generation. Plants that have a chance of surviving reproduce and repeat the above stages so that the solutions obtained in each of the iterations would be fitter. This mechanism gives plants with low proportions a chance to reproduce, and if the seeds produced by them have a good proportion in the colony, then they can survive. This algorithm stops when the number of iterations reaches the maximum number allowed. The maximum number of preset plants can be equal to the number of the initial population, in which case one of the parameters of the algorithm is removed.Fourth stage: checking the termination conditionsTypes of termination conditions in metaheuristic methods are as follows:(i)
Obtaining an acceptable minimum of solutionsIn this case, which is one of the termination conditions, assume that the total cost of a company is 1000 monetary units and the manager of the company tends to reduce this cost to 800 units, which is a desirable amount from the manager’s point of view. It is noteworthy that the value of 800 is not the optimal point, and the value of the objective function may become less than 800 units. For example, in Figure1, point A shows the desired point from the manager’s point of view, which is one of the termination conditions.(ii)
The passage of time or determined iterationIn this case, a determined iteration is introduced for the algorithm to end the algorithm. For example, in iteration 100, the algorithm ends, and there would be no need for iterations after 100, although the value of the objective function may improve after the determined iteration.(iii)
The passage of time or determined iteration without observing much improvement in the resultIf, after implementing the algorithm in successive iterations, the value of the objective function does not improve or improves insignificantly, the algorithm will stop. For example, if it does not improve from pointa to b with n iterations or improves insignificantly, the algorithm will stop. In Figure 2, the value of the objective function is constant in iteration a to iteration b that causes the algorithm to stop.Figure 1
Pseudocode of the fast nondominated sorting procedure.The general structure of MOIWO and IWO is similar. The only difference between them is that, in each iteration, the nondominated solutions are separated. This method is presented in Section4.2.1. The pseudocode of the MOIWO algorithm is shown in Figure 3.Figure 2
Comparison of solving methods based on Max Spread metric.Figure 3
Pseudocode of MOIW algorithm.


### 4.2.1. Nondominated Sorting

First, two parametersnp and Sp must be calculated; np is actually the number of solutions that have dominated the solution p, and Sp is the set of solutions that have been dominated by the solution p. Then, a set of solutions can be ranked, and they can be placed on different fronts in terms of nondomination degrees. All solutions with np=0 take the first rank and are placed in the first front. Then, for each solution p with the first rank, each member (q) from the set of Sp is met, and its nq is reduced to one unit. Thus, if nq=0, then the solution q takes the second rank and is placed in the second front. In the same way, the next fronts are formed. In the following, the fast nondominated sorting is presented in the form of pseudocode and figure. Figure 1 shows the nondominated sorting of a population of solutions in the MOIWO metaheuristic algorithm.


## 4.2.1. Nondominated Sorting

First, two parametersnp and Sp must be calculated; np is actually the number of solutions that have dominated the solution p, and Sp is the set of solutions that have been dominated by the solution p. Then, a set of solutions can be ranked, and they can be placed on different fronts in terms of nondomination degrees. All solutions with np=0 take the first rank and are placed in the first front. Then, for each solution p with the first rank, each member (q) from the set of Sp is met, and its nq is reduced to one unit. Thus, if nq=0, then the solution q takes the second rank and is placed in the second front. In the same way, the next fronts are formed. In the following, the fast nondominated sorting is presented in the form of pseudocode and figure. Figure 1 shows the nondominated sorting of a population of solutions in the MOIWO metaheuristic algorithm.


## 5. Numerical Results

### 5.1. Performance Metrics for Comparison of Metaheuristic Algorithms

In this section, quantitative and qualitative metrics are introduced. The main application of these metrics is to compare the performance of metaheuristic algorithms.


#### 5.1.1. SNS

This metric, also called the spread metric, is used to calculate the degree of variation of Pareto solutions. The priority of the algorithm depends on the value of this metric. Using the following equation, the value of this metric can be calculated.(22)SNS=∑i=1nMID−Ci2n−1.In this equation,n represents the number of nondominated solutions. Using the following equation, Ci can be calculated.(23)Ci=f12+f22+f32.In the above equation,f1i and f2iare, respectively, the values of the first and second objective functions for the nondominated solution of i.


#### 5.1.2. Max Spread Metric

It is necessary to calculate the optimal Pareto front solutions obtained by the algorithm in terms of scope. The Max Spread metric is used for this. The larger the value of this metric is, the higher the priority of the algorithm will be. The value of this metric is calculated using the following equation:(24)DM=∑i=1IMinfi−Maxfi2.In the above equation,Minfi represents the minimum value of the objective function among all nondominated solutions obtained from the algorithm, and Maxfi represents the maximum value of the objective function among all nondominated solutions obtained from the algorithm.


#### 5.1.3. MID (Mean Ideal Distance)

Using this metric, the proximity among the obtained nondominated solutions and the ideal point is measured. The priority of the algorithm is inversely related to the value of this metric. The following equation is used to calculate this metric.(25)MID=∑i=1nf1i−f1best/f1totalmax−f1totalmin2+f2i−f2best/f2totalmax−f2totalmin2n.In the above equation,f1totalmax is the largest value between nondominated solutions. Also, f2totalmin is considered as the smallest value between the nondominated solutions obtained by the algorithm.


#### 5.1.4. RAS (The Rate of Achievement to Objectives Simultaneously)

This metric seeks a balance in achieving each of the objectives and its value is calculated using the following equation:(26)RAS=∑i=1nf1ix−f1ibestx/f1ibestx+f2ix−f2ibestx/f2ibestxn.The value of this metric is inversely related to the priority of the algorithm.


#### 5.1.5. Spacing

By using this metric, the uniformity of the distribution of nondominated solutions is obtained according to the following equation. The value of this metric is inversely related to the priority of the algorithm.(27)SM=∑i=1nd¯−din−1d¯.In the above equation,n is the number of obtained nondominated solutions, di is the smallest Euclidean distance between the solution i and the solutions found in the set of obtained nondominated solutions, and d¯is the average of these distances.


#### 5.1.6. NPS

This metric is used to calculate the number of nondominated solutions. The value of this metric is directly related to the priority of the algorithm.It is worth noting that the introduced metrics are presented for a problem with three objective functions. In the case of increasing the number of objective functions, it is enough to include the values ​​of the objective functions in the formula and generalize the desired formula. It should be considered that Max Spread, NPS, Spacing, SNS, RAS, and MID metrics have been used in this study to compare the efficiency of solving methods.


### 5.2. Investigation of the Efficiency of the Proposed Solving Methods

To compare the considered methods based on the studied metrics, 10 examples in different scales have been produced. These examples have been randomly produced. The information of these examples is as shown in Table2. It should be considered that the time limit of 3600 seconds is regarded for each solving method.Table 2
Information of numerical problems.
Example numberNumber of suppliersNumber of factoriesNumber of distribution centersNumber of productsNumber of raw materialsNumber of transportation systemsNumber of time periods1524111727363217394853194115107429513612952126157141162127178161373158199181583159211020179320102312221910420Then, the values of quality metrics of the problem solution are calculated for both solving methods and reported in Tables3 and 4.Table 3
Epsilon constraint output for solved examples.
Example numberMIDMax_SpreadSMNPSRASSpacingCPU time112.52211.6417.99470.3187.8846.243213.34712.0019.927100.3919.77913.923316.21914.4349.719110.63812.03180.205416.31816.05715.639110.68217.772425.727518.03717.3420.293130.89721.792994.581617.84422.86322.848151.52823.3361272.391720.00631.39530.102172.06433.2152464.165820.40540.31134.497193.10935.1373600.299921.62942.37452.884193.61649.153600.57110000000Not solvedAverage17.36923.15722.65613.5561.47123.3441384.234Table 4
MOIWO output for solved examples.
Example numberMIDMax_SpreadSNSNPSRASSpacingCPU time118.8638.12810.08370.4467.03813.628221.19816.77812.985190.5549.12613.067314.08617.49613.515170.76717.50631.58420.84924.96524.499220.78616.4557.938525.84418.69326.859201.22424.94674.544619.88433.12624.651371.57730.37396.803720.77838.49330.473402.52946.129129.095827.45847.35649.158494.22343.533134.002925.43347.81360.994473.86553.544159.8271026.68155.86260.842486.49364.762171.294Average22.10730.87131.405930.62.24731.34188.17783


#### 5.2.1. Investigation of Solving Methods Based on MID Metric

The average of this metric is 17.5 for the Epsilon constraint and 22.10 for the MOIWO algorithm. In Figure4, the value of this metric is shown for each of the solved examples by two algorithms.Figure 4
Comparison of solving methods based on MID metric.According to the results, it is concluded that the value of MOIWO was higher in all examples because using EPC, the best values are obtained. Therefore, the short distance between the metaheuristic method and these values shows the quality of this method.


#### 5.2.2. Investigation of Solving Methods Based on Max Spread Metric

The average of this metric is 23.15 for the Epsilon constraint and 31.87 for the MOIWO algorithm. Figure2 shows the value of this metric for each of the examples solved by the two solving methods.According to Figure5, the EPC method provided a better value than MOIWO only in the first example. In other examples, the MOIWO method provided a better value. That shows the efficiency of this method to find the solutions.Figure 5
Comparison of solving methods based on NPS metric.


#### 5.2.3. Investigation of the Solving Methods Based on SNS Metric

The average of this metric is 22.65 for the Epsilon constraint and 31.40 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.Figure 6
Comparison of solving methods based on SNS metric.Regarding the SNS metric, same as the Max Spread metric, the ability of the MOIWO algorithm to find optimal solutions is obvious. Following the increase in the number of problems, the superiority of the MOIWO method in terms of the SNS metric becomes more evident.


#### 5.2.4. Investigation of the Solving Methods Based on NPS Metric

The average of this metric is 13.55 for the Epsilon constraint and 30.60 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.The NPS metric is a metric that shows the priority of the MOIWO method compared to ECP. Due to the constraints which are added to the problem, the ability of the EPC method to find different solutions reduces. MOIWO algorithm can find more optimal solutions, which is the reason that this algorithm performs better in terms of this metric.


#### 5.2.5. Investigation of the Solving Methods Based on RAS Metric

The average of this metric is 1.47 for the Epsilon constraint and 2.24 for the MOIWO algorithm. Figure7 shows the value of this metric for each of the solved examples by the two solving methods.Figure 7
Comparison of solving methods based on RAS metric.As can be seen, the MOIWO algorithm provides higher values compared to the EPC. This small difference in many examples between the MOIWO algorithm and the EPC method, which is on average about 0.6, indicates the ability of the MOIWO algorithm.


#### 5.2.6. Investigation of the Solving Methods Based on the Spacing Metric

The average of this metric is 23.34 for the Epsilon constraint and 32.34 for the MOIWO algorithm. Figure8 shows the value of this metric for each of the solved examples by the two algorithms.Figure 8
Comparison of solving methods based on SNS metric.The analysis of this metric is the same as the RAS metric. The solution obtained from the EPC method has lower values, and the MOIWO algorithm provides an SNS value close to this method. The above results show the fact that this algorithm has been able to achieve the near-ideal solution.


#### 5.2.7. Investigation of the Solving Methods Based on Solving Time

Epsilon constraint has an average solving time of 1384.23 seconds. However, this metric is 88.17 for MOIWO. A comparison of the solving times of the two methods is shown in Figure9.Figure 9
Comparison of the solving time of the two methods.The solving time of the EPC method has an ascending solution time to solve with a very steep slope. This method solved problems 7 and 8 with the maximum possible time and also could not solve problem 10. On the other hand, the solving time of the MOIWO algorithm has increased with a very low slope.As a result, it is determined that the MOIWO algorithm can provide the outputs very close to the optimal solution in much less time than EPC and even perform better than EPC in some metrics. This issue well demonstrates the ability of the MOIWO algorithm to solve the problem.


## 5.1. Performance Metrics for Comparison of Metaheuristic Algorithms

In this section, quantitative and qualitative metrics are introduced. The main application of these metrics is to compare the performance of metaheuristic algorithms.


### 5.1.1. SNS

This metric, also called the spread metric, is used to calculate the degree of variation of Pareto solutions. The priority of the algorithm depends on the value of this metric. Using the following equation, the value of this metric can be calculated.(22)SNS=∑i=1nMID−Ci2n−1.In this equation,n represents the number of nondominated solutions. Using the following equation, Ci can be calculated.(23)Ci=f12+f22+f32.In the above equation,f1i and f2iare, respectively, the values of the first and second objective functions for the nondominated solution of i.


### 5.1.2. Max Spread Metric

It is necessary to calculate the optimal Pareto front solutions obtained by the algorithm in terms of scope. The Max Spread metric is used for this. The larger the value of this metric is, the higher the priority of the algorithm will be. The value of this metric is calculated using the following equation:(24)DM=∑i=1IMinfi−Maxfi2.In the above equation,Minfi represents the minimum value of the objective function among all nondominated solutions obtained from the algorithm, and Maxfi represents the maximum value of the objective function among all nondominated solutions obtained from the algorithm.


### 5.1.3. MID (Mean Ideal Distance)

Using this metric, the proximity among the obtained nondominated solutions and the ideal point is measured. The priority of the algorithm is inversely related to the value of this metric. The following equation is used to calculate this metric.(25)MID=∑i=1nf1i−f1best/f1totalmax−f1totalmin2+f2i−f2best/f2totalmax−f2totalmin2n.In the above equation,f1totalmax is the largest value between nondominated solutions. Also, f2totalmin is considered as the smallest value between the nondominated solutions obtained by the algorithm.


### 5.1.4. RAS (The Rate of Achievement to Objectives Simultaneously)

This metric seeks a balance in achieving each of the objectives and its value is calculated using the following equation:(26)RAS=∑i=1nf1ix−f1ibestx/f1ibestx+f2ix−f2ibestx/f2ibestxn.The value of this metric is inversely related to the priority of the algorithm.


### 5.1.5. Spacing

By using this metric, the uniformity of the distribution of nondominated solutions is obtained according to the following equation. The value of this metric is inversely related to the priority of the algorithm.(27)SM=∑i=1nd¯−din−1d¯.In the above equation,n is the number of obtained nondominated solutions, di is the smallest Euclidean distance between the solution i and the solutions found in the set of obtained nondominated solutions, and d¯is the average of these distances.


### 5.1.6. NPS

This metric is used to calculate the number of nondominated solutions. The value of this metric is directly related to the priority of the algorithm.It is worth noting that the introduced metrics are presented for a problem with three objective functions. In the case of increasing the number of objective functions, it is enough to include the values ​​of the objective functions in the formula and generalize the desired formula. It should be considered that Max Spread, NPS, Spacing, SNS, RAS, and MID metrics have been used in this study to compare the efficiency of solving methods.


## 5.1.1. SNS

This metric, also called the spread metric, is used to calculate the degree of variation of Pareto solutions. The priority of the algorithm depends on the value of this metric. Using the following equation, the value of this metric can be calculated.(22)SNS=∑i=1nMID−Ci2n−1.In this equation,n represents the number of nondominated solutions. Using the following equation, Ci can be calculated.(23)Ci=f12+f22+f32.In the above equation,f1i and f2iare, respectively, the values of the first and second objective functions for the nondominated solution of i.


## 5.1.2. Max Spread Metric

It is necessary to calculate the optimal Pareto front solutions obtained by the algorithm in terms of scope. The Max Spread metric is used for this. The larger the value of this metric is, the higher the priority of the algorithm will be. The value of this metric is calculated using the following equation:(24)DM=∑i=1IMinfi−Maxfi2.In the above equation,Minfi represents the minimum value of the objective function among all nondominated solutions obtained from the algorithm, and Maxfi represents the maximum value of the objective function among all nondominated solutions obtained from the algorithm.


## 5.1.3. MID (Mean Ideal Distance)

Using this metric, the proximity among the obtained nondominated solutions and the ideal point is measured. The priority of the algorithm is inversely related to the value of this metric. The following equation is used to calculate this metric.(25)MID=∑i=1nf1i−f1best/f1totalmax−f1totalmin2+f2i−f2best/f2totalmax−f2totalmin2n.In the above equation,f1totalmax is the largest value between nondominated solutions. Also, f2totalmin is considered as the smallest value between the nondominated solutions obtained by the algorithm.


## 5.1.4. RAS (The Rate of Achievement to Objectives Simultaneously)

This metric seeks a balance in achieving each of the objectives and its value is calculated using the following equation:(26)RAS=∑i=1nf1ix−f1ibestx/f1ibestx+f2ix−f2ibestx/f2ibestxn.The value of this metric is inversely related to the priority of the algorithm.


## 5.1.5. Spacing

By using this metric, the uniformity of the distribution of nondominated solutions is obtained according to the following equation. The value of this metric is inversely related to the priority of the algorithm.(27)SM=∑i=1nd¯−din−1d¯.In the above equation,n is the number of obtained nondominated solutions, di is the smallest Euclidean distance between the solution i and the solutions found in the set of obtained nondominated solutions, and d¯is the average of these distances.


## 5.1.6. NPS

This metric is used to calculate the number of nondominated solutions. The value of this metric is directly related to the priority of the algorithm.It is worth noting that the introduced metrics are presented for a problem with three objective functions. In the case of increasing the number of objective functions, it is enough to include the values ​​of the objective functions in the formula and generalize the desired formula. It should be considered that Max Spread, NPS, Spacing, SNS, RAS, and MID metrics have been used in this study to compare the efficiency of solving methods.


## 5.2. Investigation of the Efficiency of the Proposed Solving Methods

To compare the considered methods based on the studied metrics, 10 examples in different scales have been produced. These examples have been randomly produced. The information of these examples is as shown in Table2. It should be considered that the time limit of 3600 seconds is regarded for each solving method.Table 2
Information of numerical problems.
Example numberNumber of suppliersNumber of factoriesNumber of distribution centersNumber of productsNumber of raw materialsNumber of transportation systemsNumber of time periods1524111727363217394853194115107429513612952126157141162127178161373158199181583159211020179320102312221910420Then, the values of quality metrics of the problem solution are calculated for both solving methods and reported in Tables3 and 4.Table 3
Epsilon constraint output for solved examples.
Example numberMIDMax_SpreadSMNPSRASSpacingCPU time112.52211.6417.99470.3187.8846.243213.34712.0019.927100.3919.77913.923316.21914.4349.719110.63812.03180.205416.31816.05715.639110.68217.772425.727518.03717.3420.293130.89721.792994.581617.84422.86322.848151.52823.3361272.391720.00631.39530.102172.06433.2152464.165820.40540.31134.497193.10935.1373600.299921.62942.37452.884193.61649.153600.57110000000Not solvedAverage17.36923.15722.65613.5561.47123.3441384.234Table 4
MOIWO output for solved examples.
Example numberMIDMax_SpreadSNSNPSRASSpacingCPU time118.8638.12810.08370.4467.03813.628221.19816.77812.985190.5549.12613.067314.08617.49613.515170.76717.50631.58420.84924.96524.499220.78616.4557.938525.84418.69326.859201.22424.94674.544619.88433.12624.651371.57730.37396.803720.77838.49330.473402.52946.129129.095827.45847.35649.158494.22343.533134.002925.43347.81360.994473.86553.544159.8271026.68155.86260.842486.49364.762171.294Average22.10730.87131.405930.62.24731.34188.17783


### 5.2.1. Investigation of Solving Methods Based on MID Metric

The average of this metric is 17.5 for the Epsilon constraint and 22.10 for the MOIWO algorithm. In Figure4, the value of this metric is shown for each of the solved examples by two algorithms.Figure 4
Comparison of solving methods based on MID metric.According to the results, it is concluded that the value of MOIWO was higher in all examples because using EPC, the best values are obtained. Therefore, the short distance between the metaheuristic method and these values shows the quality of this method.


### 5.2.2. Investigation of Solving Methods Based on Max Spread Metric

The average of this metric is 23.15 for the Epsilon constraint and 31.87 for the MOIWO algorithm. Figure2 shows the value of this metric for each of the examples solved by the two solving methods.According to Figure5, the EPC method provided a better value than MOIWO only in the first example. In other examples, the MOIWO method provided a better value. That shows the efficiency of this method to find the solutions.Figure 5
Comparison of solving methods based on NPS metric.


### 5.2.3. Investigation of the Solving Methods Based on SNS Metric

The average of this metric is 22.65 for the Epsilon constraint and 31.40 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.Figure 6
Comparison of solving methods based on SNS metric.Regarding the SNS metric, same as the Max Spread metric, the ability of the MOIWO algorithm to find optimal solutions is obvious. Following the increase in the number of problems, the superiority of the MOIWO method in terms of the SNS metric becomes more evident.


### 5.2.4. Investigation of the Solving Methods Based on NPS Metric

The average of this metric is 13.55 for the Epsilon constraint and 30.60 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.The NPS metric is a metric that shows the priority of the MOIWO method compared to ECP. Due to the constraints which are added to the problem, the ability of the EPC method to find different solutions reduces. MOIWO algorithm can find more optimal solutions, which is the reason that this algorithm performs better in terms of this metric.


### 5.2.5. Investigation of the Solving Methods Based on RAS Metric

The average of this metric is 1.47 for the Epsilon constraint and 2.24 for the MOIWO algorithm. Figure7 shows the value of this metric for each of the solved examples by the two solving methods.Figure 7
Comparison of solving methods based on RAS metric.As can be seen, the MOIWO algorithm provides higher values compared to the EPC. This small difference in many examples between the MOIWO algorithm and the EPC method, which is on average about 0.6, indicates the ability of the MOIWO algorithm.


### 5.2.6. Investigation of the Solving Methods Based on the Spacing Metric

The average of this metric is 23.34 for the Epsilon constraint and 32.34 for the MOIWO algorithm. Figure8 shows the value of this metric for each of the solved examples by the two algorithms.Figure 8
Comparison of solving methods based on SNS metric.The analysis of this metric is the same as the RAS metric. The solution obtained from the EPC method has lower values, and the MOIWO algorithm provides an SNS value close to this method. The above results show the fact that this algorithm has been able to achieve the near-ideal solution.


### 5.2.7. Investigation of the Solving Methods Based on Solving Time

Epsilon constraint has an average solving time of 1384.23 seconds. However, this metric is 88.17 for MOIWO. A comparison of the solving times of the two methods is shown in Figure9.Figure 9
Comparison of the solving time of the two methods.The solving time of the EPC method has an ascending solution time to solve with a very steep slope. This method solved problems 7 and 8 with the maximum possible time and also could not solve problem 10. On the other hand, the solving time of the MOIWO algorithm has increased with a very low slope.As a result, it is determined that the MOIWO algorithm can provide the outputs very close to the optimal solution in much less time than EPC and even perform better than EPC in some metrics. This issue well demonstrates the ability of the MOIWO algorithm to solve the problem.


## 5.2.1. Investigation of Solving Methods Based on MID Metric

The average of this metric is 17.5 for the Epsilon constraint and 22.10 for the MOIWO algorithm. In Figure4, the value of this metric is shown for each of the solved examples by two algorithms.Figure 4
Comparison of solving methods based on MID metric.According to the results, it is concluded that the value of MOIWO was higher in all examples because using EPC, the best values are obtained. Therefore, the short distance between the metaheuristic method and these values shows the quality of this method.


## 5.2.2. Investigation of Solving Methods Based on Max Spread Metric

The average of this metric is 23.15 for the Epsilon constraint and 31.87 for the MOIWO algorithm. Figure2 shows the value of this metric for each of the examples solved by the two solving methods.According to Figure5, the EPC method provided a better value than MOIWO only in the first example. In other examples, the MOIWO method provided a better value. That shows the efficiency of this method to find the solutions.Figure 5
Comparison of solving methods based on NPS metric.


## 5.2.3. Investigation of the Solving Methods Based on SNS Metric

The average of this metric is 22.65 for the Epsilon constraint and 31.40 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.Figure 6
Comparison of solving methods based on SNS metric.Regarding the SNS metric, same as the Max Spread metric, the ability of the MOIWO algorithm to find optimal solutions is obvious. Following the increase in the number of problems, the superiority of the MOIWO method in terms of the SNS metric becomes more evident.


## 5.2.4. Investigation of the Solving Methods Based on NPS Metric

The average of this metric is 13.55 for the Epsilon constraint and 30.60 for the MOIWO algorithm. Figure6 shows the value of this metric for each of the examples solved by the two solving methods.The NPS metric is a metric that shows the priority of the MOIWO method compared to ECP. Due to the constraints which are added to the problem, the ability of the EPC method to find different solutions reduces. MOIWO algorithm can find more optimal solutions, which is the reason that this algorithm performs better in terms of this metric.


## 5.2.5. Investigation of the Solving Methods Based on RAS Metric

The average of this metric is 1.47 for the Epsilon constraint and 2.24 for the MOIWO algorithm. Figure7 shows the value of this metric for each of the solved examples by the two solving methods.Figure 7
Comparison of solving methods based on RAS metric.As can be seen, the MOIWO algorithm provides higher values compared to the EPC. This small difference in many examples between the MOIWO algorithm and the EPC method, which is on average about 0.6, indicates the ability of the MOIWO algorithm.


## 5.2.6. Investigation of the Solving Methods Based on the Spacing Metric

The average of this metric is 23.34 for the Epsilon constraint and 32.34 for the MOIWO algorithm. Figure8 shows the value of this metric for each of the solved examples by the two algorithms.Figure 8
Comparison of solving methods based on SNS metric.The analysis of this metric is the same as the RAS metric. The solution obtained from the EPC method has lower values, and the MOIWO algorithm provides an SNS value close to this method. The above results show the fact that this algorithm has been able to achieve the near-ideal solution.


## 5.2.7. Investigation of the Solving Methods Based on Solving Time

Epsilon constraint has an average solving time of 1384.23 seconds. However, this metric is 88.17 for MOIWO. A comparison of the solving times of the two methods is shown in Figure9.Figure 9
Comparison of the solving time of the two methods.The solving time of the EPC method has an ascending solution time to solve with a very steep slope. This method solved problems 7 and 8 with the maximum possible time and also could not solve problem 10. On the other hand, the solving time of the MOIWO algorithm has increased with a very low slope.As a result, it is determined that the MOIWO algorithm can provide the outputs very close to the optimal solution in much less time than EPC and even perform better than EPC in some metrics. This issue well demonstrates the ability of the MOIWO algorithm to solve the problem.


## 6. Conclusion

In this study, a multiobjective mathematical model with a green and agile approach was presented. The integration of agility and environmental concepts in supply chain design is an important innovation of research. The use of a new multiobjective metaheuristic algorithm called MOIWO is the next innovation of this research. The results showed that economic, environmental, and agility objectives in the supply chain are in conflict with each other, and their independent optimization is not efficient enough to be used in the supply chain. Also, the MOIWO algorithm has performed well compared with the Epsilon constraint method. This algorithm has been better than the Epsilon constraint method in terms of NPS, SNS, and Max Spread metrics. In terms of the solving time, the average solving time of this algorithm was about 0.1% of the solving time of the Epsilon constraint method. Therefore, the efficiency of this multiobjective algorithm has been confirmed by various metrics, and this algorithm can be introduced as an appropriate tool for solving supply chain design problems. To develop this research and for future studies, it is suggested that uncertainty be considered in the model parameters such as customer demand, manufacturing, and distribution costs in the mathematical model. This uncertainty can be implemented with approaches such as possibilistic programming or robust optimization and the results can be examined. It is also recommended to use another metaheuristic algorithm and compare them with each other. To demonstrate the performance of the proposed mathematical model, it is suggested that the model be implemented in the real world and the results are reported.



---

*Source: 3811417-2021-11-02.xml*