Decentralized State Estimation of Combined Heat and Power System Considering Communication Packet Loss

Wenjian Zheng; Zhigang Li; Xinyu Liang; Jiehui Zheng; Q. H. Wu; Fan Hu

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Decentralized State Estimation of Combined Heat and Power System Considering Communication Packet Loss PDF

- ORCID：
Wenjian Zheng ¹
- ORCID：
Zhigang Li ¹
- ORCID：
Xinyu Liang ¹
- ORCID：
Jiehui Zheng ¹
- ORCID：
Q. H. Wu ¹
- ORCID：
Fan Hu ²

1. School of Electric Power Engineering, South China University of Technology, Guangzhou 510641, China； 2. Guangzhou Power Supply, Guangdong Power Grid Co., Ltd., Guangzhou 510000, China

Updated：2020-07-10

DOI：10.35833/MPCE.2020.000120

OUTLINE

Abstract

In order to obtain an accurate state estimation of the operation in the combined heat and power system, it is necessary to carry out state estimation. Due to the limited information sharing among various energy systems, it is practical to perform state estimation in a decentralized manner. However, the possible communication packet loss is seldomly considered among various energy systems. This paper bridges this gap by proposing a relaxed alternating direction method of multiplier algorithm. It can also improve the computation efficiency compared with the conventional alternating direction of the multiplier algorithm. Case studies of two test systems are carried out to show the validity and superiority of the proposed algorithm.

Keywords

Alternating direction of multiplier algorithm; combined heat and power system; communication packet loss; decentralized optimization; state estimation

I. Introduction

WITH the increasing pressure from energy and environmental protection, the energy utilization efficiency is improved recently, which has promoted the transformation from the conventional single energy system to the integrated energy system (IES) with the harmonization of multi-energy [

1]. One form of the IES is the combined heat and power system (CHPS). It creates the linkage between the electric power systems (EPSs) and the district heating networks (DHNs) through the coupling components such as combined heat and power (CHP) units, gas boilers and heat pumps. Due to the high fuel efficiency of CHP units, the CHPS can promote the efficient and flexible utilization of energy [2], [3].

The energy management is based on the system operation states. In order to obtain an accurate estimate of the operation states in CHPS, it is necessary to perform state estimation (SE). SE is based on the measurement data from metering devices. Due to the measurement incompleteness and the random measurement error of metering devices, SE is an important part of the monitoring and operation of the system. SE of EPSs and DHNs is widely investigated [

4], [5]. However, there is few research on CHPS. Reference [6] considers the coupling components in the combined network and realizes SE for CHPS. However, the DHN model is oversimplified and the dynamic characteristics of the water pipelines are ignored. To estimate the actual operation states more accurately, a two-stage SE approach is proposed in [7] considering the dynamic property of the pipelines. An alternating estimation strategy is adopted in [8] to cope with the quasi-dynamics of temperature in the pipelines. On this basis, considering the limited information exchange among different energy systems, [9] proposes a decentralized algorithm based on the asynchronous alternating direction method of multipliers (ADMM) to improve the computation efficiency. Reference [10] applies the ADMM algorithm in SE of IES. Besides, it develops a bilinear measurement model to avoid non-convexity in the optimization. To suppress the negative effects of bad data, a heating load pseudo-measurement model based on artificial neural network is proposed in SE of CHPS [11].

In the above literatures, SE is carried out in a centralized or a decentralized manner. However, it is noteworthy that different energy systems are operated by different energy system operators. The measurement information among different energy systems cannot be shared fully due to the privacy protection and the communication unreliability among different subsystems. Thus, the decentralized solution strategy is more practical for SE in CHPS.

The decentralized solution strategy decomposes the primal optimization problem into several subproblems. In terms of subproblem solutions, the decentralized optimization techniques can be divided into two categories [

12]. One is based on the optimality condition decomposition, which applies the Newton-Raphson method to the Karush-Kuhn-Tucker condition of each subproblem to obtain the incremental subproblem solutions [13]. The other is based on the augmented Lagrangian decomposition, which introduces the multiplier variables to relax the coupling constraints. It solves the subproblems and updates the multipliers in each iteration. As a representative, the ADMM algorithm with great performance on computation and convergence has been applied in the decentralized SE problem. In [14], the ADMM algorithm is adopted to perform SE for EPSs with linear measurements. The algorithm requires low communication burden and exhibits great convergence performance in the simulations, but it does not consider the applicability on the non-convex problems. In [15], considering the non-convex issue of nonlinear measurements, a distributed bilinear SE procedure based on the ADMM is developed to multi-area EPSs. As mentioned earlier, the ADMM algorithm is applied in SE for IES [9], [10].

Although the decentralized SE algorithms mainly focus on the practicability and validity, the possible communication failures among different subsystems are randomly considered. Recently, the communication failures have been gaining a lot of attention because of their negative effects on the operation of EPSs [

16]. The communication failures will cause measurements loss, which reduces the accuracy of SE. To address this problem, the pseudo-measurements [17] are made based on its correlation with the real-time measurements to estimate the lost measurements. In [18], a new dynamic SE approach is proposed based on the load forecasting technique and extended Kalman filter, which can predict and estimate system states with partial communication failures. The aforementioned literatures focus on the missing measurements due to the communication failures and SE is carried out in a centralized way. However, the communication failures may occur in the process of information exchange among subsystems due to the instability of the communication channels. To handle this problem, a linear matrix inequality approach is developed for distributed dynamic SE in [19]. This approach considers packet loss in the interconnected smart grid subsystems. In addition, a relaxed ADMM (R-ADMM) algorithm for distributed optimization over lossy networks is proposed in [20]. The algorithm is shown to be robust against the random packet loss in the neighboring areas. In the above literatures, few studies investigate SE of CHPS considering the communication packet loss among subsystems. To bridge this gap, the R-ADMM algorithm is further developed for decentralized SE in this paper. The contributions of this paper are summarized as follows:

1) An R-ADMM algorithm is developed for the decentralized SE of the CHPS. The algorithm shows faster convergence rates than the conventional ADMM algorithm, and it improves the commutation efficiency of the decentralized solving procedure.

2) The robust implementation of the R-ADMM algorithm is proposed considering communication failures among different subsystems. The algorithm can still converge and maintain great computation efficiency with high probability of packet loss, which reflects its robustness.

The remainder of this paper is organized as follows. Section II provides SE formulation of CHPS. Section III presents the R-ADMM algorithm. Section IV presents the case study. Section V draws the conclusions and suggests the future works.

II. Formulation of SE for CHPS

In this section, the SE problem for CHPS is formulated. The structure of an exemplary CHPS is shown in Fig. 1. The CHPS creates the linkage between the EPSs and the DHNs by the coupling units such as CHP units, gas boilers and heat pumps.

Fig. 1 Structure of an exemplary CHPS.

The formulation of SE for CHPS [

9] is adopted in this paper, which considers the quasi-dynamics of temperature in the DHN pipelines that yields more accurate estimate results.

A. SE Formulation of EPS

1)　Objective Function

Weighted least squares (WLS) is widely used for SE. The objective function based on WLS criterion is to minimize the weighted sum of squares of the measurement residuals. Therefore, the objective function of the SE in the EPS is expressed as:

$J^{E P S} (V, θ, P^{i n j}, Q^{i n j}) = \sum_{i \in 𝒩^{V}} \frac{({\tilde{V}}_{i} - V_{i})^{2}}{W_{i}^{V}} + \sum_{i \in 𝒩^{P^{i n j}}} \frac{({\tilde{P}}_{i}^{i n j} - P_{i}^{i n j})^{2}}{W_{i}^{P^{i n j}}} + \sum_{i \in 𝒩^{Q^{i n j}}} \frac{({\tilde{Q}}_{i}^{i n j} - Q_{i}^{i n j})^{2}}{W_{i}^{Q^{i n j}}}$

(1)

where $V$ is the bus voltage magnitude; $θ$ is the bus voltage phase angle; $P^{i n j}, Q^{i n j}$ are the power injections; $V_{i}$ is the voltage magnitude; $W_{i}^{V}$ is the weight coefficient of the measurement of voltage magnitude; $W_{i}^{P^{i n j}}$ , $W_{i}^{Q^{i n j}}$ are the weight coefficients of the measurement of active and reactive power injections, respectively; $𝒩^{V}$ is the index set of measurements of voltage magnitude; and $𝒩^{P^{i n j}}$ , $𝒩^{Q^{i n j}}$ are the index sets of measurements of active and reactive power injections, respectively. In this paper, the variables with a tilde above them represent measurements of those variables and the corresponding variables without a tilde above them represent the estimated values.

2)　Constraints

The constraints are the measurement equations that establish the relationship between the measurements and the state variables. For the EPS, the power injection equations are shown as:

$P_{i, t} = V_{i} \sum_{j \in ℐ^{b u s}} V_{j} (G_{i j} c o s θ_{i j} + B_{i j} s i n θ_{i j}) \forall i \in ℐ^{b u s}, t \in 𝒯$

(2)

$Q_{i, t} = V_{i} \sum_{j \in ℐ^{b u s}} V_{j} (G_{i j} s i n θ_{i j} + B_{i j} c o s θ_{i j}) \forall i \in ℐ^{b u s}, t \in 𝒯$

(3)

where $P_{i, t}, Q_{i, t}$ are the active and reactive power injections at bus i in period t, respectively; $G_{i j}$ is the $i j^{t h}$ element of the conductance matrix; $B_{i j}$ is the $i j^{t h}$ element of the susceptance matrix; $ℐ^{b u s}$ is the index set of buses in the EPS; and $𝒯$ is the index set of time periods.

B. SE Formulation of DHN

1)　Objective Function

The weighted sum of squares of the measurement residuals for DHN is minimized as:

$\begin{array}{l} J^{D H N} (p s, m, Φ, T) = \sum_{i \in 𝒩^{p s}} \frac{({\tilde{p s}}_{i} - p s_{i})^{2}}{W_{i}^{p s}} + \sum_{i \in 𝒩^{m}} \frac{({\tilde{m}}_{i} - m_{i})^{2}}{W_{i}^{m}} + \\ \sum_{i \in 𝒩^{Φ}} \frac{({\tilde{Φ}}_{i} - Φ_{i})^{2}}{W_{i}^{Φ}} + \sum_{i \in 𝒩^{T}} \frac{({\tilde{T}}_{i} - T_{i})^{2}}{W_{i}^{T}} \end{array}$

(4)

where $p s$ is the water pressure; $m$ is the mass flow rate; $Φ$ is the heat power; $T$ is the temperature; $𝒩^{p s}$ , $𝒩^{m}$ , $𝒩^{Φ}$ , $𝒩^{T}$ are the index sets of pressure measurements, mass flow rate measurements, heat power measurements, and temperature measurements, respectively; and $W_{i}^{p s},$ $W_{i}^{m},$ $W_{i}^{Φ},$ $W_{i}^{T}$ are the weight coefficients of the measurements of the pressure, mass flow rate, heat power, and temperature, respectively. $p s$ , $m$ , $Φ$ , and $T$ are the decision variables.

2)　Constraints

The constraints in DHN depict the relationship among the decision variables above. They are given below.

1) Continuity of mass flow: the mass flow that enters into a node is equal to the mass flow leaves the node. The constraint is expressed as:

$\sum_{b \in 𝒮_{i}^{P i p e +}} m_{b, t}^{P} = \sum_{b \in 𝒮_{i}^{P i p e -}} m_{b, t}^{P} \forall i \in ℐ^{n o d e}, t \in 𝒯$

(5)

where $𝒮_{i}^{P i p e +},$ $𝒮_{i}^{P i p e -}$ are the index sets of pipelines that end at the node and start from the node, respectively; $m_{b, t}^{P}$ is the mass flow rate of pipeline b in period t; and $ℐ^{n o d e}$ is the index set of the nodes in the DHN.

2) Head loss: head loss is the pressure change due to the pipeline friction. The relation between the mass flow and the head losses along each pipeline is:

$p s_{b, t}^{i n} - p s_{b, t}^{o u t} = K^{f} m_{b, t}^{P} |m_{b, t}^{P}| \forall b \in ℐ^{p i p e}, t \in 𝒯$

(6)

where $p s_{b, t}^{i n}, p s_{b, t}^{o u t}$ are the pressures at the inlet and outlet of pipeline b in period t, respectively; $K^{f}$ is the pressure loss coefficient of a pipeline; and $I^{p i p e}$ is the index set of pipelines in the DHN.

3) Heat power: the heat power at a load node or a source node is determined by the mass flow rates and the temperature changes, expressed by:

$Φ_{i, t}^{L D} = C_{p} m_{i, t}^{L D} (T_{i, t}^{L D, i n} - T_{i, t}^{L D, o u t}) \forall i \in ℐ^{l o a d}, t \in 𝒯$

(7)

$Φ_{i, t}^{C H P} = C_{p} m_{i, t}^{S} (T_{i, t}^{S, i n} - T_{i, t}^{S, o u t}) \forall i \in ℐ^{C H P}, t \in 𝒯$

(8)

where $Φ_{i, t}^{L D},$ $Φ_{i, t}^{C H P}$ are the heat power supplied to load node and heat source in period t, respectively; $C_{p}$ is the specific heat capacity of water; $m_{i, t}^{L D},$ $m_{i, t}^{S}$ are the mass flow rates entering/exiting load node i and source node i in period t, respectively; $T_{i, t}^{L D, i n},$ $T_{i, t}^{L D, o u t}$ are the temperatures of the mass flow entering load node i and exiting load node i in period t, respectively; $T_{i, t}^{S, i n}$ , $T_{i, t}^{S, o u t}$ are the temperatures of the mass flow entering source node i and exiting source node i in period t, respectively; and $I^{l o a d},$ $I^{C H P}$ are the index sets of loads and CHP units in the DHN, respectively.

4) Mixing temperature: for the temperature of the mass flow leaving a node with more than one pipeline, the mixing temperature is given as:

$\sum_{b \in 𝒮_{i}^{p i p e +}} m_{b, t}^{P} T_{b, t}^{P, o u t} = T_{i, t}^{N D} \sum_{b \in 𝒮_{i}^{p i p e +}} m_{b, t}^{P} \forall i \in ℐ^{n o d e}, t \in 𝒯$

(9)

where $T_{b, t}^{P, o u t}$ is the temperature of the mass flow at the outlet of pipeline b in period t considering the heat loss; and $T_{i, t}^{N D}$ is the mixing temperature at node i in period t.

The inlet temperature of mass flow is equal to the mixing temperature at the node as:

$T_{b, t}^{P, i n} = T_{i, t}^{N D} \forall b \in 𝒮_{i}^{p i p e -}, t \in 𝒯$

(10)

where $T_{b, t}^{P, i n}$ is the temperature of the mass flow at the inlet of pipeline b in period t considering the heat loss.

5) Heat dynamics and heat loss: in the pipeline, the outlet temperature responds to the inlet temperature with the time delays [

21]. Due to the difference between the water temperature and the ambient temperature, the temperature of the water drops as the water flows. To simulate the heat dynamics and the heat losses more accurately, we adopt the method in [22] in this paper.

As shown in Fig.2, the outlet temperature without heat loss can be denoted by the average temperature of the water masses flowing out of the pipeline during the period interval $Δ t$ .

Fig. 2 Vertical section of a pipeline in DHN.

It can be estimated by a number of historical inlet temperatures and the corresponding time-delay coefficients as follows:

$T_{b, t}^{P, o u t *} = \sum_{k = t - ϕ_{b, t}}^{t - γ_{b, t}} K_{b, t, k}^{T D} T_{b, k}^{P, i n} \forall b \in ℐ^{p i p e}, t \in 𝒯$

(11)

where $t - γ_{b, t}, t - ϕ_{b, t}$ are the indexes of the last period whose water mass flows out of the pipeline before the end of period $t$ and $t - 1$ , respectively; and $K_{b, t, k}^{T D}$ is the time-delay coefficient. These variables are determined by the pipeline parameters and the mass flow rate as:

$K_{b, t, k}^{T D} = \{\begin{array}{l} (m_{b, t}^{P} Δ t + ρ A_{b} L_{b} - S_{b, t}) / m_{b, t}^{P} & k = t - ϕ_{b, t} \\ m_{b, k}^{P} Δ t / m_{b, t}^{P} & k = t - ϕ_{b, t} + 1, \dots, t - γ_{b, t} - 1 \\ (R_{b, t} - ρ A_{b} L_{b}) / m_{b, t}^{P} & k = t - γ_{b, t} \\ 0 & o t h e r w i s e \end{array}$

(12)

$γ_{b, t} = \underset{x}{m i n} \{x : \sum_{k = t - x}^{t} m_{b, k}^{P} Δ t \geq ρ A_{b} L_{b}, x \geq 0, x \in ℤ\}$

(13)

$ϕ_{b, t} = \underset{y}{m i n} \{y : \sum_{k = t - y}^{t} m_{b, k}^{P} Δ t \geq ρ A_{b} L_{b} + m_{b, t}^{P} Δ t, x \geq 0, x \in ℤ\}$

(14)

where $ρ$ is the density of water; $R_{b, t}$ is the total mass flowing into the pipeline from the time $t - γ_{b, t}$ to $t$ , as shown in Fig. 2; $S_{b, t}$ is the total mass flowing into the pipeline from the time $t - γ_{b, t}$ to $t$ ; $∆ t$ is the time interval per period; $A_{b}$ is the cross-sectional area of the pipeline; and $L_{b}$ is the length of the pipeline. These variables can be defined as:

$R_{b, t} = \sum_{k = γ_{b, t}}^{t} m_{b, k}^{P} Δ t$

(15)

$S_{b, t} = \{\begin{array}{l} \sum_{k = t - ϕ_{b, t} + 1}^{t} m_{b, k}^{P} Δ t & ϕ_{b, t} \geq γ_{b, t} + 1 \\ R_{b, t} & o t h e r w i s e \end{array}$

(16)

Besides, the outlet temperature considering the heat loss is calculated as:

$\begin{array}{l} T_{b, t}^{P, o u t} = T_{t}^{A} + (T_{b, t}^{P, o u t *} - T_{t}^{A}) e x p (- \frac{c_{b}^{h} Δ t}{A_{b} ρ C_{p}} (γ_{b, t} + \frac{1}{2} + \frac{S_{b, t} - R_{b, t}}{m_{b, t - γ_{b, t}}^{P} Δ t})) \\ \forall b \in ℐ^{p i p e}, t \in 𝒯 \end{array}$

(17)

where $T_{t}^{A}$ is the ambient temperature in period t; and $c_{b}^{h}$ is the heat transfer coefficient of pipeline b.

C. SE Formulation of CHPS

SE formulation of CHPS consists of the EPS SE problems, DHN SE problems and the coupling constraints between EPSs and DHNs. Therefore, SE formulation of CHPS can be expressed as:

$m i n \{\sum_{i \in 𝒮^{E P S}} J_{i}^{E P S} (x_{i}^{E P S}) + \sum_{j \in 𝒮^{D H N}} J_{j}^{D H N} (x_{j}^{D H N})\}$

(18)

s.t.

$g_{i}^{E P S} (x_{i}^{E P S}) = 0 \forall i \in 𝒮^{E P S}$

(19)

$g_{j}^{D H N} (x_{j}^{D H N}) = 0 \forall j \in 𝒮^{D H N}$

(20)

$\{\begin{array}{l} E_{i} x_{i}^{E P S} = \sum_{j \in 𝒩_{i}} E D_{i j} \cdot x_{j}^{D H N} \forall i \in 𝒮^{E P S} \\ D_{j} x_{j}^{D H N} = \sum_{i \in 𝒩_{j}} D E_{j i} \cdot x_{i}^{E P S} \forall j \in 𝒮^{D H N} \end{array}$

(21)

where $x^{E P S}$ includes $V, θ, P^{i n j}, a n d Q^{i n j}$ ; $x^{D H N}$ includes $p s, m,$ $Φ, a n d T$ ; $E_{i}, E D_{i j}, a n d D E_{j i}$ are the incidence matrices, e.g., $E D_{i j}$ reflects the relation between $x_{i}^{E P S}$ and $x_{j}^{D H N}$ ; $𝒩_{i}$ , $𝒩_{j}$ are the index sets of the neighboring areas of the areas i and j, respectively; and $𝒮^{D H N},$ $𝒮^{E P S}$ are the index sets of DHN and EPS areas, respectively. Equation (19) refers to the EPS constraints (2) and (3), while (20) refers to the DHN constraints (5)-(17). In CHPS, the energy transfer is realized by the coupling components. The CHP unit is considered as the coupling component in this paper. The constraint (21) reflects the relationships of energy conversion on CHP units. If the $m^{t h}$ CHP unit in the $i^{t h}$ EPS area corresponds to the $n^{t h}$ heat source in the $j^{t h}$ DHN area, the coupling constraints can be expressed as:

$c^{m} P_{m}^{C H P} = Φ_{n}^{C H P}$

(22)

where $c^{m}$ is the heat-to-power ratio of the CHP unit.

III. Decentralized SE Using R-ADMM

The conventional ADMM algorithm has slow convergence rate, and it needs much time to converge when the communication failures occur. To address this problem, the R-ADMM algorithm is adopted in this paper. It has faster convergence rate than ADMM, and the computation efficiency is improved under the situations of communication failures.

A. R-ADMM Without Package Loss

The R-ADMM algorithm is derived using the relaxed Peaceman-Rachford splitting (R-PRS) to the Lagrangian dual problem [

23]. The detailed description of the original R-PRS is given in Appendix A.

For the two-area optimization problem, the formulation can be expressed as:

$\{\begin{array}{l} \underset{x, y}{m i n} \{f (x) + g (y)\} \\ \begin{matrix} s . t . & A x + B y = c \end{matrix} \end{array}$

(23)

By applying the R-PRS to the Lagrangian dual of (23), the variables can be computed as [

23]:

$y_{k} = \underset{y}{a r g m i n} \{g (y) - z_{k}^{T} (B y - c) + \frac{β}{2} {‖B y - c‖}^{2}\}$

(24)

$ψ_{k} = z_{k} - β (B y - c)$

(25)

$x_{k} = \underset{x}{a r g m i n} \{f (x) - (2 ψ_{k} - z_{k})^{T} (A x) + \frac{β}{2} {‖A x‖}^{2}\}$

(26)

$ξ_{k} = 2 ψ_{k} - z_{k} - β (A x_{k})$

(27)

$z_{k + 1} = z_{k} + 2 α (ξ_{k} - ψ_{k})$

(28)

where $y_{k}$ , $x_{k}$ are the solutions of the primal problem in the k^th iteration; $ψ_{k}$ , $ξ_{k}$ are the auxiliary variables which represent the solutions of the Lagrange multiplier in the k^th iteration; $z_{k}$ is the auxiliary variable and the iteration stops when $z_{k}$ reaches its limit value; $β$ is the penalty parameter; and $α \in (0,1]$ is the relaxation parameter.

For (28), the special cases with $α = 0.5$ and $α = 1$ are called the Douglas-Rachford splitting (DRS) and PRS algorithms, respectively. When $α = 0.5$ , i.e., the DRS is applied to the Lagrangian dual of the primal problem, the R-ADMM algorithm is recovered from the conventional ADMM algorithm [

23].

However, the SE formulation (18)-(21) is a general multi-area optimization problem rather than a two-area one. The R-ADMM algorithm cannot be applied in the SE problem directly. We use the method in [

20] to address this problem.

Define $x = [x_{1}^{E P S}, x_{2}^{E P S}, \dots, x_{M}^{E P S}, x_{1}^{D H N}, x_{2}^{D H N}, \dots, x_{N}^{D H N}]$ , where $M$ is the number of EPS areas and $N$ is the number of DHN areas. Define $Ω$ as the feasible region of CHPS, where $x$ satisfy the constraints (19)-(21). Introduce the bridge variable $y = [y_{1}^{E P S}, y_{2}^{E P S}, \dots, y_{M}^{E P S}, y_{1}^{D H N}, y_{2}^{D H N}, \dots, y_{N}^{D H N}]$ , where $y$ is the boundary information of each area, i.e., $y_{i}^{E P S} = E_{i} x_{i}^{E P S},$ $y_{j}^{D H N} = D_{j} x_{j}^{D H N}$ .

Therefore, SE formulation can be reformulated as:

$\{\begin{array}{l} \underset{x \in Ω}{m i n} J (x) \\ \begin{array}{l} s . t . & C x + y = 0 \\ y = P y \end{array} \end{array}$

(29)

where $C$ is the incidence matrix reflecting the relation between the boundary information $y$ and primal variable $x$ ; and $P$ is the permutation matrix, which swaps the EPS boundary information with the DHN boundary information in $y$ according to (21).

Introduce the indicator function $l_{I - P} (y)$ which is equal to 0 if $y = P y$ and $+ \infty$ otherwise. Then the SE formulation is rewritten as:

$\{\begin{array}{l} \underset{x \in Ω, y}{m i n} \{J (x) + l_{I - P} (y)\} \\ \begin{matrix} s . t . & C x + y = 0 \end{matrix} \end{array}$

(30)

The formulation (30) is similar to (23). By applying the (24)-(28) to the Lagrange dual of (30), the variable $x$ can be updated iteratively in the following way:

$x_{k} = \underset{x}{a r g m i n} \{J (x) + (P z_{k})^{T} C x + \frac{β}{2} {‖C x‖}^{2}\}$

(31)

$z_{k + 1} = (1 - α) z_{k} - α (P z_{k}) - 2 α β (C x_{k})$

(32)

Equations(31) and (32) can be partitioned into several subproblems for each area. For the $i^{t h}$ EPS area subproblem, we can obtain:

$\underset{x_{i}^{E P S}}{m i n} \{J_{i}^{E P S} x_{i}^{E P S} - \sum_{j \in 𝒩_{i}} (z_{j i, k}^{E P S})^{T} E_{i} x_{i}^{E P S} + \frac{β}{2} {‖E_{i} x_{i}^{E P S}‖}^{2} : g_{i}^{E P S} (x_{i}^{E P S}) = 0\}$

(33)

$\begin{array}{l} z_{j i, k + 1}^{E P S} = (1 - α) z_{j i, k}^{E P S} - α z_{j i, k}^{D H N} + 2 α β \cdot E D_{i j} \cdot x_{j}^{D H N} = \\ (1 - α) z_{j i, k}^{E P S} + α U_{j i, k}^{D H N} \end{array}$

(34)

$U_{j i, k}^{D H N} = - α z_{j i, k}^{D H N} + 2 α β \cdot E D_{i j} \cdot x_{j}^{D H N}$

(35)

where $z_{j i, k}^{E P S}$ is the auxiliary variable in the objective function of the $i^{t h}$ EPS area subproblem, which is updated based on the auxiliary variable $z_{j i, k}^{D H N}$ ; and $U_{j i, k}^{D H N}$ is the boundary information sending from the $j^{t h}$ DHN area to the $i^{t h}$ EPS area.

For the $j^{t h}$ DHN subproblem, we can obtain:

$\underset{x_{j}^{D H N}}{m i n} \{J_{j}^{D H N} x_{j}^{D H N} - \sum_{i \in 𝒩_{j}} (z_{i j, k}^{D H N})^{T} D_{j} x_{j}^{D H N} + \frac{β}{2} {‖D_{j} x_{j}^{D H N}‖}^{2} : g_{j}^{D H N} (x_{j}^{D H N}) = 0\}$

(36)

$z_{i j, k + 1}^{D H N} = (1 - α) z_{i j, k}^{D H N} - α z_{i j, k}^{E P S} + 2 α β \cdot D E_{j i} \cdot x_{i}^{E P S} = (1 - α) z_{i j, k}^{D H N} + α U_{j i, k}^{E P S}$

(37)

$U_{i j, k}^{E P S} = - α z_{i j, k}^{E P S} + 2 α β \cdot D E_{j i} \cdot x_{i}^{E P S}$

(38)

where $U_{i j, k}^{E P S}$ is the boundary information sending from the $i^{t h}$ EPS area to the $j^{t h}$ DHN area.

The communication procedure between the $i^{t h}$ EPS area and the $j^{t h}$ DHN area is presented in Fig. 3. During the iteration, each area updates its local auxiliary variables $z$ after receiving the boundary information $U$ from its neighboring areas.

Fig. 3 Communication procedure between neighboring areas.

The algorithm stops when the primal residual $r_{k}$ and the dual residual $s_{k}$ are smaller than their corresponding tolerances:

$\{\begin{array}{l} {‖r_{k}‖}_{\infty} \leq ε^{r e s} \\ {‖s_{k}‖}_{\infty} \leq ε^{d u a l} \end{array}$

(39)

where ${‖r_{k}‖}_{\infty}$ , ${‖s_{k}‖}_{\infty}$ are the primal residual and dual residual during the k^th R-ADMM iteration, respectively; and $ε^{r e s}$ , $ε^{d u a l}$ are the feasibility tolerances for the primal residual and dual residual, respectively. More information about these residuals and tolerances can be found in [

24].

For the convex optimization problem, the optimal solution of the Lagrangian dual problem is equal to that of the primal problem. With $α \in (0,1)$ and $β > 0$ , the proposed algorithm can converge to its global optimal solution [

20]. The special case

$α = 1$ requires more restrictive assumptions to ensure the convergence [25]. However, SE formulation is a non-convex optimization problem. The algorithm cannot guarantee to converge to its global optimal solution [26].

B. R-ADMM Considering Communication Packet Loss

Considering the instability or potential malfunction of the communication channels, one area may not receive the boundary information from its neighboring areas by chance. This situation is considered as the communication packet loss, and the probability of packet loss occurs as the packet loss probability.

Assume that the packet loss probability from area $i$ to area $j$ is equal to $p$ , we can obtain：

$\{\begin{array}{l} P (L_{i j} = 1) = p \\ P (L_{i j} = 0) = 1 - p \end{array}$

(40)

where $L_{i j} = 1$ if the packet transmitted by the area $i$ to the area $j$ is lost, and $L_{i j} = 0$ otherwise.

For the auxiliary variable $z_{i j}$ , if the area $j$ receives the boundary information $U_{i j}$ from the area $i$ , the variable $z_{i j}$ is updated using the received boundary information. Otherwise, its value is not changed. The detailed expression is:

$z_{i j, k + 1} = L_{i j, k} z_{i j, k} + (1 - L_{i j, k}) [(1 - α) z_{i j, k} + α U_{i j, k}]$

(41)

The convergence performance of the algorithm working in the lossy scenarios can be found in [

20]. For the convex optimization problem, the R-ADMM with

$α \in (0,1)$ and

$β > 0$ can still converge in the lossy scenarios.

C. R-ADMM-based Procedure of SE for CHPS

In the process of solving SE for CHPS, it is noteworthy that the heat dynamic constraints (11)-(17) cannot be handled readily. We use the alternating estimation method in [

8] to address the DHN subproblems. The procedure for solutions includes two stages. In the first stage, the complicated variables (CVs) including

$K_{b, t, k}^{T D}, γ_{b, t}, a n d ϕ_{b, t}$ are updated according to the measurement information. In the second stage, the noncomplicated variables (NCV) including

$p s_{i}, m_{b, t}^{P}, T_{i, t}^{N D},$

$ϕ_{i, t}^{L D}, a n d ϕ_{i, t}^{C H P}$ are updated based on the CVs and measurement information. The measurement information is updated with the SE results, and the two-stage problems are solved in an alternating way until the convergence criterion is met.

By combining this alternating strategy and the R-ADMM algorithm, the calculation process of SE for CHPS is presented in Fig. 4, and its detailed steps are described as follows.

Fig. 4 Procedure of R-ADMM of SE for CHPS.

Step 1: initialize the state variables. Set the index of iteration $k = 1$ .

Step 2: solve the SE problem for each EPS area (33) and each DHN area (36) using the alternating estimation strategy described above.

Step 3: calculate the boundary information according to (35) or (38) for each subsystem and transmit it to the neighboring area.

Step 4: update the auxiliary variables based on the communication situation using (41) for each subsystem.

Step 5: if the termination criterion (39) is met, stop the algorithm. Otherwise, set the number of iterations $k$ as $k + 1$ and return to the Step 2.

IV. Case Studies

Two numerical simulations are conducted to analyze the performance of the proposed algorithm. Case 1 is a CHPS with a 9-bus EPS and a 6-node DHN. EPS is connected to DHN through a CHP unit at bus 3. Case 2 is a CHPS with a 118-bus EPS, an 8-node DHN and a 45-node DHN. EPS is connected to DHNs through the CHP units at bus 12 and bus 54. The detailed information on the above cases can be found in [

9], [27].

The measurement errors are simulated as independent zero-mean Gaussian distribution with the standard deviation of 1% of the measurement ranges of each metering device. The weight coefficients used in the objective function are calculated based on the standard deviations of the corresponding measurements. The probability distribution of packet losses is according to (40). The parameters $α = 1 a n d β = 30$ are set for the R-ADMM algorithm by the default. The variables are initialized using their measurement data. According to [

20], the auxiliary variables are initialized as:

$\{\begin{array}{l} z_{j i, t, 0}^{E P S} = β \sum_{j \in 𝒩_{i}} E D_{i j} \cdot x_{j, t, 0}^{D H N} \\ z_{i j, t, 0}^{D H N} = β \sum_{i \in 𝒩_{j}} D E_{j i} \cdot x_{i, t, 0}^{E P S} \end{array}$

(42)

All tests are performed on a computer with four processors running at 1.8 GHz and 16 GB of RAM. All the optimization problems are solved with the commercial solver IPOPT 3.12.9 [

28] using MATLAB R2018a.

A. Scenarios with Low Packet Loss Probability

In this sub-section, we study the performance of the proposed algorithm with the low packet loss probability. For the reliable communication, the packet loss probability is low. The packet loss probability is set as 0.05 for each area [

29] and a SE problem with 12 periods is considered. To simulate the random measurement errors and the random packet loss among the neighboring areas, 500 Monte Carlo experiments are conducted for each period.

1)　Accuracy of SE

We analyze the accuracy of the proposed algorithm with a performance metric $S^{H} / S^{M}$ [

30], where

$S^{H}$ is the normalized average estimation error and

$S^{H}$ is the normalized average measurement error:

$S^{H} = \frac{1}{N} \sum_{n = 1}^{N} \sqrt[]{\frac{1}{M} \sum_{i = 1}^{M} \frac{(f_{i}^{e s t} - f_{i}^{t r u e})^{2}}{W_{i}}} \times 100 %$

(43)

$S^{M} = \frac{1}{N} \sum_{n = 1}^{N} \sqrt[]{\frac{1}{M} \sum_{i = 1}^{M} \frac{(f_{i}^{m e a s} - f_{i}^{t r u e})^{2}}{W_{i}}} \times 100 %$

(44)

where $f_{i}^{e s t}$ is the estimation result; $f_{i}^{m e a s}$ is the measurement data; $f_{i}^{t r u e}$ is the actual operation data; and $W_{i}$ is the weight coefficient.

A smaller value of $S^{H} / S^{M}$ indicates a greater improvement from the measurement states to the estimation states, which implies higher quality of SE result. $S^{H} / S^{M}$ values of the R-ADMM results for Case 1 and Case 2 are presented in Figs.5 and 6, respectively.

Fig. 5 $S^{H} / S^{M}$ values for R-ADMM results in Case 1.

Fig. 6 $S^{H} / S^{M}$ values for R-ADMM results in Case 2.

For the whole CHPS, the $S^{H} / S^{M}$ values of Case 1 ranges from 43.66% to 48.91%, and the $S^{H} / S^{M}$ values of Case 2 ranges from 50.58% to 56.48%. These results show that the R-ADMM can greatly improve the accuracy of the measurement data.

To further illustrate the effectiveness of the proposed algorithm, we compare the SE solutions between the R-ADMM and the centralized SE (CSE) algorithm. 500 SE tests are carried out using these algorithms on the same measurement data for each period. The mean values of the relative errors (MREs) of the R-ADMM algorithm for Case 1 and Case 2 are presented in Tables I and II, respectively. As shown in these tables, the order of MREs of the R-ADMM algorithm is less than 1×10^⁃5. This shows the rationality of the proposed algorithm.

TABLE I MREs of R-ADMM Algorithm in Case 1

Time (hour)	MRE (%)
1	2.67×10^⁃7
2	1.39×10^⁃6
3	3.21×10^⁃7
4	1.95×10^⁃7
5	6.84×10^⁃7
6	2.71×10^⁃6
7	7.80×10^⁃6
8	2.45×10^⁃7
9	6.27×10^⁃7
10	1.15×10^⁃7
11	2.31×10^⁃6
12	2.84×10^⁃7

TABLE II MREs of R-ADMM Algorithm in Case 2

Time (hour)	MRE (%)
1	1.41×10^⁃7
2	1.42×10^⁃7
3	7.83×10^⁃8
4	5.45×10^⁃8
5	2.25×10^⁃7
6	3.16×10^⁃7
7	9.31×10^⁃8
8	3.31×10^⁃7
9	2.08×10^⁃7
10	1.45×10^⁃7
11	1.37×10^⁃7
12	7.52×10^⁃8

2)　Comparison of Computation Efficiency

To study the computation efficiency of R-ADMM, an SE test is carried out using the R-ADMM with different relaxation parameters on the same measurement data. The comparisons of the results of different algorithms for Case 1 and Case 2 are listed in Tables III and IV, respectively. The fourth column of these tables shows the MREs of the state variables between the CSE algorithm and the decentralized algorithms.

TABLE III Comparison of Results of Different Algorithms in Case 1

Algorithm	CPU time (s)	$n_{i t e}$	MRE (%)	$S^{H} / S^{M}$ (%)
CSE	1.0163			48.7662
ADMM ( $α = 0.5$ )	9.3704	20	3.1736×10^⁃7	48.7662
R-ADMM ( $α = 0.6$ )	7.7554	17	2.6059×10^⁃7	48.7662
R-ADMM ( $α = 0.8$ )	5.4627	12	1.7661×10^⁃7	48.7662
R-ADMM ( $α = 1$ )	3.0831	7	1.8176×10^⁃7	48.7662

TABLE IV Comparison of Results of Different Algorithms in Case 2

Algorithm	CPU time (s)	$n_{i t e}$	MRE (%)	$S^{H} / S^{M}$ (%)
CSE	7.1014			50.891
ADMM ( $α = 0.5$ )	21.2412	19	2.2782×10^⁃8	50.891
R-ADMM ( $α = 0.6$ )	17.1962	16	1.7467×10^⁃8	50.891
R-ADMM ( $α = 0.8$ )	15.4243	14	1.6605×10^⁃8	50.891
R-ADMM ( $α = 1$ )	11.0602	10	1.5028×10^⁃8	50.891

SE problem is a non-convex optimization problem. Different values of relaxation parameter for R-ADMM result in differences in the auxiliary variables of their objective functions in each iteration, which may lead to different local optima [

31]. As shown in the fourth column of Tables III and IV, the differences of the SE solutions are minor, which shows that R-ADMM with different values of relaxation parameter converges to the same solution in the simulation tests.

Because R-ADMM needs iterations to reach the termination criterion, the computation time of R-ADMM is longer than that of the CSE algorithm. The advantage of the R-ADMM over the centralized algorithm mainly depends on its practicality. In the actual production process, different energy systems are operated by different energy system operators independently. The measurement information among different energy systems cannot be shared fully due to the privacy protection. Therefore, it is impractical to gather all measurement information in a centralized manner. While the R-ADMM can guarantee the independent operation and the coupled results in a decentralized manner in this situation.

The evolution of primal residuals and dual residuals in Cases 1 and 2 are shown in Figs.7 and 8, respectively.

Fig. 7 Evolution of residuals with different $α$ in Case 1. (a) Primal residual. (b) Dual residual.

Fig. 8 The evolution of residuals with different $α$ in Case 2. (a) Primal residual. (b) Dual residual.

When $α \in (0,1]$ [

23] shows that the primal residuals and the dual residuals hold the faster convergence rates when

$α$ increases. As shown in Figs.7 and 8, the relaxation parameter larger than 0.5 corresponds to a faster convergence rate on both primal residuals and dual residuals. For the benefit of the fast convergence rates, the algorithm needs fewer iterations and shorter computation time as the relaxation parameter increases. Note that setting

$α = 0.5$ yields the conventional ADMM. It is obvious that the R-ADMM can improve the computation efficiency, which reflects its advantages compared with the ADMM.

B. Scenarios with High Packet Loss Probability

To study the performance of the R-ADMM with high packet loss probability, we compare three scenarios with different packet loss probabilities for each case. The parameters $α = 1 a n d β = 30$ are set for the R-ADMM algorithm. The measurement data are set the same for each scenario. The 500 Monte Carlo experiments are conducted to simulate the random packet loss among the neighboring areas for each scenario. The comparison of different algorithms in different scenarios for Cases 1 and 2 are listed in Tables V and VI, respectively. The third column of these tables shows the MREs of the state variables between the scenario of the high packet loss probability ( $p = 0.50$ or $p = 0.95$ ) with that of the low packet loss probability ( $p = 0.05$ ) for the same algorithm.

TABLE V Comparison of Results of Different Algorithms in Different Scenarios in Case 1

Scenario	Algorithm	MRE (%)	CPU time (s)	$S^{H} / S^{M}$ (%)
$p = 0.05$	ADMM		21.2412	50.8910
$p = 0.05$	R-ADMM		11.0602	50.8910
$p = 0.5$ 0	ADMM	1.9740×10^⁃6	43.9914	50.8910
$p = 0.5$ 0	R-ADMM	1.2707×10^⁃6	23.2012	50.8910
$p = 0.95$	ADMM	2.8567×10^⁃5	269.0110	50.8912
$p = 0.95$	R-ADMM	1.5013×10^⁃5	102.8172	50.8912

TABLE VI Comparison of Results of Different Algorithms in Different Scenarios in Case 2

Scenario	Algorithm	MRE (%)	CPU time (s)	$S^{H} / S^{M}$ (%)
$p = 0.05$	ADMM		9.3704	48.7662
$p = 0.05$	R-ADMM		3.0831	48.7662
$p = 0.5$ 0	ADMM	3.2684×10^⁃7	19.0363	48.7662
$p = 0.5$ 0	R-ADMM	1.9231×10^⁃7	10.0696	48.7662
$p = 0.95$	ADMM	3.4064×10^⁃7	234.0327	48.7662
$p = 0.95$	R-ADMM	1.9232×10^⁃7	77.1191	48.7662

As shown in the fourth and fifth columns in Tables V and VI, a high packet loss probability will lead to a minor change of the SE result and a long computation time. When the communication packet is lost, the auxiliary variables used in the objective function during the next iteration cannot be updated. These “out-of-date” auxiliary variables will affect the accuracy of SE and reduce the computation efficiency.

As shown in the third column in Tables V and VI, the changes of the state variables among the scenarios of different packet loss probabilities with R-ADMM are smaller than those of the ADMM. As shown in the fourth column of Tables V and VI, the increased computation time under the high packet loss probability with R-ADMM are much smaller than that with ADMM. The influences of the high packet loss probability on R-ADMM is smaller than those on ADMM, which shows greater robustness of the proposed algorithm.

C. Scenario with Different Packet Loss Probabilities for Different Areas

The packet loss probabilities for different sub-systems are set to different values in Case 2 (118-bus EPS (EPS1)), as shown in Table VII. $p_{A \to B}$ represents the packet loss probability sending from area $A$ to area $B$ .

TABLE VII Different Packet Loss Probabilities

8-node DHN (DHN1)	45-node DHN (DHN2)
$p_{E P S 1 \to D H N 1} = 0.1$	$p_{E P S 1 \to D H N 2} = 0.1$
$p_{D H N 1 \to E P S 1} = 0.5$	$p_{D H N 2 \to E P S 1} = 0.3$

In the simulations, R-ADMM meets the convergence criterion with a $S^{H} / S^{M}$ value of 50.3136% in 11.8584 s, the ADMM meets the convergence criterion with a $S^{H} / S^{M}$ value of 50.3163% in 35.9279 s. The evolution of the primal residuals and the dual residuals for ADMM and R-ADMM are depicted in Figs.9 and 10, respectively. R-ADMM converges after 16 iterations while ADMM converges after 36 iterations. Due to the faster convergence rate, R-ADMM maintains greater computation efficiency compared with ADMM.

Fig. 9 Evolution of residuals with iterations of ADMM.

Fig. 10 Evolution of residuals with iterations of R-ADMM.

Ⅴ. Conclusion

This paper proposes a decentralized solution for SE of CHPS. An R-ADMM algorithm is developed based on the boundary information communication among the neighboring areas. In addition, the robust implementation of the R-ADMM algorithm considering the communication failures among different sub-systems is proposed. Case studies show that the proposed algorithm can yield reasonable SE results with great computation efficiency and robustness. In future work, the additional functionalities will be incorporated into the proposed algorithm such as observability analysis and bad data detection and identification. In addition, more differences among various energy measurement systems will be considered such as different sample rates and different sample instant.

References

D. Wang, L. Liu, H. Jia et al., “Review of key problems related to integrated energy distribution systems,” CSEE Journal of Power and Energy Systems, vol. 4, no. 2, pp. 130-145, Jun. 2018. [百度学术]

X. Liu, J. Wu, N. Jenkins et al., “Combined analysis of electricity and heat networks,” Applied Energy, vol. 162, pp. 1238-1250, Jan. 2015. [百度学术]

J. Wu, J. Yan, and H. Jia, “Integrated energy systems,” Applied Energy, vol. 167, pp. 155-157, Apr. 2016. [百度学术]

T. Fang and R. Lahdelma, “State estimation of district heating network based on customer measurements,” Applied Thermal Engineering, vol. 73, no. 1, pp. 1211-1221, Dec. 2014. [百度学术]

A. Abur and A. Gomez-Exposito. Power System State Estimation Theory and Implementation. Perak Darul Ridzuan:Universiti Teknologi Petronas, 2004. [百度学术]

J. Dong, H. Sun, Q. Guo et al., “State estimation for combined electricity and heat networks,” Power System Technology, vol. 40, no. 6, pp. 1635-1641, Jun. 2016. [百度学术]

T. Sheng, Q. Guo, H. Sun et al., “Two-stage state estimation approach for combined heat and electric networks considering the dynamic property of pipelines,” Energy Procedia, vol. 142, pp. 3014-3019, Dec. 2017. [百度学术]

T. Zhang, Z. Li, Q. H. Wu et al., “Dynamic state estimation of combined heat and power system considering quasi-dynamics of temperature in pipelines,” in Proceedings 2018 International Conference on Power System Technology, Guangzhou, China, Nov. 2018, pp. 232-237. [百度学术]

T. Zhang, Z. Li, Q. H. Wu et al., “Decentralized state estimation of combined heat and power systems using the asynchronous alternating direction method of multipliers,” Applied Energy, vol. 248, pp. 600-613, Apr. 2019. [百度学术]

Y. Du, W. Zhang, and T. Zhang, “ADMM-based distributed state estimation for integrated energy system,” CSEE Journal of Power and Energy Systems, vol. 5, no. 2, pp. 275-283, Jun. 2019. [百度学术]

H. Zang, M. Geng, M. Xue et al., “A robust state estimator for integrated electrical and heating networks,” IEEE Access, vol. 7, pp. 109990-110001, Aug. 2019. [百度学术]

D. K. Molzahnl, F. Dörfler, H. Sandberg, et al., “A survey of distributed optimization and control algorithms for electric power systems,” IEEE Transaction on Smart Grid, vol. 8, no. 6, pp. 2941-2962, Nov. 2017. [百度学术]

J. Huang, Z. Li, and Q. H. Wu, “Coordinated dispatch of electric power and district heating networks: a decentralized solution using optimality condition decomposition,” Applied Energy, vol. 206, pp. 1508-1522, Sept. 2017. [百度学术]

V. Kekatos and G. B. Giannakis, “Distributed robust power system state estimation,” IEEE Transactions on Power Systems, vol. 28, no. 2, pp. 1617-1626, May. 2013. [百度学术]

W. Zheng, W. Wu, and A. Gomez-Exposito, “Distributed robust bilinear state estimation for power systems with nonlinear measurements,” IEEE Transactions on Power Systems, vol. 32, no. 1, pp. 499-509, Jan. 2017. [百度学术]

J. Stamp, A. McIntyre, and B. Ricardson, “Reliability impacts from cyber attack on electric power systems,” in Proceedings of 2009 IEEE PES Power Systems Conference and Exposition, Seattle, USA, Apr. 2009, pp. 1-8. [百度学术]

G. N. Anandini and I. Gupta, “Estimating the lost real-time measurements under communication failure for distribution system state estimation,” in Proceedings of 2016 National Power Systems Conference, Bhubaneswar, India, Jan. 2016, pp. 1-5. [百度学术]

C. Gu and P. Jirutitijaroen, “Dynamic state estimation under communication failure using Kriging based bus load forecasting,” IEEE Transactions on Power Systems, vol. 30, no. 6, pp. 2831-2840, Nov. 2015. [百度学术]

M. M. Rana and A. Shahirinia, “Distributed dynamic state estimation considering packet losses in interconnected smart grid subsystems: linear matrix inequality approach,” IEEE Access, vol. 8, pp. 2687-2693, Jan. 2020. [百度学术]

N. Bastianello, M. Todescato, R. Carli et al., “Distributed optimization over lossy networks via relaxed Peaceman-Rachford splitting: a robust ADMM approach,” in Proceedings of 2018 European Control Conference, Limassol, Cyprus, Sept. 2018, pp. 477-482. [百度学术]

A. Benonysson, B. Bøhm, and H. F. Ravn, “Operational optimization in a district heating system,” Energy Conversion and Management, vol. 36, no. 5, pp. 297-314, May 1995. [百度学术]

Z. Li, W. Wu, J. Wang et al., “Transmission-constrained unit commitment considering combined electricity and district heating networks,” IEEE Transactions on Sustainable Energy, vol. 7, no. 2, pp. 480-492, Apr. 2016. [百度学术]

D. Davis and W. Yin, “Faster convergence rates of relaxed Peaceman-Rachford and ADMM under regularity assumptions,” Mathematics of Operations Research, vol. 42, no. 3, pp. 783-805, Feb. 2017. [百度学术]

S. Boyd, N. Parikh, E. Chu et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, Jan. 2011, pp. 1-122. [百度学术]

P. Giselsson and S. Boyd, “Linear convergence and metric selection for Douglas-Rachford splitting and ADMM,” IEEE Transactions on Automatic Control, vol. 62, no. 2, pp. 532-544, Feb. 2017. [百度学术]

K. Guo, D. Han, and T. Wu, “Convergence of alternating direction method for minimizing sum of two nonconvex functions with linear constraints,” International Journal of Computer Mathematics, vol. 94, no. 8, Sept. 2016, pp. 1653-1669. [百度学术]

Z. Li. (2015, Aug.). Test data of 6-bus system for UC-CEHN. [Online]. Available: http://motor.ece.iit.edu/data/UCCEHN_6bus.xls [百度学术]

Modern Investment Technologies, Ltd.. (2008, May). IPOPT homepage. [Online]. Available: https://projects.coin-or.org/Ipopt [百度学术]

M. Rana, L. Li, and S. Wu, “Distributed state estimation over unreliable communication networks with an application to smart grids,” IEEE Transactions on Green Communications and Networking, vol. 1, no. 1, Mar. 2017, pp. 89-96. [百度学术]

E. Yu, Power System State Estimation. Beijing: Water Resources and Electric Power Press, 1985. [百度学术]

G. Li, T. Liu, and T. Pong, “Peaceman-Rachford splitting for a class of nonconvex optimization problems,” Computational Optimization and Applications, vol. 68, no. 2, Nov. 2017, pp. 407-436. [百度学术]

H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces. New York: Springer, 2011, vol. 408. [百度学术]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher