A Stone Resource Assignment Model under the Fuzzy Environment

Abstract

This paper proposes a bilevel multiobjective optimization model with fuzzy coefficients to tackle a stone resource assignment problem with the aim of decreasing dust and waste water emissions. On the upper level, the local government wants to assign a reasonable exploitation amount to each stone plant so as to minimize total emissions and maximize employment and economic profit. On the lower level, stone plants must reasonably assign stone resources to produce different stone products under the exploitation constraint. To deal with inherent uncertainties, the object functions and constraints are defuzzified using a possibility measure. A fuzzy simulation-based improved simulated annealing algorithm (FS-ISA) is designed to search for the Pareto optimal solutions. Finally, a case study is presented to demonstrate the practicality and efficiency of the model. Results and a comparison analysis are presented to highlight the performance of the optimization method, which proves to be very efficient compared with other algorithms.

1. Introduction

The dust and the waste water from the stone industry can cause serious damage to the regional ecological environment. The overexploitation and the stone processing have resulted in the vegetation decrement, and the pollution of air and water in those areas with rich stone resources. The annual amount of waste generated include 700,000 tons of slurry waste as well as 1 million tons of solid waste. The consequent dumping of this waste in open areas has created several environmental problems and has negatively impacted agriculture, local inhabitants, and groundwater [1]. Therefore, it is urgent to normalize the quarrying and processing of the stone resource. Some technologies are introduced to save energy and reduce the emission in the stone industry by many scholars [2, 3]. Some other scholars [4–6] considered the use of the marble powder to reduce the waste, but few literatures discussed the quantitative relationship between the emission and the exploiting and processing amount. In fact, a reasonable assignment of stone resources could significantly reduce the emissions. This paper considers the government as the upper level and the stone plant as the lower level to develop a bi-level model. For the stone industry, the objectives of the government authority are to minimize environmental pollution and maximize social employment and economic revenue. This can be achieved by optimizing the amount of stone extracted and exploited between the participating plants, which are assumed to cooperate and act as a lower-level decision maker. Noting that industrial symbiosis implicitly requires the cooperative behavior of the participants [7, 8], the government can influence the plants by imposing disincentives by assigning different amounts to stone plants according to their production scale and clean technology level. The plants operate independently of each other. Each plant has its own goals, which are to maximize the profit from the sale of nano calcium carbonate, marble products, granite slabs, and man-made slabs and to minimize the emissions of stone dust and waste water.

To develop the bi-level optimization model for assigning the stone resources, some emission coefficients have to be effectively estimated. It is usually difficult to collect the exact data of emissions of stone dust and waste water when exploiting the stone mine and processing stone products. The fuzzy number is an efficient tool to describe the variables without crisp information. The membership function of fuzzy sets can be used to describe the possibility that emission coefficients take the value according to the experience of those people in the stone industry. Actually, there has been some studies describing the uncertainty by fuzzy sets. For example, Petrovic et al. [9] used fuzzy sets to describe the customer demand, supply deliveries along the SC and the external or market supply, and develop a supply chain model with fuzzy coefficients. Lee and Yao [10] fuzzify the demand quantity and the production quantity per day to solve the economic production quantity. These studies inspire us to use the fuzzy sets to interpret the vague and imprecise about the emissions of stone dust and waste water. For the fuzzy bi-level optimization problem, a satisfactory (near-optimal or “satisficing”) solution can be reached by providing tolerances in the objective functions and constraints and by defining corresponding degrees of satisfaction through membership functions to indicate the preference of the decision makers which is typical of decision making in a fuzzy environment [11]. The followers then communicate their results to the leader, who modifies his goals and control variables if the original tolerances are not met. The process continues iteratively until a solution which satisfies the goals of both leader and follower is reached.

A bi-level multiobjective model with fuzzy coefficients is always an NP hard problem, and it is especially difficult for nonlinear bi-level programming under a fuzzy environment to find a numerical solution. Some existing methods mainly focus on metaheuristics which include the genetic algorithm [12], the simulated annealing [13], and the hybrid tabu-ascent algorithm [14]. However, as these need to be designed for single-objective problems with crisp coefficients, it is difficult to find a usual or normal pattern for a bi-level model with fuzzy coefficients. This paper proposes an improved simulated annealing based on a fuzzy simulation to search for a Pareto optimal solution after a possibilistic check. The following sections of this paper are organized as follows. In Section 2, the reason a bi-level multi-objective model is used to optimize the stone industry is explained. The process of data fuzzification is introduced in detail. A possibilistic bi-level multi-objective programming model is developed. In Section 3, a fuzzy simulation-based improved simulated algorithm is proposed to solve the bi-level multi-objective programming model with fuzzy coefficients. In Section 4, a practical case is presented to show the significance of the proposed models and algorithms. Finally, conclusions are given in Section 5.

2. Mathematical Modelling

In order to develop the mathematical model, some basic background and descriptions are introduced.

2.1. Key Problems Description

For the stone industry, local government and stone plants play important roles to perform the responsibilities, respectively. Government has the authority to decide the amount that should be exploited and then needs to make a sustainable plan to avoid overexploitation and pollution. On the other hand, stone plants need to make the production plan according to the stone quota that government gives. As shown in Figure 1, the local government has environmental protection and maximum employment as its most important goals and then considers the revenue. On the lower level, the stone-material plants usually consider economic profit as their first goal. Due to a limitation on the amount that can be extracted and the environmental protection requirements, they also have to consider minimizing emissions. At the same time under a government policy and according to capacity, they also need to think about employment. In addition, it is necessary to increase the investment to improve emission reduction capacity to satisfy sustainable development requirements. Considering the above, the problem should be regarded as a bi-level optimization model in which the government authority is the upper-level decision maker and the stone plants are the lower-level decision makers. It is assumed that there is a perfect exchange of information between all the participants such that the objectives and constraints are known.

Figure 1 

Stone industry decision making hierarchy.

