Risk management for sulfur dioxide abatement under multiple uncertainties

Abstract

In this study, interval-parameter programming, two-stage stochastic programming (TSP), and conditional value-at-risk (CVaR) were incorporated into a general optimization framework, leading to an interval-parameter CVaR-based two-stage programming (ICTP) method. The ICTP method had several advantages: (i) its objective function simultaneously took expected cost and risk cost into consideration, and also used discrete random variables and discrete intervals to reflect uncertain properties; (ii) it quantitatively evaluated the right tail of distributions of random variables which could better calculate the risk of violated environmental standards; (iii) it was useful for helping decision makers to analyze the trade-offs between cost and risk; and (iv) it was effective to penalize the second-stage costs, as well as to capture the notion of risk in stochastic programming. The developed model was applied to sulfur dioxide abatement in an air quality management system. The results indicated that the ICTP method could be used for generating a series of air quality management schemes under different risk-aversion levels, for identifying desired air quality management strategies for decision makers, and for considering a proper balance between system economy and environmental quality.

Keywords

risk management / conditional value-at-risk / interval optimization / two-stage programming / uncertainty / air quality management

Introduction

Two-stage stochastic programming (TSP), as a typical stochastic programming method, can effectively tackle optimization problems where an analysis of policy scenarios is desired and the related data are mostly uncertain ( Birge and Louveaux, 1988). In applying this model, a decision is first undertaken before values of random variables are disclosed, and then, after the random events have occurred and their values are known, a recourse action is made in order to minimize “penalties” that may appear due to any infeasibilities ( Ruszczyński, 1993). The initial action is called the first-stage decision, and the corrective action is named the second-stage decision. Over the past decades, the TSP model with the expected recourse function was widely explored ( Huang and Loucks, 2000; Seifi and Hipel, 2001; Ahmed, 2004; Li et al., 2006; Cai et al., 2011a, b; Shao et al., 2011; Dong et al., 2012). However, a potential limitation of the conventional TSP model is that it can only account for the expectation as the optimization criterion without any consideration on the variability of the recourse values; hence, it is a risk neutral approach. As a result, the TSP model may become infeasible when the decision maker is risk averse under extra-high variability conditions.

A risk measurement method that could help tackle the above shortcomings is conditional value-at-risk (CVaR), which could capture the notion of risk in stochastic programming ( Dentcheva and Ruszczynski, 2006; Miller and Ruszczynski, 2008; Lim et al., 2010). This method is modified from the value-at-risk (VaR) model which has become an essential tool in financial markets when quantifying portfolio market risk ( Hsu et al., 2012). Previously, Dentcheva and Ruszczynski ( 2006) utilized a CVaR model for solving the problem of constructing a portfolio of finitely many assets . Miller and Ruszczynski ( 2008) introduced a risk-adjusted measure with CVaR to evaluate the market portfolios. These studies demonstrated that a CVaR-based model cannot only handle the expected loss under extreme conditions, but also calculate the risk through a linear model rather than nonlinear ones ( Shao et al., 2011). The idea of introducing CVaR into the objective functions in the framework of a TSP model is contained by fairly recent research topics ( Noyan, 2012). For instance, Schultz and Tiedemann ( 2006) introduced CVaR and integer programming into a TSP framework, and transformed it into an explicit mixed-integer linear programming model when the probability distribution is discrete and finite. Fabian ( 2008) proposed the decomposition frameworks for handling CVaR objectives and constraints in a TSP model. Schultz ( 2011) formulated a CVaR-based TSP method, and explored the structural properties of these optimization problems and the algorithms for the solutions. Noyan ( 2012) considered a CVaR-based TSP method, and employed it for supporting disaster management where uncertain information was addressed through a discrete scenario method. The CVaR-based TSP method can effectively help to analyze predefined policy scenarios to deal with uncertainties in stochastic formats, to reflect the expected losses from extreme events, and to balance the trade-off between system cost and risk. However, few previous studies were reported in applying inexact CVaR-based TSP methods to address multiple uncertainties that exist in air quality management systems.

In fact, in air quality management problems, various uncertainties exist in a number of system components as well as their interrelationships ( Liu et al., 2003). The random characteristics of pollutants emission (e.g., sulfur dioxide (SO₂) generation rate) caused by combustion conditions (e.g., fuel property, temperature, pressure, and air inflow), the interval characteristics of pollution control technology (e.g., treatment efficiency, and control technology capacity), the dispersion processes of pollutants (e.g., pollutant emission amount, exit velocity and temperature, wind speed and direction, and atmospheric stability), the errors in estimated modeling parameters, and the vagueness of system objectives and constraints are possible sources of uncertainties ( Byun et al., 2003; Athanasiadis and Mitkas, 2007). In general, the system objectives are often associated with a number of socio-economic and ecological factors such as economic cost, environmental protection, and ecological sustainability, while the constraints are related to pollutant discharges, capacity of pollution control technologies, environmental requirements, and policy regulations. Moreover, these uncertainties may be further exaggerated by not only interactions among various uncertain and dynamic impact factors, but also their associations with economic implications of violated environmental requirements ( Li et al., 2006). Previously, a number of inexact optimization techniques, such as interval-parameter programming, stochastic-parameter programming, and fuzzy-parameter programming were developed for dealing with uncertainties in air quality management problems ( Lejano et al., 1997; Wang and Milford, 2001; Zolezzi, 2001; Ma and Zhang, 2002; Liu et al., 2003; Haurie et al., 2004; Li et al., 2006; An and Eheart, 2007; Li et al., 2008; Lv et al., 2011). Generally, fuzzy-parameter programming can deal with imprecision or vagueness (expressed by membership functions) in decision-making problems; however, it may lead to complicated submodels that are not applicable to practical problems ( Inuiguchi and Sakawa, 1998; Liu et al., 2003; Tan et al., 2011). Although stochastic-parameter programming, including TSP and CVaR-based TSP, can directly incorporate various probabilistic uncertainties within the optimization frameworks, the increased data requirements for specifying the parameters’ probability distributions may affect its practical applicability ( Marti and Ploöchinger, 1990; Li et al., 2006). For uncertainties that cannot be quantified as membership or distribution functions, interval-parameter programming provides an alternative way to directly communicate uncertainties into the optimization process without leading to more complicated intermediate models ( Dai and Huang, 2012). Nevertheless, it may become infeasible when its right-hand-side parameters are highly uncertain ( Li et al., 2008). Moreover, the interval-parameter programming cannot reflect the trade-off between system costs and risks of violated environmental standards caused by the extra high SO₂ generated levels.

Consider an air quality management system in a region where multiple pollution-control technologies are used to reduce the air pollutants discharged by multiple emission sources. Moreover, these pollutants may disperse along the atmosphere with temporal and spatial variations, and influence the air quality of multiple receptor zones. Any change in one activity may lead to a series of environmental effects. Uncertainties exist in a number of impact factors and pollution-related processes, such as pollutant generation amounts, treatment efficiencies, emissions, and transport, as well as impacts on related environmental standards; Moreover, such uncertainties may be further multiplied by their associations with economic penalties if the promised targets are violated ( Huang and Chang, 2003; Lv et al., 2011). Furthermore, variations in SO₂ generation rates may result in extra high SO₂ generated levels even if the probability that such conditions happen is small. In fact, such extra high SO₂ generated levels would increase the risk of violated environmental standards. Therefore, incorporation of various uncertainties, economic penalties, and risks within a general mathematical programming framework is desired for supporting air quality management and planning.

