A dual-interval vertex analysis method and its application to environmental decision making under uncertainty

European Journal of Operational Research 200 (2010) 536–550

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Decision Support

A dual-interval vertex analysis method and its application to environmental decision making under uncertainty Y.P. Li a,*, G.H. Huang b,c,1, P. Guo c,2, Z.F. Yang d,3, S.L. Nie e,4 a

College of Urban and Environmental Sciences, Peking University, Beijing 100871, China Chinese Research Academy of Environmental Science, Beijing Normal University, Beijing 100012-100875, China c Faculty of Engineering, University of Regina, Regina, Saskatchewan, Canada S4S 0A2 d State Key Laboratory of Water Environment Simulation, School of Environment, Beijing Normal University, Beijing 100875, China e College of Mechanical Engineering and Applied Electronics Technology, Beijing University of Technology, Beijing 100022, China b

a r t i c l e

i n f o

Article history: Received 31 August 2007 Accepted 14 January 2009 Available online 21 January 2009 Keywords: Air quality Decision making Dual interval Environment Fuzzy programming Optimization Vertex analysis Uncertainty

a b s t r a c t In this study, a dual-interval vertex analysis (DIVA) method is developed, through incorporating the vertex method within an interval-parameter programming framework. The developed DIVA method can tackle uncertainties presented as dual intervals that exist in the objective function and the left- and right-hand sides of the modeling constraints. An interactive algorithm and a vertex analysis approach are proposed for solving the DIVA model. Solutions under an associated a-cut level can be generated by solving a series of deterministic submodels. They can help quantify relationships between the objective function value and the membership grade, which is meaningful for supporting in-depth analyses of tradeoffs between environmental and economic objectives as well as those between system optimality and reliability. A management problem in terms of regional air pollution control is studied to illustrate applicability of the proposed approach. The results indicate that useful solutions for planning the air quality management practices have been generated. They can help decision makers to identify desired pollution-abatement strategies with minimized costs and maximized environmental efficiencies. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Previously, a number of optimization methods were developed to deal with uncertainties in various environmental decision making problems; they included fuzzy, stochastic and interval mathematical programming (abbreviated as FMP, SMP and IMP) (Ellis et al., 1985, 1986; Guldman, 1988; Huang et al., 1992; Macchiato et al., 1994; Teng and Tzeng, 1994; Lejano et al., 1997; Chen and Huang, 2001; Zolezzi, 2001; Cifarelli et al., 2002; Haurie et al., 2004; Maqsood et al., 2005; Li et al., 2006a,b; Li and Huang, 2007; Qin et al., 2007). The IMP methods are effective in tackling uncertainties expressed as interval values with known lower and upper bounds but unknown distribution functions; however, they are incapable of dealing with uncertainties expressed as possibilistic and probabilistic distributions (Huang et al., 1992). Although SMP can deal with various probabilistic uncertainties, the increased data requirements for specifying the parameters’ probability distributions can affect their practical applicability. In comparison, FMP can deal with vagueness and ambiguity based on fuzzy set theory, where uncertainties are handled in a direct way without a large number of realizations (Li et al., 2006c). In general, the FMP methods are classified into three categories in view of the forms of uncertainties: (i) fuzzy flexible programming (FFP), (ii) fuzzy possibilistic programming (FPP), and (iii) robust programming (RP) (Inuiguchi and Tanino, 2000). In detail, FFP can deal with decision problems under fuzzy goal and constraints; however, it has difficulties in tackling ambiguous coefficients of the objective function and constraints (Dubois and Prade, 1988). RP improves upon FFP by allowing fuzzy parameters in the constraints to be represented by possibility distributions (El-Ghaoui and Lebert, 1997; El-Ghaoui et al., 1998; Inuiguchi and Sakawa, 1998; Dubois and Prade, 1999a,b; Ben-Tal and Nemirovski, 1999, 2002; Bertsimas and Sim, 2004; Nie et al., 2007; Li et al., 2008). However, the main limitations * Corresponding author. Tel./fax: +86 10 62756137. E-mail addresses: [email protected] (Y.P. Li), [email protected] (G.H. Huang), [email protected] (P. Guo), [email protected] (Z.F. Yang), [email protected] (S.L. Nie). 1 Tel.: +1 306 585 4095; fax: +1 306 585 4855. 2 Tel.: +1 306 585 4958; fax: +1 306 585 4855. 3 Tel.: +86 10 58807596; fax: +86 10 58800397. 4 Tel.: +86 10 62767681; fax: +86 10 62767681. 0377-2217/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2009.01.013

537

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

of this method remain in its difficulties in tackling uncertainties in a non-fuzzy decision space. In FPP, fuzzy parameters that are regarded as possibility distributions are introduced into the modeling frameworks. FPP can handle ambiguous coefficients in the left- and right-hand sides of the constraints and in the objective function. Previously, a number of FPP methods were developed (Dubois and Prade, 1980, 1986; Dong and Shah, 1987; Luhandjula et al., 1992; Tanaka and Guo, 1999; Dou et al., 1997; Chen et al., 1998; Inuiguchi and Tanino, 2000; Dubois et al., 2001; Iskander, 2005, 2006). Among them, the vertex method proposed by Dong and Shah (1987) is a useful tool for dealing with ambiguous coefficients based on fuzzy set theory and a-cut analysis. In fact, in many real-world problems, various uncertainties may exist with varied presentation formats in the system components. For example, uncertainties may be estimated as interval values; at the same time, the lower and upper bounds of these intervals may also be fuzzy in nature, leading to dual intervals. However, both IMP and FPP have difficulties in addressing such a complexity. One potential approach for accounting for such dual uncertainties is to link the vertex method with the IMP, leading to a dual-interval vertex analysis (DIVA) approach. Therefore, the objective of this study is to develop such a DIVA method to tackle the uncertainties presented as dual intervals that exist in the objective function and the left- and right-hand sides of the constraints. Then, a management problem in terms of regional air pollution control will be studied to illustrate applicability of the proposed approach. The results will help decision makers to identify desired pollution-mitigation strategies with minimized treatment cost and maximized environmental efficiency. Willingness to pay a higher cost may guarantee meeting environmental objectives; conversely, a desire to acquire a lower cost could run into a risk of violating environmental criteria. 2. Methodology 2.1. Dual-interval vertex analysis method Firstly, consider a deterministic linear programming (LP) problem as follows:

Min f ¼

n X

c j xj

ð1aÞ

j¼1

subject to : n X arj xj 6 br ;

r ¼ 1; 2; . . . ; s;

ð1bÞ

j¼1

xj P 0;

8j;

ð1cÞ

where cj, arj and br form sets of deterministic parameters. However, in many real-world problems, the quality of information that can be obtained is mostly not satisfactory enough to be known with precision. Consequently, when uncertainties exist as intervals with known lower and upper bounds, the above problem can be converted into an interval-parameter linear programming (ILP) model as follows:

Min f ¼

n X

cj xj

ð2aÞ

j¼1

subject to : n X r ; rj xj 6 b a

r ¼ 1; 2; . . . ; s;

ð2bÞ

j¼1

xj P 0;

8j;

ð2cÞ

r form sets of interval values with deterministic lower and upper bounds. Obviously, the ILP method can tackle uncertainrj and b where cj , a ties expressed as intervals with known lower and upper bounds (Huang et al., 1992). The detailed definitions related to the ILP method is presented in Appendix A to this paper. However, in many real-world problems, the lower and upper bounds of some interval parameters can rarely be acquired as deterministic values (Luciano and Peccati, 2001; Huang and Chang, 2003). Instead, they may often be given as subjective information that can only ~U are lower and upper bounds of dual-interval param~L and b be expressed as fuzzy sets; this leads to dual uncertainties (Fig. 1). In Fig. 1, b  ~ and can be quantified as membership functions. Such a dual interval cannot be addressed through the conventional ILP method. eter b, Fuzzy vertex analysis is effective for dealing with ambiguous and vague information (Dong and Shah, 1987; Zimmermann, 1995). A number

Membership Grade

1

~

bL

~

bU

Fig. 1. Interval number with fuzzy boundaries.

~

b

538

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

of definitions related to fuzzy vertex analysis will be introduced in Appendix B to this paper. As a result, fuzzy vertex method will be introduced into the ILP framework to handle such complexities; this leads to a dual-interval vertex analysis (DIVA) model as follows:

 X  ~  Min ~f ¼ ðcj þ dj Þxj n

ð3aÞ

j¼1

subject to : n X r ; rj xj 6 b a

r ¼ 1; 2; . . . ; s;

ð3bÞ

t ¼ s þ 1; s þ 2; . . . ; m;

ð3cÞ

j¼1 n X

~ ~ x 6 b a t; tj j

j¼1

8j;

xj P 0;

ð3dÞ

   r form sets of single intervals with deterministic lower and upper bounds; a ~t and d ~j form sets of dual intervals with fuzzy rj and b ~tj , b where cj , a lower and upper bounds; r is the number of single-interval constraints; t is the number of dual-interval constraints; m is the number of total constraints. Obviously, in model (3), dual uncertainties exist in coefficients of the objective function as well as left- and right-hand sides of the constraints. Fig. 2 provides a dual interval with two triangular membership functions for its lower and upper bounds. Assume that there is no intersection between the fuzzy sets at the two bounds. This is due to satisfy the definition of an interval value that its lower-bound should not be larger than its upper-bound (see Appendix A). Secondly, interval numbers are used to express uncertainties without distribution information. If the fuzzy sets of a dual-interval’s lower and upper bounds intersect, then the so-called ‘‘interval” is actually described by fuzzy membership functions, such that the interval representation becomes unnecessary (Nie et al., 2007). Thirdly, if the fuzzy sets of lower and upper bounds intersect, the interactive algorithm for solving the interval-parameter linear programming problem cannot be used for solv ~L ; b ~U  ¼ ½½bL ; b L ; ½bU ; b U , where bL is the lower boundary of the interval’s lower-bound; b L is the ~ t ¼ ½b ing such a DIVA model. Therefore, let b t t t t t t t t U is the upper boundary of the upper bound. Then, the upper boundary of the lower-bound; bUt is the lower boundary of the upper bound; b t foregoing problem [i.e., model (3)] can be reformulated as follows:

 X   X L L U x Min ~f ¼ ½½dj ; dj ; ½dUj ; d cj xj þ j j n

n

j¼1

j¼1

subject to : n X r ; rj xj 6 b a

ð4aÞ

r ¼ 1; 2; . . . ; s;

ð4bÞ

j¼1 n X

L ; ½bU ; b U ; Ltj ; ½aUtj ; a Utj xj 6 ½½bLt ; b ½½aLtj ; a t t t

t ¼ s þ 1; s þ 2; . . . ; m;

ð4cÞ

j¼1

8j:

xj P 0;

ð4dÞ