As shown in Figure 2, the granite is first exploited from the stone mine and cut into the primary products such as granite slabs and man-made slabs; then, these are processed into the floor or other products. Scrap materials are usually processed into fine powder calcium carbonate and nano calcium carbonate to meet market demand. During the complete process, a great deal of stone dust and waste water are produced. Since it is technically difficult to collect the exact data of emissions, we usually make a rough estimation by the difference of weight before and after exploiting and processing and then look for the possibility for every weight by the professional advices in the stone industry. Therefore, the fuzzy number is an efficient tool to describe this situation by its membership function. Actually, the fuzzy environment has been successfully studied and applied in many areas, such as flow shop scheduling problem [17], supply chain problem [18], and assignment problem [19]. These studies show the necessity of considering fuzzy environment in practical problems. It is also the motivation for considering fuzzy environment in the stone resources assignment problems.

Figure 2 

Process flow diagram for granite processing operations.

2.2. Assumptions and Notations

Before the data fuzzification and developing the optimization model, some assumptions should be introduced.(1)Emission of stone dust and waste water is proportional to the amount of stone processed into products.(2)Employment level is also proportional to the amount of stone processed into products.(3)The constant cost of product 𝑗 only exists when the stone-material plant produces product 𝑗.(4)Since the government endows different subsidies to plants, it is assumed that each plant has its own tax rate 𝑆_𝑖 and the tax is proportional to the turnover of all stone products.

The notations are used to describe the subsidy model in the investigation are referred to in the Abbreviations Section.

2.3. Data Fuzzification Based on Crossover Validation Test

Often there is little historical data to describe emission reduction due to the raw development of the last decade. For example, some research considers the transport cost as uncertain coefficients because of the changing weather and the unpredictable road condition [20, 21]. In this paper, the emission coefficients cannot be estimated using the statistical methods and then have to fuzzificated according to those insufficient data.

The essence of fuzzification is to find an approximate membership function to describe the fuzzy number [22]. Many scholars have described some methods of determining the membership functions that are essentially based on direct methods of inquiry made on human beings and corrected using indirect methods [23, 24]. Some other scholars propose the automatic methods to determine the membership functions when no expert is available or in the case when there are so many data [25]. In the present paper, we will propose the fuzzification methods by combining the 5-parameter membership function and crossover validation test. Taking the stone dust emission coefficients 𝐸𝑑 as an example, the process for fuzzificating can be summarized as follows.

Step 1. Split the data set 𝑆 of stone dust emission coefficients 𝐸𝑑 into a training set 𝑆_tr and a validation set 𝑆_𝑣.

Step 2. Find the smallest, middle, and largest data in 𝑆; denote them 𝐸𝑑_𝑠, 𝐸𝑑_𝑚, and 𝐸𝑑_𝑙, respectively.

Step 3. Compute the left and right slops for the data in 𝑆_tr by the following equations; respectively, 𝐸𝑑_𝛼(𝑥)=𝐸𝑑_𝑠+𝐸𝑑_𝑚−2𝑥𝐸𝑑_𝑠+𝐸𝑑_𝑚,if𝐸𝑑_𝑠≤𝑥≤𝐸𝑑_𝑠+𝐸𝑑_𝑚2,2𝑥−𝐸𝑑_𝑠+𝐸𝑑_𝑚2𝐸𝑑_𝑚,if𝐸𝑑_𝑠+𝐸𝑑_𝑚2≤𝑥≤𝐸𝑑_𝑚,𝐸𝑑_𝛽(𝑥)=𝐸𝑑_𝑙+𝐸𝑑_𝑚−2𝑥𝐸𝑑_𝑙+𝐸𝑑_𝑚,if𝐸𝑑_𝑠≤𝑥≤𝐸𝑑_𝑙+𝐸𝑑_𝑚2,2𝑥−𝐸𝑑_𝑙+𝐸𝑑_𝑚2𝐸𝑑_𝑚,if𝐸𝑑_𝑙+𝐸𝑑_𝑚2≤𝑥≤𝐸𝑑_𝑚. (2.1)(2.2) Then we get the set of left slops 𝐸𝑑^𝐿={𝛼∣𝐸𝑑_𝛼(𝑥),𝑥∈𝑆_tr,𝑥≤𝐸𝑑_𝑚} and the set of left slops 𝐸𝑑^𝑅={𝛽∣𝐸𝑑_𝛽(𝑥),𝑥∈𝑆_tr,𝑥≥𝐸𝑑_𝑚}.

Step 4. Define the membership function as follows: 𝜇_𝐸𝑑(𝑥)=0,if𝑥<𝐸𝑑_𝑠,2^𝛼−1𝑥−𝐸𝑑_𝑠𝐸𝑑_𝑚−𝐸𝑑_𝑠^𝛼,if𝐸𝑑_𝑠≤𝑥≤𝐸𝑑_𝑠+𝐸𝑑_𝑚2,1−2^𝛼−1𝐸𝑑_𝑚−𝑥𝐸𝑑_𝑚−𝐸𝑑_𝑠^𝛼,if𝐸𝑑_𝑠+𝐸𝑑_𝑚2<𝑥<𝐸𝑑_𝑚,1−2^𝛽−1𝑥−𝐸𝑑_𝑚𝐸𝑑_𝑙−𝐸𝑑_𝑚^𝛽,if𝐸𝑑_𝑚≤𝑥<𝐸𝑑_𝑙+𝐸𝑑_𝑚2,2^𝛽−1𝐸𝑑_𝑙−𝑥𝐸𝑑_𝑙−𝐸𝑑_𝑚^𝛽,if𝐸𝑑_𝑙+𝐸𝑑_𝑚2≤𝑥≤𝐸𝑑_𝑙,0,if𝑥>𝐸𝑑_𝑙, (2.3) where 𝛼∈𝐸𝑑^𝐿 and 𝛽∈𝐸𝑑^𝑅.

Step 5. Take all the data in 𝑆_𝑣 in the above equation and compute the membership 𝜇_𝐸𝑑(𝑥;𝛼,𝛽), where 𝑥∈𝑆_𝑣, 𝛼∈𝐸𝑑^𝐿 and 𝛽∈𝐸𝑑^𝑅.

