MODELING AND COMPUTATION OF ENERGY EFFICIENCY MANAGEMENT WITH EMISSION PERMITS TRADING

. In this paper, we present an optimal feedback control model to deal with the problem of energy eﬃciency management. Especially, an emission permits trading scheme is considered in our model, in which the decision maker can trade the emission permits ﬂexibly. We make use of the optimal control theory to derive a Hamilton-Jacobi-Bellman (HJB) equation satisﬁed by the value function, and then propose an upwind ﬁnite diﬀerence method to solve it. The stability of this method is demonstrated and the accuracy, as well as the usefulness, is shown by the numerical examples. The optimal management strategies, which maximize the discounted stream of the net revenue, together with the value functions, are obtained. The eﬀects of the emission permits price and other parameters in the established model on the results have been also examined. We ﬁnd that the inﬂuences of emission permits price on net revenue for the economic agents with diﬀerent initial quotas are quite diﬀerent. All the results demonstrate that the emission permits trading scheme plays an important role in the energy eﬃciency management.

1. Introduction. Although energy consumption is an essential element in economic development, we are also becoming increasingly aware of the negative impacts of energy use. The anthropogenic greenhouse gas (GHG) emission has risen dramatically during the last few decades, which mainstream researchers believe to be the main cause of climate change, especially the global warming. According to International Energy Agency (IEA) (2006) [1], for the sake of meeting the energy demand, the carbon dioxide emitted mainly by burning fossil fuels accounts for around 80% of the global emission. Therefore, it is necessary for us to realize the importance of reducing the GHG emission through improving energy use.
One way to improve energy use is increasing the primary energy efficiency. It is possible to burn less fossil fuel to obtain the same energy-consuming products by developing technologies of minimizing heat loss and maximize energy consumption efficiency. In other words, these technologies can be used to minimize the energy consumption per product in manufacturing. Higher energy efficiency can lower the operating costs and make a company more competitive in local and global markets. In addition, at a national level, the benefits from higher energy efficiency can reduce the dependence of a country on imported energy and extend the energy reserves.
Because of its importance, more and more researchers have paid their attention to energy efficiency [9,[11][12][13]20]. Especially, in [24], Worrell et al. discuss the potential contribution of industrial energy-efficiency technologies and policies to reduce energy use and greenhouse gas emissions to 2030. Patterson in [18] critically reviews the range of energy efficiency indicators that can be used, particularly at the policy-making level. The specific limitations and appropriate uses of physicalthermodynamic, economic-thermodynamic and pure economic indicators of energy efficiency are also considered. Besides, to discover the mechanism of international technology diffusion, Jin and Zhang develop an efficiency-improving vertical innovation model where energy technological progress is specified as an improvement in primary energy use efficiency in [14].
However, to our best knowledge, there are very few studies on energy efficiency management which take emission permits trading into consideration. For the past few years, the issues on climate change and emission reduction have attracted the attention of politicians and scholars from all over the world. In order to mitigate climate change and improve global environment, some emission permits trading markets have emerged and are developing prosperously. At the same time, a large number of articles have studied the emission permits prices theoretically and empirically [4,5,7,10,21]. Moreover, some problems about production decision also include the emission permits trading scheme [2,3,15,25].
Motivated by the those mentioned above, in this paper, we present an optimal feedback control model, in which the revenue from emission permits trading are included, to study the energy efficiency management. In our model, the decision maker should choose proper scale in the investment on undertaking indigenous innovation and absorbing foreign knowledge diffusion, to maximize the net revenue stream. Additionally, an initial quota from the emission permits trading scheme is given, and the emission permits can be purchased from the market compulsively if the quota is insufficient and the unused emission permits can be also sold to others.
By means of optimal control theory, we derive a Hamilton-Jacobi-Bellman (HJB) equation satisfied by the value function. Since this equation cannot be solved analytically, we propose an upwind finite difference method, which is based on an explicit finite-difference scheme, to solve it. In computational physics, an upwind scheme is a method for solving hyperbolic partial differential equations numerically. The upwind scheme attempts to discretize hyperbolic partial differential equations by using a different differencing scheme determined by the sign of the characteristic directions. It has been used in many fields and is becoming more and more popular. See, for instance, [17] for degenerate parabolic problems, [16] for hyperbolic problems, and [19] for elliptic problems.
The paper is organized as follows. In Section 2, we present the optimal feedback control model governing the energy efficiency management and make use of optimal control theory to derive an HJB equation satisfied by the value function. Then, an upwind finite difference method is proposed for the discretization of the HJB equation in Section 3 in which stability of the numerical method is also analysed. In Section 4, some numerical examples are presented to illustrate the efficiency and usefulness of the numerical method, and the effects of the parameters on the value function as well as the optimal management strategies are also showed in this section. Finally, concluding remarks are given in Section 5.
2. Energy efficiency management model. For the purpose of illustrating the energy efficiency management problem in a commitment period, we propose a finitehorizon optimal feedback control framework. Also, we assume that the decision maker does not take external economic environment into account during the decision making to maximize his own revenue. We do believe that this assumption can be expanded into a differential game framework in the future works.
Following [14], we denote the outputs of end-use by Y (t), and the secondary energy products or services produced by the economic agent with an aggregate production function is given by: where K(t) and E(t) are the inputs of capital (power generating machines and equipments) and the primary energy resource (coal, oil, natural gas) at time t, respectively, and A(t) is the efficiency of converting primary energy to end-use, secondary energy products. Without loss of generality, we assume that the production function takes the conventional Cobb-Douglas form [6]: where α and β are the elasticities of capital stocks and energy resource, respectively. We suppose that the energy efficiency A(t) is a time-varying variable and is determined by two factors, namely indigenous innovation and foreign technology diffusion. The law of motion of A(t) is given by the following controlled process: where A W T F denotes the level of energy use efficiency of the world technology frontier, λ and σ are the control variables representing the level of undertaking indigenous innovation and absorbing foreign knowledge diffusion, respectively. In (2), the term λ · A(t) implies that the energy efficiency in the next period should be more affected by the current level of energy efficiency. While the term σ · (A W T F − A(t)) implies the effect of technology spillover. Decision makers need to find the best policy λ and σ to maximize the net revenue stream. Consequently, the capital stock accumulates according to the following dynamics: where δ denotes the rate of depreciation. In the part of emission permits trading scheme, the trading volume of emission permits at price S is given by where E 0 denotes the initial quota, when τ > 0, the economic agent purchases the emission permits from the market, and τ < 0 means that the economic agent sells the unused emission permits to others.
Consequently, if we choose a logarithmic utility function, the economic agent's revenue involving emission permits trading at time t can be represented as Furthermore, we suppose that the cost of undertaking indigenous innovation and absorbing foreign knowledge diffusion are 1 2 λ 2 and 1 2 σ 2 , respectively. The current objective of economic agent is to find an optimal plan which maximizes the present flow of instantaneous net revenue. That is, the objective functional and constraints of economic agent are as follows: where r > 0 is the social risk-free discount rate, t = 0 is the initial time, T is the terminal time, and A 0 and K 0 are given initial conditions. The optimal control model (5) can be solved with the help of HJB equation. According to [8], we can derive an HJB equation satisfied by the value function with the terminal condition V (A, K, T ) = 0.
3. Numerical method. In this section, we will present an upwind finite difference method to discretize (6) since it is difficult to be solved analytically. We will also prove that the numerical method developed is stable. We first denote the optimal management strategies by λ * , σ * and E * . From the first-order optimality condition, we know that (6) can be split into the following coupled equations: and E * solves the following nonlinear equation where respectively.
3.1. Upwind finite difference method. To implement the upwind finite difference method, we first truncate the whole space into the finite computational domain , and then uniformly divide the intervals into N A and N K sub-intervals, respectively: for l = 0, · · · , N t , is the partitions along the t axe and let ∆t , whose mesh lines are perpendicular to one of the axes, is defined.
Using the above mesh, we can approximate (8)-(11) by the following difference equation: and (E * ) l−1 i,j solves the nonlinear equation for i = 1, 2, · · · , N A − 1, j = 1, 2, · · · , N K − 1, and l = 1, · · · , N t , where sign z denotes the sign of z, and . This is the upwind finite-difference scheme, which reduces to the forward-difference scheme when the coefficients are positive, and to the backward-difference scheme when the coefficients are negative. Then, we propose the following algorithm based on the explicit scheme (12) for the numerical solution of (6)- (7).
i,j and (E * ) 0 i,j for i = 1, 2, · · · , N A − 1 and j = 1, 2, · · · , N K − 1 using the terminal condition (7); • For l = 1 to N t compute V l i,j , (λ * ) l i,j , (σ * ) l i,j and (E * ) l i,j as follows. For i = 1, 2, · · · , N A − 1 and j = 1, 2, · · · , N K − 1, set and (E * ) l i,j solves the nonlinear equation Remark 1. According to [23], a hyperpyramidal region, which contains the original hypercubic one as a subset, is needed to solve the problem due to the explicit nature of the scheme (16). In the numerical experiments, we will implement the above algorithm in a hyperpyramidal region.

Remark 2.
To solve the nonlinear equation (19), we use the standard Newtonian method in which (E * ) l−1 i,j is chosen as the initial guess.

3.2.
Stability of the discretization. Since the numerical scheme (12) is explicit in the time-direction, we need to discuss its stability. The stability of (12) is established in the following theorem.
Proof. To begin with, we can know from [22] that (12) is stable only if (12) with F = 0 is stable. In the following discussions, we will derive the sufficient condition for the stability of the numerical scheme.
If we set F = 0, then (12) becomes . Squaring both sides and summing over i and j, we obtain So, a sufficient condition is The stability of the scheme under the condition (20) guarantees that the error due to the discretization at each time step will not grow in the time-stepping procedure. Thus, the error is dominated by the approximation of the derivatives in space. In the next section, we will numerically demonstrate the convergence of the scheme.
4.1. The efficiency of the numerical method. First we consider the convergence rate of our discretization method. Since the closed-form solution of the HJB equation (6) cannot be found, we use the numerical solution on the uniform mesh with mesh nodes N A = 512 = N K and N t = 128 as the "exact solution". We compute the errors in the discrete L ∞ -norm at the computational final time step t = 0 on a sequence of spatial meshes with N A = N K = 2 n+2 and N t = 2 n for a positive integer n from n = 1 to a maximum n = 6. The discrete L ∞ -norm is defined as: where V h denotes the numerical solution. The log-log plots of the computed maximum errors, along with the linear fitting, are depicted in Figure 1. From the figure we can see that the rate of convergence of V h in the discrete L ∞ norm is of the order O(h 1.1398 ), where h = max(∆A, ∆K). Additionally, this demonstrates numerically that our numerical method for (6) governing energy efficiency problem is useful and efficient. A rigorous mathematical proof of this rate of convergence will be presented in a future paper.

4.2.
The solution of the model. In this and next parts, we will illustrate the results by tables and figures. We plot the value function V , indigenous innovation level λ, absorbed knowledge level σ and emission level E in Figure 2 and list some of their values in Table 1. In what follows, we will explain our results, particularly the trends of the computed results in energy efficiency A and capital stock K, qualitatively, rather than quantitatively. From the third row of Table 1 and Figure 2(a) we see that the value function is an increasing function of both the energy efficiency and the capital stock. This implies that a higher energy efficiency, as well as a more capital stock, can result in a higher net revenue of the economic agent. This is a reasonable phenomenon. On the one hand, the economic agent can save more fossil fuels in the production process when the energy efficiency is higher. On the other hand, more capital stocks can results in more output, and thus more revenues.
Then, the results about the indigenous innovation λ and the absorbed knowledge σ, which are plotted in Figure 2(b) and Figure 2(c) and listed in Table 1's fourth row and fifth row, show that with the increase of energy efficiency, the economic agent would like to increase the level of undertaking indigenous innovation rather than the level of absorbing foreign knowledge diffusion. That is to say, when the energy efficiency is higher and approaching the level of world technology frontier, it would be more difficult to absorb foreign knowledge, and the economic agent would be more likely to invest on indigenous innovation. Moreover, the capital stock seems to have few effects on the decision about undertaking indigenous innovation and absorbing foreign knowledge diffusion.
In addition, we can clearly see from Figure 2(d) and Table 1's last row that the emission level increases as energy efficiency increases and decreases as capital stock increases.. This implies that the economic agent takes advantage of the higher energy efficiency to increase the emission level for more production revenues. In another aspect, increasing capital stock, combined with decreasing emission level, results in an invariant production revenue as well as more gains from emission permits trading.
Last but not least, we plot a phase portrait of energy efficiency A v.s. capital stock K in Figure 3. From the trend of horizontal axis about vertical one, we can conclude that the optimal capital stock increases with the increase in optimal energy efficiency, which implies that the two state variables in the optimal control system are positively correlated in the time horizon.  4.3.1. The effects of parameters S and E 0 . In this and next two subsections, we will show the effects of parameters on the results by presenting some figures, in which t = 0 and K = 0.5. Since our main work in this paper is to integrate emission permits trading into energy efficiency management, we first examine the effects of emission permits price on the results by Figure 4. Note that S is set to be 0.3, 0.5, and 0.7 in each figure.
We can see from Figure 4 that a higher emission permit price results in more revenue, a lower indigenous innovation level as well as a lower emission level. Furthermore, the absorbed knowledge level is not sensitive to the emission permit price. In this example, the initial quota E 0 is set to be larger, and the sum of revenues from emission permits trading and saved cost from indigenous innovation is more than the loss in production due to the emission reduction.
However, if we set the initial quota E 0 to be 0.2, the effects of S on the results do not change besides of the one on value function. We plot the value function for the different values of S when E 0 = 0.2 in Figure 5. On this occasion, the revenue reduces with the increase of emission permits price. That is to say, the sum of revenues from emission permits trading and saved cost from indigenous    innovation is less than the loss in production due to the emission reduction. So, we can conclude that the economic agent with higher initial emission quota prefers to higher emission permits prices.  In addition, the influence of initial quotas, which are set to be 0.5, 0.7 and 0.9, on the results are also plotted in Figure 6. The initial quota can be regarded as the inherent revenue and its allocation rule should be based in part on historical data and current actual capabilities. We find that the initial quota only influences the net revenue without any changes in other three control variables, and the more the initial quota is, the greater the net revenue is.

4.3.2.
The effects of parameters α and β. In this part, the effects of the elasticities α and β on the results are studied. In Figure 7, α is set to be 0.6, 0.8 and 1, and in Figure 8, β is set to be 0.1, 0.2 and 0.3. According to the production function theory, if α+β = 1, the production function has constant returns to scale, which means that doubling the usage of capital and energy will also double output. If α+β < 1 returns to scale are decreasing, and if α + β > 1 returns to scale are increasing. Following the increase of α and β, the returns to scale are also change from decreasing to increasing, so the indigenous innovation level, absorbed knowledge level and the emission level all increase to pursue more revenues as depicted in Figures 7(b), 7(c), 7(d) and in Figures 8(b), 8(c), 8(d).
However, in the case that β increases, the gains from the expanded production can not offset the added costs in increasing the level of undertaking indigenous innovation and absorbing foreign knowledge diffusion and the losses in emission permits trading, which leads to the difference in the changed trends of net revenue with respect to α and β.   Figure 9 shows the effects of depreciation rate of capital stock on the results. In the figure, δ is set to be 0.3, 0.5 and 0.7. Similar to initial quota E 0 , depreciation rate δ only has effect on the net revenue. Moreover, the larger the depreciation rate is, the more the capital invested by the economic agent is. Therefore, the net revenue should be decreasing with increasing the depreciation rate.

4.3.4.
The effects of parameters A 0 and K 0 . Since the initial values of energy efficiency A 0 and capital stock K 0 do not influence the optimal revenues and optimal control variables, in this part, we only examine their effects on the phase portrait of A v.s. K. In Figure 10(a), A 0 is set to be 0.3, 0.4 and 0.5, and in Figure 10(b), K 0 is set to be 0.3, 0.4 and 0.5, respectively. From the results, we can see that the relative trend of state variables A and K will not change with the initial values.
5. Concluding remarks. In this paper, we present an optimal feedback control model, in which an emission permits trading scheme is involved, to deal with the problem of energy efficiency management. We make use of the optimal control theory to derive a Hamilton-Jacobi-Bellman (HJB) equation for the value function, and then propose an upwind finite difference method to solve it. The stability of the proposed numerical method has also been proved. Extensive numerical experiments have been performed to show the accuracy and usefulness of both the model and numerical method. Using the computed optimal solutions with various value of   system parameters, we have examined the effects of the emission permits price and other parameters in the established model on the optimal strategies. We find that the influences of emission permits price on the net revenue for the economic agents with different initial quota are quite different. Future works can be performed from two aspects. Mathematically, we will design a new efficient numerical scheme, such as discontinuous Galerkin finite element methods or finite volume methods, to solve the derived HJB equation. From the viewpoint of economics and management science, our model can be extented to a differential game framework, in which the decision maker should take external economic environment into account. Furthermore, the influence of pollution resulting from emission can be also considered.