# Solving a Two-Stage Stochastic Capacitated Location-Allocation Problem with an Improved PSO in Emergency Logistics

Authors: Ye Deng; Wanhong Zhu; Jian Tang; Jianfei Qin
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

A stochastic expected value model and its deterministic conversion are developed to formulate a two-stage stochastic capacitated location-allocation (LA) problem in emergency logistics; that is, the number and capacities of supply centers are both decision variables. To solve these models, an improved particle swarm optimization algorithm with the Gaussian cloud operator, the Restart strategy, and the adaptive parameter strategy is developed. The algorithm is integrated with the interior point method to solve the second-stage model. The numerical example proves the effectiveness and efficiency of the conversion method for the stochastic model and the proposed strategies that improve the algorithm.


---

## Body

## 1. Introduction

A location-allocation (LA) problem, also known as facility location problem (FLP), involves locating a number of facilities to which customers are allocated to minimize the cost of satisfying customer demands. It is an important problem in supply chain or logistics management and greatly affects long-term transportation and storage decisions. Many enterprises and government departments focus on this problem to reduce cost and to improve efficiency, especially in emergency logistics given the frequent occurrence of disasters, epidemic, security incidents, and other emergencies nowadays. This problem has complex and uncertain features, such as the changing demands, allocations, and locations of customers or facilities.Many studies have been conducted on the basic LA problem since the proposal of Cooper [1] in 1963. However, most of the studies were for deterministic cases and not for uncertain occurrences. In the recent two decades, several models for uncertain occurrences have been proposed and solved by different algorithms. Logendran and Terrell [2] first considered a stochastic uncapacitated LA problem and proposed an expected value model (EVM) to maximize the net profits. Carrizosa et al. [3, 4] proposed a LA problem that considers the locations of both customers and facilities, which may be regions that have several probability distributions. Liu [5, 6] contributed to the uncertainty theory by proposing three stochastic models [7] and three fuzzy programing models [8] for the capacitated LA problem. A hybrid intelligent algorithm that consists of a network simplex algorithm, a simulation, and a genetic algorithm was developed to solve the stochastic and fuzzy models above. Silva and De La Figuera [9] studied the capacitated facility location problem with constrained backlogging probabilities and solved it using a heuristic method based on a reactive greedy adaptive search procedure. Wang and Shi-Wei [10] proposed a robust optimization model for a logistics center LA problem and compared it with stochastic and deterministic optimization models. Two algorithms, namely, an enumeration method and a genetic algorithm, were adopted to solve the problem. Yao et al. [11] considered a joint facility LA and inventory problem with stochastic demands. The problem involves identifying the best locations of warehouses and the inventory levels and allocating customers. A heuristic integrated approximation and transformation technique was developed to solve the problem. Wen and Kang [12] considered a facility LA problem with random fuzzy demands. They proposed a hybrid intelligent algorithm, similar to the method of Zhou and Liu [8], which consists of the simplex algorithm, a random fuzzy simulation, and a genetic algorithm. A similar method was also adopted by Mousavi and Niaki [13] in solving a LA problem with a fuzzy variable and customer location demands, which were normal random variables. Vidyarthi and Jayaswal [14] proposed a nonlinear integer programing model to solve a LA problem with immobile servers, stochastic demands, and congestions. Pereira et al. [15] presented a probabilistic maximal covering LA problem and proposed a hybrid algorithm to solve it. They formulated a linear programing model to efficiently solve small and medium problems and a flexible adaptive large neighborhood search heuristic to solve large problems. Alizadeh et al. [16] considered a capacitated multifacility LA problem with stochastic demands. The capacitated subsources of each facility could be utilized when the number of demand points, that is, the planned total requirements, are exceeded. Alizadeh et al. [17] transformed the mixed-integer nonlinear programing model to a simple formulation model and proposed a genetic algorithm (GA) and a colonial competitive algorithm (CCA) to solve medium and large problems.In this study, we consider a two-stage stochastic capacitated LA problem (SCLAP) in the context of emergency logistics. In the management of emergency logistics, the core problem is utilizing the limited relief supplies rapidly and efficiently. Hence, predesigned supply centers, that is, emergency logistics distribution centers, are important. The “appropriate” number, size, and location of supply centers have become a comprehensive decision problem that should not be addressed separately. The demands of different customers are mostly uncertain and depend largely on different scenarios. Hence, we focus on the uncertainty in demand quantity, which is assumed as the only independent stochastic variable in this paper that follows a given regular stochastic distribution. To solve SCLAP two conditions must be satisfied: the constraint of stochastic quantity demands and the minimization of the generalized cost. The generalized cost consists of two parts, namely, the sum of the costs of building and maintaining supply centers and the stochastic costs of transportation to each customer from each supply center. Therefore, we determine the appropriate number, capacities, and locations of supply centers in this paper.In a traditional capacitated LA problem (CLAP), each customer can be supplied by existing supply centers and can be supplied by more than one center at the same time. Hence, the problem becomes NP-hard and difficult to solve [7]. Moreover, in the two-stage SCLAP, the number, the capacities, and the locations of supply centers are considered decision variables. The research on this model is extremely weak, and few relevant papers have been found [18–23]. Furthermore, most of the research considered the number and the capacities of supply centers separately, and no research that considered both variables together has been found.To efficiently solve this model, an improved particle swarm optimization (PSO) algorithm is proposed. The algorithm consists of three improvement strategies: the Gaussian cloud operator, the Restart strategy, and the adaptive parameter strategy. The second stage of the problem is modeled as a linear program. Hence, we adopt the interior point method of the time-consuming simplex method. We convert the initial stochastic programing model to a crisp model, thus reducing the computing time dramatically based on the assumption of the demands with independent regular distributions and the uncertainty theory proposed by Liu [6].The remainder of this paper is organized as follows. Section2 describes the random EVM and the crisp model for the two-stage SCLAP. Section 3 presents the details of the hybrid algorithm solution to the model. Section 4 introduces a case study of the new model and verifies the algorithm efficiency with the improvement strategies. Section 5 concludes with the contributions and innovations of this paper and presents the future research directions.


## 2. Model Formulation

### 2.1. Problem Description and Theoretical Foundation

To model the two-stage SCLAP, the following assumptions should be considered: the graph of all the nodes is a complete graph, each customer node can be connected with all supply nodes but cannot be connected with another customer node, the weight of each edge between two nodes is measured by the Euclidian distance plus the transportation volume, the locations of customer nodes are fixed and the demand quantities are stochastic, and the capacity constraint is only imposed on supply nodes. The notation and variables used in the following formulations are defined in Descriptions of Notations and Variables.To model the two-stage SCLAP, we first apply the EVM introduced by Zhou and Liu [7] to the SCLAP. Then, we extend the classic one-stage EVM to a two-stage model and provide a deterministic equivalent form. We introduce several basic definitions, theorems of probability, and uncertainty theories.Definition 1 (see [5]).
LetΩ be a nonempty set and A the σ-algebra of the subsets (called events) of Ω. The set function Pr is called a probability measure if it satisfies the following conditions.Axiom 1 (normality).
(1)
Pr
⁡
Ω
=
1
.Axiom 2 (nonnegativity).
(2)
Pr
⁡
A
≥
0
f
o
r
a
n
y
A
∈
A
.Axiom 3 (countable additivity).
For every countable sequence of mutually disjoint events{Ai}, we obtain(3)Pr⁡⋃i=1∞Ai=∑i=1∞Pr⁡Ai.Definition 2 (see [5]).
LetΩ be a nonempty set, A the σ-algebra of the subsets of Ω, and Pr the probability measure. Then, the triplet (Ω,A,Pr) is called a probability space.Definition 3 (see [5]).
A random variable is a measurable function from the probability space(Ω,A,Pr) to the set of real numbers; that is, for any Borel set B of real numbers, the set is(4)ξ∈B=ω∈Ω∣ξω∈B.Liu [6] also proposed the definitions formeasure inversion theorem,regular uncertainty distributions, andinverse uncertainty distribution as supplements to the uncertainty theory. The uncertainty distributions are specifically presented as stochastic distributions, and the uncertainty measure M can also be replaced by the probability measure Pr. He proposed several theorems that can help in the conversion of the stochastic model into a crisp model based on these definitions. The related definitions and theorems are listed as follows.Definition 4 (see [6] measure inversion theorem).
Letξ be an uncertain variable with an uncertainty distribution Φ. Then, for any real number x, we have(5)Mξ≤x=Φx,Mξ>x=1-Φx.Definition 5 (see [6] regular uncertainty distribution).
An uncertainty distributionΦ(x) is regular if it is a continuous and strictly increasing function with respect to x at which 0<Φ(x)<1 and (6)limx→-∞⁡Φx=0,limx→+∞⁡Φx=1.Definition 6 (see [6] inverse uncertainty distribution).
Letξ be an uncertain variable with a regular uncertainty distribution Φ(x). Then, the inverse function Φ-1(x) is the inverse uncertainty distribution of ξ.Theorem 7 (see [6]).
Letξ1,ξ2,…,ξn be independent uncertain variables with regular uncertainty distributions Φ1,Φ2,…,Φn, respectively. If f is a strictly increasing function, then ξ=f(ξ1,ξ2,…,ξn) has an inverse uncertainty distribution(7)Ψ-1α=fΦ1-1α,Φ2-1α,…,Φn-1α.Theorem 8 (see [6]).
Letξ1 and ξ2 be independent normal uncertain variables N(e1,σ1) and N(e2,σ2), respectively. Then, the sum of ξ1+ξ2 is also a normal uncertain variable N(e1+e2,σ1+σ2); that is,(8)Ne1,σ1+Ne2,σ2=Ne1+e2,σ1+σ2.


### 2.2. Initial Expected Value Model

The mathematical formulation, which is an expected value model (EVM), proposed by Zhou and Liu for the initial SCLAP can be defined as follows [7]:(9)Fs=minx,y⁡Eminz∈Zω⁡∑i=1n∑j=1mzijxi-aj2+yi-bj2,which is subject to(10)∑i=1nzij=ξjω,∀j∈M,ω∈Ω;(11)∑j=1mzij≤si,∀i∈N;(12)zij≥0,∀i∈N,j∈M;(13)xi,yi∈R2,∀i∈N;(14)aj,bj∈R2,∀j∈M,where E(X(ω)) denotes the expected value of uncertain variable X. If (15)Cx,y∣ω=minz∈Zω⁡∑i=1n∑j=1mzijxi-aj2+yi-bj2,(9) will be reformulated as follows:(16)Fs=minx,y⁡∫0∞Prω∈Ω∣Cx,y∣ω≥rdr.In this mathematical model, the decision variable is the location of the supply centers(x,y), and the objective function is the minimum value of the sum of the demand and distance from each supply center to the customer in (9) or (16). The first constraint (see (10)) states that a customer’s stochastic demand must be satisfied. The second constraint (see (11)) states that the capacity of each supply center must be sufficient to supply all customers with demands. Equations (12), (13), and (14) are used to keep the variables within a certain range.According to Zhou and Liu [7], the EVM can be easily solved by a hybrid intelligent algorithm that consists of the network simplex algorithm, stochastic simulations, and the GA. However, in emergency logistics management, in addition to making a decision on supply center locations, the number and capacities of supply centers should also be considered. Therefore, the SCLAP should be extended into a two-stage model. Variables n and si are both fixed in the initial EVM, whereas, in the two-stage EVM, these variables are decision variables.


### 2.3. Two-Stage Expected Value Model

