Urban per capita carbon emission is a core indicator of the carbon emission level of a city, and CD electricity-related carbon emission contributes 65% of it (Zhang and Lu, 2021). Obviously, reductions in electricity-related carbon emissions of CDs are the key for the Chinese government to pave the path to the carbon peak and carbon neutrality goal. In 2021, the National Development and Reform Commission of China issued a notice on “Further Deepening the Market-oriented Reform of Feed-in Tariff for Coal-fired Power Generation,” requiring the power company to undertake CD power supply and maintenance tasks and charge electricity bills based on tiered electricity prices (Li et al., 2021). Obviously, the power company has become the best candidate for administrating the electricity-related carbon emission quota trading of CDs. Therefore, in this paper, the proposed electricity-related carbon emission quota trading market is administrated by a power company that is entrusted by the municipal government. In summary, the government devises the carbon emission policy. The power company is the executor of the electricity carbon emission policy entrusted by the government. City dwellers are regulated by the power company and are responsible for their carbon emissions. At the same time, CDs become participants in electricity-related carbon emission reduction.

The electricity-related carbon emission quota trading model is designed via drawing on the experiences of how household-based garbage removal fees are collected and the manner of tiered electricity pricing that increases prices after a given quota has been exceeded. In the trading model, first, the power company analyzes the carbon emission levels of individual CDs and formulates a total quota on the carbon emission over all CDs. Then, the power company purchases the calculated amount of carbon emission quota in the provincial carbon emission credit spot market and approves the basic carbon emission quotas of individual CDs according to the total urban resident population. For CDs with family size exceeding corresponding to the basic quota, their quotas will be adjusted according to their actual family sizes. The 2019 edition of the National Greenhouse Gas Inventory Guidelines issued by the Intergovernmental Panel on Climate Change (IPCC) is the basis for approving CDs’ electricity-related carbon emission quotas (Viola, 2020; McLean and Gibbs, 2022).

In fact, carbon emissions are not generated by the electricity consumption of CDs but from the thermal units that burn coal for generating the consumed electricity. In addition to thermal units, supply resources of CDs typically include green energy resources and energy storage systems. Among them, the carbon emission intensities of thermal units are determined according to unit types; the carbon emission intensities of green energy resources, such as hydropower, solar energy, and wind energy, are set as 0. Energy storage systems do not directly generate carbon emissions, but their stored electricity could partly come from thermal units, therefore implicitly introducing carbon emissions.

We use n to denote the number of resources. All resources are considered current resources. After the network power flow snapshot is selected, all the ground branches at the CD connection node b can be considered equivalent impedances connected to node b , denoted as r . The current injected to node b can be represented as the sum of the currents from all resources in the network, and the power injected to node b can be calculated as in (1), where I j A l l is the current from resource j . In fact, P b A l l consists of three parts: thermal power, green power, and power from energy storage systems that are correspondingly denoted as P b T , P b G , and P b S . As an example, P b D represents the power drawn by node b from resource c . P b A l l = ∑ j = 1 n r I j A l l 2 , (1) P c T = P c D P b T + P b G + P b S P b D . (2)

In Eq. 1, the line losses of the distribution area caused by supplying CDs are not considered. However, the municipal government may require including the line losses in the calculation of carbon emissions in order to accurately capture the actual amounts. We use t to denote the total number of time intervals in a statistical period and τ j to denote the electricity-related carbon emission of different resources. The carbon emission of green energy resources is 0, and the carbon emission of the energy storage is determined by its stored energy from thermal units. P i A t represents the total power injection to node b ; P i L o s s is the line loss at time interval j . The total carbon emissions of node b in a statistical period can be calculated as in Eq. 3. C c T = ∑ j = 1 , i = 1 t , n P i A T τ j × 1 + P i L o s s . (3)

Using Eq. 3, the total carbon emissions of all CDs in the city can be obtained. In order to avoid the increase of alternative energy usage that leads to extra carbon emissions due to enforced reduction in the proportion of CDs’ electricity-related carbon emissions, a proportion factor is adopted to approve the basic quota of CDs’ electricity-related carbon emissions. The proportion factor is usually issued by carbon emission regulatory agencies based on their statistics and monitoring data. In our case, the proportion factor was issued by the SACE.

