The key step of Unscented Kalman Filter is to generate Sigma points through unscented transformation based on the probability distribution of the state variable. The specific generation method is shown in Equations 1, 2: χ k − 1 , j = x k − 1 , j = 0 x k − 1 + n + λ P k j , j = 1 , . . . , n x k − 1 − n + λ P k j , j = n + 1 , . . . , 2 n (1) λ = α 2 n + κ − n (2) where χ k − 1 , j represents the Sigma point sets, and n is the dimension of the state variable. α is the scale correction factor, and Κ is a free parameter to ensure that the matrix n + λ P k is semi-definite. n + λ P k i is column i of the matrix.

Based on the Unscented Kalman Filter algorithm, the state estimation steps of the distribution system are as follows:

1) The Sigma point set { χ k − 1 } is generated through unscented transformation.

2) State prediction: In terms of prediction, machine learning and exponential smoothing method are adopted in abundant research. Although machine learning is better performed in many aspects, its training relies on reliable data sets and is time-consuming. Therefore, in this paper, the prediction of the next moment is carried out separately for each Sigma point using exponential smoothing method based on the previous state of these points, and the mean and covariance of the points are calculated in Equations 3–7: ξ k | k − 1 , j = f χ k − 1 , j (3) x ¯ k | k − 1 = ∑ j = 0 2 n W j , m ξ k | k − 1 , j (4) P k | k − 1 = ∑ j = 0 2 n W j , c ξ k | k − 1 , j − x ¯ k | k − 1 ξ k | k − 1 , j − x ¯ k | k − 1 T + Q k (5) W j , m = λ n + λ , j = 0 1 2 n + λ , j = 1 , . . . , 2 n (6) W j , c = λ n + λ + 1 − α 2 + β , j = 0 1 2 n + λ , j = 1 , . . . , 2 n (7) where Q k is the process noise covariance matrix, W j , m is the weight of the mean of the Sigma point set, and W j , c is the weight of the covariance of the Sigma point set.

3) Measurement prediction: Based on the measurement function and the state variable, the prediction of measurement point set z k | k − 1 , j and its mean z ¯ k | k − 1 are obtained in Equations 8, 9. The covariance matrix P k , z z of the measurement prediction, the covariance matrix P k , x z of the state variable prediction and the measurement prediction and the Kalman gain K k are calculated respectively in Equations 10–12. z k | k − 1 , j = h ξ k | k − 1 , j (8) z ¯ k | k − 1 = ∑ j = 0 2 n W j , m z k | k − 1 , j (9) P k , z z = ∑ j = 0 2 n W j , c z k | k − 1 , j − z ¯ k | k − 1 z k | k − 1 , j − z ¯ k | k − 1 T + R k (10) P k , x z = ∑ j = 0 2 n W j , c ξ k | k − 1 , j − x ¯ k | k − 1 z k | k − 1 , j − z ¯ k | k − 1 T (11) K k = P k , x z P k , z z − 1 (12) where R k is the covariance matrix of the measurement noise.

4) State update: After obtaining a new round of measurement z k , the state estimation x k   and its covariance matrix P k are updated according to Equations 13, 14. x k = x ¯ k | k − 1 + K k z k − z ¯ k | k − 1 (13) P k = P k | k − 1 − K k P k , z z K k T (14)

In traditional risk assessment, a sampling process based on the probability distribution of uncertainty factors is often required before conducting risk assessment with the scenarios it obtained. However, in real-time operation, it is difficult to adjust the probability distribution of uncertainty factors simultaneously according to the real measurement information, which means that the traditional risk assessment methods are lack of timeliness and not suitable for real-time dynamic risk assessment of the system.

