(Article reproduction) A preliminary study on the calculation method of carbon emission flow in power system (including matlab code)

References: Zhou Tianrui, Kang Chongqing, Xu Qianyao, Chen Qixin. Preliminary Study on the Calculation Method of Carbon Emission Flow in Power System [J]. Automation of Electric Power System, 2012,36(11):44-49.

The dual carbon target is still very popular in the power system field recently It’s hot, and the calculation of carbon emission flow is also a hot topic. This paper can be said to be the most basic model, and the number of citations is also very high. I will write a blog to simply reproduce this paper. Friends need it can refer to.

1. Fundamentals

Taking the IEEE14-node system as an example, the basic principles and steps of power system carbon emission flow analysis are explained below. Suppose there are N nodes in the system, K nodes have generator power injection, and M nodes have loads.

1.1 Power flow calculation

The carbon emission flow and the power flow of the power system are closely related but different. Overall, the flow calculation is the basis for the calculation of the carbon emission flow. Therefore, before calculating the carbon emission flow, the power flow calculation of the system should be carried out first. There are many methods of power flow calculation, such as the most commonly used cow-pull method, PQ decomposition method, or calling matpower in matlab to get the results directly. I won’t introduce too much here. For specific steps, please refer to Newton’s method and P-Q decomposition method IEEE14 system Power flow calculation (with matlab code)

After obtaining the results of power flow calculation, data such as branch power flow, generator injected power, node injected power, and load size need to be used for the calculation of carbon emission flow.

1.2 Introduction to Correlation Matrix

(1) Branch power flow distribution matrix

The branch power flow distribution matrix is an N-order square matrix, represented by $P_B$ . The purpose of defining this matrix is to describe the active power flow distribution of the power system, and the calculation formula is as follows:

$P_B_{ij}=\left\{\begin{matrix}p_{ij},p_{ij}>0\ 0,p_{ij}\leq0\end{matrix}\right.” class =”mathcode” src=”//i2.wp.com/latex.csdn.net/eq?P_B_{ij}=\left\{\begin{matrix}p_{ij},p_{ij}>0\ 0,p_{ij}\leq0\end{matrix}\right.”></p> <p>In the formula, <img alt=$ indicates the flow through branch ij active power.

%% branch flow distribution matrix PB
for k=1:length(branch(:,1))
    branch(k,14)=PB(branch(k,1),branch(k,2));
    branch(k,15)=PB(branch(k,2),branch(k,1));
end
for k=1:14
    for kk=1:14
        if PB(k,kk)<0
           PB(k,kk)=0;
        end
    end
end

(2) Unit injection distribution matrix

The unit injection distribution matrix is a matrix of K X N order, represented by $P_G$ . The purpose of defining this matrix is to describe the active power injected into the system by all generator sets, and the calculation formula is as follows:

$P_G_{kj}=\left\{\begin{matrix} P_k,k,j\in G_{kj}\ 0,k,j\ otin G_{kj} \end{matrix}\right .$

In the formula, $G_{kj}$ represents a collection of generators. The meaning of this formula is that if the kth generator set is connected to node j, then $P_{Gkj}$ Equal to generator injected power, otherwise zero.

%% unit injection distribution matrix
PG=zeros(K,N);
for k=1:length(gen(:,1))
    PG(k,gen(k,1))=gen(k,2);
end

(3) Load distribution matrix

The load distribution matrix is an order MXN matrix, represented by $P_L$ . The purpose of defining this matrix is to describe the connection relationship between all electric loads and the power system and the amount of active load. The calculation formula is as follows:

$P_{Lmj}=\left\{\begin{matrix} pload_{mj},m,j \in M_{mj} \ 0,else \end{matrix}\right.$

In the formula, $M_{mj}$ represents the load set. The meaning of this formula is that if the mth load is connected to node j, then $P_{Lmj}$ Equal to the magnitude of the active load, otherwise zero.

PL=zeros(M,N);
bus_copy=bus;
for k=1:M
    for kk=1:N
        if bus_copy(kk,3)>0
            PL(k,kk)=bus_copy(kk,3);
            bus_copy(kk,3)=0;
            break
        end
    end
end

(4) Node flux matrix

The node active flux matrix is an N-order diagonal matrix, represented by $P_N$ , calculated The formula is:

$P_Z=[P_B\ P_G]^T \ P_N=diag(\xi _{N + K}P_Z)$

In the formula, $P_Z$ is the auxiliary variable introduced, $\xi _{N + K}$ is an element all A row vector of 1.

%% Auxiliary variable PZ
xigama=ones(1,K+N);
PZ=[PB;PG];
%% node active flux matrix
PN=diag(xigama*PZ);

1.3 Calculation of carbon emission flow

(1) Carbon emission intensity vector of generating units

The carbon emission intensity vector of the generator set is a K-dimensional column vector, represented by $E_G$ , are generally regarded as known conditions.

(2) Node carbon potential vector

The node carbon potential vector is an N-dimensional column vector, represented by $E_N$ , calculated The formula is as follows:

$E_N=(P_N-P_B^T)^{-1}P_G^TE_G$

%% node carbon potential distribution vector
EN=(PN-PB')^(-1)*PG'*EG;

(2) Branch carbon flow rate distribution matrix

After the node carbon potential vector is calculated, similar to the power flow calculation, the carbon flow rate of each branch of the system can be further obtained. Thus, the branch carbon flow rate distribution matrix is defined as an N-th order square matrix, using $R_B$ indicates that the calculation formula is:

$R_B=P_Bdiag(E_N)$

%% branch carbon sulfur ratio distribution matrix
RB=PB*diag(EN);

(3) Loaded carbon flow rate vector

After the node carbon potential vector is calculated, the carbon emission intensity of electricity consumption at the node load is equal to the node carbon potential. Combined with the load distribution matrix, the carbon flow rate corresponding to all loads can be obtained. The physical meaning is the carbon emissions generated by the power generation side for the supply node load per unit time. Then the load carbon flow rate vector is represented by $R_L$ , and the calculation formula is:

$R_L=P_LE_N$

%% load carbon flow rate vector
RL=PL*EN;

2. Comparison of running results and original literature

Since the branch loss is not considered in the original literature, the DC power flow method is used to calculate the power flow distribution of the system. In the code, I use the cow-pulling method in polar coordinates, so the results will have some errors, but they are basically the same.

2.1 Node active flux and node carbon potential

(1) Original results

(2) Calculation results

2.2 Loaded Carbon Flow Rate

(1) Original results

(2) Calculation results

PS: You can see the complete code by poking this link:

Matlab code reproduction of “Preliminary Study on the Calculation Method of Carbon Emission Flow in Power System”