Gathered Village Location Optimization for Chinese Sustainable Urbanization Using an Integrated MODM Approach under Bi-Uncertain Environment

Abstract

Urbanization has become a main challenge all over developing countries in the 21st Century. However, decision making should take into account the different national situations with their complex factors to achieve sustainable development. As standards of living have risen in urban areas, local/neighbor urbanization has become a coming trend in China. With this in mind, the paper focuses on the optimization of nearby gathered village locations in Population Migration (PM) with consideration of both qualitative and quantitative criteria. Therefore, an integrated multiple objective decision making approach (MODM) under a bi-uncertain environment is proposed to solve this problem, which is based on the comprehensive Economy-Society-Ecology-Resource-Religion (ESERR) urbanization concept. The first step is to establish a bi-uncertain multiple objective programming model orienting the problem. Secondly, the model process is composed of fuzzy random variable transformation and the expected value model based on a new fuzzy measure, which is given accordingly to obtain the equivalent model. Thirdly, in order to describe the model efficiently, the Multi-Objective Adaptive Global Local Neighbor Particle Swarm Optimization (MOAGLNPSO) with three-dimensional Pareto optimal judgment criteria is designed. Finally, a case study is tested to validate the effectiveness and to illustrate the advantages of the whole approach. This novel approach can help optimize sustainable urbanization strategies and ensure their realistic application.

1. Introduction

Urbanization has been brought into focus in most developing countries. However, China has its own characteristic urbanization road, since its history of development is different from other developing countries and developed countries. Tan et al. [1] say that much of the rural population continues to migrate to urban areas during the process of China’s urbanization. According to the Department of Economic and Social Affairs of the United Nations [2], more and more Chinese will live in cities in 2050. This change in population distribution resulting from the movement of people from rural to urban areas matches the definition of urbanization [3]. At present, urbanization is one of the most prominent trends. As a result, Population Migration (PM) is one of the most common methods of sustainable urbanization in today’s China.

Lang et al. [4] state that China’s urbanization rate has grown year-by-year, quickly growing the urban economy and bringing social change. However, a new threat for China’s sustainable urbanization progress has emerged: China’s urban living standards are getting higher and higher. Along with the new situation, local/neighbor urbanization will be the next trend in China. Lin and Meulder [5] demonstrate that nearby PM has a positive impact on rapid urbanization by supplying adequate housing, public education services and job opportunities for migrants. Besides, Xie et al. [6] demonstrate that PM will cause the loss of labor force, leading to a large number of abandoned farmland. Moreover, Gao et al. [7] hold that rural housing abandonment will be happening along with PM. Thus, PM must consider the intricate local situation, such as diverse economy, society, ecology, resource and religion conditions in an uncertain environment. Furthermore, China is a country with a complex climatic condition, fragile ecological environment and frequent disasters. For example: (1) strong dust attacked 13 provinces in northern China in 2006; (2) the Wenchuan Earthquake caused great damage in 2008; (3) the Maoxian landslide caused serious losses in 2017; and (4) dense smog appeared in many places across the country recently. Those disasters caused tremendous economic losses and social chaos. Consequently, PM will not only narrow the income gap, but will also avoid those serious disaster security risks. Hence, the process of population growth and migration is a complex and systematic project that involves multiple areas, various participants and conflicting objectives. Under such challenges, scientifically-gathered area/village location is the precondition for solving the problem.

Many studies have investigated the Village location Problem (VLP) in the past few years (e.g., [8,9,10,11,12,13,14,15,16,17]). Research in VLP has referred to the orientation, theory, methodologies, etc., of this domain. Research topics such as economic impact [8], social relationships [9], energy-saving trends [10] and comprehensive aims [11], mathematic analysis [12] and solution approaches [13,14] and the utilization of Geographic Information System (GIS) [15,16] and Information System (IS) [17] have often been the study focus in VLP. However, little scholarly work has been dedicated to the multiple development goals in disaster risk. As there is higher demand for urbanization, traditional VLP has shifted its focus to a comprehensive goal under the disaster risk threat. Thus, it is very important to obtain scientifically-gathered village location information.

Local/neighbor urbanization forces PM progress to face the new challenges. Furthermore, the demand for coordinating conflicting-objectives and reducing cost with avoiding geological disaster risks in PM has forced researchers to focus on more effective gathered village location information. This study investigates how these challenges can be overcome by an integrated Multiple Objective Decision Making (MODM) approach under a bi-uncertain (i.e., recombination of two types of uncertainties: fuzzy and random) environment to optimize the gathered village location problem. The main objectives the present study endeavors to achieve are:

• To explore the optimal strategies of the nearby gathered village location in PM progress for local/neighbor urbanization with Chinese characteristics;
• To make a comprehensive balance of Economy-Society-Ecology-Resource-Religion (ESERR) aspects under a bi-uncertain environment in China’s urbanization;
• To take into full consideration climate-induced geological disaster risk, which might trigger migration;
• To look to the future to improve the key study from the hukou system, pollution and industrialization influences.

In particular, the integrated MODM approach is composed of a multiple objective programming model, a bi-uncertain parameter transformation process and a Multi-Objective Adaptive Global Local Neighbor Particle Swarm Optimization (MOAGLNPSO). The model is established for local governments to pursue a comprehensive balance aim. The decision maker’s first objective is the distance from the gathered village locations to the urban area. Afterwards, the moving resettlement cost is another necessary objective. Furthermore, the integrated urbanization level (i.e., involving economic, social and ecology) is the most important objective for local governments to determine. Meanwhile, the decision must be satisfied with the security constraint (i.e., to avoid locating in a high-risk geological disaster area), development constraint (i.e., to avoid locating in the restricted development area for eco-environmental protection) and logical constraint (i.e., to avoid a negative variable-value). Considering the bi-uncertain rebuilding cost and climate-induced dangerousness weight, the transformation process has a theorem [18], and this theorem could transform fuzzy random variables into trapezoidal fuzzy numbers first. However, the decision maker will have different optimistic-pessimistic attitudes in the realistic uncertain decision making process [19]. To avoid extreme attitudes, it is necessary to utilize a more flexible measure (i.e., measure $Me$) to establish the fuzzy Expected Value Model (EVM), which can finally calculate the fuzzy random variables. To solve the proposed model efficiently, the MOAGLNPSO is a combination of a Pareto Archived Evolution Strategy (PAES) [20] and an Adaptive Global Local Neighbor Particle Swarm Optimization (AGLNPSO), which is developed by incorporating an Adaptive Particle Swarm Optimization (APSO) [21] with a Global Local Neighbor Particle Swarm Optimization (GNLPSO) [22] and a Multi-Objective Particle Swarm Optimization (MOPSO) [23]. In order to simplify the solution verification and describe the Pareto optimal front more intuitively, a judgment criterion is given for the more complex three-dimensional (3D) Pareto optimal solution.

This study explores VLP for PM progress in sustainable urbanization through comprehensive consideration of the ESERR goal with climate-induced geological disaster risk. The proposed model could enrich the orientation of traditional VIP and the existing methodologies. To ensure the convenience in practical applications, the crisp equivalence model and computer algorithm can be used to help researchers and practitioner get a scientific and effective conclusion.

The remainder of this paper is as follows. Firstly, the background of the study is given in Section 2. Secondly, the multiple objective programming model is established in Section 3. Thirdly, the bi-uncertain parameter transformation is processed in Section 4. Furthermore, the MOAGLNPSO is developed in Section 5. In addition, a case study is presented in Section 5. Last, but not least, the advantages, limitations and possible future extensions of this work are referred to in Section 6.

2. Background of the Study

With the high living density in China’s urban areas, local/neighbor urbanization is the most effective approach for sustainable development. This is a way for the rural population not to migrate to large or medium-sized cities, but to nearby small cities and towns. The lifestyle of the migrant has been improved with the productivity improvement, income growth and life quality promotion. Under these circumstances, confirming the gathered area locations of PM seems much more important.

2.1. The Related Literature

Decision making for the gathered area locations of PM in sustainable urbanization involves bi-uncertainty factors. In the past several years, numerous efforts have been made to promote the development of the related research issues. The literature we have studied can be divided into four aspects: (1) sustainable urbanization; (2) VLP; (3) MODM methods; (4) uncertainty and solution algorithm.

2.1.1. Sustainable Urbanization

According to the statement in [1], China’s urbanization road has been unique because of its own national situation. Thus, it can be concluded that population migration will deserve more attention in sub-urbanization or the new style local/neighbor urbanization through the study. Wegren [24] discovers that the fast-growing industrialization, irrational site selection and transformation of land function have impeded rural sustainable development. In addition, Lang et al. [4] point out that disordered development has been hindering China’s urbanization for a long time. Consequently, they study how to form urban and rural communities that meet people’s social demand for sustainable urbanization. Besides, many scholars hold that environmental management is the most influential factor in sustainable development [25,26,27,28,29,30]. Therefore, it is necessary to consider the social, economic, environmental and other sustainable factors in the urbanization process.

2.1.2. VLP

