Optimal Cluster Scheduling of Active/Reactive Power for Distribution Network Considering Aggerated Flexibility of Heterogeneous Building-Integrated DERs

1. Introduction

Construction sector is the third largest energy consumer, accounting for 33% of the global energy demand and 16% of global CO₂ emissions in 2022 [1]. Considering the growing tendency of urbanization progress and the continuous development of the power industry, the proportion of construction energy consumption will be increased by 40% in 2030 and even 50% in 2050 [2]. The aggregation of flexible DERs in the buildings is recognized as a promising alternative method to reduce CO₂ emissions and energy consumption, promoting the profound integration of distributed energy resources (DERs) and building energy supply [3]. The smart building system equipped with power access points will achieve the flexible resources aggregation and exploit the regulation capacity of these resources sufficiently, intensifying the interdependency between the buildings and electricity network [3,4]. Nevertheless, the location of building-integrated photovoltaic (BIPV), building-integrated electric vehicle (BIEV), building-integrated energy storage (BIES), building-integrated temperature-controlled load (BITC) and other comprehensive energy production and storage equipment in the smart building system is generally scattered, resulting in the high complexity of centralized optimization scheduling due to the lack of information exchange and energy interaction [4,5,6]. Consequently, the optimization scheduling of distribution network with cluster division for flexible resources integrated into the smart building system will be ideal for the high-penetrated deployment of DERs to meet the carbon-neutral targets.

One typical method studied in [7] for the cluster optimization scheduling is to use the second-order neighborhood similarity matrix and modular matrix to extract relatively more comprehensive complex network feature information for network division. A chromosome encoding method considering modular index based on the electrical designed to satisfy the demand of distribution network programmer [8]. A K-means clustering algorithm based on the variation sensitivity of voltage-to-power between nodes to solve the voltage overlimit and voltage fluctuations [9]. In addition, in order to tackle the difficulty in the uncertainty and centralized control, a multi-time rolling dynamic cluster division method with the equilibrium degree of reactive power for active distribution network [10]. However, the cluster division in these studies only involve a single index, without considering the comprehensive indexes including the flexibility regulation capacity and electrical connection.

Integration and quantification for the regulation capacity of flexible resources is regarded as an effective means to take advantage of the flexibility of the smart building system [4,11,12,13,14,15]. Lots of inspiring works have been reported on the virtual battery model [11], virtual synchronous machine model [12], flexible supply and demand balance model [13], and node power model [14] for the aggregation of flexible resources. A generic and scalable approach considering zonotopic sets for flexible energy systems is studied in [15] to quantitatively describe their flexibility. A constraint space superposition method based on Minkowski Sum is designed to seek the maximal inner cube constraints for the virtual synchronous machine of the heat pump cluster [16]. Moreover, a day-ahead optimal scheduling model based on the generalized virtual battery model is developed to determine the energy schedule and reserve capacity of electric vehicle and air conditioning [17]. After the flexible resources aggregation, the coordination and optimization of distribution network with cluster division based on the flexible resources in the smart building system can be achieved. Nevertheless, the potential active and reactive power support capacitates from the clusters of the building-integrated flexible DERs for distribution networks are not involved.

In this paper, the customized virtual battery models are proposed for building-integrated flexible DERs to aggregate the power characteristics of flexible resources while quantifying the regulation capacity of flexible resources with time-shifting power and energy boundaries. Moreover, a cluster division algorithm considering the structural and functional cluster division indexes for building-integrated flexible DERs is developed to automatically divide the distribution network with high proportion renewables and flexible DERs into the optimal clusters. Furthermore, an optimal active-reactive power collaborative scheduling model for the clusters of the building-integrated flexible DERs is formulated to minimize the operational cost, the power loss cost, and penalty cost for flexibility deficiency while handling the uncertainties from the renewable energy in the distribution network.

2. Flexible regulation capacity of flexible DERs in smart building system

2.1. Smart building system

Figure 1 depicts the typical modality of smart building system among the high renewable penetrated distribution network. The smart building system is generally equipped with cogeneration systems, distributed photovoltaics, electric vehicle, heat pumps, energy storage to satisfy the demand of electricity, cold and heat load in the building [5]. In order to fully utilize the space and time flexibility of the building-integrated flexible DERs, the establishment of virtual battery model for flexible DERs is necessary, which can quantify and integrate the flexibility and regulation capacity of these resources.

2.2. Battery model of ES in the smart building system

BIES can transfer the electricity generated by BIPVs to balance the energy of the power system, improving the reliability of system operation and power supply quality. A battery model of ES in the smart building system (1), consisting of charge power constraints, discharge power constraints, and energy state constraints is developed to mathematically describe the operating characteristic,

$$P_{\mathrm{ch},i,t,\min }^{\mathrm{ES}}\le P_{\mathrm{ch},i,t}^{\mathrm{ES}}\le P_{\mathrm{ch},i,t,\max }^{\mathrm{ES}}$$

(1a)

$$P_{\mathrm{dis},i,t,\min }^{\mathrm{ES}}\le P_{\mathrm{dis},i,t}^{\mathrm{ES}}\le P_{\mathrm{dis},i,t,\max }^{\mathrm{ES}}$$

(1b)

$$E_{i,t,\min }^{\mathrm{ES}}\le E_{i,t}^{\mathrm{ES}}\le E_{i,t,\max }^{\mathrm{ES}}$$

(1c)

$$E_{i,t+1}^{\mathrm{ES}}=E_{i,t}^{\mathrm{ES}}+P_{i,t}^{\mathrm{ES}}\Delta t=E_{i,t}^{\mathrm{ES}}+\left({\eta }_{\mathrm{ch}}^{\mathrm{ES}}{P}_{\mathrm{ch},i,t}^{\mathrm{ES}}-\frac{P_{\mathrm{dis},i,t}^{\mathrm{ES}}}{{\eta }_{\mathrm{dis}}^{\mathrm{ES}}}\right)\Delta t$$

(1d)

where $P_{\mathrm{c}\mathrm{h},i,t}^{\mathrm{E}\mathrm{S}}$and $P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t}^{\mathrm{E}\mathrm{S}}$ are the charging and discharging power of BIES connected at node $i$ and time $t$; $P_{\mathrm{c}\mathrm{h},i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{S}}$, $P_{\mathrm{c}\mathrm{h},i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{S}}$, $P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{S}}$ and $P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{S}}$ denote the lower and upper bounds of the charging and discharging power of BIES; $E_{i,t}^{\mathrm{E}\mathrm{S}}$and $P_{i,t}^{\mathrm{E}\mathrm{S}}$ are the operation capacity and actual power of BIES; $E_{i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{S}}$ and $E_{i,t,\mathrm{m}\mathrm{i}\mathrm{m}}^{\mathrm{E}\mathrm{S}}$ are the threshold of the operation capacity of BIES; ${\eta }_{\mathrm{c}\mathrm{h}}^{\mathrm{E}\mathrm{S}}$ and ${\eta }_{\mathrm{d}\mathrm{i}\mathrm{s}}^{\mathrm{E}\mathrm{S}}$ are the charge and discharge efficiency coefficient of BIES, and $∆t$ is the time interval.

The flexible regulation capacity of energy storage in the smart building system with energy time-shift capacity and power regulation capacity can be expressed as,

$$F_{i,t}^{\mathrm{ES},\mathrm{s},\mathrm{up}}=P_{i,t,\max }^{\mathrm{ES}}-P_{i,t}^{\mathrm{ES}}$$

(2a)

$$F_{i,t}^{\mathrm{ES},\mathrm{s},\mathrm{dn}}=P_{i,t}^{\mathrm{ES}}-P_{i,t,\min }^{\mathrm{ES}}$$

(2b)

where $F_{i,t}^{\mathrm{E}\mathrm{S},\mathrm{s},\mathrm{u}\mathrm{p}}$ and $F_{i,t}^{\mathrm{E}\mathrm{S},\mathrm{s},\mathrm{d}\mathrm{n}}$ denote the upgraded flexible energy supply and downgraded flexible energy supply of BIES connected at node $i$ and time $t$, deriving from the lower and upper bounds of the actual power $P_{i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{S}}$ and $P_{i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{S}}$ of BIES.

2.2. Virtual battery (VB) model in the smart building system

The VB model can accurately describe the operation and regulation characteristics of flexible DERs by constraining the energy state, energy boundary and power boundary, which can be expressed as,

$$P_{i,t,\min }^G\le P_{i,t}^G\le P_{i,t,\max }^G$$

(3a)

$$E_{i,t,\min }^G\le E_{i,t}^G\le E_{i,t,\max }^G$$

(3b)

$$E_{i,t+1}^G=E_{i,t}^G+P_{i,t}^G\Delta t+\Delta E_{i,t}^G$$

(3c)

where $P_{i,t}^G$ and $E_{i,t}^G$ are the regulation power and flexibility reserve energy of the flexible DERs in the smart building system connected at node $i$ and time $t$; $P_{i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^G$, $P_{i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^G$, $E_{i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^G$, $E_{i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^G$ denote the lower and upper bounds of the regulation power and flexibility reserve energy; (3c) represents the state of the flexibility reserve energy of the flexible DERs and $∆E_{i,t}^G$ is the effect of other factors on the electric energy of the VB model.

2.3. VB model of EV in the smart building system

In general, the power battery of BIEV can achieve power interaction with the grid through flexibly switching between charging, placing, and discharging states. The grid-connected operation status of BIEV is illustrated in Figure 2, where BIEV with an initial operating capacity of $E_{i,t_0}^{\mathrm{E}\mathrm{V}}$ are connected to the grid with the maximum charge and discharge power of $P_{\mathrm{c}\mathrm{h},i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{V}}$and$P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{V}}$ at time $t_0$, reaching the expected operating capacity of $E_{i,t}^{\mathrm{E}\mathrm{V}}$ when off the grid at time $t_1$. Thus, the virtual battery model of BIEV can be formulated as follows,

$$P_{\mathrm{ch},i,t,\min }^{\mathrm{EV}}\le P_{\mathrm{ch},i,t}^{\mathrm{EV}}\le P_{\mathrm{ch},i,t,\max }^{\mathrm{EV}}$$

(4a)

$$P_{\mathrm{di}s,i,t,\min }^{\mathrm{EV}}\le P_{\mathrm{di}s,i,t}^{\mathrm{EV}}\le P_{\mathrm{di}s,i,t,\max }^{\mathrm{EV}}$$

(4b)

$$E_{i,t_0}^{\mathrm{EV}}-P_{\mathrm{dis},i,t,\max }^{\mathrm{EV}}(t-t_0)\le E_{i,t}^{\mathrm{EV}}\le E_{i,t_0}^{\mathrm{EV}}+P_{\mathrm{ch},i,t,\max }^{\mathrm{EV}}(t-t_0)$$

(4c)

$$E_{i,t+1}^{\mathrm{EV}}=E_{i,t}^{\mathrm{EV}}+P_{i,t}^{\mathrm{EV}}\Delta t=E_{i,t}^{\mathrm{EV}}+\left({\eta }_{\mathrm{ch}}^{\mathrm{EV}}{P}_{\mathrm{ch},i,t}^{\mathrm{EV}}-\frac{P_{\mathrm{dis},i,t}^{\mathrm{EV}}}{{\eta }_{\mathrm{dis}}^{\mathrm{EV}}}\right)\Delta t$$

(4d)

where $P_{\mathrm{c}\mathrm{h},i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{V}}$ and $P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{V}}$ are the minimum threshold the charging and discharging power of BIEV connected at node $i$ and time $t$; $P_{\mathrm{c}\mathrm{h},i,t}^{\mathrm{E}\mathrm{V}}$, $P_{\mathrm{d}\mathrm{i}\mathrm{s},i,t}^{\mathrm{E}\mathrm{V}}$,$P_{i,t}^{\mathrm{E}\mathrm{V}}$, $E_{i,t}^{\mathrm{E}\mathrm{V}}$ denote the the charging power, discharging power, actual power, and operation capacity; ${\eta }_{\mathrm{c}\mathrm{h}}^{\mathrm{E}\mathrm{V}}$ and ${\eta }_{\mathrm{d}\mathrm{i}\mathrm{s}}^{\mathrm{E}\mathrm{V}}$ are the charge and discharge efficiency coefficient of BIEV, (4c) and (4d) represent the boundary of operation capacity and operation state.

Then, considering the power boundary $P_{i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{E}\mathrm{V}}$ and$P_{i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{E}\mathrm{V}}$of BIEV is affected by the operation capacity, the upgraded flexible energy supply $F_{i,t}^{\mathrm{E}\mathrm{V},\mathrm{s},\mathrm{u}\mathrm{p}}$ and downgraded flexible energy supply $F_{i,t}^{\mathrm{E}\mathrm{V},\mathrm{s},\mathrm{d}\mathrm{n}}$ of EV in the smart building system can be expressed as,

$$F_{i,t}^{\mathrm{EV},\mathrm{s},\mathrm{up}}=\min \{P_{i,t,\max }^{\mathrm{EV}}-P_{i,t}^{\mathrm{EV}},\frac{E_{i,t,\max }^{\mathrm{EV}}-E_{i,t}^{\mathrm{EV}}}{\Delta t}\}$$

(5a)

$$F_{i,t}^{\mathrm{EV},\mathrm{s},\mathrm{dn}}=\min \{P_{i,t}^{\mathrm{EV}}-P_{i,t,\min }^{\mathrm{EV}}P,\frac{E_{i,t}^{\mathrm{EV}}-E_{i,t,\min }^{\mathrm{EV}}}{\Delta t}\}$$

(5b)

2.4. VB model of TCL in the smart building system

The equivalent thermal parameter model is adopted to describe the thermal dynamics characteristics of TCL, which equivalents internal environment, external environment, and heating (cooling) capacity to circuit devices to analyze the dynamic change relationship of indoor temperature and power [18,19]. As illustrated in Figure 3, the heating capacity $Q_{i,t}^1$of the temperature control load can be defined as the sum of the heat absorbed by the room $Q_{i,t}^2$and the air convection heat $Q_{i,t}^3$. When the ambient temperature $T_{\mathrm{o}\mathrm{u}\mathrm{t},i,t}^{\mathrm{T}\mathrm{C}\mathrm{L}}$ is higher than the indoor temperature $T_{\mathrm{i}\mathrm{n},i,t}^{\mathrm{T}\mathrm{C}\mathrm{L}}$, the temperature of the heating equipment $T_{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t},i,t}^{\mathrm{T}\mathrm{C}\mathrm{L}}$ and the ambient temperature combine to raise the indoor temperature, increasing the energy of the VB model, as follows,

$$\{\begin{array}{l}{Q}_{i,t}^1=Q_{i,t}^2+Q_{i,t}^3\\ Q_{i,t}^1=A_1{\eta }^{\mathrm{TCL}}(T_{\mathrm{heat},i,t}^{\mathrm{TCL}}-T_{\mathrm{in},i,t}^{\mathrm{TCL}})\\ Q_{i,t}^2=\rho \mathrm{c}V(T_{\mathrm{in},i,t+1}^{\mathrm{TCL}}-T_{\mathrm{in},i,t}^{\mathrm{TCL}})\\ Q_{i,t}^3=A_2{\eta }^{\mathrm{Wall}}(T_{\mathrm{in},i,t}^{\mathrm{TCL}}-T_{\mathrm{out},i,t}^{\mathrm{TCL}})\\ {\eta }^{\mathrm{COP}}=\frac{P_{i,t}^{\mathrm{TCL}}\Delta t}{Q_{i,t}^1}\end{array}$$

(6a)

$$P_{i,t,\min }^{\mathrm{TCL}}\le P_{i,t}^{\mathrm{TCL}}\le P_{i,t,\max }^{\mathrm{TCL}}$$

(6b)

$$T^{\mathrm{base}}-\sigma \le T_{\mathrm{in},i,t}^{\mathrm{TCL}}\le T^{\mathrm{base}}+\sigma$$

(6c)

$$T_{\mathrm{in},i,t+1}^{\mathrm{TCL}}=(1-\frac{A_2{\eta }^{\mathrm{Wall}}}{\rho \mathrm{c}V})T_{\mathrm{in},i,t}^{\mathrm{TCL}}+\frac{A_2{\eta }^{\mathrm{Wall}}}{\rho \mathrm{c}V}{T}_{\mathrm{out},i,t}^{\mathrm{TCL}}+\frac{P_{i,t}^{\mathrm{TCL}}}{\rho \mathrm{c}V{\eta }^{\mathrm{COP}}}\Delta t$$

(6d)

where $A_1$ and $A_2$ are the heat dissipation area of BITC connected at node $i$ and time $t$ and the wall; ${\eta }^{\mathrm{T}\mathrm{C}\mathrm{L}}$, ${\eta }^{\mathrm{W}\mathrm{a}\mathrm{l}\mathrm{l}}$, and ${\eta }^{\mathrm{C}\mathrm{O}\mathrm{P}}$ denote the heating efficiency, wall surface thermal radiation rate, and energy conversion efficiency; $\rho$, $c$, $V$ are the indoor air density, air heat capacity, and room volume. (6b) and (6c) represent the boundary of actual power and the indoor temperature, where $T^{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}$ and $\sigma$ are the indoor temperature baseline and dead zone value.

The upgraded flexible energy supply $F_{i,t}^{\mathrm{T}\mathrm{C}\mathrm{L},\mathrm{s},\mathrm{u}\mathrm{p}}$ and downgraded flexible energy supply $F_{i,t}^{\mathrm{T}\mathrm{C}\mathrm{L},\mathrm{s},\mathrm{d}\mathrm{n}}$ of TCL in the smart building system can be expressed as,

$$F_{i,t}^{\mathrm{TCL},\mathrm{s},\mathrm{up}}=P_{i,t,\max }^{\mathrm{TCL}}-P_{i,t}^{\mathrm{TCL}}$$

(7a)

$$F_{i,t}^{\mathrm{TCL},\mathrm{s},\mathrm{dn}}=P_{i,t}^{\mathrm{TCL}}-P_{i,t,\min }^{\mathrm{TCL}}$$

(7b)

3. Active-reactive power collaborative optimization scheduling of distribution network cluster with the aggregation of flexible DERs in smart building system

3.1. Cluster division index of building-integrated flexible DERs

Structural and functional cluster division indexes of the building-integrated flexible DERs are designed for distribution network with high penetration of renewable energy to achieve independent autonomy within clusters and coordinated interaction among clusters. In order to ensure the close electrical connection of the nodes within the clusters, the module degree index based on the electrical distance is adopted to describe the connection strength of the nodes in structure [8]. Functional index refers to the flexible balance contribution index, which can evaluate the regulation capacity of the building-integrated flexible DERs to ensure the power balance under the high penetration rate of distributed energy access.

3.1.1. Module degree index

Generally, electrical distance $L_{ij}$ is exploited to measure the tightness of electrical coupling between nodes in the distribution network based on the sensitivity relationship of node power variation and node voltage [8], and thus the module degree index $\xi$ for the cluster division structure strength of complex networks is formulated as follows,

$$\xi =\frac{1}{2m}{\displaystyle \sum _i{\displaystyle \sum _j\left({\mu }_{ij}-\frac{r_i{r}_j}{2m}\right)\delta (i,j)}}$$

(8a)

$$\{\begin{array}{l}{\mu }_{ij}=1-\frac{L_{ij}}{L_{ij,\max }}\\ m={\displaystyle \sum _i{\displaystyle \sum _j{\mu }_{ij}}}\\ r_i={\displaystyle \sum _i{\mu }_{ij}}\end{array}$$

(8b)

$$L_{ij}=\sqrt{{(d_{i1}-d_{j1})}^2+{(d_{i2}-d_{j2})}^2+\mathrm{...}+{(d_{in}-d_{jn})}^2}$$

(8c)

$$d_{ij}=(S_{ii}^P-S_{ij}^P)+(S_{ii}^Q-S_{ij}^Q)$$

(8d)

$$\Delta U=S^P\Delta P+S^Q\Delta Q$$

(8e)

where ${\mu }_{ij}$, $m$, and $r_i$ denote the edge weight connecting node $i$ and $j$, sum of the edge weight in the network, edge weight connecting node $i$; $\delta \left(i,j\right)$, a 0-1 variable, equals to 1 if node 1 and node 2 are in the same cluster. Since there is a coupling relationship between the two nodes and the surrounding nodes, $L_{ij}$ can be formulated with the combined effect $d_{ij}$ of the power change of node $i$ and $j$ as shown in (8c) and (8d); $L_{ij,\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum threshold of the electrical distance; (8e) indicates the sensitivity relationship between node voltage and node injected power; $∆U$, $∆P$, and $∆Q$are the variation of the voltage amplitude, active power and reactive power; $S^P$ and $S^Q$ are the sensitivity matrix of voltage-active power and voltage-reactive power.

3.1.2. Flexible balance contribution index

The regulation capacity of the building-integrated flexible DERs refers to its flexible energy supply capacity to satisfy the flexible demands of net load in the smart building system, which can be divided into upgraded flexible energy supply capacity and downgraded flexible energy supply capacity. Flexible balance contribution index $\phi$ can measure the matching degree between distributed power consumption demands and the supply of the building-integrated flexible DERs through calculating the equilibrium proportion of the flexible DERs bearing the node flexibility demands at a certain moment,

$$\phi =\frac{1}{n}({\displaystyle \sum _{i=1}^n{\alpha }_{i,t}^{G,\mathrm{up}}}+{\displaystyle \sum _{i=1}^n{\alpha }_{i,t}^{G,\mathrm{dn}}})$$

(9a)

$$\{\begin{array}{l}{\alpha }_{i,t}^{G,\mathrm{up}}=\frac{F_{i,t}^{G,\mathrm{s},\mathrm{up}}}{F_{i,t}^{\mathrm{dem}}}\\ {\alpha }_{i,t}^{G,\mathrm{dn}}=\frac{F_{i,t}^{G,\mathrm{s},\mathrm{dn}}}{F_{i,t}^{\mathrm{dem}}}\end{array}\mathrm{G}\in \left\{\mathrm{ES},\mathrm{EV},\mathrm{TCL}\right\}$$

(9b)

$$F_{i,t}^{\mathrm{dem}}=P_{i,t}^{\mathrm{Load}}-P_{i,t}^{\mathrm{PV}}$$

(9c)

where ${\alpha }_{i,t}^{G,\mathrm{u}\mathrm{p}}$ and ${\alpha }_{i,t}^{G,\mathrm{d}\mathrm{n}}$ denote the upgraded and downgraded flexible balance contribution degree; (9c) represents the flexible demands of net load in the smart building system, deriving from the difference between the load demands $P_{i,t}^{\mathrm{L}\mathrm{o}\mathrm{a}\mathrm{d}}$and BIPV generation.

Then, the cluster division index $\gamma$ of the building-integrated flexible DERs considering the module degree index and the flexible balance contribution index with different weights $a_1$ and $a_2$ can be calculated from formula (8a) and formula (9a),

$$\gamma =a_1\xi +a_2\phi$$

(10)

The critical procedures of the cluster division process of the building-integrated flexible DERs are summarized as follows,

• Each node in the distribution network is regarded as a separate cluster initially to calculate the cluster division index $\gamma$ of the building-integrated flexible DERs based on formula (10);
• Node $j$ is randomly selected from the remaining other nodes to merge with node $i$ forming a new cluster $k(i,j)$, and then the variation of the cluster division index $∆\gamma {=\gamma }^{\text{'}}-\gamma$can be derived from the division index of cluster $k(i,j)$. When $∆\gamma$ reaches the maximum positive value, the two nodes can be divided into the same cluster;
• The new cluster $k(i,j)$ will be regarded as a new node to repeating the second procedure. The division procedures will be stopped until the nodes in the network cannot merge and the cluster division index $\gamma$ of the building-integrated flexible DERs reaches the maximum.

3.2. Active-reactive power collaborative optimization scheduling within clusters

3.2.1. Optimization objective

The active-reactive power collaborative optimization scheduling within clusters can be formulated as a multi-object optimization model based on the flexibility of the building-integrated flexible DERs to handle the uncertainties from the renewable energy. Specifically, each cluster $k$ aims to minimize the operational cost $f_k^1$, power loss cost $f_k^2$, and penalty cost for flexibility deficiency $f_k^3$, as shown in (11a), (12a), and (13a). $f_k^1$ is composed of the maintenance cost $C_{k,i,t}^{\mathrm{P}\mathrm{V}}$ of BIPVs due to the renewable generation curtailment, lifetime degradation cost $C_{k,i,t}^{\mathrm{E}\mathrm{S}}$ of BIESs, compensation cost $C_{k,i,t}^{\mathrm{C}\mathrm{L}}$ of BIEVs and BITCs stemmed from involving power regulation, and the operation cost $C_{k,i,t}^{\mathrm{S}\mathrm{V}\mathrm{C}}$ of static var compensator (SVC) from the uncertainties considering the renewable generation. $f_k^2$ can be derived from the power flow calculation during the scheduling circulation, where $U_{k,i,t}$, $G_{k,i,t}$, $B_{k,i,t}$, ${\delta }_{k,ij}$ represent the nodal voltage magnitudes, conductance, susceptance, and phase angle of line $ij$ at time $t$ within cluster $k$; ${\lambda }_k^{\mathrm{L}\mathrm{o}\mathrm{s}\mathrm{s}}$ is the compensation price of network loss. The penalty cost $f_k^3$ for flexibility deficiency is introduced when the flexibility supply cannot satisfy the load demands within the clusters, including the upgraded flexibility deficiency $F_{k,t}^{\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{i},\mathrm{u}\mathrm{p}}$ and downgraded flexibility deficiency $F_{k,t}^{\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{i},\mathrm{d}\mathrm{n}}$, as follows,

$$f_k^1={\displaystyle \sum _{t\in \Gamma }\left({\displaystyle \sum _{i\in {\Omega }^{\mathrm{PV}}}^{N_k}{C}_{k,i,t}^{\mathrm{PV}}}+{\displaystyle \sum _{i\in {\Omega }^{\mathrm{ES}}}^{N_k}{C}_{k,i,t}^{\mathrm{ES}}}+{\displaystyle \sum _{i\in {\Omega }^{\mathrm{CL}}}^{N_k}{C}_{k,i,t}^{\mathrm{CL}}}+{\displaystyle \sum _{i\in {\Omega }^{\mathrm{SVC}}}^{N_k}{C}_{k,i,t}^{\mathrm{SVC}}}\right)}k=1,2,\mathrm{3...},K$$

(11a)

$$C_{k,i,t}^{\mathrm{PV}}={\lambda }_k^{\mathrm{PV}}{P}_{k,i,t}^{\mathrm{PV}}={\lambda }_k^{\mathrm{PV}}\left(P_{k,i,t}^{\mathrm{PV}0}-\Delta P_{k,i,t}^{\mathrm{PV}}\right)$$

(11b)

$$C_{k,i,t}^{\mathrm{ES}}={\lambda }_k^{\mathrm{ES}}\left|P_{k,i,t}^{\mathrm{ES}}\right|$$

(11c)

$$C_{k,i,t}^{\mathrm{CL}}={\lambda }_k^{\mathrm{EV}}\left|P_{k,i,t}^{\mathrm{EV}}\right|+{\lambda }_k^{\mathrm{TCL}}\left|P_{k,i,t}^{\mathrm{TCL}}\right|$$

(11d)

$$C_{k,i,t}^{\mathrm{SVC}}={\lambda }_k^{\mathrm{SVC}}\left|Q_{k,i,t}^{\mathrm{SVC}}\right|$$

(11e)

$$f_k^2={\displaystyle \sum _{t\in \Gamma }\left({\lambda }_k^{\mathrm{Loss}}{\displaystyle \sum _{(i,j)\in N_k}{G}_{k,ij}\left({\left(U_{k,i,t}\right)}^2+{\left(U_{k,j,t}\right)}^2-2U_{k,i,t}{U}_{k,j,t}\cos {\delta }_{k,ij}\right)}\right)}k=1,2,\mathrm{3...},K$$

(12a)

$$f_k^3={\displaystyle \sum _{t\in \Gamma }\left({\lambda }_k^{\mathrm{defi}}(F_{k,t}^{\mathrm{defi},\mathrm{up}}+F_{k,t}^{\mathrm{defi},\mathrm{dn}})\right)}k=1,2,\mathrm{3...},K$$

(13a)

$$F_{k,t}^{\mathrm{defi},\mathrm{up}}={\displaystyle \sum _{i\in N_k}\Delta }{F}_{k,i,t}^{\mathrm{up}}=\{\begin{array}{l}{\displaystyle \sum _{i\in N_k}\left|F_{k,i,t}^{G,\mathrm{s},\mathrm{up}}-F_{k,i,t}^{\mathrm{dem}}\right|}\Delta F_{k,i,t}^{\mathrm{up}}\le 0\\ 0\Delta F_{k,i,t}^{\mathrm{up}}\ge 0\end{array}$$

(13b)

$$F_{k,t}^{\mathrm{defi},\mathrm{dn}}={\displaystyle \sum _{i\in N_k}\Delta }{F}_{k,i,t}^{\mathrm{dn}}=\{\begin{array}{l}{\displaystyle \sum _{i\in N_k}\left|F_{k,i,t}^{G,\mathrm{s},\mathrm{dn}}-F_{k,i,t}^{\mathrm{dem}}\right|}\Delta F_{k,i,t}^{\mathrm{dn}}\le 0\\ 0\Delta F_{k,i,t}^{\mathrm{dn}}\ge 0\end{array}$$

(13c)

where ${\mathsf{\Omega }}^{\mathrm{P}\mathrm{V}}$, ${\mathsf{\Omega }}^{\mathrm{E}\mathrm{S}}$, ${\mathsf{\Omega }}^{\mathrm{C}\mathrm{L}}$, and ${\mathsf{\Omega }}^{\mathrm{S}\mathrm{V}\mathrm{G}}$ are the node sets of BIPV, BIES, BIEV, BITC, and SVC; $N_k$ and $K$ indicate the $k$-th distributed resource cluster and total number of divided clusters; ${\lambda }_k^{\mathrm{P}\mathrm{V}}$, ${\lambda }_k^{\mathrm{E}\mathrm{S}}$, ${\lambda }_k^{\mathrm{E}\mathrm{V}}$, ${\lambda }_k^{\mathrm{T}\mathrm{C}\mathrm{L}}$, and ${\lambda }_k^{\mathrm{S}\mathrm{V}\mathrm{C}}$ represent unit operation prices of BIPV, BIES, BIEV, BITC, and SVC connected at node $i$ and time $t$ within cluster $k$; $P_{k,i,t}^{\mathrm{P}\mathrm{V}}$, $P_{k,i,t}^{\mathrm{P}\mathrm{V}0}$, and $\mathsf{\Delta }{P}_{k,i,t}^{\mathrm{P}\mathrm{V}}$ denote the actual generation, predicted generation, and generation curtailment of BIPV. Due to the frequent charging/discharging behaviors of the building-integrated flexible DERs, the aging of these resources is inevitable, where $P_{k,i,t}^{\mathrm{E}\mathrm{S}}$, $P_{k,i,t}^{\mathrm{E}\mathrm{V}}$, $P_{k,i,t}^{\mathrm{T}\mathrm{C}\mathrm{L}}$ represent the actual power of BIES, BIEV, BITC connected at node $i$ and time $t$ within cluster $k$; $Q_{k,i,t}^{\mathrm{S}\mathrm{V}\mathrm{C}}$ is the reactive power output of SVC. (13b) and (13c) indicate the upgraded and downgraded flexibility deficiency, deriving from the flexibility margin $\mathsf{\Delta }{F}_{k,i,t}^{\mathrm{u}\mathrm{p}}$ and $\mathsf{\Delta }{F}_{k,i,t}^{\mathrm{d}\mathrm{n}}$ of the negative value, and ${\lambda }_k^{\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{i}}$ is the penalty factor for the flexibility deficiency.

Considering the operational cost, power loss cost, and penalty cost for flexibility deficiency comprehensively, the multi-objective function with different weights ${\lambda }_1$, ${\lambda }_2$ and ${\lambda }_3$ can be represented as,

$$\min f_k={\lambda }_1{f}_k^1+{\lambda }_2{f}_k^2+{\lambda }_3{f}_k^3$$

(14)

3.2.2. Operation constraints

Newton-Laphson method is adopted to describe the power flows within the clusters, as shown in (15a). Constraint (15b) represents the active and reactive power balance at each bus, and constraint (15c) imposes the upper and lower limits on the nodal voltage magnitudes, $U_{k,i,t,\mathrm{m}\mathrm{a}\mathrm{x}}$ and $U_{k,i,t,\mathrm{m}\mathrm{i}\mathrm{n}}$.

$$\{\begin{array}{l}\Delta P_{k,i,t}-U_{k,i,t}{\displaystyle \sum _{(i,j)\in N_k}{U}_{k,j,t}\left(G_{k,ij}\cos {\delta }_{k,ij}+B_{k,ij}\sin {\delta }_{k,ij}\right)}=0\\ \Delta Q_{k,i,t}-U_{k,i,t}{\displaystyle \sum _{(i,j)\in N_k}{U}_{k,j,t}\left(G_{k,ij}\sin {\delta }_{k,ij}-B_{k,ij}\cos {\delta }_{k,ij}\right)}=0\end{array}$$

(15a)

$$\{\begin{array}{l}\Delta P_{k,i,t}=P_{k,,i,t}^{\mathrm{PV}}+P_{k,i,t}^{\mathrm{ES}}-P_{k,i,t}^{\mathrm{EV}}-P_{k,i,t}^{\mathrm{TCL}}-P_{k,i,t}^{\mathrm{Load}}\\ \Delta Q_{k,i,t}=Q_{k,i,t}^{\mathrm{PV}}+Q_{k,i,t}^{\mathrm{ES}}+Q_{k,i,t}^{\mathrm{SVC}}\end{array}$$

(15b)

$$U_{k,i,t,\min }\le U_{k,i,t}\le U_{k,i,t,\max }$$

(15c)

The power output curtailment of BIPV is supposed to be lower than the specified maximum thresholds of the active and reactive power, $\mathsf{\Delta }{P}_{k,i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{P}\mathrm{V}}$ and $\mathsf{\Delta }{Q}_{k,i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{P}\mathrm{V}}$. Moreover, sufficient power factor $\mathrm{c}\mathrm{o}\mathrm{s}\theta$ for the operation of photovoltaic inverter is needed to limit the reactive power circulation, where the amount of reactive power injected or absorbed is constrained by the bound (16a) and (16b),

$$（P_{k,i,t}^{\mathrm{PV}}{）}^2+（Q_{k,i,t}^{\mathrm{PV}}{）}^2\le （S_{k,i,t}^{\mathrm{PV}}{）}^2$$

(16a)

$$\cos \theta \sqrt{（P_{k,i,t}^{\mathrm{PV}}{）}^2+（Q_{k,i,t}^{\mathrm{PV}}{）}^2}\le P_{k,i,t}^{\mathrm{PV}}$$

(16b)

$$0\le \Delta P_{k,i,t}^{\mathrm{PV}}\le \Delta P_{k,i,t,\max }^{\mathrm{PV}}$$

(16c)

$$0\le \Delta Q_{k,i,t}^{\mathrm{PV}}\le \Delta Q_{k,i,t,\max }^{\mathrm{PV}}$$

(16d)

Also, the state of reactive power of SVC should be constrained within their lower and upper bounds $Q_{k,i,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{S}\mathrm{V}\mathrm{C}}$ and $Q_{k,i,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{S}\mathrm{V}\mathrm{C}}$ , as follows,

$$Q_{k,i,t,\min }^{\mathrm{SVC}}\le Q_{k,i,t}^{\mathrm{SVC}}\le Q_{k,i,t,\max }^{\mathrm{SVC}}$$

(17)

In addition, the sum of the regulation capacity of the building-integrated flexible DERs based on the constrained space superposition method with Minkowski sum $⨁$ can be expressed as formula (18a), and the state of the total sum of the regulation capacity $P_{k,t}^{\mathrm{F}\mathrm{l}\mathrm{e}\mathrm{x}}$ should always be limited within the allowable lower and upper bounds $P_{k,t,\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{F}\mathrm{l}\mathrm{e}\mathrm{x}}$ and $P_{k,t,\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{F}\mathrm{l}\mathrm{e}\mathrm{x}}$,

$$P_{k,t}^{\mathrm{Flex}}=P_{k,t}^{\mathrm{ES}}\oplus P_{k,t}^{\mathrm{EV}}\oplus P_{k,t}^{\mathrm{TCL}}={\displaystyle \sum _{i\in N_k}{P}_{k,i,t}^{\mathrm{ES}}}+{\displaystyle \sum _{i\in N_k}{P}_{k,i,t}^{\mathrm{EV}}}+{\displaystyle \sum _{i\in N_k}{P}_{k,i,t}^{\mathrm{TCL}}}$$

(18a)

$$P_{k,t,\min }^{\mathrm{Flex}}\le P_{k,t}^{\mathrm{Flex}}\le P_{k,t,\max }^{\mathrm{Flex}}$$

(18b)

$$\{\begin{array}{l}{P}_{k,t,\max }^{\mathrm{Flex}}={\displaystyle \sum _{i\in N_k}{P}_{i,t,\max }^{\mathrm{ES}}}+{\displaystyle \sum _{i\in N_k}{P}_{i,t,\max }^{\mathrm{EV}}}+{\displaystyle \sum _{i\in N_k}{P}_{i,t,\max }^{\mathrm{TCL}}}\\ P_{k,t,\min }^{\mathrm{Flex}}={\displaystyle \sum _{i\in N_k}{P}_{i,t,\min }^{\mathrm{ES}}}+{\displaystyle \sum _{i\in N_k}{P}_{i,t,\min }^{\mathrm{EV}}}+{\displaystyle \sum _{i\in N_k}{P}_{i,t,\min }^{\mathrm{TCL}}}\end{array}$$

(18c)

3.2.3. Solution methodology

Due to the nonlinearity and nonconvexity of the power flow equations increasing the heavy computation burden of optimization model, the second-order cone-based branch flow method is exploited to decompose this challenging problem into the linearized cone model by introducing intermediate variable $X_{k,i,t}$, $Y_{k,i,t}$, $Z_{k,i,t}$ and cone constraints, as shown in (19),

$$X_{k,i,t}={(U_{k,i,t})}^2$$

(19a)

$$Y_{k,i,t}=U_{k,i,t}{U}_{k,j,t}\cos {\delta }_{k,ij}$$

(19b)

$$Z_{k,i,t}=U_{k,i,t}{U}_{k,j,t}\sin {\delta }_{k,ij}$$

(19c)

$$\sqrt{{(Y_{k,i,t})}^2+{(Z_{k,i,t})}^2+{\left(\frac{X_{k,i,t}-X_{k,j,t}}{2}\right)}^2}\le \frac{X_{k,i,t}+X_{k,j,t}}{4}$$

(19d)

Then the power loss cost function (12a), the power flow equations (15a), and the nodal voltage magnitudes constraint (15c) can be reformulated as follows,

$$f_k^2={\displaystyle \sum _{t\in \Gamma }\left({\rho }_k^{\mathrm{Loss}}{\displaystyle \sum _{(i,j)\in N_k}{G}_{k,ij}\left(X_{k,i,t}+X_{k,j,t}-2Y_{k,i,t}\right)}\right)}k=1,2,\mathrm{3...},K$$

(20)

$$\{\begin{array}{l}\Delta P_{k,i,t}-\left(G_{k.ii}{X}_{k,i,t}+{\displaystyle \sum _{\begin{array}{l}(i,j)\in N_k\\ j\ne 1\end{array}}^n\left(G_{k,ij}{Y}_{k,i,t}+B_{k,ij}{Z}_{k,i,t}\right)}\right)=0\\ \Delta Q_{k,i,t}-\left(-B_{k,ii}{X}_{k,i,t}-{\displaystyle \sum _{\begin{array}{l}(i,j)\in N_k\\ j\ne 1\end{array}}^n\left(B_{k,ij}{Y}_{k,i,t}-G_{k,ij}{Z}_{k,i,t}\right)}\right)=0\end{array}$$

(21)

$${(U_{k,i,t,\min })}^2\le X_{k,i,t}\le {(U_{k,i,t,\max })}^2$$

(22)

4. Discussion

4.1. System data

The proposed active-reactive power collaborative optimization scheduling model with cluster division is tested on a modified IEEE 33-bus distribution system, as depicted in Figure 4. The installed capacities of BIPVs in nodes 1, 6, 13, 17, 18, 24, 29 are set as 0.5MW, 1.0MW, 0.2MW, 1.5MW, 0.8MW, 1.2MW, 0.4MW to simulate the diverse PV endowments. Figure 5 illustrates the predicted energy outputs of BIPVs derived from the actual values in [6,20,21] and the estimated 24-h loads with a system standard load of 3.175+j2.3MVA [22,23]. The prediction errors of BIPVs generation and node loads follow Gaussian distributions with zero means and 10% standard deviations of the prediction values [21,22,24]. The detailed technical specifications of building-integrated equipment and distribution network are summarized in Table 1 [3,4,24].

The appropriate weights for the cluster division index of the building-integrated DERs and the multi-object function are set as 0.41, 0.29; 0.63, 0.48, and 0.25 after a large number of simulation tests, respectively. The optimization scheduling of distribution network with cluster division for flexible DERs integrated into the smart building system is implemented on an hourly basis over a 24-h horizon, and solved by the CPLEX Optimizer on MATLAB platform.

4.2. Comparative Results and Analysis

Three comparative schemes are performed for in depth analysis on the effectiveness of the proposed methodology.

• Scheme 1 performs the proposed optimization scheduling of distribution network with cluster division for flexible DERs integrated into the smart building system in Section 3.
• Scheme 2 adopts the centralized optimization scheduling of distribution network without considering the cluster division.
• Scheme 3 is the initial distribution network before the optimization scheduling.

4.2.1. Cluster division of the building-integrated flexible DERs

Consider the differences in node voltage amplitude, power injection, and load demands at different times, the detailed data of each node at 12:00 noon is exploited for cluster division. Figure 6 depicts the function curves for the cluster division index of the building-integrated DERs corresponding to the number of different clusters. It can be found that the index function achieves the maximum value $\gamma =0.92$ when the number of clusters is divided into 6. Therefore, the result of the optimal cluster division at $K=6$ is shown in Figure 7,

4.2.2. Flexible regulation capacity results of the building-integrated flexible DERs

Figure 8 illustrates the flexible regulation capacity of the building-integrated flexible DERs based on the flexibility constraint space superimposed by Minkowski sum. It can be observed from Figure 8 that the self-consistent rate is allowed to reflect the contribution of different flexibility resources to the flexibility balance of the power grid. For instance, the flexibility power region of building-integrated flexible DERs in cluster 2 is larger than that in cluster 1. In addition, the sum of the total amount of BIPV generation and the regulation capacity of the building-integrated flexible DERs can satisfy the total load demands within the cluster 2, indicating that cluster 2 achieves the complete self-consistency of the active power. On the contrary, the total load demands within the cluster 1 is greater than the sum of the total amount of BIPV generation and the regulation capacity of the building-integrated flexible DERs at time $t=4$ and $t=20$, thus give rise to the flexibility power deficiency.

4.2.3. Optimization scheduling results with scheme 1-3

The power loss and the node voltage amplitude with schemes 1-3 are illustrated in Figure 9, where the node voltage at 12:00 noon is selected for analysis since the generation of BIPVs reaches the peak at that time. It can be found that the optimized system power loss with scheme 1 is lower than that with scheme 2 due to the comprehensive scheduling of the BIPVs generation and the regulation capacity of building-integrated flexible DERs within clusters. In addition, the voltage regulation effect of the centralized active power-reactive power cooperative optimization of the distribution network is unsatisfactory, remaining the node voltage overlimit. However, since the building-integrated flexible DERs are involved in the cluster power regulation, the node voltage amplitude with scheme 1 is within the specified range. Table 2 summarizes the statistical data of comparative economic performance results with schemes1 and 2, including the total cost, operational cost, power loss cost, and flexibility deficiency cost. Compared to scheme 2, scheme 1 considering flexibility of building-integrated flexible DERs can better estimate the flexibility regulation capacity, and thus the operational cost and flexibility deficiency cost can be reduced by 11.69% and 9.21%, respectively. Consequently, scheme 1 outperforms schemes 2 and schemes 3 on the economic merits and voltage regulation voltage regulation.

5. Conclusions

In this paper, an optimal active-reactive power collaborative scheduling model for the building flexible DER clusters is proposed to resolve the uncertainties from the renewable energy in distribution network. The key findings of this study are as follows: (1) The developed cluster division algorithm considering structural and functional indexes can achieve independent autonomy within clusters and coordinated interaction among clusters, contributing to decreasing the whole network loss; (2) Since the customized virtual battery models integrate and quantify the regulation capacity of flexible resources, the flexibility deficiency cost can be reduced by 9.21%; (3) The developed optimal active-reactive power collaborative scheduling model can achieve a better performance in economic merits and voltage regulation voltage regulation than the centralized optimization scheduling of distribution network without considering the cluster division, whose operational cost can be decreased by 11.69%. In this article, the building-integrated flexible DERs including BIES, BIEV, and BITC were acknowledged as the resources in demand side, further on research will focus on the collaborate flexibility regulation capacity of building-integrated flexible DERs in the supply side and demand side.