Step 6. Carry out the crossover validation test proposed by Kohavi [26]. Compute the memberships of 𝑥_𝑖∈𝑆_𝑣 for any combination (𝛼,𝛽)∈(𝐸𝑑^𝐿,𝐸𝑑^𝑅). Then compute the percentage of correct results by the following equation: PCC_𝑣1=100×𝑁𝑉_{𝑥𝑖,𝜇𝐸𝑑𝑥𝑖;𝛼,𝛽∈𝑆𝑣}𝛿𝑥_𝑖,𝜇_𝐸𝑑𝑥_𝑖;𝛼,𝛽, (2.4) where PCC_𝑣 denotes the percentage of correct results over the validation set 𝑆_𝑣, 𝑁𝑉 is the number of data points in validation set 𝑆_𝑣, and 𝛿(𝜇_𝐸𝑑(𝑥_𝑖;𝛼₁,𝛽₁),𝜇_𝐸𝑑(𝑥_𝑖;𝛼₂,𝛽₂))=1 if 𝜇_𝐸𝑑(𝑥_𝑖;𝛼₁,𝛽₁)=𝜇_𝐸𝑑(𝑥_𝑖;𝛼₂,𝛽₂), while 𝛿(𝜇_𝐸𝑑(𝑥_𝑖;𝛼₁,𝛽₁),𝜇_𝐸𝑑(𝑥_𝑖;𝛼₂,𝛽₂))=0 if 𝜇_𝐸𝑑(𝑥_𝑖;𝛼₁,𝛽₁)≠𝜇_𝐸𝑑(𝑥_𝑖;𝛼₂,𝛽₂).

Step 7. Find the combination (𝛼,𝛽) by which the largest percentage of correct results can be obtained when carrying out the crossover validation test with each other. Then we get the membership function.

2.4. Model Formulation

The bi-level multiobjective optimization model under a fuzzy environment for assigning stone resources can be mathematically formulated as follows.

2.4.1. Government Model (Upper Level)

As the upper level, the government has the obligation to protect the local environment, solve employment issues, and promote economic revenue. Generally, the following goals are usually considered by the government.

To achieve minimum emissions, including the stone dust (∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖) when all plants exploit the stone mine, the stone dust (∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗) when plants produce stone products, and total waste water (∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗) when all the plants exploit the stone mine and produce those stone products is the first objective. Since 𝐸𝑑_𝑖, 𝑒𝑑_𝑖𝑗, and 𝑒𝑤_𝑖𝑗 are all fuzzy numbers which are obtained by fuzzification due to insufficient historical data, it is usually difficult to derive precise minimum emissions, and decision makers only require a minimum objective (𝐹₁) under some possibilistic level (𝛿^𝑈₁) [27]. Hence the following possibilistic objective function and constraint are derived: min𝐹₁ (2.5) subject to (s.t.) Pos^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐹₁≥𝛿^𝑈₁, (2.6) where Posis the possibility measure proposed by Dubois and Prade [28] and 𝛿^𝑈 is the possibilistic level representing the possibility that decision makers achieve the minimum objective. All fuzzy arithmetic in (2.6) and the following equations come from the operation proposed by Kaufmann and Gupta [29].

To achieve maximum employment 𝐹₂ which consisted of constant workers (𝑃_𝑖) and variable workers (𝑝_𝑖𝑗𝑋_𝑖𝑗), the following objective function is obtained: max𝐹₂=^𝑚_𝑖=1^𝑛_𝑗=1𝑝_𝑖𝑗𝑋_𝑖𝑗+𝑃_𝑖. (2.7)

To achieve the maximum economic output which can be obtained by multiplying unit amount (𝑐_𝑗), conversion rate (𝜃_𝑖𝑗), and amount of stone (𝑋_𝑖𝑗), the following objective function is obtained: max𝐹₃=^𝑚_𝑖=1𝑆_𝑖^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗. (2.8)

Generally, some mandatory conditions must be satisfied when the government makes a decision. These are listed as follows.

The total exploration quantity (∑^𝑚_𝑖=1𝑌_𝑖) cannot exceed the upper limitation (𝑅^𝑈) of the total stone resources in the region: ^𝑚_𝑖=1𝑌_𝑖≤𝑅^𝑈. (2.9)

The stone dust from exploiting (∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖) and producing (∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗) and the waste water (∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗) should be less than the predetermined levels (𝐸𝐷^𝑈 and 𝐸𝑊^𝑈) in order to guarantee air and water quality. Two constraints are derived under the possibilistic levels (𝛿^𝑈₂ and 𝛿^𝑈₃): Pos^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝐷^𝑈≥𝛿^𝑈₂,Pos^𝑚𝑛_𝑖=1𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝑊^𝑈≥𝛿^𝑈₃. (2.10)(2.11)

The output of some products (∑^𝑚_𝑖=1𝜃_𝑖𝑗𝑋_𝑖𝑗) should meet the market demand (𝐷^𝐿_𝑗). For example, the nano calcium carbonate is very popular in many areas, so the stone plants should provide enough output to meet the demand: ^𝑚_𝑖=1𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝐷^𝐿_𝑗∀𝑗. (2.12)

2.4.2. Plant Model (Lower Level)

On the lower level, the stone plants usually pursue maximum profit and then try to reduce the emissions. Thus, the following two objectives are introduced.

Each plant wishes to achieve maximum profit which consisted of total sales (∑^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗) minus the production cost (𝑓(𝑋_𝑖𝑗)) and the inventory cost (ℎ_𝑖(𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗)); then the following objective function is determined: max𝐻¹_𝑖=^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗−^𝑛_𝑗=1𝑓𝑋_𝑖𝑗−ℎ_𝑖𝑌_𝑖−^𝑛_𝑗=1𝑋_𝑖𝑗, (2.13) where 𝑓(𝑋_𝑖𝑗) is the production-cost function as follows [27]: 𝑓𝑋_𝑖𝑗=𝑡_𝑖𝑗𝑋_𝑖𝑗+𝐶_𝑖𝑗,if𝑋_𝑖𝑗>0,0,if𝑋_𝑖𝑗=0. (2.14) Every plant also wishes to achieve minimum emissions. However, since the emissions 𝑒𝑑_𝑖𝑗 and 𝑒𝑤_𝑖𝑗 are fuzzy numbers, it is usually difficult to determine the precise minimum emissions, and decision makers only require a minimum objective (𝐻²_𝑖) under some possibilistic level (𝜎^𝐿_𝑖). Hence, the possibilistic constraint is as follows: min𝐻²_𝑖 (2.15) subject to Pos^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐻²_𝑖≥𝜎^𝐿_𝑖, (2.16) where 𝜎^𝐿_𝑖 is the possibilistic level under which decision makers require the minimum objective.