According to the forecasted total electricity-related carbon emissions of CDs, the power company purchases an equivalent amount of carbon emissions quota C c A l l from the carbon emission credit spot market. We denote the number of CDs in the city as m . When the proportion of all CDs’ electricity-related carbon emission amount fluctuates beyond the threshold by ∆ φ , the quota will be deducted by ∆ σ . The fluctuation threshold is set referring to the typical value issued by the SACE.

We use C c P and C c C to, respectively, denote the electricity-related carbon emission of CD b in the previous and current statistical periods and C c A P and C c A C to, respectively, denote the total electricity electricity-related carbon emission of the city in the previous and current statistical periods. On this basis, a CD’s basic electricity-related carbon emission quota can be calculated as in Eq. 4, where C c P C c A P and C c C C c A C satisfy Eq. 5. (5) represents the simple fact that the share of the electricity-related carbon emission of a CD over the total electricity-related carbon emission of the city must be less than 1 in both the previous and current statistical periods. Based on Eq. 4, the power company charges CDs the basic carbon emission fees from their deposits. C c F = C c A l l m × 1 − C c P C c A P − C c C C c A C × Δ σ , (4) 0 ≤ C c P C c A P ≤ 1 0 ≤ C c C C c A C ≤ 1 . (5)

The purpose of establishing the electricity-related carbon emission quota trading market is to help CDs recognize the limitations on carbon emissions and thereby motivate them to actively reduce their carbon emissions. In the electricity-related carbon emission quota trading market, CDs with different incomes would show different degrees of sensitivity to the carbon emission fee. Intuitively, high-income CDs would be less sensitive to it and might be willing to pay extra beyond the bill of the basic quota; by contrast, low-income CDs would be highly sensitive to the carbon emission fee and might be willing to sell their remaining quotas for offsetting part of their electricity bills. To this end, in this paper, a clustering analysis based on CDs’ electricity usage is conducted to obtain their sensitivity level classification.

Density-based spatial clustering of applications with noise (DBSCAN) is a clustering method based on density detection and with noise canceling ability. This method has the ability to form clusters from high-density regions with limited features (Juan et al., 2021; He et al., 2022) and does not require pre-specifying the number of clusters. However, in the case of large CD samples, its convergence efficiency would be of concern. In other words, when facing a large number of CD samples, the convergence time of the DBSCAN algorithm may become unacceptable and divergence may happen. To this end, in this paper, the DBSCAN algorithm is improved through limiting the size of searching neighbors, achieving the processing capability of a large number of CD samples.

Considering the Mahalanobis distance involves covariance of distances and can comprehensively consider the characteristic relationship between CD samples (Chen et al., 2019; Cui and Xia, 2022), in the DBSCAN algorithm, it is used to measure the distances between CD samples.

Using e i and e k to represent the feature vectors of CD samples i and k , respectively, the Mahalanobis distance of the two CD samples can be calculated as in Eq. 6, where n c is the number of features and O a is the covariance matrix of the Mahalanobis distances. d i k = e i − e k n c O − 1 e i − e k . (6)

A neighborhood parameter, ϖ , is used by the DBSCAN algorithm to depict the tightness of the CD cluster distribution. For e i , taking it as the center, the circle with radius of ϖ forms its neighborhood V ϖ , as represented in Eq. 7, where G a is the total number of samples. V ϖ = e i ∈ G a d i k ≤ ϖ . (7)

In the DBSCAN algorithm, if the neighborhood of e i contains multiple CD samples, it becomes the core object, and these CD samples in its neighborhood are directly density-reachable to e i . For any e i and e j , if they are not directly density-reachable, but if there exists a path of samples that starts from e i and ends at e j , while any two consecutive DC samples are directly density-reachable, e i and e j are density-reachable and form a directly density-reachable closure. As an example, considering e i , e k , e k + 1 , and e j , if the pairs of e i and e k , e k and e k + 1 , and e k + 1 and e j are all directly density-reachable, e i and e j are density-reachable. This is further illustrated as in Figure 2. For e i and e k , if there exists e s , that is, density-reachable to e i and e k , e i and e k are density-connected, as illustrated in Figure 3. Thereafter, based on the regional scale of the urban power supply grid, the searching threshold is set as Δ λ .

From Figures 2, 3, it can be seen that the DBSCAN algorithm explores density-connected CDs from density-reachable CDs and clusters those CDs into one cluster, finally forming clusters of CDs. This process is repeated until all samples are traversed.

The proposed credit system is essentially a user incentive model. By setting different preferential quota prices for CDs with different credit ratings, CDs could be motivated to participate in carbon emission reduction. The power company establishes and operates the credit system for CDs, including updating their credit scores and credit ratings and adjusting the floating carbon emission quota prices.