In the two-stage EVM for SCLAP, the decision variables extend to the numberp, the capacities S, and the locations of supply centers (x,y). The first stage aims to determine the minimum value of the generalized cost and generates the number and capacity values for the second stage. The second stage utilizes the number and capacity values from the first stage and determines the minimum value of the stochastic demand cost of transportation. In addition, to provide a comprehensive description of the optimization objective in emergency logistics, we introduce a nonlinear function based on the objective function in the initial model to calculate the generalized cost. The two-stage EVM for SCLAP can be formulated as follows.First stage(17)minp,S,x,y⁡E∑i=1pA+δ·si+Fsp,Sis subject to(18)Smin≤si≤Smax,∀i∈N;(19)∑i=1psi≥E∑j=1mξjω,∀p∈N.Second stage(20)Fsp,S=minz∈Zω∑i=1p∑j=1mzijxi-aj2+yi-bj2is subject to(21)∑i=1pzij=ξjω,∀j∈M,ω∈Ω;(22)∑j=1mzij≤si,∀i∈N,si∈S;(23)zij≥0,∀i∈N,j∈M;(24)xi,yi∈R2,∀i∈N;(25)aj,bj∈R2,∀j∈M,where A (A>0) denotes the fixed cost of constructing a supply center and δ(0≤δ≤1) denotes the variable cost coefficient to maintaining an operating supply center. Equation (17) represents the expected generalized cost. The first portion, ∑i=1pA+δ·si, is the sum of the cost of construction and the cost of maintaining the supply centers, and the latter portion, Fsp,S, is the sum of the costs of the demands for transportation from each supply center to the customer, similar to the initial EVM. The first constraint in (18) ensures that the capacity of each supply center is limited within a reasonable range from Smin to Smax. Equation (19) stipulates that the total capacity of supply centers must exceed the total expected demands of customers regardless of the value of N or the number of supply centers p. The formulations in the second stage denote the same items as those in the initial model, except that fixed constants n and si become variables.


### 2.4. Deterministic Equivalent Models

Evidently, the two-stage EVM model is much more complicated than the initial model, and the stochastic simulation method is extraordinarily time-consuming because the number and the capacities of supply centers are both unknown. However, on the basis of the definitions and theorems in Section2.1, if we assume that the demand variables are nonnegative normal distributions, that is, ∀j∈M,ξj~N(ej,σj), we can convert the two-stage EVM of the SCLAP into a deterministic two-stage model. The following equivalent formula of the first stage model can be expressed as (26)minp,S,x,y⁡pA+δ∑i=1psi+EFsp,Sand subject to(27)Smin≤si≤Smax,∀i∈N,∑i=1psi≥∑j=1mej,∀p∈N.Proof.
Customer demandξ follows the normalized distribution (Theorem 8), which belongs to theregular uncertainty distribution, and f=∑j=1mξj(ω) is a strictly increasing function. Hence, one has the following.
(1) Objective conversion is as follows:(28)E∑i=1pA+δ·si+Fsp,S=E∑i=1pA+δ·si+EFsp,S=pA+δ∑i=1psi+EFsp,S.
(2) Constraint conversion is as follows:∵
ξj~N(ej,σj) and according to Theorem 8∴
∑j=1mξj(ω)~N∑j=1mej,∑j=1mσj∴
∑i=1psi≥E∑j=1mξj(ω)⇒∑i=1psi≥∑j=1mEξj(ω)=∑j=1mejThe following equivalent formula of the second-stage model can be established as (29)EFsp,S=minz∈EZω∑i=1p∑j=1mzijxi-aj2+yi-bj2and subject to(30)∑i=1pzij=ej,∀j∈M,ω∈Ω;∑j=1mzij≤si,∀i∈N,si∈S;zij≥0,∀i∈N,j∈M;xi,yi∈R2,∀i∈N;aj,bj∈R2,∀j∈M.Proof.
Customer demandξj=∑i=1pzij follows the normalized distribution (Theorem 8), which belongs to theregular uncertainty distribution. Hence, one has the following.
(1) Objective conversion is as follows:(31)Eminz∈Zω⁡∑i=1p∑j=1mzijxi-aj2+yi-bj2=minz∈Zω⁡∑j=1mE∑i=1pzijxi-aj2+yi-bj2=minz∈Zω⁡∑j=1m∑i=1pEzijxi-aj2+yi-bj2=minz∈EZω⁡∑i=1p∑j=1mzijxi-aj2+yi-bj2.
(2) Constraint conversion is as follows:(32)ξj=∑i=1pzij⟹Eξj=E∑i=1pzij=∑i=1pEzij=∑i=1pzij′=ej,zij′∈EZω.
After the conversion of the two-stage EVM, we can obtain a two-stage deterministic model, which can be solved much easily. Similarly, we can also convert the chance-constrained programing (CCP) model into a deterministic model. Assuming that the supply centers satisfy the demand of customerj with a probability βj, the total supply capacity satisfies the total demand of customers with probability α, and the demand variable ξj(∀j∈M) complies with the same kind of regular stochastic distribution Φj, the CCP of the two-stage SCLAP can be formulated as follows (only the values different from the EVM are listed for simplicity).
Equation (19) is reformulated as(33)Pr⁡∑i=1psi≥∑j=1mξjω≥α,∀p∈N;and (21) is reformulated as (34)Pr⁡∑i=1pzij≥ξjω≥βj,∀j∈M,ω∈ΩAccording to Theorem 7, the stochastic function f=∑j=1mξj(ω) is strictly increasing. Therefore, one has the following.
Equation (33) can be converted into(35)∑i=1psi≥∑j=1mΦj-1α,∀p∈N,and (34) can be converted into(36)∑i=1pzij≥Φj-1βj,∀j∈M.However, if the SCLAP is a dependent-chance programing (DCP) model (only the values different from the EVM are listed for simplicity), (17) is reformulated as(37)minp,S,x,y⁡Pr⁡∑i=1pA+δ·si+Fsp,S≤Cand (19) is reformulated as(38)∑i=1psi≥∑j=1mξjω,∀p∈N.
The objective is to maximize the probability, which cannot be converted into a deterministic model. Hence, the stochastic simulation process must be activated to calculate the probability. In this paper, we only focus on the EVM and develop an improved PSO (IPSO) algorithm to solve this problem given that the two-stage SCLAP is newly proposed.


## 2.1. Problem Description and Theoretical Foundation

To model the two-stage SCLAP, the following assumptions should be considered: the graph of all the nodes is a complete graph, each customer node can be connected with all supply nodes but cannot be connected with another customer node, the weight of each edge between two nodes is measured by the Euclidian distance plus the transportation volume, the locations of customer nodes are fixed and the demand quantities are stochastic, and the capacity constraint is only imposed on supply nodes. The notation and variables used in the following formulations are defined in Descriptions of Notations and Variables.To model the two-stage SCLAP, we first apply the EVM introduced by Zhou and Liu [7] to the SCLAP. Then, we extend the classic one-stage EVM to a two-stage model and provide a deterministic equivalent form. We introduce several basic definitions, theorems of probability, and uncertainty theories.Definition 1 (see [5]).
LetΩ be a nonempty set and A the σ-algebra of the subsets (called events) of Ω. The set function Pr is called a probability measure if it satisfies the following conditions.Axiom 1 (normality).
(1)
Pr
⁡
Ω
=
1
.Axiom 2 (nonnegativity).
(2)
Pr
⁡
A
≥
0
f
o
r
a
n
y
A
∈
A
.Axiom 3 (countable additivity).
For every countable sequence of mutually disjoint events{Ai}, we obtain(3)Pr⁡⋃i=1∞Ai=∑i=1∞Pr⁡Ai.Definition 2 (see [5]).
LetΩ be a nonempty set, A the σ-algebra of the subsets of Ω, and Pr the probability measure. Then, the triplet (Ω,A,Pr) is called a probability space.Definition 3 (see [5]).
A random variable is a measurable function from the probability space(Ω,A,Pr) to the set of real numbers; that is, for any Borel set B of real numbers, the set is(4)ξ∈B=ω∈Ω∣ξω∈B.Liu [6] also proposed the definitions formeasure inversion theorem,regular uncertainty distributions, andinverse uncertainty distribution as supplements to the uncertainty theory. The uncertainty distributions are specifically presented as stochastic distributions, and the uncertainty measure M can also be replaced by the probability measure Pr. He proposed several theorems that can help in the conversion of the stochastic model into a crisp model based on these definitions. The related definitions and theorems are listed as follows.Definition 4 (see [6] measure inversion theorem).
Letξ be an uncertain variable with an uncertainty distribution Φ. Then, for any real number x, we have(5)Mξ≤x=Φx,Mξ>x=1-Φx.Definition 5 (see [6] regular uncertainty distribution).
An uncertainty distributionΦ(x) is regular if it is a continuous and strictly increasing function with respect to x at which 0<Φ(x)<1 and (6)limx→-∞⁡Φx=0,limx→+∞⁡Φx=1.Definition 6 (see [6] inverse uncertainty distribution).
Letξ be an uncertain variable with a regular uncertainty distribution Φ(x). Then, the inverse function Φ-1(x) is the inverse uncertainty distribution of ξ.Theorem 7 (see [6]).
Letξ1,ξ2,…,ξn be independent uncertain variables with regular uncertainty distributions Φ1,Φ2,…,Φn, respectively. If f is a strictly increasing function, then ξ=f(ξ1,ξ2,…,ξn) has an inverse uncertainty distribution(7)Ψ-1α=fΦ1-1α,Φ2-1α,…,Φn-1α.Theorem 8 (see [6]).
Letξ1 and ξ2 be independent normal uncertain variables N(e1,σ1) and N(e2,σ2), respectively. Then, the sum of ξ1+ξ2 is also a normal uncertain variable N(e1+e2,σ1+σ2); that is,(8)Ne1,σ1+Ne2,σ2=Ne1+e2,σ1+σ2.


## 2.2. Initial Expected Value Model

The mathematical formulation, which is an expected value model (EVM), proposed by Zhou and Liu for the initial SCLAP can be defined as follows [7]:(9)Fs=minx,y⁡Eminz∈Zω⁡∑i=1n∑j=1mzijxi-aj2+yi-bj2,which is subject to(10)∑i=1nzij=ξjω,∀j∈M,ω∈Ω;(11)∑j=1mzij≤si,∀i∈N;(12)zij≥0,∀i∈N,j∈M;(13)xi,yi∈R2,∀i∈N;(14)aj,bj∈R2,∀j∈M,where E(X(ω)) denotes the expected value of uncertain variable X. If (15)Cx,y∣ω=minz∈Zω⁡∑i=1n∑j=1mzijxi-aj2+yi-bj2,(9) will be reformulated as follows:(16)Fs=minx,y⁡∫0∞Prω∈Ω∣Cx,y∣ω≥rdr.In this mathematical model, the decision variable is the location of the supply centers(x,y), and the objective function is the minimum value of the sum of the demand and distance from each supply center to the customer in (9) or (16). The first constraint (see (10)) states that a customer’s stochastic demand must be satisfied. The second constraint (see (11)) states that the capacity of each supply center must be sufficient to supply all customers with demands. Equations (12), (13), and (14) are used to keep the variables within a certain range.According to Zhou and Liu [7], the EVM can be easily solved by a hybrid intelligent algorithm that consists of the network simplex algorithm, stochastic simulations, and the GA. However, in emergency logistics management, in addition to making a decision on supply center locations, the number and capacities of supply centers should also be considered. Therefore, the SCLAP should be extended into a two-stage model. Variables n and si are both fixed in the initial EVM, whereas, in the two-stage EVM, these variables are decision variables.


## 2.3. Two-Stage Expected Value Model

In the two-stage EVM for SCLAP, the decision variables extend to the numberp, the capacities S, and the locations of supply centers (x,y). The first stage aims to determine the minimum value of the generalized cost and generates the number and capacity values for the second stage. The second stage utilizes the number and capacity values from the first stage and determines the minimum value of the stochastic demand cost of transportation. In addition, to provide a comprehensive description of the optimization objective in emergency logistics, we introduce a nonlinear function based on the objective function in the initial model to calculate the generalized cost. The two-stage EVM for SCLAP can be formulated as follows.First stage(17)minp,S,x,y⁡E∑i=1pA+δ·si+Fsp,Sis subject to(18)Smin≤si≤Smax,∀i∈N;(19)∑i=1psi≥E∑j=1mξjω,∀p∈N.Second stage(20)Fsp,S=minz∈Zω∑i=1p∑j=1mzijxi-aj2+yi-bj2is subject to(21)∑i=1pzij=ξjω,∀j∈M,ω∈Ω;(22)∑j=1mzij≤si,∀i∈N,si∈S;(23)zij≥0,∀i∈N,j∈M;(24)xi,yi∈R2,∀i∈N;(25)aj,bj∈R2,∀j∈M,where A (A>0) denotes the fixed cost of constructing a supply center and δ(0≤δ≤1) denotes the variable cost coefficient to maintaining an operating supply center. Equation (17) represents the expected generalized cost. The first portion, ∑i=1pA+δ·si, is the sum of the cost of construction and the cost of maintaining the supply centers, and the latter portion, Fsp,S, is the sum of the costs of the demands for transportation from each supply center to the customer, similar to the initial EVM. The first constraint in (18) ensures that the capacity of each supply center is limited within a reasonable range from Smin to Smax. Equation (19) stipulates that the total capacity of supply centers must exceed the total expected demands of customers regardless of the value of N or the number of supply centers p. The formulations in the second stage denote the same items as those in the initial model, except that fixed constants n and si become variables.


## 2.4. Deterministic Equivalent Models

Evidently, the two-stage EVM model is much more complicated than the initial model, and the stochastic simulation method is extraordinarily time-consuming because the number and the capacities of supply centers are both unknown. However, on the basis of the definitions and theorems in Section2.1, if we assume that the demand variables are nonnegative normal distributions, that is, ∀j∈M,ξj~N(ej,σj), we can convert the two-stage EVM of the SCLAP into a deterministic two-stage model. The following equivalent formula of the first stage model can be expressed as (26)minp,S,x,y⁡pA+δ∑i=1psi+EFsp,Sand subject to(27)Smin≤si≤Smax,∀i∈N,∑i=1psi≥∑j=1mej,∀p∈N.Proof.
Customer demandξ follows the normalized distribution (Theorem 8), which belongs to theregular uncertainty distribution, and f=∑j=1mξj(ω) is a strictly increasing function. Hence, one has the following.
(1) Objective conversion is as follows:(28)E∑i=1pA+δ·si+Fsp,S=E∑i=1pA+δ·si+EFsp,S=pA+δ∑i=1psi+EFsp,S.
(2) Constraint conversion is as follows:∵
ξj~N(ej,σj) and according to Theorem 8∴
∑j=1mξj(ω)~N∑j=1mej,∑j=1mσj∴
∑i=1psi≥E∑j=1mξj(ω)⇒∑i=1psi≥∑j=1mEξj(ω)=∑j=1mejThe following equivalent formula of the second-stage model can be established as (29)EFsp,S=minz∈EZω∑i=1p∑j=1mzijxi-aj2+yi-bj2and subject to(30)∑i=1pzij=ej,∀j∈M,ω∈Ω;∑j=1mzij≤si,∀i∈N,si∈S;zij≥0,∀i∈N,j∈M;xi,yi∈R2,∀i∈N;aj,bj∈R2,∀j∈M.Proof.
Customer demandξj=∑i=1pzij follows the normalized distribution (Theorem 8), which belongs to theregular uncertainty distribution. Hence, one has the following.
(1) Objective conversion is as follows:(31)Eminz∈Zω⁡∑i=1p∑j=1mzijxi-aj2+yi-bj2=minz∈Zω⁡∑j=1mE∑i=1pzijxi-aj2+yi-bj2=minz∈Zω⁡∑j=1m∑i=1pEzijxi-aj2+yi-bj2=minz∈EZω⁡∑i=1p∑j=1mzijxi-aj2+yi-bj2.
(2) Constraint conversion is as follows:(32)ξj=∑i=1pzij⟹Eξj=E∑i=1pzij=∑i=1pEzij=∑i=1pzij′=ej,zij′∈EZω.
After the conversion of the two-stage EVM, we can obtain a two-stage deterministic model, which can be solved much easily. Similarly, we can also convert the chance-constrained programing (CCP) model into a deterministic model. Assuming that the supply centers satisfy the demand of customerj with a probability βj, the total supply capacity satisfies the total demand of customers with probability α, and the demand variable ξj(∀j∈M) complies with the same kind of regular stochastic distribution Φj, the CCP of the two-stage SCLAP can be formulated as follows (only the values different from the EVM are listed for simplicity).
Equation (19) is reformulated as(33)Pr⁡∑i=1psi≥∑j=1mξjω≥α,∀p∈N;and (21) is reformulated as (34)Pr⁡∑i=1pzij≥ξjω≥βj,∀j∈M,ω∈ΩAccording to Theorem 7, the stochastic function f=∑j=1mξj(ω) is strictly increasing. Therefore, one has the following.
Equation (33) can be converted into(35)∑i=1psi≥∑j=1mΦj-1α,∀p∈N,and (34) can be converted into(36)∑i=1pzij≥Φj-1βj,∀j∈M.However, if the SCLAP is a dependent-chance programing (DCP) model (only the values different from the EVM are listed for simplicity), (17) is reformulated as(37)minp,S,x,y⁡Pr⁡∑i=1pA+δ·si+Fsp,S≤Cand (19) is reformulated as(38)∑i=1psi≥∑j=1mξjω,∀p∈N.
The objective is to maximize the probability, which cannot be converted into a deterministic model. Hence, the stochastic simulation process must be activated to calculate the probability. In this paper, we only focus on the EVM and develop an improved PSO (IPSO) algorithm to solve this problem given that the two-stage SCLAP is newly proposed.


## 3. Improved IPSO

The LAP has been proven to be NP-hard [24]. Hence, it is a two-stage SCLAP. Heuristic methods are the best at handling large NP-hard problems, especially under uncertain environments. We develop an improved IPSO combined with the interior point method [25] to solve the deterministic equivalent conversion of the two-stage EVM.


### 3.1. Design of the IPSO

In 1995, Kennedy and Eberhart [26] first proposed a swarm intelligence optimization algorithm, which was inspired by the flocking of birds and schooling of fish against predation, called the particle swarm optimization (PSO). PSO has been successfully applied to parameter optimization [27], combinatorial optimization [28], pattern recognition [29, 30], data mining [31], and other fields [32, 33]. However, local optimal solution is often obtained when general PSO algorithm is applied to solve the continuous location of distribution centers. Therefore, scholars have proposed many improvement strategies for a much efficient and effective PSO [34–37]. Recently, Zhan et al. [38] proposed an improved PSO based on neighbor heuristic and Gaussian cloud learning, and the results proved its superiority over many PSO variants. We also adopt a Gaussian cloud operator combined with an adaptive parameter strategy and a Restart strategy to improve the general PSO algorithm of the two-stage SCLAP. For a detailed description of the IPSO algorithm, an algorithm flow chart is provided as a supplement in Figure 1.Figure 1
IPSO algorithm flow chart.In Figure1, the loop program is repeated for IterMax times, and three improvement strategies, which are shown in the curved boxes, are adopted. If no better solution is found during ΔT(ΔT=ΔIter), theRestart strategy will be activated; otherwise, theGaussian cloud operator will be activated. In addition, several IPSO parameters will adaptively change along with the iteration (theadaptive parameter strategy).Based on (27), the range of supply center number p can be determined by (39)p≥∑j=1mejSmax+1,where [x] denotes taking the integer portion of x. The capacity of each supply center must be within [Smin,Smax], which means that if a candidate solution of si exceeds Smax, a penalty will be given in the general cost function, and if si is less than Smin, it should be reset to Smin. Therefore, the minimum number p exists but the maximum number does not. The cost of building and maintaining supply centers will be greatly wasted if the total capacity of supply centers exceeds the actual expected demands of customers, especially in the EVM and its crisp variant model. Fixed value A, which denotes building a supply center, is the decision factor in choosing the best number p for the optimal solution (discussed in Section 4). For each p, the main procedure of IPSO can be described as follows.Step 1.
Initialize each individual of the particle swarm population with a given sizepop_size as follows:(40)Indk=xk,yk,Sk1×3p=x1k,x2k,…,xpk,y1k,y2k,…,ypk,s1k,s2k,…,spk,k=1,2,…,pop_size,where Indk must satisfy the constraints (introduced in Section 3.3), which are generated randomly as in Algorithm 1. TD denotes the total demand of customers andE(TD) denotes the expected value equal to ∑j=1mej.Algorithm 1:The initialization algorithm.
Note: Initialize individuals with pop_size, p, Smin, Smax, and TD
(1) for all  k∈pop_sizedo
(2) forj=1→2pdo
(3)  ifj<=pthen
(4)   x←random(0,1)
(5)   NIndkj←x
(6)  else ifj>p&j<=2pthen
(7)   y←random(0,1)
(8)   NIndkj←y
(9)  end if
(10) end for
(11) whilesum(S)<E(TD)do
(12)  forj=1→pdo
(13)   Sj←random(Smin,Smax)
(14)  end for
(15) end while
//Normalize S and complete the construction of NInd
(16) S′←(S-Smin)/(Smax-Smin)
(17) forj=2p+1→3pdo
(18)  NIndkj←Sj-2p′
(19) end for
(20) end for
Output: NIndStep 2.
Encode the initial individuals into normalized values to update the particle swarm. The normalization formula is as follows:(41)NIndk=NormIndk=xk-xminxmax-xmin,yk-yminymax-ymin,Sk-SminSmax-Smin1×3p=x′k,y′k,S′k1×3p,where xmax and ymax denote the upper-bounds of the location coordinates and xmin and ymin denote the lower-bounds of the location coordinates. Specifically, we set the position coordinate bounds as follows:(42)xmin=minj∈M⁡aj,xmax=minj∈M⁡aj;ymin=minj∈M⁡bj,ymax=minj∈M⁡bj.Then, the particle individuals are normalized, and NIndk is equal to the position of particle k in the continuous search space.Step 3.
Set the parameters for the IPSO algorithm operating at each iteration. These parameters are listed in Table1.
w
,
c
1, and c2 change adaptively at each iteration. This phenomenon is called adaptive unfixed phenomenon, which is discussed in Section 3.2. The positions are set between 0 and 1 because the individuals have been normalized. The speed limitation can be set to any reasonable symmetric interval. It is set at [–0.3,0.3]. The dimensions of the position and the speed vectors are both equal to 3p.Table 1
Description of parameters in the IPSO algorithm.
Notation
Name
Range
State
w
Inertia weight
0
<
w
<
1
Adaptive unfixed
c
1
,
c
2
Acceleration factor
c
1
,
c
2
≥
0
Adaptive unfixed
[
X
m
i
n
,
X
m
a
x
]
D
=
3
p
Position limitation
X
m
i
n
=
0, Xmax=1
Fixed
[
V
m
i
n
,
V
m
a
x
]
D
=
3
p
Speed limitation
V
m
i
n
=
-
0.3, Vmax=0.3
FixedStep 4.
Check the constraints for each individual. If feasible, calculate the fitness values of the second and first stage objective functions. The second-stage programing is linear programing, which is solved by the interior point method and discussed in Section3.3. Otherwise, regenerate the individual using the initialization process in Algorithm 1.Step 5.
Update the position and speed of particles using different strategies. If the global optimum does not improve within a given cumulative iterationΔT=ΔIter, theRestart strategy is activated to reinitialize all individuals based on the global best-so-far solution and the Gaussian cloud model. Otherwise, theGaussian cloud operator is activated to update the particles based on the global best-so-far solution and the individual best-so-far solution. The two updating strategies are discussed in Section 3.2.Step 6.
Repeat Steps3–5 for IterMax.Step 7.
Decode the global best-so-far solution into the initial range and generate the global optimum. The inverse normalization formula is as follows:(43)Indk=Norm-1NIndk=x′k·xmax-xmin+xmin,y′k·ymax-ymin+ymin,S′k·Smax-Smin+Smin1×3p=xk,yk,Sk1×3p.