In order to realize the sustainable goal of local population migration, the primary task is to determine the reasonable gathered area/village locations. Accordingly, many scholars have been studying VLP for a long time. Wegren et al. [8] focus on economic impact in VLP. Shi et al. [9] work out that the sustainable urbanization mode, supporting policy, national standards and assessment tools, and sustainable planning are the most important factors in VLP, while Prinsloo et al. [10] propose that VLP should follow the popular energy-saving trend. Liu et al. [11] indicate that de-industrialization will cause brownfield redevelopment, so VLP should consider the environmental, economic and social factors. In problem solution methods, Tang et al. [12] present a hierarchical simulation model to investigate complex rural settlement. Further, Liu et al. [13] analyze the major factors that influence the villages’ locations and emphasize the need to solve the problem dynamically. Some other approaches such as GIS, IS and related methods, have been discussed in VLP [15,16,17]. In particular, Trivedi and Singh [14] address the selecting of location process with multiple objectives, which is similar to the problem of gathered village locations.

2.1.3. MODM Methods

In reality, the gathered village locations for local population migration face various factors, such as distance, cost, risk, and so on. Therefore, MODM is an effective framework for the decision maker to evaluate location rankings. Until now, many methods have been utilized to solve the MODM problems [31,32,33,34]. For example, Zhang et al. [31] study the city sustainability evaluation problem by using objective weights approach; Govindan and Sivakumar [32] propose a fuzzy Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) method to solve the multi-objective linear programming method in green and low-carbon development; Gutjahr and Pichler [33] use non-scalarizing methods to optimize stochastic multi-objective decision making; Delgoda et al. [34] design a novel generic optimization method for irrigation scheduling under multiple objectives; Deng et al. [35] propose an improved APSO algorithm to solve the multi-objective optimization model. In particular, MODM [36] and multi-criteria [37] have been used to assess people’s life and environment for sustainable urbanization.

2.1.4. Uncertainty and Solution Algorithm

Since the decision making process usually involves with many natural and artificial uncertainties, it should consider using the MODM combined with uncertain theory to solve VLP. Furthermore, Shapiro [38] hold that randomness and fuzziness are complementary. Uno et al. [39] use fuzzy random programming to solve facility location problems. Therefore, this kind of combination has similar advantages in VLP. To study the transforming process for fuzzy random variables [18], many scholars also use various computational algorithms to process fuzzy random numbers and obtain problem solutions. Zhong et al. [40] utilize the genetic algorithm to calculate fuzzy random programming models. Similarly, Wang et al. [41] take advantage of the Particle Swarm Optimization algorithm (PSO) to search for approximate optima.

2.2. Problem Description

Based on the perspective of local/neighborhood urbanization, the gathered area should be located nearby towns, which will be conducive to the development of economic production and social life for a long time. As preliminary work, Wang and Gan [42] have already established a village evaluation indicator system based on the fuzzy Data Envelopment Analysis (DEA) method. Furthermore, this research has been applied in the Northwest Sichuan Tibetan Region (NSTR), which is under the threat of serious geological disasters due to the complex local situations. Moreover, the central village is in rapid development, while the weak villages are in recession. According to the paper, the selected central village is close to the town. Consequently, the gathered area is preferred to be located around the central village and gradually developed into a village.

Urbanization is a complex system affected by various factors [43]: for example, labor transfer [44,45], land use/cover [46,47], detrimental impact on ecology and environment [48,49]. Liu et al. [50] have summarized that urbanization is made up of many intertwined and interrelated aspects, such as regional economies, social life, ecological environment, and so on. Accordingly, VLP for PM must coordinate conflicting objectives, which are already referred to above. In addition, religious belief and resource scarcity should also be taken into consideration. Therefore, this study focuses on the comprehensive ESERR goal in sustainable urbanization.

Appropriate solutions in PM will make people have their own stable habitation and job. To solve PM, it is necessary to rebuild houses. On the one hand, the rebuilding price usually combines the depreciation level with the construction-installation expense in the finance subsidy of governments. However, the construction-installation expense itself is various for different house types. On the other hand, the depreciation level is influenced deeply by human subjective judgment. For example, the typical house falls into three categories, and they have discrete probabilities. In the meantime, the depreciation level is vague and uncertain for different house types. Thus, the rebuilding price needs to be qualified with fuzziness and randomness. Hence, due to the complexity of rebuilding price, it is subject to bi-uncertainty with fuzzy random variables.

A similar nature also belongs to the climate-induced dangerousness weight. For example, debris flow is under the influence of stochastic rainfall with fuzzy expression by people. In the NSTR case, this study considers debris flow as the most frequent and serious geological disaster. The type of precipitation amount is divided into six categories, and there is a set of discrete probabilities associated with them. Further, the possible result description of dangerousness level for the different precipitation amounts is vague and uncertain. Hence, the climate-induced dangerousness weight is subject to bi-uncertainty with both fuzziness and randomness.

2.3. Methods Description

In order to solve the aforesaid problem, a novel integrated MODM approach is developed to accomplish the aims, which is combined in the programming model under a bi-uncertain environment and parameter transforming process MOAGLNPSO in the paper.

In reality, people always need to balance the conflictive objectives according to the evaluation criteria. Thus, MODM is a suitable method that can help choose a satisfying solution from several feasible schemes and then achieve the multiple evaluation targets [51]. Based on the characteristics of the proposed problem, the MODM programming model is superior to other models and provides some advantages in this study: (1) representing the complex and conflictive objectives of economy, society and ecology for the target system; (2) ensuring that the decision results conform to the security, development and logical constraints; (3) combining the model with a bi-uncertain environment conveniently.

On the other hand, the model needs to be further processed and then transformed into a solvable model with mathematical meaning as it has the bi-uncertain (i.e., fuzzy random) parameter, which reflects the true reality in practice. Through transforming fuzzy random variables into trapezoidal fuzzy numbers and finally establishing fuzzy EVM, the transforming process can be realized. In addition, MOAGLNPSO has been adopted to solve the multiple objective optimization problem, which is combined wit the Pareto optimal strategy, the adaptive mechanism and a wider range of search space.

2.4. Proposed Research Framework

Above all, the framework of gathered village location optimization for China’s sustainable urbanization can be expressed as in Figure 1 including three research phases given in the following:

Phase I: To describe the problem effectively through the multiple objective gathered village location program model with the fuzzy random rebuilding price and climate-induced dangerousness.

Phase II: To obtain the crisp equivalent model in the bi-uncertain parameter transforming process.

Phase III: To work out the optimal gathered village locations by using MOAGLNPSO.

3. Modeling

Research Phase I: This paper utilizes a bi-uncertain multiple objective programming model for optimal gathered village locations based on the NSTR case via three steps. Problem analysis and model establishment of other areas in China can be obtained according to this case similarly. The following notations are used.

3.1. Notations

Indices:

e	East longitude;
n	North longitude;
r	Religious pilgrimage area, r = 1,2, …,R;
w	Weak village, w = 1,2, …,W;
h	Hospital, h = 1,2, …,H;
s	School, s = 1,2, …,S;
$da$	Disaster high-risk area, da = 1,2, …,DA;
m	Micro element, m = 1,2, …,M;
Variables
$F_1:D$	Distance: consideration of central village and religious pilgrimage area;
$F_2:O$	Cost: consideration of reconstruction and relocation compensation;
$F_3:U$	Urbanization level: combined consideration of economic, society and ecology;
$PS$	Subsidy of price gap;
$PL$	Land price;
$PT$	Transport price;
$AS$	Area of structure;
$AC$	Area of covered;
$F{H}^w$	Farmer household in weak village;
$GD{P}^p$	Per capita GDP;
$R{S}_m^{da}$	Circle radius of micro element for disaster high-risk area;
$\left(R{P}_e^r,R{P}_n^r\right)$	Religious pilgrimage area;
$\left(V_e^w,V_n^w\right)$	Weak village;
$\left(F_e^h,F_n^h\right)$	Hospital ;
$\left(F_e^s,F_n^s\right)$	School;
$\left(C_{me}^{da},C_{mn}^{da}\right)$	Circle center position of micro element for disaster high-risk area;
$\left(V_e^c,V_n^c\right)$	Central village ;
$\left(F{A}_e,F{A}_n\right)$	Farming area;
$\left(T{W}_e,T{W}_n\right)$	Town;
$\left(R{D}_e,R{D}_n\right)$	Restricted development area;
Decision variables
$x_e$,$x_n$	Gathered village;
Bi-uncertain parameter
$\tilde{\overline{\varsigma }}$	Rebuilding price, fuzzy random variable;
$\tilde{\overline{\zeta }}$	Climate-induced dangerousness weight, fuzzy random variable;
Weight
${\alpha }_1$,${\alpha }_2$	Weight for distance;
${\beta }_1$,${\beta }_2$,${\beta }_3$	Weight for urbanization;
${\rho }_h$	Weight for hospital;
${\rho }_s$	Weight for school;
${\varpi }^{da}$	Geography factor dangerousness weight.

3.2. Gathered Village Location Model