The designed credit system adopts a hundred-mark system. The range of 60–100 points indicates qualified credit. The 60th point corresponds to the base price of electricity-related carbon emission quota. With credit scores between 61 and 100 points for CDs who purchase and sell quotas, the prices will be respectively lowered and increased to different degrees on the basis of the base price. The range of 0–59 points indicates unqualified credit. For CDs in this area, on the contrary, on the basis of the base price, their purchase and sell prices will be respectively lowered and increased to different degrees. Table 1 shows the carbon emission credit ratings and the corresponding floating prices.

A CD’s initial credit score is set by the power company based on the standards issued by the electricity–carbon industry association. For simplicity, in this paper, the credit score initialization is based on the punctuality of the CDs’ electricity bill payment. Considering the due time of an electricity bill is T e i and the actual payment time is T h i , if T h i is earlier than T e i , a normal payment will be recorded; otherwise, it will be recorded as a default.

In a statistical period Z a , the number of normal payments is denoted as f a and the total number of bill cycles is denoted as F a . Then, the normal payment ratio is equal to f a F a . If there is no default, namely, delayed payments, this ratio is equal to 1. The initial credit score of a CD is determined by the normal payment ratio and the delayed time in the presence of defaults. Assuming the basic credit score is G b , the initial credit score of a CD, namely, G c , can be calculated as in Eq. 8, where G c ∈ 0 , 100 . G c = G b + f a F a − ∑ i = 1 n T T i R − T i E T i E . (8)

A CD’s electricity-related carbon emission quota trading records, as well as violations on delivering traded quotas, will be stored, and at the same time, its credit scores will be dynamically updated based on these records. The power company supervises the carbon emission trading between CDs and rewards or penalizes CDs based on the rules in Table 1. Before the final delivery of traded quotas, all deals are first checked by the power company and then reviewed by the municipal government to verify the fairness and protect the interests of all involved parties.

On the basis of the initial credit score, the power company enforces the rules of credit score updating, as designed in Table 2. The rules are designed referring the standards of the electricity–carbon industry association.

The purpose of analyzing the total selling and purchasing offers is to judge whether the total carbon emission quota of CDs has been reduced based on the trading information. If the total amount of purchasing offers of CDs is greater than that of the selling offers, which means the total amount of carbon emissions will increase, the government will apply the carbon emission protection mechanism. Carbon emission quota trading will be settled with the government-approved ceiling price. For the part that exceeds the supply, the power company will purchase more green electricity on behalf of CDs so that CDs can offset their excess. The corresponding price will increase on the basis of the government-approved ceiling price. If the total amount of purchasing offers of CDs is less than that of the selling offers, a Stackelberg game between CD buyers, CD sellers, and the power company will be carried out to determine the price.

The total amount of CDs’ electricity-related carbon emissions, C c C A , can be calculated as in Eq. 9, where C c F represents the basic quota of a CD; C c S A represents the total sold quota; and C c B A represents the total purchased quota. C c C A = C c F × m + C c S A − C c B A . (9)

If C c C A increases, the protection price of the carbon emission quota price will be implemented. For a CD who purchases carbon emission quota, its electricity carbon emission expenditure can be calculated as in Eq. 1, where V F represents the base price; C c R e a l i t y represents the actual carbon emission quota; V Max represents the maximum price; C c B u y represents quota purchased from other CDs; and Z E is the power company’s green power purchase price. Z E = C c F × V F + C c B u y × V Max + C c R e a l i t y − C c F − C c B u y × V F + V Max . (10)

For a CD who sells carbon emission quota, the income of the CD, Z S , can be calculated as in Eq. 11, where C l represents quota sold to other CDs. It is worthwhile to mention that a fee taking s of the income will be charged by the power company. Z S = C c B u y × V Max × 1 − s . (11)

The carbon emission pricing analysis is used to determine the base price. On the premise that the price in carbon emission quota trading is higher than that in the carbon emission credit spot market, CD sellers, CD buyers, and the power company will play a Stackelberg game, seeking an equilibrium price for the three parties. This equilibrium price will be set as the settlement price.

The Stackelberg game is a leader–follower game model. In this model, the leader with advantages first gives its competitive strategy, and then the followers give their own strategies responding to the strategy of the leader (Mukaidani et al., 2021; Bouchaib et al., 2023). In the proposed carbon emission quota trading market, the competitive positions of the power company and CDs are obviously not equal. Intuitively, the power company will be the leader in the game and CDs, either as buyers or sellers, are the followers in the game.

In the game, we use q to denote the power company, u c to denote CD sellers, and u d to denote CD buyers. They are collected by O as in Eq. 12. O = u c , u d , q . (12)

In the game, considering the power company as the leader bids θ prices, denoted by l , the leader’s bidding price vector can be written as Eq. 13. L = l 1 , l 2 , … l θ . (13)

Considering the CD buyers, as the follower, bid ε quantities, denoted by y , the buying quantity vector can be written as Eq. 14. Y = y 1 , y 2 , … y ε . (14)

Considering the CD sellers, as the follower, bid ψ quantities, denoted by x , the selling quantity vector can be written as Eq. 15. X = x 1 , x 2 , … x ψ . (15)

In the Stackelberg game, the expenditure of CD buyers is denoted as Z i B u y ; the income of CD sellers is denoted as Z i S ; and the profit of the power company is denoted as Z i S T . The Stackelberg game model can be expressed as in Eq. 16. W = O ; L ; Y , X ; Z i B u y , Z i S , Z i S T . (16)

In the gaming process represented by Eq. 16, the CD sellers pursue higher incomes and credit scores; and the CD buyers pursue lower expenditures and higher credit scores. The power company pursues the highest trading agency fee under the condition that the overall electricity-related carbon emissions can be reduced. All the three parties use benefit-maximizing strategies in the game until they reach an equilibrium.

In the game, in a trading cycle k , the quota selling price is V i S C ; the sold quota is C i S C ; the cost for the basic quota is V i Q ; and the credit score loss cost is denoted as C i C . The income of a CD seller can be calculated as in Eq. 17. C i S R = ∑ i = 1 k V i S C × C i S C − V i Q − C i C . (17)

The CD seller needs to satisfy Eq. 18, where C i S M a x is the basic quota. C i S M i n ≤ C i S R ≤ C i S M a x . (18)

In a trading cycle, the quota buying price is V j B u y ; the purchased quota is C j B u y ; the cost for the basic quota is C j B a s e ; and the credit score loss cost is denoted as C j C . The expenditure of a CD buyer can be calculated as in Eq. 19. C j E = ∑ j = 1 k V j B u y × C j B u y + C j B a s e − C j C . (19)

CDs are connected to distribution transformers through 0.4 kV distribution lines. The upper bound of their carbon emissions is indeed limited by the rated capacities of the distribution transformers that they are connected to. Thus, the power input to CD buyers needs to satisfy constraint Eq. 20. The upper bound is determined by the rated capacity of the distribution transformer, C j R a t e d , and the total load of CDs that connect to the transformer, C j S u m . Distribution transformer load capacity (kW) = the rated capacity (kVA) × the power factor. The power factor is set as 0.8 according to the industry standard of China. The corresponding description has been revised in the paper. C i S M i n = C j B M i n C i S M a x ≤ C j B M a x C j B M a x = C j R a t e d − C j S u m . (20)

The power supply company only administrates the trading market and charges agency fees according to the traded amount of carbon emission quotas. The benefit of the power company is, in fact, the collected agency fees as in Eq. 11. The benefit of the power company can be calculated as in Eq. 21, where n j is the number of CD sellers and n k is the number of CD buyers. S c i and S q j are the income of CD seller i and the expenditure of CD buyer j , respectively. C j R = s ∑ j = 1 e C i E . (21)

We use β i , β j , and β k strategy spaces of CD sellers, CD buyers, and the power company, respectively, and then the game between them can be represented as Eq. 22. max S c s . t . c m 1 , c m 2 , … , c m t b ∈ β i max S q s . t . c n 1 , c n 2 , … , c n t b ∈ β j max S d s . t . v h 1 , v h 1 , … v h θ ∈ β k . (22)

With a Nash equilibrium solution, neither the power supply company nor any CD buyers and sellers can unilaterally change the strategy to achieve more profit. To achieve Nash equilibrium in the game, the optimal conditions of all participants should be satisfied. Under this premise, CD sellers and the power company look for the best quota buying price from high to low for profit maximization, while CD buyers look for the best quota selling price from low to high for cost minimization. In the process of solving the Nash equilibrium of the Stackelberg game, CD sellers, the power company, and CD buyers offer the initial prices. Then, in each iteration, the price change gradients of the three participants are calculated and their offered prices are updated. The iterative process continues until their current price offers are equal to the maximum values of their strict convex functions so as to obtain the Nash equilibrium. The detailed solving process is available in Jin et al. (2021) and Veeramsetty (2021) and thus is not described redundantly in this paper.

In the settlement, the credit scores can be used to exchange for preferential treatment. When the exchange amount is greater than ∆ u , the credit score is deducted by ∆ g . Thereafter, the credit score will be updated following the observation in Table 2. The credit score is updated in every trading cycle. In fact, in the settlement, the credit value deduction mechanism for CDs that have enjoyed preferential treatment is to avoid CDs becoming inert after reaching high credit scores and credit ratings, thereby repeatedly incentivizing CDs to participate in carbon emission quota trading.

For a CD seller, if his/her actual delivered quota changes due to the increase or decrease in his/her own power consumption, there will be a deviation from its offered amount. The ratio of deviation can be calculated as in Eq. 23, where F S is the ratio of deviation; C j R S is the offered quota amount; and C j R is the actual sold quota amount. When F S is negtive, i.e., the actual sold is less than offered, default happens. F S = C j R − C j R S C j R S , (23) F B = C j B − C j R B C j R B . (24)

For a CD buyer, if his/her actual bought quota decreases due to the increase or decrease in his/her own power consumption, there will be a deviation from its asked amount. The ratio of deviation can be calculated as in Eq. 24, where F B is the ratio of deviation; C j R B is the actual bought quota amount; and C j B is the asked quota amount. When F B is positive, i.e., the actual bought amount is less than asked, default happens.

Using Q a and Q b to denote the initial credit scores of CD buyers and CD sellers, and v c 1 and v c 2 to denote carbon emission quota pruchasing and selling prices, respectively, credit scores, Q x and Q y , after updating can be calculated as Eqs 25, 26. Q x = Q a + − F S Δ k 1 Δ f 1 + C j R Δ k 2 Δ f 2 − C j R × V j R E 1 Δ μ Δ g Δ f 3 , (25) Q y = Q b + − F B Δ h 1 Δ j 1 + C j B Δ h 2 Δ j 2 − C j B × V j R E 2 Δ μ Δ g Δ j 3 . (26)

An actual community in a city in China is selected to validate the proposed electricity-related carbon emission quota trading market model. This community is equipped with 24 distribution transformers and has 10,762 CDs. The resident population of this community is 22,736. A single CD has only one house. The carbon emission quota trading is carried out monthly from January to December 2021. We consider that each CD corresponds to one unit (apartment) in this community. The typical values from the SACE are shown in Tables 1, 2. Other detailed parameters are shown in Tables 3, 4.

The Hungarian algorithm in Du et al. (2020) is compared with the proposed trading mechanism. The Hungarian algorithm has the characteristics of accurate results and high analysis efficiency and is widely used in carbon emission credit trading markets.

All numerical simulations are conducted on a server with Intel Xeon E5-2650 CPU (8 cores), 32 GB RAM, and 10 TB SSD. The operating system is Windows Server 2019.

The electricity-related carbon emission quota trading market proposed in this paper is cleared once a month. Considering China’s carbon emission credit price will gradually be in line with that of Europe and the United States, the initial carbon emission quota price is set at 0.26¥/kg, referring to the price in Europe (Frank, 2021).

In this paper, the tiered electricity prices of CDs have three tiers based on the electricity consumption: 1) within 100 kWh, charged at benchmark electricity price, 2) 101–200 kWh, charged at 120% of the benchmark electricity price, and 3) above 200 kWh, charged at 150% of the benchmark electricity price. Only CDs with electricity consumption in the first tier sell their remaining carbon emission quotas, while CDs in the other two tiers need to purchase quotas. For credit score and rating initialization, the electricity bill payment records of 60 months, from 2015 to 2020, are used. The power and capacity of equipped energy storage systems are 200 and 500 KWh, respectively. The capacity of distribution transformers supplying this community is 2000 KVA. The maximum purchase of carbon emission quota will be restricted by the capacity of the distribution transformers.

The electricity consumption of this community in January 2021 was 1872.6 MWh, and the line loss rate of its distribution area was 4.28%. The daily carbon emissions of this community in January 2021 are shown in Figure 4. In the basic carbon emission quota approval, the proportion of its electricity-related carbon emission in its total carbon emission amount increase does not exceed the threshold of Δ φ . As shown in Figure 4, the carbon emission quota for this community is 164.27 kg.