Weighing Efficiency-Robustness in Supply Chain Disruption by Multi-Objective Firefly Algorithm

Abstract

This paper investigates various supply chain disruptions in terms of scenario planning, including node disruption and chain disruption; namely, disruptions in distribution centers and disruptions between manufacturing centers and distribution centers. Meanwhile, it also focuses on the simultaneous disruption on one node or a number of nodes, simultaneous disruption in one chain or a number of chains and the corresponding mathematical models and exemplification in relation to numerous manufacturing centers and diverse products. Robustness of the design of the supply chain network is examined by weighing efficiency against robustness during supply chain disruptions. Efficiency is represented by operating cost; robustness is indicated by the expected disruption cost and the weighing issue is calculated by the multi-objective firefly algorithm for consistency in the results. It has been shown that the total cost achieved by the optimal target function is lower than that at the most effective time of supply chains. In other words, the decrease of expected disruption cost by improving robustness in supply chains is greater than the increase of operating cost by reducing efficiency, thus leading to cost advantage. Consequently, by approximating the Pareto Front Chart of weighing between efficiency and robustness, enterprises can choose appropriate efficiency and robustness for their longer-term development.

1. Introduction

2014 witnessed the epidemic Ebola, affecting Cote d’Ivoire, the main supplier of cocoa to the Nestle Corporation. A crisis of raw materials occurred in supply chains, leading to a huge loss for the Nestle Corporation. The National Seismic Center in Chile announced that an earthquake of 8.2 on the Richter scale took place in the northern Chile and the price of bronze increased dramatically, which brought substantial pressure related to the cost of materials to producers of sanitary appliances and metals. On 27 January 2015, heavy snowfall struck Huoshan County in China and trees fell down on the high-voltage power lines, resulting in blackout in 45 base stations and in the majority of places in the county.

Extant literature has described huge damages caused by supply chain disruptions. For instance, Li and Liu (2014) [1] proposes an analytic algorithm to calculate the purchasing time and inventory based on the sensitivity model, which does not consider the uncertainty of demand and supply among supply chain members. Mizgier et al. (2015) [2] investigate whether or not capital charge could be combined with process improvement, an approach predominantly employed for managing high-frequency, low-impact operational disruptions. They note that process improvement can be a first line of defense to manage certain types of operational risk events, but their study is not without limitations. Their categorization of event types is not developed in the literature but through regulation, and it is actively followed in practice. In addition, the event types are broad and can be made granular. Constrained by the order quantity, Gorji et al. (2014) [3] discuss the optimal inventory decision-making related to two-echelon supply chains, assuming that clients would wait for new purchase in shortage, but they do not note that suppliers have limited transportation facilities. Mizgier et al. (2015) [4] propose a model for the calculation of the loss distribution from disruptions in a SCN. This gives supply chain risk managers an efficient tool for the quantification of supply chain risks intrinsically embedded in the existing network design, but they concentrate only on the downstream risk arising from supplier disruptions. Wang et al. (2013) [5] present a systematic cost model to facilitate manufacturers to evaluate and increase the minimal distribution cost of the supplier quality, whereas they do not take the interaction between suppliers into account. When demand is disrupted, Cao, et al. (2013) [6] develop a coordinated mechanism for a one manufacturer and n Gounod competitive retailers supply chain; in reality, however, there are multiple manufacturers simultaneously. Hishamuddin et al. (2014) [7] consider a recovery mechanism for disruptions in two-echelon supply chains, which only involves two-echelon supply chain model with one manufacturer and one retailer. Shu et al. (2014) [8] investigate risk control of production disruption in supply chains on the basis of Generic Bill of Materials (GBOM), and study the strategies related to the overall optimal profits achieved by enterprises when production is uncertain. They make a hypothesis regarding corporate efficiency and the disruptions of market demand, and it is difficult to obtain the real and accurate disruption probability as a result of the uncertainty of supply chain disruptions. The studies mentioned above do not consider the robustness of supply chains; however, efficiency can be increased in addressing supply chain disruption if robustness is introduced in the design of supply chain network. Here, the robustness of supply chains refers to the robust supply chain. The more robust the supply chain is, the less likely the supply chain tends to be disrupted. As such, robustness of supply chains is added as a key consideration here. The introduction of robustness of supply chains involves weighing efficiency and robustness of enterprises, to optimize multiple objectives. Rosenberg (1967) [9] applies a genetic algorithm to address the problem of multiple objectives and since then swarm intelligence algorithm has been widely used to tackle multi-objective problems. Mizgier et al. (2015) [10] formulate a multi-objective optimization model for the capital allocation problem and apply it to a real-world case of two financial conglomerates. They extend the optimum capital allocation problem to a multi-objective optimization problem by incorporating the cost of capital. Nevertheless, their choice of parameters has only partially been validated with real data. Yang (2013) [11] proposes the firefly algorithm and since then a number of scholars have improved it and applied it in all fields. The firefly algorithm mimics some characteristics of tropic firefly swarms and their flashing behavior. Sayadi et al. (2010) [12] adopt the firefly algorithm to study the minimal completion time for the scheduling problem in assembly lines, but the setup of parameters is not discussed sufficiently. Yang (2008) [13] extends the firefly algorithm to address multi-objective problems and tests the validity of the algorithm by functions and applies it for designing benchmark optimization. In addition, they apply the firefly algorithm to discreet optimization of manufacturing units (Sayadi et al., 2013; Sayadi et al., 2010) [12,14], involving optimization with a single objective. Chandrasekaran et al. (2013) [15] solve the problem of unit commitment by means of the firefly algorithm with a real-number encoding binary system, without considering the optimization of multiple conflicting objectives. Talatahari et al. (2014) [16] optimize a tower structure with the firefly algorithm, whereas they fail to conceive of the case of discreteness. Marichelvam et al. (2014) [17] address the scheduling problem of multi-objective mixed assembly lines by the discreet firefly algorithm; however, they do not consider that machines may be irrelevant or different in each stage.

The majority of studies on supply chains have considered only one product and one manufacturing center, which is ideal. In reality, most enterprises involve more than one product and have a number of manufacturing centers in different regions. Our study, therefore, suits the supply chain with multiple products and manufacturing centers. The probability of supply chain disruption is difficult to quantify and most studies consider strategies for addressing the aftermath of the disruption, which is a passive approach. This paper considers the robustness of supply chains at the supply chain design stage itself; at the same time, node disruption and link disruption of supply chains are also taken into account. In addition, efficiency and robustness of supply chains are weighed, providing evidence for enterprises to choose the appropriate efficiency and robustness of supply chains. The multi-objective firefly algorithm involves maximal efficiency and robustness of supply chains simultaneously, which is different from the discreet points from single-objective algorithms. The multi-objective firefly algorithm produces the weighing line graph for approximating the Pareto Front, whose validity has been proved and applied in each field.

This paper, for the first time, introduces the issue of weighing efficiency and robustness of supply chains disrupted by the multi-objective firefly algorithm, which is deployed to solve the weighting problem. The result is the continuous Pareto curve instead of a number of independent points, and the robustness and cost of supply chains corresponding to their efficiency can be derived. The results show that the total cost with the optimal target function is lower than that with the most effective supply chain. In other words, the reduction of anticipated disruption cost caused by the increase of supply chain robustness is greater than the rise of operation cost caused by the drop of the efficiency of supply chains, thus creating the cost advantage. This has both theoretical and practical significance in this competitive age. Therefore, enterprises can choose appropriate efficiency and robustness to promote their long-term development based on the approximate Pareto front figure of weighing efficiency and robustness.

The organization of this paper is as follows: the supply chain model is first constructed; the firefly algorithm and the multi-objective firefly algorithm are explicated in detail; the validity of the multi-objective firefly algorithm is tested; the multi-objective firefly algorithm is exemplified and simulated; node disruption and link disruption are considered at the same time; and the simulated results are finally analyzed.

2. Building Models

The previous research on supply chains has mostly considered only one product and one manufacturing center, as Shukla et al. (2011) [18] do. This design does not fit the reality, and thus this paper considers supply chains with multiple products and manufacturing centers. The majority of enterprises have manufacturing centers for their products, delivering goods to their distribution centers and client regions where goods are in demand, and the three-echelon model has been adopted by many researchers (Shukla et al., 2011; Meepetchdee and Shah, 2007; Peng et al., 2011) [18,19,20]. As a result, this paper designs the three-echelon model with multiple products and manufacturing centers. The supply chain comprises the fixed manufacturing centers a with multiple products, the potential distribution centers b and fixed client zones c, as shown in Figure 1. In this supply chain, the varieties and quantities of products in multiple manufacturing centers rest on the model; the numbers, locations, varieties and quantities of delivered goods depend on the model; and the construction of the distribution center leads to the corresponding cost of the infrastructure. The client zone represents the demand for one or multiple products. The operational cost includes infrastructure cost, manufacturing cost and material carrying cost in distribution centers and transport cost. The scenario planning is adopted to calculate and analyze the anticipated disruption cost in different cases of supply chain disruption. The scenarios can be the node disruption in the manufacturing center or distribution center, or the link disruption from manufacturing center to the distribution center or from distribution center to the client zone. To ensure stable operations in the overall supply chain system, once the distribution center is constructed, it must serve some client zones. Each distribution center and each client zone can be supplied by only one manufacturing center and one distribution center. There is no inventory accumulation or loss. The demand of each client zone can be satisfied. The data given consist of the demand of client zones, cost, manufacturing center and the distance between distribution centers and client zones, the probability of disruption in each scenario and the quantity of disrupted products. The aim is to maximize the two conflicting targets of efficiency and robustness of supply chains, and weigh and analyze efficiency and robustness and compare the total cost at the same time.

Tsiakis et al. (2001) [21], Meepetchdee and Shah (2007) [19], Shukla et al. (2011) [18], Peng et al. (2011) [20] and Kalaitzidou et al. (2014) [22] present mathematical models of robustness of supply chains in their designs. Mixed integer linear planning is most common in the design of supply chains, for addressing node and link disruptions of supply chains. Parameter setting contains index, decision-making variables, demand parameters, efficiency parameters, cost parameters, distance parameters, probability parameters and disrupted product quantity parameters. In addition, the decision-making variables can be divided into two categories. In the first category, binary system variables define whether the manufacturing center, the distribution center and client zones are connected. In the second category, the integer variable defines the product quantity.

2.1. Parameters Setting

• Index

$a$: manufacturing centers; $b$: distribution centers; $c$: client zones; $N$: product types; $s$: scenario collections.

• Decision-making variable

Define whether manufacturing centers supply distribution centers as the binary system variables.

$$V\left(a,b\right)\{\begin{array}{c}1, manufacturing centers a supply distribution centers b\\ 0,manufacturing centers a do not supply distribution centers b \end{array}$$

Define whether distribution centers supply client zones as the binary system variables.

$$V\left(b,c\right)\{\begin{array}{c}1, distribution centers b supply client zones c\\ 0,distribution centers b do not supply client zones c \end{array}$$

Define the description of multi-echelon supply chains as integer variables.

$R\left(N,a\right)$: quantity of product $N$ by the manufacturing center $a$; $Q\left(N,a,b\right)$: quantity of product $N$ from the manufacturing center $a$ to the distribution cente $b$; $Q\left(N,b,c\right)$: quantity of product from the distribution center $b$ to the client zone $c$.

• Demand parameter

$D\left(N,c\right)$: annual demand of the product $N$ in client zone $c$.

• Efficiency parameter

${\eta }_e$: efficiency of supply chains; ${\eta }_r$: robustness of supply chains

• Cost parameter

$c\left(b,f\right)$: the fixed cost amortized annually when distribution center $b$ is constructed; $C\left(N,b,v\right)$: the unit changeable cost of product $N$ amortized annually when the distribution center $b$ is constructed; $C\left(N,a\right)$: the production cost of the unit product $N$ in the manufacturing center $a$; $C\left(N,b,h\right)$: unit carrying cost of the product $N$ in the distribution center $b$; $C\left(N,a,b\right)$: unit shipment cost of the unit product $N$ from the manufacturing center $a$ to the distribution center $b$; $C\left(N,b,c\right)$: the unit shipment cost of the unit product $N$ from the distribution center $b$ to the client zone $c$; $u\left(N\right)$: opportunity cost, i.e., unit marginal profit of Product $N$.

• Distance parameter

$G\left(a,b\right)$: the distance from the manufacturing center $a$ to the distribution center $b$; $\mathrm{H}\left(b,c\right)$: the distance from the distribution center $b$ to the client zone $c$.

• Probability parameter

$p\left(s\right)$: probability of scenario $s$

• Disrupted product quantity parameter

$q\left(N,s\right)$: quantity of disrupted product $N$ in scenario $s$

2.2. Constraints

• Network structure constraint

All the relevant network structure constraints in the manufacturing center, the distribution center and the client zone can be summarized as follows.

$$V\left(a,b\right)\le {\displaystyle \sum }_{c}V\left(b,c\right) \forall a,b$$

(1)

$$V\left(b,c\right)\le {\displaystyle \sum }_{a}V\left(a,b\right) \forall a,b,c$$

(2)

$$\mathrm{Q}\left(\mathrm{N},\mathrm{a},\mathrm{b}\right)\le \mathrm{k}\times \mathrm{V}\left(\mathrm{a},\mathrm{b}\right)\forall \mathrm{N},\mathrm{a},\mathrm{b}$$

(3)

$$Q\left(N,b,c\right)\le k\times V\left(b,c\right) \forall N,b,c$$

(4)

$${\displaystyle \sum }_{a}V\left(a,b\right)\le 1 \forall b$$

(5)

$${\displaystyle \sum }_{b}V\left(b,c\right)=1 \forall c$$

(6)

Formula (1) shows that if the manufacturing center $a$ serves distribution center $b$, then the distribution center $b$ supplies at least some client zone. Formula (2) shows that if the distribution center $b$ is built, then the client zone $c$ may or may not be supplied by distribution center $b$. As long as manufacturing center $a$ supplies distribution center $b$, manufacturing center $a$ can provide distribution center $b$ with product $N$. Thus, constraint (3) is formed, and $k$ is the appropriate large number, with k = 1,000,000,000. The same constraint can be applied to distribution center $b$ and the client zone $c$, as shown in Formula (4). Formulas (5) and (6) are constraints with a single source to ensure that each distribution center and each client zone can be supplied by one manufacturing center and one distribution center.

• Material balance constraint

If there is no inventory accumulation and loss, the material balance constraint can be summarized as follows.

$$R\left(N,a\right)={\displaystyle \sum }_{b}Q\left(N,a,b\right) \forall N,a$$

(7)

$${\displaystyle \sum }_{a}Q\left(N,a,b\right)={\displaystyle \sum }_{c}Q\left(N,b,c\right) \forall N,b$$

(8)

$${\displaystyle \sum }_{b}Q\left(N,b,c\right)=D\left(N,c\right) \forall N,c$$

(9)

Formula (7) shows that the quantity of product $N$ from manufacturing center $a$ to the distribution center $b$ amounts to the quantity of product $N$ in manufacturing center $a$. Likewise, the quantity of product $N$ from manufacturing center $a$ to distribution center $b$ amounts to the quantity of product $N$ from distribution center $b$ to client zone $c$, as shown in Formula (8). Formula (9) ensures that the demand of each client zone can be satisfied.

• Non-negativity constraints

All the consecutive variables must be non-negativity.

$$R\left(N,a\right)\ge 0 \forall N,a$$

(10)

$$Q\left(N,a,b\right)\ge 0 \forall N,a,b$$

(11)

$$Q\left(N,b,c\right)\ge 0 \forall N,b,c$$

(12)

To reduce the search space efficiently, efficiency and robustness of supply chains must be non-negative.

$${\eta }_e\ge 0$$

(13)

$${\eta }_r\ge 0$$

(14)

2.3. Objective Functions

The construction of supply chain considers efficiency and robustness, and thus the objective targets are defined as the two conflicting targets of the maximal efficiency and robustness. The efficiency of supply chains is expounded in terms of operational cost, whereas robustness of supply chains is explicated in terms of the anticipated disruption cost.

$$max{\eta }_e=\frac{c{\left(O\right)}_{max}-c\left(O\right)}{c{\left(O\right)}_{max}-c{\left(O\right)}_{min}}$$

(15)

$$max{\eta }_r=\frac{c{\left(E\right)}_{max}-c\left(E\right)}{c{\left(E\right)}_{max}-c{\left(E\right)}_{min}}$$

(16)

$c{\left(O\right)}_{max}$: the operations cost in most robust supply chains; $c{\left(O\right)}_{min}$: the operations cost in the most effective supply chains; $c{\left(E\right)}_{min}$: the anticipated disruption cost in the most robust supply chains; $c{\left(E\right)}_{max}$: the anticipated disruption cost in the most effective supply chains.

The operations cost of the objective functions include the infrastructure cost, production cost, the material carrying cost and shipment cost in the distribution center.

$${\displaystyle \sum }_{a,b}V\left(a,b\right)\times c\left(b,f\right)+{\displaystyle \sum }_{b,N}(C\left(N,b,v\right)\times {\displaystyle \sum }_{a}Q\left(N,a,b\right))$$

(17)

$${\displaystyle \sum }_{N,a}R\left(N,a\right)\times C\left(N,a\right)$$

(18)

$${\displaystyle \sum }_{N,b,c}Q\left(N,b,c\right)\times C\left(N,b,h\right)$$

(19)

$${\displaystyle \sum }_{N,a,b}Q\left(N,a,b\right)\times G\left(a,b\right)\times C\left(N,a,b\right)$$

(20)

$${\displaystyle \sum }_{N,b,c}Q\left(N,b,c\right)\times H\left(b,c\right)\times C\left(N,b,c\right)$$

(21)

Formula (17) is the infrastructure cost incurred for the construction of the distribution center, which is a fixed cost. However, the changeable cost is the annual amortization, related to the unit changeable cost of product $N$ in the distribution center, multiplied by the quantity. If the production cost in the manufacturing center is in proportion to the quantity of products, the total production cost in the manufacturing center is as shown in Formula (18). Formula (19) shows that the material carrying cost in the distribution center is in proportion to the total handling capacity of the distribution center. Formulas (20) and (21) show the shipment cost from the manufacturing center to the distribution center, and from the distribution center to the client zone, respectively. The shipment cost is the function of product quantity, distance and unit shipment cost. Generally, trucks tend to carry the full cargo and therefore the economy of scale effects on the shipment cost are neglected here.

As discussed above, the operation cost can be expressed as follows.

$$c(O)={\displaystyle \sum }_{a,b}V\left(a,b\right)\times c\left(b,f\right)+{\displaystyle \sum }_{b,N}(C(N,b,v)\times {\displaystyle \sum }_{a}Q\left(N,a,b\right))+{\displaystyle \sum }_{N,a}R\left(N,a\right)\times C\left(N,a\right)+{\displaystyle \sum }_{N,b,c}Q\left(N,b,c\right)\times C\left(N,b,h\right)+{\displaystyle \sum }_{N,a,b}Q\left(N,a,b\right)\times G\left(a,b\right)\times C(N,a,b)+{\displaystyle \sum }_{N,b,c}Q\left(N,b,c\right)\times H\left(b,c\right)\times C\left(N,b,c\right)$$

(22)

The anticipated disruption cost of objective functions can be defined by the scenario methods, which is an ancient concept. In the earliest records, people were already interested in the future and considered scenario methods as the indirect approach to explore the future society and the system (Bradfielda et al., 2005) [23]. Greiner et al. (2014) [24] investigate the support of industrial strategic decisions by means of scenario planning, but they do not subdivide the market. Kang et al. (2014) [25] consider the infrastructure of water supply by scenario planning; however, they do not consider the environmental and social impacts. Therefore, these factors should be included in the setup of the control variables. Menezes et al. (2014) [26] design the production plans for petroleum refining by scenario planning, whereas they do not note that the government may delay the investment. Tsiakis et al. (2001) [27] examine the supply chain network in uncertain circumstances for minimizing the infrastructure and operation cost by scenario planning, while they only investigate the case of demand disruption.

This paper calculates and analyzes the anticipated disruption cost in different disruption scenarios of supply chains, and the scenarios are concerned with node disruption in the manufacturing center and the distribution center or the link disruption from the manufacturing center to the distribution center or from the distribution center to the client zone. The anticipated disruption cost can be expressed by multiplication of the probability in the scenario $s$, the disrupted quantity of the product $N$, and the unit marginal profit of the product $N$.

$$c\left(E\right)={\displaystyle \sum }_{s,N}p\left(s\right)\times q\left(N,s\right)\times u\left(N\right)$$

(23)

Here, $p\left(s\right)$ indicates the probability of scenario $s$; $u\left(N\right)$ indicates the opportunity cost, that is, the unit marginal profit of product $N$; $q\left(N,s\right)$ shows the quantity of product $N$ disrupted in scenario $s$.

$$q\left(N,s\right)={\displaystyle \sum }_{N,c}\left[D\left(N,c\right)-{\displaystyle \sum }_{b\notin s}Q\left(N,b,c\right)\right]$$

(24)

Therefore, substituting Formula (24) into Formula (23), the expression of the anticipated disruption cost is as follows.

$$\begin{array}{c}c\left(E\right)={\displaystyle \sum }_{s,N,c}p\left(s\right)\times \left[D\left(N,c\right)-{\displaystyle \sum }_{b\notin s}Q\left(N,b,c\right)\right]\times u\left(N\right)\hfill \\ ={\displaystyle \sum }_{s}p\left(s\right)\times \left\{{\displaystyle \sum }_{N,c}\left[D\left(N,c\right)-{\displaystyle \sum }_{b\notin s}Q\left(N,b,c\right)\right]\times u\left(N\right)\right\}\hfill \\ ={\displaystyle \sum }_{s}p\left(s\right)\times {\displaystyle \sum }_{N,c}D\left(N,c\right)\times u\left(N\right)-{\displaystyle \sum }_{s}p\left(s\right)\times {\displaystyle \sum }_{N,c,b\notin s}Q\left(N,b,c\right)\times u\left(N\right)\hfill \end{array}$$

(25)

3. Multi-objective Firefly Algorithm

The heuristic algorithm is quite powerful in solving optimization problems, which has been well researched in numerous papers and monographs (Coello, 1999; Deb, 1999; Geem, 2009; Talbi, 2009) [28,29,30,31,32]. Most heuristic algorithms are based on the so-called swarm intelligence. A typical example is Particle Swarm Optimization (PSO), which imitates some features of swarms of birds and fish. In contrast, the firefly algorithm is different from PSO and can have two advantages: local attractions and automatic regrouping. This latter advantage makes it particularly suitable for multimodal global optimization problems (Yang, 2009; 2010) [33,34]. In our study, the multi-objective firefly algorithm considers the two conflicting targets of the maximal efficiency and robustness of supply chains, which is different from the discreet points achieved by the single-objective algorithm which avoids local optimization.

3.1. Firefly Algorithm

Yang (2008–2014) [11,13,17,33,34,35,36,37,38] proposed the firefly algorithm which was later improved, on the basis of the behavioral features of ideal glittering fireflies. Firstly, fireflies can attract other fireflies, male or female. Secondly, the degrees of attraction are in proportion to the brightness which declines as the distance increases. As a result, between two glittering fireflies, the less bright firefly moves to the brighter one. Any given firefly, without a brighter one, moves randomly. Finally, the degrees of brightness depend on the value of the objective functions (Marichelvam et al., 2014; Yang, 2008, 2009, 2010, 2012, 2013; Yang and Deb, 2010; Yang et al., 2012) [11,13,17,33,34,35,36,37,38].

In terms of maximization, the brightness of fireflies can be simply defined as proportionate with the value of objective functions. In the firefly algorithm, the relative brightness and attraction of fireflies affect their movements, and thus it is necessary to make a definition. For simplicity, if the attraction of fireflies rests on their brightness, the objective functions determine their brightness.

Conversely, the attraction of fireflies is connected with the distance between fireflies. As the distance varies, the attraction of fireflies can be defined as follows.

$$\beta ={\beta }_0{e}^{-\gamma r_{ij}^2}$$

(26)

Here, ${\beta }_0$ shows the attraction at $r=0$, that is, the maximal attraction at the light source; $\gamma$ shows the absorption coefficient of light intensity, simulating the feature of light damping in the air, which can be constant; $r_{ij}$ shows the Descartes distance between any two fireflies $i$ and $j$ at $x_i$ and $x_j$; $r_{ij}=\left|\right|x_i-x_j\left|\right|$. It is important to note that $r$ defined above does not confine to Euclidean distance; any measurement solving optimization effectively can be viewed as $r$.

The location where firefly $i$ is attracted by the brighter firefly $j$ and moves to $j$ can be re-defined as follows.

$$x_i^{t+1}=x_i^t+{\beta }_0{e}^{-\gamma r_{ij}^2}\left(x_j^t-x_i^t\right)+\alpha \left(rand-1/2\right)$$

(27)

In the Formula, $\alpha$ is the step-length factor, which is constant; $rand$ is the random factor equally distributed at [0,1]. The second part of the Formula is produced by attraction, and the third part is the random disturbance term to avoid plunging into the local optimization in advance.

The process of optimizing algorithm is as follows. The firefly groups are randomly distributed in the search space. With different locations, each firefly emits different lights. Comparing the values of the objective functions, the less bright fireflies are attracted by and moves towards the brighter ones. The distance depends on the degree of attraction shown in Formula (26). Then, the location is adjusted according to Formula (27). After many moves, all the fireflies tend to gather in the brightest location and thus achieve optimization.

3.2. Multi-Objective Optimization

In terms of optimization in the management science and engineering, the targets of optimization are more than one object, and hence the multi-objective functions are not uncommon. In the multi-objective optimization process, each sub-objective tends to be in conflict. If the performance of one sub-objective increases, another sub-objective’s performance might decrease. Consequently, it is impossible for all the sub-objectives to achieve optimization. There is compromise and coordination among these sub-objectives. Different from the optimization in the single-objective algorithm, the results might be a collection of compromised solutions, which are called Pareto optimized solutions (Coello, 1999) [28,29]. The multi-objective optimization can be described in terms of the following definitions (Pareto, 1897) [39].

Definition 1, Pareto dominance: if the conditions $\forall x\in \left\{1, 2,\dots ,m\right\}$, $u_x\le v_x$, and $\exists y\in \left\{1, 2, \dots , m\right\}$, $u_y<v_y$, are satisfied, then the vector quantity $\overrightarrow{v}=\left(v_1,v_2,\dots ,v_m\right)$ Pareto dominates another vector $\overrightarrow{u}=\left(u_1,u_2,\dots ,u_m\right)$, marked as

$$\overrightarrow{v}\prec \overrightarrow{u}$$

(28)

Definition 2, Pareto optimized solution: if and only if the conditions $\exists {\overrightarrow{x}}_u\in X$, $\overrightarrow{u}\prec \overrightarrow{v}$; here, $\overrightarrow{v}=\overrightarrow{f}\left({\overrightarrow{x}}_v\right)=\left(v_1,v_2,\dots ,v_m\right)$; and $\overrightarrow{u}=\overrightarrow{f}\left({\overrightarrow{x}}_u\right)=\left(u_1,u_2,\dots ,u_m\right)$ are satisfied, then ${\overrightarrow{x}}_v\in X$ is the Pareto optimized solution on $X$.

Definition 3, collection of Pareto optimized solutions: in terms of the given multi-objective optimization, the collection of Pareto optimized solutions can be defined as

$$\rho f^{*}=\left\{\overrightarrow{v}=\overrightarrow{f}\left(\overrightarrow{x}\right)|\overrightarrow{x}\in {\rho }^{*}\right\}$$

(29)

Here, the collection of the optimal solutions in the objective space is called the Pareto Front. It is important that the controllable symbol $\prec$ above can be changed into $\succ$ in terms of minimization.

To achieve a better Pareto Front, various solutions with effective strategies are proposed (Burachik et al., 2014; Campigotto et al., 2014; Chen and Zou, 2014; Khorram et al., 2014) [40,41,42,43].

3.3. Multi-objective Firefly Algorithm

One way of multi-objective optimization is to group all objectives into one single objective and thus the single-objective algorithm is not greatly modified. For instance, Apostolopoulos and Vlachos (2011) [44] investigate multi-objective optimization by the firefly algorithm in this way. Another method is to produce Pareto Front by extending the firefly algorithm; by expanding the firefly algorithm, Yang (2013) [11] proposes the multi-objective algorithm.

The process of optimization of the multi-objective algorithm is as follows. Firstly, the objective functions are defined. Secondly, the firefly groups are initialized and fireflies are distributed equally in the search space, which can be achieved by sampling. The tolerable errors and maximal iteration are defined, and then the brightness of fireflies can be evaluated and each firefly can be compared by the value of objective functions. If firefly $j$ dominates firefly $i$, firefly $i$ moves to firefly $j$, according to Formula (27). After moving, if $i$ cannot satisfy the constraint, a new firefly is produced. If one firefly is not dominated by any firefly, then the firefly is placed in the Pareto Front and a random vector is generated (the sum is one), and thus the optimal collective solution ${\rho }_{*}^t$ is achieved. Subsequently, the non-dominant solution collection is transmitted to the next iteration. After multiple iterations, the maximal iteration is achieved, and $n$ non-dominating solutions collections of the approximate Pareto can be generally obtained, and thus optimization is achieved.

For more effective random movement, the current optimized solution $g_{*}^t$ can be achieved by weighted sum of the minimal objective functions. Here,

$$\psi \left(x\right)={\sum }_{k=1}^k{w}_k{f}_k, {\sum }_{k=1}^k{w}_k=1$$

(30)

Here, $w_k=\raisebox{1ex}{$p_k$}\!\left/ \!\raisebox{-1ex}{$k$}\right.$, $p_k$ are random figures equally distributed from [0,1]. To ensure $\sum _k{w}_k=1$, $k$ figures equally distributed are generated and re-operated by zooming. In each iteration, $w_k$ should be chosen randomly, and hence the non-dominant solution can be sampled from the diversity of the Pareto optimized fronts.

From the perspective of the Pareto Front, if one firefly is not dominated by other fireflies, then the firefly is re-located as

$$x_i^{t+1}=g_{*}^t+\alpha \left(rand-1/2\right)$$

(31)

In the Formula, $g_{*}^t$ is the current optimized solution in terms of one group of a given random weight.

As discussed above, the procedure of the multi-objective firefly algorithm in weighing efficiency and robustness in disrupted supply chains is as follows.

Step 1: The objective functions are defined and the firefly groups $x_i \left(i=1, 2,\dots ,n\right)$ are initialized.

Step 2: Comparing $\rho f_i^{*}$ and $f_j^{*}$ $\left(i, j=1, 2, \dots , n, i\ne j\right)$, if $\rho f_j^{*}$ dominates $\rho f_i^{*}$, firefly $i$ moves to firefly $j$ according to Formula (27). If $\rho f_i^{*}=\rho f_j^{*}$ or both are equal, there is no non-dominant solutions and thus the random weight $w_k \left(k=1, 2, \dots , k\right)$ is produced. The optimized solution $g_{*}^t$ is found in all fireflies. $\psi$ is minimized in Formula (30), and the firefly moves randomly according to Formula (31).

Step 3: The Pareto Front is updated by non-dominant solutions and the number of optimized solutions is recorded. All the non-dominant solutions are transmitted to the next iteration and the brightness and locations of fireflies are updated.

Step 4: Step 2 is repeated till the maximal iteration is achieved, and then all the Pareto optimized solutions are obtained. The current optimized Pareto Front approximated has been found.

4. Performance Test of Multi-objective Firefly Algorithm

In order to test the validity of the multi-objective algorithm, MATLAB was used in a PC with Intel(R) Core(TM) i3-2310M, 2.1GHz, with a memory of 4G. The test functions in

Table 1. 5 testing functions.

Problems	Dimensions	Range	Objective functions
Zhang et al. (2004) [45]	1	[−5,7]	$f_1\left(x\right)=x^2$, $f_2\left(x\right)={\left(x-2\right)}^2$
E. Wang (2004) [46]	2	[−5,10]	$f_1\left(x,y\right)={\left(x^2+y^2\right)}^{\frac{1}{8}}$
			$f_2\left(x,y\right)={({\left(x-0.5\right)}^2+{\left(y-0.5\right)}^2)}^{\frac{1}{4}}$
FON (Deb et al., 2000) [47]	3	[−4,4]	$f_1\left(x\right)=1-\text{exp}(-\sum _{i=1}^3{(x_i-\frac{1}{\sqrt{3}})}^2)$
			$f_2\left(x\right)=1-\text{exp}(-\sum _{i=1}^3{(x_i+\frac{1}{\sqrt{3}})}^2)$
KUR (Deb et al., 2000) [47]	3	[−5,5]	$f_1\left(x\right)=\sum _{i=1}^{n-1}(-10\text{exp}(-0.2\sqrt{x_i^2+x_{i+1}^2}))$
			$f_2\left(x\right)=\sum _{i=1}^n({\left|x_i\right|}^{0.8}+5\sin x_i^3)$
ZDT3 (Deb et al., 2000) [47]	30	[0,1]	$f_1\left(x\right)=x_1$
			$f_2\left(x\right)=g\left(x\right)\left[1-\sqrt{x_1/g\left(x\right)}-\frac{x_1}{g\left(x\right)}\sin \left(10\pi x_1\right)\right]$
			$g\left(x\right)=1+9(\sum _{1=2}^n{x}_i)/\left(n-1\right)$

are considered.

The parameters of the multi-objective firefly algorithm are defined as follows: the number of fireflies is $n=500$, $\alpha =1.5$, ${\beta }_0=1$, $\gamma =1$, and the number of iterations is 200. The results of test functions are shown in Figure 2, Figure 4, Figure 6, Figure 8 and Figure 10, and the true value of the test functions is shown in Figure 3, Figure 5, Figure 9 and Figure 11.

From Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, the horizontal axis is $f_1$, and the vertical axis is $f_2$. The simulated results show that all the selected test functions are calculated by the multi-objective firefly algorithm. After the iteration, the Pareto Front is reached, satisfying the evaluation standards of the optimized solutions. Each test function achieves optimization by the multi-objective firefly algorithm.

5. Exemplifications

To illustrate the applicability of the multi-objective firefly algorithm to the three-echelon supply chain with multiple products and manufacturing centers in the mixed integer and linear planning model, a computer manufacturing company in China is examined. The company deals in a number of products with many manufacturing centers and distribution centers, which suits the three-echelon supply chain model with multiple products and manufacturing centers. As a result, the example can test the applicability of the multi-objective firefly algorithm. Two manufacturing centers of the company are located in two developed cities: Beijing and Shanghai. The company chooses 7 potential distribution centers in the most convenient and advanced provinces and the municipal cities in East China, South China, Central China, North China, Northwest, Southwest and Northeast China, namely, Beijing, Liaoning, Shanghai, Hubei, Guangdong, Sichuan and Shaanxi. There are 33 client zones in total, including 33 provincial regions like Hong Kong and Macau.

5.1. Elaboration

Two manufacturing centers can satisfy all requirements and produce two types of computers. One is the functional product: ordinary computers with marginal profit of 200 RMB. The other is an innovative product: a new-style computer with marginal profit of 1300 RMB. Fisher (1997) [48] separates products into two types in terms of the demand model, namely, the functional product and the innovative product. The demand for the functional product can be anticipated, with long product lifecycles but low marginal contribution, whereas the demand for the innovative product cannot be predicted, with short product lifecycles but high marginal contribution. For instance, salt, tissue and toothbrushes are examples of functional products. The launch of new cars and fashionable bags are instances of innovative products. The distribution centers can be located in 7 potential centers and the optimized distribution center relies on the model. The distances between the manufacturing center, the distribution center and the client zone can be recommended by means of the Baidu map. The focus of this study is node and link disruptions of supply chains and therefore demand in the client zone is definite and proportional to the population in the client zone. The population in the provincial regions is drawn from “China Statistics Yearbook in 2013” (China Statistic Yearbook, 2013; 2014) [49].

5.2. Demand in Client Area

The annual demand for the two products in 33 client zones is shown in

Table 2. Functional products in client area and demand for innovative products.

Client Area	Provincial Administrative Region	Demand for Ordinary Computers $N_1$ (Sets)	Demand for New Computers $N_1$ (Sets)
1	Beijing (BJ)	2069	207
2	Tianjin (TJ)	1413	141
3	Hebei (HEB)	7288	729
4	Shanxi (SAX)	3611	361
5	Neimenggu (NMG)	2490	249
6	Liaoning (LN)	4389	439
7	Jiling (JL)	2750	275
8	Heilongjiang (HLJ)	3834	383
9	Shanghai (SH)	2380	238
10	Jiangsu (JS)	7920	792
11	Zhejiang(ZJ)	5477	548
12	Anhui(AH)	5988	599
13	Fujian (FJ)	3748	375
14	Jiangxi (JX)	4504	450
15	Shandong (SD)	9685	969
16	Henan (HEN)	9406	941
17	Hubei (HUB)	5779	578
18	Hunan (HUN)	6639	664
19	Guangdong (GD)	10,594	1059
20	Guangxi (GX)	4682	468
21	Hainan (HAN)	887	89
22	Chongqing (CQ)	2945	295
23	Sichuan (SC)	8076	808
24	Guizhou (GZ)	3484	348
25	Yunnan (YN)	4659	466
26	Tibet (TB)	308	31
27	Shaanxi (SHX)	3753	375
28	Gansu (GS)	2578	258
29	Qinghai (QH)	573	57
30	Ningxia (NX)	647	65
31	Xinjiang (XJ)	2233	223
32	Hongkong (HK)	716	72
33	Macao (MAC)	57	6

. The annual demand for ordinary computers is 135,562, whereas the annual demand for the new-style computers is 13,558.

5.3. Distance among Manufacturing Center, Distribution Center and Client Areas

The distances between the manufacturing center, the distribution center and the client zone are calculated by the Baidu maps. The distance $G \left(a,b\right)$ from the manufacturing center $a$ to the distribution center $b$ is shown in

Table 3. Distance between manufacturing center and distribution center (kilometers).

	$b$	BJ	LN	SH	HUB	GD	SC	SHX
$a$
BJ		0	695	1261	1160	2117	1800	1088
SH		1257	1728	0	837	1475	1962	1375

and the distance $H \left(b,c\right)$ from the distribution center $b$ to the client zone $c$ is shown in

Table 4. Distance between distribution center and client areas (kilometers).

	$b$	BJ	LN	SH	HUB	GD	SC	SHX
$a$
BJ		0	695	1257	1162	2114	1809	1103
TJ		137	669	1086	1141	2094	1823	1145
HEB		315	1016	1138	969	1922	1501	794
SAX		521	1210	1356	949	1866	1317	610
NMG		485	1180	1730	1384	2301	1687	977
LN		695	0	1728	1819	2771	2473	1767
JL		1011	333	2034	2132	3085	2787	2080
HLJ		1251	573	2274	2372	3325	3027	2320
SH		1261	1720	0	843	1465	1965	1380
JS		1051	1526	303	539	1359	1662	1077
ZJ		1320	1785	178	721	1276	1856	1318
AH		1059	1622	465	390	1214	1513	928
FJ		1937	2402	773	910	935	2045	1656
JX		1431	2039	729	358	787	1493	1104
SD		446	1009	859	849	1817	1603	897
HEN		713	1369	943	516	1453	1193	485
HUB		1160	1816	837	0	981	1144	745
HUN		1487	2136	1086	354	651	1211	1007
GD		2117	2790	1475	984	0	1729	1636
GX		2340	2995	1902	1207	560	1215	1628
HAN		2759	3354	2043	1626	592	1700	2248
CQ		1764	2497	1685	871	1424	318	685
SC		1800	2530	1962	1146	1727	0	712
GZ		2148	2803	1843	1149	1096	663	1068
YN		2662	3315	2355	1661	1350	899	1569
TB		3636	4323	4196	3578	3606	2095	2833
SHX		1088	1821	1375	746	1636	712	0
GS		1485	2174	2010	1392	2282	857	646
QH		1686	2375	2247	1629	2519	1072	883
NX		1174	1843	1950	1461	2350	1429	736
XJ		3161	3850	3899	3270	4160	2792	2534
HK		2204	2936	1526	1110	180	1894	1799
MAC		2272	2916	1614	1139	139	1861	1792

.

5.4. Probability of Supply Chain Disruptions

Supply chain disruption may be caused by natural or human factors. Although it is difficult to quantify the probability of natural disasters, data recorded in history can help predict the probability. Li et al. (2013) [50] predict the probability of occurrence of natural disasters in the future on the basis of the data in the history. Here, it is assumed that most supply chain disruptions are caused by natural disasters in China while political situation is stable. The relative probability of occurrence of natural disasters at provincial level is calculated based on the loss caused by natural disasters in different regions in “Statistic Yearbook in China in 2013” (China Statistic Yearbook, 2013; 2014) [49]. The probability of disruption of the manufacturing center and the potential disruption is shown in

Table 5. Disruption probability of manufacturing center and distribution center.

Provincial Administrative Regions	BJ	LN	SH	HUB	GD	SC	SHX
Disruption probability	0.041	0.049	0.001	0.031	0.018	0.096	0.021

.

5.5. Relevant Costs

It is assumed that the cost of building each distribution center is 100,000,000 RMB. The lifespan of a distribution center is assumed to be 20 years and thus the fixed cost to be amortized is 500,000 RMB each year. When the distribution center is built, amortization cost assigned to the normal computers and new-type computers is 200 RMB and 800 RMB respectively. The production cost for each normal computer and new-type computer is 2500 RMB and 6000 RMB, respectively. The unit carrying cost for the normal computer and the new-type computer is 50 RBM and 100 RMB, respectively, in each distribution center. The unit transportation cost of each normal computer and new-type computer from the manufacturing center to the distribution center and also from the distribution center to the client zone is 40 RMB. The multi-objective firefly algorithm is adopted in the programming and operation of the corresponding exemplification in Matlab R2014b.

6. Results and Discussions

The supply chain disruption might involve node disruption, i.e., disruption at the distribution center, or disruption in the link between the manufacturing center and the distribution center. In what follows, these two circumstances are considered.

6.1. Disruption of Distribution Center

As the distance between distribution centers is large, it is assumed that the disruption in each distribution center is separate, and disruptions can occur simultaneously in many distribution centers. Each distribution center can have only one of the two possible statuses: normal or disrupted. Assume the distribution center will lose all its capacity in disruption. The probability of the disruption occurrence depends on the provincial administrative regions (Table 5). A disruption may occur in the distribution center, or in many centers simultaneously. Here, three simultaneous disruptions are considered, because the probability of disruption occurrence is small in four or more than four distribution centers. The model involves 1586 constraints, 739 variables and $63 \left(c_7^1+c_7^2+c_7^3\right)$ scenarios.

Before the model is operated, it is necessary to calculate the upper and lower boundary of $c\left(O\right)$ and $c\left(E\right)$. $c{\left(O\right)}_{min}$ is achieved by minimizing $c\left(O\right)$, and $c\left(E\right)$ is maximal at this point. $c{\left(E\right)}_{min}$ is directly achieved by minimizing $c\left(E\right)$, and the inventory selected is the most robust here. $c\left(O\right)$ is maximal through the inventory selected. In the disrupted distribution center, the boundary value between $c\left(O\right)$ and $c\left(E\right)$ is shown in

Table 6. Boundary value of $c\left(O\right)$ and $c\left(E\right)$ in disruption of distribution center.

$c{\left(O\right)}_{min}$	$c{\left(E\right)}_{max}$	$c{\left(E\right)}_{min}$	$c{\left(O\right)}_{max}$
7,383,873,820	10,410,177	218,552	15,235,138,900

.

The multi-objective firefly algorithm is simulated by 200 iterations and the most effective and robust supply chain networks with disrupted distribution center are shown in Figure 12 and Figure 13.

The total cost in the most efficient and the most robust supply chains and the minimal total cost with the optimal objective functions are shown in

Table 7. Comparing each cost in the disrupted distribution center.

Exemplification	Total Cost in the Most Efficient Supply Chains (RMB)	Total Cost in the Most Robust Supply Chains (RMB)	Total Cost of the Optimal Objective Functions (RMB)
Disrupted distribution center	7,394,283,997	15,235,357,452	7,386,193,817

. However, it can be seen from Table 7 that the total cost with the optimal objective function is lower than that in the most efficient supply chains because the operations cost increases with the optimal objective functions, while the reduction of the anticipated disruption cost is more than the increase in operations cost. In the increasingly competitive markets nowadays, the cost reduction may bring competitive edge to enterprises. When the optimal total cost is minimal with the objective function, the supply chain network can be as shown in Figure 14. Comparing Figure 12 and Figure 14, a portion of products is transmitted from the Beijing distribution center to Hubei distribution center. Hubei distribution center delivers the goods due to the fact that the disruption probability of Hubei distribution center is lower than that of the Beijing distribution center.

The multi-objective firefly algorithm is used and after 200 iterations, the approximate Pareto front is achieved (Figure 15), with horizontal axis as efficiency of supply chain and vertical axis as robustness of the supply chain. From Figure 15, it can be seen that when nodes or links are disrupted, the robustness decreases as the efficiency of the supply chain increases. It is a trading-off process. Meanwhile, Table 7 shows that the best solution makes the target function reach maximum. By multiple-objective algorithm, the Paredo front achieved is continuous, and thus efficiency and robustness can be derived at any point. Enterprises can choose the appropriate efficiency and robustness of supply chain in terms of specific features of their respective businesses and the prevailing status. For instance, businesses frequently affected by natural disasters may increase robustness of supply chain whereas businesses which are not sensitive to natural disasters can attach greater importance to efficiency.

6.2. Linkage Disruption between Manufacturing Center and Distribution Center

The link involves the two nodes and thus it is assumed that the probability of link disruption is the mean of the probability of disruption in all the linking nodes, and the probability of disruption in the link between the manufacturing center and the distribution center is shown in

Table 8. Link disruption probability between manufacturing center and distribution center.

Link Disruption Probability	BJ	LN	SH	HUB	GD	SC	SHX
Beijing	0.041	0.045	0.021	0.036	0.0295	0.0685	0.031
Shanghai	0.021	0.025	0.001	0.016	0.0095	0.0485	0.011

. A disruption may occur between the manufacturing center and the distribution center, or many simultaneous disruptions may occur, because the probability of disruptions is small in three or more than three links. To avoid complexity, at the most two simultaneous disruptions are considered. The model involves 1586 constraints, 739 variables and $105 \left(c_{14}^1+c_{14}^2\right)$ scenarios.

Before the operation of the model, the upper and lower boundary of $c\left(O\right)$ and $c\left(E\right)$ must be calculated. $c{\left(O\right)}_{min}$ can be achieved by minimizing $c\left(O\right)$ directly and $c\left(E\right)$ is maximal at this point. $c{\left(E\right)}_{min}$ can be achieved by minimizing $c\left(E\right)$ directly and $c\left(O\right)$ is maximal now. In the case of disruption in the link between the manufacturing center and the distribution center, boundary values of $c\left(O\right)$ and $c\left(E\right)$ are shown in

Table 9. Boundary values of $c\left(O\right)$ and $c\left(E\right)$ in Link disruption between manufacturing center and distribution center.

$c{\left(O\right)}_{min}$	$c{\left(E\right)}_{max}$	$c{\left(E\right)}_{min}$	$c{\left(O\right)}_{max}$
7,383,873,820	15,091,028	1,147,609	1,523,5138,900

.

The multi-objective firefly algorithm is simulated and after 200 iterations, the most effective supply chain network with disruption at the link between the manufacturing center and the distribution center and the supply chain network with the minimal optimized total cost by the objective functions are shown in Figure 16 and Figure 17. Comparing Figure 16 and Figure 17, a portion of products is transmitted from Beijing distribution center to Shanghai and Shaanxi distribution centers. Shanghai and Shaanxi distribution centers deliver the goods, due to the fact that the disruption probability of Shanghai and Shaanxi distribution centers is lower than that in the Beijing distribution center.

The multi-objective firefly algorithm is used and after 200 iterations, the approximate Pareto front is achieved (Figure 18), with horizontal axis as efficiency and vertical axis as robustness of supply chain. From Figure 15 and Figure 18, it can be seen that when nodes or links are disrupted, the robustness decreases as efficiency increases. Based on Figure 18, enterprises can weigh efficiency and robustness of supply chain and determine the appropriate efficiency and robustness in terms of the nature of their industry and their own circumstances. For instance, those likely to be affected by natural disasters would increase their robustness to some extent; those insensitive to the natural disasters would attach great importance to efficiency.

7. Conclusions

The design of supply chain is extremely demanding, due to the internal complexity and uncertainty of external factors. In the past, the complexity centered on each component of the supply chain, whereas the attention has recently been shifted to the design, performance and analysis of the overall supply chain. This paper proposes a mathematical model of mixed integer linear planning with an aim to address the complexity of the above issue. Robustness is considered in the design of supply chains with multiple products and manufacturing centers, including node and link disruptions. The mathematical model of mixed integer linear planning concerns optimization of multiple objectives, focusing on weighing efficiency and robustness of supply chain. Efficiency and robustness can be in conflict and contradiction and optimization of one objective is at the cost of the other and hence it is difficult to achieve optimization. It is insufficient to solve the problem by the traditional mathematical means. For instance, it is required that each objective function and constraint should be differentiable with low efficiency and sophisticated calculation. This paper, for the first time, introduces the multi-objective firefly algorithm based on Pareto optimization to address the problem of supply chain disruption. The simulation is used to test its validity. The Pareto lines are provided to weigh efficiency and robustness in disrupted supply chains and partial optimization can also be avoided, which may extend the research and application of the firefly algorithm deployed to solve the weighing problem. The result is the continuous Pareto curve instead of a number of independent points, and the robustness and cost of supply chains corresponding to their efficiency can be derived.

From the exemplification of the mixed integer linear planning of the three-echelon supply chain with multiple products and manufacturing centers, it can be seen that the robustness of supply chain reduces as efficiency increases with node disruption or link disruption, whereas enterprises expect to increase their efficiency and also improve their robustness. These are inversely proportionate and therefore it is important to weigh efficiency and robustness in disrupted supply chains to achieve the optimal objective functions. The simulation shows that the total cost with the optimal objective functions is lower than that with the most efficient supply chain. In other words, the reduction of anticipated disruption cost caused by the increased robustness is more than the increase of the operations cost caused by decreasing efficiency of supply chains. As a result, the Pareto front approximated, when the objective function is optimal, provides the evidence for enterprises to opt for appropriate efficiency and robustness of supply chains. Overall, this approach provides a tool for quantifying the decision of planning supply chains. In addition, it increases robustness of supply chains, and reduces the total cost of supply chains.

The three-echelon supply chain model involves multiple products and manufacturing centers, and other supply chain models are analyzed on this basis. To avoid complexity, the mathematical model of mixed integer linear planning is adopted, and the non-linear mathematical model can also be considered in the future. In the model, it is assumed that once the distribution center is built, it has to serve the client zone. Each distribution center and client zone can only be supplied by one manufacturing center and distribution center, respectively. There is no inventory accumulation and loss, and the demand of each client zone can be satisfied. In the future studies, these assumptions can be relaxed.