Aimed at the balance of society, economy and ecology, the local governments need to minimize the distance/cost and maximize the combination urbanization level through the optimal gathered village locations firstly. Meanwhile, they also need to keep the gathered village locations away from the high-risk and restricted development area secondly. Finally, the whole programming model needs to be proposed.

Step 1. Objectives:

Urbanization is the process of the rural population gathering around the urban areas. Hence, it is necessary to think about a shorter distance from the gathered village to the existing central village around the town. In addition, it is also reasonable to consider the adjacent religious pilgrimage area. However, various religions have diverse culture, such as Buddhism, Daoist, Confucianist, Muslims, etc. Each religion has its own respected creeds, eating habits, ancient laws, living regulations and so on; thus it cannot treat all groups as a whole, or it will dissatisfy its believers. In NSTR, the most popular religious belief is Tibetan Buddhism. To avoid making a mistake about the local culture, this paper only considers the influences of Tibetan Buddhism instead of other religions. Besides, the state has invested much manpower, many material resources and many financial resources to help protect and carry forward the fine traditional Tibetan culture. Therefore, being close to the religious pilgrimage area cannot be ignored in PM.

$$\begin{array}{cc}min{F}_1:D=\hfill & {\alpha }_1\sqrt{{\left(x_e-V_e^c\right)}^2+{\left(x_n-V_n^c\right)}^2}+\hfill \\ & {\alpha }_2{\displaystyle \sum _{r=1}^R}\frac{\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}{{\displaystyle \sum _{r=1}^R}\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}\sqrt{{\left(x_e-R{P}_e^r\right)}^2+{\left(x_n-R{P}_n^r\right)}^2}\hfill \end{array}$$

(1)

where $\frac{\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}{{\displaystyle \sum _{r=1}^R}\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}$ is the weight for the distance of distinguishing the religious pilgrimage area, which is based on the consideration of the distance from the town.

In the population gathering process, a low moving resettlement cost is an important factor. The total cost is consist with rebuilding cost, land cost and transportation cost. $\tilde{\overline{\varsigma }}$ is a bi-uncertain variable combined with the construction-installation expense and depreciation level.

$$\begin{array}{cc}min{F}_2:O=\hfill & \left(\tilde{\overline{\varsigma }}+PS\right)AS+PL·AC+PT·F{H}^w{\displaystyle \sum _{w=1}^W}\sqrt{{\left(x_e-V_e^w\right)}^2+{\left(x_n-V_n^w\right)}^2}\hfill \end{array}$$

(2)

In order to realize the goal of sustainable urbanization, it is necessary to pursue an aim combining the economy, society and ecology [50]. First, the $GD{P}_p$ in each distance reflects the economy aim. Secondly, the distance of the gathered village location to the infrastructures reflects the society aim. Finally, the distance of the gathered village location to the farming area is based on ecology protection. Notably, the hukou system plays a significant role in the process of urbanization in China. However, according to the field research, the local data from the Public Security Bureau, Civil Affairs Bureau and Village Committees show that this system has little impact on local migration. Therefore, the third objective does not consider the influence of the hukou system.

$$\begin{array}{cc}max{F}_3:U=\hfill & {\beta }_1\frac{GD{P}^p}{\sqrt{{\left(x_e-C{Y}_e\right)}^2+{\left(x_n-C{Y}_n\right)}^2}}+{\beta }_2\left[{\displaystyle \sum _{h=1}^H}\frac{{\rho }_h}{\sqrt{{\left(x_e-F_e^h\right)}^2+{\left(x_n-F_n^h\right)}^2}}+\right.\hfill \\ & \left.{\displaystyle \sum _{s=1}^S}\frac{{\rho }_s}{\sqrt{{\left(x_e-F_e^s\right)}^2+{\left(x_n-F_n^s\right)}^2}}\right]+{\beta }_3\frac{1}{\sqrt{{\left(x_e-F{A}_e\right)}^2+{\left(x_n-F{A}_n\right)}^2}}\hfill \end{array}$$

(3)

Step 2. Constraints:

Security constraint: The gathered village should not be located in the geological disaster high-risk area for personal and property security.

The occurrence of geological disaster is affected by various factors. In particular, it is greatly impacted by the weather. In order to express the security constraint effectively, this study tries to: (1) describe the fuzzy random climate-induced dangerousness weight; (2) consider the climate, topography and geomorphology, geological structure and lithology type (except climate, others are called by a joint name: geography factor) to get comprehensive dangerousness weights in different disaster risk districts; (3) process weights’ relative comparison to obtain the larger/smaller proportion in different disaster risk districts; (4) use the micro element method to develop the location constraint. The details are shown as below.

$${\left(x_e-C_{me}^{da}\right)}^2+{\left(x_n-C_{mn}^{da}\right)}^2\ge \tilde{\overline{\zeta }}{\varpi }^{da}{\left(R{S}_m^{da}\right)}^2;\forall da\in DA;\forall m\in M$$

(4)

Development constraint: the gathered village locations should keep away from the restricted development area in order to preserve the ecological environment.

$$R{D}_e^{min}\le x_e\le R{D}_e^{max};R{D}_n^{min}\le x_n\le R{D}_n^{max}$$

(5)

Logical constraint: the longitude coordinate of the gathered village locations should not be negative.

$$x_e\in R^{+};x_n\in R^{+}$$

(6)

Step 3. Programming model:

In summary, the bi-uncertain gathered village location programming model, which has multiple objectives, is as below.

$$\begin{array}{cc}min{F}_1:D=\hfill & {\alpha }_1{\displaystyle {\displaystyle \sum _{r=1}^R}}\frac{\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}{{\displaystyle \sum _{r=1}^R}\sqrt{{\left(R{P}_e^r-C{Y}_e\right)}^2+{\left(R{P}_n^r-C{Y}_n\right)}^2}}\sqrt{{\left(x_e-R{P}_e^r\right)}^2+{\left(x_n-R{P}_n^r\right)}^2}+\hfill \\ & {\alpha }_2\sqrt{{\left(x_e-V_e^c\right)}^2+{\left(x_n-V_n^c\right)}^2}\hfill \\ min{F}_2:O=\hfill & \left(\tilde{\overline{\varsigma }}+PS\right)AS+PL·AC+PT·F{H}^w{\displaystyle \sum _{w=1}^W}\sqrt{{\left(x_e-V_e^w\right)}^2+{\left(x_n-V_n^w\right)}^2}\hfill \\ max{F}_3:U=\hfill & {\beta }_1\frac{GD{P}^p}{\sqrt{{\left(x_e-C{Y}_e\right)}^2+{\left(x_n-C{Y}_n\right)}^2}}+{\beta }_2\left[{\displaystyle \sum _{h=1}^H}\frac{{\rho }_h}{\sqrt{{\left(x_e-F_e^h\right)}^2+{\left(x_n-F_n^h\right)}^2}}+\right.\hfill \\ & \left.{\displaystyle \sum _{s=1}^S}\frac{{\rho }_s}{\sqrt{{\left(x_e-F_e^s\right)}^2+{\left(x_n-F_n^s\right)}^2}}\right]+{\beta }_3\frac{1}{\sqrt{{\left(x_e-F{A}_e\right)}^2+{\left(x_n-F{A}_n\right)}^2}}\hfill \\ s.t.\hfill & \left\{\begin{array}{c}{\left(x_e-C_{me}^{da}\right)}^2+{\left(x_n-C_{mn}^{da}\right)}^2\ge \tilde{\overline{\zeta }}{\varpi }^{da}{\left(R{S}_m^{da}\right)}^2;\forall da\in DA;\forall m\in M\hfill \\ R{D}_e^{min}\le x_e\le R{D}_e^{max}\hfill \\ R{D}_n^{min}\le x_n\le R{D}_n^{max}\hfill \\ x_e\in R^{+}\hfill \\ x_n\in R^{+}\hfill \end{array}\right.\hfill \end{array}$$

(7)

4. Model Process

Research Phase II: As the model has bi-uncertain parameters (i.e., fuzzy random variables: rebuilding price $\tilde{\overline{\varsigma }}$ and climate-induced dangerousness weight $\tilde{\overline{\zeta }}$), it requires further treatment and needs to be transformed into a solvable model with mathematical meaning. In order to describe the transforming process in detail, some basic knowledge is stated below.

Definition 1.

Let there be a domain U. Let $\tilde{A}$ be a fuzzy set, which is defined on U. If α is the possibility level and $0\le \alpha \le 1$, ${\tilde{A}}_{\alpha }$ consists of all elements whose degrees of membership in $\tilde{A}$ are greater than or equal to α.

$${\tilde{A}}_{\alpha }=\{x\in U|{\mu }_{\tilde{A}}\left(x\right)\ge \alpha \}$$

(8)

then ${\tilde{A}}_{\alpha }$ is called the α-level set of fuzzy set $\tilde{A}$.

Definition 2.

Let ε be a discrete random variable defined on a probability space $(\Omega ,\mathcal{A},Pr)$ with the discrete distribution $P_{\epsilon }\left(x\right)=P\{x=x_n\},n=1,2,\dots$, and θ be any given probability level and $0\le \theta \le max{P}_{\epsilon }\left(x\right)$. ${\epsilon }_{\theta }$ consists of all elements whose values of $P_{\epsilon }\left(x\right)$ for ε are greater than or equal to θ.

$${\epsilon }_{\theta }=\{x\in \mathbb{R}|P_{\epsilon }\left(x\right)\ge \theta \}$$

(9)

then ${\epsilon }_{\theta }$ is called the θ-level set of random variable ε.

Definition 3.

Let $(\Theta ,Pos)$ be a possibility space and A be a set in $P(\Theta )$. Then, the fuzzy measure of A is:

$$Me\left\{A\right\}=Nec\left\{A\right\}+\lambda \left(Pos\right\{A\}-Nec\{Ag\left\}\right)$$

(10)

where $\lambda (0\le \lambda \le 1)$ is the optimistic-pessimistic parameter to determine the combined attitude of a DM. $Pos\left\{A\right\}$ and $Nec\left\{A\right\}$ are the possibility and necessity fuzzy measures proposed by Dubois [52].

Theorem 1.

Let $\tilde{\overline{\xi }}=\left\{\begin{array}{cc}\left(a_{1L},a_{1C},a_{1R}\right)\hfill & withprobability{p}_1\hfill \\ ⋮\hfill & ⋮\hfill \\ \left(a_{iL},a_{iC},a_{iR}\right)\hfill & withprobability{p}_i\hfill \\ ⋮\hfill & ⋮\hfill \\ \left(a_{IL},a_{IC},a_{IR}\right)\hfill & withprobability{p}_I\hfill \end{array}\right.$ be a fuzzy random variable, which has a discrete random distribution with fluctuating lower, central and upper parameters for the fuzzy property. The discrete distribution is $P_{\psi }\left(x\right)$. δ is any given probability level of a random variable, and η is any given possibility level of fuzzy variable, then the fuzzy random variable can be transformed into a $(\delta ,\eta )$-level trapezoidal fuzzy variable $\tilde{\xi }$.

Proposition 1.

Let $\tilde{\xi }=\left(r_1,r_2,r_3,r_4\right)$ be a trapezoidal fuzzy variable. Then, its expected value is:

$$E^{Me}\left[\tilde{\xi }\right]=\left\{\begin{array}{cc}\frac{\lambda }{2}\left(r_1+r_2\right)+\frac{1-\lambda }{2}\left(r_3+r_4\right)\hfill & if{r}_4\le 0\hfill \\ \frac{\lambda }{2}\left(r_1+r_2\right)+\frac{\lambda r_4^2-\left(1-\lambda \right)r_3^2}{2\left(r_3+r_4\right)}\hfill & if{r}_3\le 0\le r_4\hfill \\ \frac{\lambda }{2}\left(r_1+r_2+r_3+r_4\right)\hfill & if{r}_2\le 0\le r_3\hfill \\ \frac{(1-\lambda )r_2^2-\lambda r_1^2}{2\left(r_1+r_2\right)}+\frac{\lambda }{2}\left(r_3+r_4\right)\hfill & if{r}_1\le 0\le r_2\hfill \\ \frac{1-\lambda }{2}\left(r_1+r_2\right)+\frac{\lambda }{2}\left(r_3+r_4\right)\hfill & if0\le r_1\hfill \end{array}\right.$$

(11)

Although there are many properties and transformation approaches for fuzzy random variables, Gan and Xu [18] propose a theorem that could transform fuzzy random variables into fuzzy variables and similar to trapezoidal fuzzy numbers for greater convenience. Meanwhile, Xu and Zhou [19] give a fuzzy measure $Me$, which is suitable for the realistic uncertain decision making process. In this paper, the transforming process involving two steps with Theorem 1 and the fuzzy EVM based on measure $Me$ is used to deal with problems in the bi-uncertain model.

Step 1: Transformation of fuzzy random variable $\tilde{\overline{\xi }}$.

Through Theorem 1, the fuzzy random rebuilding price, viz. $\tilde{\overline{\xi }}$, can be transformed into $(\delta ,\eta )$-level trapezoidal fuzzy variable ${\tilde{\overline{\xi }}}_{(\delta ,\eta )}$ as shown in Figure 2.

According to Definition 2, the $\delta$-level sets (or $\delta$-cuts) of the discrete random variable $\psi$ can be denoted as follows:

$${\psi }_{\delta }=\left[{\psi }_{\delta }^L,{\psi }_{\delta }^R\right]=\{x\in \mathbb{R}|P_{\psi }\left(x\right)\ge \delta \}$$

(12)

Here, ${\psi }_{\delta }^L=min\{x\in \mathbb{R}|P_{\psi }\left(x\right)\ge \delta \}$ and ${\psi }_{\delta }^R=max\{x\in \mathbb{R}|P_{\psi }\left(x\right)\ge \delta \}$. The parameter $\delta \in \left[0,max{P}_{\psi }\left(x\right)\right]$ here reflects the optimism degree for the decision maker. These intervals indicate where the range of the data lies at the probability level $\delta$. Note that ${\psi }_{\delta }$ is a crisp set.

Let $X=\{x_{\omega }=\psi \left(\omega \right)\in \mathbb{R}|P_{\psi }\left(\psi \left(\omega \right)\right)\ge \delta ,\omega \in \Omega \}$; it is not hard to prove that $X=\left[{\psi }_{\delta }^L,{\psi }_{\delta }^R\right]={\psi }_{\delta }$, viz. $minX={\psi }_{\delta }^L$ and $maxX={\psi }_{\delta }^R$. In other words, ${\psi }_{\delta }^L$ is the minimum value that $\psi$ achieves with probability $\delta$; ${\psi }_{\delta }^R$ is the maximum value that $\psi$ achieves with probability $\delta$. Therefore, the $\delta$-level fuzzy random variable ${\tilde{\overline{\xi }}}_{\delta }$ can be defined as ${\tilde{\overline{\xi }}}_{\delta }=\left\{\begin{array}{cc}{\psi }_{\delta }^L=\left(a_{(\delta ,L)}^L,a_{(\delta ,C)}^L,a_{(\delta ,R)}^L\right)& \mathrm{with}\mathrm{probability}{p}_{\delta }^L\\ ⋮& ⋮\\ {\psi }_{\delta }^R=\left(a_{(\delta ,L)}^R,a_{(\delta ,C)}^R,a_{(\delta ,R)}^R\right)& \mathrm{with}\mathrm{probability}{p}_{\delta }^R\end{array}\right.$.

It can also be denoted as follows:

$${\tilde{\overline{\xi }}}_{\delta }=\{{\tilde{\xi }}_{\delta }\left(\omega \right)=(a_{(\delta ,L)}\left(\omega \right),a_{(\delta ,C)}\left(\omega \right),a_{(\delta ,R)}\left(\omega \right))\mathrm{with}\mathrm{probability}p\left(\omega \right)|x_{\omega }\in X,\omega \in \Omega \}$$

(13)

where ${\tilde{\xi }}_{\delta }\left(\omega \right)$ is a fuzzy variable. The variable ${\tilde{\overline{\xi }}}_{\delta }$ can be expressed in another form as ${\tilde{\overline{\xi }}}_{\delta }={\bigcup }_{\omega \in \Omega }{\tilde{\xi }}_{\delta }\left(\omega \right)={\tilde{\xi }}_{\delta }(\Omega )$; here, ${\tilde{\xi }}_{\delta }\left(\omega \right)(\omega \in \Omega )$ are fuzzy variables. Therefore, the fuzzy random variable $\tilde{\overline{\xi }}$ is transformed into a group of fuzzy variables ${\tilde{\xi }}_{\delta }\left(\omega \right)(\omega \in \Omega )$, which is denoted as ${\tilde{\xi }}_{\delta }(\Omega )$. On the basis of the concept of fuzzy variable $\eta$-level sets (or $\eta$-cuts) (see Definition 2), for the parameter $0\le \eta \le 1$, let:

$${\tilde{\xi }}_{(\delta ,\eta )}\left(\omega \right)=\left[{\xi }_{(\delta ,\eta )}^L\left(\omega \right),{\xi }_{(\delta ,\eta )}^R\left(\omega \right)\right]=\{x\in U|{\mu }_{{\tilde{\xi }}_{\delta }\left(\omega \right)}\left(x\right)\ge \eta \}$$

(14)

then the $\eta$-level sets (or $\eta$ cuts) of ${\tilde{\xi }}_{\delta }(\Omega )$ are defined as follows:

$${\tilde{\xi }}_{(\delta ,\eta )}(\Omega )=\left\{{\tilde{\xi }}_{(\delta ,\eta )}\left(\omega \right)=\left[{\xi }_{(\delta ,\eta )}^L\left(\omega \right),{\xi }_{(\delta ,\eta )}^R\left(\omega \right)\right]|\omega \in \Omega \right\}$$

(15)

here, ${\xi }_{(\delta ,\eta )}^L\left(\omega \right)=inf{\mu }_{{\tilde{\xi }}_{\delta }\left(\omega \right)}^{-1}\left(\eta \right)$, ${\xi }_{(\delta ,\eta )}^R\left(\omega \right)=sup{\mu }_{{\tilde{\xi }}_{\delta }\left(\omega \right)}^{-1}\left(\eta \right)$, $\omega \in \Omega$. Inspired by the fuzzy expected value of the fuzzy random variable proposed by [38], it can be obtained as follows:

$$\begin{array}{cc}{a}_{(\delta ,L)}={\displaystyle \sum _{\omega }}p\left(\omega \right)a_{(\delta ,L)}\left(\omega \right);\hfill & a_{(\delta ,R)}={\displaystyle \sum _{\omega }}p\left(\omega \right)a_{(\delta ,R)}\left(\omega \right)\hfill \\ {\xi }_{(\delta ,\eta )}^L={\displaystyle \sum _{\omega }}p\left(\omega \right){\xi }_{(\delta ,\eta )}^L\left(\omega \right);\hfill & {\xi }_{(\delta ,\eta )}^R={\displaystyle \sum _{\omega }}p\left(\omega \right){\xi }_{(\delta ,\eta )}^R\left(\omega \right)\hfill \end{array}$$

(16)

Consequently, $\tilde{\overline{\xi }}$ can be transformed into ${\tilde{\xi }}_{(\delta ,\eta )}$ by the $\delta$-cuts and $\eta$-cuts.

Where $0\le \eta \le 1$ and $\delta \in [0,max{P}_{\psi }\left(x\right)]$, let $a_{(\delta ,L)}={\left[s\right]}_L$, $a_{(\delta ,R)}={\left[s\right]}_R$ ${\xi }_{(\delta ,\eta )}^L={\displaystyle \underline{s}}$ and ${\xi }_{(\delta ,\eta )}^R=\overline{s}$, then the fuzzy random variable $\tilde{\overline{\xi }}$ can be transformed into the $(\delta ,\eta )$-level trapezoidal fuzzy variable ${\tilde{\xi }}_{(\delta ,\eta )}$ by the following equation:

$$\tilde{\overline{\xi }}⟶{\tilde{\xi }}_{(\delta ,\eta )}=\left({\left[s\right]}_L,{\displaystyle \underline{s}},\overline{s},{\left[s\right]}_R\right)$$

(17)

The parameters $\delta$ and $\eta$ reflect the optimism degree of the decision maker. Thus, the fuzzy random variable $\tilde{\overline{\xi }}$ is transformed into a fuzzy variable, which is a trapezoidal fuzzy number with the membership function ${\mu }_{{\tilde{\xi }}_{(\delta ,\eta )}\left(x\right)}$. The value of ${\mu }_{{\tilde{\xi }}_{(\delta ,\eta )}\left(x\right)}$ at $x\in \left[{\left[s\right]}_L,{\left[s\right]}_R\right]$ is considered subjectively to be one as below:

$${\mu }_{{\tilde{\xi }}_{(\delta ,\eta )\left(x\right)}}=\left\{\begin{array}{ccc}1\hfill & \mathrm{if}\hfill & {\displaystyle \underline{s}}\le x<\overline{s}\hfill \\ \frac{x-{\left[s\right]}_L}{{\displaystyle \underline{s}}-{\left[m\right]}_L}\hfill & \mathrm{if}\hfill & {\left[s\right]}_L\le x<{\displaystyle \underline{s}}\hfill \\ \frac{{\left[s\right]}_R-x}{{\left[s\right]}_R-\overline{s}}\hfill & \mathrm{if}\hfill & \overline{s}\le x<{\left[s\right]}_R\hfill \\ 0\hfill & \mathrm{if}\hfill & x<{\left[s\right]}_L,x>{\left[s\right]}_R\hfill \end{array}\right.$$

(18)

Therefore, the fuzzy random objective Equation (2) and constraint Equation (4) can be transformed into Equations (19) and (20) with a fuzzy parameter as follows.

$$\begin{array}{cc}min{F}_2:O=\hfill & ({\tilde{\varsigma }}_{(\delta ,\eta )}+PS)AS+PL·AC+PT·F{H}^w{\displaystyle \sum _{w=1}^W}\sqrt{{(x_e-V_e^w)}^2+{(x_n-V_n^w)}^2}\hfill \end{array}$$

(19)

$${\left(x_e-C_{me}^{da}\right)}^2+{\left(x_n-C_{mn}^{da}\right)}^2\ge {\tilde{\zeta }}_{(\delta ,\eta )}{\varpi }^{da}{\left(R{S}_m^{da}\right)}^2;\forall da\in DA;\forall m\in M$$

(20)

Step 2: EVM based on measure $Me$ of fuzzy variable ${\tilde{\overline{\xi }}}_{(\delta ,\eta )}$.

According to Proposition 1, the expected value of trapezoidal fuzzy variable ${\tilde{\xi }}_{(\delta ,\eta )}=\left({\left[s\right]}_L,{\displaystyle \underline{s}},\overline{s},{\left[s\right]}_R\right)$ is as below.

$$E^{Me}\left[{\tilde{\xi }}_{(\delta ,\eta )}\right]=\frac{1-\lambda }{2}\left({\left[s\right]}_L+{\displaystyle \underline{s}}\right)+\frac{\lambda }{2}\left(\overline{s}+{\left[s\right]}_R\right)$$

(21)

where $\left({\left[s\right]}_L,{\displaystyle \underline{s}},\overline{s},{\left[s\right]}_R\right)\ge 0$.

When $\lambda =0.5$, it is a special case of $Me$. This means that the DM takes a compromise attitude, then:

$$E^{Me}\left[{\tilde{\xi }}_{(\delta ,\eta )}\right]=\frac{\left({\left[s\right]}_L+{\displaystyle \underline{s}}+\overline{s}+{\left[s\right]}_R\right)}{4}$$

(22)

Then, the fuzzy objective Equation (19) and constraint Equation (20) can be transformed into Equations (23) and (24).

$$\begin{array}{cc}min{F}_2:O=\hfill & (E^{Me}\left[{\tilde{\varsigma }}_{(\delta ,\eta )}\right]+PS)AS+PL·AC+PT·F{H}^w{\displaystyle \sum _{w=1}^W}\sqrt{{(x_e-V_e^w)}^2+{(x_n-V_n^w)}^2}\hfill \end{array}$$

(23)

$${\left(x_e-C_{me}^{da}\right)}^2+{\left(x_n-C_{mn}^{da}\right)}^2\ge E^{Me}\left[{\tilde{\zeta }}_{(\delta ,\eta )}\right]{\varpi }^{da}{\left(R{S}_m^{da}\right)}^2;\forall da\in DA;\forall m\in M$$

(24)

where $E^{Me}\left[{\tilde{\varsigma }}_{(\delta ,\eta )}\right]$ and $E^{Me}\left[{\tilde{\zeta }}_{(\delta ,\eta )}\right]$ are the expected values of the random variables.

5. Solution Method

Research Phase III: PSO has been adopted for dealing with multiple objective optimization problems and has been found to be very successful in heuristics [23]. Thus, PSO is adopted in this study based on this consideration. This paper proposes the MOAGLINPSO algorithm, which is made up of PAES [20], AGLNPSO [21,22] and MOPSO [23]. Of course, this proposed algorithm may not be the best. However, it can assist in obtaining an effective solution, which has been demonstrated in the analysis of the case. In the future, in order to get better solutions more effectively, alternative approaches and algorithms (e.g., other exact approaches, (meta-)heuristics, evolutionary algorithms, etc.) will be compared.

5.1. Overall Procedure for the Proposed Algorithm

The flowchart of the proposed algorithm is shown in Figure 3 including seven steps.

Step 1.

Initialize the parameters: swarm_size, iteration_max, the range of velocity and position for the variables, the personal best position acceleration constant, the global best position acceleration constant, the local best position acceleration constant, the near neighbor best acceleration constant, and the inertia weight_max/weight_min. Then, initialize the velocities and positions of the particle-represented solutions.

Step 2.

Check the feasibility and decode the particles.

Step 3.

Calculate the three objectives to evaluate every particle.

Step 4.

Calculate the $pbest$, $gbest$, $lbest$ using the multi-objective method and $nbest$. Restore the non-dominating solutions (i.e., the (global) elite individuals) and objective values.

Step 5.

Update the inertia weight for each iteration.

Step 6.

Update the velocity and position of each particle.

Step 7.

Check the MOAGLNPSO termination: If the stopping criterion (i.e., iteration_max) is met, then end the MOAGLNPSO procedure to obtain the optimal solution. and it terminates. Otherwise, go to Step 2.

The details of the MOAGLNPSO are described as follows, and the notations used are shown:

$s:$	Particle index, $s=1,\dots S$
$\tau :$	Iteration index, $\tau =1,\dots T$
$u_r:$	Uniform random number in the interval $[0,1]$
$w\left(\tau \right):$	Inertia weight in the $\tau$-th iteration
$w^{max}:$	Maximum inertia weight value
$w^{min}:$	Minimum inertia weight value
${\omega }_s\left(\tau \right)$:	Velocity of the s-th particle in the $\tau$-th iteration
${\theta }_s\left(\tau \right):$	Position of the s-th particle in the $\tau$-th iteration
${\psi }_s^p$:	Personal best position of the s-th particle
${\psi }_s^g$:	Global best position of the s-th particle
${\psi }_s^L$:	Local best position of the s-th particle
${\psi }_s^N:$	Near neighbor best position of the s-th particle
$c_p:$	Personal best position acceleration constant
$c_g:$	Global best position acceleration constant
$c_l:$	Local best position acceleration constant
$c_n:$	Near neighbor best position acceleration constant
${\omega }^{max}:$	Maximum velocity value
${\omega }^{min}:$	Minimum velocity value
${\theta }^{max}:$	Maximum position value
${\theta }^{min}:$	Minimum position value
$R_s:$	The s-th set of solutions
$c:$	The current solution randomly selected from the
	non-dominated solutions
$c^N:$	New generated solution

Step 1: Solution representation and particle swarm initialization.

In this paper, the particle-represented solutions are $x_e\mathrm{and}{x}_n$ (i.e., gathered village location), which are the coordinates produced by the east-north longitude.

Initialize S particles as a swarm; generate the s-th particle with random position ${\theta }_s$ in the coordinate range of the considered area scope. Randomly generate velocity for each particle in the range ${\omega }^{min}\le \omega \le {\omega }^{max}$. Set the iteration $\tau =1$. Set swarm_size S, iteration_max T, personal best position acceleration constant $c_p$, global best position acceleration constant $c_g$, local best position acceleration constant $c_l$ and near neighbor best position acceleration constant $c_n$.

Step 2: Feasibility checking and decoding method.

Since the gathered village location should satisfy the security, development and logical constraints, checking and abandoning the infeasible particles are needed. Then, the particle-represented solution can be decoded into a solution as the east-north longitude through a common way for the problem.

Step 3: Particle evaluation.

For $s=1,\dots ,S$, set ${\theta }_s\left(\tau \right)$ into the solution $R_s$, that is $x_e,x_n$ (i.e., gathered village location), and put it into the the optimal objectives $F_1:D,F_2:O,F_3:U$ and calculate them, respectively.

Step 4: Multi-objective method.

Procedure: PAES

generate a new solution $c^N$

if (c dominates $c^N$)

discard $c^N$

else if ($c^N$ dominates c)

replace c with $c^N$ and add $c^N$ to the archive

else if ($c^N$ is dominated by any member of the archive)

discard $c^N$

else if ($c^N$ dominates any member of the archive)

replace it with $c^N$ and add $c^N$ to the archive and discard all other members dominated

by $c^N$

else

apply test procedure to c, $c^N$, determine which to become the new current solution

and whether to add $c^N$ to the archive

until a termination criterion has been reached, return to the beginning

Procedure: test

if the archive is not full

add $c^N$ to the archive

if ($c^N$ is in a less crowded region of the archive than c)

accept $c^N$ as the new current solution

else

maintain c as the current solution

else if ($c^N$ is in a less crowded region of the archive than any other member on the archive)

add $c^N$ to the archive and remove a member of the archive from the most crowded region

if ($c^N$ is in a less crowded region of the archive than c)

accept $c^N$ as the new current solution

else

maintain c as the current solution

else

do not add $c^N$ to the archive

Selection: (1) divide 10 by the number of particles in each hypercube to get its score; (2) apply roulette wheel selection to the hypercube according to their scores and select a hypercube; (3) uniformly choose a member of that hypercube.

The multi-objective method consists of the PAES procedure and the test procedure, and the selection is introduced to calculate $pbest$, $gbest$ and $lbest$. This method uses a truncated archive to store the elite individuals (i.e., non-dominated solutions), which is used to separate the objective function space into hypercubes, each of which has a score based on its density. The best selection is based on a roulette wheel to select the most suitable hypercube first and then uniformly choose a solution. Note that the initialized solution is regarded as the $pbest$ and the non-dominated solution of each particle at the first iteration. When the iteration updates, the updated solution and the non-dominated solutions are used to calculate the $pbest$ by the method. After the $pbest$ has been confirmed at each iteration, the $pbest$ non-dominated solutions for all particles are considered with the $gbest$ non-dominated solutions (i.e., there is no $gbest$ non-dominated solution at initialization) to calculate the $gbest$ by the method. Similar to the $gbest$, among all the $pbest$ non-dominated solutions from K neighbors of the s-th particle and $lbest$ non-dominated solutions, $lbest$ is also set using this method. For each particle, set ${\psi }_s^N={\psi }_o^N$ maximizing $\frac{Z\left({\theta }_s\right)-Z\left({\psi }_o\right)}{{\theta }_s-{\psi }_o}$ to get $nbest$, $o\in S\backslash s$. Here, Z refers to the objective functions, and the difference position value considers the distance. The details for the PAES procedure, test procedure and selection procedure are outlined similarly for the $pbest$, $gbest$, $lbest$ as above, where c is the current solution randomly selected from the non-dominated solutions. Note that c is randomly selected from the $pbest$ non-dominated solutions to calculate the $gbest$, etc., at the first iteration. Therefore, the $gbest$ non-dominated solutions at the T-th particle is the final solutions to the problem.

Step 5: Inertia weight updating.

Update the inertia weight for iteration $\tau$ by using the equations:

$$\begin{array}{c}\overline{\omega }=\frac{{\displaystyle \sum _{s=1}^S}{\displaystyle \sum _{h=1}^H}\left|{\omega }_{sh}\right|}{S·H}\hfill \\ {\omega }^{*}=\left\{\begin{array}{cc}\left(1-\frac{1.8\tau }{T}\right){\omega }^{max},\hfill & 0\le \tau \le T/2\hfill \\ \left(0.2-\frac{0.2\tau }{T}\right){\omega }^{max},\hfill & T/2\le \tau \le T\hfill \end{array}\right.\hfill \\ ▵w=\frac{({\omega }^{*}-\overline{\omega })}{{\omega }^{max}}(w^{max}-w^{min})\hfill \\ w=w+▵w\hfill \\ w=w^{max}ifw>w^{max}\hfill \\ w=w^{min}ifw>w^{min}\hfill \end{array}$$

(25)

Step 6: Velocity and position updating.

Update the velocity and the position of each s-th particle by using the equations:

$$\begin{array}{c}{\omega }_s(\tau +1)=w\left(\tau \right){\omega }_s\left(\tau \right)+c_p{u}_r({\psi }_s-{\theta }_s\left(\tau \right))+c_g{u}_r({\psi }_g-{\theta }_s\left(\tau \right))+c_l{u}_r({\psi }_s^L-{\theta }_s\left(\tau \right))\hfill \\ {\theta }_s(\tau +1)={\theta }_s\left(\tau \right)+{\omega }_s(\tau +1)\hfill \\ \mathrm{If}{\theta }_s(\tau +1)>{\theta }^{max},\mathrm{then}\mathrm{set}{\theta }_s(\tau +1)={\theta }^{max}{\omega }_s(\tau +1)=0\hfill \\ \mathrm{If}{\theta }_s(\tau +1)<{\theta }^{min},\mathrm{then}\mathrm{set}{\theta }_s(\tau +1)={\theta }^{min}{\omega }_s(\tau +1)=0\hfill \end{array}$$

(26)

Step 7: Check the algorithm termination.

If the stopping criterion is met (i.e., iteration_max), end the MOAGLNPSO procedure to obtain the optimal solution and terminate it. Otherwise, the algorithm needs to continue.

5.2. Non-Dominating Solution Evaluation

In order to evaluate the quality of the non-dominating solution set, this paper gives four suitable indicators based on the study of Zitzler [53].

(a) The average distance ${\Phi }_1(\Theta )$ gives the average distance to the non-dominating solution set.

$${\Phi }_1(\Theta ):=\frac{1}{|\Theta |=I}{\displaystyle \sum _{i=1}^I}\sqrt{{(d_i-\overline{d})}^2}$$

(27)

where $\overline{d}:=\frac{1}{|\Theta |}{\displaystyle \sum _{i=1}^I}{d}_i$, $d_i=min\{{\displaystyle \sum _{q=1}^Q}\sqrt{{(F_q\left({\theta }_i\right)-F_q\left({\theta }_j\right))}^2};{\theta }_i,{\theta }_j\in \Theta ,i\ne j\}$ and $|·|$ denote the number of the set’s elements.

(b) The distribution ${\Phi }_2(\Theta )$ takes the distribution in combination with the non-dominating solution set.

$${\Phi }_2(\Theta ):=\frac{1}{C_I^2}{\displaystyle \sum _{\theta \in \Theta }}\left|\right\{\Vert {\theta }_i-{\theta }_j\Vert <\sigma ;{\theta }_i,{\theta }_j\in \Theta ,i\ne j\}$$

(28)

where $\Vert ·\Vert$ denotes the elements’ distance.

(c) The extent ${\Phi }_3(\Theta )$ considers the extent of the Pareto optimal front.

$${\Phi }_3(\Theta ):=\sqrt{{\displaystyle \sum _{q=1}^Q}{(F_q\left({\theta }_i\right)-F_q\left({\theta }_j\right))}^2};{\theta }_i,{\theta }_j\in \Theta ,i\ne j$$

(29)

(d) The set convergence $\varpi$ reflects the stabilization of the non-dominating solution.

$$\begin{array}{c}m\in M,n\in N,\chi =0.\hfill \\ \mathrm{If}\mathrm{traversal}M\mathrm{for}\mathrm{any}m\mathrm{there}\mathrm{is}n=m,n\in N\mathrm{then}\chi =\chi +1\hfill \\ \varpi =\frac{\chi }{\left|M\right|}\hfill \end{array}$$

(30)

That is to say, if $\varpi \ge ϵ$, the non-dominating solution set is stable and the approximation termination is achieved.

6. A Case Study

In this section, computational experiments were carried out in Batang town, which is located in the west of Ganzi. As Figure 4 shows, the location is in the hill plateau mountainous area of NSTR. The MODM approach is validated, and the efficiency of the algorithm is tested through the illustrative example on the dataset adopted from the case.

6.1. Presentation of the Case Problem

Through the study of Wang and Gan [42], Xiasangka has been selected as the central village with its own development advantage. In contrast, the other seven villages (i.e., Renai, Yudi, etc.) have been regarded as the weak villages, which need out-migration. That is to say that, the rural population would migrate from these weak villages to the gathered village around Xiasangka as shown below in Figure 5A. Besides, there are seven hospitals and 10 schools in Batang town. The religious pilgrimage area, farming area, climate-induced high-risk geological disaster area (divided into four districts according to the dangerousness level) and restricted development area locations in the town are described in the following Figure 5B.

The fuzzy random rebuilding price and climate-induced dangerousness weight are obtained through the Analytic Hierarchy Process (AHP) [54] method by an expert based on the investigation and historical data from Batang town, which is given in the following.

$$\tilde{\overline{\varsigma }}=\left\{\begin{array}{cc}(490,580,667)\hfill & \mathrm{with}\mathrm{probability}0.56\hfill \\ (510,600,720)\hfill & \mathrm{with}\mathrm{probability}0.17\hfill \\ (500,593,700)\hfill & \mathrm{with}\mathrm{probability}0.27\hfill \end{array}\right.$$

The calculation result of expected value is 589.61.

$$\tilde{\overline{\zeta }}=\left\{\begin{array}{cc}(0.1,0.15,0.2)\hfill & \mathrm{with}\mathrm{probability}0.07\hfill \\ (0.1,0.2,0.3)\hfill & \mathrm{with}\mathrm{probability}0.11\hfill \\ (0.3,0.35,0.4)\hfill & \mathrm{with}\mathrm{probability}0.15\hfill \\ (0.5,0.55,0.6)\hfill & \mathrm{with}\mathrm{probability}0.20\hfill \\ (0.6,0.7,0.8)\hfill & \mathrm{with}\mathrm{probability}0.21\hfill \\ (0.8,0.85,0.9)\hfill & \mathrm{with}\mathrm{probability}0.26\hfill \end{array}\right.$$

The calculation result of expected value is 0.56.

In a similar way, geography factor dangerousness weights for the four high-risk geological disaster districts can be obtained from

Table 1. The dangerousness weights for geological disaster high-risk districts.

Weight Type	Deda Country District	Lieyi Country District	Batang Town District	Changbo Country District
Geography factor dangerousness weights	0.2338	0.2727	0.1429	0.3506
Comprehensive dangerousness weights	0.1328	0.1550	0.0081	0.1992
Normalization weights	0.9305	1.0910	0.5714	1.4026

. Therefore, the final comprehensive dangerousness weights are also presented in Table 1.

The weights in the proposed Model 7 are presented in

Table 2. The weights in the proposed model.

Distance:
${\alpha }_1$	${\alpha }_2$
0.41	0.59
Urbanization Level:
${\beta }_1$	${\beta }_2$	${\beta }_3$
0.38	0.46	0.16
Hospital:
${\rho }_1$	${\rho }_2$	${\rho }_3$	${\rho }_4$	${\rho }_5$	${\rho }_6$	${\rho }_7$
0.13	0.09	0.21	0.17	0.18	0.15	0.07
School:
${\varrho }_1$	${\varrho }_2$	${\varrho }_3$	${\varrho }_4$	${\varrho }_5$	${\varrho }_6$	${\varrho }_7$	${\varrho }_8$	${\varrho }_9$	${\varrho }_{10}$
0.08	0.05	0.1	0.12	0.15	0.11	0.1	0.08	0.15	0.04

, which can be gained through a similar way as above.

6.2. Case Solution

The developed algorithm was run using software MATLAB 7.0 on an Inter Core 2, 2.00-GHz clock pulse with 8192 MB memory. The algorithmic parameters for the case problem were set as follows: swarm_size $S=30$, iteration_max $T=150$, inertia weight_max $w^{max}=0.9$, inertia weight_min $w^{min}=0.1$, personal best position acceleration constant $c_p=0.5$, global best position acceleration constant $c_g=0.5$, local best position acceleration constant $c_l=0.2$ and near neighbor best acceleration constant $c_n=0.2$.

After 150 iterations of MOAGLNPSO, the algorithm termination was achieved within 6 min on an average of 10 runs. Thus, the time is acceptable. The optimal solutions (i.e., the optimal gathered village locations) are shown in

Table 3. The non-dominating gathered village location.

No.	1	2	3	4	5	6
E	99${}^{\circ }$11${}^{\prime }$17.59${}^{″}$	99${}^{\circ }$10${}^{\prime }$57.15${}^{″}$	99${}^{\circ }$10${}^{\prime }$54.57${}^{″}$	99${}^{\circ }$19${}^{\prime }$50.82${}^{″}$	99${}^{\circ }$11${}^{\prime }$10.13${}^{″}$	99${}^{\circ }$25${}^{\prime }$43.57${}^{″}$
N	30${}^{\circ }$06${}^{\prime }$9.45${}^{″}$	30${}^{\circ }$06${}^{\prime }$7.01${}^{″}$	30${}^{\circ }$05${}^{\prime }$59.96${}^{″}$	30${}^{\circ }$14${}^{\prime }$11.36${}^{″}$	30${}^{\circ }$05${}^{\prime }$14.79${}^{″}$	30${}^{\circ }$19${}^{\prime }$24.16${}^{″}$
No.	7	8	9	10	11	12
E	99${}^{\circ }$10${}^{\prime }$33.40${}^{″}$	99${}^{\circ }$13${}^{\prime }$35.10${}^{″}$	99${}^{\circ }$10${}^{\prime }$34.40${}^{″}$	99${}^{\circ }$11${}^{\prime }$40.18${}^{″}$	99${}^{\circ }$10${}^{\prime }$48.31${}^{″}$	99${}^{\circ }$12${}^{\prime }$13.86${}^{″}$
N	30${}^{\circ }$06${}^{\prime }$19.17${}^{″}$	30${}^{\circ }$04${}^{\prime }$20.63${}^{″}$	30${}^{\circ }$06${}^{\prime }$20.24${}^{″}$	30${}^{\circ }$06${}^{\prime }$9.18${}^{″}$	30${}^{\circ }$05${}^{\prime }$57.15${}^{″}$	30${}^{\circ }$07${}^{\prime }$20.28${}^{″}$
No.	13	14	15	16	17	18
E	99${}^{\circ }$06${}^{\prime }$54.44${}^{″}$	99${}^{\circ }$10${}^{\prime }$39.39${}^{″}$	99${}^{\circ }$20${}^{\prime }$52.12${}^{″}$	99${}^{\circ }$24${}^{\prime }$4.93${}^{″}$	99${}^{\circ }$11${}^{\prime }$32.53${}^{″}$	99${}^{\circ }$14${}^{\prime }$11.82${}^{″}$
N	30${}^{\circ }$01${}^{\prime }$5.31${}^{″}$	30${}^{\circ }$02${}^{\prime }$59.74${}^{″}$	30${}^{\circ }$12${}^{\prime }$43.65${}^{″}$	30${}^{\circ }$17${}^{\prime }$38.32${}^{″}$	30${}^{\circ }$06${}^{\prime }$13.32${}^{″}$	30${}^{\circ }$08${}^{\prime }$4.38${}^{″}$
No.	19	20	21	22	23	24
E	99${}^{\circ }$28${}^{\prime }$6.17${}^{″}$	99${}^{\circ }$11${}^{\prime }$35.93${}^{″}$	99${}^{\circ }$21${}^{\prime }$5.56${}^{″}$	99${}^{\circ }$23${}^{\prime }$59.52${}^{″}$	99${}^{\circ }$16${}^{\prime }$6.03${}^{″}$	99${}^{\circ }$10${}^{\prime }$45.64${}^{″}$
N	30${}^{\circ }$18${}^{\prime }$43.55${}^{″}$	30${}^{\circ }$06${}^{\prime }$10.72${}^{″}$	30${}^{\circ }$13${}^{\prime }$38.86${}^{″}$	30${}^{\circ }$15${}^{\prime }$3.48${}^{″}$	30${}^{\circ }$12${}^{\prime }$17.22${}^{″}$	30${}^{\circ }$05${}^{\prime }$57.08${}^{″}$
No.	25	26	27
E	99${}^{\circ }$11${}^{\prime }$59.57${}^{″}$	99${}^{\circ }$17${}^{\prime }$59.70${}^{″}$	99${}^{\circ }$17${}^{\prime }$59.70${}^{″}$
N	30${}^{\circ }$06${}^{\prime }$47.81${}^{″}$	30${}^{\circ }$15${}^{\prime }$43.98${}^{″}$	30${}^{\circ }$16${}^{\prime }$20.28${}^{″}$

including all the non-dominating solutions. Table 3 shows an optimal gathered village location set with 27 solutions expressing in east-north longitude. The local governments as the decision makers are able to choose their preferred plan from the set. If they prefer to pursue the shorter distance $F_1:D$ as more important, they would choose the minimum distance plan, and vice versa. Although, there are fuzzy numbers in Model 7, they are easily transformed into equivalent crisp forms by many fuzzy theories. This will not influence the decision result.

Since the proposed model has three conflicting objectives, judging its non-dominating solution is more complex compared with the two-dimensional situation. In order to simplify the solution and describe the Pareto optimal front more intuitively, this study is given a judgment rule as follows.

Let the two minimum objective functions in the each iteration be the $x,y$ axis. Meanwhile, let the maximum objective function be the z axis. Therefore, a 3D Pareto optimal solution can be presented as a cuboid in the three-dimensional stereogram, which is shown in Figure 6A. If we transform it into the plane projection, the tendency of the optimal solution can be presented as in Figure 6B.

Through Figure 6B, it can be concluded that if one solution of a 3D Pareto optimal solution is close to $f\left(x\right)=\frac{1}{x}(x>0)$, it is more optimal. With the characteristic of the PSO algorithm, if the solution amount increases, the obtained 3D Pareto optimal solution is a better choice.

Based on the above judgment rule, the below Figure 7 describes the iterative progress of the Pareto optimal solutions.

Above all, the local governments in NSTR can obtain the optimal gathered village locations. Moreover, in order to realize the different aims of urbanization, the decision maker can choose the corresponding different schemes according to the development situations and change with the non-dominating solutions.

6.3. Analytic Results of the Proposed Approach

(1) Fast-growing economy:

The gathered villages are in close proximity to the 318 National Highway. Therefore, the convenient transportation will raise the productivity of Chinese medicinal herbs, such as Cordyceps sinensis, Notopterygium root, Rheum officinale, Chinese rhubarb, etc. Furthermore, animal husbandry can be promoted, such as Tibetan pork, yak meat, and so on. Moreover, it is convenient for local villages to make a pilgrimage to develop a sustainable culture.

(2) Social services improvement:

Figure 8 shows that the optimal gathered villages are relatively near the central village, so the level of social services will be improved. The optimal villages can construct more efficient and sustainable infrastructures, including high-tech hospitals, high-level schools, and so on. Meanwhile, the local governments play an important role in urban development, so they will have dedicated efforts to formulating and implementing policies to promote local/neighborhood urbanization.

(3) Eco-environmental protection:

Sustainable use of natural resources plays a significant role in the process of local sustainable urbanization. The gathered villages are farther away from restricted development than the weak villages, so the ecological environment can be well protected. Moreover, the local governments can conduct sustainable planning to repair the deteriorated environment, and thus, the villagers can improve the living conditions and feel a greater level of happiness.

(4) Disaster risk mitigation:

The picture demonstrates that keeping the optimal villages away from the climate-induced high-risk geologic disaster areas is very important. Thus, it can avoid the injury and distribution of livelihood. Furthermore, the geologic hazards, such as earthquake, debris flow, mountain collapse and so on, will lead to the complete destruction of work places. The low-risk areas can protect people against joblessness and minimize the loss of industries and agriculture.

6.4. Comparative Analysis and Discussion

(1) Worthiness of modeling and solutions:

The proposed model embodies the comprehensive ESERR goal by integrating the objectives and constraints. Meanwhile, the MODM method is introduced to determine the Pareto optimal solution and provides more effective and non-dominated alternatives for the decision maker. Compared with the traditional weight-sum method of the multi-objective, the solutions in this paper have more values and reflect the users’ preference requirements. Therefore, it is more worthwhile.

The fuzzy random programming approach explicitly considers the entire range of bi-uncertain scenarios; thus, it conforms more to reality. Although it increases the complexity of modeling, the model is well brought into life. Therefore, the extra effort on modeling and solving fuzzy random programming is worthwhile.

(2) Efficiency of the algorithm:

For multiple objective optimization, the definition of quality is substantially more complex than for single objective optimization. For further expression for the efficiency of the non-dominating solutions, the solution amounts and four performance metrics are studied.

Table 4. Evaluation of the non-dominating solution.

Iteration	The Solution Amount	The Average Distance	The Distribution	The Extent	The Set Convergence
1	5	3.0988	0.6000	27.6237	0.6000
5	11	0.7681	0.6727	28.2246	0.9091
10	13	2.0580	0.6154	31.7128	0.9231
11	15	1.5690	0.6095	31.7128	0.8667
15	15	1.5694	0.6286	31.7128	0.9333
20	16	3.3987	0.8750	67.0213	0.9375
30	17	2.9558	0.8824	67.0213	0.9412
75	22	1.4833	0.9091	76.1165	0.9545
87	26	0.4224	0.9231	76.0754	0.9615
88	26	0.4225	0.9231	76.0754	0.9615
96	26	0.4225	0.9231	76.0754	0.9615
102	26	0.4224	0.9231	76.0754	0.9615
146	27	0.2035	0.9259	76.0754	0.9630
150	27	0.2035	0.9259	76.0754	1

shows the metrics of non-dominating solutions proposed above after 10 runs.

This paper compares basic multi-objective PSO and the developed algorithm. In the developed algorithm, the particle-represented solutions tie the particles of PSO with the problem’s solution. The hybrid particle-updating mechanism (i.e., more particle-updating orientations viz. $lbest$ and $nbest$ and inertia weight updating) successfully enhances the searching capability. As shown in

Table 5. Performance comparison of MOAGLNPSO and basic MOPSO.

Algorithm Type	Iteration	The Solution Amount	The Average Distance	The Distribution	The Extent	The Set Convergence
MOAGLNPSO	150	27	0.2035	0.9529	76.0754	1
MOPSO	150	24	0.2108	0.7246	65.6142	1

, the developed algorithm is a useful tool for problem solution by comparing with the basic one.

7. Conclusions

This paper studies VLP for PM of local/neighborhood urbanization in China for sustainable development. Under the comprehensive ESERR goal orientation, an integrated MODM approach is proposed to obtain the optimal gathered village locations. The fuzzy random multiple objective programming model is established primarily. Further, the transformation process mentioned in this study is first to obtain the equivalent fuzzy bi-level programming model. Second, a fuzzy EVM based on the fuzzy measure $Me$ is utilized suitably for the realistic uncertain decision making process. Afterwards, MOAGLNPSO is developed to solve the problem. Meanwhile, a Pareto optimal solution judgment criterion is proposed for the convenience of discussion. Finally, a case study is presented as an illustrative example of this problem. The results validate the worthiness of modeling and solutions and test the efficiency of the algorithm and parameters.

The contributions of this paper to the literature are: (1) This study discusses the emerging new challenge in China’s urbanization progress. In order to promote sustainable urbanization, the optimal gathered village location programming provides a more reasonable and practical description of the problem. (2) Although there are many works about climate-induced geological disaster, few papers consider it in the VLP model. Thus, this paper increases the awareness of the problem. (3) This paper uses fuzzy random house rebuilding price to describe the bi-uncertain situation. To the best of our knowledge, it has never been done before. (4) MOAGLNPSO with a suitable judgment criterion is developed as one of the useful tools to solve such a problem.

This work is original. However, there are still three areas suggested for future research. Firstly, a more detailed description of the objective functions needs to be investigated and developed further. In particular, considering the key study, a wider scope is necessary, such as: (1) how to apply the proposed approach to the Tibet region even in the Muslim-dominated Xinjiang region; (2) how to explore the role of the hukou system in China’s urbanization; (3) how to describe the influences of industrialization and environmental pollution in sustainable urbanization. Due to the diversity of different religious cultures, the research methods proposed in this paper can provide the theoretical basis for religious groups. Beside, this paper has been researching nearby migration, which does not involve inter-provincial migration, so the hukou system has little influence on local/neighborhood urbanization. For future research, it is necessary to consider the hukou system, which could slow down the migration in the cities and regulate urbanization. Moreover, to obtain better and more effective solutions with less memory and computing time, alternative approaches and algorithms (e.g., other exact approaches, (meta-)heuristics, evolutionary algorithms, etc.) could be explored. Finally, the follow-up task also should be considered, such as the reasonable infrastructure layouts and so on. All these areas are very important and equally worthy of concern.