As an extension of previous efforts, an interval-parameter CVaR-based two-stage programming (ICTP) method will be developed for air quality management, where techniques of interval-parameter programming, two-stage stochastic programming (TSP), and conditional value-at-risk (CVaR) will be incorporated within a general framework. The developed ICTP will be used for handling uncertainties expressed as probability distributions and interval values, and also for effectively analyzing predefined policy scenarios. Moreover, the ICTP method will help control the system risks of violated environmental standards, and avoid the serious problems related to excessive emissions of air pollutants. A case study will then be provided for demonstrating applicability of the developed methodology. The results can help decision makers to not only manage air quality problems, but also gain insight into the tradeoffs between environmental and economic objectives.

Methodology

When uncertainties of the right sides of the model are expressed as probability distribution functions, and decisions need to be made periodically over time (i.e., SO₂ generation rates are random), the problem can be formulated as a TSP model ( Huang and Loucks, 2000). A general form of the TSP model can be formulated as follows ( Birge and Louveaux, 1988):

(1a) \min \mathfrak{E}\left[f\left(\mathit{x},\omega \right)\right]= \min \left\{{\mathit{c}}^{\mathrm{T}}\mathit{x}+{\mathfrak{E}}_{\omega \in \mathrm{\Omega }}\left[Q\left(\mathit{x},\omega \right)\right]:\mathit{x}\in X\right\},

with

(1b) Q\left(\mathrm{x},\omega \right)=\min \left\{\mathit{q}{\left(\omega \right)}^{\mathrm{T}}\mathit{y}\left(\omega \right):\mathit{y}\in Y\right\},

subject to

(1c) \mathit{D}\left(\omega \right)\mathit{y}\left(\omega \right)+\mathit{T}\left(\omega \right)\mathit{x}\ge \mathit{h}\left(\omega \right),

where c\in {\mathrm{ℝ}}^{n_{\mathrm{1}}} , X\subseteq {\mathrm{ℝ}}^{n_{\mathrm{1}}} , and Y\subseteq {\mathrm{ℝ}}^{n_{\mathrm{2}}} ; \omega is a random variable from probability space (\mathrm{\Omega },F,P) with \mathrm{\Omega }\to {\mathrm{ℝ}}^k , q:\mathrm{\Omega }\to {\mathrm{ℝ}}^{n_{\mathrm{2}}} , h:\mathrm{\Omega }\to {\mathrm{ℝ}}^{m_{\mathrm{2}}} , D:\mathrm{\Omega }\to {\mathrm{ℝ}}^{m_{\mathrm{2}}\times n_{\mathrm{2}}} , and T:\mathrm{\Omega }\to {\mathrm{ℝ}}^{m_{\mathrm{2}}\times n_{\mathrm{1}}} ; \mathbb{E} denotes the expectation operates; x and y are the vectors of first-stage and second-stage decision variables, respectively. Moreover, f(x,\omega )=c^{T}x+Q(x,\omega ) is the total cost function of the first-stage problem [i.e. model (1a)] and Q(x,\omega ) , known as the recourse function, is the objective function of the second-stage problem [i.e., model (1b)]. Model (1) is risk-neutral two-stage stochastic programming model in the sense that it focuses on the optimization of an expectation objective and ignores the risk associated with the variability of random variables. Conditional value-at-risk (CVaR) as a risk measurement method can capture the notion of risk in stochastic programming. Therefore, we consider the CVaR-based TSP method as follows:

(2) \min \left\{\mathfrak{E}\left[f\left(\mathit{x},\omega \right)\right]+\lambda \mathrm{C}\mathrm{V}\mathrm{a}{\mathrm{R}}_{\beta }\left[f\left(\mathit{x},\omega \right)\right]:\mathit{x}\in \mathit{X}\right\},

where {\mathrm{CVaR}}_{\beta } denotes the conditional value-at-risk at level \beta ; \lambda is a nonnegative weight to trade-off expected cost with risk, which is specified by decision makers according to their risk preferences. Birbil et al. ( 2009) indicated that {\mathrm{CVaR}}_{\beta }(Z+a)={\mathrm{CVaR}}_{\beta }(Z)+a for a\in \mathrm{ℝ} and Z\in  . Then, we have {\mathrm{CVaR}}_{\beta }[f(x,\omega )]=c^{T}x+{\mathrm{CVaR}}_{\beta }[Q(x,\omega )] , and

(3a) \mathfrak{E}\left[f\left(x,\omega \right)\right]+\lambda {\mathrm{CVaR}}_{\beta }\left[f\left(x,\omega \right)\right]

(3b) =\left(1+\lambda \right)c^{\mathrm{T}}x+\mathfrak{E}\left[Q\left(x,\omega \right)\right]+\lambda {\mathrm{CVaR}}_{\beta }\left[Q\left(x,\omega \right)\right] .

Suppose the random vector \omega ( \omega \in {\mathbf{ℝ}}^m ) has a probability density function p(\cdot ) , and Q(x,\omega ) is continuous in x and measurable in \omega where for each x\in X . Given the decision variable x , the probability of Q(x,\omega ) not exceeding a threshold of \alpha is represented as \psi (x,\alpha )={\displaystyle {\int }_{Q(x,\omega )\le \alpha }p(\omega )d\omega } . Given a confidence level \beta and a fixed x , the value-at-risk (VaR) is defined by {\mathrm{VaR}}_{\beta }(x)=\min \{\alpha \in \mathbf{ℝ}:\psi (x,\alpha )\ge \beta \} . And then, {\mathrm{CVaR}}_{\beta }(x) defined as the expected value of the loss that exceeds {\mathrm{VaR}}_{\beta }(x) , can be written as follows,

(4) {\mathrm{CVaR}}_{\beta }(x)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\displaystyle {\int }_{Q(x,\omega )\ge {\mathrm{VaR}}_{\beta }(x)}Q(x,\omega )p(\omega )d\omega }.

Since model (4) is convoluted and implicit, Rockafellar and Uryasev ( 2000) developed an equivalent function to solve model (4) as follows:

(5) F_{\beta }\left(\mathit{x},\alpha \right)= \alpha +{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\displaystyle {\int }_{\omega \in {\mathbf{ℝ}}^m}\left\{\max \left\{\mathrm{0}, \left[Q\left(\mathit{x},\omega \right)-\alpha \right]\right\}\right\}p\left(\omega \right)\mathrm{d}\omega },

where F_{\beta }(x,\alpha ) is shown to be convex and continuously differentiable with respect to \alpha , and \underset{x}{\min } \{{\mathrm{CVaR}}_{\beta }(x)\} is equals to \underset{\alpha }{\min }\{F_{\beta }(x,\alpha )\} . By introducing an auxiliary variable z , the minimization of {\mathrm{CVaR}}_{\beta }[Q(x,\omega )] is equivalent to the following model ( Tong et al., 2010):

(6a) \underset{\left(x,\alpha ,z\right)}{\min }\left\{\alpha +{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\mathfrak{E}}_{\omega \in {\mathbf{ℝ}}^m}\left[z\left(\omega \right)\right]\right\},

subject to

(6b) \mathit{x}\in X,

(6c) z\left(\omega \right)\ge Q\left(\mathit{x},\omega \right)-\alpha ,

(6d) z\left(\omega \right)\ge \mathrm{0.}

Accordingly, the CVaR-based TSP method [model (2)] can be transformed into the following form:

(7a) \min f=(\mathrm{1}+\lambda )c^{T}x+\mathfrak{E}[Q(x,\omega )]+\lambda \left\{\alpha +{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\mathfrak{E}}_{\omega \in {\mathbf{ℝ}}^m}[z(\omega )]\right\},

subject to

(7b) z(\omega )\ge Q(x,\omega )-\alpha ,

(7c) z(\omega )\ge \mathrm{0},

(7d) \alpha \ge \mathrm{0},

(7e) x\in X,

with

(7f) Q(x,\omega )=\min q{(\omega )}^{T}y(\omega ),

subject to

(7g) D(\omega )y(\omega )+T(\omega )x\ge h(\omega ),

(7h) y\in Y.

By letting random variables (i.e., ω) take the discrete value ω_s with the probability level p_s (s=1,2,...,v and {\displaystyle \sum p_s}=\mathrm{1} ), model (7) can be equivalently formulated as a following linear programming model ( Noyan, 2012):

(8a) \min f=(1+\lambda )c^{T}x+{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{(q_s)}^T{y}_s}+\lambda \left(\alpha +{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{z}_s}\right),

subject to

(8b) T_{s}x+D_s{y}_s\ge h_s,s=\mathrm{1},\cdots ,v,

(8c) x\in X,

(8d) y_s\ge \mathrm{0},s=\mathrm{1},\cdots ,v,

(8e) z_s\ge {(q_s)}^T{y}_s-\alpha ,s=\mathrm{1},\cdots ,v,

(8f) \alpha \in \mathrm{ℝ},z_s\ge \mathrm{0},s=\mathrm{1},\cdots ,v.

Model (8) can deal with uncertainties in the right-hand sides of constraints presented as probability distributions when coefficients in the left-hand sides of constraints and in the objective function are deterministic. However, in real-world optimization problems, the quality of available information is mostly not satisfactory enough to be presented as probabilities ( Li et al., 2006). Such uncertainties cannot be solved through model (8). The interval-parameter programming is effective in tackling uncertainties expressed as interval values with known lower and upper bounds but unknown distribution functions ( Huang et al., 1994). Therefore, through incorporating interval-parameter programming and CVaR-based TSP method within a general optimization framework, an interval-parameter CVaR-based TSP (ICTP) method can be formulated as follows:

(9a) \min f^{\pm }=(1+\lambda ){(c^{\pm })}^T{x}^{\pm }+{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{(q_s^{\pm })}^T{y}_s^{\pm }}+\lambda \left({\alpha }^{\pm }+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{z}_s}\right),

subject to

(9b) T_s^{\pm }{x}^{\pm }+D_s^{\pm }{y}_s^{\pm }\ge h_s^{\pm },s=\mathrm{1},\cdots ,v,

(9c) z_s\ge {(q_s^{\pm })}^T{y}_s^{\pm }-{\alpha }^{\pm },s=\mathrm{1},\cdots ,v,

(9d) {\alpha }^{\pm }\in \{{\mathrm{ℝ}}^{\pm }\}, z_s\ge \mathrm{0},s=\mathrm{1},\cdots ,v,

(9e) x\in X,

(9f) y_s\ge \mathrm{0},s=\mathrm{1},\cdots ,v.

Model (9) can be transformed into two deterministic submodels corresponding to the lower and upper bounds of the objective function. This transformation process is based on an interactive algorithm, which is different from the best/worst case analysis ( Huang et al., 1994). The interval solutions associated with varying levels of constraint–violation risk can then be obtained through solving the two submodels sequentially. The submodel corresponding to the lower-bound objective function ( f^{-} ) can be first formulated as follows:

(10a) \min f^{-}=(1+\lambda )\left({\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{1}}}{c}_j^{-}{x}_j^{-}}+{\displaystyle \sum _{j=k_{\mathrm{1}}+\mathrm{1}}^{n_{\mathrm{1}}}{c}_j^{-}{x}_j^{+}}\right)+{\displaystyle \sum _{s=\mathrm{1}}^v\left({\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{p}_s{q}_{js}^{-}{y}_{js}^{-}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{p}_s{q}_{js}^{-}{y}_{js}^{+}}\right)} +\lambda \left({\eta }^{-}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\alpha }}{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{z}_s}\right),

subject to

(10b) {\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{1}}}{\left|t_{js}\right|}^{+}\mathrm{Sign}(t_{js}^{+})x_j^{-}}+{\displaystyle \sum _{j=k_{\mathrm{1}}+\mathrm{1}}^{n_{\mathrm{1}}}{\left|t_{js}\right|}^{-}\mathrm{Sign}(t_{js}^{-})x_j^{+}}+ {\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{\left|d_{js}\right|}^{+}\mathrm{Sign}(d_{js}^{+})y_{js}^{-}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{\left|d_{js}\right|}^{-}\mathrm{Sign}(d_{js}^{-})y_{js}^{+}}\ge h_s^{-},\forall s,

(10c) {\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{\left|q_{js}\right|}^{-}\mathrm{Sign}(q_{js}^{-})y_{js}^{-}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{\left|q_{js}\right|}^{-}\mathrm{Sign}(q_{js}^{-})y_{js}^{+}}\le z_s+{\eta }^{-},\forall s,

(10d) z_s\ge \mathrm{0},\forall s,

(10e) {\eta }^{-}\ge \mathrm{0},

(10f) x_j^{-}\ge \mathrm{0},j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{1}},

(10g) x_j^{+}\ge \mathrm{0},j=k_{\mathrm{1}}+\mathrm{1},k_{\mathrm{1}}+\mathrm{2},\cdots ,n_{\mathrm{1}},

(10h) y_{js}^{-}\ge \mathrm{0},\forall s,j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{2}},

(10i) y_{js}^{+}\ge \mathrm{0},\forall s,j=k_{\mathrm{2}}+\mathrm{1},k_{\mathrm{2}}+\mathrm{2},\cdots ,n_{\mathrm{2}},

where x_j^{\pm },j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{1}} , are interval variables with positive coefficients in the objective function; x_j^{\pm },j=k_{\mathrm{1}}+\mathrm{1},k_{\mathrm{1}}+\mathrm{2},\cdots ,n_{\mathrm{1}} , are interval variables with negative coefficients; y_{js}^{\pm }, j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{2}} and s=\mathrm{1},\mathrm{2},\cdots ,v , are interval random variables with positive coefficients in the objective function; y_{js}^{\pm }\ge \mathrm{0}, j=k_{\mathrm{2}}+\mathrm{1},k_{\mathrm{2}}+\mathrm{2},\cdots ,n_{\mathrm{2}} and s=\mathrm{1},\mathrm{2},\cdots ,v , are interval random variables with negative coefficients. Solutions of {\eta }_{\mathrm{opt}}^{-} , x_{j,\mathrm{opt}}^{-} ( j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{1}} ), x_{j,\mathrm{opt}}^{+} ( j=k_{\mathrm{1}}+\mathrm{1},k_{\mathrm{1}}+\mathrm{2},\cdots ,n_{\mathrm{1}} ), y_{js,\mathrm{opt}}^{-} ( j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{2}} ), and y_{js,\mathrm{opt}}^{+} ( j=k_{\mathrm{2}}+\mathrm{1},k_{\mathrm{2}}+\mathrm{2},\cdots ,n_{\mathrm{2}} ) can be obtained through submodel (10). Based on the above solutions, the second submodel for f^{+} can be formulated as follows:

(11a) \min f^{+}=(1+\lambda )\left({\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{1}}}{c}_j^{+}{x}_j^{+}}+{\displaystyle \sum _{j=k_{\mathrm{1}}+\mathrm{1}}^{n_{\mathrm{1}}}{c}_j^{+}{x}_j^{-}}\right)+{\displaystyle \sum _{s=\mathrm{1}}^v\left({\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{p}_s{q}_{js}^{+}{y}_{js}^{+}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{p}_s{q}_{js}^{+}{y}_{js}^{-}}\right)} +\lambda \left({\eta }^{+}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\alpha }}{\displaystyle \sum _{s=\mathrm{1}}^v{p}_s{z}_s}\right),

subject to

(11b) {\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{1}}}{\left|t_{js}\right|}^{-}\mathrm{Sign}(t_{js}^{-})x_j^{+}}+{\displaystyle \sum _{j=k_{\mathrm{1}}+\mathrm{1}}^{n_{\mathrm{1}}}{\left|t_{js}\right|}^{+}\mathrm{Sign}(t_{js}^{+})x_j^{-}} +{\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{\left|d_{js}\right|}^{-}\mathrm{Sign}(d_{js}^{-})y_{js}^{+}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{\left|d_{js}\right|}^{+}\mathrm{Sign}(d_{js}^{+})y_{js}^{-}}\ge h_s^{+},\forall s,

(11c) {\displaystyle \sum _{j=\mathrm{1}}^{k_{\mathrm{2}}}{\left|q_{js}\right|}^{+}\mathrm{Sign}(q_{js}^{+})y_{js}^{+}}+{\displaystyle \sum _{j=k_{\mathrm{2}}+\mathrm{1}}^{n_{\mathrm{2}}}{\left|q_{js}\right|}^{+}\mathrm{Sign}(q_{js}^{+})y_{js}^{-}}\le z_s+{\eta }^{+},\forall s,

(11d) z_s\ge \mathrm{0},\forall s,

(11e) {\eta }^{+}\ge {\eta }_{\mathrm{opt}}^{-},

(11f) x_j^{+}\ge x_{j,\mathrm{opt}}^{-},j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{1}},

(11g) \mathrm{0}\le x_j^{-}\le x_{j,\mathrm{opt}}^{+},j=k_{\mathrm{1}}+\mathrm{1},k_{\mathrm{1}}+\mathrm{2},\cdots ,n_{\mathrm{1}},

(11h) y_{js}^{+}\ge y_{js, opt}^{-},\forall s,j=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{2}},

(11i) \mathrm{0}\le y_{js}^{-}\le y_{js,\mathrm{opt}}^{+},\forall s,j=k_{\mathrm{2}}+\mathrm{1},k_{\mathrm{2}}+\mathrm{2},\cdots ,n_{\mathrm{2}},

Solutions of {\eta }_{\mathrm{opt}}^{+} , x_{j,\mathrm{opt}}^{+} ( j=\mathrm{1},\mathrm{2},\mathrm{...},k_{\mathrm{1}} ), x_{j,\mathrm{opt}}^{-} ( j=k_{\mathrm{1}}+\mathrm{1}, k_{\mathrm{1}}+\mathrm{2},\mathrm{...},n_{\mathrm{1}} ), y_{js,\mathrm{opt}}^{+} ( j=\mathrm{1},\mathrm{2},\mathrm{...},k_{\mathrm{2}} ), and y_{js,\mathrm{opt}}^{-} ( j=k_{\mathrm{2}}+\mathrm{1}, k_{\mathrm{2}}+\mathrm{2},\mathrm{...},n_{\mathrm{2}} ) can be obtained through submodel (11). Through integrating solutions of submodels (10) and (11), interval solution for model (9) under a set of p_s ( s=\mathrm{1},\mathrm{2},\mathrm{...},v ) levels can be obtained.

Case study

Air pollution issues have been of substantial concern because they are related to not only a variety of human activities and economic implications, but also to multiple human health risks ( Liu et al., 2003). Hence, effective planning for air quality management is important for facilitating sustainable socioeconomic development. In air quality management systems, uncertainties may exist in a variety of impact factors and pollution-related processes, such as pollutant characteristics, emission rates, mitigation measures, pollutant transport, 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., 2006). For example, the operating cost of SO₂ control technology may be imprecise, which can be represented as an interval parameter (with lower- and upper- bounds). Moreover, SO₂ generation rates may vary with the operating process and adopted fuel, which can be expressed as a random variable (with a probability distribution). In general, an air quality management system can be characterized by one or several emission sources that generate negative impacts on a given receptor zone ( Li et al., 2006; Lv et al., 2011). The allowable pollutant emission amounts promised for each source are critical in ensuring a healthy air quality level. Moreover, if the emission allowance is not exceeded, it will cost a standard fee for pollution control. Otherwise, if the allowance is exceeded, the surplus (excess) emission will be subject to the following economic penalties ( Lv et al., 2011): 1) operating cost for mitigating the excess pollution, 2) expenditure for buying additional allowances, and 3) a fine from the government. When penalties occur, the final pollutant emission amount will be the summation of both fixed-allowable and probabilistic surplus ones. Furthermore, the extra high SO₂ generated levels may occur with a small probability due to the variation characteristic of random variables (i.e., SO₂ generation rates), while such conditions would increase the risk of violated environmental standards. The risk cost calculated by CVaR, is used to quantify the risk of violated environmental standards. Therefore, the problem under consideration is how to minimize the system cost, including expected cost and risk cost, for pollution abatement while satisfying the overall environmental goal. Specifically, formulation of desired plans for air pollution is critical for: 1) satisfying the environmental requirements, 2) reducing the risk of violated environmental standards, 3) minimizing the expected cost (i.e., the sum of regular and penalty costs) and risk cost, and 4) identifying effective management policies.

In this study, two power plants (i = 1, 2), a chemical industry plant (i = 3), a petroleum refinery (i = 4), and a steel mill (i = 5) are considered as SO₂ emission sources. Five receptor zones including one residential zone (p = 1), one scenic zone (p = 2), two agricultural zones (p = 3 and p = 5), and one industrial zone (p = 4) are affected by the SO₂ emission from five emission sources. Figure 1 shows their locations in the hypothetical system. Investigations on the local observed meteorological data show that the stability of class D and east wind have the highest occurrence frequency. The effective stack height of five emission sources are estimated at 252 m, 265 m, 150 m, 110 m, and 100 m, respectively ( Liu et al., 2003). The receptor zones are adversely affected by the SO₂ generated from the five emission sources. Because the sulfur content of the coals used varied and the related combustion conditions are different, the relevant SO₂ emission rates can be expressed as random variables ( Li et al., 2006). Meanwhile, the stochastic variability of random SO₂ generation rates may result in the extreme expected loss. The data of SO₂ generation rates within the planning period are acquired through survey and forecast measures. 1,000 specimens of SO₂ generation rates are selected for each emission source. Thus, 5,000 specimens would be needed in the study. Kolmogorov–Smirnov tests are then used to identify the probability distribution function (PDF) of SO₂ generation rates. Figure 2 shows the frequency of the SO₂ generation rate for each emission source, which indicates they follow the normal distribution. Next, the normal probability distribution functions can be assumed to be discrete with 20 scenarios (i.e., 20 probability levels from extra low to extra high). Table 1 presents the discrete results, which show the interval SO₂ generation rates and the associated probabilities of occurrence. To keep the environmental quality at a more healthy level, the allowable emission amounts for five sources will be limited to a relatively low level which must be less than [40.0, 44.4], [26.0, 29.3], [3.7, 4.1], [7.5, 9.2] and [0.9, 1.5] tonnes/day, respectively. Moreover, the SO₂ loading capacities in the five receptor zones must be limited to under 0.15 mg/m³, 0.05 mg/m³, 0.15 mg/m³, 0.25 mg/m³, and 0.15 mg/m³, respectively.

For the study system, even the use of low-sulfur coal (e.g., 0.2%) can hardly satisfy the SO₂ emission standard. Nevertheless, coals being mined naturally contain [0.2, 5]% of sulfur in China. Thus, to meet the environmental requirement, each emission source has to adopt pollution control technologies, including oxidation packed absorption (OPA), alkali absorption (AA), and limestone wet scrubbing (LWS), to mitigate the SO₂ emission and to avoid penalties from the government. Table 2 presents the parameters for different control technologies (i.e., efficiency, regular operation cost, and penalty operation cost)( Liu et al., 2003; Li et al., 2006; Lv et al., 2011). The efficiency of each pollution-control technology may vary with different operating conditions such as reagent ratios, temperature, and inlet SO₂ concentrations. Moreover, an elevated SO₂ level in the ambient air may cause a wide variety of human health threats and environmental impacts. According to the standard of the World Health Organization, the SO₂ loading level does not exceed the threshold of 50 μg/(m³·year) in this study. The regular operation costs involve the capital and operating costs for mitigating SO₂ emissions through three pollution-control techniques. The penalty rates are significantly higher than the regular costs. The AA technology has the highest SO₂ removal efficiency, but it also has the highest operating cost due to the highest reagent cost.

This study focuses on the mitigation of SO₂ emission from multiple sources through different control measures. The pre-set policy formulated by local authorities is considered when a decision analysis is undertaken. Uncertainties that exist in the system components may be presented as PDFs and discrete intervals. Moreover, the risks from the variability of random parameters are quantified as CVaR. Therefore, the proposed ICTP method is suitable for tackling this type of management problem. Thus, we have:

\min f^{\pm }=\left(\mathrm{1}+\lambda \right){\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{5}}{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{3}}{X}_{ij}^{\pm }A{O}_j^{\pm }}}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{5}}{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{3}}{\displaystyle \sum _{s=\mathrm{1}}^{\mathrm{20}}{p}_{is}{Y}_{ijs}^{\pm }E{O}_j^{\pm }}}}+\lambda \left({\alpha }^{\pm }+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\beta }}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{5}}{\displaystyle \sum _{s=\mathrm{1}}^{\mathrm{20}}{p}_{is}{z}_{is}}}\right)

subject to

(12a) z_{is}\ge {\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{3}}E{O}_j^{\pm }{Y}_{ijs}^{\pm }}-{\alpha }^{\pm },\forall i,s

[Constraint of auxiliary variable]

(12b) {\displaystyle \sum _{i=1}^5\left(X_{ij}^{\pm }+Y_{ijs}^{\pm }\right)}\le T{C}_j^{\pm },\forall j,s

[Constraint of capacity for each technology]

(12c) {\displaystyle \sum _{j=1}^3\left(X_{ij}^{\pm }+Y_{ijs}^{\pm }\right)}\ge G_{is}^{\pm },\forall i,s

[Constraint of SO₂ mitigation requirement]

(12d) {\displaystyle \sum _{j=1}^3\left(1-{\eta }_j^{\pm }\right)\left(X_{ij}^{\pm }+Y_{ijs}^{\pm }\right)}\le e_i^{\pm },\forall i,s

[Constraint of SO₂ emission allowance]

(12e) {\displaystyle \sum _{i=1}^5{\displaystyle \sum _{j=1}^3{t}_{ip}u\left(1-{\eta }_j^{\pm }\right)\left(X_{ij}^{\pm }+Y_{ijs}^{\pm }\right)}}\le a_p^{\pm },\forall p,s

[Constraint of ambient air quality requirement]

(12f) {\alpha }^{\pm }\ge 0

(12g) z_{is}\ge \mathrm{0},\forall i,s

(12h) \mathrm{0}\le Y_{ijs}^{\pm }\le X_{ij}^{\pm }\le X_{ij,\max },\forall i,j,s

[Non-negativity constraints]

where f^{\pm } is the mean-risk function value ($); i is the name of SO₂ emission source, i = 1, 2,…, 5; j is the type of SO₂ control technology, j = 1, 2, 3, where j = 1 for the OPA technology, j = 2 for the AA technology, and j = 3 for the LWS technology; s is the scenario of SO₂ emission rates, s = 1, 2,…, 20; p is the name of receptor zone, p = 1, 2,…, 5; λ and β are risk parameters; z_{is} is the auxiliary variable; A{O}_j^{\pm } is the operating cost of control technology j for allowable SO₂ emissions ($/tonne) (the first-stage cost parameter); E{O}_j^{\pm } is the operating and penalty costs of control technology j for excess SO₂ emissions ($/tonne) (the second-stage cost parameter), where E{O}_j^{\pm }>A{O}_j^{\pm } ; T{C}_j^{\pm } is the maximum allowable capacity of control technology j (tonne/day); e_i^{\pm } is the SO₂ emission allowance for emission source i (tonne/day); G_{is}^{\pm } is the amount of SO₂ generated from source i under scenario s of SO₂ generation rates (tonne/day); {\eta }_j^{\pm } is the efficiency of control technology j; u is the unit conversion factor changing “tonne/day” into “mg/s,” where the work time of emission sources varies between 8 and 10 h and the value of u equals [2.78, 3.47] × 10⁴; X_{ij}^{\pm } is the amount of SO₂ generated from source i to be mitigated by control technology j, which can be derived based on the emission standard as pre-regulated by the authorities (the first-stage decision variable) (tonne/day); Y_{ijs}^{\pm } is the amount by which the SO₂ emission allowance ( X_{ij}^{\pm } ) is exceeded when the SO₂ generation rate is G_{is}^{\pm } with the probability p_{is} (tonne/day) (the second-stage decision variable); a_p^{\pm } is the environmental loading capacity of receptor zone p (mg/m³); t_{ip} is the contribution of unit pollutant emission rate at emission source i to receptor zone p. t_{ip} represents the mapping relationship between SO₂ emission source i and receptor zone p, which can be calculated based on the Gaussian dispersion model.

In the Gaussian dispersion model, the ground-level concentration of SO₂ at each arbitrary downwind location (x, y) can be estimated through the following equation:

(13) C\left(x,y\right)={\displaystyle \frac{Q}{\mathrm{\pi }\overline{u}{\sigma }_y{\sigma }_z}}\exp \left[\left({\displaystyle \frac{-y^{\mathrm{2}}}{\mathrm{2}{\sigma }_y^{\mathrm{2}}}}\right)-\left({\displaystyle \frac{H^{\mathrm{2}}}{\mathrm{2}{\sigma }_z^{\mathrm{2}}}}\right)\right],

where H represents average effective stack height (m); Q denotes pollutant emission rate (mg/s); \overline{u} is wind velocity (m/s); {\sigma }_y and {\sigma }_z represent standard deviations of the plume at y and z directions (m), respectively. Through fitting the Pasquill-Gifford curves, {\sigma }_y and {\sigma }_z can be calculated as follows ( Turner, 1970):

(14a) {\sigma }_y=a{x}^b,

(14b) {\sigma }_z=c{x}^d,

where the coefficients a, b, c and d depend on Pasquill Stability Category. Let x_{ip} denote the distance between the i^th SO₂ emission source and p^th receptor zone at x directions (km). Meanwhile, let y_{ip} denote the distance between the i^th SO₂ emission source and p^th receptor zone at y directions (km). Therefore, the concentration contribution of the i^th emission source to the ground concentration in the p^th receptor zone can be estimated as follows:

(15) C\left(x_{ip},y_{ip}\right)={\displaystyle \frac{Q_i}{\mathrm{\pi }\overline{u}ac{x}_{ip}^{b+d}}}\exp \left[\left({\displaystyle \frac{-y_{ip}^{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}{x}_{ip}^{\mathrm{2}b}}}\right)-\left({\displaystyle \frac{H^{\mathrm{2}}}{\mathrm{2}{c}^{\mathrm{2}}{x}_{ip}^{\mathrm{2}d}}}\right)\right],

where Q_i denotes the pollutant emission rate of the i^th SO₂ emission source (mg/s). Consequently, the transfer factor t_{ip} can be defined as follows:

(16) t_{ip}={\displaystyle \frac{x_{ip}^{-b-d}}{\mathrm{\pi }\overline{u}ac}}\exp \left[\left({\displaystyle \frac{-y_{ip}^{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}{x}_{ip}^{\mathrm{2}b}}}\right)-\left({\displaystyle \frac{H^{\mathrm{2}}}{\mathrm{2}{c}^{\mathrm{2}}{x}_{ip}^{\mathrm{2}d}}}\right)\right].

In model (12), the objective is to minimize the sum of system expected cost for SO₂ abatement and the CVaR quantifying the risk of violated environmental standard; the constraints include the SO₂ emission standards (or regulated allowances) and the allowable pollutant-loading levels. Figure 3 illustrates the general framework of the ICTP method model (12)], which is based on the interval-parameter programming, CVaR and TSP technologies. The solution algorithm of model (12) can then be summarized through using the following pseudo-code:

Step 1: Sample the SO₂ concentration for each emission source and obtain the related PDF;

Step 2: Employ the Gaussian dispersion model to calculate the transfer factor ( t_{ip} );

Step 3: Formulate the ICTP model [model (12)] for risk management of sulfur dioxide abatement;

Step 4: Introduce the discretized PDF and the transfer factor ( t_{ip} ) into model (12);

Step 5: Transform model (12) into two submodels, which correspond to f^{-} and f^{+} respectively, based on the methodology developed in Section 2;

Step 6: Solve the two submodels, respectively, through the Lingo 9.0 optimization software;

Step 7: Obtain the solution of the ICTP model [model (12)] and get the optimal SO₂ emission allocation pattern under different risk levels;

Step 8: Stop.

Result analysis

In this study, two risk parameters, five levels of λ (i.e., 0.5, 1, 2, 5, 10) and three levels of β (i.e., 0.6, 0.8 and 0.99), were considered. Table 3 presents the solutions obtained through the ICTP method under λ = 1 and β = 0.99, for abating the allowable and excess SO₂ emissions (i.e., the first- and second-stage variables). The solutions presented as interval numbers indicate that the related decisions should be sensitive to the change of input parameters. An excess SO₂ emission will be generated if the SO₂ emission allowance is exceeded (i.e., excess SO₂ emission= SO₂ generation rate − SO₂ emission allowance). If the SO₂ generation rate is higher, a control technology with a higher efficiency should be installed to abate the excess emissions for meeting the air quality standard and reducing the economic penalties. Therefore, the amount treated by the AA technology should be assigned first in case of excess SO₂ emissions. For example, in power plant 1, the amounts of allowable SO₂ emission to be treated by OPA, AA, and LWS technology would be 30 tonnes/day, 30 tonnes/day, and 24 tonnes/day, respectively. In addition, some excess SO₂ emissions should be treated under varied SO₂ generation rates (with different probability levels). The amounts of excess emission mitigated by OPA, AA, and LWS technology would respectively be 10.3 tonnes/day, 18.7 tonnes/day, and [19.8, 23.4] tonnes/day under the highest SO₂ generation rate (in scenario 20 with the probability of 0.1%), while no excess emission was allotted under the lowest SO₂ generation rate (in scenario 1 with the probability of 0.2%).

Table 4 illustrates the optimized allocation pattern for the total SO₂ emission from the five sources to the OPA, AA, and LWS technologies under different risk levels. The results indicate that the optimized scheme for SO₂ emission allocation would be related to the risk levels (i.e., λ and β). For example, for β = 0.6, the amounts of total SO₂ emission abated by the OPA technology should be [34.2, 34.7] tonnes/day under λ = 0.5, [32.3, 32.5] tonnes/day under λ = 1, [31.5, 31.7] tonnes/day under λ = 2, [31, 31.7] tonnes/day under λ = 5, [30.9, 31.7] tonnes/day under λ = 10; while for λ = 1, the amounts of total SO₂ emission mitigated by the OPA technology should be [32.3, 32.5] tonnes/day under β = 0.6, [31.8, 32.1] tonnes/day under β = 0.8, [31.4, 31.7] tonnes/day under β = 0.99. Figure 4 shows the changing trend of the excess SO₂ emissions treated by the OPA, AA, and LWS technologies under different λ and β levels. It is indicated that the excess SO₂ emissions are affected by the risk level as well. For instance, the excess SO₂ emissions to any control technologies would decrease when λ is fixed and β changes from 0.6 to 0.99 (or β is fixed and λ changes from 0.5 to 10). Therefore, lower λ (and/or β) means a higher level of system risk, which may lead to more excess SO₂ emissions; conversely, higher λ (and/or β) means a lower level of system risk, which may lead to less excess SO₂ emissions.

Figure 5(a) shows the mean-risk function value of f^{\pm } at different risk levels. The result indicated that the optimal value of the mean-risk function would increase when the risk parameters (λ and/or β) increase. For example, the mean-risk function value would increase from $[14.9, 17.6] × 10³ to $[123.8, 145.8] × 10³ when β is fixed at 0.6 and λ increases from 0.5 to 10. In addition, the mean-risk function value would increase from $[20.7, 24.3] × 10³ to $[23.2, 27.1] × 10³ when λ is fixed at 1 and β increases from 0.6 to 0.99. Figure 5 (b) illustrates the total expected cost under different λ and β levels. It is indicated that the total expected cost would increase due to higher risk parameters (λ and/or β). This implies that a higher system cost would guarantee a lower system risk; conversely, if the decision maker aims toward an economical plan, a higher expected loss could be confronted. In the real world application, the decision makers may need to choose between a riskier solution with a lower system cost and a more conservative solution with a higher system cost.

Figure 6 presents the recourse cost under different λ and β levels. It is indicated that risk parameters can affect the recourse cost. When fixing λ at 0.5 and increasing β from 0.6 to 0.99, the recourse cost would decrease from $[1.82, 2.51] × 10³ to $[1.57, 2.20] × 10³. Besides, the recourse cost would decrease from $[1.57, 2.20] × 10³ to $[1.00, 1.50] × 10³ when β is fixed at 0.8 and λ changes from 0.5 to 10. According to these results, increasing λ may lead to a more risk-averse policy with lower recourse costs in general. Therefore, a more risk-averse policy keeps lower excess SO₂ emissions and results in a lower recourse cost. Figure 7 presents the values of CVaR under different λ and β levels. It appears that CVaR increases as β increases due to the definition of CVaR. For example, when λ is fixed at 0.5 and β changes from 0.6 to 0.99, the value of CVaR would increase from $[11.7, 13.7] × 10³ to $[14.2, 16.4] × 10³. CVaR_β quantifies the expectation of the worst (1‒β) % of the total costs, and a specified β level indicates the confidence level in percentage terms. At a larger value of β, the corresponding VaR may increase, and CVaR_β would account for the risk of larger realizations. Therefore, larger β values would lead to more conservative policies, which give more weight to worse scenarios. In addition, CVaR decreases as λ increases. For example, when β is fixed at 0.8 and λ changes from 0.5 to 10, the value of CVaR would decrease from $[12.2, 14.2] × 10³ to $[11.9, 14.0] × 10³. This is because of the changing trade-off between the expectation total cost and the CVaR criterion. A higher λ would increase the relative importance of the risk term, leading to more risk-averse policies. Therefore, increasing the parameter λ (and/or β) indicates a higher level of risk aversion.

Discussion

To better reflect the advantages of the proposed ICTP method, an interval-parameter TSP (ITSP) model is applied to the study case for comparison purposes. The total expected cost solved by ITSP model is $[9.1, 10.7] × 10³, which is lower than that solved by ICTP method under any risk parameters. However the recourse cost, which equals $[1.9, 2.6] × 10³, solved by ITSP model is higher than those solved by ICTP method under any risk parameters. Table 5 presents the solutions of ITSP for decision variables. Different from the ICTP method, the optimal solution in the ITSP model is to obtain the minimum cost regardless of risk aversion. Consequently, it is incapable of analyzing the trade-off between the system cost and risk of violated environmental standards. For example, the lower levels of allowable SO₂ emission (i.e., X_{ik,opt}^{\pm } ), higher levels of excess SO₂ emission (i.e., Y_{isk,opt}^{\pm } ), and the lower system cost would be generated by the ITSP than the ICTP. This implies that the system cost largely relies on the SO₂ generation condition, and tends to fluctuate more intensively with the change of SO₂ generation rates when the CVaR is preferred with more weight. Furthermore, the solutions in the ITSP seem to be exactly the same as those in ICTP when λ is zero, which indicates that the ITSP can be considered a special or simplified case of the ICTP. Therefore, the risk-averse method (i.e., ICTP) that considers the effects of random variables would provide more robust solutions compared to the risk-neutral approach (i.e., ITSP).

Without consideration of interval-parameter programming, the CVaR-based TSP method can be used to solve the study case through letting all interval parameters be equal to their mid-value. The solution is a set of deterministic values. It represents a decision under an input scenario (mid-values for all parameters), and is actually one of many alternatives from the ICTP. In such a case, the decision alternative would be restricted to a single solution. Thus, the effectiveness and flexibility of the alternative would be reduced since impact of uncertainties is not considered. Although further sensitivity analysis can be undertaken for the CVaR-based TSP solution, a number of possibilities exist when many inputs are uncertain. It can increase the computational burdens. For each possibility, the sensitivity analysis can provide an individual response to variations of the uncertain inputs; however, it can hardly reflect interactions among these uncertainties. In ICTP, the interactive algorithm is used to handle the interval parameters. As a result, the solution corresponding to f^{-} can be first solved (when the objective is to be minimized), and the relevant solution corresponding to f^{+} is to be feasible as one of the two bounds of the interval solution, which leads to a set of optimal and stable interval solutions. In comparison, the best/worst case analysis ( Liu et al., 2011), which can also address the interval parameters, would obtain the solution results under two extreme scenarios (i.e., best and worst conditions). It is useful for judging the system’s capability to realize the desired goal; but will not necessarily construct a set of stable intervals for decision variables.

The robust risk analysis method (RRAM) was introduced in Chen et al. (2013a, b), which incorporated interval-parameter programming, two-stage stochastic programming, and variance within a general modeling framework. Thus, the main difference between ICTP and RRAM is the quantitative risk measure for the second-stage random variable, where the former uses the CVaR and the latter uses the variance. Moreover, compared with ICTP, the disadvantage of RRAM is that it considers the under-and-over-performances of the probability distribution equally. For example, the high SO₂ generation rate may result in the system risk, while the low waste generation rate is safe to the management system. Therefore, quantitatively evaluating the right tail of distributions of the SO₂ generation rate is necessary. However, consideration of the left tail of distributions must not only waste the expected cost, but also be of no use for the system stability. Moreover, Shao et al. (2011) developed a risk-averse method through incorporating the CVaR into the constraint of the optimization model. This method also can analyze the trade-offs between system cost and risk of violated environmental standard. However, it may have computational difficulties in estimating the maximum acceptable loss (i.e., the limit value of right hand in the constraint containing CVaR).

Conclusions

An interval-parameter CVaR-based two-stage programming (ICTP) method was developed by incorporating the interval-parameter programming, two-stage stochastic programming (TSP) and conditional value at risk (CVaR) within a general optimization framework. In the developed ICTP, uncertainties could be presented as probabilistic distribution functions and discrete intervals. Also, CVaR as a kind of risk measurement method was successfully incorporated within the objective function of ICTP to describe the expected losses under extreme conditions. In general, the proposed ICTP method had several advantages: (i) its objective function simultaneously took expected cost and risk cost into consideration, and also used discrete random variables and discrete intervals to reflect uncertain properties. Therefore, the solution results could possess characteristics of CVaR, interval-parameter programming, and the TSP model; (ii) it quantitatively evaluated the right tail of distributions of random variables, which could better calculate the risk of violated environmental standards compared to the robust risk analysis method; (iii) it was useful for helping decision makers to analyze the trade-offs between cost and risk, and identified desired air quality management strategies under complex uncertainties; and (iv) it was effective to penalize the second-stage costs, as well as to capture the notion of risk in stochastic programming.

An air quality management case was used to demonstrate the applicability of the ICTP method. The obtained decision variables values can be used to reflect decision makers’ opinions on the SO₂ abatement planning; and the risk acceptance λ and the percentile β can reflect decision maker’s preference toward the system cost and risk control. As multiple scenarios of solutions are normally available, it is suggested that the decision makers should use multi-criteria decision analysis (MCDA) or data envelopment analysis (DEA) technologies for evaluating and ranking various possible alternatives. If the decision makers need only one concrete planning pattern, they should clearly identify their preferences on system economy and risk, and specify the corresponding SO₂ abatement alternative for meeting their requirement. Thus, the ICTP method is valuable for supporting the adjustment or justification of the existing SO₂ emission allocation patterns. On the other hand, to reflect the pollutants dispersion in the atmospheric environment, the Pasquill-Gifford curves were used to calculate the related parameters for the Gaussian dispersion model. Although somewhat simplified, the curves are representative and straightforward to understand. In fact, the variation of the atmosphere stability degree and the dispersion parameters is continuous due to the continuous changes of the atmosphere turbulence. However, the turbulence-related parameters in the Gaussian dispersion model are treated as discrete values, which may lead to the deviations of solutions in the air quality management model. The California puff (CALPUFF) modeling system can predict the ground-level concentration of SO₂ with a higher accuracy, through the comprehensive consideration of boundary layer parameters such as Monin-Obukhov length, convection velocity, boundary height friction velocity, and surface roughness length. However, it requires intensive data inputs. Moreover, incorporating the outputs of CALPUFF into the optimization models is a complex problem and deserves further investigations.

References

[1] Ahmed S (2004). Mean-risk objectives in stochastic programming. Technical report, Georgia Institute of Technology, available at 