There are two core steps in the process of state estimation of the distribution system based on Unscented Kalman Filter. The first is to predict the state variable through Sigma points, and the second is to revise the prediction of state variable according to new measurement data. The sampling that is completed in the prediction step of the Unscented Kalman Filter algorithm can not only be combined with the measurement information of the next moment to estimate the operational state of the distribution system, but also be used for the early evaluation of the operational risk before the measurement information of the next moment is obtained. Therefore, this paper takes the output of the prediction step of the Unscented Kalman Filter as the input of the risk assessment model, and the complete process is shown in Figure 1. Based on the various Sigma points generated by the Unscented Kalman Filter, the system’s operational state at the next moment, namely the voltage amplitude and phase of all buses, can be forecasted, and the power flow can be calculated immediately based on the bus voltage. Then, the operational risk of the distribution system is fully evaluated with the obtained operational state and the established risk indices. By making full use of the state estimation results, the difficulty of simultaneously updating the probability distribution is avoided and the calculation efficiency is significantly improved, so as to enable real-time operational risk assessment and provide the corresponding data basis for the risk control, which contributes to better operation of the distribution system.

The overvoltage risk index reflects the risk when the voltage of each bus of the distribution system exceeds the safe voltage range, as is calculated in Equation 15: R i v h = ∑ s ∈ D τ s B i , s v h V i , s S e v v h V i , s (15) where R i v h is the overvoltage risk of bus i, and τ s is the probability of the scene. B i , s v h is a flag which indicates whether there is risk on bus i. V i , s is the voltage of bus i in scene s, and S e v v h is a function describing the severity of the risk, as is shown in Equation 16: S e v v h V i , s = e V i , s − V max − 1 e − 1 (16) where V max = 1.05 p . u .

The low voltage risk index reflects the risk when the bus voltage falls below the lower voltage limit and is calculated in Equation 17: R i v l = ∑ s ∈ D τ s B i , s v l V i , s S e v v l V i , s (17) where R i v l is the low voltage risk of bus i, and B i , s v l is a flag which indicates whether there is risk on bus i. S e v v l is a function describing the severity of the risk, as is shown in Equation 18: S e v v l V i , s = e V min − V i , s − 1 e − 1 (18) where V min = 0.95 p . u .

The line overload risk index reflects the risk when the active power flow of the line exceeds the maximum limit, as is calculated in Equation 19: R j p = ∑ s ∈ D τ s B j , s P j , s S e v P P j , s (19) where R j p is the line overload risk of branch j; B j , s marks whether branch j is overloaded in scene s; P j , s is the active power flow of line j in scene s; S e v P is a function that describes the severity of the line overload risk, as in Equation 20: S e v P P j , s = P j , s − P j , ⁡ max (20) where P j , ⁡ max is the rated capacity of branch j.

Unlike line overload risk which results from high load demand at buses, reverse line overload is caused by high renewable energy penetration in the system. When there is a large amount of renewable power in the system, the residual power will be transmitted backwards through the line from the installation point of the renewable power to the head of the line. The reverse line overload risk index reflects the risk that the reverse active power carried on the line exceeds the maximum limit, (the risk of bus voltage increase caused by incomplete absorption of the renewable power has been included in the overvoltage operational risk mentioned previously, so only the risk caused by its impact on the power flow is considered here), and it is calculated in Equation 21: R j pr = ∑ s ∈ D τ s B j , s r P j , s S e v Pr P j , s (21) where R j pr is the risk of reverse line overload at branch j; B j , s r marks whether branch j is reversely overloaded in scene s; S e v Pr is a function describing the severity of the risk of reverse line overload, as is calculated in Equation 22: S e v Pr P j , s = − P j , s − P j , ⁡ max (22)

Because of the different location of each bus in the distribution system and their different amount of load demand, the voltage risk of each bus has varying degrees of impact, so this paper assigns different weights to each bus, and summarizes voltage risks of buses to the system level using the bus weight, as is shown in Equations 23, 24. R s y s v h = ∑ i = 1 N b w i R i v h (23) R s y s v l = ∑ i = 1 N b w i R i v l (24) where R s y s v h and R s y s v l are the risks of overvoltage and low voltage in the distribution system, respectively; N b is the total number of buses in the distribution system; w i is the weight coefficient of each bus, and the calculation is shown in Equation 25: w i = L i E P i ∑ E P i (25) where E P i is the active power of bus i; L i is the number of branches connected to the bus.

The overload and reversing overload operational risks of distribution lines are the cumulative sum of the overload operational risks of each branch in the distribution system, as is shown in Equations 26, 27: R s y s p = ∑ j = 1 N l R j , i p (26) R s y s pr = ∑ j = 1 N l R j , i pr (27) where R s y s p and R s y s pr are the risks of overload and reverse line overload in the distribution system, respectively; N l is the total number of branches in the distribution system.

When there is an excessive amount renewable energy that cannot be fully absorbed, it can lead to serious overvoltage and reverse line overload risks in the distribution system. In this case, the impact of renewable energy can be diminished by curtailing power, but it can also result in lower utilization rates of renewable energy. By using the prediction of the system operational state, namely the bus voltage, during the state estimation process, the power flow on each branch can be quickly calculated, and the injection power of each bus can be obtained. Based on this, this section constructs an optimization model to optimize the amount of renewable energy curtailment required to avoid serious voltage and power flow risks before calculating the risk of renewable energy curtailment. The objective function of the optimization model is described in Equation 28: f C r , U o , P o = C r + M U o + P o (28) where C r represents the amount of curtailed renewable energy power, M is a relatively large adjustment coefficient that ensures that the voltage and current can be controlled within the specified limits through renewable energy curtailment, and U 0 and P 0 represent the degree to which the voltage and the power flow exceed the specified limits. If M is too large, it may lead to computational difficulties. However, if M is too small, it may result in insufficient renewable energy curtailment. For a specific distribution system, the value of M can be figured out through testing to confirm its appropriateness. Specifically, if the optimization results show that the voltage and power flow are not controlled within the limits and there is room for further curtailment of renewable energy, it indicates that the value of M needs to be increased. The triggering point of curtailment is determined based on the notion that it only occurs when the risks reach a certain level, and the appropriate points can be set based on varying conditions in different models. For bus voltage beyond the limit, this optimization model sets the action of renewable energy curtailment when the bus voltage amplitude reaches below 0.9 times or above 1.1 times the standard value per unit. For line power flow exceeding the limit, the renewable energy curtailment is triggered when the power flow or reverse power flow exceeds 1.2 times the standard value per unit. The expressions for U 0 and P 0 are described in Equation 29. U o = V i , s − 1.1 p . u . , 0 , 0.9 p . u . − V i , s , i f   V i , s ≥ 1.1 p . u . i f   0.9 p . u . ≤ V i , s ≤ 1.1 p . u . i f   V i , s ≤ 0.9 p . u . , P o = P j , s − 1.2 p . u . , 0 , i f i f   P j , s ≥ 1.2 p . u . P j , s < 1.2 p . u . (29)

In addition, the optimization model of renewable energy curtailment is also constrained by power flow and power balance. As optimization problems based on accurate power flow constraints cannot be solved, this paper adopts the following linear approximation power flow constraints (Yang et al., 2018; Yang et al., 2019). P i j = v i − v j g i j − b i j θ i j + Δ P i j (30) Q i j = − 1 + 2 v i b i j , 0 − v i − v j b i j − g i j θ i j + Δ Q i j (31)

In the Equations 30, 31, P i j and Q i j represent the active and reactive power flowing from bus i to bus j respectively; Δ P i j and Δ Q i j represent active and reactive power losses respectively; g i j and b i j represent the conductance and susceptance of the branch from bus i to bus j respectively; θ i j represents the voltage phase difference between bus i and bus j.

Given the complicated operational requirements of distribution systems and the requirement for timeliness of its real-time risk assessment, the simplification of calculation needs to be considered. Although approximating network losses may sacrifice a bit of accuracy, it is an effective way to decouple voltage amplitude v and phase angle θ , as well as greatly reduce computational burden. In order to obtain approximate expressions for the network losses Δ P i j and Δ Q i j , this paper uses a quadratic Taylor expansion of the cosine function of the phase angle difference, θ i j ≈ 1 − θ i j 2 / 2 , and ignores higher-order terms given that the voltage phase angle difference is small. In addition, considering that the bus voltage amplitude is close to 1p.u., an approximate value of V i V j θ i j 2 ≈ θ i j 2 is used for the product of voltage phase angle difference and amplitude. Afterwards, first-order Taylor expansions are applied to the square of the voltage phase angle difference and the square of the amplitude difference, respectively, and high-order terms are ignored, reaching an approximate expression of the losses shown in Equations 32, 33. Δ P i j = g i j v i 2 + v j 2 − 2 v i v j g i j ⁡ cos ⁡ θ i j ≈ g i j v i − v j 2 + g i j θ i j 2 ≈ 2 g i j v i , 0 − v j , 0 v i − v j − g i j v i , 0 − v j , 0 2 + 2 g i j θ i j , 0 θ i j − g i j θ i j , 0 2 (32) Δ Q i j = − b i j v i 2 + v j 2 + 2 v i v j b i j ⁡ cos ⁡ θ i j ≈ − b i j v i − v j 2 − b i j θ i j 2 ≈ − 2 b i j v i , 0 − v j , 0 v i − v j + b i j v i , 0 − v j , 0 2 − 2 b i j θ i j , 0 θ i j + b i j θ i j , 0 2 (33)

The constraints for bus injection power are described in Equations 34, 35. P i = ∑ i , j ∈ Κ P i j + ∑ j = 1 N G i j v i 2 (34) Q i = ∑ i , j ∈ Κ Q i j + ∑ j = 1 N − B i j v i 2 (35) where P i and Q i represent the active and reactive injection power at bus i respectively. G i j and B i j are the elements in the bus admittance matrix.

The constraints on nodal power balance are shown in Equations 36, 37. P r e s , i − P r e s , i q u i t − P d , i = P i (36) − Q d , i = Q i (37) where P d , i and Q d , i represent the active and reactive power of the load at bus i respectively. P r e s , i and P r e s , i q u i t represent the total renewable energy power and curtailed power at bus i respectively.

The constraint for renewable energy curtailment is described in Equation 38. P r e s , i q u i t ≤ P r e s , i q , ⁡ max (38) where P r e s , i q , ⁡ max is the maximum amount of curtailed renewable power at bus i.

The renewable energy curtailment won’t happen unless the bus voltage or branch power flow reach the triggering point. When they reach the threshold for triggering curtailment, the optimization model reduces the part of bus voltage and power flow exceeding the threshold limit by curtailing part of the renewable energy, thereby reducing the overall operational risk of the distribution system. In the cases where renewable energy curtailment is called for, the distribution system status after curtailment will be calculated through power flow calculation before calculating the voltage and power flow related risks.

The various operational risk indices of the distribution system reflect different kinds of operational risks. In order to evaluate the overall risk of the system, appropriate weights need to be assigned to each risk index to form a comprehensive risk index system. In order to achieve reasonable risk assessment, this paper combines the Analytic Hierarchy Process and Entropy Weight Method. The Analytic Hierarchy Process determines index weights subjectively based on past knowledge (Saaty, 1990), while the Entropy Weight Method determines index weights objectively based on the amount of information contained in the data (Chen, 2020). The objective weights calculated by the Entropy Weight Method are objective and unbiased for it eliminates the different dimensions of various indices before comparing them. However, the method is based solely on the characteristics of data without considering the actual meaning of it, which may lead to results that go against common sense. Thus, in this paper, subjective weights with experts’ intelligence and experience contained are combined to objective weights to obtain comprehensive weights of the risk indices.

1) Subjective weights. The Analytic Hierarchy Process decomposes complex decision-making problems into a series of sub-problems, and then further divides sub-problems to build a hierarchical evaluation index system, which fully considers the relationship between various influencing factors. When the Analytic Hierarchy Process is used to determine the subjective weights, a three-layer structure of “target layer - criterion layer - scheme layer” is constructed, and the number of elements in each layer usually does not exceed nine. The structure is shown in Figure 2.