### 3.2. Improvement Strategies

#### 3.2.1. Gaussian Cloud Operator

TheGaussian cloud operator is an algorithm application of the membership clouds. The membership cloud was first proposed by Deyi et al. [39], a member of the Chinese Academy of Sciences. The method bridges the gap between quantitative methodology and qualitative methodology based on the fuzzy set theory and has been successfully applied in algorithm improvement [38, 40]. A classic Gaussian cloud is described in Figure 2.Figure 2
Classic Gaussian cloud with expectation Ex, entropy En, and hyperentropy He.Assuming∀k∈pop_size, we have Xk=NIndk, and the particle population is X=X1,X2,…,Xpop_sizeT. Each particle position vector Xk=xk1,xk2,…,xkD in a D-dimension solution space (D=3p) represents a potential solution, which is accordingly associated with a speed vector Vk=vk1,vk2,…,vkD. The individual best-so-far solution of each particle can be expressed as Pk=pk1,pk2,…,pkD and the global best-so-far solution can be expressed as Pg=pg1,pg2,…,pgD based on the previous fitness value of each particle. In a traditional PSO updating process, Pk and Pg are used as heuristic factors in (44)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2pgdt-xkdt,(45)xkdt+1=xkdt+vkdt+1,where d=1,2,…,D, and other notations are shown in Table 1. However, with theGaussian cloud operator, (39) will be modified as (46)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2gcpgdt-xkdt,where gcpgd=GC(Ex,En,He), Ex=pgd, and x=GC(Ex,En,He) are the cloud drop generator functions, which can be described as follows:Step 1.
Generate a normal random numberσ based on the expectation value En, and the standard deviation value is He, that is, σ~N(En,He).Step 2.
Generate a normal random numberx based on the expectation value Ex, and the standard deviation value is σ; that is, x~N(Ex,σ).We do not replace each dimension of initialPg with cloud drop gcpgd because it may induce excessive noise and disturbance. The reasonable approach is to replace one random dth-dimension of the initial Pg at each iteration.


#### 3.2.2. Restart Strategy

TheRestart strategyis another method to keep the algorithm from premature convergence. If the global optimum Pg does not improve within ΔT time (ΔT=ΔIter), that is, PgT=PgT+ΔT, the particle population is disrupted and reinitialized based on Pg. In ∀k∈pop_size, we also use theGaussian cloud operator to reinitialize the different input parameters as follows:(47)xkd=GCpgd,En′,He′,where En′=0.2 and He′=0.02. Subsequently, the position vector should be checked and modified, as in Algorithm 1, to satisfy the limitation constraint.


#### 3.2.3. Adaptive Parameter Strategy

Parameter setting always significantly affects the intelligent algorithm performance. The main related parameters of the IPSO algorithm are as follows: the number of iterationIterMax, the population size (pop_size), the inertia weight (w), the acceleration factors (c1 and c2), the entropy (En), and the hyperentropy (He) in theGaussian cloud operator. IterMax and pop_size are relatively fixed and have strong effects on computing time. Hence, we only utilize the adaptive parameter strategies on w, c1, c2, En, and He.The inertia weight (w) represents the particle inheritability of the previous speed, that is, a larger w means a stronger ability in global searching, whereas a smaller w means a stronger ability in local searching. A linear decreasing inertia weight (LDIW), which was first proposed by Shi and Eberhart [41], is adopted in (48)wIter=wmax-wmax-wminIterIterMax,where wmin is set to 0.4 and wmax is set to 0.9.The acceleration factors (c1 and c2) represent the particle weights on the attractiveness of Pk and Pg, respectively, which also reflect the balance between global and local searching. A typical setting is c1=c2=2. Several studies also proposed other effective settings, such as c1=c2=1.49 and c1=2.8, c2=1.3 [42]. We propose a linear adaptive parameter strategy that can be formulated as shown in (49)c1=2.5-IterIterMax,c2=3.5-c1.As a result, the algorithm will have a larger c1 and a smaller c2 at the beginning, which benefit global searching, and a smaller c1 and a larger c2 in the latter stages, which benefit local searching.The two important parameters in theGaussian cloud operator, which are called the entropy En and the hyperentropy He, could also be changed adaptively. The mutation of a random dimension value is more beneficial at the beginning of the iteration than at the later time. Hence, we propose a nonlinear adaptive parameter strategy for En and He as follows:(50)En=En01-IterIterMax2,He=0.1En,where En0 is set to 0.05 based on multiple tests.


### 3.3. Interior Point Method

To rapidly solve the second-stage programing ((29) and (30)), we adopt the interior point method instead of the simplex method in this paper. The interior point method was first proposed by Karmarkar [25] in 1984. It is a polynomial algorithm for linear programing that requires O(n3.5L) arithmetic operations on O(L) bit numbers in the worst case, where n is the number of variables and L is the number of bits in the input. Compared with the simplex method, it is very suitable for solving large-scale problems because the computing time does not increase with the size of the problem. Moreover, it is not strict with the initial point requirements of quadratic convergence and robustness. The improvement of the algorithm efficiency is verified in Section 4.3.


## 3.1. Design of the IPSO

In 1995, Kennedy and Eberhart [26] first proposed a swarm intelligence optimization algorithm, which was inspired by the flocking of birds and schooling of fish against predation, called the particle swarm optimization (PSO). PSO has been successfully applied to parameter optimization [27], combinatorial optimization [28], pattern recognition [29, 30], data mining [31], and other fields [32, 33]. However, local optimal solution is often obtained when general PSO algorithm is applied to solve the continuous location of distribution centers. Therefore, scholars have proposed many improvement strategies for a much efficient and effective PSO [34–37]. Recently, Zhan et al. [38] proposed an improved PSO based on neighbor heuristic and Gaussian cloud learning, and the results proved its superiority over many PSO variants. We also adopt a Gaussian cloud operator combined with an adaptive parameter strategy and a Restart strategy to improve the general PSO algorithm of the two-stage SCLAP. For a detailed description of the IPSO algorithm, an algorithm flow chart is provided as a supplement in Figure 1.Figure 1
IPSO algorithm flow chart.In Figure1, the loop program is repeated for IterMax times, and three improvement strategies, which are shown in the curved boxes, are adopted. If no better solution is found during ΔT(ΔT=ΔIter), theRestart strategy will be activated; otherwise, theGaussian cloud operator will be activated. In addition, several IPSO parameters will adaptively change along with the iteration (theadaptive parameter strategy).Based on (27), the range of supply center number p can be determined by (39)p≥∑j=1mejSmax+1,where [x] denotes taking the integer portion of x. The capacity of each supply center must be within [Smin,Smax], which means that if a candidate solution of si exceeds Smax, a penalty will be given in the general cost function, and if si is less than Smin, it should be reset to Smin. Therefore, the minimum number p exists but the maximum number does not. The cost of building and maintaining supply centers will be greatly wasted if the total capacity of supply centers exceeds the actual expected demands of customers, especially in the EVM and its crisp variant model. Fixed value A, which denotes building a supply center, is the decision factor in choosing the best number p for the optimal solution (discussed in Section 4). For each p, the main procedure of IPSO can be described as follows.Step 1.
Initialize each individual of the particle swarm population with a given sizepop_size as follows:(40)Indk=xk,yk,Sk1×3p=x1k,x2k,…,xpk,y1k,y2k,…,ypk,s1k,s2k,…,spk,k=1,2,…,pop_size,where Indk must satisfy the constraints (introduced in Section 3.3), which are generated randomly as in Algorithm 1. TD denotes the total demand of customers andE(TD) denotes the expected value equal to ∑j=1mej.Algorithm 1:The initialization algorithm.
Note: Initialize individuals with pop_size, p, Smin, Smax, and TD
(1) for all  k∈pop_sizedo
(2) forj=1→2pdo
(3)  ifj<=pthen
(4)   x←random(0,1)
(5)   NIndkj←x
(6)  else ifj>p&j<=2pthen
(7)   y←random(0,1)
(8)   NIndkj←y
(9)  end if
(10) end for
(11) whilesum(S)<E(TD)do
(12)  forj=1→pdo
(13)   Sj←random(Smin,Smax)
(14)  end for
(15) end while
//Normalize S and complete the construction of NInd
(16) S′←(S-Smin)/(Smax-Smin)
(17) forj=2p+1→3pdo
(18)  NIndkj←Sj-2p′
(19) end for
(20) end for
Output: NIndStep 2.
Encode the initial individuals into normalized values to update the particle swarm. The normalization formula is as follows:(41)NIndk=NormIndk=xk-xminxmax-xmin,yk-yminymax-ymin,Sk-SminSmax-Smin1×3p=x′k,y′k,S′k1×3p,where xmax and ymax denote the upper-bounds of the location coordinates and xmin and ymin denote the lower-bounds of the location coordinates. Specifically, we set the position coordinate bounds as follows:(42)xmin=minj∈M⁡aj,xmax=minj∈M⁡aj;ymin=minj∈M⁡bj,ymax=minj∈M⁡bj.Then, the particle individuals are normalized, and NIndk is equal to the position of particle k in the continuous search space.Step 3.
Set the parameters for the IPSO algorithm operating at each iteration. These parameters are listed in Table1.
w
,
c
1, and c2 change adaptively at each iteration. This phenomenon is called adaptive unfixed phenomenon, which is discussed in Section 3.2. The positions are set between 0 and 1 because the individuals have been normalized. The speed limitation can be set to any reasonable symmetric interval. It is set at [–0.3,0.3]. The dimensions of the position and the speed vectors are both equal to 3p.Table 1
Description of parameters in the IPSO algorithm.
Notation
Name
Range
State
w
Inertia weight
0
<
w
<
1
Adaptive unfixed
c
1
,
c
2
Acceleration factor
c
1
,
c
2
≥
0
Adaptive unfixed
[
X
m
i
n
,
X
m
a
x
]
D
=
3
p
Position limitation
X
m
i
n
=
0, Xmax=1
Fixed
[
V
m
i
n
,
V
m
a
x
]
D
=
3
p
Speed limitation
V
m
i
n
=
-
0.3, Vmax=0.3
FixedStep 4.
Check the constraints for each individual. If feasible, calculate the fitness values of the second and first stage objective functions. The second-stage programing is linear programing, which is solved by the interior point method and discussed in Section3.3. Otherwise, regenerate the individual using the initialization process in Algorithm 1.Step 5.
Update the position and speed of particles using different strategies. If the global optimum does not improve within a given cumulative iterationΔT=ΔIter, theRestart strategy is activated to reinitialize all individuals based on the global best-so-far solution and the Gaussian cloud model. Otherwise, theGaussian cloud operator is activated to update the particles based on the global best-so-far solution and the individual best-so-far solution. The two updating strategies are discussed in Section 3.2.Step 6.
Repeat Steps3–5 for IterMax.Step 7.
Decode the global best-so-far solution into the initial range and generate the global optimum. The inverse normalization formula is as follows:(43)Indk=Norm-1NIndk=x′k·xmax-xmin+xmin,y′k·ymax-ymin+ymin,S′k·Smax-Smin+Smin1×3p=xk,yk,Sk1×3p.