Since production in all the plants is influenced by government policy and market demand, there are some conditions that need to be satisfied.

The amount used for production (∑^𝑛_𝑗=1𝑋_𝑖𝑗) should not exceed the total limitation (𝑌_𝑖): ^𝑛_𝑗=1𝑋_𝑖𝑗≤𝑌_𝑖. (2.17)

The inventory amount (𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗) should not exceed the maximum limitation (𝐼𝑉^𝑈_𝑖): 𝑌_𝑖−^𝑛_𝑗=1𝑋_𝑖𝑗≤𝐼𝑉^𝑈_𝑖. (2.18)

The production cost which consisted of two parts including product cost (∑^𝑛_𝑗=1𝑓(𝑋_𝑖𝑗)) and total inventory cost (∑^𝑛_𝑗=1ℎ_𝑖(𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗)) should not exceed the predetermined level (𝑃𝐶^𝑈_𝑖): ^𝑛_𝑗=1𝑓𝑋_𝑖𝑗+ℎ_𝑖𝑌_𝑖−^𝑛_𝑗=1𝑋_𝑖𝑗≤𝑃𝐶^𝑈_𝑖. (2.19)

Some products (𝜃_𝑖𝑗𝑋_𝑖𝑗) should not be less than the lowest production level (𝑃^𝐿_𝑖𝑗) in plant 𝑖: 𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝑃^𝐿_𝑖𝑗. (2.20)

2.4.3. Bilevel Model

In such a complicated system, both the leader and the followers should simultaneously consider the objectives and constraints and then make the decision. Therefore, from (2.5)∼(2.20), the complete bi-level multiobjective optimization model under a fuzzy environment is as follows: min𝐹₁max𝐹₂=∑^𝑚∑_𝑖=1^𝑛_𝑗=1𝑝_𝑖𝑗𝑋_𝑖𝑗+𝑃_𝑖max𝐹₃=∑^𝑚_𝑖=1𝑆∑_𝑖^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗s.t.Pos∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐹₁≥𝛿^𝑈₁Pos∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝐷^𝑈≥𝛿^𝑈₂Pos∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝑊^𝑈≥𝛿^𝑈_3𝑚𝑖∑₌₁𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝐷^𝐿_𝑗∀𝑗max𝐻¹_𝑖=∑^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗−∑^𝑛_𝑗=1𝑓𝑋_𝑖𝑗−ℎ_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗min𝐻²_𝑖s.t.Pos∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐻²_𝑖≥𝜎^𝐿_𝑖𝑛𝑗∑₌₁𝑋_𝑖𝑗≤𝑌_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗≤𝐼𝑉^𝑈_𝑖𝑛𝑗∑₌₁𝑓𝑋_𝑖𝑗+ℎ_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗≤𝑃𝐶^𝑈_𝑖𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝑃^𝐿_𝑖𝑗. (2.21)

3. Solution Approach

Generally, bi-level programming is an NP-hard problem, and it is difficult to determine an optimal solution [30–32]. In the proposed model, decision makers on the upper and lower levels have to face more than two conflicting objectives and then make a decision under a fuzzy environment. This significantly increases the difficulty of finding an optimal strategy for both the upper and lower levels. Therefore, the fuzzy simulation-based improved simulated annealing (FS-ISA) is designed to solve the bi-level optimization model with fuzzy coefficients.

3.1. Fuzzy Simulation for Possibilistic Constraints

Fuzzy simulation is usually proposed to approximate the possibility measure according to the membership function of a fuzzy number [27]. Taking the constraint (2.5) as an example, we will introduce the key principle of the fuzzy simulation and find the minimum 𝐹₁ such that the constraint holds.

Let 𝑌^∗_𝑖 and 𝑋^∗_𝑖𝑗 be predetermined feasible solutions for 𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛, which will be regarded as input variables. Firstly, set 𝐹₁=𝑀, where 𝑀 is a sufficiently large number. Secondly, randomly generate 𝜆_𝑖, 𝜂_𝑖𝑗, and 𝜅_𝑖𝑗 from the 𝛿^𝑈₁-level set of the fuzzy numbers 𝐸𝑑_𝑖, 𝑒𝑑_𝑖𝑗, and 𝑒𝑤_𝑖𝑗, respectively. Thirdly, compute the value 𝑓=∑^𝑚_𝑖=1𝜆_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1(𝜂_𝑖𝑗+𝜅_𝑖𝑗)𝑋_𝑖𝑗. If 𝐹₁>𝑓, replace it with 𝑓. Finally, repeat this process for 𝑁 times. The value 𝐹₁ is regarded as the estimation. Then the simulation process can be summarized in Procedure 1.

Input: Decision variables 𝑌𝑖 and 𝑋𝑖𝑗 Output: The minimum 𝐹1 Step 1. Set 𝐹1=𝑀, where 𝑀 is sufficiently large number; Step 2. Randomly generate 𝜆𝑖, 𝜂𝑖𝑗 and 𝜅𝑖𝑗 from the 𝛿𝑈1-level set of the fuzzy numbers 𝐸𝑑𝑖, 𝑒𝑑𝑖𝑗 and 𝑒𝑤𝑖𝑗 respectively; Step 3. Compute 𝑓=∑𝑚𝑖=1𝜆𝑖𝑌𝑖+∑𝑚𝑖=1∑𝑛𝑗=1(𝜂𝑖𝑗+𝜅𝑖𝑗)𝑋𝑖𝑗 and replace 𝐹1 with 𝑓 provided that 𝐹1>𝑓; Step 4. Repeat the second and third steps 𝑁 times; Step 5. Return 𝐹1.

Algorithm 1 

Fuzzy simulation for possibilistic constraints.