Assuming that the criterion layer has m elements, the scheme layer has n elements, and the elements in the criterion layer and the scheme layer are compared in pair-to-pair to form a judgment matrix A m × m and m n-order judgment matrices B 1 n × n ∼ B m n × n . Take matrix A as an example, and its expression is shown in Equation 39: A = a 11 … a 1 m ⋮ a i j ⋮ a m 1 ⋯ a m m (39)

In the formula, a i j represents the evaluation value obtained by comparing elements a i and a j , and the value of a i j is determined by inviting experts according to the scale shown in Table 1. If i = j, a i j = 1; If i≠j, a i j = 1/ a j i .

Consistency test was performed on the judgment matrix, and the consistency index CI is calculated in Equation 40: C I = λ max − m m − 1 (40) where λ max is the largest eigenvalue of matrix A.

The random consistency index RI is obtained according to Table 2, and then the consistency ratio CR is calculated in Equation 41: C R = C I R I (41)

If CR < 0.1, matrix A is considered to pass the consistency test, and then the weights corresponding to each judgment matrix A = a i j m × m are calculated in Equation 42. ω i = ∏ j = 1 m a i j 1 m ∑ k = 1 m ∏ j = 1 m a k j 1 m , i = 1 , 2 , . . . , m (42)

Record the weight set of matrix A as vector w A = ω 1 A , … , ω m A .

Likewise, the consistency tests of matrix B 1 n × n ∼ B m n × n are carried out respectively, and the weight vectors w B 1 ∼ w B m are calculated after passing the test. Ultimately, the subjective weight set is calculated in Equation 43. W 1 × n = ω 1 C , … , ω n C 1 × n = w A 1 × m w B 1 ⋮ w B m m × n = ω 1 A , … , ω n A ω 1 B 1 … ω n B 1 ⋮ ⋱ ⋮ ω 1 B m ⋯ ω n B m (43)

2) Objective weights. Assuming that there are a total of s risk assessment indices, and each risk assessment index has n samples. The objective weight evaluation matrix J is established in Equation 44: J = b 11 … b 1 t ⋮ b i j ⋮ b s 1 ⋯ b s t s × t (44) where element b i j represents the jth sample risk value of the ith risk assessment index. The standardized matrix C = c i j s × t is obtained by standardizing matrix J in Equation 45. c i j = b i j − min b i max b i − min b i , the index is positive max b i − b i j max b i − min b i , the index is negative (45) where max   b i represents the maximum value of the element in row i of matrix J, and min   b i represents the minimum value in row i of matrix j.

The entropy calculation formula of the risk assessment index is shown in Equation 46: E i = − 1 ln ⁡ t ∑ j = 1 t c i j ∑ j = 1 t c i j · ⁡ ln c i j ∑ j = 1 t c i j , i = 1 , 2 , . . . , s , j = 1 , 2 , . . . , t (46)

The information utility value is defined as 1 − E i and it can be normalized to obtain the objective weight set Z of s risk assessment indicators in Equation 47. Z i = 1 − E i s − ∑ i = 1 s E i , i = 1 , 2 , . . . , s (47)

After determining subjective and objective weights, they are integrated to form the comprehensive risk index weights, as is shown in Equation 48. w i = μ Z i + 1 − μ W i (48) where W is the subjective weight set obtained through the Analytic Hierarchy Process, Z is the objective weight set calculated by the Entropy Weight Method, and μ is the weight coefficient used to combine subjective and objective weights. The coefficient μ is determined based on the quality of historical data in the Entropy Weight Method and the authority of expert information in the Analytic Hierarchy Process.

After determining the weight of risk indices, the comprehensive operational risk of the distribution system can be calculated by weighted sum of overvoltage risk, low voltage risk, line overload, reverse line overload risk, and renewable energy curtailment risk, as is shown in Equation 49. C R I = w 1 R s y s v h + w 2 R s y s v l + w 3 R s y s p + w 4 R s y s pr + w 5 R s y s q u i t (49) where w 1 , w 2 , w 3 , w 4 , w 5 are the weights of each kind of operational risk respectively, and they sum up to 1.

The operational risk of the distribution system is mainly determined by the connected load and renewable energy, in addition to its own topology and network parameters. This section will discuss their impact on the operational risk of the distribution system from the perspectives of load power, renewable energy power, and the location of renewable energy. The results are shown in Figures 3–5. It can be seen that high load and renewable energy injection will bring about different kind of risks. Along with the increase of load and renewable energy, the operational risk significantly increases, and the growth is not linear. The growth rate increases with the increase of load and renewable energy power. In addition, as is shown in Figure 5, the location of renewable energy also has a significant impact on the risk of the distribution systems. When renewable energy is connected to the head of the line, its power flows to the transmission network through fewer branches, and the bus voltage in the distribution system gradually decreases from the first bus to the end without any occurrence of overvoltage. When renewable energy is connected to the middle section of the line, its power flows through more branches, resulting in a gradual increase in bus voltage from the head to the middle section of the distribution system, and overvoltage risk occurs at buses from the middle to the end. When renewable energy is connected to the end of the line, its output is transmitted backwards through the entire distribution system, posing a serious risk of overvoltage in the distribution system.

The entropy weight method is used to determine the weight of risk indices objectively based on the amount of information reflected in historical data, as a supplement to subjective weights. For different systems and when different operating modes are adopted, the risks faced by the distribution system vary, and the objective weights of the risk indices calculated by the entropy weight method also differs. This section discusses the entropy weight method from the perspective of whether renewable energy curtailment is carried out, and it can be seen from the results shown in Figure 6 that after considering the curtailment of renewable energy, the weights of the remaining four risk indices have decreased due to the addition of one risk index. However, compared to low voltage risk and overload risk, the objective weight of overvoltage risk and reverse line overload risk has decreased more significantly. This is due to the decrease of overvoltage and reverse line overload risk resulting from renewable energy curtailment. According to the Entropy Weight Method, the information contained in a risk index is positively correlated with its degree of dispersion. The wider its variation range, the more information it contains. After conducting the curtailment operations, the overvoltage risk and reverse line overload risk are reduced to the threshold limit, so the originally different severity of overvoltage and reverse line overload risks become similar, which means their fluctuation ranges are narrowed and the information contained in these two indices is reduced, thus resulting in a significant decrease in the objective weights assigned to them.

Curtailing renewable energy sources can alter the operational status of the distribution system, thereby affecting its risk level. This section discusses the effect and the results are shown in Table 5. It can be seen that after the curtailment, the overvoltage risk and reverse line overload risk caused by renewable energy have decreased. Although the renewable energy curtailment risk has significantly increased, due to the fact that it only occurs when the voltage and overload risks reach a certain level and that its value is less dispersed than the other risks, both the subjective and objective weights of renewable energy curtailment risk we obtained are relatively lower than voltage and overload risks, so the weighted total risk of the distribution system is been reduced after curtailment.

This paper utilizes the results of the prediction step of the Unscented Kalman Filter state estimation as the input for risk assessment, so the risk assessment is closely related to the accuracy of the state estimation. In order to demonstrate the relationship between state estimation and risk assessment in detail, this paper conducted comparative experiments by adjusting the variance of measurement noise. In order to simulate the measurement information of SCADA data containing noise, this paper randomly generates noise based on the proposed measurement noise variance, and adds it to the accurate information from power flow calculation to substitute the measurement data used in state estimation. Table 6 records the three noise variances used in this section, and demonstrates the impact of state estimation accuracy on risk assessment through the differences in numerical values of these three variances. In practical applications, noise variances need to be selected based on the characteristics of the measuring device. Due to the randomness of noise, the experimental results also have a certain degree of randomness. In order to better compare the results under different noises, this section conducted 10 experiments for each group, and took the average value for comparison between groups. According to the results recorded in Table 6, as the variance of measurement noise increases, the risk of the distribution system also increases. This is due to deteriorated accuracy of state estimation and the greater variance of the results of state estimation under larger measurement noise, which also means that the possible distribution system states fluctuate in a wider range, and that there is a threat of greater risks, ultimately leading to an increase in the value of comprehensive risk.