In this article, the system is assumed in a quasi-steady state operation, and a widely used dynamic model is adopted to represent the equation of state: x k + 1 = A k x k + B k + ω k (4) where A _k is the state transfer matrix, which represents the speed of the state transfer; B _k is the input matrix, reflecting the changing trend in the state trajectory. The Holt-Winters two-parameter exponential smoothing method is employed to compute A _k and B _k online. A k = α 1 + β I n B k = 1 + β 1 − α x ̂ k | k − 1 − β c k − 1 + 1 − β b k − 1 (5) subject to c k = α x ̂ k + 1 − α x ̂ k | k − 1 b k = β c k − c k − 1 + 1 − β b k − 1 (6) in which, α and β indicate smoothing parameters, taking the values between 0 and 1; x ̂ k | k − 1 is the state prediction value at k − 1; and x ̂ k is the state estimate value at k.

At present, hybrid measurements with DRTUs and PMUs are adopted in the PDNs to coordinate the economy of power network construction and the accuracy of state estimation. The use of various measurement devices increases the measurement redundancy and then provides reasonable measurement data for the PDNs state estimation. The DRTU measurements are expressed as U i p = U i r e , p 2 + U i i m , p 2 , (7) P i − j p = ∑ q ∈ a , b , c G i − j p q U i r e , p U ⃗ i − j r e , q + U i i m , p U ⃗ i − j i m , q + B i − j p q U i i m , p U ⃗ i − j r e , q − U i r e , p U ⃗ i − j i m , q Q i − j p = − b i − j p 2 U i r e , p 2 + U i i m , p 2 + ∑ q ∈ a , b , c G i − j p q U i i m , p U ⃗ i − j r e , q − U i r e , p U ⃗ i − j i m , q − B i − j p q U i i m , p U ⃗ i − j i m , q + U i r e , p U ⃗ i − j r e , q , (8) P i p = ∑ j = 1 , j ≠ i j = N γ j − i P j − i p Q i p = ∑ j = 1 , j ≠ i j = N γ j − i Q j − i p (9) where U ⃗ i − j r e , q = U i r e , q − U j r e , q , U ⃗ i − j i m , q = U i i m , q − U j i m , q , P i − j p and Q i − j p are the p-phase active and reactive power flow measurements at the i − j branch, respectively. P i p and Q i p are the p-phase active and reactive power injection measurements at the ith bus, respectively. b i − j p represents p-phase ground guide at the i − j branch. γ _j−i is a scalar describing the mutual coupling between the network nodes. If Θ represents the set of neighbor nodes of node i, then γ j − i = 1 , j ∈ Θ 0 , j ∉ Θ (10)

Besides, the PMU measurements are expressed as U i r e , p = U i r e , p U i i m , p = U i i m , p , (11) I i − j r e , p = − b i − j p 2 U i i m , p + ∑ q ∈ a , b , c G i − j p q U ⃗ i − j r e , q − B i − j p q U ⃗ i − j i m , q I i − j i m , p = b i − j p 2 U i r e , p + ∑ q ∈ a , b , c G i − j p q U ⃗ i − j i m , q + B i − j p q U ⃗ i − j r e , q (12) where I i − j r e , p and I i − j i m , p represent the p-phase real and imaginary currents at the i − j branch, respectively.

Figure 1 shows the three phase unbalanced line model of the PDNs, where the missing phase line model is similar to the three-phase model and is not fully shown here.

Actually, the linearization error boundary determined by the interval analysis method is too conservative, and its cumulative effect will also affect the accuracy and even lead to filter divergence. In addition, the Lagrange remainder of the linearized nonlinear function can be used to represent the linearization error, but the process of optimizing the Hessian matrix is more complex. In this article, the linearization error is limited by solving the dual form of semi-infinite programming, which not only avoids the solution of Hessian matrix but also overcomes the defects of the interval analysis method, and a more compact linearization error boundary can be obtained.

A lemma will first be introduced for the formula derivation in this section.

Lemma 1

Wang et al. (2021) The optimization problem of the ellipsoid can be expressed as min P , u ̂ log ⁡ det P s . t . u − u ̂ T P − 1 u − u ̂ ⩽ 1 , ∀ u ∈ C (22) with C ≜ u : u = s x , x ∈ ϕ ϕ ≜ x : x − x ̂ T P ̃ − 1 x − x ̂ ⩽ 1 . (23)

Then, the optimal solution of Eq. 22 is P * = ∫ C u u T d τ * − u * u * T (24) where u ̂ * = ∫ C u d τ * is the optimal solution to the following optimization problem: max τ i ⩾ 0 log ⁡ det ∫ C × 1 u ⌣ u ⌣ T d τ s . t . ∫ C × 1 d τ = 1 (25) with u ⌣ = [ u T , 1 ] T , and Eq. 25 is the dual problem of Eq. 22. For the proof of this lemma, please refer to Wang et al. (2021).

Next, we introduce the specific application of this lemma in this article. First, the Taylor expansion of the measurement equation is performed around x ̂ k + 1 | k . g ( x k + 1 ) = g ( x ̂ k + 1 | k ) + G k + 1 ( x k + 1 − x ̂ k + 1 | k ) + Δ g ( x k + 1 ) (26)

According to the properties of the ellipsoid, Eq. 26 is equivalent to: g x ̂ k + 1 | k + E k + 1 | k z k + 1 | k = g x ̂ k + 1 | k + G k + 1 E k + 1 | k z k + 1 | k + Δ g z k + 1 | k (27) where z k + 1 | k ≤ 1 , P k + 1 | k = E k + 1 | k E k + 1 | k T . Thus, the higher-order linearization error can be expressed as Δ g z k + 1 | k = g x k + 1 − g x ̂ k + 1 | k − G k + 1 x k + 1 − x ̂ k + 1 | k = g x ̂ k + 1 | k + E k + 1 | k z k + 1 | k − g x ̂ k + 1 | k − G k + 1 E k + 1 | k z k + 1 | k . (28)

As shown in Figure 3, we try to find a compact ellipsoid to cover the higher-order remainder term, assuming that the outsourced ellipsoid of the higher-order remainder term is expressed as ϕ _g,k+1. This article uses the nonlinear programming method to optimize the boundary of linearization error, which can get a tight boundary estimate. In general, it can be obtained from Lemma 1 that the linear programming Eq. 25 is the dual form of the nonlinear programming Eq. 22. Therefore, in order to achieve a balance between the estimation accuracy and the calculation amount, the continuous system is discretized by the Monte Carlo sampling method. In this regard, the nonlinear programming is converted into a linear programming problem. The following expression is established: Δ g z k + 1 | k ≜ u | u − u ̂ g , k + 1 T P g , k + 1 − 1 u − u ̂ g , k + 1 ≤ 1 ≜ u | u = u ̂ g , k + 1 + E g , k + 1 | k z g , k + 1 | k , P g , k + 1 = E g , k + 1 | k E g , k + 1 | k T , z g , k + 1 | k ≤ 1 ∈ ϕ g , k + 1 (29) where u ̂ g , k + 1 and P _g,k+1 denote the center and shape matrix of the ellipsoid ϕ _g,k+1, respectively.

A nonlinear programming method can be used to optimize the linearized error boundary, so the following nonlinear programming problem can be constructed: min P g , k + 1 , u ̂ g , k + 1 log ⁡ det P g , k + 1 s . t . u − u ̂ g , k + 1 T P g , k + 1 − 1 u − u ̂ g , k + 1 ⩽ 1 , ∀ u ∈ C (30) where C is defined as C ≜ u : u = h x k + 1 , x k + 1 ∈ ϕ k + 1 | k ϕ k + 1 | k ≜ x k + 1 : x k + 1 − x ̂ k + 1 | k T P k + 1 | k − 1 x k + 1 − x ̂ k + 1 | k ⩽ 1 (31) in which, h(⋅) is an arbitrary nonlinear continuous function, expressed here as a higher-order remainder term function.

The minimal volume of the linear error outsourcing ellipsoid can be obtained by solving the nonlinear programming problem. In order to make it more convenient, the continuous system is discretized by the Monte Carlo sampling method, and the discrete expression of the center and shape matrix of the ellipsoid are obtained by Eqs. 32, 33, which corresponds to Eq. 24 under continuous form. Moreover, the approximate expression of Eq. 25 is obtained, as shown in Eq. 34. Thus, an outsourcing compact ellipsoid with the linearization error is developed. The results obtained in this section will be used in the design of the dual set-membership filter. The specific process will be shown in detail in Section 4.2. P g , k + 1 = ∑ i = 1 m τ i u i u i T − u ̂ g , k + 1 u ̂ g , k + 1 T (32) with u ̂ g , k + 1 = ∑ i = 1 m τ i u i (33) where coefficients τ _i (i = 1, … , m) are obtained by solving the linear programming shown in the solution Eq. 34 max τ i ⩾ 0 log ⁡ det ∑ i m τ i u ⌣ i u ⌣ i T s . t . ∑ i m τ i = 1 . (34)

By Lemma 1, we obtain the dual form of nonlinear programming and the continuous expressions of the center and shape matrix of the linearized error ellipsoid. Then, the continuous system is discretized by the Monte Carlo sampling method, and the semi-infinite programming problem is transformed into a finite programming problem, so the discrete expressions of the ellipsoid center and shape matrix are obtained.

Remark 2

Different from the UKF, the approximate probability distribution of the nonlinear function cannot be obtained directly, since we assume that the true state is restricted to the ellipsoid, rather than obeying a particular Gaussian distribution. Wang et al. (2021) avoids the process of Taylor expansion by directly solving the approximate distribution of the nonlinear measurement function values. But this will still inevitably produce large errors for strong nonlinear systems. Therefore, the method in this article strikes a balance between the computational quantity and the estimation accuracy.

In this section, a set-membership filtering algorithm based on the event-triggered mechanism is proposed to solve the state estimation problem of the PDNs. And the real state is always confined to the set of state estimation. In order to deal with quadratic constraints such as event-triggered mechanism and external noises, RLMI technology is used to transform the externally bounded ellipsoid filtering problem into a convex optimization problem with linear matrix inequality constraints. The sufficient conditions for the existence of the estimated state set are strictly derived, and the filter gain parameters are easily obtained by solving the convex optimization problem(Boyd et al., 1994).

We begin by introducing two lemmas, which will be used in the proof of this section.

Lemma 2

(S-Procedure) (Bai et al., 2021a): Let ℏ ₀(⋅), ℏ ₁(⋅), …, ℏ _p (⋅) be the quadratic function of the variable χ ∈ R n , ℏ _i (χ) = χ ^T M _i χ(i = 0, 1, … , p), with M i = M i T , if there exist λ ₁⩾0, … , λ _p ⩾0 such that M 0 − ∑ i = 1 p λ i M i ⩽ 0 , then the following is true: ℏ ₁(χ)⩽0, … , ℏ _p (χ)⩽0 → ℏ ₀(χ)⩽0.

Lemma 3

(Schur Complement Equivalence) (Bai et al., 2021a): For matrices H ₁, H ₂, H ₃, where S 1 = S 1 T ⩾ 0 , H 2 = H 2 T ⩾ 0 , H 3 T H 2 − 1 H 3 + H 1 ⩽ 0 , if and only if − H 2 H 3 H 3 T H 1 < 0 or H 1 H 3 T H 3 − H 2 < 0 .

One-step state prediction value x ̂ k + 1 | k is approximated by the Holt-Winters two-parameter exponential smoothing method. Based on the state estimate at k, one-step predicted ellipsoid ϕ _k+1|k is obtained by solving the semi-definite programming problem Eq. 36. ϕ k | k + 1 = x k + 1 : x k + 1 − x ̂ k + 1 | k T P k + 1 | k − 1 x k + 1 − x ̂ k + 1 | k ⩽ 1 (35) where P _k+1|k is the desired optimization variable.

It should be noted that calculating the volume of the ellipsoidal set is more difficult than the trace of the matrix P _k+1. The trace of P _k+1 is the sum of the semi-axes of the ellipsoid, by which the size of the ellipsoidal set can be reflected effectively. Therefore, the objective function Eq. 36 is the trace function. min P k + 1 | k , λ 1 ⩾ 0 , λ 2 ⩾ 0 f P k + 1 | k s . t . Ξ k + 1 | k Φ k + 1 | k T Φ k + 1 | k − P k + 1 | k < 0 (36) where Φ k + 1 | k = 0 , A k E k , I , (37) Ξ k + 1 | k = diag λ 1 + λ 2 − 1 , − λ 1 I , − λ 2 Q k − 1 . (38)

Proof 1

In view of Equation 4, we can get one-step prediction state error e k + 1 | k = x k + 1 − x ̂ k + 1 | k = A k E k z k + w k = Φ k + 1 | k ζ . (39)

Let, ζ = [ 1 , z k T , w k T ] T . Thus, x _k+1 ∈ ϕ _k+1|k is equivalent to ζ T Φ k + 1 | k T P k + 1 | k − 1 Φ k + 1 | k ζ ≤ 1 . (40)

The following unequal conditions exist naturally: z k ≤ 1 w k T Q k − 1 w k ≤ 1 . (41)

Eq. 41 is equivalent to ζ T diag − 1 , I , 0 ζ ≤ 0 ζ T diag − 1,0 , − Q k − 1 ζ ≤ 0 . (42)

By means of Lemma 1, ζ T [ Φ k + 1 | k T P k + 1 | k − 1 Φ k + 1 | k ] ζ ≤ 1 holds, if positive scalar λ ₁, λ ₂ exists such that Φ k + 1 | k T P k + 1 | k − 1 Φ k + 1 | k − diag 1,0,0 − λ 1 diag − 1 , I , 0 − λ 2 diag − 1,0 , − Q k − 1 ≤ 0 . (43)

Equation 43 can be rewritten as Φ k + 1 | k T P k + 1 | k − 1 Φ k + 1 | k + Ξ k + 1 | k ≤ 0 (44) where Ξ k + 1 | k = diag λ 1 + λ 2 − 1 , − λ 1 I , − λ 2 Q k − 1 . (45)

By Lemma 2, Equation 45 can be converted to the following formula: Ξ k + 1 | k Φ k + 1 | k T Φ k + 1 | k − P k + 1 | k < 0 . (46)

The proof is completed.

On the other hand, based on the proposed event-triggered mechanism, the estimates can be expressed by: x ̂ k + 1 = A k x ̂ k + B k + K k + 1 y ̃ k + 1 − g x ̂ k + 1 | k (47) where K _k+1 is the filter gain matrix to be solved. In addition, the estimated ellipsoid ϕ _k+1 is obtained by solving the semi-definite programming problem, based on the one-step state prediction value at k. ϕ k + 1 = x k + 1 : x k + 1 − x ̂ k + 1 T P k + 1 − 1 x k + 1 − x ̂ k + 1 ≤ 1 (48)

The optimal variable P _k+1 and the gain matrixes K _k+1 can be obtained by solving the following optimization problem: min P k + 1 , K k + 1 , λ 3 ⩾ 0 , … , λ 7 ⩾ 0 f P k + 1 s . t . Ξ k + 1 Φ k + 1 T Φ k + 1 − P k + 1 < 0 (49) where Φ k + 1 = − K k + 1 e g , k + 1 , I − K k + 1 G k + 1 A k E k , I − K k + 1 G k + 1 , − K k + 1 E g , k + 1 , − K k + 1 , − K k + 1 (50) Ξ k + 1 = diag λ 3 + λ 4 + λ 5 + λ 6 + λ 7 δ − 1 , − λ 3 I , − λ 4 Q k − 1 , − λ 5 I , − λ 6 R k + 1 − 1 , − λ 7 I (51)

Proof 2

According to Eq. 47, e k + 1 = x k + 1 − x ̂ k + 1 = x k + 1 − x ̂ k + 1 | k − K k + 1 G k + 1 x k + 1 − x ̂ k + 1 | k + e g , k + 1 + E g , k + 1 z g , k + 1 + v k + 1 + σ k + 1 = I − K k + 1 G k + 1 A k E k z k + w k − K k + 1 e g , k + 1 − K k + 1 E g , k + 1 z g , k + 1 − K k + 1 v k + 1 − K k + 1 σ k + 1 = − K k + 1 e g , k + 1 + I − K k + 1 G k + 1 A k E k z k + I − K k + 1 G k + 1 w k − K k + 1 e g , k + 1 − K k + 1 E g , k + 1 z g , k + 1 − K k + 1 v k + 1 − K k + 1 σ k + 1 = Φ k + 1 η (52) with η = 1 , z k T , w k T , z g , k + 1 T , v k + 1 T , σ k + 1 T T . (53)

Thus, x _k+1 ∈ ϕ _k+1 is equivalent to η T [ Φ k + 1 T P k + 1 − 1 Φ k + 1 ] η ⩽ 1 . The following conditions are naturally satisfied: z k ⩽ 1 w k T Q k − 1 w k ⩽ 1 z g , k + 1 ⩽ 1 v k + 1 T R k + 1 − 1 v k + 1 ⩽ 1 σ k + 1 T σ k + 1 ⩽ δ (54) which can be rewritten as follows: η T diag − 1 , I , 0,0,0,0 η ⩽ 0 η T diag − 1,0 , Q k − 1 , 0,0,0 η ⩽ 0 η T diag − 1,0,0 , I , 0,0 η ⩽ 0 η T diag − 1,0,0,0 , R k + 1 − 1 , 0 η ⩽ 0 η T diag − δ , 0,0,0,0 , I η ⩽ 0 . (55)

It can be seen from Lemma 1 that Eq. 48 is established if positive scalar λ ₃, λ ₄, λ ₅, λ ₆, λ ₇ exists to make the following inequality hold: Φ k + 1 T P k + 1 − 1 Φ k + 1 − diag 1,0,0,0,0,0 − λ 3 diag − 1 , I , 0,0,0,0 − λ 4 diag − 1,0 , Q k − 1 , 0,0,0 − λ 5 diag − 1,0,0 , I , 0,0 − λ 6 diag − 1,0,0,0 , R k + 1 − 1 , 0 − λ 7 diag − δ , 0,0,0,0 , I ⩽ 0 . (56)

Equation 45 is equivalent to Φ k + 1 T P k + 1 − 1 Φ k + 1 + Ξ k + 1 ⩽ 0 (57) where Ξ k + 1 = diag λ 3 + λ 4 + λ 5 + λ 6 + λ 7 δ − 1 , − λ 3 I , − λ 4 Q k − 1 , − λ 5 I , − λ 6 R k + 1 − 1 , − λ 7 I . (58)

Using Lemma 2, the following inequality relationship holds: Ξ k + 1 Φ k + 1 T Φ k + 1 − P k + 1 < 0 . (59)

The proof is completed.

In order to more clearly show the filtering algorithm in this article, the ET-DSMF algorithm is summarized as in Algorithm 1.

Remark 3

It is worth mentioning that the filtering algorithms under the Kalman framework are studied based on the assumption that the noises obey Gaussian distribution (Yang and Li, 2009). However, because the distribution network is located at the end of the power system and there is serious electromagnetic interference, the above assumption is not in line with the actual project. In fact, compared with the statistical characteristics of noises, it is easier for us to know its boundaries. Fortunately, set-membership filtering can effectively deal with this unknown but bounded noises and provide an estimation set that always surrounds the real states.

Remark 4

In order to alleviate communication pressure and improve the induction phenomenon of measurement data in the transmission process, the event-triggered mechanism is integrated into the design of the set-membership filter. To tackle the effects of σ _k , we employ the prescribed triggering condition which can be rewritten as a quadratic constraint. Subsequently, by using S-procedure, this quadratic constraint can be easily dealt with and reflected in the obtained RLMIs. And the gain matrix and the filter of the optimal estimated ellipsoid are obtained by solving the convex optimization problem.

Figures 5–7 show the performance of the DSMF and ESMF algorithms, taking the simulation results of a-phase at bus 4 as an example. It can be seen that both algorithms can achieve good estimation results. However, the MSE generated by the method proposed in this article is smaller. As shown in Figure 7, although the estimated boundaries obtained by the two algorithms can contain the real state, the proposed algorithm is less conservative owing to considering the compact linearization error ellipsoid. It is obvious that the boundary obtained by the DSMF algorithm is more compact than that obtained by the ESMF.

Firstly, the event-triggering threshold is set as ɛ _i = 5 × 10^–6. Figure 8 shows the state estimation ellipsoid center and boundary of the ET-DSMF algorithm, where κ = 83%. Although the estimation performance of the proposed algorithm degrades at a lower transmission rate, it ensures that the real state always resides in the estimation ellipsoid because the effective information is transmitted to the estimation center. In order to show the filtering results more clearly, Figure 9 shows the MSEs generated by the proposed algorithm. It can be seen that the estimation error is always below the minimum upper bound and close to the upper bound, which further shows the low conservatism of the proposed algorithm.

In order to verify the influence of the different triggering thresholds on the performance of the ET-DSMF algorithm, different triggering threshold parameters are selected, i.e., ET-DSMF1 with ɛ _i = 4 × 10^–6, ET-DSMF2 with ɛ _i = 6 × 10^–6, and ET-DSMF3 with ɛ _i = 8 × 10^–6. Figure 10 and Figure 11 show the filtering curves and MSEs of c-phase at bus 4 under three different threshold parameters. In addition, to more clearly show the triggering time, the triggering sequences is shown in Figure 12. With the increase of the triggering threshold, the MSE of the ET-DSMF algorithm increases and the data transmission data decreases accordingly. Therefore, selecting the appropriate triggering threshold parameters can effectively alleviate the pressure of network communication and reduce the probability of network-induced phenomenon.