Sometimes, we need to check whether a solution satisfies the possibilistic constraint. This means that we need to compute the possibility and compare it with the predetermined possibilistic level. Then another simulation is applied to check the constraint. Taking the constraint (2.10) as an example, we will introduce how to simulate the possibility 𝐿=Pos{∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝑊^𝑈}.

Let 𝑋^∗_𝑖𝑗 be predetermined solution for 𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛, which will be regarded as input variables. Give a lower estimation of the possibility 𝐿, denoted by 𝛿. Then we randomly generate 𝜅_𝑖𝑗 from the 𝛿-level set of the fuzzy numbers 𝑒𝑤_𝑖𝑗. If the 𝛿-level set is not easy for a computer to describe, we can give a larger region, for example, a hypercube containing the 𝛿-level set. Certainly, the smaller the region, the more effective the fuzzy simulation. Now we set 𝜇=max𝜇𝑒𝑤_𝑖𝑗,𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛. (3.1) If ∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝜅_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝑊^𝑈 and 𝐿<𝜇, then we set 𝐿=𝜇. Repeat this process 𝑁 times. The value 𝐿 is regarded as an estimation of the possibility. Then the process for constrain check can be summarized in Procedure 2.

Input: Decision variables 𝑋𝑖𝑗 Output: The possibility Step 1. Set 𝐿=𝛼 as a lower estimation Step 2. Randomly generate 𝜅𝑖𝑗 from the 𝛿-level set of the fuzzy numbers 𝑒𝑤𝑖𝑗 Step 3. Set 𝜇=𝑚𝑎𝑥{𝜇𝑒𝑤𝑖𝑗,𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛} Step 4. If ∑𝑚𝑖=1∑𝑛𝑗=1𝜅𝑖𝑗𝑋𝑖𝑗≤𝐸𝑊𝑈 and 𝐿<𝜇, set 𝐿=𝜇 Step 5. Repeat the second and third steps 𝑁 times Step 6. Return 𝐿.

Algorithm 2 

Possibilistic constraint check.

3.2. Fuzzy Simulation-Based Improved Simulated Annealing Algorithm

Simulated annealing algorithm (SA) is proposed for the problem of finding, numerically, a point of the global optimization of a function defined on a subset of a 𝑛-dimensional Euclidean space [33–35]. Many fruitful results are obtained in the past decades. Steel [36, 37] calls simulated annealing the most exciting algorithmic development of the decade. For the multiobjective optimization problems, some scholars have introduced many progressive simulated annealing algorithms to solve them. Especially, Suppapitnarm et al. [38] designed a simulated annealing algorithm along with archiving the Pareto optimal solutions coupled with return to base strategy (SMOSA) to explore the trade-off between multiple objectives in optimization problems. Suman and Kumar [39, 40] introduced four simulated annealing algorithms including SMOSA, UMOSA, PSA, and WMOSA to solve multiobjective optimization of constrained problems with varying degree of complexity and then proposed a new algorithm PDMOSA. Sanghamitra et al. [41] proposed a simulated annealing-based multiobjective optimization algorithm (AMOSA) that incorporates the concept of archive in order to provide a set of trade-off solutions for the problem under consideration.

In the following part, we will incorporate the fuzzy simulation into the SMOSA algorithm proposed by Suppapitnarm and Parks [16] and use the interactive method to search the Pareto optimal solution for the bi-level multiobjective optimization with fuzzy possibilistic constraints. Take the problem (A.10) as an example and denote 𝐗_𝑖=(𝑋_𝑖1,𝑋_𝑖2,…,𝑋_𝑖𝑛) and 𝐘=(𝑌₁,𝑌₂,…,𝑌_𝑚), the process of FS-ISA can be summarized in Procedure 3.

Input: The initial temperature 𝑡0 Output: Pareto-solution 𝑌∗𝑖 and 𝑋∗𝑖𝑗 for all 𝑖 and 𝑗 Step 1. Randomly generate a feasible solution 𝑌 according to the fuzzy simulation for possibilistic constraints and take it as the initial parameter for the lower level; Step 2. Solve all the multiobjective optimization problems on the lower level by SMOSA based on the fuzzy simulation and we obtain the Pareto optimal solution 𝑋𝑖 for all 𝑖. Put 𝐺=(𝑋𝑇1,𝑋𝑇2,…, 𝑋𝑇𝑚,𝑌𝑇) into a Pareto set of solutions and compute all objective values of the upper and lower levels; Step 3. Generate a new solution 𝐺1=(𝑋11𝑇,𝑋21𝑇,…,𝑋𝑚1𝑇,𝑌1𝑇) in the neighborhood of 𝐺 by the random perturbation; Step 4. Check the feasibility by fuzzy simulation according to all the constraints on both levels. If not, return to Step 3; Step 5. Compute the objective values on both level, respectively. Compare the generated solution with all solutions in the Pareto set and update the Pareto set if necessary; Step 6. Replace the current solution 𝐺 with the generated solution 𝐺1 if 𝐺1 is archived and go to Step 7; Step 7. Accept the generated solution 𝑌1 as the input solution for the lower level if it is not archived with the probability: 𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦(𝑝)=𝑚𝑖𝑛(1,𝑒𝑥𝑝{−Δ𝑠𝑖/𝑡𝑖}), where Δ𝑠𝑖=𝐹∗1(𝐺)−𝐹∗1(𝐺1)+𝐹2(𝐺1)−𝐹2(𝐺)+𝐹3(𝐺1)−𝐹3(𝐺). If the generated solution is accepted, take it into the lower level and solve them. Then we get a new solution 𝐺1∗=(𝑋11𝑇∗,𝑋21𝑇∗,…,𝑋𝑚1𝑇∗,𝑌1𝑇) and put it into the Pareto set. If not, go to Step 9; Step 8. Compare 𝐺1∗ and 𝐺1 according to the evaluation function based on the compromise approach proposed by Xu and Li [15]. If 𝐺1∗ is more optimal than 𝐺1, let 𝐺=𝐺1∗. If not 𝐺=𝐺1 Step 9. Periodically, restart with a randomly selected solution from the Pareto set. While periodically restarting with the archived solutions, Suppapitnarm et al. [16] have recommended biasing towards the extreme ends of the trade-off surface; Step 10. Periodically reduce the temperature by using a problem-dependent annealing schedule Step 11. Repeat steps 2–10, until a predefined number of iterations is carried out.

Algorithm 3 

FS-ISA algorithm for bi-level multi-objective programming.

Above all, the whole procedure of FS-ISA for bi-level multiobjective optimization problems with fuzzy coefficients is described in Figure 3.

Figure 3 

Flow chart for FS-ISA.

4. A Case Study

In the following, a practical example in China is introduced to demonstrate the complete modelling and algorithm process.

4.1. Data and Computation

Yingjing County is a famous county in China for its rich mineral products. The granite in this area has stable physical and chemical properties so that it can be processed into many useful stone products, mainly including granite slabs, man-made composite slabs, granite sands, and nano calcium carbonates (see Figure 4). Figure 5 shows the actual stone industry process from exploitation to production. Due to the vegetation deterioration, air and water pollution, and an aggravation of the ecological environment caused by the disordered exploitation and production manner, it is urgent for both the Yingjing government and the stone plants to optimize the assignment strategy.

Figure 4 

Products from the granite.

Figure 5 

Basic flow chart of stone industry.

Up to now, around 1 billion m³ of granite available is being exploited in Yingjing County according to the investigation. At present, only 7 stone plants have been built in this county, but the government plans to extend this to 10 stone plants in 2013, with all ten plants sharing the granite resource. From the historical data, stone dust and waste water emission coefficients are fuzzificated and crossover validation tested. The test demonstrates that when the membership function is triangular, the percentage of correct results is the largest 92.32%. Therefore, the emission coefficients are regarded as fuzzy numbers in Tables 1 and 3. According to the environmental sector in this county, stone dust emissions should not exceed 2500 tonnes and waste water emission should not exceed 2500 tonnes. Although it is difficult to satisfy the constrained index in a short time due to uncertainty, the possibility of holding the two constraints should not be less than 0.9 which indicates that the possibilistic levels 𝛿^𝑈₂ and 𝛿^𝑈₃ for the government should also be 0.9. For total emissions, the environmental sector requires the minimum objective to be under the possibilistic level 𝛿^𝑈₁=0.85. As the demand and the price of the four stone products sharply increase, the government requires that their output from all the plants should at least satisfy the basic market demand 𝐷^𝐿_𝑗(𝑗=1,…,4) as in Table 2. Each stone products’ unit price is in Table 2. For the 10 stone plants, the inventory and the production upper limitations for each plant are listed in Table 1. The possibilistic level 𝛿^𝐿_𝑖 that plant 𝑖 needs to obtain for minimum emissions is in Table 1. Since every plant has a different capacity for controlling emissions, the fixed and unit variable cost, emission coefficients, and constant costs are different as outlined in Table 3. The transformation rate 𝜃_𝑖𝑗 and the lower limitation of product 𝑗 in plant 𝑖 are also listed in Table 3.

Table 1 

Parameters for every stone plant.

Table 2 

Parameters of each product.

Table 3 

Parameters for product 𝑗 produced by plant 𝑖.

Taking all the numerical values into (2.21) and setting the initial temperature 𝑇₀=500, the last temperature is 0 and the cooling method is 1 decrement once. The neighbourhood can be developed as 𝑌¹_𝑖=𝑌⁰_𝑖+𝑟ℎ and 𝑋¹_𝑖𝑗=𝑋⁰_𝑖𝑗+𝑟ℎ, where 𝑟 is a random number in (−1,1) and ℎ is the step length (here ℎ=2.0). After a simulation of many cycles, the Pareto optimal solution and theobjective value are determined as shown in Tables 4 and 5. The results illustrate that although some plants have the highest productive efficiency, their high emission coefficient will result in the low exploiting quotas such as Kai Quan, Guo Jian, and De Sheng. On the other hand, stone plants will tend to produce the high value-added but low emission products due to the environmental pressure and the limitation of exploiting quotas, such as nano calcium carbonates and man-made composite slabs. However, stone plants will abundantly produce the traditional products such as granite slabs because of the huge cost of those new products.

Table 4 

Assignment results for different products.

Table 5 

Objectives for both the upper and lower levels.

4.2. Sensitivity Analysis

In fact, the decision maker is able to adjust the parameter to obtain different level solutions. From theoretical deduction, it is known that the possibilistic level is a key factor impacting the results. If the accuracy of 𝛿^𝑈_𝑖 and 𝛿^𝐿_𝑖 decreases, the feasible set is expanded and then a better Pareto optimal solution and a better Pareto optimal point are determined. From Table 5, it can be seen that the emissions increase and the economic profit and the employment decrease as the possibilistic level 𝛿^𝑈_𝑖(𝑖=1,2,3) decreases indicating that the government requirements are less strict which results in the stone plants pursuing the economic profit and neglecting the emissions and the employment objectives. Finally, the total emissions increase and the government tax revenue decreases. On the other hand, if the possibilistic level 𝛿^𝑈_𝑖(𝑖=1,2,3) increases, the government requirements are more strict and hence the total emissions decrease and the government tax revenue increases.

Similarly, for the following level, if the possibilistic levels 𝛿^𝐿_𝑖(𝑖=1,2,…,10) decrease, the plants pay less attention to the stone dust and waste water emissions resulting in an increase in profit and consequently more emission.

4.3. Comparison Analysis

For the proposed case, all the emission coefficients including 𝐸𝑑_𝑖, 𝑒𝑑_𝑖𝑗, and 𝑒𝑤_𝑖𝑗 are fuzzificated as triangular fuzzy numbers according to the real-life situation. Because all the equations in the model are linear, it actually can be easily converted into a crisp model without uncertain coefficients by the possibility measure. Lemma A.1 is given to show the process in the Appendix section, and we can get the crisp model according to (A.10). Taking all the numerical values into (A.10) and setting the same parameters for ISA, we can easily get the optimal solutions. The error analysis and computation time are listed in Table 6 and Table 7, respectively. It is obvious that the results from solving the crisp equivalent model are close to the results from simulating the model. It shows that the the fuzzy simulation technique is reasonable and efficient to solve the model for bi-level multiobjective optimization problems with fuzzy coefficients. At the same time, it is found from Table 7 that the average computational time by ISA is less than the time by FS-ISA. It is also reasonable because the process of fuzzy simulation for possibilistic constraint will spend much time to get the approximate value. However, not all possibilistic constraints can be directly converted into crisp ones. Lemma A.1 is efficient only for the special membership functions such as the triangular and trapezoidal fuzzy numbers.

Table 6 

Errors analysis by solving crisp equivalent model.

Table 7 

Computing time and memory by ISA and GA.

To illustrate that FS-ISA is suitable for this kind of fuzzy bi-level model, the results are compared with a genetic algorithm (GA). GA is one of the most popular algorithms. Many scholars also made the comparison between SA and GA in solving bi-level optimization problems [12, 42, 43]. They regard that different data scales will result in huge differences on the computational efficiency. To ensure the fairness, we also design the GA based on the fuzzy simulation for the bi-level multiobjective optimization with fuzzy coefficients. We set the chromosome number 20, the crossover rate 0.6, the mutation rate 0.8, and the iterative number 500. The average computing time and memory are listed in Table 7. Experiments show that the similar optimal results can be obtained by both FS-ISA and FS-GA, but the computational efficiency is different when the number of stone resources and stone plants changes. It is found that when the number of stone resources and stone plants is small, FS-ISA is more efficient than GA in solving the bi-level multiobjective optimization and much more computational effort is needed for FS-GA to achieve the same optimal solution as FS-ISA. However, when the data scale is large, FS-GA can reach a more optimal solution at the expense of more computation time. The result is in accordance with the findings by Xu et al. [43]. Of course, if the fuzzy bi-level multi-objective optimization model can be easily converted into the crisp model, we can obtain a more accurate solution and spend less time by ISA than that by FS-ISA.

5. Conclusions

In this paper, we have developed a bi-level multi-objective optimization model with possibilistic constraints under the fuzzy environment. In the model, the government is considered as the leader level for minimizing the emissions of the stone dust and the waste water and maximizing the employment and economic growth, and then stone plants are considered as the follower level for maximizing the profit and minimizing the emissions. Then we propose an algorithm FS-ISA to solve the model. Finally, a practical case proves that the proposed model and algorithm are efficient.

Although the model proposed in this paper should be helpful for solving some real-world problems, it only dealt with by the possibilistic constraints. If DM has different purposes such as maximizing the possibility that the predetermined goals are achieved, we can apply dependent-chance constraint to deal with it. In further research to be undertaken, a detailed analysis will be given.

Appendix

Lemma A.1. Assume that 𝐸𝑑_𝑖, 𝑒𝑑_𝑖𝑗, and 𝑒𝑤_𝑖𝑗 (𝑖=1,2,…,𝑚;𝑗=1,2,…,𝑛) are L-R fuzzy numbers with the following membership functions: 𝜇Ed_𝑖(𝑡)=𝐿𝐸𝑑_𝑖−𝑡𝛼^Ed_𝑖,𝑡<𝐸𝑑_𝑖,𝛼_𝑖^𝐸𝑑𝑅>0,𝑡−𝐸𝑑_𝑖𝛽_𝑖^𝐸𝑑,𝑡≥𝐸𝑑_𝑖,𝛽_𝑖^𝐸𝑑𝜇>0,𝑒𝑑_𝑖𝑗(𝑡)=𝐿𝑒𝑑_𝑖𝑗−𝑡𝛼^𝑒𝑑_𝑖𝑗,𝑡<𝑒𝑑_𝑖𝑗,𝛼^𝑒𝑑_𝑖𝑗𝑅>0,𝑡−𝑒𝑑_𝑖𝑗𝛽^𝑒𝑑_𝑖𝑗,𝑡≥𝑒𝑑_𝑖𝑗,𝛽^𝑒𝑑_𝑖𝑗𝜇>0,𝑒𝑤_𝑖𝑗(𝑡)=𝐿𝑒𝑤_𝑖𝑗−𝑡𝛼^𝑒𝑤_𝑖𝑗,𝑡<𝑒𝑤_𝑖𝑗,𝛼^𝑒𝑤_𝑖𝑗𝑅>0𝑡−𝑒𝑤_𝑖𝑗𝛽^𝑒𝑤_𝑖𝑗,𝑡≥𝑒𝑤_𝑖𝑗,𝛽^𝑒𝑤_𝑖𝑗>0, (A.3) where 𝛼_𝑖^𝐸𝑑, 𝛽_𝑖^𝐸𝑑 are positive numbers expressing the left and right spreads of 𝐸𝑑, 𝛼^𝑒𝑑_𝑖𝑗, 𝛽^𝑒𝑑_𝑖𝑗 are positive numbers expressing the left and right spreads of  𝑒𝑑, and 𝛼^𝑒𝑤_𝑖𝑗, 𝛽^𝑒𝑤_𝑖𝑗 are positive numbers expressing the left and right spreads of  𝑒𝑤, 𝑖=1,2,…,𝑚, 𝑗=1,2,…,𝑛. Reference functions 𝐿,𝑅∶[0,1]→[0,1] with 𝐿(1)=𝑅(1)=0 and 𝐿(0)=𝑅(0)=1 are nonincreasing, continuous functions. Then one has Pos{∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1(𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗)≤𝐹₁}≥𝛿^𝑈₁ if and only if 𝐹₁≥^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿^𝑈₁^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗. (A.4)

Proof. Let 𝜔∈[0,1] be any positive real number and 𝐿((𝐸𝑑_𝑖−𝑥)/𝛼_𝑖^𝐸𝑑)=𝐿((𝑒𝑑_𝑖𝑗−𝑦)/𝛼^𝑒𝑑_𝑖𝑗) =𝐿((𝑒𝑤_𝑖𝑗−𝑧)/𝛼^𝑒𝑤_𝑖𝑗)= 𝜔, then from (A.3) we have 𝑥=𝐸𝑑_𝑖−𝛼_𝑖^𝐸𝑑𝐿⁻¹(𝜔),𝑦=𝑒𝑑_𝑖𝑗−𝛼^𝑒𝑑_𝑖𝑗𝐿⁻¹(𝜔),𝑧=𝑒𝑤_𝑖𝑗−𝛼^𝑒𝑤_𝑖𝑗𝐿⁻¹(𝜔). (A.5) For any 𝑌_𝑗,𝑋_𝑖𝑗≥0(𝑖=1,2,…,𝑚;𝑗=1,2,…,𝑛), it easily follows that 𝑡=^𝑚_𝑖=1𝑥𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑦𝑋_𝑖𝑗+𝑧𝑋_𝑖𝑗=^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗𝐿⁻¹(𝜔). (A.6) Therefore, we have 𝐿∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝑡∑^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗=𝜔. (A.7) It is also proved by similar technique that 𝑅𝑡−∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗∑^𝑚_𝑖=1𝛽_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝛽^𝑒𝑑_𝑖𝑗+𝛽^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗=𝜔. (A.8) Hence, it is easily found that ∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1(𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗) is also a L-R fuzzy number with the left spread ∑^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1(𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_{𝑖𝑗𝑐𝑒})𝑋_𝑖𝑗 and the right spread ∑^𝑚_𝑖=1𝛽_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1(𝛽^𝑒𝑑_𝑖𝑗+𝛽^𝑒𝑤_𝑖𝑗)𝑋_𝑖𝑗. According to the definition of possibility measure proposed by Dubois and Prade [28], it can be obtained as follows: Pos^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐹₁≥𝛿^𝑈₁∑⟺𝐿^𝑚𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐹₁∑^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≥𝛿^𝑈₁⟺∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐹₁∑^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚_𝑖=1∑^𝑛_𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐿⁻¹𝛿^𝑈₁⟺^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿^𝑈₁^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+^𝑚𝑛_𝑖=1𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐹₁. (A.9) This completes the proof.

From Lemma A.1, the model (2.21) is equivalent to the following bi-level multi-objective programming problem: min𝐹^∗₁=∑^𝑚_𝑖=1𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿∑^𝑈₁^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗max𝐹₂=∑^𝑚∑_𝑖=1^𝑛_𝑗=1𝑝_𝑖𝑗𝑋_𝑖𝑗+𝑃_𝑖max𝐹₃=∑^𝑚_𝑖=1𝑆∑_𝑖^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗s.t.^𝑚_𝑖∑₌₁𝐸𝑑_𝑖𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿∑^𝑈₂^𝑚_𝑖=1𝛼_𝑖^𝐸𝑑𝑌_𝑖+∑^𝑚𝑛_𝑖=1∑_𝑗=1𝛼^𝑒𝑑_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝐷^𝑈_𝑚^𝑛∑_𝑖=1∑_𝑗=1𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿∑^𝑈₃^𝑚𝑛_𝑖=1∑_𝑗=1𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗≤𝐸𝑊^𝑈_𝑚𝑖∑₌₁𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝐷^𝐿_𝑗∀𝑗max𝐻¹_𝑖=∑^𝑛_𝑗=1𝑐_𝑗𝜃_𝑖𝑗𝑋_𝑖𝑗−∑^𝑛_𝑗=1𝑓𝑋_𝑖𝑗−ℎ_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗min𝐻_𝑖^2∗=∑^𝑛_𝑗=1𝑒𝑑_𝑖𝑗+𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗−𝐿⁻¹𝛿∑^𝑈₃^𝑛_𝑗=1𝛼^𝑒𝑑_𝑖𝑗+𝛼^𝑒𝑤_𝑖𝑗𝑋_𝑖𝑗s.t.^𝑛_𝑗∑₌₁𝑋_𝑖𝑗≤𝑌_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗≤𝐼𝑉^𝑈_𝑖𝑛𝑗∑₌₁𝑓𝑋_𝑖𝑗+ℎ_𝑖𝑌_𝑖−∑^𝑛_𝑗=1𝑋_𝑖𝑗≤𝑃𝐶^𝑈_𝑖𝜃_𝑖𝑗𝑋_𝑖𝑗≥𝑃^𝐿_𝑖𝑗. (A.10)

Abbreviations

Indices

𝑖∶	Index of stone-material plants, 𝑖=1,2,…,𝑚
𝑗∶	Index of stone products, 𝑗=1,2,…,𝑛.

Parameters

𝐸𝑑𝑖∶	Stone dust emissions coefficient when plant 𝑖 exploits
𝑒𝑑𝑖𝑗∶	Stone dust emissions coefficient when plant 𝑖 produces product 𝑗
𝑒𝑤𝑖𝑗∶	Waste water emissions coefficient when that plant 𝑖 produces product 𝑗
𝑝𝑖𝑗∶	Employment coefficient that plant 𝑖 produces product 𝑗
𝑃𝑖∶	Basic employment that plant 𝑖 needs
𝑆𝑖∶	Unit tax rate that plant 𝑖 pays to the government
𝑐𝑗∶	Unit price of product 𝑗
𝑡𝑖𝑗∶	Unit variable cost when plant 𝑖 produces product 𝑗
ℎ𝑖∶	Unit cost when plant 𝑖 holds remnant stone materials
𝜃𝑖𝑗∶	Transformation rate when plant 𝑖 produces product 𝑗
𝐶𝑖𝑗∶	Constant cost if plant 𝑖 produces product 𝑗
𝑅𝑈∶	Total stone resources upper limitation in the region
𝐷𝐿𝑗∶	Lower limitation of product 𝑗 demand
𝐸𝐷𝑈∶	Stone dust total emissions upper limitation in the region
𝐸𝑊𝑈∶	Waste water total emissions upper limitation in the region
𝐼𝑉𝑈𝑖∶	Inventory upper limitation for plant 𝑖
𝑃𝐶𝑈𝑖∶	Production cost upper limitation for plant 𝑖
𝑃𝐿𝑖𝑗∶	Lower limitation for product 𝑗 in plant 𝑖.

Decision variables

𝑌𝑖∶	Amount that the government allows plant 𝑖 to exploit
𝑋𝑖𝑗∶	Amount that plant 𝑖 uses to produce the product 𝑗.

Acknowledgments

The work is supported by the Key Program of National Natural Science Foundation of China (Grant no. 70831005) and also supported by the “985” Program of Sichuan University (Innovative Research Base for Economic Development and Management).