Taxing Strategies for Carbon Emissions: a Bilevel Optimization Approach

This paper presents a quantitative and computational method to determine the optimal tax rate among generating units. To strike a balance between the reduction of carbon emission and the profit of energy sectors, the proposed bilevel optimization model can be regarded as a Stackelberg game between the government agency and the generation companies. The upper-level, which represents the government agency, aims to limit total carbon emissions within a certain level by setting optimal tax rates among generators according to their emission performances. The lower-level, which represents decision behaviors of the grid operator, tries to minimize the total production cost under the tax rates set by the government. The bilevel optimization model is finally reformulated into a mixed integer linear program (MILP) which can be solved by off-the-shelf MILP solvers. Case studies on a 10-unit system as well as a provincial power grid in China demonstrate the validity of the proposed method and its capability in practical applications.


Introduction
Global climate change is mostly attributed to the excessive emission of greenhouse gases, especially carbon dioxide (CO 2 ). According to [1], the concentration of CO 2 in the Earth's atmosphere has grown by 35% since the Industrial Revolution. The emissions of CO 2 will still increase by 1.9% annually if no effective action is taken. The strong requirement for a cleaner environment and sustainable development brings new challenges to all emission intensive sectors.
In order to combat global warming issues, many countries have launched their future plans in emission reduction and devoted substantial efforts to achieve their goals through some market-based policies. Such policies regulate customers' behavior through market forces. Briefly speaking, the emission of CO 2 would be reduced if the price of emission intensive products increases and consequently the demands shift to other alternative environmentally-friendly goods. From an economic point of view, market-based policies effectively achieve the target of emission reduction at the lowest cost to society [2]. Among various market based policies, the cap-and-trade (C&T) program [3] and carbon tax scheme [4] are the most widely implemented regulations.
In the C&T program, a government agency specifies a limit or cap on the total amount of carbon emissions. Then the cap is allocated to entities under regulation in the form of carbon credits. Generation companies (GENCOs) whose carbon credits are more than their actual emissions are allowed to trade them [3]. In the trading market, the buyer is charged for polluting the environment, while the seller is rewarded for protecting the environment. In practice, the European Union Emissions Trading System, known as the first and largest emissions trading system in the world [5], has become the key tool for reducing greenhouse gas emissions in Europe. The workings of spot, forward, future and option markets are studied in [6], the economic analysis is carried out in [7] and examined among 12 European countries. In the United States, a national market for SO 2 trading and several regional markets for nitrogen oxides trading have been established since 1990. Other countries such as New Zealand and Japan have also set up their own CO 2 trading programs.
The carbon tax is a direct charge imposed on carbon emissions. GENCOs are obligated to pay fees proportional to the quantity of their emissions. In the early 1990s, Sweden became the first country attempting to manage the carbon emission with the carbon tax scheme. In the following 20 years, a number of countries have tried to implement carbon tax or energy tax schemes that aim at reducing carbon emissions [8][9][10] and obtained positive effects. In 2009, the Chinese Ministry of Finance, together with the National Development and Reform Commission, promulgated a research report on the feasibility of implementing carbon tax policy in china, and announced future plans in the coming decade.
There has been a hot debate between C&T and carbon tax scheme [11][12][13][14]. Some argued that the carbon tax is easy to carry out since it does not require trading platforms. Moreover, it provides a price signal that stimulates emitters to seek long-term solutions such as investing on renewable generation sources [12]. But others insists that C&T has some distinctive advantages over a carbon tax, such as international fairness, broad participation, and most important, direct regulation on the total quantity of emissions [11]. It's worth mentioning that a comprehensive bilevel optimization model was proposed in [14] to investigate and compare the C&T and carbon tax policies in terms of effectiveness and efficiency from a generation expansion planning perspective. Relative advantages and disadvantages of both policies with respect to different criteria were revealed.
Despite the debate about C&T and carbon tax policy, a couple of strategic approaches have been studied to reduce carbon emission from energy and manufacturing industry. An inexact mix-integer two-stage program model was proposed in [15] for long term decision making problems under uncertainty, such as capacity expansion, facility improvement as well as coal inventory planning. A low-carbon dispatch model of wind-incorporated power systems is proposed in [16], in which the bi-objective short term dispatch problem is converted into a single nonlinear program using fuzzy satisfaction-maximizing approach. In [17], decision architecture for tri-level production process is proposed for implementing "environmental-benign" and low emission manufacturing industry. Optimization based approaches for low-carbon energy systems planning under uncertainty are comprehensively reviewed in [18]. In traditional single level optimization methods, emission constraints are incorporated in the economic dispatch (ED) model, so that the actual emission corresponding to the optimal solution will not violate the pre-specified value. However, due to lack of legal means to regulate the ED decisions, speculative operators might relax the emission constraints. As a result, the total emission cannot achieve the desired level. The carbon tax is a direct legal mean to supervise the operator's behaviors. A method is proposed in [19] to design the tax rate based on game theory. It was pointed out in [19] that allowing the tax rate to vary among manufacturers according to their corresponding energy efficiency is more effective than adopting a uniform tax rate.
Inspired by the work in [14,19], in this paper a bilevel optimization model to determine the optimal carbon tax rate is proposed. In order to capture the hierarchical relationship between the target of emission reduction and the profit of energy industry, the model can be viewed as a Stackelberg game whose players include the government and the GENCOs. The government sets tax rate in the upper level while considering reactions from the grid operator of the lower level. The emission constraint is explicitly considered in the upper level. Different from [14], the objective of the upper level authority is to minimize the actual total tax levied on GENCOs, which is a nonlinear function and finally linearized by using primal and dual variables of the lower level problem. Different from [19], the power system dispatch decision is explicitly modeled in our formulation. Note that in our settings, each generator has a different tax rate. This differs from the usual practice of setting a single tax rate or carbon price across all generators or even across all industrial sectors or emitters which is the first best economically efficient policy. The reason for this difference is that, as explained above, the objective of the regulator is to minimize the total tax revenue given the emissions constraint, which differs from the usual objective of minimizing economic cost given the emissions constraint.
The reminder of this paper is organized as follows: the mathematical formulation is provided in Section 2. The solution method is developed in Section 3. Test results on a 10-unit system as well as a provincial power gird in China are given in Section 4. Finally, conclusions are given in Section 5.

Mathematical Formulation
In this paper, the problem is established against the following background. A certain number of generators are operated by a control center, say an Independent System Operator (ISO), or a regional dispatch center in China. The carbon emission of each generator is assumed to be proportional to its output. Aiming at limiting the carbon emission within a pre-specified level in a target year, the government determines a tax rate for carbon emission on each generator. In response to the tax signal, the ISO executes an ED program by taking into account the carbon tax cost in its objective function and determines the electricity production of each generator.
Different from C&T programs directly setting the quantity of carbon emission, the carbon tax scheme indirectly controls the emission level through prices. Therefore, the most critical problem confronting the tax decision maker is to determine an appropriate tax rate that achieves the target of emission reduction in a desired period. If the carbon tax rate is too high, the GENCOs may suffer from financial burdens, and the economic growth will be hindered. On contrary, if the carbon tax rate is too low, the actual emission cannot be reduced to the desired level. Therefore, it is important for the tax decision maker to strike a proper balance between the environment protection and economy growth.
A distinguished feature of the proposed model is it captures the reactions from the ISO to the carbon tax signal from a lower level optimization problem, so the quantity of emission corresponding to the ISO's dispatch decision can be easily estimated by the tax maker. This feature overcomes the traditional difficulty of the carbon tax policy. In the following subsections, a load model is firstly introduced to simplify the ED problem of ISO, and then the bilevel optimization model of the taxing problem is established.

Load Modeling
Similar to [20], the yearly load curve is approximated through a piecewise constant duration curve containing several demand blocks and shown in Figure 1. This manifestation can be regarded as rearranging the yearly time sequence load curve into a quantity descending sequence, and merging periods with similar load levels into a single demand block. This approach is appropriate for a long term policy making problem. However, this modeling usually does not allow accurately representing time coupled operating constraints such as minimum up/down time constraints and ramping constraints of conventional units.

Bilevel Optimization Model
The parameters and variables used in the model are defined as follows: The tax policy making problem can be formulated as a bilevel optimization, and regarded as a Stackelberg game, in which the government is the leader, whose strategy is the tax rate, and ISO is the follower, whose strategy is the output of generators. The upper-level (UL, the leader's problem) represents the tax rate decision process of the government agency with the target of minimizing the total levied tax subject to restricting the emission within a certain level. The lower-level (LL, the follower's problem) represents the ISO's ED problem subject to operating constraints with the tax rate fixed. The bilevel optimization model is provided as follows: Note that in the lower-level problems (2), dual variables are provided at the right side of their corresponding constraints following a colon, they are also decision variables following the paradigm of complementarity modeling approach [21].
The UL problem represents the total tax minimization for the government, constrained by: (i) tax rate bound (Equation (1b); (ii) carbon emission limit (Equation (1c)); (iii) a set of LL problems (Equations (2a)-(2c)). The tax rates determined in the UL problem will affect the power production of each unit in the LL problems. The ED problem minimizes the production cost of each demand block, subject to generation limits (Equation (2b)) and power balance condition (Equation (2c)). It decides the electricity generation of each unit, which in turn influences the quantity of total emission concerned by the government. The interaction between the government and the GENCOs constitutes a Stackelberg competition, in which the government is the leader, and the ISO is the follower. At the optimal solution, for the government, the target of emission reduction can be fulfilled with least additional tax.
For the energy industry, the production cost is minimized subject to the given tax rate. From a game theoretical point of view, the optimal solution {d * i , P * ij , ∀i,j} renders a Stackelberg equilibrium, i.e.: It should be pointed out that: 1. In the UL problem, the tax rate may vary among different units. This setting provides more flexibility and economic superiority than using an equal tax rate on all generators. To see this, adding the constraint d 1 = d 2 = … = d n to constraint (1b), the tax rate will be equal at the optimal solution. Consequently, because the feasible region becomes smaller, so the optimal value will be higher. In short, allowing the tax rate to vary according to different emission performance provides flexible and economical means to control emission effectively. 2. The value of E p is released by the government. Lower E p represents more strict restriction on carbon emission. Higher E p relaxes the emission limitations and gives priority to the economic targets. 3. Because the time scale involved in the ED problem is one year, so system upgrading is not directly considered in the current formulation. To pay the tax is the only choice of generator owners. However, if long term decisions, such as the generation expansion planning decisions, are incorporated, our modeling framework is able to model long term resorts to reduce emissions under the taxation policy in the time scale of several decades rather than paying the tax, such as upgrading the generation equipment, investment on new technologies and renewable energies can be incorporated. This will lead to a different formulation of the decision problem, but the solution method directly applies. It is a very interesting topic that we are still working on. 4. Direct Current (DC) power flow constraints also can be imposed in the LL ED problem to prevent transmission congestions. 5. As for the objective function (2a), because the prices of the electricity is constant in China, so the income of the energy industry will also be a constant provided the total demand is fixed, namely maximizing profits is equivalent to minimizing costs. In a competitive market environment, the profit, which will be the difference of income and cost, can be used as the objective function of ED. 6. It's important to restrict the cost coefficient (b i + d i e i ) in the objective (2a) different for each generator in order to avoid a degenerate ED problem in which the linear independent constraint qualification is not satisfied (the optimal solution is not unique). This can be implemented via several means. The simplest way is to add small disturbance on each (b i + d i e i ) in the ED problem, or impose the "alldifferent" constraint provided by YALMIP [22] on {b i + d i e i }, ∀i. It should be pointed out that the "alldifferent" command produces additional binary variables in the final model, and may increase the computational complexity.

Equivalent Single Level Problem
As analyzed before, the UL problem (1) and the collection of LL problems (2) are interacted with each other, they should be solved simultaneously. To use off-the-shelf solvers, it's necessary to convert the bilevel model (1)-(2) into a single level optimization problem by replacing the LL problem with its first order optimality condition. Because the LL problem appears to be a linear program (LP) when the tax rate is fixed, two options are available for this task: the KKT formulation and the primal-dual formulation [23,24]. In this paper, the former is adopted because the strong duality condition will be used to linearize the nonlinear objective function (1a) following the method presented in [21,25,26].
Note that for fixed tax rates, the LL problems (2) are LPs and decoupled with respect to block j, their corresponding KKT conditions are as follows: where the expression " a b ⊥ " means at most one of a and b can take a strictly positive value, while the other being at 0. So problem (1)-(2) becomes: However, problem (4) is instinctively hard to solve because the Mangasarian-Fromovitz constraint qualification fails at every feasible point [27] due to the presence of the complementarity and slackness conditions in Equations (4f) and (4g). Moreover, the objective function in Equation (4a) is bilinear and thus non-convex, it's difficult to retrieve a global optimal solution. In the next section, problem Equation (4) will be transformed into an equivalent MILP whose global optimum can be solved efficiently by state-of-the-art MILP solvers.

Solution of the Bilevel Model
The single level optimization problem (4) includes two kinds of nonlinearities: (1) the complementarity and slackness conditions in (4f) and (4g); (2) the bilinearity in the objective function (4a); They are subsequently linearized in the following two subsections.

Linearizing the Complementarity Constraints
The complementarity and slackness conditions in (4f) and (4g) can be linearized by introducing auxiliary binary variables l ij z and u ij z , and replacing them with following disjunctive constraints [28]: where M is a conceptual large enough constant. If M is properly chosen, it can be seen from Equations (5a)-(5b) that at most one of {P ij − P l i , l ij η } can take a strictly positive value, because l ij z is binary. It should be noticed that the value of M must be selected carefully. If M is too small, constraint (5) will impose smaller bounds to the original model [29]. In contrast, if M is too large, the resulting problem will be ill-condition and difficult to converge [24]. There is no universal method to select an appropriate value of M for an arbitrary model. This paper gives a heuristic way to determine proper bounds of dual variables based on the physical interpretation of Lagrange multipliers.
Variables u ij η can be regarded as the sensitivity of the objective function (2a) with respect to the right hand side term u i P . Suppose that u i P is changed to u u i i P P + Δ , the original optimal value f of LL problem (2) is consequently changed to f + Δf. The incremental cost Δf can be interpreted as follows: the output of unit i increases u i P Δ , according to the power balance condition (2c), there must be some other unit decreases its output by u i P Δ , so the change in the objective function Δf satisfies: As a result: The same procedure can be applied to drive the bound of l ij η . So the parameter M in Equations (5b) and (5d) can be selected as: The value of M for Equations (5a) and (5c) with primal variables can be simply selected as: To sum up, M should be selected as:

Linearizing the Objective Function
The bilinearity in the objective function (4a) can be linearized following the method presented in [21,25,26]. According to the strong duality theorem of LP, the objective functions of a primal LP and its corresponding dual LP takes the same value at the optimum. Thus the following equality holds for LL problem (2): Thus: Substituting Equation (10) into Equation (4a) yields: The right hand side of Equation (11) is linear in variables l ij η , u ij η , μ j and P ij .

Case Study
To validate the effectiveness of the proposed model and algorithm, numeric experiments on a 10-unit system and a provincial power grid in China are carried out. YALMIP is used to formulate the model. CPLEX 12.2 is used to solve related MILP problems.

10-Unit System
The parameters of generators and demand blocks of the 10-unit system are given in Tables 1 and 2, respectively. In the 10-unit system, the parameters of generators are collected from typical generators in China. Each generator has a coal consumption rate c i measured from operating data. The carbon emission coefficient e i = 44c i /12, because the molecular weight of carbon and CO 2 is 12 and 44, respectively, this means burning 1kg coal produces 44/12 kg CO 2 . As for the parameters of load blocks, the demand levels are assumed according to the capacity of generators, the time durations are the same as that in [20]. Note that the parameter b is the equivalent production cost of generators, which depends on multiple factors, such as the pure generation cost, CO 2 /SO 2 capture cost, maintenance fee, repayment for the construction investment and so on. Usually, the larger plants have lower unit production costs, but higher other costs than the smaller ones, especially the repayment during the first decade since it is built. So in a period of a few years, b i may be higher for larger generators. Nevertheless, this is not a universal law. In the next case of the Guangdong power grid, some large generators have lower b i than smaller ones.
In this case, the parameter E p is determined through two optimization problems. First, the optimal production cost problem (13) is solved. The optimal value is C min , and the corresponding carbon emission is E max : . . , , Then the optimal emission problem (14) is solved. The optimal value is E min and the corresponding production cost is C max . The results of problem (13) and (14) with respect to the considered system are shown in Table 3: Table 3. Optimal cost and emission.
Optimal production cost problem Optimal emission problem C min (CNY) E max (kg) C max (CNY) E min (kg) 1.6352 × 10 10 3.9941 × 10 10 1.8149 × 10 10 3.8775 × 10 10 The value of E p is chosen between E min and E max as follows: The optimal tax rate under different permitted emission level α is computed from MILP (12) and shown in Table 4.   Table 4 we can see, in general, generators with larger emission rate will be levied higher tax, which coincides with our intuition, but it's not necessary that high emission generators should always be taxed with high priority because their capacity may be small thus has less impact on the quantity of emission.
With the optimal tax rates fixed, the LL ED problem is solved. The production cost (∑ j ∑ i ΔT j b i P ij ), the total tax cost (∑ j ∑ i ΔT j d i e i P ij ), the permitted amount of carbon emission E p (defined in Equation (15)) and the expected emission (∑ j ∑ i ΔT j e i P ij ) are shown in Table 5. The total electricity generated by each unit is shown in Table 6. From Table 5, it can be seen that the expected emission corresponding to the optimal dispatch decision of the ISO is always less than that it is allowed. This important feature demonstrates the bilevel model endows the carbon tax scheme an indirect but effective control ability on the quantity of emission. Meanwhile, the tax cost is 1.3%-4% of the production cost when α varies from 0.2 to 0.6, which is relative reasonable, namely the profit of GENCOs will not be dramatically influenced even the selling price does not change. From Table 6, it's easily concluded that low emission generators produce more electricity with E p decreasing (α increasing).

Guangdong Power Grid of China
The realistic Guangdong power grid of China is studied in this case. One hundred and seventy four (174) generating units with a total 58,744 MW capacity are available in this system. The topology of the 500 kV main transmission network is shown in Figure 2. The parameters of demand blocks are given in Table 7.  In this system, it's found that some generators in the same power plant have identical parameters. To simplify the bilevel model (1)-(2) and reduce binary variables in MILP (12), generators with the same cost coefficient b i and carbon emission coefficient e i are aggregated into a single equivalent one whose minimal output l i P /maximal output u i P is the sum of all generators to be aggregated, while i b and e i remain the same. For example, there are three generators whose production cost coefficients b = 500 CNY/MWh, carbon emission coefficient e = 1050 kg/MWh, the capacity parameter is It's also found that the equivalent production costs of most new-built generators with CO 2 /SO 2 capture facilities are usually higher than the old ones without advanced environmental protection equipment. Therefore, a traditional cost-minimized energy scheduling scheme will result in higher emissions. If carbon tax is imposed on Guangdong power grid, bilevel model (1)-(2) is used to study how to determine an appropriate tax rate among generating units. The results are shown through Figures 3-5. The permitted amount of carbon emission E p is chosen using the same method proposed in the previous subsection, and characterized by the parameter α. In practice, the National Development and Reform Commission of China, or the Southern Power Grid Company of China is responsible to release the parameter E p for Guangdong power grid.   According to our modeling method, each aggregated unit has a corresponding optimal tax rate for a fixed permitted emission level. To illustrate them clearly, the 52 aggregated generators are divided into 6 groups. Each group contains a number of aggregated generators whose coal consumption rate is within a certain range. The coal consumption interval of each group is given in Table 8. The carbon emission coefficient e i is 44/12 times the corresponding coal consumption rate.  Figures 3 and 4 show the weighted average tax rate and total electricity production of each group under different permitted emission level, respectively. They are calculated as follows: , Electricity From Figures 3 and 4 we can see, to reduce carbon emission, a carbon tax is usually priorly levied on the generators with high coal consumption level. As a result, the electricity produced by such units is reduced. To meet the load demand, low-emission units generate more electricity, so the total emission decreases, at the cost of increasing the funds to pay.
The permitted amount of carbon emission E p and expected emission (∑ j ∑ i ΔT j e i P ij ) is compared in Table 9. It is verified that the expected emission is always no more than it is allowed. In other words, the pre-specified emission gives a valid upper bound if the production cost is minimized at the lower level. Because the optimal solution of ED happens at one vertex of its feasible polyhedron, the actual emission is not continuous. Therefore, in Table 9 (as well as in Table 5) there is slight difference between the permitted emission and expected emission. Finally, Figure 5 shows the total tax (∑ j ∑ i ΔT j d i e i P ij ) levied from Guangdong power grid under different permitted emission levels. If a "unit cost of carbon emission" is defined as follows: where d * j and P * ij is the optimal tax rate and energy offer under the permitted emission level. It's found that C U increases with the permitted emission level α increasing.

Conclusions
This paper presents a bilevel optimization model for the carbon emission taxing problem that can strike a proper balance between the target of emission reduction and the profit of emission sectors. The model is finally converted into a tractable MILP without any approximation and can be efficiently solved by commercial solvers. Our model possesses two distinguished features: (1) The Stackelberg competition between the tax maker and emission sectors under regulation is explicitly modeled; (2) The actual emission corresponding to the optimal ED solution is no more than the pre-specified value.
It's found in the case studies that: (1) according to the design parameters and operating data of typical generators in China, generators with higher emission coefficients will generally have higher tax rates; (2) the unit cost of carbon emission increases with the decreasing permitted emissions, but the total tax is still within a reasonable range when α varies from 0.2 to 0.6. These results provide practical references for relevant government agencies.
Further work will be focused on extending the proposed formulation to a generation expansion planning framework, thus upgrading the generation equipment, investment on new technologies and renewable energies could be modeled. If the uncertainty resulting from long term decisions and renewable generation is considered, the problem will be more challenging and worthy of studying.