A two-step solution method is proposed for solving model (4). In the first step, based on the interactive algorithm, the DIVA model can be transformed into two sets of submodels for each membership grade (a-cut). The submodels corresponding to ~f L can be first formulated  ~t > 0, and  ~f > 0Þ: r > 0, b when the objective function is to be minimized. Thus, we have (assume that b

Min ~f L ¼

k1 n X X L xL þ L xU ½cLj þ ½dLj ; d ½cLj þ ½dLj ; d j j j j j¼1

ð5aÞ

j¼k1 þ1

subject to : k1 X

jarj jU SignðaUrj ÞxLj þ

j¼1 k1 X

n X

U

jarj jL SignðaLrj ÞxUj 6 br ;

r ¼ 1; 2; . . . ; s;

ð5bÞ

j¼k1 þ1

tj jU Signða Utj ÞxLj ½jatj jU SignðaUtj Þ; ja

j¼1

þ

n X

U ; tj jL Signða Ltj ÞxUj 6 ½bUt ; b ½jatj jL SignðaLtj Þ; ja t

t ¼ s þ 1; s þ 2; . . . ; m;

ð5cÞ

j¼k1 þ1

xLj P 0;

j ¼ 1; 2; . . . ; k1 ;

xUj

j ¼ k1 þ 1; k1 þ 2; . . . ; n;

P 0;

xLj

ð5dÞ ð5eÞ xUj

where (j = 1, 2, . . . , k1) are variables with positive coefficients in the objective function; (j = k1 + 1, k1 + 2, . . . , n) are variables with negL , ative coefficients in the objective function. Based on vertex analysis method, take one end point from each of the fuzzy intervals (i.e., ½dLj ; d j U Þ; then, the obtained end points can be combined into an n-array, leading to 2n combinations for n fuzzy sets. Through Utj , and ½bUt ; b ½aUtj ; a t solving 2n problems, a set of lower-bound objective function values ðf1L ; f2L ; . . . ; f2Ln Þ can be obtained. The optimized lower-bound of the interval for the objective function value (under an a-cut level) can be identified as follows:

½f Lopt ; f Lopt a ¼ ½minðf1L ; f2L ; . . . ; f2Ln Þ; maxðf1L ; f2L ; . . . ; f2Ln Þa :

ð6Þ

539

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

~

μ Ã (bt ) 1 α3 α2 α1 0

b tL

bt L

b Ut

bt U

~

bt

Fig. 2. Fuzzy interval with triangular membership functions.

In the second step, another set of submodels corresponding to ~f U can be formulated: k1 n X X U xU þ U xL ½cUj þ ½dUj ; d ½cUj þ ½dUj ; d j j j j

Min ~f U ¼

j¼1

ð7aÞ

j¼k1 þ1

subject to : k1 X

jarj jL SignðaLrj ÞxUj þ

j¼1 k1 X

n X

tj jL Signða Ltj ÞxUj þ ½jatj jL SignðaLtj Þ; ja

j¼1

r ¼ 1; 2; . . . ; s;

ð7bÞ

06

n X

L ; tj jU Signða Utj ÞxLj 6 ½bLt ; b ½jatj jU SignðaUtj Þ; ja t

t ¼ s þ 1; s þ 2; . . . ; m;

ð7cÞ

j¼k1 þ1

xUj P xLjopt ;

xLjopt

L

jarj jU SignðaUrj ÞxLj 6 br ;

j¼k1 þ1

xLj

6

j ¼ 1; 2; . . . ; k1 ;

xUjopt ;

ð7dÞ

j ¼ k1 þ 1; k1 þ 2; . . . ; n;

ð7eÞ

xUjopt

where (j = 1, 2, . . . , k1) and (j = k1 + 1, k1 + 2, . . . , n) are solutions corresponding to the upper boundary of lower interval objective value U , ½aL ; a L Þ, leading to 2n combinations for n fuzzy sets Ltj , and ½bLt ; b (i.e., f Lopt Þ. Similarly, take one end point from each of the fuzzy sets (i.e., ½dUj ; d t j tj under each a-cut level. Through solving 2n deterministic problems, a set of upper-bound objective function values ðf1U ; f2U ; . . . ; f2Un Þ can be obtained. The optimized upper-bound interval for the objective function value (under an a-cut level) can be identified as follows:

½f Uopt ; f Uopt a ¼ ½minðf1U ; f2U ; . . . ; f2Un Þ; maxðf1U ; f2U ; . . . ; f2Un Þa ;

ð8Þ

f U is the upper where is lower boundary of the interval’s upper bound (for the objective function value) with ¼ opt U U U U  boundary of the upper-bound with f opt ¼ maxðf1 ; f2 ; . . . ; f2n Þ. Then, through integrating the solutions of the two sets of submodels, the solution of dual interval for the objective function value (under an a-cut level) can be obtained as follows: f Uopt

U fopt

minðf1U ; f2U ; . . . ; f2Un Þ;

 ð~f opt Þa ¼ ½½f Lopt ; f Lopt a ; ½f Uopt ; f Uopt a :

ð9Þ

Additional dual-interval solutions based on the inputs of the fuzzy sets can be obtained through repeating the process with changed a-cut levels. Then, the final dual-interval solutions for the objective function value under various a-cut levels can be obtained. Obviously, in DIVA, uncertain parameters in the constraints’ left- and right-hand sides and the objective function can be expressed as dual intervals, such that the robustness of optimization effort can be enhanced (Here, ‘‘robustness” means the developed DIVA possesses enhanced capacities in reflecting complexities of system uncertainties). The detailed solution processes for the DIVA model can be described as follows: Step 1: Formulate the DIVA model. Step 2: Transform the DIVA model into two submodels, where the lower bound submodel corresponding to ~f L is firstly desired when the  objective is to minimize ~f . Step 3: Discretize the range of membership grade [0, 1] into a finite number of a-cut levels. For each membership grade aj, find the corresponding intervals for fuzzy sets. Step 4: Take one end point from each of the intervals (under an a-cut level), there are 2n combinations for n fuzzy sets. Step 5: Solve 2n ~f L submodels, a set of lower-bound objective function values ðf1L ; f2L ; . . . ; f2Ln Þ can be obtained. Step 6: The desired lower-bound of the interval for the objective function value can be identified as: ½f Lopt ; f Lopt a ¼ ½minðf1L ; f2L ; . . . ; f2Ln Þ; maxðf1L ; f2L ; . . . ; f2Ln Þa . Step 7: Formulate and solve the ~f U submodel (under the same a-cut level) and obtain a set of values ðf1U ; f2U ; . . . ; f2Un Þ. Step 8: The desired upper-bound of the interval for the objective function value can be identified as: ½f Uopt ; f Uopt a ¼ ½minðf1U ; f2U ; . . . ; f2Un Þ; maxðf1U ; f2U ; . . . ; f2Un Þa . Step 9: Integrate the solutions of the two sets of submodels, the solution of dual interval for the objective function value (under an a-cut level) can be obtained as

 ð~f opt Þa ¼ ½½f Lopt ; f Lopt a ; ½f Uopt ; f Uopt a : Step 10: Repeat the processes (Step 4 to Step 9) at every a-cut level; then we can obtain the final dual-interval solutions.

540

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

2.2. Illustrative example The following illustrative example can be formulated to demonstrate the applicability of the proposed method. First, consider an ILP model as follows (Huang et al., 1998):

Max f ¼ ½26; 30x1  ½5:5; 6:0x2 subject to : ½8; 10x1  ½12; 14x2 6 ½3:8; 4:2; ½2:4; 2:8x1 þ ½3; 4x2 6 ½6:0; 6:5; x1 ; x2 P 0:

ð10aÞ ð10bÞ ð10cÞ ð10dÞ

According to an interactive algorithm, the above model can be converted into the following two submodels (Huang et al., 1998):

Max f U ¼ 30xU1  5:5xL2

ð11aÞ

subject to : 8xU1  14xL2 6 4:2;

ð11bÞ

2:4xU1 xU1 ; xL2

ð11dÞ

þ

4xL2

6 6:5;

P0

ð11cÞ

and

Max f L ¼ 26xL1  6:0xU2

ð12aÞ

subject to : 10xL1  12xU2 6 3:8;

ð12bÞ

2:8xU1 þ 3xL2 6 6:0;

ð12cÞ

06 xU2 xU1opt

xL1

P

6

xU1opt ;

ð12dÞ

xL2opt ;

ð12eÞ

xL2opt

where and are solutions from submodel (11). The solutions for the ILP model are x1opt ¼ ½1:32; 1:62, x2opt ¼ ½0:64; 0:73, and f opt ¼ ½29:94; 45:08. When dual intervals exist in model (10), a DIVA model can be formulated as follows:

 Max ~f ¼ ½½26; 27; ½29; 30x1  ½5:5; 6:0x2 subject to : ½8; 10x1  ½12; 14x2 6 ½3:8; 4:2; ½½2:4; 2:5; ½2:7; 2:8x1 þ ½3; 4x2 6 ½½6:0; 6:2; ½6:3; 6:5; x1 ; x2 P 0:

ð13aÞ ð13bÞ ð13cÞ ð13dÞ ð13eÞ

The DIVA model can be solved through a two-step solution method. First, based on the interactive algorithm, model (13) can be transferred into the following two submodels that correspond to the lower and upper bounds of objective function value.

Max ~f U ¼ ½29; 30xU1  5:5xL2

ð14aÞ

subject to : 8xU1  14xL2 6 4:2;

ð14bÞ

½2:4; 2:5xU1 xU1 ; xL2 P 0

ð14dÞ

þ

4xL2

6 ½6:3; 6:5;

ð14cÞ

and

Max ~f L ¼ ½26; 27xL1  6:0xU2

ð15aÞ

subject to : 10xL1  12xU2 6 3:8;

ð15bÞ

½2:7; 2:8xL1 þ 3xU2 0 6 xL1 6 xU1opt ;

ð15dÞ

xU2 P xL2opt :

6 ½6:0; 6:2;

ð15cÞ ð15eÞ

Secondly, according to the vertex analysis method, under each a-cut level, take one end point from each of the fuzzy intervals; then, each submodel (under an a-cut level) can be converted into a number of deterministic problems. Table 1 presents the solutions for submodels (14) and (15) under a = 0.5. The results indicate that solutions for ~f U are (43.053, 43.493, 43.602, 43.846, 44.047, 44.294, 44.405, 44.859), and ½f Uopt ; f Uopt 0:5 ¼ ½43:053; 44:859; solutions for ~f L are (30.099, 30.369, 30.502, 30.763, 30.776, 31.039, 31.175, 31.455), and  ½f Lopt ; f L opt 0:5 ¼ ½30:099; 31:455. Thus, we have ð~f opt Þ0:5 ¼ ½½30:099; 31:455; ½43:053; 44:859 and, correspondingly, x0:5 1opt ¼ ð1:327; 1:345; 0:5 ¼ ð0:606; 0:628; 0:789; 0:816Þ. Additional solutions based on the input fuzzy sets can be obtained through repeating 1:586; 1:624Þ and x2opt the process with changed a-cut levels. Then, the final solutions for the objective function value and decision variables under multiple a-cut levels can be obtained.

541

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550 Table 1 Solutions for submodels (14) and (15) when a = 0.5. ~f U xU1 xL2

43.052 1.586 0.606

43.493 1.603 0.616

43.602 1.607 0.618

43.846 1.586 0.606

44.047 1.624 0.628

44.294 1.603 0.616

44.405 1.607 0.618

44.859 1.624 0.628

~f L xL1 xU2

30.099 1.327 0.789

30.369 1.340 0.780

30.502 1.346 0.805

30.763 1.327 0.789

30.776 1.359 0.816

31.039 1.340 0.800

31.175 1.346 0.805

31.455 1.359 0.816

3. Application to air quality management 3.1. Overview of the study system For decades, issues of air pollution control have been of substantial concerns since the increasing pollutant concentrations in the ambient environment have caused adverse effects on crops, trees, materials and human health (Li et al., 2006a).The impacts of air pollution are being observed across many regions and socio-economic sectors. For example, air pollution may be linked to a sequence of events such as pollutant generation, emission, minimization, mitigation and transport, as well as impacts on various receptors (Flagan and Seinfeld, 1988). Furthermore, pollution-mitigation processes are related to many factors such as production scale, plant size and geographical location (Wark et al., 1998). These processes and the related factors are associated with a variety of uncertainties. Therefore, development of effective optimization methodologies, which can tackle various uncertainties and complexities and thus provide bases for supporting air quality management and pollution-control planning, is desired. An air pollution control system can be characterized by one or several sources that generate adverse impacts on a given receptor. In North America, over 65% of SO2 (or more than 13 million tons per year) released to the atmosphere comes from electric utilities that burn coal (USEPA, 2004). Since it is economically infeasible or technically impossible to design processes with zero emission of air pollutants, it is desired that authorities seek to control the emissions to levels at which the effects are minimized. In the study system, two power plants (i.e., coal-fired units) are considered as major SO2 emission sources within a planning horizon of three 5-year periods. The residential zone is adversely affected by the emitted SO2 from the two sources (Fig. 3). The problem under consideration is how to minimize the total cost for pollution abatement while satisfying the overall environmental regulation. The Gaussian dispersion model can be employed to predict SO2 concentrations at different locations of the receptor zone under different meteorological conditions. If a ‘‘well-behaved” emission plume is assumed, air pollutant could have a Gaussian distribution about the plume centerline, which is located a distance H above the ground. For a steady rate of emission, the ground-level pollutant concentration at arbitrary downwind location (u, w) can be estimated as follows (Kumar, 1977; Haith, 1982; Nevers, 2000):

Cðu; wÞ ¼

q

puru rz

" exp

w2 2r2u

! 

H2 2r2z

!# ð16Þ

;

where C(u, w) is ground-level pollutant concentration (mg/m3) at point (u,w); u is downwind distance from source (m); w is crosswind dis is mean wind velocity (m/ tance from source (m); H represents average effective stack height (m); q denotes pollutant emission rate (mg/s); u s); rw and rz represent standard deviations of the plume at w and z directions (m), respectively. The rw and rz levels depend on wind speed, solar heating, and local turbulence, and they are functions of both meteorological conditions and downwind distance. Correspondingly, to use above equation for calculating the air-pollutant concentrations, meteorological conditions must be known. Thus, the centerline ground-level concentration can be given as follows (Kumar, 1977):

Cðu; wÞ ¼

q

purw rz

exp 

! H2 : 2r2z

ð17Þ

Farmland zone 1 Power plant 1 Residential zone

Steel mill Oil refinery

Power plant 2

N

Farmland zone 2 Wind direction Fig. 3. The study system.

W

E S

542

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

Assume that rw and rz may be approximated by (Kumar, 1977)

rz ¼ cub andrw ¼ eud ;

ð18Þ

where c, b, e and d are constants. Thus, the maximum ground-level concentration can be represented as follows (Kumar, 1977; Kumar and Djurfors, 1978):

C max

  q exp  bþd 2b ¼ : h 2 ibþd 2b puce c2bH ðbþdÞ

ð19Þ

Consequently, a transfer factor (representing the worst case of maximum ground-level concentration) can be defined as follows:

  exp  bþd 2b nip ¼ ; h 2 ibþd 2b puce c2bH ðbþdÞ

ð20Þ

where nip represents the contribution of each source (i) to the ground concentration in each receptor zone (p). More extensive discussions of air-pollutant dispersion models are provided in Kumar (1977), Haith (1982), Wark et al. (1998) and Nevers (2000). Therefore, systems analysis techniques could be used for managing the regional air quality in more efficient and environmentally benign ways, which may be helpful for generating a desired compromise between required environmental objectives and minimized system costs. The objective is to minimize the cost for SO2 abatement; the constraints include the SO2 emission standards, pollution-control capacities, and allowable pollutant-loading levels. Thus, we have

Min f ¼ Lk

2 X 2 X 3 X i¼1

j¼1 k¼1

X ijk P Sik ;

8i; k;

Djk X ijk

ð21aÞ

subject to : 2 X

ð21bÞ

j¼1

½constraints of SO2 -mitigation requirement X ijk 6 FC ;ijk

8i; j; k;

ð21cÞ

½constraints of SO2 -mitigation capacity 2 X

ð1  gj ÞX ijk 6 e;ik

8i; k;

ð21dÞ

j¼1

½constraints of SO2 -emission allowance 2 X 2 X i¼1

nip ð1  gj ÞX ijk 6 h;k

8k;

ð21eÞ

j¼1

½constraints of ambient air quality requirement X ijk P 0;

8i; j; k;

ð21fÞ

½non-negativity constraints where hk denotes the total allowable SO2-loading level at the two emission sources during period k (lg/m3); i is the name of SO2 emission source; j is the type of SO2 control measure; k means the planning period; Lk denotes the length of period k (day); gj is the efficiency of control measure j; Djk is the operating cost of control measure j for SO2 mitigation during period k ($/tonne); eik denotes the SO2-emission allowance for source i during period k (tonne/day); FCijk is the SO2-mitigation capacity of measure j at source i in period k (tonne/day); Sik denotes the amount of SO2 generated from source i during period k (tonne/day); nip is the transfer factor of pollution contribution from emission source i to receptor zone p (day/m3); Xijk denotes the amount of SO2 (generated from source i) to be mitigated by control measure j during period k (tonne/day). Uncertainties exist in a variety of impact factors and pollution-related processes, such as pollutant characteristics, emission rates, mitigation measures, and pollution impacts. These uncertainties would inevitably affect the efforts in modeling pollutant transport in the atmosphere, which is important for relating source conditions to ambient air quality (Li et al., 2006a). Many system components such as the related cost coefficients, the SO2-treatment capacities, the SO2-emission allowance, and the allowable SO2-loading level may not be available as deterministic values; they may be presented in terms of interval values. For example, the cost of SAS measure for SO2-treatment is approximately 60.5 to 76.7 dollars per tonne; power plant 1 has a capacity to treat 260–300 tonnes of SO2 per day; SO2-emission allowance for power plant 1 is in the range of 23.5 and 28.9 tonnes per day. Moreover, the coals used by the electricity-generation facilities have various sulfur contents (low-sulfur coals have less than 1% sulfur, while high-sulfur ones contain 2–5% sulfur); this may lead to a varying SO2-generation rate. Furthermore, efficiencies of the control measures to mitigate SO2 emissions may also vary with operating conditions (e.g., reagent ratio, temperature, and inlet SO2 concentration). These uncertainties can be estimated as interval values; at the same time, the lower and upper bounds of these intervals may also be fuzzy in nature, leading to dual uncertainties. Due to the existence of dual

543

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

uncertainties expressed as fuzzy sets and interval values, the developed DIVA method is considered suitable for tackling the above complexities. Therefore, model (21) can be reformulated as follows:

Min f ¼ Lk

2 X 2 X 3 X i¼1

 jk X  ijk D

ð22aÞ

j¼1 k¼1

subject to : 2 X

  ijk P e S ik ; X

8i; k;

ð22bÞ

j¼1

 ijk 6 FC ijk ; X 2 X

8i; j; k;

~ j ÞX  ijk 6 eik ; ð1  g

ð22cÞ

8i; k;

ð22dÞ

j¼1 2 X 2 X i¼1

~ j ÞX  ijk 6 hk ; nip ð1  g

8k;

ð22eÞ

j¼1

 ijk P 0; X

8i; j; k;

ð22fÞ

  ~ j Þ are presented as dual intervals with their bounds being prewhere the SO2-generation rates ðe S ik Þ and the pollution-control measures ðg sented as triangular membership functions; this results in multiple levels of interactive relationships. Model (22) can be solved based on the solution method as described in Section 2. Table 2 shows the related technical and economic data including SO2-generation rates, SO2-emission allowances, allowable SO2-loading levels, and SO2-mitigation costs. The SO2-generation rates vary with the type of coal and the related combustion conditions (low-sulfur coals have less than 1% sulfur, while high-sulfur ones contain 2–5% sulfur). Moreover, from a long-term planning point of view, economic development within the region can result in an increased demand for electricity, leading to a continuous increase of SO2 generation from each source. To meet the environmental requirement, each power plant has to adopt control measures to mitigate the SO2 emission and to avoid penalties from the government. Table 3 shows the efficiencies and capacities of different control measures. The measure of wet soda ash scrubber (SAS) has a higher acid-removal efficiency than the lime spray dryer (LSD); meanwhile, this approach also has a relatively high operating cost (due to the high reagent cost). The representative cost and technical data were investigated based on governmental reports and the related literature (Flagan and Seinfeld, 1988; Wark et al., 1998; Nevers, 2000; Liu et al., 2003; USEPA, 2004; Li et al., 2006a). 3.2. Result analysis In this study, five a-cut levels (0, 0.25, 0.50, 0.75 and 1) were considered; for each a-cut level, sixteen scenarios were examined based on different combinations of the resulting dual intervals. Fig. 4 presents the solutions of total cost from the DIVA model under a-cut levels of 0, 0.25, 0.5 and 0.75. The results indicate that different combinative considerations on the uncertain inputs would lead to varied objective function values. Here, the lower-bound cost corresponds to advantageous conditions, while the upper-bound one is associated with more demanding conditions. For example, when a = 0.25, the lower-bound costs would be 119.1, 122.7, 124.1, 126.9, 127.4, 127.6, 130.2, 130.8, 132.4, 132.5, 135.6, 135.7, 136.3, 139.2, 141.9 and 144.6 million dollars. Among them, the minimum and maximum values would respec-

Table 2 Technical and economic data. Planning period k=1

k=2

k=3

~sik (tonne/day) SO2-generation rate,  Power plant 1 Power plant 2

[[230, 250], [320, 350]] [[138, 150], [215, 230]]

[[246, 266], [342, 372]] [[150, 162], [235, 250]]

[[265, 285], [370, 400]] [[165, 177], [254, 270]]

eik (tonne/day) SO2-emission allowance,  Power plant 1 Power plant 2

[36.3, 47.6] [26.0, 33.2]

[35.4, 45.9] [25.1, 32.0]

[34.5, 44.5] [24.2, 30.6]

 jk Regular cost for treating allowable SO2 emission ($/tonne), D Soda ash scrubber (SAS) [60.5, 76.7] Lime spray dryer (LSD) [32.0, 38.4]

[70.0, 84.5] [35.0, 42.3]

[81.0, 97.2] [40.0, 48.6]

Allowable SO2-loading level,  hk (lg/m3)

[18.7, 23.1]

[17.8, 21.8]

[19.6, 24.2]

Table 3 SO2 treatment efficiencies and capacities. Soda ash scrubber (SAS)

Lime spray dryer (LSD)

[[0.85, 0.88], [0.95, 0.99]]

[[0.72, 0.74], [0.83, 0.85]]

[260, 300] [170, 200]

[160, 185] [150, 170]

 ~ j (%) SO2-treatment efficiency, g SO2-treatment capacity, FC ijk (tonne/day) Power plant 1 Power plant 2

544

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

160

(a) Lower bound System cost ($10 6)

150

α=0

α = 0.25

α = 0.5

α = 0.75

140 130 120 110 1

16

Scenario

240

System cost ($10 6)

(b) Upper bound α=0

α = 0.25

α = 0.5

α = 0.75

230

220

210 1

Scenario

16

Fig. 4. Total costs under different a-cut levels.

tively be 119.1 and 144.6 million dollars; they would form the lower and upper boundaries for the lower-bound cost (i.e., ~f L 0:25 ¼ $½119:1; 144:6  106 ). Similarly, under a = 0.25, the upper-bound cost would be 213.9, 213.9, 214.4, 217.5, 217.6, 218.4, 220.2, opt 220.9, 220.9, 221.4, 223.9, 224.58, 225.0, 225.4, 227.2 and 231.5 million dollars. Among them, the minimum and maximum values would be 213.9 and 231.5 million dollars; correspondingly, the lower and upper boundaries for the upper-bound cost would respectively be 6 $213.9  106 and $231.5  106 (i.e., ~f U0:25 opt ¼ $½213:9; 231:5  10 ). Through integration of the lower- and upper-bound values, the total cost 6 would be $[[119.1, 144.6], [213.9, 231.5]]  10 under a = 0.25, which is a dual interval. Similarly, when a = 0, the minimum and maximum values for the lower-bound cost would be $115.3  106 and $149.3  106, while the minimum and maximum values for the upper-bound cost would be $216.7  106 and $237.1  106; they constitute the total cost of  ~f 0 ¼ $½½115:3; 149:3; ½216:7; 237:1  106 . When a = 0.5, the minimum and maximum values for the lower-bound cost would be opt $123.2  106 and $140.2  106, while the minimum and maximum values of upper-bound cost would be $212.4  106 and $226.0  106; correspondingly, the total cost under a = 0.5 would be $[[123.2, 140.2], [212.4, 226.0]]  106. When a = 0.75, the minimum and maximum values for the lower-bound cost would be $127.2  106 and $135.7  106, while the minimum and maximum values of upper-bound cost would be $212.7  106 and $220.1  106; correspondingly, the total cost under a = 0.75 would be $[[127.2, 135.7], [212.7, 220.1]]  106. The total cost would be [131.5, 215.6]  106 when a = 1, presenting in a single interval; this is because the triangular membership functions for fuzzy sets are used in this study. Each a-cut level denotes a subset of elements that belong to a fuzzy set at least to a membership degree of a (also named degree of plausibility). The results indicate that different a-cut levels correspond to different costs. Under a lower degree of plausibility (i.e., lower a-cut level), the intervals are wider; conversely, a higher degree of plausibility would lead to narrower intervals. For example, when a = 0 (with the lowest plausibility degree), the total cost would be in the range of $115.3  106 to 237.1  106; in comparison, under a = 1 (with the highest plausibility degree), the total cost would be $131.5  106 to 215.6  106 (i.e., a much narrower interval). The results also indicate that the marginal cost would be decreased as the a-cut level is raised (i.e., the marginal variation of total cost with a-cut level). For example, an average cost [i.e., ðf Lopt þ f Lopt þ f Uopt þ f Uopt Þ=4] could be calculated under each a-cut level. Correspondingly, the average costs would be $179.6  106, $177.3  106, $175.5  106, $173.9  106 and $173.6  106 under a-cut levels of 0, 0.25, 0.5, 0.75 and 1, respectively. Letting the average cost under a = 0 be a reference one; thus, the variations of average cost (Dfavg) would be $2.3  106 under Da = 0.25, $4.1  106 under Da = 0.5, $5.7  106 under Da = 0.75, and $6.0  106 under Da = 1. The marginal variation of cost with a-cut level (i.e., Dfavg/Da) would be 9.2 when a = 0.25, 8.2 when a = 0.5, 7.6 when a = 0.75, and 6.0 when a = 1, demonstrating that the marginal cost would be decreased as the a-cut level is raised. Table 4 presents the solutions for mitigating SO2 emissions under a-cut levels of 0, 0.25, 0.50 and 0.75. For power plant 1 during period 1, the amounts (of SO2 emission) treated by soda ash scrubber (SAS) would be 96.3, 160.3, 163.8 and 186.3 tonne/day (i.e., X aijk , X bijk , X cijk and X dijk Þ when a = 0.25; there are four mitigation options, where the first option (i.e., 96.3 tonne/day) corresponds to the lower boundary of the interval’s lower-bound ðf Lopt Þ, and 160.3, 163.8 and 186.3 tonne/day correspond to f Lopt , f Uopt and f Uopt , respectively. Besides, for power plant 1, the amounts of SO2 treated by lime spray dryer (LSD) would be 136.2, 87.2, 160.0 and 160.0 tonne/day, respectively corresponding to the

545

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550 Table 4 Solutions of treated SO2 amounts under different a-cut levels. Treated SO2 amount ðX aijk ; X bijk ; X cijk ; X dijk Þ (tonne/day)

Power plant

Period

Control measure

a=0

a = 0.25

a = 0.50

a = 0.75

1 1 2 2

1 1 1 1

SAS LSD SAS LSD

(87.1, 172.3, 172.3, 193.3) (142.9, 77.7, 147.7, 156.7) (19.1, 67.7, 67.7, 109.2) (118.9, 82.3, 147.3, 120.8)

(96.3, 160.3, 163.8, 186.3) (136.2, 87.2, 160.0, 160.0) (24.2, 60.7, 66.9, 99.2) (115.3, 87.8, 150.0, 129.0)

(106.7, 149.2, 167.5, 182.5) (128.3, 95.8, 160.0, 160.0) (30.3, 54.5, 68.8, 94.6) (110.7, 92.5, 150.0, 134.6)

(116.5, 137.8, 171.3, 178.8) (121.0, 104.7, 160.0, 160.0) (35.7, 47.8, 70.6, 81.5) (106.8, 97.7, 150.0, 142.9)

1 1 2 2

2 2 2 2

SAS LSD SAS LSD

(129.0, 219.8, 219.8, 232.0) (117.0, 46.2, 122.2, 140.0) (50.0, 102.8, 102.8, 145.0) (100.0, 59.2, 132.2, 105.0)

(138.9, 207.2, 207.2, 213.9) (109.6, 56.3, 138.6, 154.3) (55.6, 95.3, 95.3, 133.9) (95.9, 65.2, 141.6, 114.2)

(149.9, 195.3, 195.3, 204.5) (101.1, 65.7, 154.2, 160.0) (62.1, 88.5, 88.8, 124.3) (90.9, 70.5, 150.0, 122.0)

(160.5, 183.2, 193.3, 200.8) (93.0, 75.3, 160.0, 160.0) (68.2, 81.4, 96.0, 114.1) (86.3, 76.1, 144.6, 130.3)

1 1 2 2

3 3 3 3

SAS LSD SAS LSD

(174.3, 271.5, 271.5, 271.7) (90.7, 13.5, 98.5, 128.3) (94.6, 153.1, 153.1, 166.1) (70.4, 23.9, 100.9, 103.9)

(185.1, 258.1, 258.1, 259.6) (82.6, 24.4, 115.6, 136.6) (101.1, 145.0, 145.0, 168.5) (65.4, 30.5, 111.0, 99.5)

(196.7, 245.3, 245.3, 245.3) (73.3, 34.7, 132.2, 147.2) (108.1, 137.3, 137.3, 157.5) (59.9, 36.7, 120.7, 108.5)

(208.2, 232.5, 232.5, 232.5) (64.3, 45.0, 148.8, 156.3) (114.9, 129.6, 129.6, 146.2) (54.6, 42.9, 130.4, 117.8)

Note: X aijk , X bijk , X cijk and X dijk are amounts of SO2 (generated from source i) to be treated by control measure j during period k (tonne/day); they correspond to f Lopt , f Lopt , f Uopt and f U , respectively. opt

Lower-bound amount of SO 2 treated by SAS

Treated SO2 amount (tonne/day)

Upper-bound amount of SO 2 treated by SAS Lower-bound amount of SO 2 treated by LSD Upper-bound amount of SO 2 treated by LSD

250 200 150 100 50 0

Plant 1

Plant 2

Period 1

Plant 1

Plant 2

Period 2

Plant 1

Plant 2

Period 3

Fig. 5. Amounts of SO2 treated by SAS and LSD under a = 1.

objective function values of f Lopt , f Lopt , f Uopt and f Uopt . For power plant 2 during period 1, the amounts of SO2 treated by SAS and LSD would be (24.2, 60.7, 66.9, 99.2) and (115.3, 87.8, 150.0, 129.0) tonne/day when a = 0.25, respectively. The total mitigated SO2 is the sum of SO2 amounts treated by SAS and LSD. The SAS has a higher efficiency and a higher cost in treating SO2 than the LSD. Solutions with more SAS involvement would correspond to a higher cost, which could guarantee that air quality criteria be met; however, a decision aiming towards less use of SAS would result in a lower cost but a higher risk of violating environmental requirements. The results for other acut levels and time periods can be similarly interpreted based on the solutions presented in Table 4. Fig. 5 shows the amounts of SO2 (generated from the two power plants) treated by SAS and LSD under a = 1. The results are expressed as   ~ j are deterministic values when a = 1. The amounts of SO2 treated by SAS single intervals since the boundaries of fuzzy intervals e S ik and g would keep rising along with the increasing SO2-generation rates and environmental requirements. This is because the relevant pollutionabatement efficiencies need to be continuously improved to satisfy the increasing requirements of ambient air quality. Consequently, the results provide useful bases for generating desired technology combinations that would lead to satisfied environmental quality as well as least-possible pollution-abatement costs. 4. Comparison with the conventional ILP methods The study problem could also be solved through an interval-parameter linear programming (ILP) method by simplifying the fuzzy lower  and upper bounds into deterministic values (Huang et al., 1992). This implies that, in ILP, SO2-generation rates ðe S ik Þ and pollution-control  ~ j Þ are considered to be single intervals with deterministic lower and upper bounds. For example, in DIVA, the amount of SO2 measures ðg generation (for power plant 1 in period 1) is [[230, 250], [320, 350]] tonne/day with their two bounds being triangular membership functions; in comparison, in ILP, a single interval of [230, 350] tonne/day is simply adopted as SO2-generation rate for power plant 1 in period 1. Table 5 presents the solutions obtained from the ILP model. It is indicated that the ILP solutions are significantly different from the DIVA ones. The ILP solutions also lack confidence-degree information as defined by membership grade (a) in the DIVA solutions. The total cost obtained from ILP would be $[129.9, 223.6]  106 and the corresponding uncertainty degree (UD) be 26.5%. The concept of UD is useful for quantitatively evaluating levels of uncertainty for input parameters (Huang and Moore, 1993). For DIVA, let the solutions of the total costs (i.e., dual intervals) under different a-cut levels be equal to their mid values. For example, when a = 0, the average lower- and upper-bound costs would be $132.3  106 [i.e., $(115.3 + 149.3)  106/2] and $226.9  106 [i.e., $(216.7 + 237.1)  106/2], respectively; the

546

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

Table 5 Solutions obtained from the ILP method. Treated SO2 (tonne/day)

Power plant i

Control measure j

Period k

Solution

X 111 X 121 X 211 X 221 X 112 X 122 X 212 X 222 X 113 X 123 X 213 X 223 f ($106)

1 1 2 2 1 1 2 2 1 1 2 2

SAS LSD SAS LSD SAS LSD SAS LSD SAS LSD SAS LSD

1 1 1 1 2 2 2 2 3 3 3 3

[129.0, 190.0] [100.8, 160.0] [41.8, 80.0] [96.2, 150.0] [176.8, 212.0] [69.2, 160.0] [76.9, 100.0] [73.1, 150.0] [228.5, 240.0] [36.5, 160.0] 127.3 [37.7, 142.7] [129.9, 223.6]

corresponding UD level would be 26.3% under a = 0. Similarly, the average costs would be [131.8, 222.7]  106 with a UD level of 25.6% when a = 0.25, $[131.7, 219.2]  106 with UD = 24.9% when a = 0.5, and $[131.5, 216.4]  106 with UD = 24.2% when a = 0.75. When a = 1, the total cost would be $[131.5, 215.6]  106 and the corresponding UD level be 24.2%. The results indicate that the total cost from ILP has a higher UD level than that from DIVA under each a-cut level. The lower level of uncertainty in the DIVA solution demonstrates that this proposed method can more effectively specify the variety of uncertainties based on the a-cut analysis technique. On the other hand, DIVA can directly incorporate more uncertainties expressed as intervals and fuzzy sets within its optimization framework; in comparison, ILP is a simplified approach in dealing with the dual uncertainties. Moreover, Nie et al. (2007) proposed an interval-parameter fuzzy robust programming (IFRP) method. The IFRP incorporated robust programming within an ILP framework, where uncertain parameters could be represented as interval numbers and/or fuzzy membership functions. However, IFRP had difficulties in dealing with fuzzy information in the objective function. Moreover, in its solution process, the IFRP model was converted into a deterministic version through transforming m imprecise constraints into 2 km precise inclusive ones that correspond to k a-cut levels; thus the derived solution was a set of single intervals. This approach led to a lack of possibility information in the obtained solutions (as defined by a-cut levels). In comparison, DIVA could handle dual uncertainties expressed as both fuzzy sets and intervals that exist in the constraints’ left- and right-hand sides and the objective function. Dual-interval solutions with an associated a-cut level could be generated by solving two 2n deterministic submodels; besides, the information of possibility (as represented by a-cut levels) for the objective function value would be available in the DIVA solutions. Therefore, compared with the conventional ILP methods, DIVA is more effective in reflecting system uncertainties. 5. Conclusions In this study, a dual-interval vertex analysis (DIVA) method has been developed through incorporating the vertex method within an interval-parameter programming framework. The DIVA can handle uncertainties presented as dual intervals that exist in the objective function and the left- and right-hand sides of the modeling constraints. Based on an interactive algorithm and a vertex analysis approach, solutions under an associated a-cut level can be generated by solving two 2n deterministic submodels. The related fuzzy information can also be robustly reflected in the solutions for the objective function value and decision variables. The solutions can help quantify relationships between the cost and the membership grade, which is meaningful for supporting in-depth analyses of tradeoffs between environmental and economic objectives as well as those between system optimality and reliability. A management problem in terms of regional air pollution control has been studied to illustrate applicability of the proposed approach. Results of the case study indicate that useful solutions for planning the air quality management practices have been generated. They can help decision makers to identify desired pollution-mitigation strategies with minimized total costs and maximized environmental efficiencies. Willingness to pay a higher cost may guarantee meeting environmental objectives; conversely, a desire to acquire a lower cost could run into a risk of violating environmental criteria. The major limitation of the DIVA method is its capability of handling large-scale problems. This is because two 2n deterministic submodels under each a-cut level have to be solved for generating optimal solution when n fuzzy sets exist. Particularly, when n is large, a large number of possibilistic combinations will have to be examined, resulting in a relatively high computational requirement. This will lead to difficulties in the method’s practical application. Therefore, one potential extension of this research is about the development of more effective algorithm for solving the large-scale DIVA problems. Moreover, DIVA can tackle uncertainties presented as intervals and fuzzy sets; however, it has difficulties in dealing with uncertainties expressed as random variables. For example, the SO2-generation rate may vary with coal types at power plants, which could be expressed as a random variable (with a probability distribution). Stochastic programming with recourse is effective for problems where an analysis of policy scenarios is desired and coefficients are random with known probability distributions. Therefore, DIVA can be integrated with stochastic programming methods to enhance its capacities in dealing with uncertainties expressed in multiple formats. Acknowledgements This research was supported by the Major State Basic Research Development Program of MOST (2005CB724201, 2005CB724207 and 2006CB403303), the special fund of 08ESPCT-K, and the Natural Sciences and Engineering Research Council of Canada. The authors are very grateful to the editors and the anonymous reviewers for their insightful comments and suggestions.

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

547

Appendix A. Interval-parameter programming Let x denote a closed and bounded set of real numbers. An interval number with known lower and upper bounds but unknown distribution can be defined as follows:

x ¼ ½xL ; xU  ¼ ft 2 xjxL 6 t 6 xU g;

ða1Þ

where xL and xU represent the lower and upper bounds of  x, respectively. When xL = xU, x becomes a deterministic number, i.e., x ¼ xL ¼ xU . For x, its Signð xÞ can be defined as follows:

SignðxÞ ¼



ðx P 0Þ; 1 ðx < 0Þ;

1

ða2Þ

where x P 0 if xU P 0 and xL P 0;

ða3Þ

x < 0 if xU < 0 and xL < 0:

ða4Þ

Therefore, when interval number  x is greater than zero, its sign is 1; when  x is less than zero, its sign is 1. For x, its absolute value jxj can be defined as follows:

jxj ¼



x

if x P 0; x if x < 0;

ða5Þ

where

jxjL ¼



xL xU

if x P 0; if x < 0

ða6Þ

if x P 0; if x < 0:

ða7Þ

and

jxjU ¼



xU xL

Let R denote a set of interval numbers. An interval row-vector X is a tuple of interval numbers

X 2 fRg1n ;

X ¼ fxi ¼ ½xLi ; xUi j8ig;

ða8Þ

where lower/upper bounds of interval vector X can be expressed as X L ¼ fxLi j8ig and X U ¼ fxUi j8ig: Then, for an interval vectors (or matrix), we have

8i; j;

X P 0 iff xij P 0;

8i; j;

X 6 0 iff xij 6 0;

X 2 fRgmn ; X 2 fRgmn ;

m P 1;

ða9Þ

m P 1:

ða10Þ

, we have Let * 2 {+, , , } be a binary operation on interval numbers. For two interval numbers x and y

x  y  ¼ ½minfx  yg; maxfx  yg;

xL 6 x 6 xU ;

y L 6 y 6 yU :

ða11Þ

 does not contain a zero. Hence, we have In case of division, it is assumed that y

x þ y  ¼ ½xL þ yL ; xU þ yU ; L

U

U

L

x  y  ¼ ½x  y ; x  y ; x  y  ¼ ½minfx  yg; maxfx  yg; x  y  ¼ ½minfx  yg; maxfx  yg:

ða12Þ ða13Þ ða14Þ ða15Þ

Uncertain degree (UD) of an interval number ð xÞ can be defined as follows (Huang and Moore, 1993):

UD ðxÞ ¼

x U  xL  100%: x U þ xL

ða16Þ

When  x is highly uncertain, UD ð xÞ would become close to 100%; when  x turns towards deterministic, UD ð xÞ would move towards 0. An interval-parameter linear programming (ILP) model can be defined as follows (Huang et al., 1992):

Min f ¼ CX

ða17Þ

subject to :  AX 6 B;

ða18Þ

X P 0;

ða19Þ

mn

 2 fRgm1 , C 2 fRg1n and X 2 fRgn1 ðR denote a set of interval values). An interactive solution algorithm was proposed to where A 2 fRg ,B solve the above problem through analyses of the interrelationships between the parameters and the variables and between the objective function and the constraints. According to Huang et al. (1992), when the objective is to be minimized, a submodel corresponding to fL  > 0, and f > 0Þ: can be firstly formulated as follows (assume that b i

548

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

Min f L ¼

k1 X

n X

cLj xLj þ

j¼1

cLj xUj

ða20Þ

j¼k1 þ1

subject to : k1 X

n X

jaij jU SignðaUij ÞxLj þ

j¼1

U

jaij jL SignðaLij ÞxUj 6 bi ;

8i;

ða21Þ

j¼k1 þ1

xLj P 0;

j ¼ 1; 2; . . . ; k1 ;

ða22Þ

xUj

j ¼ k1 þ 1; k1 þ 2; . . . ; n;

ða23Þ

P 0;

where xLj (j = 1, 2, . . . , k1) are lower bounds of interval variables with positive coefficients in the objective function; xUj (j = k1 + 1, k1 + 2, . . . , n) L can be are upper bounds of interval variables with negative coefficients. Solutions of xLjopt (j = 1, 2, . . . , k1), xUjopt (j = k1 + 1, k1 + 2, . . . , n), and fopt U obtained from submodel (a20–a23). Thus, the submodel corresponding to f can be formulated as follows: k1 X

Min f U ¼

n X

cUj xUj þ

j¼1

cUj xLj

ða24Þ

j¼k1 þ1

subject to : k1 X

jaij jL SignðaLij ÞxUj þ

j¼1

L

jaij jU SignðaUij ÞxLj 6 bi ;

8i;

ða25Þ

j¼k1 þ1

xUj P xLjopt ; 06

n X

xLj

6

j ¼ 1; 2; . . . ; k1 ;

xUjopt ;

j ¼ k1 þ 1; k1 þ 2; . . . ; n:

ða26Þ ða27Þ

Then, through integration of the solutions from the two submodels, solutions for model (a17–a19) presented as sets of intervals can be obtained:

xjopt ¼ ½xLjopt ; xUjopt ; f opt ¼ ½f L ; f U : opt opt

8j;

ða28Þ ða29Þ

Generally, ILP method has the following advantages: (i) it allows uncertainties to be directly communicated into the optimization process, (ii) it does not lead to more complicated intermediate models, and thus has a relatively low computational requirement, and (iii) it does not require distributional information for model parameters, which is particularly meaningful for practical applications because it is typically much more difficult for planners/engineers to specify distributions than to define fluctuation intervals (Huang et al., 1992). Appendix B. Fuzzy vertex analysis method Fuzzy sets are defined as sets whose members are vague objects. The general notation of a fuzzy set can be presented as follows (Zimmermann, 1995):

e AðxÞ ¼ fx; leðxÞjx 2 X and le ðxÞ 2 ½0; 1g; A A

ðb1Þ

e where X = {x} is a universe set of elements; AðxÞ is a fuzzy set of X; le ðxÞ is the membership function or grade of membership. The le ðxÞ value A A ranges from 0 to 1, where 1 represents full membership and 0 denotes non-membership. The closer le ðxÞ is to 1, the more likely it is that an A e e e element x belongs to A; conversely, the closer le ðxÞ is to 0, the less likely it is that x belongs to A. Moreover, a fuzzy set A is normal if and only A e is convex if a < c < b for every set of real numbers (a, b, if _x le ðxÞ ¼ 1, which means that at least one x exists such that leðxÞ ¼ 1. Fuzzy set A A A c 2 X), then le ðcÞ P min½leðaÞ; le ðbÞ. This means that the membership function consists of increasing and decreasing parts. A A A e at least to the degree a, and this degree is also called the degree of An a-cut is defined as the set of elements that belong to fuzzy set A confidence or the degree of plausibility. The a-cut can be described as follows (Zimmermann, 1995):

e a ¼ fxjl ðxÞ P a; x 2 Xg; A e A

ðb2Þ

e a is the a-cut level of A, e and it consists of all components of X whose membership grade is greater than or equal to a. The support of where A e is defined by the classical set as follows: fuzzy set A

e ¼ fxjl ðxÞ > 0g: suppð AÞ e A

ðb3Þ

e may conveniently be exThe convexity condition ensures that the support is in an interval. The membership function of any fuzzy set A pressed for all x 2 X in canonical form (Dubois and Prade, 1986; Klir, 1997)

8 f ðxÞ > > > eA > <1 leA ðxÞ ¼ > ge ðxÞ > > > : A 0

when x 2 ½a; bÞ; when x 2 ½b; c; when x 2 ðc; d;

ðb4Þ

otherwise;

where a, b, c, d 2 X and a 6 b 6 c 6 d, fe is a real-valued function that is increasing and right-continuous, and ge is a real-valued function that A A is decreasing and left-continuous. There are many approaches for generating membership functions. However, the triangular membership function is employed by many researchers for the sake of tractability, simplicity, and unimodality (Civanlar and Trussel, 1986; Dou et al.,

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

549

1997; Li, 2003). A triangular membership function can be easily defined by specifying three numbers: (1) the most credible value, (2) the lowest possible value, and (3) the highest possible value. The most credible value is assigned a membership grade of 1, and any number that falls short of the lowest possible value or exceeds the highest possible value will get a membership grade of 0. The intermediate membership grades can be obtained by linear interpolation. The algebraic operations on real numbers can be extended to fuzzy sets by the use of extension principle (Zadeh, 1978; Dubois and e1, A e2, . . . , A e n be n fuzzy subsets Prade, 1980; Zimmermann, 1995). Let X be a Cartesian product of universe, X = X1  X2      Xn, and A e e e of X1, X2, . . . , Xn, respectively. The Cartesian product of A 1 , A 2 , . . . , A n is defined as follows (Dubois and Prade, 1980; Kaufmann and Gupta, 1991):

e1      A en ¼ A

Z X 1 X n

min½le ðx1 Þ; . . . ; le ðxn Þ=ðx1 ; . . . ; xn Þ: A1

ðb5Þ

An

Let f be a mapping from X1  X2      Xn to a universe Y such that y = f(x1, x2, . . . , xn), where x1 2 X1, . . . , xn 2 Xn. The extension principle can e2, . . . , A e n such that (Lai and Hwang, 1992; Zimmermann, 1995): e1, A e on Y through f from n fuzzy sets A induce a fuzzy subset B

leB ðyÞ ¼

sup x1 ;...;xn y¼f ðx1 ;...;xn Þ

fmin½le ðx1 Þ; le ðx2 Þ; . . . ; le ðxn Þg; A1

A2

ðb6Þ

An

leB ðyÞ ¼ 0 if f 1 ðyÞ ¼ ;;

ðb7Þ

where f1(y) is the inverse image of y. leðyÞ is the highest value among membership grades le e ðx1 ; x2 ; . . . ; xn Þ of the realization of y B A 1  A 2 e An using an n-array (x1, x2, . . . , xn). Application of the extension principle to fuzzy sets can be viewed as its extension to a-cuts when the membership functions are continuous (Dou et al., 1997; Li, 2003). The vertex method based on a-cut analysis is useful for dealing with fuzzy sets (Dong and Shah, 1987; Dou et al., 1997). Using the a-cut concept, each fuzzy variable characterized by a convex membership function can be converted into a group of intervals with various a levels. Then, intervals with the same a-cut level from all fuzzy variables can be processed through interval analysis, resulting in an interval function associated with an a-cut level. At an a-cut level, the interval function can be denoted as (Chen et al., 1998)

ya ¼ f ðxa1 ; xa2 ; . . . ; xan Þ; ya

½ya ; ya ; xa

ðb8Þ

a ½aa ; b ; xa

a ½aa ; b ; . . . ; xa

a ½aa ; b .

aa

a

where ¼ L R 1¼ 1 1 2¼ 2 2 n ¼ n n When i ¼ bi , the interval reduces to a point (a deterministic value). The interval analysis is equivalent to solving a minimization problem for the lower-bound and a maximization problem for the upper-bound

yaL ¼ min f ðx1 ; x2 ; . . . ; xn Þ and yaR ¼ max f ðx1 ; x2 ; . . . ; xn Þ; a

a

ðb9Þ

a

where x1 2 ½aa1 ; b1 ; x2 2 ½aa2 ; b2 ; . . . ; xn 2 ½aan ; bn . When the function is continuous, the desired value of interval function can be obtained as follows:

ya ¼ fmin½f ðcj Þ; f ðdi Þ; max½f ðcj Þ; f ðdi Þ; j ¼ 1; 2; . . . ; 2n ;

8ig;

ðb10Þ

e can be obtained based on analysis of lower where cj is the ordinate of jth vertex and di is an extreme point. Final solution to the fuzzy set B and upper bounds of the output values at different a-cut levels. References A.Ben-Tal, A.NemirovskiRobust solutions to uncertain linear programsOperational Research Letters 25, 1–13. A.Ben-Tal, A.NemirovskiRobust optimization – methodology and applicationsMathematic Program Search B 92, 453–480. D.Bertsimas, M.SimThe price of robustnessOperations Research 52 (1), 35–53. M.J.Chen, G.H.HuangA derivative algorithm for inexact quadratic program-application to environmental decision-making under uncertaintyEuropean Journal of Operational Research 128, 570–586. H.K.Chen, W.K.Hsu, W.L.ChiangA comparison of vertex method with JHE methodFuzzy Sets and Systems 95, 201–214. M.D.Cifarelli, D.Masciandaro, L.Peccati, S.Salsa, A.TaglianiSuccess or failure of a firm under different financing policies: A dynamic stochastic modelEuropean Journal of Operational Research 136 (3), 471–482. M.R.Civanlar, H.J.TrusselConstructing membership functions using statistical dataFuzzy Sets and Systems 18, 1–13. W.Dong, H.C.ShahVertex method for computing functions of fuzzy variablesFuzzy Sets and Systems 24, 65–78. C.Dou, W.Woldt, I.Bogardi, M.DahabNumerical solute transport simulation using fuzzy sets approachJournal of Contaminant Hydrology 27 (1–2), 107–126. D.Dubois, H.PradeSystems of linear fuzzy constraintsFuzzy Sets and Systems 3, 37–48. D.Dubois, H.PradeFuzzy sets and statistical dataEuropean Journal of Operational Research 25, 345–356. D.Dubois, H.PradePossibility Theory – An Approach to the Computerized Processing of Uncertainty. Plenum Press, New York. D.Dubois, H.PradeA synthetic view of belief revision with uncertain inputs in the framework of possibility theoryInternational Journal of Approximate Reasoning 17 (2–3), 295–306. D.Dubois, H.PradeOn fuzzy interpolationInternational Journal of General Systems 28 (2), 103–112. D.Dubois, H.Prade, R.SabbadinDecision-theoretic foundations of qualitative possibility theoryEuropean Journal of Operational Research 128, 459–478. L.El-Ghaoui, H.LebertRobust solutions to least-square problems to uncertain data matricesSIAM Journal of Matrix Annual Application 18, 1035–1064. L.El-Ghaoui, F.Oustry, H.LebretRobust solutions to uncertain semidefinite programsSIAM Journal on Optimization 9, 33–52. J.H.Ellis, E.A.McBean, G.J.FarquharChance-constrained stochastic linear programming model for acid rain abatement – I. Complete colinearity and noncolinearityAtmospheric Environment 19, 925–937. J.H.Ellis, E.A.McBean, G.J.FarquharChance-constrained stochastic linear programming model for acid rain abatement – II. Limited colinearityAtmospheric Environment 20, 501–511. R.C.Flagan, J.H.SeinfeldFundamentals of Air Pollution Engineering. Prentice-Hall Inc., Englewood Cliffs, New Jersey. J.-M.GuldmanChance-constrained dynamic model of air quality managementFuzzy Sets and Systems 114 (5), 1116–1126. A.D.HaithEnvironmental Systems Optimization. John & Sons, Inc., New York. A.Haurie, J.J.Kübler, A.Clappier, H.Bergh, A.van denMetamodeling approach for integrated assessment of air quality policiesEnvironmental Modeling and Assessment 9, 1–12. G.H.Huang, N.B.ChangThe perspectives of environmental informatics and systems analysisJournal of Environmental Informatics 1, 1–6. G.H.Huang, R.D.MooreGrey linear programming, its solving approach, and its applicationInternational Journal of Systems Science 24, 159–172. G.H.Huang, B.W.Baetz, G.G.PatryA grey linear programming approach for municipal solid waste management planning under uncertaintyCivil Engineering Systems 9, 319– 335. G.H.Huang, B.W.Baetz, G.G.PatryTrash-flow allocation: Planning under uncertaintyInterfaces 28 (6), 36–55. M.Inuiguchi, M.SakawaRobust optimization under softness in a fuzzy linear programming problemInternational Journal of Approximate Reasoning 18, 21–34.

550

Y.P. Li et al. / European Journal of Operational Research 200 (2010) 536–550

M.Inuiguchi, T.TaninoPortfolio selection under independent possibilistic informationFuzzy Sets and Systems 115, 83–92. M.G.IskanderA suggested approach for possibility and necessity dominance indices in stochastic fuzzy linear programmingApplied Mathematics Letters 18, 395–399. M.G.IskanderExponential membership function in stochastic fuzzy goal programmingApplied Mathematics and Computation 173, 782–791. A.Kaufmann, M.M.GuptaIntroduction a Fuzzy Arithmetic: Theory and Applications. Van Nostrand Reinhold, New York. G.J.KlirThe role of constrained fuzzy arithmetic in engineering. In: B.M.Ayyub, M.M.GuptaUncertainty Analysis in Engineering and Sciences: Fuzzy, Logic Statistics, and Neural Network Approach. Kluwer Academic Publishers, Norwell, MA, USA, pp. 1–19. Kumar, A., 1977. Pollutant Dispersion in the Planetary Boundary Layer, Ph.D. Thesis. University of Waterloo. Kumar, A., Djurfors, S.G., 1978. Maximum and critical ground level concentrations from a synthetic crude oil plant, paper #78-48.5. In: Proceedings of 1978 Annual Meeting of the Air Pollution Control Association, Houston, Texas. Y.J.Lai, C.L.HwangFuzzy Mathematical Programming Methods and Applications. Springer-Verlag, Berlin. R.P.Lejano, P.M.Ayala, E.A.GonzalesOptimizing the mix of strategies for control of vehicular emissionsEnvironmental Management 21 (1), 79–87. Li, J.B., 2003. Development of an Inexact Environmental Modeling System for the Management of Petroleum-Contaminated Sites, Ph.D. Dissertation, University of Regina, Regina, Saskatchewan, Canada. Y.P.Li, G.H.HuangFuzzy two-stage quadratic programming for planning solid waste management under uncertaintyInternational Journal of Systems Science 38 (3), 219–233. Y.P.Li, G.H.Huang, A.Veawab, X.H.Nie, L.LiuTwo-stage fuzzy-stochastic robust programming: A hybrid model for regional air quality managementJournal of the Air and Waste Management Association 56, 1070–1082. Y.P.Li, G.H.Huang, S.L.NieAn interval-parameter multistage stochastic programming model for water resources management under uncertaintyAdvances in Water Resources 29, 776–789. Y.P.Li, G.H.Huang, S.L.Nie, Y.F.HuangIFTSIP: Interval fuzzy two-stage stochastic mixed-integer programming: A case study for environmental management and planningCivil Engineering and Environmental Systems 23 (2), 73–99. Y.P.Li, G.H.Huang, X.H.Nie, S.L.NieA two-stage fuzzy robust integer programming approach for capacity planning of environmental management systemsEuropean Journal of Operational Research 189, 399–420. L.Liu, G.H.Huang, Y.Liu, G.A.FullerA fuzzy-stochastic robust programming model for regional air quality management under uncertaintyEngineering Optimization 35 (2), 177– 199. E.Luciano, L.PeccatiCycles optimization: The equivalent annuity and the NPV approachesInternational Journal of Production Economics 69 (1), 65–83. M.K.Luhandjula, H.Ichihashi, M.InuiguchiFuzzy and semi-infinite mathematical programmingInformation Science 61, 233–250. M.F.Macchiato, C.Cosmi, M.Ragosta, G.TosatoAtmospheric emission reduction and abatement costs in regional environmental planningJournal of Environmental Management 41, 141–156. I.Maqsood, G.H.Huang, J.S.YeomansAn interval-parameter fuzzy two-stage stochastic program for water resources management under uncertaintyEuropean Journal of Operational Research 167, 208–225. N.de NeversAir Pollution Control Engineering, second ed. McGraw-Hill Companies, Inc.. X.H.Nie, G.H.Huang, Y.P.Li, L.LiuIFRP: A hybrid interval-parameter fuzzy robust programming approach for municipal solid waste management planning under uncertaintyJournal of Environmental Management 84, 1–11. X.S.Qin, G.H.Huang, G.M.Zeng, A.Chakma, Y.F.HuangAn interval-parameter fuzzy nonlinear optimization model for stream water quality management under uncertaintyEuropean Journal of Operational Research 180, 1331–1357. H.Tanaka, P.J.GuoPortfolio selection based on upper and lower exponential possibility distributionsEuropean Journal of Operational Research 114, 115–126. J.Teng, G.TzengMulticriteria evaluation for strategies of improving and controlling air quality in the super city: A case study of Taipei CityJournal of Environmental Management 40, 213–229. USEPA (United States Environmental Protection Agency), SO2 Allowance Market Analysis: 2004 Update, Washington. K.Wark, C.F.Warner, W.T.DavisAir Pollution: Its Origin and Control, third ed. Harper Collins Publishers, New York. L.A.ZadehFuzzy sets as a basis for a theory of possibilityFuzzy sets and systems 1, 3–28. H.-J.ZimmermannFuzzy Set Theory and its Applications, third ed. Kluwer Academic Publishers. T.ZolezziWell-posedness and optimization under perturbationsAnnals of Operations Research 101 (1), 351–361.