## 3.2. Improvement Strategies

### 3.2.1. Gaussian Cloud Operator

TheGaussian cloud operator is an algorithm application of the membership clouds. The membership cloud was first proposed by Deyi et al. [39], a member of the Chinese Academy of Sciences. The method bridges the gap between quantitative methodology and qualitative methodology based on the fuzzy set theory and has been successfully applied in algorithm improvement [38, 40]. A classic Gaussian cloud is described in Figure 2.Figure 2
Classic Gaussian cloud with expectation Ex, entropy En, and hyperentropy He.Assuming∀k∈pop_size, we have Xk=NIndk, and the particle population is X=X1,X2,…,Xpop_sizeT. Each particle position vector Xk=xk1,xk2,…,xkD in a D-dimension solution space (D=3p) represents a potential solution, which is accordingly associated with a speed vector Vk=vk1,vk2,…,vkD. The individual best-so-far solution of each particle can be expressed as Pk=pk1,pk2,…,pkD and the global best-so-far solution can be expressed as Pg=pg1,pg2,…,pgD based on the previous fitness value of each particle. In a traditional PSO updating process, Pk and Pg are used as heuristic factors in (44)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2pgdt-xkdt,(45)xkdt+1=xkdt+vkdt+1,where d=1,2,…,D, and other notations are shown in Table 1. However, with theGaussian cloud operator, (39) will be modified as (46)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2gcpgdt-xkdt,where gcpgd=GC(Ex,En,He), Ex=pgd, and x=GC(Ex,En,He) are the cloud drop generator functions, which can be described as follows:Step 1.
Generate a normal random numberσ based on the expectation value En, and the standard deviation value is He, that is, σ~N(En,He).Step 2.
Generate a normal random numberx based on the expectation value Ex, and the standard deviation value is σ; that is, x~N(Ex,σ).We do not replace each dimension of initialPg with cloud drop gcpgd because it may induce excessive noise and disturbance. The reasonable approach is to replace one random dth-dimension of the initial Pg at each iteration.


### 3.2.2. Restart Strategy

TheRestart strategyis another method to keep the algorithm from premature convergence. If the global optimum Pg does not improve within ΔT time (ΔT=ΔIter), that is, PgT=PgT+ΔT, the particle population is disrupted and reinitialized based on Pg. In ∀k∈pop_size, we also use theGaussian cloud operator to reinitialize the different input parameters as follows:(47)xkd=GCpgd,En′,He′,where En′=0.2 and He′=0.02. Subsequently, the position vector should be checked and modified, as in Algorithm 1, to satisfy the limitation constraint.


### 3.2.3. Adaptive Parameter Strategy

Parameter setting always significantly affects the intelligent algorithm performance. The main related parameters of the IPSO algorithm are as follows: the number of iterationIterMax, the population size (pop_size), the inertia weight (w), the acceleration factors (c1 and c2), the entropy (En), and the hyperentropy (He) in theGaussian cloud operator. IterMax and pop_size are relatively fixed and have strong effects on computing time. Hence, we only utilize the adaptive parameter strategies on w, c1, c2, En, and He.The inertia weight (w) represents the particle inheritability of the previous speed, that is, a larger w means a stronger ability in global searching, whereas a smaller w means a stronger ability in local searching. A linear decreasing inertia weight (LDIW), which was first proposed by Shi and Eberhart [41], is adopted in (48)wIter=wmax-wmax-wminIterIterMax,where wmin is set to 0.4 and wmax is set to 0.9.The acceleration factors (c1 and c2) represent the particle weights on the attractiveness of Pk and Pg, respectively, which also reflect the balance between global and local searching. A typical setting is c1=c2=2. Several studies also proposed other effective settings, such as c1=c2=1.49 and c1=2.8, c2=1.3 [42]. We propose a linear adaptive parameter strategy that can be formulated as shown in (49)c1=2.5-IterIterMax,c2=3.5-c1.As a result, the algorithm will have a larger c1 and a smaller c2 at the beginning, which benefit global searching, and a smaller c1 and a larger c2 in the latter stages, which benefit local searching.The two important parameters in theGaussian cloud operator, which are called the entropy En and the hyperentropy He, could also be changed adaptively. The mutation of a random dimension value is more beneficial at the beginning of the iteration than at the later time. Hence, we propose a nonlinear adaptive parameter strategy for En and He as follows:(50)En=En01-IterIterMax2,He=0.1En,where En0 is set to 0.05 based on multiple tests.


## 3.2.1. Gaussian Cloud Operator

TheGaussian cloud operator is an algorithm application of the membership clouds. The membership cloud was first proposed by Deyi et al. [39], a member of the Chinese Academy of Sciences. The method bridges the gap between quantitative methodology and qualitative methodology based on the fuzzy set theory and has been successfully applied in algorithm improvement [38, 40]. A classic Gaussian cloud is described in Figure 2.Figure 2
Classic Gaussian cloud with expectation Ex, entropy En, and hyperentropy He.Assuming∀k∈pop_size, we have Xk=NIndk, and the particle population is X=X1,X2,…,Xpop_sizeT. Each particle position vector Xk=xk1,xk2,…,xkD in a D-dimension solution space (D=3p) represents a potential solution, which is accordingly associated with a speed vector Vk=vk1,vk2,…,vkD. The individual best-so-far solution of each particle can be expressed as Pk=pk1,pk2,…,pkD and the global best-so-far solution can be expressed as Pg=pg1,pg2,…,pgD based on the previous fitness value of each particle. In a traditional PSO updating process, Pk and Pg are used as heuristic factors in (44)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2pgdt-xkdt,(45)xkdt+1=xkdt+vkdt+1,where d=1,2,…,D, and other notations are shown in Table 1. However, with theGaussian cloud operator, (39) will be modified as (46)vkdt+1=wvkdt+c1r1pkdt-xkdt+c2r2gcpgdt-xkdt,where gcpgd=GC(Ex,En,He), Ex=pgd, and x=GC(Ex,En,He) are the cloud drop generator functions, which can be described as follows:Step 1.
Generate a normal random numberσ based on the expectation value En, and the standard deviation value is He, that is, σ~N(En,He).Step 2.
Generate a normal random numberx based on the expectation value Ex, and the standard deviation value is σ; that is, x~N(Ex,σ).We do not replace each dimension of initialPg with cloud drop gcpgd because it may induce excessive noise and disturbance. The reasonable approach is to replace one random dth-dimension of the initial Pg at each iteration.


## 3.2.2. Restart Strategy

TheRestart strategyis another method to keep the algorithm from premature convergence. If the global optimum Pg does not improve within ΔT time (ΔT=ΔIter), that is, PgT=PgT+ΔT, the particle population is disrupted and reinitialized based on Pg. In ∀k∈pop_size, we also use theGaussian cloud operator to reinitialize the different input parameters as follows:(47)xkd=GCpgd,En′,He′,where En′=0.2 and He′=0.02. Subsequently, the position vector should be checked and modified, as in Algorithm 1, to satisfy the limitation constraint.


## 3.2.3. Adaptive Parameter Strategy

Parameter setting always significantly affects the intelligent algorithm performance. The main related parameters of the IPSO algorithm are as follows: the number of iterationIterMax, the population size (pop_size), the inertia weight (w), the acceleration factors (c1 and c2), the entropy (En), and the hyperentropy (He) in theGaussian cloud operator. IterMax and pop_size are relatively fixed and have strong effects on computing time. Hence, we only utilize the adaptive parameter strategies on w, c1, c2, En, and He.The inertia weight (w) represents the particle inheritability of the previous speed, that is, a larger w means a stronger ability in global searching, whereas a smaller w means a stronger ability in local searching. A linear decreasing inertia weight (LDIW), which was first proposed by Shi and Eberhart [41], is adopted in (48)wIter=wmax-wmax-wminIterIterMax,where wmin is set to 0.4 and wmax is set to 0.9.The acceleration factors (c1 and c2) represent the particle weights on the attractiveness of Pk and Pg, respectively, which also reflect the balance between global and local searching. A typical setting is c1=c2=2. Several studies also proposed other effective settings, such as c1=c2=1.49 and c1=2.8, c2=1.3 [42]. We propose a linear adaptive parameter strategy that can be formulated as shown in (49)c1=2.5-IterIterMax,c2=3.5-c1.As a result, the algorithm will have a larger c1 and a smaller c2 at the beginning, which benefit global searching, and a smaller c1 and a larger c2 in the latter stages, which benefit local searching.The two important parameters in theGaussian cloud operator, which are called the entropy En and the hyperentropy He, could also be changed adaptively. The mutation of a random dimension value is more beneficial at the beginning of the iteration than at the later time. Hence, we propose a nonlinear adaptive parameter strategy for En and He as follows:(50)En=En01-IterIterMax2,He=0.1En,where En0 is set to 0.05 based on multiple tests.


## 3.3. Interior Point Method

To rapidly solve the second-stage programing ((29) and (30)), we adopt the interior point method instead of the simplex method in this paper. The interior point method was first proposed by Karmarkar [25] in 1984. It is a polynomial algorithm for linear programing that requires O(n3.5L) arithmetic operations on O(L) bit numbers in the worst case, where n is the number of variables and L is the number of bits in the input. Compared with the simplex method, it is very suitable for solving large-scale problems because the computing time does not increase with the size of the problem. Moreover, it is not strict with the initial point requirements of quadratic convergence and robustness. The improvement of the algorithm efficiency is verified in Section 4.3.


## 4. Numerical Example

### 4.1. Data and Implementation