[2] An H, Eheart J W (2007). A screening technique for joint chance-constrained programming for air-quality management. Oper Res, 55(4): 792–798

[3] Athanasiadis I N, Mitkas P A (2007). Knowledge discovery for operational decision support in air quality management. Journal of Environmental Informatics, 9(2): 100–107

[4] Birbil S I, Frenk J, Kaynar B, Noyan N (2009). Risk measures and their applications in asset management. In: Gregoriou G N, ed. The VaR Implementation Handbook. New York: The McGraw-Hill Companies, 311–337

[5] Birge J R, Louveaux F V (1988). A multicut algorithm for two-stage stochastic linear programs. Eur J Oper Res, 34(3): 384–392

[6] Byun D W, Kim S T, Chen F Y, Kim S B, Cuclis A, Moon N K (2003). Information infrastructure for air quality modeling and analysis: application to the houston-galveston ozone non-attainment area. Journal of Environmental Informatics, 2(2): 38–57

[7] Cai Y P, Huang G H, Tan Q, Yang Z F (2011a). An integrated approach for climate-change impact analysis and adaptation planning under multi-level uncertainties. Part I: methodology. Renewable & Sustainable Energy Reviews, 15: 2779–2790

[8] Cai Y P, Huang G H, Tan Q, Liu L (2011b). An integrated approach for climate-change impact analysis and adaptation planning under multi-level uncertainties. Part II. case study. Renewable & Sustainable Energy Reviews, 15: 3051–3073

[9] Chen C, Huang G H, Li Y, Li M (2013a). Model of risk analysis on site selection of biomass power plant based on stochastic robust interval method. Transactions of the Chinese Society of Agricultural Engineering, 29(20): 206–213

[10] Chen C, Huang G H, Li Y P, Zhou Y (2013b). A robust risk analysis method for water resources allocation under uncertainty. Stochastic Environ Res Risk Assess, 27(3): 713–723

[11] Dai C, Li Y, Huang G (2012). An interval-parameter chance-constrained dynamic programming approach for capacity planning under uncertainty. Resour Conserv Recycling, 62: 37–50

[12] Dentcheva D, Ruszczynski A (2006). Portfolio optimization with stochastic dominance constraints. J Bank Finance, 30(2): 433–451

[13] Dong C, Huang G H, Cai Y P, Liu Y (2012). An inexact optimization modeling approach for supporting energy systems planning and air pollution mitigation in Beijing City. Energy, 37(1): 673–688

[14] Fábián C I (2008). Handling CVaR objectives and constraints in two-stage stochastic models. Eur J Oper Res, 191(3): 888–911

[15] Haurie A, Kübler J, Clappier A, Van Den Bergh H (2004). A metamodeling approach for integrated assessment of air quality policies. Environ Model Assess, 9(1): 1–12

[16] Hsu C P, Huang C W, Paul Chiou W J P (2012). Effectiveness of copula-extreme value theory in estimating value-at-risk: empirical evidence from Asian emerging markets. Rev Quant Finance Account, 39(4): 447–468

[17] Huang G, Loucks D (2000). An inexact two-stage stochastic programming model for water resources management under uncertainty. Civ Eng Environ Syst, 17(2): 95–118

[18] Huang G H, Baetz B W, Patry G G (1994). Grey dynamic programming for waste-management planning under uncertainty. J Urban Plann Dev, 120(3): 132–156

[19] Huang G H, Chang N B (2003). Perspectives of environmental informatics and systems analysis. J Environ Inform, 1(1): 1–6

[20] Inuiguchi M, Sakawa M (1998). Robust optimization under softness in a fuzzy linear programming problem. Int J Approx Reason, 18(1‒2): 21–34

[21] Lejano R P, Ayala P M, Gonzales E A (1997). RESEARCH: optimizing the mix of strategies for control of vehicular emissions. Environ Manage, 21(1): 79–87

[22] Li Y, Huang G, Nie X, Nie S (2008). An inexact fuzzy-robust two-stage programming model for managing sulfur dioxide abatement under uncertainty. Environ Model Assess, 13(1): 77–91

[23] Li Y, Huang G H, Veawab A, Nie X, Liu L (2006). Two-stage fuzzy-stochastic robust programming: a hybrid model for regional air quality management. J Air Waste Manag Assoc, 56(8): 1070–1082

[24] Lim C, Sherali H D, Uryasev S (2010). Portfolio optimization by minimizing conditional value-at-risk via nondifferentiable optimization. Comput Optim Appl, 46(3): 391–415

[25] Liu L, Huang G, Liu Y, Fuller G, Zeng G (2003). A fuzzy-stochastic robust programming model for regional air quality management under uncertainty. Eng Optim, 35(2): 177–199

[26] Liu Y, Zou R, Guo H (2011). Risk explicit interval linear programming model for uncertainty-based nutrient-reduction optimization for the Lake Qionghai Watershed. J Water Resour Plann Manage, 137(1), 83–91

[27] Lv Y, Huang G H, Li Y P, Yang Z F, Sun W (2011). A two-stage inexact joint-probabilistic programming method for air quality management under uncertainty. J Environ Manage, 92(3): 813–826

[28] Ma X, Zhang F (2002). A genetic algorithm based stochastic programming model for air quality management. J Environ Sci (China), 14(3): 367–374

[29] Marti K, Ploöchinger E (1990). Optimal step sizes in semi-stochastic approximation procedures. I. Optimization, 21(1): 123–153

[30] Miller N, Ruszczynski A (2008). Risk-adjusted probability measures in portfolio optimization with coherent measures of risk. Eur J Oper Res, 191(1): 193–206

[31] Noyan N (2012). Risk-averse two-stage stochastic programming with an application to disaster management. Comput Oper Res, 39(3): 541–559

[32] Rockafellar R T, Uryasev S (2000). Optimization of conditional value-at-risk. Journal of Risk, 2: 21–42

[33] Ruszczyński A (1993). Parallel decomposition of multistage stochastic programming problems. Math Program, 58(1–3): 201–228

[34] Schultz R (2011). Risk aversion in two-stage stochastic integer programming. Stochastic Programming: 165–187

[35] Schultz R, Tiedemann S (2006). Conditional value-at-risk in stochastic programs with mixed-integer recourse. Math Program, 105(2‒3): 365–386

[36] Seifi A, Hipel K (2001). Interior-point method for reservoir operation with stochastic inflows. J Water Resour Plan Manage, 127(1): 48–57

[37] Shao L G, Qin X S, Xu Y (2011). A conditional value-at-risk based inexact water allocation model. Water Resour Manage, 25(9): 2125–2145

[38] Tan Q, Huang G H, Wu C Z, Cai Y P (2011). IF-EM: an interval-parameter fuzzy linear programming model for environment-oriented evacuation planning under uncertainty. J Adv Transport, 45(4): 286–303

[39] Tong X J, Qi L Q, Wu F, Zhou H (2010). A smoothing method for solving portfolio optimization with CVaR and applications in allocation of generation asset. Appl Math Comput, 216(6): 1723–1740

[40] Turner D B (1970). Workbook of atmospheric dispersion estimates. Washington D C: Office of Air Programs, Publication No. AP-26, U.S. Environmental Protection Agency

[41] Wang L, Milford J B (2001). Reliability of optimal control strategies for photochemical air pollution. Environ Sci Technol, 35(6): 1173–1180

[42] Zolezzi T (2001). Well-posedness and optimization under perturbations. Ann Oper Res, 101(1/4): 351–361

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg