{'fraction_non_alphanumeric': 0.05591866376180102, 'fraction_numerical': 0.016864358912289195, 'mean_word_length': 4.2172174279099135, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pollution control for switching diffusion models: Approximation methods and numerical results

This work focuses on optimal pollution controls. The main effort is devoted to obtaining approximation methods for optimal pollution control. To take into consideration of random environment and other random factors, the control system is formulated as a controlled switching diffusion. Markov chain approximation techniques are used to design the computational schemes. Convergence of the algorithms are obtained. To demonstrate, numerical experimental results are presented. A particular feature is that computation using real data sets is provided.

1. Introduction. Professor Peter E. Kloeden has made fundamental contributions to stochastic differential equations, numerical methods, and applications. The book by Kloeden and Platen [8] has been used training a generation of mathematicians and scientists. An updated treatments can be found in Han and Kloeden [3] and references therein. Devoted to numerical solutions of pollution management problems, on the happy occasion of celebrating his 70th birthday, we dedicate this paper to Peter.
Environment pollution was identified as responsible for 9 million premature deaths in 2015, or 16% of all deaths worldwide. The Lancet Commission on Pollution and Health estimated that in 2016, ambient air pollution alone costs the global economy 5.7 trillion dollars, or 4.4 percent of global GDP. Consequently, it has been widely recognized that pollution management is of vital importance for environmental health, social welfare, and world economy. A major issue of pollution management is concerned with the tradeoff of pollution accumulation and consumption.
Following the seminal paper of Keeler et al. [6], much work has been devoted to the study of optimal control of the related dynamic systems or economic systems. In [5], Kawaguchi and Morimoto treated a pollution accumulation problem of maximizing the long-run average welfare using a controlled diffusion model. Assume that an economy consumes some goods and also generates certain amount of pollution. The pollution stock is gradually degraded and its instantaneous growth rate encounters randomness that can be approximately represented by a random disturbance with mean zero and constant variance. The social welfare is then defined as a suitable combination of the utility of the consumption and the disutility of pollution. The basic optimal control problem is to find the optimal consumption strategies for the society in the long-run average sense.
Departing from the existing frameworks, we consider a general and yet more realistic model. Real-world economic systems are highly complex. It is difficult to capture the underlying structural and dynamic changes by merely the traditional diffusion models. To accommodate high level uncertainty, we introduce a discreteevent component. The new model is of hybrid in nature in the sense that it involves both continuous dynamics and discrete events. The continuous and discrete events coexist and interact. Some of the recent formulations and related work can be found in Yin and Zhu [16]; see also Jasso-Fuentes and Yin [4]. In this paper, we use a continuous-time Markov chain to represent the discrete event process. In our model, the Markov chain serves as a modulating force. The systems that we consider in this paper are controlled switching diffusions.
Within this formulation, stochastic control techniques can be employed to find optimal controls and establish their essential properties. At present, closed-form solutions to controlled switching diffusions are difficult to find. In order to implement the control policy in practice, one has to resort to numerical solutions. This paper is devoted to computational methods and numerical solutions for the class of controlled switching diffusion arising from pollution control problems. However, computable solutions are highly appreciated.
Our contributions of this paper can be summarized as follows.
(i) We propose a new model to address issues arising in environment management. Switching processes are incorporated into the existing diffusion models to capture random structural and dynamic changes in pollution dynamics that were not considered in the existing literature. (ii) We develop new feasible numerical methods for the resulting stochastic control problem. By employing the key idea of Markov chain approximation techniques, we construct discrete-time controlled Markov chains, and show that the controlled chains are consistent with the original continuous-time controlled switching diffusions. Furthermore, we prove that our numerical schemes are convergent by means of the weak convergence method. (iii) We verify our theoretical results by a number of numerical case studies. (iv) We further validate our results by using real-world data. One of the interesting aspects is that this paper might be one of the first to use numerical stochastic control techniques to deal with real data problem arising from pollution con-trols. Naturally, in analyzing real data, there are unknown parameters and also the switching states may be hidden. We apply statistical methods from machine learning to estimate the parameters and hidden state. In addition, the principle component analysis is used to resolve the curse of dimensionality issues when various sources of pollution are involved. Such an attempt paves a way for a wide range of applications not only in pollution control but also other environment management problems in the real-world scenarios.
The rest of the paper is arranged as follows. The problem formulation and controlled systems are presented in Section 2. Section 3 discusses main properties of the associated value functions. Numerical algorithms are constructed to treat the problem in Section 4. The convergence of the numerical algorithm is demonstrated in Section 5. To illustrate how the numerical algorithms can be used, we provide a number of case studies in Section 6. This section also provides computational results using real data. Finally, Section 7 concludes the paper with some further remarks.
2. Problem formulation. Suppose that d, d 1 , d 2 are positive integers, and that M = {1, 2, . . . , m} is a finite set. Let C(t) = (c 1 (t), . . . , c d1 (t)) denote the consumption rate of goods, G ∈ R d×d1 denote the generation of polluting materials when consumption plan implemented, and X(t) = (x 1 (t), . . . , x d (t)) be the amount of which represents the self growth rate and interaction between the pollution materials, where ξ j (i) are row vectors for j = 1, · · · , d, σ(x, i) is an R d×d2 matrix, α(t) is a continuous-time Markov chain with state space M, and W (·) is an R d2 -dimensional vector-valued standard Brownian motion independent of α(t).
The amount of pollution materials X(t) satisfies the system equation The corresponding objective function is given by where U (·) is a utility function satisfying the following conditions: (H1) Given x and i, suppose that U is a smooth function w.r.t. the control C in that U (x, i, ·) ∈ C 2 [0, ∞).
We say that the consumption rate is admissible if it is F t -progressively measurable, where F t = {(X(s), α(s)) : s ≤ t} such that 0 ≤ |C(t)| ≤ K 0 for some K 0 > 0. Assume that the pollution level α(·) is a continuous-time Markov chain with a finite state space M = {1, . . . , m}. As another generalization of Morimoto [5], we assume that σ to be dependent of (X(t), α(t)).
Example 2.1. A particular, one-dimensional example is the following pollution control model. Let c(t) ∈ R, t ≥ 0, denote the consumption rate. Here c(t) represents how fast the underlying economy consumes some goods and in the meantime generates pollution. The rate c(t) is real valued also called the flow of pollution and is the control process. Let X(t) be a real-valued process denoting the stock of pollution at time t and let α(t) be a finite-state continuous-time Markov chain representing the environment mode. The pollution decay rate is given by ρ(α(t)), a function of α(t). These variables satisfy the following stochastic differentiable equation where σ 0 is a function of (x, α) and W (t) is a standard Brownian motion. In this paper, the environment mode process α(t) represents the states of macro economic factors such as the level of resources available for environment protection, government regulatory policies on pollution control, etc. We assume that α(t) takes values in M = {1, 2, . . . m}.
The overall objective is the consumption of goods at a level for maintaining good economic activities and in the meantime keeping pollution in check. To quantify, let the social utility function be U (c) and the social disutility of pollution be D(x). The goal is to choose c(t) over time to maximize the expected welfare where λ > 0 is a discount factor. For example, we can take U (c) = c β with 0 < β < 1 and D(x) = x 2 .
3. HJB equations. Consider the objective (2) with constraint (1). Denote the corresponding value function by where x ∈ R d , α ∈ {1, 2, . . . }. We obtain the following system of dynamic programming equations where the operator L is defined as: (7) Assume that for each α, v(·, α) has continuous partial derivatives up to the second order, that Q = (q ij ) is the generator of the Markov chain α(·), and that x is a d-dimension vector with x = (x 1 , x 2 , . . . , x d ) , and z denotes the transpose of z ∈ R d . 4. Numerical method. To solve equation (6), a closed-form solution is difficult or virtually impossible to obtain due to the Markov switching. We need to deal with a system of equations instead of a single equation. As a viable alternative, we aim to obtain numerical solutions.
First, using the finite difference method with step size h > 0, we approximate the derivatives as: is the standard basis in R d . It follows that we can rewrite (6) as an iterative algorithm: and a + and a − denote the positive and negative parts of a, respectively. We denote the operator on the right-hand side of equation (11) as T .
It can be proved that the operator T is a contraction mapping. We refer the reader to Barles and Jakobsen [1] and Lu et al. [12] for details. Then by the Banach fixed-point theorem, we conclude that there exists a unique V h satisfying V h = T V h .

5.
Convergence of Markov chain approximation. By Theorem 4.1, we need only demonstrate that this limit gives a good approximation to the value function of the original problem (5). In this section, we first that show the limit in theorem 4.1 coincides with the value function for a constructed controlled Markov chain, denoted by V h M C . Then we demonstrate that V h M C approximates the value function V in equation (5). 5.1. Markov chain approximation. We propose to simulate the controlled stochastic process (1) by a controlled Markov chain and study underlying relation between our numerical approximation and the optimization problem under the constructed Markov chain.
To begin, we consider a discrete Markov chain with state space . . } and denote the temporal difference as ∆t. Furthermore, we let C h n be a sequence of control actions, which is adapted to the smallest σ-algebra generated by To retain the local behavior of the original process (similar to the consideration of local consistency defined in Kushner and Dupuis [10, p.71]), we define the transition probabilities of the Markov chain {η h i , α h i } given control action C by: We proceed to verify that our constructed approximating Markov chain is locally like the switching diffusion (1). Due to the presence of α(t), we need to modify the local consistency definition in [10]. The definition is given as follows.
Definition 5.1. The sequence (η h n , α h n ), with temporal difference ∆t(x, i, C), is said to be locally consistent with equation (1) if The definition above reflects the idea that in a small interval ∆t(x, i, C), the constructed approximating Markov chain is in line with the switching diffusion in that it has the same local properties as (1). Using the above definition, direct computation leads to the local consistency of our approximating Markov chain. For simplicity, the detailed calculation is omitted.
and use continuous-time interpolations Using the interpolation of the controlled Markov chain {(η h (·), α h (·)), C h (·)}, we consider an optimal control problem with the objective function as (2). The corresponding problem is formulated as: which can be regarded as a Markov decision process (MDP) Khalvati et al. [7]. Through dynamic programming principle, we can obtain the iterative form of the value function as Note that the operator on the right-hand side agrees with that in (11). By the uniqueness of fixed point, we conclude that for any x > 0, . This shows that the value function obtained by our numerical algorithm can be considered as the value function of a related problem by solving an Markov decision process based on a simulated process (controlled Markov chain).

5.2.
Weak convergence to controlled switching diffusions. This section demonstrates that appropriately interpolated process converges weakly to the controlled switching diffusion. First, the process η n may be written as We denote which is a martingale w.r.t. the smallest σ-algebra generated by ). Using the local consistence property, we can interpret its continuous-time interpolation as a stochastic differential equation with driving force W h (·), and o h (t) → 0 in probability as h → 0.
To proceed, it is more convenient that we use a relaxed control representation [10, Section 9.5], denoted by m h (·). Using the relaxed control to represent the control actions C h (·), leads to the use of a nice compactness property in the analysis. To proceed, we denote the set of all admissible controls by C. We map the set C into a larger space with m h (·, ·) defined on B(C × [0, ∞]), where B denotes the Borel sets, which satisfies: where F t is the smallest σ-algebra generated by {η h (s), α h (s), s ≤ t}. Moreover, we can further rewrite the process η h (·) as where we have used m h s (dc), the "derivative" (see [10, p.267] and also Song et al. [13]).
Proof. We define F h t as the smallest σ algebra generated by Denote the corresponding conditional expectation on F h t by E h t . For any ε > 0, and t, s > 0 satisfying s ≤ ε, it can be shown that (see, for example, using similar techniques in [14, pp. 218-220]) for some γ(h, ε) > 0, The process α h (·) is thus tight by Kushner [9, Theorem 3, p.47].
Thus we obtain the tightness of {W h (·)}. Likewise, we can show η h (·) is tight. Next, note that the measures m h (·) live in a compact set, thus it is tight. Putting the estimates together, the desired result of the lemma follows.
Proof. By the tightness and the Prohorov theorem, there is a convergent subse- For notational simplicity, we still index the subsequence by h with limit denoted by {η(·), α(·), m(·), W (·)}. We first show that W h (·) converges weakly to W (·), which is a standard Brownian motion. For ε > 0 and any process y(·), define the process y ε (·) by y ε (t) = y(iε), t ∈ [iε, iε + ε). We note that for any t, s > 0, where E|e h ε (t)| → 0 as ε → 0 uniform in t in any bounded time interval. To characterize W (·), let t > 0, δ > 0, κ 1 , κ 2 , {t k : k ≤ κ 1 } be given such that t k ≤ t ≤ t + r for all k ≤ κ 1 , φ j (·) for j ≤ κ 2 is real-valued and continuous functions on U × [0, ∞) with compact support for all j ≤ κ 2 . Define Let H(·) be a real-valued and continuous function of its arguments with compact support. Note that W h (·) is an F h t -martingale. Thus we obtain Using the Skorohod representation and the dominant convergence theorem, passing to the limits as h → 0, we have Since W (·) has continuous sample paths, W (·) is a continuous F t -martingale. Note that By the weak convergence and the Skorohod representation theorem, we have The quadratic variation of the martingale W (t) is t, which implies W (·) is an F t -Wiener process. Then, by the tightness of {η h (·), α h (·)}, we can further get from (21),  E|e h ε (t)| → 0, uniformly on any finite time interval.
Letting h → 0, the convergence with probability one (using Skorohod representation) yields uniformly in t. In addition, {m h (·)} converges in the "compact weak" topology. For any bounded and continuous function φ(·) with compact support, The weak convergence and the Skorohod representation imply that [Gc + Ξ(α(s))X(s)]m s (dr)ds, uniformly in t on any bounded interval with probability one (using Skorohod representation). Because η h ε (·) and α h ε (·) are piecewise constant functions, it follows from the probability one convergence, Combining the estimates obtained thus far, [Gc + Ξ(α(s))X(s)]m s (dc)ds where lim ε→0 E|e ε (t)| = 0 uniformly on any bounded interval. Letting ε → 0 in the above leads to the desired result.

Convergence of objective and value functions.
It remains to demonstrate that V h M C (x, α) approximates V (x, α), the value function of the original problem.
, C h (·)) → J(X(·), α(·), C(·)). Proof. In view of the weak convergence and the Skorohod representation theorem, the desired result follows by noting the continuity and hence the boundedness of the objective function.
follows directly from the weak convergence in Theorem 5.4. For any ε > 0, We need to prove is the optimal control for with trajectory X ε (·) and relaxed control representation m op (·). By the chattering lemma [10], there exists an ordinary control C δ op (·), which is constant on [kδ, (k + 1)δ], and only takes finite values such that Then we consider optimal control problem of the process X ε,∆t (t) = X ε (k∆t) for t ∈ [k∆t, (k + 1)∆t], where k = 0, 1, 2, . . . , which can be considered as a piecewise constant approximation of the process X ε (·). Denote the optimal control asĈ(·), with trajectoryX ε (·) and relaxed controlm(·). It follows: By the chattering lemma again, there exists ∆t > 0 and an ordinary control C ∆t (·), which takes finite values and which is constant on [k∆t, (k + 1)∆t], for k = 0, 1, 2, . . . such that It followsĈ ∆t (·) is admissible with respect to the controlled Markov chain defined in Section 5.1. Using lemma 5.5 we got The proof is complete. 6. Numerical examples. In this section, we use several examples to illustrate our numerical approximation methods. Example 6.1. We first consider a one-dimension case given in (3) and (4). In particular, we take D(x) = 100x,Ũ (c) = 2 √ c, and σ 0 (x, i) = σ(i)x. Suppose that the Markov chain has two states, denoted by 1 and 2. The generator of the Markov chain is given by

CAOJIN ZHANG, GEORGE YIN, QING ZHANG AND LE YI WANG
Let λ = 0.05 with pollution generation/decay rates ρ(1) = 0.7, ρ(2) = 0.3, the noise "volatility" σ(1) = 0.1, σ(2) = 2.3, and consumption rate c(t) ∈ [0, 3]. We use our numerical scheme to obtain the optimal control actions and corresponding value functions. The results are depicted in Figures 1 and 2, respectively. As it can be seen that the optimal controls are of the threshold type. In state 1, if the accumulation of pollution is less than x = 63, we can perform full consumption rate. If the stock exceeds x = 63, we should stop consuming. In state 2, the threshold value is x = 80. This value is greater than that of state 1. It is due to the greater decay rate ρ(2) in state 2. We obtain the following result.  Fig. 3.
The thresholds are smaller than those in Example 6.1. The reason is that a quadratic disutility function has a much larger rate of increasing with respect to x, the stock of pollution. With control actions, although it will bring the benefits from utility function U (·), the disutility part increases more dramatically with the accumulation of pollution. As a tradeoff, the optimal decision is not to implement the control when x > 41 in state 1 and x > 69 in state 2. Example 6.3. Next, we try to find the optimal strategy from a real-data scenario. We use real air pollution data provided in De Vito et al. [2] to calibrate our regimeswitching model. In particular, it records the accumulation of 13 materials that   Fig. 4. For experimental purpose, we split the data set into two parts. The first part contains the data from 03/10/2004, 18:00 to 11/15/2004, 17:00 and is used for model calibration and for derivation of the optimal strategy, denoted as the training set D train . The second part is from 11/15/2004, 18:00 to 04/04/2005, 14:00, denoted as the testing set D test . We assume D train is historical data and D test is data collected in future which is not observed by the model. We calibrate the model using D train and apply the model implied consumption strategy on D test to test the model's generalization ability on new samples.  In this example, we analyze the accumulation process of NO 2 . Using an Expectation Maximization (EM) based algorithm combined with a Bayesian Inference method as Zhang et al. [18], we fit the hourly observations of NO 2 in D train into a one-dimensional regime-switching diffusion model (1) with 3 underlying states.
We obtain ρ(0) = 0.0078, σ 2 (0) = 0.043 for the 0th hidden state, ρ(1) = −8.609E, Eσ 2 (1) = 0.566 for the first hidden state, and ρ(2) = 8.422, σ 2 (2) = 0.640 for the second hidden state. Also, the estimated generator is Given the utility function: with the control set c ∈ [0, 1]. We consider optimal consumption plan as to the maximize the utility function under the pollution accumulation process. Using our numerical method, we demonstrate the value functions and corresponding optimal controls in Figure 5, respectively. Using our numerical method, the optimal controls are obtained as Given the same utility function, we obtain the optimal controls as  Example 6.5. We consider a two-dimensional system. As can be seen that the processes NOx and NO 2 are highly correlated. It is natural to integrate them as a 2-dimension correlated process of (1) with a full covariance matrix. Consider D train for N O 2 and N O x : {x 1 1 , . . . , x 1 N }, {x 2 1 , . . . , x 2 N }, respectively. It is known that combustion of natural gas and coal are two important industrial sources to produce pollution materials N O x and N O 2 . We consider control variable C = (c 1 , c 2 ), corresponding to the consumption of natural gas and coal, respectively. We assume that the consumption of per unit natural gas generates 0. We further assume the pollution materials can be described by which is known as the system of mixture regression model, where ∆t is the length of each time interval, and 1 , 2 ∼ N (0, 1) are two independent standard Gaussian noises.
Using the EM based algorithm, we calibrate the mixture regression model (24) with data from D train and obtain parameter estimates as Based on these relations, we obtain the optimal control policy given in Fig 8 (b), (d), and (f). Consider the utility function given by where x 1 , x 2 represent the state variables of NOx and NO 2 concentration, respectively, and c 1 , c 2 ∈ [0, 1] are the corresponding control variables. We solve the maximization problem and obtain 4 optimal control actions: (c 1 , c 2 ) = {(0, 0), (1, 0), (0, 1), (1, 1)} indexed as 0,1,2,3. We then obtain the value functions and optimal actions in Next, we examine the performance of our optimal strategy on the testing data. Using a K-Nearest-Neighbor (KNN) Larose [11] algorithm, we estimate the optimal control actions on the testing data set and calculate the accumulative utility: −114.25. Then we use the dimension reduction by X = XH to obtain a covariate matrix R 6000×2 and fit it into a 2-dimension regime switching diffusion model as (24). We obtain the parameter estimates as: We consider the objective function To proceed, we apply our numerical method to derive the optimal controls C * ( X) = (c 1 , c 2 ) with respect to the state space X. There are 4 optimal control actions: (c 1 , c 2 ) = {(0, 0), (1, 0), (0, 1), (1, 1)} indexed as 0,1,2,3. We demonstrate the value functions and optimal actions with respect to the state space X as in Figure 9.
We compute the optimal controls on D test , by implementing the data transformation X = XH and exploiting our optimal controls C * ( X). The optimal control on D test is demonstrated in Fig. 10. The accumulative utility is −415.37, which is much better than the loss −1503.33 by using a random policy.
7. Further remarks. This paper is devoted to numerical solutions of a pollution control problem. While this study is limited to a simplified model structure, more general model structures can be considered such as nonlinear models in the continuous state x system. In addition, to balance the future and current generations, more advanced criteria [4] can be treated. These will be pursued in our future studies.