In this section, a 20-customer network of allocation and transportation in emergency logistics is designed, and the decision objective is to determine the number, capacities, and locations of supply centers for material support. The customer locations and their stochastic demands are listed in Table2.Table 2
Customer locations and stochastic demand distributions.
Number
Locations
Demand distributions
(1)
(28, 42)
N
(
2,1
)
(2)
(18, 50)
N
(
4,1
)
(3)
(74, 34)
N
(
7,1
)
(4)
(74, 6)
N
(
6,1
)
(5)
(70, 18)
N
(
4,1
)
(6)
(72, 98)
N
(
10,1
)
(7)
(60, 50)
N
(
4,1
)
(8)
(36, 40)
N
(
2,1
)
(9)
(12, 4)
N
(
5,1
)
(10)
(18, 20)
N
(
10,1
)
(11)
(14, 78)
N
(
2,1
)
(12)
(90, 36)
N
(
6,1
)
(13)
(78, 20)
N
(
7,1
)
(14)
(24, 52)
N
(
10,1
)
(15)
(54, 6)
N
(
8,1
)
(16)
(62, 60)
N
(
7,1
)
(17)
(98, 14)
N
(
2,1
)
(18)
(36, 58)
N
(
6,1
)
(19)
(38, 88)
N
(
9,1
)
(20)
(32, 54)
N
(
7,1
)We adopt the stochastic simulation method [7] to solve the two-stage SCLAP, and the simulation number of time R is set to 100 for a fair comparison with the traditional EVM. The generalized cost is significantly affected by fixed cost A and variable cost coefficient δ. Hence, we test eight pairs of different values for comparison and analysis. The eight pairs are divided into two classes using the controlling variable method for testing. The values are listed in Table 3.Table 3
Fixed costA and variable cost coefficient δ.
Class 1
A
200
300
500
1000
δ
0.1
0.1
0.1
0.1
Class 2
A
400
400
400
400
δ
0.1
0.2
0.4
0.7The experiments are conducted on a laptop with Intel® Core™ i7-6700HQ 2.60 GHz quad-core processors, and the algorithms are coded in MATLAB R2016a environment. For convenience, we adopt the built-in functionlinprogin MATLAB to achieve linear programing with the simplex method or the interior point method for the second-stage optimization problem.


### 4.2. Computational Results

#### 4.2.1. Results of the Deterministic Model

We can achieve an optimal solution for the deterministic model from the initial EVM using the IPSO.pop_size is set to 10 and IterMax is set to 500 for the IPSO algorithm. The capacity limitation of supply centers is set as Smin=30 and Smax=100. The fixed cost and the variable cost coefficient for the generalized cost function are set as A = 500 and δ=0.1, respectively. The results are listed in Table 4 where the number of supply centers p will accordingly change from two to four. Cap. denotes the capacity and Pos. denotes the position of each supply center. GC denotes the generalized cost of each solution with different numbers of centers and CT denotes the computing time of each situation.Table 4
Results of the deterministic model.
p
Cap.
Pos.
GC
CT
2
78.000
(32.002, 54.000)
3775.291
29.57 s
40.014
(76.318, 19.881)
3
30.000
(57.580, 86.522)
3643.147
31.33 s
44.000
(75.068, 20.762)
48.000
(25.065, 49.983)
4
30.000
(57.521, 86.513)
3655.400
33.63 s
40.000
(76.379, 19.922)
30.000
(18.000, 20.000)
37.000
(31.922, 53.975)The best decision solution is to build three supply centers, and the optimum generalized cost is 3643.147. We also present the convergence curves of differentp numbers in Figure 3 and illustrate the solutions with LA maps in Figures 4–6. In addition, the amount of traffic is simply represented with different shades of color.Figure 3
Convergence curves with different number of supply centers.Figure 4
LA map withp = 2.Figure 5
LA map withp = 3.Figure 6
LA map withp = 4.In the traditional one-stage problem, the number and capacity of a supply center are fixed, whereas the number (p) and capacity (Cap.) are decision variables in this problem. Moreover, the capacity has a range. The results of the two-stage SCLAP can provide significant practical conclusions for actual emergency logistics because the scarcity of supply materials is extremely serious and urgent such that the number and capacity of supply centers must be determined accurately based on practical situations. Therefore, we must achieve a balance between minimizing the number of centers and minimizing the distance of transport.


#### 4.2.2. Results of the EVM with Stochastic Simulation

Stochastic simulation is time-consuming and extremely sensitive to the number of simulation iterations (R). Only when R is sufficiently large can the solution become stable. As a result, R is set to 100 and a parallel strategy, that is, aparfor loop with the MATLAB software, is adopted. For each simulation, the random demands of all customers are generated, and the programing model is solved. The IPSO obtains a stable optimal solution based on the expected value of the generalized cost at each iteration. The results of the EVM with stochastic simulation are listed in Table 5.Table 5
Results of the EVM with stochastic simulation.
p
Cap.
Pos.
GC
CT
2
40.123
(76.615, 19.827)
3785.601
1010.467 s
78.093
(32.578, 55.244)
3
30.000
(55.178, 86.894)
3660.836
1035.310 s
40.034
(76.309, 19.186)
52.326
(25.479, 48.314)
4
37.128
(29.474, 52.955)
3674.496
1055.373 s
30.000
(56.231, 89.416)
30.000
(18.218, 18.944)
40.071
(76.446, 18.986)Table5 shows that the best decision solution is also to build three supply centers, and the optimum generalized cost is 3660.836. The GC values are slightly larger than those in Table 4 because the simulation inevitably induces several random factors. Thus, obtaining the theoretical optimal solution is impossible. Cap. and Pos. are also different than those in Table 4 for the same reason. However, the CT values are significantly larger than those in the deterministic model, which implies that the deterministic model can dramatically improve the computational efficiency.


#### 4.2.3. Parameter Analysis of Generalized Cost

The generalized cost of the optimal solution is dramatically affected by fixed costA and variable cost coefficient δ. We design two experiments based on the data in Table 3 to demonstrate this phenomenon, and the number of supply centers (p) is set from 2 to 8 for each class testing. The experiments are carried out by the IPSO algorithm with the same parameters for 10 times. The best results are recorded in Table 6 and shown in Figures 7 and 8.Table 6
Best results of the 10 tests with differentA or δ.
p
Class 1 (δ=0.1)
Class 2 (A=400)
A
=
200
A
=
300
A
=
500
A
=
1000
δ
=
0.1
δ
=
0.2
δ
=
0.4
δ
=
0.7
2
3180.616
3378.905
3775.291
4768.959
3582.563
3583.987
3613.299
3660.975
3
2744.854
3045.227
3643.147
5147.107
3344.534
3361.372
3382.718
3421.571
4
2460.016
2861.945
3655.400
5658.357
3259.020
3274.130
3302.576
3345.515
5
2353.784
2858.423
3854.877
6352.455
3349.796
3370.884
3408.197
3449.036
6
2290.738
2921.248
4107.793
7174.363
3451.416
3466.693
3505.663
3558.175
7
2313.417
3052.983
4516.247
7971.316
3671.614
3719.131
3775.696
3838.061
8
2424.430
3214.099
4881.651
8877.735
4031.114
4031.241
4062.632
4141.791Figure 7
Changes in GC with different fixed cost valuesA (δ=0.1).Figure 8
Changes in GC with variable cost coefficientδ (A = 400).In Table6, the bold values represent the minimum GCs under different p values. The most suitable number of supply centers should be built according to the given fixed cost value A and variable cost coefficient δ. Parameter A plays an important role in determining the best p. The best p increases when A decreases, as shown in Figure 7. In addition, the GC value nonlinearly increases with p. Therefore, if the number of supply centers exceeds the actual required number based on the actual total demand, a large cost will be wasted. With variable cost coefficient δ, the GC value increases linearly with p and has no effect on the best p decision, as shown in Table 6 and Figure 8. As a result, the most logical and economical approach is to reduce the fixed cost value of building each supply center.


### 4.3. Comparison and Analysis of the Algorithms

To provide a fair comparison of the IPSO and PSO algorithms, we use the same data in Table2 and the common parameter settings in Section 4.2.1 to solve the two-stage EVM. However, in the basic PSO algorithm, the fixed parameters are set as follows: w=0.9 and c1=c2=2. The test is repeated 10 times, and the statistical results are presented in Table 7. Avg GC denotes the average GC value of the 10 tests, and the Best GC denotes the best value of these tests. Avg Time denotes the average running time. Table 7 shows that the IPSO algorithm is better than the basic PSO algorithm in both convergence accuracy and efficiency, which means that the improvement strategies are effective.Table 7
Statistical results of the different algorithms used to solve the two-stage EVM.
Algorithm
Avg GC
Best GC
Avg Time
PSO
3664.726
3651.823
30.60 s
IPSO
3646.752
3643.318
30.33 s
I
P
S
O
∗
3648.000
3639.652
83.91 s
∗ denotes the use of the simplex method.The interior point method is adopted in the IPSO algorithm. However, its advantage over the simplex method must also be verified. As a result, we test the IPSO algorithm with the simplex method in solving the second-stage programing of the deterministic model. The results are listed in row 3 of Table7 with an asterisk. The Avg GC and Best GC results obtained by the simplex method are close to the results of the interior point method. However, Avg Time shows that the interior point method is much better than the simplex method in computing efficiency.To further illustrate the applicability and superiority of the proposed IPSO algorithm, additional tests, based on the classic one-stage EVM, are performed. Solving the classic one-stage EVM is the same as solving the second stage of a two-stage EVM with knownp and S. As a result, the initial individuals are recoded and Indk=xk,yk,Sk1×3p is changed to Indk=xk,yk1×3p∗, where p∗ is the given number of supply centers, which is 4, and S becomes a known variable [40,50,60,70]. Zhou and Liu [7] developed a hybrid intelligent algorithm (HIA) that consists of a network simplex algorithm, a simulation, and a genetic algorithm to solve this model. HIA is used to compare the proposed IPSO algorithm and the classic PSO and IPSO∗ algorithms. The PSO, IPSO, and IPSO∗ algorithms share the same parameter setting with the previous test. HIA parameters are set as in [7]; that is, the crossover probability is 0.3, the mutation probability is 0.2, the population size is 10, the iteration number is 500, and the simulation number is 500. The numbers of customers and locations in the test data are completely consistent with the data in Table 2. However, the demand distributions are different. The expected demand for the first 10 customers is 5, that for the last 10 customers is 10, and the customer demand variance is 1. The test is repeated 10 times, and the results are recorded in Table 8. Then, some statistical results are calculated and presented in Table 9, where Best GC denotes the minimum value of general cost as shown in (20), Avg GCE denotes the average percent error of the general cost, Best CT denotes the minimum value of computing time, and Avg CTE denotes the average percent error of computing time. The percent error is calculated as (actual value − optimal value)/optimal value × 100%, which can be used to show the deviation from the optimal solution.Table 8
Results of the 10 tests of different algorithms solving the classic one-stage EVM.
PSO
IPSO
I
P
S
O
∗
HIA
GC
CT
GC
CT
GC
CT
GC
CT
1
2202
31.8 s
2151
31.5 s
2136
81.4 s
2194
66 min
2
2191
31.5 s
2151
30.2 s
2138
87.4 s
2185
65 min
3
2132
31.7 s
2151
31.1 s
2132
79.1 s
2192
68 min
4
2202
32.4 s
2132
31.9 s
2151
80.5 s
2158
71 min
5
2132
31.1 s
2132
32.2 s
2151
92.1 s
2239
65 min
6
2132
31.2 s
2151
30.9 s
2132
82.0 s
2173
69 min
7
2292
31.4 s
2151
31.3 s
2151
83.1 s
2204
66 min
8
2217
31.8 s
2132
31.3 s
2202
83.0 s
2153
69 min
9
2247
32.1 s
2132
32.5 s
2132
81.8 s
2214
68 min
10
2221
31.8 s
2132
31.4 s
2152
81.6 s
2169
69 min
∗ denotes the use of the simplex method.Table 9
Statistical results of different algorithms solving the classic one-stage EVM.
PSO
IPSO
I
P
S
O
∗
HIA
Best GC
2132
2132
2132
2153
Avg GCE (%)
3.04
0.45
0.74
1.63
Best CT
31.1 s
30.2 s
79.1 s
65 min
Avg CTE (%)
1.86
4.07
5.18
4.00
∗ denotes the use of the simplex method.The results in Tables8 and 9 show that the performance of the IPSO algorithm is the most superior in terms of convergence accuracy and convergence speed. In addition, on the basis of the Avg GCE and the Avg CTE, IPSO is the most stable algorithm among the four in terms of accuracy. It also performs stably in terms of computing time. Adopting the simulation method for HIA results in spending excessive time in the evaluation cycle. Hence, the computing time is much longer than those of the other three algorithms.In addition, the two tests also illustrate the difference between the two EVMs. The proposed two-stage EVM assumes that the supply center number, capacities, and locations are decision variables. This approach is applicable in emergency logistics because it avoids wasting emergency resources for inappropriate number of supply centers or capacities. However, in the classic one-stage EVM, the numbers of centers and capacities are fixed. If the fixed costA is set at 500 and the variable cost coefficient δ is set at 0.1, we can also obtain the same general cost with the two-stage EVM. The first stage objective in the two-stage EVM is 4147 when the optimal solution has a one-stage EVM general cost of 2132, as shown in Table 9. However, if we initially use the two-stage EVM, we can obtain the first stage objective optimal solution of 4076. The optimal decision is to set three supply centers with capacities of 30, 65, and 55, and the second-stage objective is 2561. Evidently, the total cost obtained by the two-stage EVM is lower and much reasonable than the classic one-stage EVM.


## 4.1. Data and Implementation

In this section, a 20-customer network of allocation and transportation in emergency logistics is designed, and the decision objective is to determine the number, capacities, and locations of supply centers for material support. The customer locations and their stochastic demands are listed in Table2.Table 2
Customer locations and stochastic demand distributions.
Number
Locations
Demand distributions
(1)
(28, 42)
N
(
2,1
)
(2)
(18, 50)
N
(
4,1
)
(3)
(74, 34)
N
(
7,1
)
(4)
(74, 6)
N
(
6,1
)
(5)
(70, 18)
N
(
4,1
)
(6)
(72, 98)
N
(
10,1
)
(7)
(60, 50)
N
(
4,1
)
(8)
(36, 40)
N
(
2,1
)
(9)
(12, 4)
N
(
5,1
)
(10)
(18, 20)
N
(
10,1
)
(11)
(14, 78)
N
(
2,1
)
(12)
(90, 36)
N
(
6,1
)
(13)
(78, 20)
N
(
7,1
)
(14)
(24, 52)
N
(
10,1
)
(15)
(54, 6)
N
(
8,1
)
(16)
(62, 60)
N
(
7,1
)
(17)
(98, 14)
N
(
2,1
)
(18)
(36, 58)
N
(
6,1
)
(19)
(38, 88)
N
(
9,1
)
(20)
(32, 54)
N
(
7,1
)We adopt the stochastic simulation method [7] to solve the two-stage SCLAP, and the simulation number of time R is set to 100 for a fair comparison with the traditional EVM. The generalized cost is significantly affected by fixed cost A and variable cost coefficient δ. Hence, we test eight pairs of different values for comparison and analysis. The eight pairs are divided into two classes using the controlling variable method for testing. The values are listed in Table 3.Table 3
Fixed costA and variable cost coefficient δ.
Class 1
A
200
300
500
1000
δ
0.1
0.1
0.1
0.1
Class 2
A
400
400
400
400
δ
0.1
0.2
0.4
0.7The experiments are conducted on a laptop with Intel® Core™ i7-6700HQ 2.60 GHz quad-core processors, and the algorithms are coded in MATLAB R2016a environment. For convenience, we adopt the built-in functionlinprogin MATLAB to achieve linear programing with the simplex method or the interior point method for the second-stage optimization problem.


## 4.2. Computational Results

### 4.2.1. Results of the Deterministic Model

We can achieve an optimal solution for the deterministic model from the initial EVM using the IPSO.pop_size is set to 10 and IterMax is set to 500 for the IPSO algorithm. The capacity limitation of supply centers is set as Smin=30 and Smax=100. The fixed cost and the variable cost coefficient for the generalized cost function are set as A = 500 and δ=0.1, respectively. The results are listed in Table 4 where the number of supply centers p will accordingly change from two to four. Cap. denotes the capacity and Pos. denotes the position of each supply center. GC denotes the generalized cost of each solution with different numbers of centers and CT denotes the computing time of each situation.Table 4
Results of the deterministic model.
p
Cap.
Pos.
GC
CT
2
78.000
(32.002, 54.000)
3775.291
29.57 s
40.014
(76.318, 19.881)
3
30.000
(57.580, 86.522)
3643.147
31.33 s
44.000
(75.068, 20.762)
48.000
(25.065, 49.983)
4
30.000
(57.521, 86.513)
3655.400
33.63 s
40.000
(76.379, 19.922)
30.000
(18.000, 20.000)
37.000
(31.922, 53.975)The best decision solution is to build three supply centers, and the optimum generalized cost is 3643.147. We also present the convergence curves of differentp numbers in Figure 3 and illustrate the solutions with LA maps in Figures 4–6. In addition, the amount of traffic is simply represented with different shades of color.Figure 3
Convergence curves with different number of supply centers.Figure 4
LA map withp = 2.Figure 5
LA map withp = 3.Figure 6
LA map withp = 4.In the traditional one-stage problem, the number and capacity of a supply center are fixed, whereas the number (p) and capacity (Cap.) are decision variables in this problem. Moreover, the capacity has a range. The results of the two-stage SCLAP can provide significant practical conclusions for actual emergency logistics because the scarcity of supply materials is extremely serious and urgent such that the number and capacity of supply centers must be determined accurately based on practical situations. Therefore, we must achieve a balance between minimizing the number of centers and minimizing the distance of transport.


### 4.2.2. Results of the EVM with Stochastic Simulation

Stochastic simulation is time-consuming and extremely sensitive to the number of simulation iterations (R). Only when R is sufficiently large can the solution become stable. As a result, R is set to 100 and a parallel strategy, that is, aparfor loop with the MATLAB software, is adopted. For each simulation, the random demands of all customers are generated, and the programing model is solved. The IPSO obtains a stable optimal solution based on the expected value of the generalized cost at each iteration. The results of the EVM with stochastic simulation are listed in Table 5.Table 5
Results of the EVM with stochastic simulation.
p
Cap.
Pos.
GC
CT
2
40.123
(76.615, 19.827)
3785.601
1010.467 s
78.093
(32.578, 55.244)
3
30.000
(55.178, 86.894)
3660.836
1035.310 s
40.034
(76.309, 19.186)
52.326
(25.479, 48.314)
4
37.128
(29.474, 52.955)
3674.496
1055.373 s
30.000
(56.231, 89.416)
30.000
(18.218, 18.944)
40.071
(76.446, 18.986)Table5 shows that the best decision solution is also to build three supply centers, and the optimum generalized cost is 3660.836. The GC values are slightly larger than those in Table 4 because the simulation inevitably induces several random factors. Thus, obtaining the theoretical optimal solution is impossible. Cap. and Pos. are also different than those in Table 4 for the same reason. However, the CT values are significantly larger than those in the deterministic model, which implies that the deterministic model can dramatically improve the computational efficiency.


### 4.2.3. Parameter Analysis of Generalized Cost

The generalized cost of the optimal solution is dramatically affected by fixed costA and variable cost coefficient δ. We design two experiments based on the data in Table 3 to demonstrate this phenomenon, and the number of supply centers (p) is set from 2 to 8 for each class testing. The experiments are carried out by the IPSO algorithm with the same parameters for 10 times. The best results are recorded in Table 6 and shown in Figures 7 and 8.Table 6
Best results of the 10 tests with differentA or δ.
p
Class 1 (δ=0.1)
Class 2 (A=400)
A
=
200
A
=
300
A
=
500
A
=
1000
δ
=
0.1
δ
=
0.2
δ
=
0.4
δ
=
0.7
2
3180.616
3378.905
3775.291
4768.959
3582.563
3583.987
3613.299
3660.975
3
2744.854
3045.227
3643.147
5147.107
3344.534
3361.372
3382.718
3421.571
4
2460.016
2861.945
3655.400
5658.357
3259.020
3274.130
3302.576
3345.515
5
2353.784
2858.423
3854.877
6352.455
3349.796
3370.884
3408.197
3449.036
6
2290.738
2921.248
4107.793
7174.363
3451.416
3466.693
3505.663
3558.175
7
2313.417
3052.983
4516.247
7971.316
3671.614
3719.131
3775.696
3838.061
8
2424.430
3214.099
4881.651
8877.735
4031.114
4031.241
4062.632
4141.791Figure 7
Changes in GC with different fixed cost valuesA (δ=0.1).Figure 8
Changes in GC with variable cost coefficientδ (A = 400).In Table6, the bold values represent the minimum GCs under different p values. The most suitable number of supply centers should be built according to the given fixed cost value A and variable cost coefficient δ. Parameter A plays an important role in determining the best p. The best p increases when A decreases, as shown in Figure 7. In addition, the GC value nonlinearly increases with p. Therefore, if the number of supply centers exceeds the actual required number based on the actual total demand, a large cost will be wasted. With variable cost coefficient δ, the GC value increases linearly with p and has no effect on the best p decision, as shown in Table 6 and Figure 8. As a result, the most logical and economical approach is to reduce the fixed cost value of building each supply center.


## 4.2.1. Results of the Deterministic Model

We can achieve an optimal solution for the deterministic model from the initial EVM using the IPSO.pop_size is set to 10 and IterMax is set to 500 for the IPSO algorithm. The capacity limitation of supply centers is set as Smin=30 and Smax=100. The fixed cost and the variable cost coefficient for the generalized cost function are set as A = 500 and δ=0.1, respectively. The results are listed in Table 4 where the number of supply centers p will accordingly change from two to four. Cap. denotes the capacity and Pos. denotes the position of each supply center. GC denotes the generalized cost of each solution with different numbers of centers and CT denotes the computing time of each situation.Table 4
Results of the deterministic model.
p
Cap.
Pos.
GC
CT
2
78.000
(32.002, 54.000)
3775.291
29.57 s
40.014
(76.318, 19.881)
3
30.000
(57.580, 86.522)
3643.147
31.33 s
44.000
(75.068, 20.762)
48.000
(25.065, 49.983)
4
30.000
(57.521, 86.513)
3655.400
33.63 s
40.000
(76.379, 19.922)
30.000
(18.000, 20.000)
37.000
(31.922, 53.975)The best decision solution is to build three supply centers, and the optimum generalized cost is 3643.147. We also present the convergence curves of differentp numbers in Figure 3 and illustrate the solutions with LA maps in Figures 4–6. In addition, the amount of traffic is simply represented with different shades of color.Figure 3
Convergence curves with different number of supply centers.Figure 4
LA map withp = 2.Figure 5
LA map withp = 3.Figure 6
LA map withp = 4.In the traditional one-stage problem, the number and capacity of a supply center are fixed, whereas the number (p) and capacity (Cap.) are decision variables in this problem. Moreover, the capacity has a range. The results of the two-stage SCLAP can provide significant practical conclusions for actual emergency logistics because the scarcity of supply materials is extremely serious and urgent such that the number and capacity of supply centers must be determined accurately based on practical situations. Therefore, we must achieve a balance between minimizing the number of centers and minimizing the distance of transport.


## 4.2.2. Results of the EVM with Stochastic Simulation

Stochastic simulation is time-consuming and extremely sensitive to the number of simulation iterations (R). Only when R is sufficiently large can the solution become stable. As a result, R is set to 100 and a parallel strategy, that is, aparfor loop with the MATLAB software, is adopted. For each simulation, the random demands of all customers are generated, and the programing model is solved. The IPSO obtains a stable optimal solution based on the expected value of the generalized cost at each iteration. The results of the EVM with stochastic simulation are listed in Table 5.Table 5
Results of the EVM with stochastic simulation.
p
Cap.
Pos.
GC
CT
2
40.123
(76.615, 19.827)
3785.601
1010.467 s
78.093
(32.578, 55.244)
3
30.000
(55.178, 86.894)
3660.836
1035.310 s
40.034
(76.309, 19.186)
52.326
(25.479, 48.314)
4
37.128
(29.474, 52.955)
3674.496
1055.373 s
30.000
(56.231, 89.416)
30.000
(18.218, 18.944)
40.071
(76.446, 18.986)Table5 shows that the best decision solution is also to build three supply centers, and the optimum generalized cost is 3660.836. The GC values are slightly larger than those in Table 4 because the simulation inevitably induces several random factors. Thus, obtaining the theoretical optimal solution is impossible. Cap. and Pos. are also different than those in Table 4 for the same reason. However, the CT values are significantly larger than those in the deterministic model, which implies that the deterministic model can dramatically improve the computational efficiency.


## 4.2.3. Parameter Analysis of Generalized Cost

The generalized cost of the optimal solution is dramatically affected by fixed costA and variable cost coefficient δ. We design two experiments based on the data in Table 3 to demonstrate this phenomenon, and the number of supply centers (p) is set from 2 to 8 for each class testing. The experiments are carried out by the IPSO algorithm with the same parameters for 10 times. The best results are recorded in Table 6 and shown in Figures 7 and 8.Table 6
Best results of the 10 tests with differentA or δ.
p
Class 1 (δ=0.1)
Class 2 (A=400)
A
=
200
A
=
300
A
=
500
A
=
1000
δ
=
0.1
δ
=
0.2
δ
=
0.4
δ
=
0.7
2
3180.616
3378.905
3775.291
4768.959
3582.563
3583.987
3613.299
3660.975
3
2744.854
3045.227
3643.147
5147.107
3344.534
3361.372
3382.718
3421.571
4
2460.016
2861.945
3655.400
5658.357
3259.020
3274.130
3302.576
3345.515
5
2353.784
2858.423
3854.877
6352.455
3349.796
3370.884
3408.197
3449.036
6
2290.738
2921.248
4107.793
7174.363
3451.416
3466.693
3505.663
3558.175
7
2313.417
3052.983
4516.247
7971.316
3671.614
3719.131
3775.696
3838.061
8
2424.430
3214.099
4881.651
8877.735
4031.114
4031.241
4062.632
4141.791Figure 7
Changes in GC with different fixed cost valuesA (δ=0.1).Figure 8
Changes in GC with variable cost coefficientδ (A = 400).In Table6, the bold values represent the minimum GCs under different p values. The most suitable number of supply centers should be built according to the given fixed cost value A and variable cost coefficient δ. Parameter A plays an important role in determining the best p. The best p increases when A decreases, as shown in Figure 7. In addition, the GC value nonlinearly increases with p. Therefore, if the number of supply centers exceeds the actual required number based on the actual total demand, a large cost will be wasted. With variable cost coefficient δ, the GC value increases linearly with p and has no effect on the best p decision, as shown in Table 6 and Figure 8. As a result, the most logical and economical approach is to reduce the fixed cost value of building each supply center.


## 4.3. Comparison and Analysis of the Algorithms

To provide a fair comparison of the IPSO and PSO algorithms, we use the same data in Table2 and the common parameter settings in Section 4.2.1 to solve the two-stage EVM. However, in the basic PSO algorithm, the fixed parameters are set as follows: w=0.9 and c1=c2=2. The test is repeated 10 times, and the statistical results are presented in Table 7. Avg GC denotes the average GC value of the 10 tests, and the Best GC denotes the best value of these tests. Avg Time denotes the average running time. Table 7 shows that the IPSO algorithm is better than the basic PSO algorithm in both convergence accuracy and efficiency, which means that the improvement strategies are effective.Table 7
Statistical results of the different algorithms used to solve the two-stage EVM.
Algorithm
Avg GC
Best GC
Avg Time
PSO
3664.726
3651.823
30.60 s
IPSO
3646.752
3643.318
30.33 s
I
P
S
O
∗
3648.000
3639.652
83.91 s
∗ denotes the use of the simplex method.The interior point method is adopted in the IPSO algorithm. However, its advantage over the simplex method must also be verified. As a result, we test the IPSO algorithm with the simplex method in solving the second-stage programing of the deterministic model. The results are listed in row 3 of Table7 with an asterisk. The Avg GC and Best GC results obtained by the simplex method are close to the results of the interior point method. However, Avg Time shows that the interior point method is much better than the simplex method in computing efficiency.To further illustrate the applicability and superiority of the proposed IPSO algorithm, additional tests, based on the classic one-stage EVM, are performed. Solving the classic one-stage EVM is the same as solving the second stage of a two-stage EVM with knownp and S. As a result, the initial individuals are recoded and Indk=xk,yk,Sk1×3p is changed to Indk=xk,yk1×3p∗, where p∗ is the given number of supply centers, which is 4, and S becomes a known variable [40,50,60,70]. Zhou and Liu [7] developed a hybrid intelligent algorithm (HIA) that consists of a network simplex algorithm, a simulation, and a genetic algorithm to solve this model. HIA is used to compare the proposed IPSO algorithm and the classic PSO and IPSO∗ algorithms. The PSO, IPSO, and IPSO∗ algorithms share the same parameter setting with the previous test. HIA parameters are set as in [7]; that is, the crossover probability is 0.3, the mutation probability is 0.2, the population size is 10, the iteration number is 500, and the simulation number is 500. The numbers of customers and locations in the test data are completely consistent with the data in Table 2. However, the demand distributions are different. The expected demand for the first 10 customers is 5, that for the last 10 customers is 10, and the customer demand variance is 1. The test is repeated 10 times, and the results are recorded in Table 8. Then, some statistical results are calculated and presented in Table 9, where Best GC denotes the minimum value of general cost as shown in (20), Avg GCE denotes the average percent error of the general cost, Best CT denotes the minimum value of computing time, and Avg CTE denotes the average percent error of computing time. The percent error is calculated as (actual value − optimal value)/optimal value × 100%, which can be used to show the deviation from the optimal solution.Table 8
Results of the 10 tests of different algorithms solving the classic one-stage EVM.
PSO
IPSO
I
P
S
O
∗
HIA
GC
CT
GC
CT
GC
CT
GC
CT
1
2202
31.8 s
2151
31.5 s
2136
81.4 s
2194
66 min
2
2191
31.5 s
2151
30.2 s
2138
87.4 s
2185
65 min
3
2132
31.7 s
2151
31.1 s
2132
79.1 s
2192
68 min
4
2202
32.4 s
2132
31.9 s
2151
80.5 s
2158
71 min
5
2132
31.1 s
2132
32.2 s
2151
92.1 s
2239
65 min
6
2132
31.2 s
2151
30.9 s
2132
82.0 s
2173
69 min
7
2292
31.4 s
2151
31.3 s
2151
83.1 s
2204
66 min
8
2217
31.8 s
2132
31.3 s
2202
83.0 s
2153
69 min
9
2247
32.1 s
2132
32.5 s
2132
81.8 s
2214
68 min
10
2221
31.8 s
2132
31.4 s
2152
81.6 s
2169
69 min
∗ denotes the use of the simplex method.Table 9
Statistical results of different algorithms solving the classic one-stage EVM.
PSO
IPSO
I
P
S
O
∗
HIA
Best GC
2132
2132
2132
2153
Avg GCE (%)
3.04
0.45
0.74
1.63
Best CT
31.1 s
30.2 s
79.1 s
65 min
Avg CTE (%)
1.86
4.07
5.18
4.00
∗ denotes the use of the simplex method.The results in Tables8 and 9 show that the performance of the IPSO algorithm is the most superior in terms of convergence accuracy and convergence speed. In addition, on the basis of the Avg GCE and the Avg CTE, IPSO is the most stable algorithm among the four in terms of accuracy. It also performs stably in terms of computing time. Adopting the simulation method for HIA results in spending excessive time in the evaluation cycle. Hence, the computing time is much longer than those of the other three algorithms.In addition, the two tests also illustrate the difference between the two EVMs. The proposed two-stage EVM assumes that the supply center number, capacities, and locations are decision variables. This approach is applicable in emergency logistics because it avoids wasting emergency resources for inappropriate number of supply centers or capacities. However, in the classic one-stage EVM, the numbers of centers and capacities are fixed. If the fixed costA is set at 500 and the variable cost coefficient δ is set at 0.1, we can also obtain the same general cost with the two-stage EVM. The first stage objective in the two-stage EVM is 4147 when the optimal solution has a one-stage EVM general cost of 2132, as shown in Table 9. However, if we initially use the two-stage EVM, we can obtain the first stage objective optimal solution of 4076. The optimal decision is to set three supply centers with capacities of 30, 65, and 55, and the second-stage objective is 2561. Evidently, the total cost obtained by the two-stage EVM is lower and much reasonable than the classic one-stage EVM.


## 5. Conclusion

In this paper, a two-stage SCLAP in the context of emergency logistics is considered. In this problem, the number and capacities of supply centers are uncertain and must be determined. This condition is practical in emergency logistics compared to general logistics, which assumes that the number and capacities of supply centers are explicit and fixed. To solve this problem, a two-stage EVM and a generalized cost function are proposed. However, we convert the initial model into a deterministic model based on the stochastic theory and the uncertain theory proposed by Liu [5, 6] given that the traditional stochastic programing with a simulation method is time-consuming and unstable. To solve this model, we develop an IPSO algorithm and test it based on well-known 20-customer data. In this IPSO algorithm, three improvement strategies are introduced, namely, the Gaussian cloud operator, the Restart strategy, and the adaptive parameter strategy. In addition, we adopt the interior point method instead of the simplex method [7] to solve the second stage of programing. The comparison tests prove that the methods can improve the precision and convergence rates dramatically. The tests reveal the difference between the classic one-stage EVM and the two-stage EVM proposed in this paper and verify that the latter is much reasonable for emergency logistics. We also perform an analysis of the parameters in the proposed generalized cost function, drawing a conclusion that the fixed cost value has an important influence on the decision on the best number of supply centers.However, this study has the following limitations. (1) The conversion of the deterministic model is only possible when the demand of the customer follows a normal distribution. (2) In this paper, only the EVM is converted and studied. According to different optimization problems, stochastic programing can also be constructed as other models, such as the CCP model and the DCP model. However, only the CCP model can be converted into a deterministic model. (3) In actual emergency logistics, the transportation cost between two nodes is probably uncertain and can also be assumed as a stochastic variable. However, it is taken as a deterministic variable in this paper and can be a significant research direction in the future.In addition, several interesting subjects should be considered in future research. For example, a new SCLAP that builds new centers and opens old centers must be considered. The objective must be to choose the best locations to build new centers and find the most appropriate old centers to close simultaneously. Other models with multiple objectives, complex constraints, or fuzzy variables can also be considered. The uncertain CLAP can also be solved with robust optimization, which is becoming a popular new research topic.



---

*Source: 6710929-2017-05-31.xml*