Adaptive robust economic dispatch and real-time control of distribution system considering controllable inverter air-conditioner clusters

With a tremendous number of renewable energy sources (RES) integrated into the distribution system, the inherent uncertainty of RES power generation brings about significant challenges in distribution and power balance within the distribution system. This article proposes an adaptive robust economic dispatch (ARED) model and a real-time control strategy for distribution systems as countermeasures, which make full use of the adjustable capabilities of controllable inverter air-conditioner (IAC) clusters. Firstly, the concept of the adjustable capacity curve (ACC) is developed to accurately quantify the adjustable capacity of an IAC cluster. Afterward, a two-stage adaptive robust optimization is formulated for ARED, which comprehensively takes the adjustable capacity of the IAC cluster and the uncertainty of RES into consideration. Meanwhile, the solution methodology of ARED is also designed based on the column and constraint generation (C&CG) algorithm, where the master problem is quadratic programming with quadratic constraints (QCQP), and the max-min sub-problem is reformulated to a mixed integer linear programming (MILP) form by taking advantage of linear duality theory and big-M method. Finally, a novel real-time decentralized control strategy for IAC clusters is also proposed for purpose of hedging against stochastic RES power fluctuation after every round of ARED decisions. The results of the case study validate the effectiveness of ARED model and real-time control strategy under different uncertainty scenarios of RES power generation.

1 Introduction

In recent years, in order to suppress fossil energy depletion and extreme climate change, innovative technologies and applications of energy and power in industry are going through a highly active period (Dai et al., 2017). Clean and low-carbon energyutilization technologies have attracted broad attention from all countries, emphasizing the dominant position of renewable energy sources (RES) in future power system (Lu et al., 2017).

However, with the growing proportion of RES engaging in distribution system operation, significant challenges are brought about to distribution system operators (DSOs). On the one hand, during the power dispatching stage, traditional prediction value-based deterministic economic dispatch may suffer from suboptimality (Li et al., 2015) due to the inevitable prediction errors for RES power generation. Therefore, an uncertainty-based economic dispatch model is required, aiming to determine a relatively cost-efficient reference point for all dispatchable units. Stochastic optimization (SO) (Chen et al., 2018; Ye et al., 2018; Roldán et al., 2019; Shuai et al., 2019), chance-constrained programming (Zhou et al., 2018; Yang et al., 2021) and robust optimization (RO) (Gao et al., 2018; Liu et al., 2018; Chen et al., 2021) are three leading paradigms adopted widely in previous studies (Ning and You, 2018). Among these paradigms, SO is mainly formulated based on finite typical scenarios and corresponding discrete probability distribution after scenario generation and reduction. In the work of Ye et al. (2018), authors propose an optimal dispatch method for a power system with high penetration of wind power. Typical scenarios are generated by assuming the prediction error obeys the Gaussian probability density function. In the work of Chen et al. (2018), a generative adversarial network based data-driven approach is adopted to generate typical scenarios of RES power generation. A SO-based optimization model usually contains a relatively large number of constraints, and the applicability of the solution depends on the accuracy of the probability distribution. Distinct from SO, the two-stage adaptive RO does not require accurate probability distribution of uncertainty variables, thus is more practicable for systems with limited knowledge of stochastic RES characteristics (Zhang et al., 2019). Liu et al. (2018) develop a two-stage RO model to minimize the operation cost of the microgrid, in which uncertainties of RES and load are modeled by box constraints. In the study of Chen et al. (2021), a robust dynamic economic dispatching model of integrated transmission and distribution system is established, in which uncertainty set is decided according to conditional value-at-risk rather than manually assigned. Gao et al. (2018) explore the coordinated energy management among distribution systems and networked microgrids, in which a tractable typical budget set is adopted. To improve the conservativeness of RO’s solution, in these works, different uncertainty sets with adjustable uncertainty budgets are employed, such as polyhedron uncertainty set (Gao et al., 2018; Liu et al., 2018; Chen et al., 2021), ellipsoid uncertainty set (Roldán et al., 2019). In addition, drawing advantages of both RO and SO, distributionally robust optimization (DRO) also attracts growing attention in recent years, which builds an ambiguity set of uncertainty variables based on statistical moments information and optimizes the system’s expected cost under the worst distribution. In the work of Ruan et al. (2019), the ambiguity set is developed by 1-norm and ∞-norm constraints of the distance between stochastic probability distribution and historical statistical distribution. Zhou et al. (2018) develop a distributionally robust dispatch model considering generalized moments uncertainty of wind power and equivalently transform it to deterministic quadratically constrained quadratic programming (QCQP) problem.

On the other hand, during the real-time operation stage, actual RES power generation deviates from the predicted value with stochastic amplitude and high variation frequency, which may lead to fierce fluctuation of exchange power between the distribution system and transmission system. Due to this reason, more ancillary resources are required for the transmission system to ensure real-time power balance (Kalantar-Neyestanaki et al., 2020). However, under severe operation scenarios such as a contingency state, the regulatable capacity of the transmission system may become extremely limited to meet the requirements, which greatly stimulates the need for ancillary resources from the distribution system side (Tan et al., 2020). To release the power balance burden of the transmission system, demand-side controllable load clusters have been considered to be important regulatory resources with great potential. Among different types of controllable loads, a thermostatically controllable load such as heat pumps and inverter air-conditioners (IAC) are the most characteristic types, due to their large scale, low cost, and fast response speed (Yao and Zhang, 2018). In the work of Wei et al. (2016), a hierarchical and distributed control strategy of heat pumps is proposed to balance the fluctuations caused by RES, in which each pump’s target power is determined by solving the optimal setting temperature. Ding et al. (2021) also propose a game-theoretic demand side management strategy, which is capable of guiding users to make optimal power consumption schedules and smooth the tie-line power of microgrids simultaneously. Hui et al. (2019) propose equivalent modeling of IAC to provide frequency regulation service and verify the aggregation of IACs can be controlled in a similar way to traditional generators. Based on networked information interaction between each IAC aggregator, Jiang and Wei (2018) presents adistributed cooperation model and optimal control strategy for an IAC cluster and the power grid to avoid the occurrence of a new load peak. Based on a transactive control framework, Yao and Zhang (2018) propose a decentralized control strategy of the IAC cluster to provide peak shaving service, and furtherly explore the coordination of heterogeneous thermostatically controlled loads to provide real-time ancillary services in their follow-up study (Yao and Zhang, 2020).

Actually, the power-adjustable capability of thermostatic load and user’s comfort experience is highly related to its temperature state, hence very sensitive to power adjustment instructions of both economic dispatch stage and real-time control stage. In most of the previous studies of economic dispatch, controllable loads are scheduled like normal generators according to their power-adjustable capacity, without exploring the effect on users’ comfort experience and variation of adjustable capability. In addition, the issues of how to reduce communication costs and protect users’ private information should also be subtly considered when aggregating large-scale controllable loads. To address these issues, this paper first develops the concept of an adjustable capacity curve, which is capable of quantifying the exact adjustable capacity of each IAC cluster within a certain time period. Then, for the power dispatch stage, we propose an adaptive robust economic dispatch (ARED) model to achieve cost-efficient dispatching of all generation units by taking IAC clusters’ adjustable capability and RES uncertainty into consideration. Afterward, for the real-time operation stage, a real-time control strategy is designed, where IAC clusters function as ancillary resources to hedge against the actual power fluctuation of RES based on decentralized control. Variation of power exchange between the transmission system and distribution system is effectively eliminated. The main contributions of this paper are summarized by the following aspects:

1) Sophisticated IAC model is first proposed in this paper, which not only contains common operation constraints in existing works (Yao and Zhang, 2018; Hui et al., 2019) but also takes the nonlinear relationship of IAC’s compressor operation frequency and IAC’s electrical/cooling power into consideration. More accurate quantification of IAC’s power adjustable capacity and tractable decentralized control can be realized based on this model.

2) A two-stage adaptive robust optimization model is formulated for ARED in this paper by drawing the advantages of RO, where the uncertainty of RES is addressed by adjustable box constraints rather than by discrete scenarios in the work of Ye et al. (2018). ARED model also considers the IAC cluster as a dispatchable unit whose adjustable capability is dynamically updated along with real-time control processes of the IAC cluster. In addition, different from most literature such as Chen et al. (2021); Gao et al. (2018) where upward and downward power adjustment are modeled by two independent positive-define variables, we propose an equivalent formation by only one variable, which can simplify the expressions of objective function and constraints.

3) A decentralized real-time control strategy capable of aggregating and cooperatively controlling IACs is proposed for the purpose of compensating power deviation caused by stochastic RES power fluctuation. IAC clusters’ power adjustable capacity can be quantified and updated rapidly during real-time control, which also functions as a constraints boundary in ARED model. Most of the computation is conducted locally, and less remote communication cost is required compared with the work of Wei et al. (2016); Ding et al. (2021). Privacy protection of IACs’ parameters and users’ preference information can also be guaranteed under the designed strategy.

The remaining content of this paper is organized as follows. Section 2 introduces the model of IAC and the quantification method of the IAC cluster’s adjustable capability. In Section 3, a two-stage adaptive robust optimization model for ARED considering IAC adjustable capacity is formulated and a C&CG algorithm-based solution methodology of ARED is elaborated as well. In Section 4, the real-time control strategy of the IAC cluster is proposed to hedge against RES power fluctuation during the real-time operation stage of the distribution system. In Section 5, a case study is explored to verify the effectiveness of ARED model and real-time control strategy of the IAC cluster. Section 6 concludes this paper.

2 Inverter air-conditioner modeling and adjustable capability quantification

2.1 IAC model

For a room equipped with an IAC, the first-order equivalent thermal parameter (ETP) dynamic model of room temperature T_a in continuous form can be expressed by (1):

$\frac{d{T}_{\text{a}}}{dt}=-\frac{U_{\mathrm{o}\mathrm{a}}{A}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}+c_{\text{a}}{\rho }_{\text{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}\xi }{c_{\text{a}}{\rho }_{\text{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}}\left(T_{\text{a}}-T_{\mathrm{o}\mathrm{u}\mathrm{t}}\right)-\frac{1}{c_{\text{a}}{\rho }_{\text{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}}{q}_{\mathrm{a}\mathrm{c}}(1)$

In (1), c_a and ρ_a stand for heat capacity and density of air; A_room and V_room stand for area and volume of the room; U_oa is heat transfer coefficient of the room; ξ is the number of air exchange times per hour; T_out is outdoor ambient temperature; q_ac is cooling power of IAC. For specific sampling and control period τ, the discrete-time formulation of (1) is derived and applied in this paper, where we introduce equivalent heat dissipation coefficient $\sigma =\left(U_{\mathrm{o}\mathrm{a}}{A}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}+c_{\mathrm{a}}{\rho }_{\mathrm{a}}V\xi \right)/c_{\mathrm{a}}{\rho }_{\mathrm{a}}V$ for simplicity of expression. The room temperature at discrete time point kτ = 0, τ, 2τ, … can be calculated by:

$T_{\mathrm{a}}\left(k+1\right)=e^{-\sigma \tau }{T}_{\mathrm{a}}\left(k\right)+\frac{e^{-\sigma \tau }-1}{c_{\mathrm{a}}{\rho }_{\mathrm{a}}V\sigma }{q}_{\mathrm{a}\mathrm{c}}\left(k\right)+\left(1-e^{-\sigma \tau }\right)T_{\mathrm{o}\mathrm{u}\mathrm{t}}(2)$

For one IAC, the consumed electrical power p_ac and produced cooling power q_ac are actually adjusted by operation frequency f_ac of IAC compressor. However, the relationship between p_ac, q_ac and frequency f_ac are usually nonlinear, which means the adjustment sensitivity ∂p_ac/∂f_ac, ∂q_ac/∂f_ac are diverse at different operation frequency point. In this paper, we adopt piecewise linearization technology to approximate the original nonlinear functions, as shown in Figure 1. Suppose the operation frequency f_ac is segmented into N intervals according to (N+1) segmentation points (f¹, … , f^N+1), and so are p_ac and q_ac. Then, the slope vectors ${\mathbf{k}}_{\mathrm{p}\mathrm{s}\mathrm{e}\mathrm{g}}={(k_{\mathrm{p}}^i)}_{\mathrm{N}\times \mathrm{1}},{\mathbf{k}}_{\mathrm{q}\mathrm{s}\mathrm{e}\mathrm{g}}={(k_{\mathrm{q}}^i)}_{\mathrm{N}\times \mathrm{1}}$ concerning segment i can be calculated. Meanwhile, let ${\mathbf{f}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}={(f_{\mathrm{a}\mathrm{c}}^1,\dots ,f_{\mathrm{a}\mathrm{c}}^N)}^{\mathrm{T}}$ be the vector of the forward N segmentation points of f_ac, and ${\mathbf{f}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{B}}={(f_{\mathrm{a}\mathrm{c}}^{\mathrm{2}},\dots ,f_{\mathrm{a}\mathrm{c}}^{N+1})}^{\mathrm{T}}$ be the backward N segmentation points. In a similar way, we can as well define p_segF and q_segF. Drawing support from these definition of necessary coefficient vectors, the relationship between p_ac, q_ac and f_ac can be formulated in a compact form by a set of equalities and inequalities:

$q_{\mathrm{a}\mathrm{c}}\left(k\right)={{\mathbf{q}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}}^{\mathrm{T}}\omega \left(k\right)+{\mathbf{k}}_{\mathrm{q}\mathrm{s}\mathrm{e}\mathrm{g}}^{\mathrm{T}}\lambda \left(k\right)(3)$

$p_{\mathrm{a}\mathrm{c}}\left(k\right)={{\mathbf{p}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}}^{\mathrm{T}}\omega \left(k\right)+{\mathbf{k}}_{\mathrm{p}\mathrm{s}\mathrm{e}\mathrm{g}}^{\mathrm{T}}\lambda \left(k\right)(4)$

$0\le \lambda \left(k\right)\le \left({\mathbf{\Lambda }}_{\mathrm{B}}-{\mathbf{\Lambda }}_{\mathrm{F}}\right)\omega \left(k\right)(5)$

$1^{\mathrm{T}}\omega \left(k\right)=1(6)$

$f_{\mathrm{a}\mathrm{c}}\left(k\right)={\mathbf{f}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}^{\mathrm{T}}\omega \left(k\right)+{\mathrm{1}}^{\mathrm{T}}\lambda \left(k\right)(7)$

$f_{\mathrm{a}\mathrm{c}}^{\min }\le f_{\mathrm{a}\mathrm{c}}\left(k\right)\le f_{\mathrm{a}\mathrm{c}}^{\max }(8)$

$\left|f_{\mathrm{a}\mathrm{c}}\left(k+1\right)-f_{\mathrm{a}\mathrm{c}}\left(k\right)\right|\le r(9)$

FIGURE 1

FIGURE 1. Piecewise linearization of nonlinear p_ac—f_ac function.

In (3)-(7), binary vector $\omega (k)={({\omega }_i(k))}_{\mathrm{N}\times \mathrm{1}}$ and continuous auxiliary vector $\lambda (k)={({\lambda }^i(k))}_{\mathrm{N}\times \mathrm{1}}$ are introduced to transform nonlinear equalities p_ac = g_p (f_ac), q_ac = g_q (f_ac) to mixed integer linear equalities and inequalities. ω(k) contains only one element 1, whose position indicates the segment where the operation frequency $f_{\mathrm{a}\mathrm{c}}\left(k\right)$ locates. Λ_F = diag (f_segF) and Λ_B = diag (f_segB) are two diagonal matrix, and $\left({\mathbf{\Lambda }}_{\mathrm{B}}-{\mathbf{\Lambda }}_{\mathrm{F}}\right)\omega \left(k\right)$ determines the upper bound of auxiliary vector λ(k). Constraint (8) represents the upper and lower bound of $f_{\mathrm{a}\mathrm{c}}\left(k\right)$, and constraint (9) represents frequency regulation limitations during discrete control interval.

2.2 Adjustable capability of one IAC

For each room equipped with an IAC, suppose there always exists an acceptable room temperature range [T_min, T_max] and a most satisfactory temperature T_fit according to user’s personal preference. So, electrical power p_ac of IAC can be regulated within a certain range without violating the user’s comfort limitations, which provides substantial flexible power adjustment capability to the distribution system after aggregating.

To evaluate the power adjustment capability of one IAC and an IACs cluster, we firstly introduce a piecewise linear mapping $\phi :\mathbb{R}\mapsto \mathbb{R}$ of room temperature T_a:

$\phi \left(T_{\mathrm{a}}\right)=\left\{\begin{array}{c}\left(T_{\mathrm{a}}-T_{\mathrm{fi}\mathrm{t}}\right)/\left(T_{\max }-T_{\mathrm{fi}\mathrm{t}}\right),if{T}_{\mathrm{a}}\ge T_{\mathrm{fi}\mathrm{t}}\hfill \\ \left(T_{\mathrm{a}}-T_{\mathrm{fi}\mathrm{t}}\right)/\left(T_{\mathrm{fi}\mathrm{t}}-T_{\min }\right),if{T}_{\mathrm{a}}<T_{\mathrm{fi}\mathrm{t}}\hfill \end{array}\right.(10)$

Thus, the mapping result $s=\phi \left(T_{\mathrm{a}}\right)$ will be a scalar within the range of [ − 1, 1], which will be positive if room temperature T_a is above T_fit, and be negative if T_a is below T_fit. Meanwhile, smaller absolute value of s indicates that T_a is closer to T_fit, hence $\left|s\right|$ can be regarded as an index to quantify user’s “unsatisfactory degree” [15] of current room temperature T_a.

Then, we define a monotonically decreasing function p_ac = d_ac(s) termed adjustment capability curve (ACC) of one IAC according to user’s preference information, which is illustrated in Figure 2. The horizontal axis represents the target “unsatisfactory degree” s (k + 1) after current round of control, and vertical axis represents the electrical power p_ac required to transfer s(k) to s (k + 1). Intuitively, image of ACC is uniquely determined according to 4 characteristic points (F, A, B, C) and 2 horizontal power boundary values $p_{\mathrm{a}\mathrm{c}}^{\max \text{\_}\mathrm{o}\mathrm{p}\mathrm{t}},p_{\mathrm{a}\mathrm{c}}^{\min \text{\_}\mathrm{o}\mathrm{p}\mathrm{t}}$. For point F, its abscissa value indicates target “unsatisfactory degree” s (k + 1) is fixed to current s(k), and its ordinate value $p_{\mathrm{a}\mathrm{c}}^{\mathrm{fi}\mathrm{x}}$ represents the required electrical power to guarantee s (k + 1) = s(k) after control period τ. Similarly, for point A, it indicates the electrical power $p_{\mathrm{a}\mathrm{c}}^{\max \text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}$ required to ensure s (k + 1) = −1; for point B, it indicates the electrical power $p_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{i}\mathrm{n}\text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}$ to satisfy s (k + 1) = +1; for pointC, it indicates the electrical power $p_{\mathrm{a}\mathrm{c}}^{\mathrm{s}\mathrm{a}\mathrm{t}}$ to guarantee s (k + 1) = 0. Except for these 4 characteristic electrical power values $\{p_{\mathrm{a}\mathrm{c}}^{\mathrm{fi}\mathrm{x}},p_{\mathrm{a}\mathrm{c}}^{\max \text{\_}\mathrm{c}\mathrm{f}\mathrm{t}},p_{\mathrm{a}\mathrm{c}}^{\min \text{\_}\mathrm{c}\mathrm{f}\mathrm{t}},p_{\mathrm{a}\mathrm{c}}^{\mathrm{s}\mathrm{a}\mathrm{t}}\}$, p_ac is also restricted by IAC’s upper and lower bound of electrical power $\{p_{\mathrm{a}\mathrm{c}}^{\max \text{\_}\mathrm{o}\mathrm{p}\mathrm{t}},p_{\mathrm{a}\mathrm{c}}^{\min \text{\_}\mathrm{o}\mathrm{p}\mathrm{t}}\}$ due to physical operation limitations. Consequently, p_ac = d_ac(s) presents to be a monotonically decreasing piecewise linear function.

FIGURE 2

FIGURE 2. Adjustment capability curve (ACC) of one IAC.

To obtain adjustable capacity of one IAC, firstly we should calculate corresponding characteristic cooling power values $\{q_{\mathrm{a}\mathrm{c}}^{\mathrm{fi}\mathrm{x}},q_{\mathrm{a}\mathrm{c}}^{\max \text{\_}\mathrm{c}\mathrm{f}\mathrm{t}},q_{\mathrm{a}\mathrm{c}}^{\min \text{\_}\mathrm{c}\mathrm{f}\mathrm{t}},q_{\mathrm{a}\mathrm{c}}^{\mathrm{s}\mathrm{a}\mathrm{t}}\}$ according to discrete ETP Equation 2:

$q_{\mathrm{a}\mathrm{c}}^{\mathrm{fi}\mathrm{x}}\left(k\right)=c_{\mathrm{a}}{\rho }_{\mathrm{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}\sigma \left(T_{\mathrm{o}\mathrm{u}\mathrm{t}}-T_{\mathrm{a}}\left(k\right)\right)(11)$

$q_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}\text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}\left(k\right)=c_{\mathrm{a}}{\rho }_{\mathrm{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}\sigma \left(T_{\mathrm{o}\mathrm{u}\mathrm{t}}-\frac{T^{\min }-e^{-\sigma \tau }{T}_{\mathrm{a}}\left(k\right)}{1-e^{-\sigma \tau }}\right)(12)$

$q_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{i}\mathrm{n}\text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}\left(k\right)=c_{\mathrm{a}}{\rho }_{\mathrm{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}\sigma \left(T_{\mathrm{o}\mathrm{u}\mathrm{t}}-\frac{T^{\max }-e^{-\sigma \tau }{T}_{\mathrm{a}}\left(k\right)}{1-e^{-\sigma \tau }}\right)(13)$

$q_{\mathrm{a}\mathrm{c}}^{\mathrm{s}\mathrm{a}\mathrm{t}}\left(k\right)=c_{\mathrm{a}}{\rho }_{\mathrm{a}}{V}_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{m}}\sigma \left(T_{\mathrm{o}\mathrm{u}\mathrm{t}}-\frac{T^{\mathrm{fi}\mathrm{t}}-e^{-\sigma \tau }{T}_{\mathrm{a}}\left(k\right)}{1-e^{-\sigma \tau }}\right)(14)$

Then, by referring to aforementioned segmentation of IAC power and operation frequency, binary vector ω(k) and continuous auxiliary vector λ(k) can be determined according to current cooling power q_ac(k) via constraints (3)–(7). Afterward, characteristic power $p_{\mathrm{a}\mathrm{c}}^{•}\left(k\right)$ (• stands for {fix, max_cft, min_cft, sat}) can be calculated:

$p_{\mathrm{a}\mathrm{c}}^{•}\left(k\right)={{\mathbf{p}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}}^{\mathrm{T}}\omega \left(k\right)+\frac{{\mathbf{k}}_{\mathrm{p}\mathrm{s}\mathrm{e}\mathrm{g}}^{\mathrm{T}}\omega \left(k\right)}{{\mathbf{k}}_{\mathrm{q}\mathrm{s}\mathrm{e}\mathrm{g}}^{\mathrm{T}}\omega \left(k\right)}\left[q_{\mathrm{a}\mathrm{c}}^{•}\left(k\right)-{{\mathbf{q}}_{\mathrm{s}\mathrm{e}\mathrm{g}\mathrm{F}}}^{\mathrm{T}}\omega \left(k\right)\right](15)$

Finally, adjustable electrical power $p_{\mathrm{a}\mathrm{c}}^{\mathrm{\Delta }}(k)$ of one IAC is quantified by:

$\left\{\begin{array}{c}{p}_{\mathrm{a}\mathrm{c}}^{\mathrm{\Delta }}\left(k\right)\ge \max \left\{p_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{i}\mathrm{n}\text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}\left(k\right),p_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{i}\mathrm{n}\text{\_}\mathrm{o}\mathrm{p}\mathrm{t}}\left(k\right)\right\}-p_{\mathrm{a}\mathrm{c}}\left(k\right)\hfill \\ p_{\mathrm{a}\mathrm{c}}^{\mathrm{\Delta }}\left(k\right)\le \min \left\{p_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}\text{\_}\mathrm{c}\mathrm{f}\mathrm{t}}\left(k\right),p_{\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}\text{\_}\mathrm{o}\mathrm{p}\mathrm{t}}\left(k\right)\right\}-p_{\mathrm{a}\mathrm{c}}\left(k\right)\hfill \end{array}\right.(16)$

2.3 Adjustable capability of IAC cluster

During system-level dynamic economic dispatching and real-time control processes, the IAC cluster is the basic unit participating in power balance. So, an IAC cluster is working like a virtual power plant (VPP) to some extent, which is ought to provide its power-adjustable capacity and current operation power to DSO.

To achieve this goal, in this paper, we require all IACs belonging to the same cluster to follow a unified control signal s^∗(k + 1). The advantages are twofold: (1) Unified “unsatisfactory degree” ensures fairness of different users’ comfort experience without violating their personal preference, because s^∗(k + 1) ∈ [ − 1, 1]; (2) This requirement is sufficiently simple and tractable, which brings great convenience to decision-making process even if the number of IACs in one cluster is very large. Once the requirement is satisfied, the adjustment capability of the IAC cluster can be easily computed by simply aggregating (adding up) each IAC’s ACC as an equivalent ACC of the IAC cluster:

$p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}=d_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}\left(s\right)\stackrel{\mathrm{\Delta }}{=}\sum _{i\in \mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}{d}_{\mathrm{a}\mathrm{c},i}\left(s\right)(17)$

So, (17) indicates that p_cluster is also a monotonically decreasing piecewise linear function concerning cluster’s unified control signal s. Suppose $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}\in [p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\min },p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\max }]$, the adjustable capacity $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta }}$ of IAC cluster can be quantified similarly with (16):

$\left\{\begin{array}{c}{p}_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\min }\le p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta }}\le p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\max }\hfill \\ p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\min }=\min \left\{0,p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\min }-p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}\left(k\right)\right\}\hfill \\ p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\max }=\max \left\{0,p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\max }-p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}\left(k\right)\right\}\hfill \end{array}\right.(18)$

where $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}(k)=\sum _i{p}_{\mathrm{a}\mathrm{c},i}(k)$ is the electrical power demand of IAC cluster at current time point kτ. That’s to say, once p_cluster(k) and cluster’s ACC are determined, adjustable capacity $[p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\min },p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\max }]$ of IAC cluster can be quickly determined according to (17) (18).

3 Adaptive robust economic dispatch optimization and solution methodology

Due to the lack of uncertainty modeling of RES, under some extreme operation scenarios, the decision results of a traditional deterministic economic dispatching model based on predicted values may lead to high adjustment costs during real-time power balance. To remedy this defect, in this section, we propose a two-stage adaptive robust economic dispatch optimization model ARED in (19), which takes both uncertainty of RES and adjustable capability of IAC clusters into consideration:

$\underset{y^{\mathrm{f}}}{\min }{Q}^{\mathrm{f}}+\underset{u}{\max }\underset{y^{\mathrm{s}}}{\min }{Q}^{\mathrm{s}}(19)$

The object function of ARED is composed of two parts: normal-state operation cost Q^f (first stage decision) and re-dispatch cost Q^s under worst operation scenario (second stage decision). The decision variables of these two stages are defined by y^f, y^s respectively, and u denotes uncertainty variables, which is modeled within an adjustable uncertainty set.

3.1 First stage decision

3.1.1 Decision variables definition

In the first stage decision, the decision variables are:

$y^{\mathrm{f}}=\left\{\begin{array}{c}{v}_j|\forall j\in V\hfill \\ p_{ij};q_{ij};l_{ij}|\forall \left(i,j\right)\in \xi \hfill \\ p_{\mathrm{g},j};q_{\mathrm{g},j}|\forall j\in G\hfill \end{array}\right\}(20)$

where V stands for the set of all busbars, ξ stands for the set of all distribution lines, and G represents the set of all power generation units (generators and wind plants); Decision variable v is the square of busbar voltage, {p_ij; q_ij; l_ij} stands for active power, reactive power and square of current magnitude of directed distribution line (i, j), and p_g is active power generation of power generation unit.

3.1.2 Objective function

The object of the first stage decision is to minimize power generation cost and wind curtailment punishment according to predicted power generation of wind plants $w_{\mathrm{g}}^{\mathrm{e}}$:

$\min Q^{\mathrm{f}}=\sum _{j\in \left\{G\backslash W\right\}}\left(c_{2,j}{p}_{\mathrm{g},j}^2+c_{1,j}{p}_{\mathrm{g},j}+c_{0,j}\right)+\sum _{j\in W}{\rho }_{\mathrm{c}\mathrm{r}}\left(w_{\mathrm{g},j}^{\mathrm{e}}-p_{\mathrm{g},j}\right)(21)$

where c₂, c₁, c₀ are coefficient of generator’s quadratic generation cost function; ρ_cr is punishment coefficient of wind curtailment; W is the set of wind plants, which is a subset of G.

3.1.3 Constraints

$p_{\mathrm{g},j}-p_{\mathrm{d},j}+\sum _{\left(i,j\right)\in \xi }\left(p_{ij}-r_{ij}{l}_{ij}\right)=\sum _{\left(j,k\right)\in \xi }{p}_{jk}+g_j{v}_j(22)$

$q_{\mathrm{g},j}-q_{\mathrm{d},j}+\sum _{\left(i,j\right)\in \xi }\left(q_{ij}-x_{ij}{l}_{ij}\right)=\sum _{\left(j,k\right)\in \xi }{q}_{jk}+b_j{v}_j(23)$

$v_j=v_i+\left(r_{ij}^2+x_{ij}^2\right)l_{ij}-2\left(r_{ij}{p}_{ij}+x_{ij}{q}_{ij}\right)(24)$

${‖2p_{ij};2q_{ij};l_{ij}-v_i‖}_2\le l_{ij}+v_i(25)$

${\left(V_j^{\min }\right)}^2\le v_j\le {\left(V_j^{\max }\right)}^2(26)$

$\begin{array}{c}\left\{\begin{array}{c}{p}_{\mathrm{g},j}^{\min }\le p_{\mathrm{g},j}\le p_{\mathrm{g},j}^{\max }\hfill \\ q_{\mathrm{g},j}^{\min }\le q_{\mathrm{g},j}\le q_{\mathrm{g},j}^{\max }\hfill \end{array}\right.j\in G\hfill \\ \left\{\begin{array}{c}0\le p_{\mathrm{g},j}\le w_{\mathrm{g},j}^{\mathrm{e}}\hfill \\ p_{\mathrm{g},j}\tan {\varphi }^{\min }\le q_{\mathrm{g},j}\le p_{\mathrm{g},j}\tan {\varphi }^{\max }\hfill \end{array}\right.j\in W\hfill \end{array}(27)$

$\left\{\begin{array}{c}-p_{ij}^{\max }\le p_{ij}\le p_{ij}^{\max }\hfill \\ -p_{ij}^{\max }\le p_{ij}-l_{ij}{r}_{ij}\le p_{ij}^{\max }\hfill \\ -q_{ij}^{\max }\le q_{ij}\le q_{ij}^{\max }\hfill \\ -q_{ij}^{\max }\le q_{ij}-l_{ij}{x}_{ij}\le q_{ij}^{\max }\hfill \\ l_{ij}\ge 0\hfill \end{array}\right.(28)$

For ∀j ∈ V, constraints (22) (23) stand for nodal active power and reactive power balance equations respectively; For ∀(i, j) ∈ ξ, we specially use second-order cone relaxation (SOCR) inequality (25) to substitute original quadratic equation $v_i{l}_{ij}={(p_{ij})}^2+{(q_{ij})}^2$, which is widely used in modeling radial distribution network (Farivar and Low, 2013); Constraint (26) restricts upper and lower bound of busbar voltage; Boundary constraints of all power generation units and power factor constraints of wind plants are presented in (27), where ϕ ∈ [ϕ^min,ϕ^max] is allowable power factor angle. Constraints (28) represents bidirectional power flow limitations.

3.2 Second stage decision

3.2.1 Decision variables definition

In the second stage decision, to ensure power balance, re-dispatch of all generation units, controllable load clusters (i.e., IAC clusters in this paper), wind power curtailment, and load shedding are deployed according to actual RES power generation during real-time operation. This stage is regarded as a “corrective stage” based on the first stage decision results. The decision variables are listed in (29), where we use superscript Δ to indicate the variation of the corresponding decision variable, which can be positive or negative. In addition, two non-negative variables $p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$ and $p_{d,j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$ are introduced, which equal the absolute values of $p_{\mathrm{g},j}^{\mathrm{\Delta }}$ and $p_{\mathrm{d},j}^{\mathrm{\Delta }}$.

$y^{\mathrm{s}}=\left\{\begin{array}{c}{v}_j^{\mathrm{\Delta }};|\forall j\in V\hfill \\ p_{ij}^{\mathrm{\Delta }};q_{ij}^{\mathrm{\Delta }};l_{ij}^{\mathrm{\Delta }}|\forall \left(i,j\right)\in \xi \hfill \\ p_{\mathrm{g},j}^{\mathrm{\Delta }};q_{\mathrm{g},j}^{\mathrm{\Delta }};p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}|\forall j\in G\hfill \\ w_{g,j}^{\mathrm{c}\mathrm{r}}|\forall j\in W\hfill \\ p_{\mathrm{d},j}^{\mathrm{\Delta }};q_{\mathrm{d},j}^{\mathrm{\Delta }};p_{d,j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}|\forall j\in L\hfill \end{array}\right\}(29)$

3.2.2 Objective function

$\begin{array}{c}\underset{u}{\max }\underset{y^{\mathrm{s}}}{\min }{Q}^{\mathrm{s}}=\sum _{j\in \left\{G\backslash W\right\}}{\rho }_{\mathrm{a}\mathrm{d},j}{p}_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}+\sum _{j\in W}{\rho }_{\mathrm{c}\mathrm{r},j}{w}_{\mathrm{g},j}^{\mathrm{c}\mathrm{r}}\hfill \\ +\sum _{k\in CL}{\rho }_{\mathrm{d},j}{p}_{\mathrm{d},k}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}+\sum _{k\in \left\{L\backslash CL\right\}}{\rho }_{\mathrm{l}\mathrm{s},j}{p}_{\mathrm{d},k}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}\hfill \end{array}(30)$

The object function is to minimize the total cost of corrective actions under the worst uncertainty scenario u ∈ U, where U is the adjustable uncertainty set. CL stands for the set of all controllable IAC clusters. Coefficient {ρ_ad, ρ_cr, ρ_d, ρ_ls} represents power adjustment cost of generation units, curtailment cost of wind plants, power adjustment cost of controllable IAC clusters and load shedding cost of normal loads.

3.2.3 Constraints

$p_{\mathrm{g},j}^{\mathrm{\Delta }}+p_{\mathrm{d},j}^{\mathrm{\Delta }}+\sum _{\left(i,j\right)\in \xi }\left(p_{ij}^{\mathrm{\Delta }}-r_{ij}{l}_{ij}^{\mathrm{\Delta }}\right)=\sum _{\left(j,k\right)\in \xi }{p}_{jk}^{\mathrm{\Delta }}+g_j{v}_j^{\mathrm{\Delta }}(31)$

$q_{\mathrm{g},j}^{\mathrm{\Delta }}+q_{\mathrm{d},j}^{\mathrm{\Delta }}+\sum _{\left(i,j\right)\in \xi }\left(q_{ij}^{\mathrm{\Delta }}-x_{ij}{l}_{ij}^{\mathrm{\Delta }}\right)=\sum _{\left(j,k\right)\in \xi }{q}_{jk}^{\mathrm{\Delta }}+b_j{v}_j^{\mathrm{\Delta }}(32)$

$v_j^{\mathrm{\Delta }}=v_i^{\mathrm{\Delta }}+\left(r_{ij}^2+x_{ij}^2\right)l_{ij}^{\mathrm{\Delta }}-2\left(r_{ij}{p}_{ij}^{\mathrm{\Delta }}+x_{ij}{q}_{ij}^{\mathrm{\Delta }}\right)(33)$

$2p_{ij}{p}_{ij}^{\mathrm{\Delta }}+2q_{ij}{q}_{ij}^{\mathrm{\Delta }}-l_{ij}{v}_i^{\mathrm{\Delta }}-v_i{l}_{ij}^{\mathrm{\Delta }}=0(34)$

$\begin{array}{c}\left\{\begin{array}{c}{p}_{\mathrm{g},j}^{\min }\le p_{\mathrm{g},j}+p_{\mathrm{g},j}^{\mathrm{\Delta }}\le p_{\mathrm{g},j}^{\max }\hfill \\ q_{\mathrm{g},j}^{\min }\le q_{\mathrm{g},j}+q_{\mathrm{g},j}^{\mathrm{\Delta }}\le q_{\mathrm{g},j}^{\max }\hfill \end{array}\right.j\in G\hfill \\ \left\{\begin{array}{c}{p}_{\mathrm{g},j}+p_{\mathrm{g},j}^{\mathrm{\Delta }}=w_{\mathrm{g},j}^{\mathrm{g}}-w_{\mathrm{g},j}^{\mathrm{c}\mathrm{r}}\hfill \\ w_{\mathrm{g},j}^{\mathrm{c}\mathrm{r}}\ge 0\hfill \\ \left(p_{\mathrm{g},j}+p_{\mathrm{g},j}^{\mathrm{\Delta }}\right)\tan {\varphi }^{\min }\le q_{g,j}+q_{\mathrm{g},j}^{\mathrm{\Delta }}\hfill \\ q_{g,j}+q_{\mathrm{g},j}^{\mathrm{\Delta }}\le \left(p_{\mathrm{g},j}+p_{\mathrm{g},j}^{\mathrm{\Delta }}\right)\tan {\varphi }^{\max }\hfill \end{array}\right.j\in W\hfill \end{array}(35)$

$\begin{array}{c}{\left(V_j^{\min }\right)}^2\le v_j+v_j^{\mathrm{\Delta }}\le {\left(V_j^{\max }\right)}^{2}j\in V\hfill \\ \left\{\begin{array}{c}-p_{ij}^{\max }\le \left(p_{ij}+p_{ij}^{\mathrm{\Delta }}\right)\le p_{ij}^{\max }\hfill \\ -p_{ij}^{\max }\le \left(p_{ij}+p_{ij}^{\mathrm{\Delta }}\right)-\left(l_{ij}+l_{ij}^{\mathrm{\Delta }}\right)r_{ij}\le p_{ij}^{\max }\hfill \\ -q_{ij}^{\max }\le \left(q_{ij}+q_{ij}^{\mathrm{\Delta }}\right)\le q_{ij}^{\max }\hfill \\ -q_{ij}^{\max }\le \left(q_{ij}+q_{ij}^{\mathrm{\Delta }}\right)-\left(l_{ij}+l_{ij}^{\mathrm{\Delta }}\right)x_{ij}\le q_{ij}^{\max }\hfill \\ l_{ij}+l_{ij}^{\mathrm{\Delta }}\ge 0\hfill \end{array}\right.\left(i,j\right)\in \xi \hfill \end{array}(36)$

$\left\{\begin{array}{c}{\beta }^{\min }{p}_{\mathrm{d},j}\le p_{\mathrm{d},j}^{\mathrm{\Delta }}\le {\beta }^{\max }{p}_{\mathrm{d},j}\hfill \\ {\beta }^{\min }{q}_{\mathrm{d},j}\le q_{\mathrm{d},j}^{\mathrm{\Delta }}\le {\beta }^{\max }{q}_{\mathrm{d},j}\hfill \end{array}\right.(37)$

$\left\{\begin{array}{c}{w}_{\mathrm{g},j}^{\mathrm{g}}=w_{\mathrm{g},j}^{\mathrm{e}}+\left({\theta }_{\mathrm{g},j}^{+}-{\theta }_{\mathrm{g},j}^{-}\right)\mathrm{\Delta }{w}_{\mathrm{g},j}\hfill \\ 0\le {\theta }_{\mathrm{g},j}^{+},{\theta }_{\mathrm{g},j}^{-}\le 1\hfill \end{array}\right.(38)$

$\left\{\begin{array}{c}-p_{\mathrm{d},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}\le p_{\mathrm{d},j}^{\mathrm{\Delta }}\le p_{\mathrm{d},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}\hfill \\ -p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}\le p_{\mathrm{g},j}^{\mathrm{\Delta }}\le p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}\hfill \end{array}\right.(39)$

In these constraints, {p_g, q_g, p_ij, q_ij, l_ij,v} are derived by the solutions of first stage decision model, and should be considered as constant known values during second stage decision. Constraints (31)–(36) basically inherit first stage constraints (22)–(28), but are presented in form of incremental expressions. (34) is derived through deploying first-order Taylor expansion on original constraints $v_i{l}_{ij}={(p_{ij})}^2+{(q_{ij})}^2$ at the operation point determined by first stage decision, which becomes a linear equation. (37) represents the adjustable capacity of IAC clusters or load shedding capacity of normal loads, where β^min, β^max are two scale factors representing adjustable capacity (18) with their absolute value smaller than 1. For IAC clusters, they can not only reduce power demand, but also increase power demand, so β^min ≤ 0, β^max ≥ 0; For normal loads, they can only reduce power demand (load shedding), hence β^min = 0, β^max ≥ 0. (38) establishes a box-style uncertainty set to depict the relationship between actual power generation of wind plant $w_{\mathrm{g}}^{\mathrm{g}}$ and predicted value $w_{\mathrm{g}}^{\mathrm{e}}$. The maximum deviation Δw_g is adjustable, which can be used to obtain and compare solutions under different conservative levels. Auxiliary constraint (39) implicitly establishes non-negative restriction of variables $p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$, $p_{\mathrm{d},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$, and guarantees $p_{\mathrm{g},j}^{\mathrm{\Delta }}=p_{\mathrm{g},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$,$p_{\mathrm{d},j}^{\mathrm{\Delta }}=p_{\mathrm{d},j}^{\mathrm{\Delta },\mathrm{a}\mathrm{b}\mathrm{s}}$ at optimal solution.

Different from the model formulation in previous work (Gao et al., 2018), the model in this paper not only requires fewer decision variables but also eliminates integer variables indicating upward or downward power adjustment states, which contributes to the simplification of ARED solution.

3.3 C&CG algorithm based ARED solution method

In our two-stage ARED optimization problem, the first stage decision is a QCQP problem due to SOCR of branch power flow equations. The second stage decision is a max-min problem, which can be regarded as a linear programming problem (LP) concerning decision variable y^s for given y^f and uncertainty scenario u. Considering the problem structure of ARED, we draw experience from column and constraints generation (C&CG) algorithm (Zeng and Zhao, 2013) to solve it, which iteratively solves the master problem (MP) and sub-problem (SP) to approximate the optimal solution. The compact form of MP and SP are:

$\begin{array}{ll}\hfill & \begin{array}{c}\text{MP}:\underset{y^{\mathrm{f}},\eta ,y^{\mathrm{s},l}}{\min }{Q}^{\mathrm{f}}+\eta \hfill \\ \mathrm{s}.\mathrm{t}\text{.\hspace{0.17em}}\left\{\begin{array}{c}\begin{array}{c}\hfill \eta \ge b^{\mathrm{T}}{y}^{\mathrm{s},l}\hfill \\ \hfill \begin{array}{c}{E}_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{f}}+G_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{s},l}\le h_{\mathrm{i}\mathrm{e}}-M_{\mathrm{i}\mathrm{e}}{u}^l\hfill \\ E_{\mathrm{e}}{y}^{\mathrm{f}}+G_{\mathrm{e}}{y}^{\mathrm{s},l}=h_{\mathrm{e}}-M_{\mathrm{e}}{u}^l\hfill \end{array}\hfill \end{array}\text{\hspace{0.17em},\hspace{0.17em}}l=1,\dots ,k\hfill \\ y^{\mathrm{f}}\in {\mathrm{\Omega }}^{\mathrm{f}}\hfill \end{array}\right.\hfill \end{array}\hfill \end{array}(40)$

$\begin{array}{c}\text{SP}:\underset{u}{\max }\underset{y^{\mathrm{s}}}{\min }{Q}^{\mathrm{s}}=b^{\mathrm{T}}{y}^{\mathrm{s}}\hfill \\ \mathrm{s}.\mathrm{t}\text{.\hspace{0.17em}}\left\{\begin{array}{c}\hfill Qu\le q\hfill \\ \hfill \begin{array}{c}{G}_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{s}}\le h_{\mathrm{i}\mathrm{e}}-E_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{f}}-M_{\mathrm{i}\mathrm{e}}u\hfill \\ G_{\mathrm{e}}{y}^{\mathrm{s}}=h_{\mathrm{e}}-E_{\mathrm{e}}{y}^{\mathrm{f}}-M_{\mathrm{e}}u\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}(41)$

In MP, y^f ∈ Ω^f stands for feasible region of y^f generated by constraints (22)–(28); The first three constraints are gradually appended to MP after solving SP during each iteration, and k is current iteration times. In SP, $u\stackrel{\mathrm{\Delta }}{=}{({({\theta }^{+})}^{\mathrm{T}},{({\theta }^{+})}^{\mathrm{T}})}^{\mathrm{T}}$, and constraint Qu ≤ q stands for adjustable uncertainty set (38), and two remained equality and inequality correspond to constraints (31)–(37), (39).

To solve max-min problem SP, we firstly reformulate it to a single max-problem by deriving the dual problem of inner min-problem, which is named by DSP1:

$\begin{array}{c}\text{DSP}1:\underset{u,\pi ,\mu }{\max }\hfill \\ J={\left(h_{\mathrm{i}\mathrm{e}}-E_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{f}}-M_{\mathrm{i}\mathrm{e}}u\right)}^{\mathrm{T}}\pi +{\left(h_{\mathrm{e}}-E_{\mathrm{e}}{y}^{\mathrm{f}}-M_{\mathrm{e}}u\right)}^{\mathrm{T}}\mu \hfill \\ \mathrm{s}.\mathrm{t}\text{.\hspace{0.17em}}\left\{\begin{array}{c}\hfill Qu\le q,\pi \le 0\hfill \\ \hfill {\left[\begin{array}{c}\hfill G_{\mathrm{i}\mathrm{e}}\hfill \\ \hfill G_{\mathrm{e}}\hfill \end{array}\right]}^{\mathrm{T}}\left(\begin{array}{c}\hfill \pi \hfill \\ \hfill \mu \hfill \end{array}\right)=b\hfill \end{array}\right.\hfill \end{array}(42)$

where π, μ are Lagrange multiplicators corresponding to inequality and equality constraints in (41). Notice that bilinear terms $u^{\mathrm{T}}{M_{\mathrm{e}}}^{\mathrm{T}}\mu$ and $u^{\mathrm{T}}{M_{\mathrm{i}\mathrm{e}}}^{\mathrm{T}}\pi$ are incorporated in objective function of DSP1, which hinder DSP1 from being directly solved due to nonconvexity. Fortunately, by checking ARED model, it can be deduced that:

$\begin{array}{ll}\hfill u^{\mathrm{T}}{M_{\mathrm{e}}}^{\mathrm{T}}\mu & =\left(\begin{array}{cc}\hfill {\left({\theta }^{+}\right)}^{\mathrm{T}}\hfill & \hfill {\left({\theta }^{-}\right)}^{\mathrm{T}}\hfill \end{array}\right)\left(\begin{array}{cc}\hfill O\hfill & \hfill -{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}\hfill \\ \hfill O\hfill & \hfill {\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\mu }_{\neg \mathrm{w}\mathrm{g}}\hfill \\ \hfill {\mu }_{\mathrm{w}\mathrm{g}}\hfill \end{array}\right)\hfill \\ \hfill & =-{\left({\theta }^{+}\right)}^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{\mu }_{\mathrm{w}\mathrm{g}}+{\left({\theta }^{-}\right)}^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{\mu }_{\mathrm{w}\mathrm{g}}\hfill \end{array}(43)$

$u^{\mathrm{T}}{M_{\mathrm{i}\mathrm{e}}}^{\mathrm{T}}\pi =0(44)$

where Λ_wg = diag (Δw_g), μ_wg is the multiplicator related with equality constraint (45) synthetized from equations in (35) and 38 by eliminating intermediate variable $w_{\mathrm{g}}^{\mathrm{g}}$, and μ_¬wg is the rest multiplicator of μ.

$p_{\mathrm{g},j}+p_{\mathrm{g},j}^{\mathrm{\Delta }}=w_{\mathrm{g},j}^{\mathrm{e}}+\left({\theta }_{\mathrm{g},j}^{+}-{\theta }_{\mathrm{g},j}^{-}\right)\mathrm{\Delta }{w}_{\mathrm{g},j}-w_{\mathrm{g},j}^{\mathrm{c}\mathrm{r}}j\in W(45)$

According to the work of Shao et al. (2017), max-min RO always derives its optimal solution at one vertex of the uncertainty set, so u can be regarded as a binary vector without loss of optimality. Considering this, we can furtherly eliminate bilinear terms in (43) by introducing auxiliary variable z⁺ = θ⁺ ◦ μ_wg, z⁻ = θ⁻ ◦ μ_wg and applying big-M method:

$\left\{\begin{array}{c}{u}^{\mathrm{T}}{M_{\mathrm{e}}}^{\mathrm{T}}\mu =-1^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{z}^{+}+1^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{z}^{-}\hfill \\ {\mu }_{\mathrm{w}\mathrm{g}}+M\left({\theta }^{+}-1\right)\le z^{+}\le M{\theta }^{+}\hfill \\ -M{\theta }^{+}\le z^{+}\le {\mu }_{\mathrm{w}\mathrm{g}}-M\left({\theta }^{+}-1\right)\hfill \\ {\mu }_{\mathrm{w}\mathrm{g}}+M\left({\theta }^{-}-1\right)\le z^{-}\le M{\theta }^{-}\hfill \\ -M{\theta }^{-}\le z^{-}\le {\mu }_{\mathrm{w}\mathrm{g}}-M\left({\theta }^{-}-1\right)\hfill \end{array}\right.(46)$

where M is a big constant, ◦ stands for the Schur-Hadamard product of two matrices sharing the same size. Finally, DSP1 can be exactly linearized to an equivalent mixed integer linear programming (MILP) problem DSP2, which can be directly solved by commercial solvers:

$\begin{array}{c}\text{DSP}2:\underset{u,\pi ,\mu ,z^{+},z^{-}}{\max }J=\hfill \\ {\left(h_{\mathrm{i}\mathrm{e}}-E_{\mathrm{i}\mathrm{e}}{y}^{\mathrm{f}}\right)}^{\mathrm{T}}\pi +{\left(h_{\mathrm{e}}-E_{\mathrm{e}}{y}^{\mathrm{f}}\right)}^{\mathrm{T}}\mu -1^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{z}^{+}+1^{\mathrm{T}}{\mathrm{\Lambda }}_{\mathrm{w}\mathrm{g}}{z}^{-}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}Qu\le q,\pi \le 0\hfill \\ {\left[\begin{array}{c}\hfill G_{\mathrm{i}\mathrm{e}}\hfill \\ \hfill G_{\mathrm{e}}\hfill \end{array}\right]}^{\mathrm{T}}\left(\begin{array}{c}\hfill \pi \hfill \\ \hfill \mu \hfill \end{array}\right)=b\hfill \\ {\mu }_{\mathrm{w}\mathrm{g}}+M\left({\theta }^{+}-1\right)\le z^{+}\le M{\theta }^{+}\hfill \\ -M{\theta }^{+}\le z^{+}\le {\mu }_{\mathrm{w}\mathrm{g}}-M\left({\theta }^{+}-1\right)\hfill \\ {\mu }_{\mathrm{w}\mathrm{g}}+M\left({\theta }^{-}-1\right)\le z^{-}\le M{\theta }^{-}\hfill \\ -M{\theta }^{-}\le z^{-}\le {\mu }_{\mathrm{w}\mathrm{g}}-M\left({\theta }^{-}-1\right)\hfill \end{array}\right.\hfill \end{array}(47)$

To sum up, we use pseudo code to demonstrate the detailed solution process of ARED based on the C&CG algorithm.

4 Real-time control for compensation of res power fluctuation

The actual power generation of RES is revealed during real-time operation, which is very likely to deviate from the predicted value due to RES power fluctuation. Hence, DSO distributes the power deviation to each generation unit and IAC cluster according to a proportion of their adjustable capacity. However, unlike distributed generators, the IAC cluster is unable to constantly generate and regulate a certain amount of power by itself, and it can only provide power adjustment by affecting the user’s comfort experience (namely, by changing room temperature). Therefore, in this section, we propose a novel control strategy for IAC clusters to achieve power adjustment for compensation of RES power fluctuation without violating IAC users’ comfort preferences.

4.1 Real-time control of IAC cluster

Suppose one IAC cluster’s target power adjustment is $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\ast }(k)$. Hence, the immediate work of IAC cluster is to decompose $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\ast }(k)$ to every IAC, which turns out to find out the unified target s^∗(k + 1) for all IACs to follow. Firstly, we calculate the target (reference) power of the IAC cluster by:

$p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\ast }\left(k+1\right)=p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{\Delta },\ast }+p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}\left(k\right)(48)$

Next, according to IAC cluster’s ACC, the unified control target s^∗(k + 1) can be determined by:

$s^{\ast }\left(k+1\right)=d_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{-1}\left(p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\ast }\left(k+1\right)\right)(49)$

where $d_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{-1}$ is the inverse function of d_cluster. Intuitively, s^∗(k + 1) can derived by the intersection point of function curve p_cluster = d_cluster(s) and horizontal line $p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}=p_{\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}}^{\ast }(k+1)$.

Finally, each IAC local controller conduct autonomous control of IAC to follow the target signal s^∗(k + 1). Thus, this process can be carried out in a decentralized style. It is noteworthy that during local IAC control, to ensure accuracy and exactness of control, the sampling and control period is usually chosen to be smaller than that of determining the target signal. Therefore, a local model predictive control (MPC) optimization is developed to model and optimize the local process of tracking s^∗(k + 1) with discrete time micro step τ = τ_mpc:

$\begin{array}{c}\min J_{\mathrm{a}\mathrm{c}}={\left({\mathbf{T}}_{\mathrm{a}}-{\mathbf{T}}_{\mathrm{a}}^{\ast }\right)}^{\mathrm{T}}{\mathbf{\Lambda }}_T\left({\mathbf{T}}_{\mathrm{a}}-{\mathbf{T}}_{\mathrm{a}}^{\ast }\right)\hfill \\ +{\left({\mathbf{p}}_{\mathrm{a}\mathrm{c}}-{\mathbf{p}}_{\mathrm{a}\mathrm{c}}^{\ast }\right)}^{\mathrm{T}}{\mathbf{\Lambda }}_p\left({\mathbf{p}}_{\mathrm{a}\mathrm{c}}-{\mathbf{p}}_{\mathrm{a}\mathrm{c}}^{\ast }\right)\hfill \\ s.t.\left\{\begin{array}{c}{\mathbf{T}}_{\mathrm{a}}^{\ast }={\phi }^{-1}\left(s^{\ast }\left(k+1\right)\right)\times 1_{N_{\mathrm{m}\mathrm{p}\mathrm{c}}\times 1}\hfill \\ {\mathbf{p}}_{\mathrm{a}\mathrm{c}}^{\ast }=d_{\mathrm{a}\mathrm{c}}\left(s^{\ast }\left(k+1\right)\right)\times 1_{N_{\mathrm{m}\mathrm{p}\mathrm{c}}\times 1}\hfill \\ \text{constraints\hspace{0.17em}}\left(2\right)-\left(9\right)\text{\hspace{0.17em}hold}\hfill \end{array}\right.\hfill \end{array}(50)$

where N_mpc is the number of micro steps when applying MPC, which is also termed “prediction horizon” in some existed literature;${\mathbf{T}}_{\mathrm{a}}={[T_{\mathrm{a}}(1),T_{\mathrm{a}}(2),..T_{\mathrm{a}}(N_{\mathrm{m}\mathrm{p}\mathrm{c}})]}^{\mathrm{T}}\in {\mathbb{R}}^{N_{\mathrm{m}\mathrm{p}\mathrm{c}}}$ is the room temperature at the end of each MPC micro step, and similarly ${\mathbf{p}}_{\mathrm{a}\mathrm{c}}\in {\mathbb{R}}^{N_{\mathrm{m}\mathrm{p}\mathrm{c}}}$is IAC electrical power of each micro step; Λ_T and Λ_p are two diagonal coefficient matrix; Constraints (2)–(9) are also required to be satisfied for each micro step. During the time period before receiving the next target signal, each IAC will repeatedly execute local MPC in a rolling manner, as shown in Figure 3. After each round of MPC, only the first micro step’s control variable (namely, f_ac (1)) result is adopted by IAC controller (Chen and Liu, 2021).

FIGURE 3

FIGURE 3. Local rolling MPC of IAC controller.

4.2 Overall processes of ARED and real-time control

Now, we integrate ARED and real-time control together. The entire processes of conducting ARED and real-time control of IAC clusters are illustrated in Figure 4. The processes colored green, gray, and yellow are carried out by local IAC controllers, IAC load aggregators, and DSO respectively.

FIGURE 4

FIGURE 4. Schematic of ARED and real-time control processes.

Firstly, each IAC aggregator computes the ACC of the IAC cluster and uploads its adjustable capacity to DSO every τ_ARED = 5 min, which determines the lower and upper bound of load power adjustment in constraint (37). Based on IAC cluster adjustable capacity and predicted power generation of RES, DSO solves ARED optimization model and determines the optimal setting point for each power generation unit, which minimizes the total operation cost under the worst uncertainty scenario.

Afterward, real-time control processes are executed with a time cycle of τ_RC = 1 min to ensure the minute-level power balance of the distribution system against RES power fluctuation. Concretely speaking, IAC aggregator determine its target “unsatisfactory degree” signal s^∗ according to (49), and broadcast s^∗ to its downstream IACs.

Finally, to track the target power indicated by s^∗, each IAC will then execute rolling MPC locally and autonomously, which is conducted every τ_mpc = 15 s with coefficient Λ_T set to zero matrix and Λ_p set to identity matrix. Therefore, to minimize cost function of MPC model (50), p_ac is forced to approach target value $p_{\mathrm{a}\mathrm{c}}^{\ast }$ as quickly as possible during investigated prediction horizon, which ensures the accuracy of tracking target value. Although temperature T_a is not explicitly forced to track the target value $T_{\mathrm{a}}^{\ast }$ in objective function, its value is still located in the range of [T_min, T_max] due to s ∈ [ − 1, 1], which guarantees good comfort experience of IAC user.

5 Study case

5.1 Case description

A modified IEEE 33-busbar distribution network is adopted in this paper, the topology of which is shown in Figure 5. Two distributed generators and two wind plants are added to busbars 9, 19, 24, 27, and four IAC clusters with 15, 15, 30, and 15 IACs individually are also integrated into busbars 7, 8, 25, and 30. In this study case, wind plants account for 50% of the total installed capacity of the explored distribution system, and can undertake 43% of the total load at maximum power generation.

FIGURE 5

FIGURE 5. Case study: a modified IEEE 33-busbar system.

For the first stage decision of ARED, parameters of all dispatchable power generation units are listed in Table 1. For second stage decision of ARED, cost coefficients of wind curtailment and load shedding {ρ_cr, ρ_ls} are set to 40, 60$/MW respectively. Power adjustment cost coefficient ρ_ad of the transmission system is set to 45$/MW, and that of distributed generators is set to 35$/MW. IAC clusters in this paper are assumed to participate in power adjustment by incentive contracts, and the contract price (i.e., cost coefficient ρ_d) is set to 35$/MW.

TABLE 1

TABLE 1. Parameters of power generation units.

The parameters of IAC are inherited from the literature [18], and for each room equipped with IAC, the initial room temperature, initial IAC compressor operation frequency, and ambient temperature are randomly generated according to distributions in Table 2 where N, U stands for Gaussian distribution and uniform distribution. We also simplify each room’s structure to be a cuboid with length 12.5m, width 8m, and height 2.5 m.

TABLE 2

TABLE 2. Random initialization of iac operation condition and user preference parameters.

5.2 Anaysis of ARED results

5.2.1 Results of ARED

Without losing generality, we firstly demonstrate ARED results with uncertainty budget set to Δw_g,j = 0.2 MW, which are listed in Table 3. For the results of the first stage decision, the operation power of wind plants at busbar 9 and 19 are scheduled to the same values as the predicted ones, which maximize the use of wind power. The second stage decision reveals the worst uncertainty scenario: the power generation of wind plants 9 and 19 simultaneously deviate -0.2 MW from predicted values, which causes an 0.4 MW power deficiency in total. To compensate for this deviation, distributed generators 24 and 27 increase their power generation by 0.274 and 0.055 MW, and controllable IAC clusters at busbar 7, 8, 25, and 30 decrease their power demand by 0.066 MW in total as well, which undertakes 16.5% of required power adjustment capacity. Meanwhile, the 0 MW power adjustment of the substation indicates that uncertain power fluctuation of RES is completely compensated within the distribution system itself, which stabilized power exchange between the distribution system and transmission system.

TABLE 3

TABLE 3. Results of ared’s first and second stage decision.

5.2.2 Comparison of ARED and traditional ED

In traditional ED, only the first stage decision is carried out to minimize total operation cost under predicted RES power generation, which is a deterministic optimization problem. As shown in Table 4. The cost of the first stage decision ($536.82/h) is lower than that of ARED ($538.14/h). However, under the worst uncertainty scenario, ED’s cost of corrective regulation in the second stage ($16.78/h) is contrarily higher than ARED ($13.84/h), which simultaneously causes a higher total cost ($553.60/h). The results also demonstrate the mechanism of ARED: by deviating the solution of the first stage decision from that of a deterministic model (which actually causes a suboptimal solution of the first stage), a smaller cost of the second stage under the worst operation scenario can be achieved, which leads to a lower total cost as well.

TABLE 4

TABLE 4. Per-hour cost comparison of ed and ared.

5.2.3 Influence of uncertainty budget

By adjusting uncertainty budget Δw_g,j from 0 to 0.3 MW (interval 0.05 MW), ARED with different uncertainty budget are also solved to derive different robust optimal solutions under distinct conservative level. The results are shown in Table 5.

TABLE 5

TABLE 5. Cost under different uncertainty magnitude.

When Δw_g,j = 0, two-stage ARED model becomes a deterministic one with optimal first stage cost and 0 s stage cost. It is noteworthy that the optimal solution of the deterministic ARED model must be a feasible solution to other general ARED models with nonzero. Therefore, the optimal values of the objective function of other general ARED models are sure to be no less than that of the deterministic ARED model. Additionally, with the increase in the uncertainty budget, the cost of both two stages tends to increase simultaneously. This indicates that the more conservative the ARED model is, the less cost-effectiveness the results will be. So, the tradeoff between cost-effectiveness and robustness of ARED should be considered in practical situations by selecting an appropriate uncertainty budget.

5.3 Analysis of real-time control results

During real-time control, the deviation between RES actual and predicted power are distributed to dispatchable units for power compensation. Within these units, the IAC cluster is the most distinct because it works more like energy storage that is incapable of constantly generating and regulating a certain amount of power by itself. Hence, in this part, we focus on making an analysis of real-time control of IAC clusters.

5.3.1 Compensation of wind power fluctuation

Figure 6 demonstrates the power compensation effect of controllable IAC cluster at busbar 25 under the aforementioned control strategy. The yellow dotted line illustrates the fluctuation of the actual power of RES, which is sampled every 4 micro steps (τ_RC = 1 min), and indicates the latest power adjustment IAC cluster should undertake. By applying our control strategy, the target power of the IAC cluster is correspondingly calculated according to (48), as shown by the blue dotted line. Then, IACs execute autonomous and local MPC control to track a unified target “unsatisfactory degree” signal translated from target power according to (49). The actual power response of IACs is shown by red line, which is updated every micro step (τ_mpc = 15 s). It is clearly shown that the actual power of the IAC cluster perfectly follows the target power. Once the target power is updated, the actual power will start to vary towards the target value.

FIGURE 6

FIGURE 6. Compensation of wind power fluctuation.

For a more intuitive presentation, we define the residual between target power and actual power by “compensation residual”, which is presented by the purple line. The residual signal actually indicates the “unrealized” power adjustment to be undertaken by the transmission system. It can be found that the residual is a nonzero value at the beginning 1 micro-step of each τ_RC control interval (4 micro steps) due to updating of target power, but is gradually eliminated during the sequential 3 micro steps. This means power compensation can be fully achieved within the distribution system, and no more ancillary power adjustment is required in the transmission system.

5.3.2 Evolution of room temperature

The scaled adjustable capability of the IAC cluster is based on each IAC’s minor contribution of power adjustment, which changes the temperature of a room. Figure 7 demonstrates room temperature evolution of the 30 controllable IACs at busbar 25. In Figure 7A, we can find that initial values of room temperature are scattered in a wide range due to our random initialization. But these temperature values gradually converge to a narrow range after around 20 micro-steps and evolve following a similar tendency. This is the effect of adopting a unified target signal s∗, as shown in Figure 7C. It can be found that the “unsatisfied degree” of each room gradually converges to the same value. Thus, high similarity can be observed concerning the evolution of different rooms’ target temperature and actual temperature.

FIGURE 7

FIGURE 7. Evoluton of (A), (B) room temperature and (C) unsatisfactory degree.

As real-time control proceeds, the room temperature is still within the range determined by each user’s preference. For ease of explanation, we randomly choose 6 rooms and plot their individual temperature curve together with preference temperature levels in Figure 7B. It is clearly shown that after the initial 20 micro-steps, all rooms’ temperature gradually locate within their individual allowable interval indicated by T_max (red horizontal line) and T_min (yellow horizontal line), and varies around T_fit (purple horizontal line). In fact, target signal s∗ is bounded to [-1,1], as shown in Figure 7C, so room temperature will not violate user’s preference constraints.

5.3.3 Adjustable capacity of IAC cluster

ACC of IAC cluster is updated for calculation of IAC cluster’s adjustable capacity before conducting ARED, so its updating period is τ_ARED = 5 min, which equals 20 micro-steps. Figure 8 demonstrates adjustable capacity evolution of the IAC cluster at busbar 25. The upward adjustable capacity of the IAC cluster can be intuitively comprehended by the gap between maximum allowable power and actual power of IAC, and so is downward adjustable capacity similarly.

FIGURE 8

FIGURE 8. Adjustable capability of IAC cluster.

It can be observed that the maximum and minimum allowable power of IAC are time-variant, and this is the most prominent feature that distinguishes controllable IAC clusters from other dispatchable units such as distributed generators. The mechanism of this feature can be explained through ACC we proposed in Section 2: Firstly, higher room temperature will lead to larger maximum and minimum allowable power of IAC according to (12) and (13), which can be comprehended by shifting ACC curve upward along the vertical p_ac -axis; Meanwhile, “unsatisfactory degree” s is also monotonically increasing with respect to room temperature according to (10). Consequently, we can deduce that a bigger s corresponds to a higher value of the maximum and minimum allowable power of the IAC cluster. Referring to evolution of “unsatisfactory degree” s in Figure 7C, s = −0.354, 0.116. −0.163, −0.421, -0.535 at the 20th, 40th, 60th, 80th, and 100th micro-steps. Hence, the variation rule of the IAC cluster’s maximum and minimum allowable power can be explained.

5.4 Influencing factors of real-time control performance

In this part, we discuss influencing factors of real-time control from two aspects: bias of prediction error and fluctuation amplitude of RES stochastic deviation from the predicted value.

5.4.1 Bias of prediction deviation

Under normal circumstances, the deviation between the predicted value and actual value of RES power generation obeys Gaussian distribution with a mean value of 0. The aforementioned results of real-time control are all obtained based on this assumption. However, positive-definite or negative-definite deviation may occur due to incorrect RES power prediction algorithm or other system failures, which is termed “bias” in this part.

Figure 9 demonstrates real-time control results of the IAC cluster at busbar 25 with a negative-definite bias of -4.5 kW superposed to its power adjustment signal, which requires the IAC cluster to increase power demand by an extra 4.5 kW. It can be noticed that curves of “unsatisfactory degree” are shifted downwards compared with those in Figure 7C. This means the power adjustment is achieved with lower room temperature, which affects the user’s comfort experience. Meanwhile, the actual power of the IAC cluster is shifted upward and reaches the maximum allowable power boundary more frequently. On this occasion, the upward capacity of the IAC cluster is squeezed to 0, and the power adjustment signal may not be fully compensated, as shown at 60th and 90th micro-step in Figure 9B.

FIGURE 9

FIGURE 9. Results of real-time control with negative biased RES prediction error. (A) Adjustable capacity evolution. (B) Power compensation results. (C) Unsatisfactory evolution.

Consequently, the bias of prediction deviation will lead to orientated increasing or decreasing of IAC cluster power, which also shifts the variation range of room temperature. The bias should be avoided as far as possible by improving the RES power prediction algorithm or eliminating system faults.

5.4.2 Amplitude of fluctuation

In this part, we amplify the power adjustment signal by 1.88 times to simulate violent RES power fluctuation. The results are calculated and shown in Figure 10.

FIGURE 10

FIGURE 10. Results of real-time control with large power adjustment signal. (A) Adjustable capacity evolution. (B) Power compensation results. (C) Unsatisfactory evolution.

Comparing Figure 10A with Figure 8, we can find that the actual power of IAC will reach maximum or minimum allowable boundary much more frequently under large fluctuation. When this happens, the upward capacity or downward capacity of the IAC cluster will be alternately squeezed to 0, and the power adjustment signal may not be fully compensated, as shown at the 55th, 70th, and 90th micro-steps in Figure 10B. For room temperature or “unsatisfactory degree” in Figure 10C, it is not shifted like the results under the biased situation of RES prediction, but the variation range becomes relatively wider compared with Figure 7C.

In fact, the actual power fluctuation of RES is usually bounded within a maximum range of 20% deviation from the predicted value. Therefore, the power adjustment signal can be fully compensated by the IAC cluster by aggregating a sufficient number of IACs or selecting the appropriate coefficient of distribution proportion.

6 Conclusion

This work proposes an adaptive robust economic dispatch (ARED) and real-time control scheme of the distribution system, which comprehensively takes the uncertainty of RES and flexibility of controllable IAC clusters into consideration. Firstly, the concept of the adjustable capacity curve (ACC) is developed to quantify the exact adjustable capacity of the IAC cluster within a certain time period and guide each IAC to carry out decentralized local model predictive control. Afterward, ARED model for the power dispatch stage is formulated, which incorporates the adjustable capacity of IAC and uncertainty set of RES power generation in the constraints. An iterative solution methodology of ARED is also designed by drawing the experience of the C&CG algorithm, of which the master problem is formulated as a QCQP and the max-min sub-problem is reformulated to a MILP via linear duality theory and big-M method. Finally, to hedge against stochastic RES power fluctuation, a real-time control strategy is also proposed, where IAC clusters provide ancillary power adjustment service based on decentralized control.

A case study on a modified IEEE 33-bus system verifies that ARED model can better guarantee the economic efficiency of the distribution system against the worst uncertainty scenario under different uncertainty budgets. Meanwhile, the real-time control strategy of the IAC cluster is also proved to be effective in eliminating RES power fluctuation without violating the comfort preferences of IAC users. In addition, the influence of biased prediction error and fluctuation amplitude is also analyzed in detail, which can provide a reference for enhancing the operation stability of the IAC cluster under our control strategy. In the future, we will furtherly investigate adjustable capability quantification and real-time control of IAC clusters based on price-based demand response.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

GC: modeling, simulation, and paper writing; DL: conceptualization and draft-reviewing.

Funding

This work was supported in part by the National Natural Science Foundation of China—Key Program of Joint Fund in Smart Grid (U2166210) and the National Key Research and Development Program of China (Basic Research Class 2017YFB0903000).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Chen, G., and Liu, D. (2021). Hybrid flow model of cyber physical distribution network and an instantiated decentralized control application. CSEE J. Power Energy Syst. (Early Access). doi:10.17775/CSEEJPES.2020.05770

CrossRef Full Text | Google Scholar

Chen, Y., Wang, Y., Kirschen, D., and Zhang, B. (2018). Model-free renewable scenario generation using generative adversarial networks. IEEE Trans. Power Syst. 33, 3265–3275. doi:10.1109/TPWRS.2018.2794541

CrossRef Full Text | Google Scholar

Chen, Z., Guo, C., Dong, S., Ding, Y., and Mao, H. (2021). Distributed robust dynamic economic dispatch of integrated transmission and distribution systems. IEEE Trans. Ind. Appl. 57, 4500–4512. doi:10.1109/TIA.2021.3091663

CrossRef Full Text | Google Scholar

Dai, Y., Kang, Y., and Xiong, X. (2017). Energy development and carbon emission scenarios towards 2050 energy transition and low carbon development roadmap for China. Beijing: China Environmental Publishing House.

Google Scholar

Ding, Y., Xie, D., Hui, H., Xu, Y., and Siano, P. (2021). Game-theoretic demand side management of thermostatically controlled loads for smoothing tie-line power of microgrids. IEEE Trans. Power Syst. 36, 4089–4101. doi:10.1109/TPWRS.2021.3065097

CrossRef Full Text | Google Scholar

Farivar, M., and Low, S. H. (2013). Branch flow model: Relaxations and convexification—Part i. IEEE Trans. Power Syst. 28, 2554–2564. doi:10.1109/TPWRS.2013.2255317

CrossRef Full Text | Google Scholar

Gao, H., Liu, J., Wang, L., and Wei, Z. (2018). Decentralized energy management for networked microgrids in future distribution systems. IEEE Trans. Power Syst. 33, 3599–3610. doi:10.1109/TPWRS.2017.2773070

CrossRef Full Text | Google Scholar

Hui, H. X., Ding, Y., and Zheng, M. L. (2019). Equivalent modeling of inverter air conditioners for providing frequency regulation service. IEEE Trans. Ind. Electron. 66, 1413–1423. doi:10.1109/tie.2018.2831192

CrossRef Full Text | Google Scholar

Jiang, A., and Wei, H. (2018). Distributed cooperation model and optimal control strategy for interaction between large-scale air conditioning and power grid based on communication. Proc. Chin. Soc. Electr. Eng. 38, 6276–6283. doi:10.13334/j.0258-8013.pcsee.172036

CrossRef Full Text | Google Scholar

Kalantar-Neyestanaki, M., Sossan, F., Bozorg, M., and Cherkaoui, R. (2020). Characterizing the reserve provision capability area of active distribution networks: A linear robust optimization method. IEEE Trans. Smart Grid 11, 2464–2475. doi:10.1109/TSG.2019.2956152

CrossRef Full Text | Google Scholar

Li, Z., Shahidehpour, M., Wu, W., Zeng, B., Zhang, B., and Zheng, W. (2015). Decentralized multiarea robust generation unit and tie-line scheduling under wind power uncertainty. IEEE Trans. Sustain. Energy 6, 1377–1388. doi:10.1109/TSTE.2015.2437273

CrossRef Full Text | Google Scholar

Liu, Y., Guo, L., and Wang, C. (2018). Economic dispatch of microgrid based on two stage robust optimization. Proc. CSEE 38, 4013–4022. doi:10.13334/j.0258-8013.pcsee.170500

CrossRef Full Text | Google Scholar

Lu, Z., Huang, H., Shan, B., Wang, Y., Du, S., and Li, J. (2017). Morphological evolution model and power forecasting prospect of future electric power systems with high proportion of renewable energy. Automation Electr. Power Syst. 41, 12–18. doi:10.7500/AEPS20170109001

CrossRef Full Text | Google Scholar

Ning, C., and You, F. (2018). Data-driven stochastic robust optimization: General computational framework and algorithm leveraging machine learning for optimization under uncertainty in the big data era. Comput. Chem. Eng. 111, 115–133. doi:10.1016/j.compchemeng.2017.12.015

CrossRef Full Text | Google Scholar

Roldán, C., Mínguez, R., García-Bertrand, R., and Arroyo, J. M. (2019). Robust transmission network expansion planning under correlated uncertainty. IEEE Trans. Power Syst. 34, 2071–2082. doi:10.1109/TPWRS.2018.2889032

CrossRef Full Text | Google Scholar

Ruan, H., Gao, H., Liu, J., and Huang, Z. (2019). A distributionally robust reactive power optimization model for active distribution network considering reactive power support of dg and switch reconfiguration. Proceedings Chin. Soc. Electr. Eng. 39, 685–695. doi:10.13334/j.0258-8013.pcsee.180241

CrossRef Full Text | Google Scholar

Shao, C., Wang, X., Shahidehpour, M., Wang, X., and Wang, B. (2017). Security-constrained unit commitment with flexible uncertainty set for variable wind power. IEEE Trans. Sustain. Energy 8, 1237–1246. doi:10.1109/TSTE.2017.2673120

CrossRef Full Text | Google Scholar

Shuai, H., Fang, J., Ai, X., Tang, Y., Wen, J., and He, H. (2019). Stochastic optimization of economic dispatch for microgrid based on approximate dynamic programming. IEEE Trans. Smart Grid 10, 2440–2452. doi:10.1109/TSG.2018.2798039

CrossRef Full Text | Google Scholar

Tan, Z., Zhong, H., Xia, Q., Kang, C., Wang, X. S., and Tang, H. (2020). Estimating the robust p-q capability of a technical virtual power plant under uncertainties. IEEE Trans. Power Syst. 35, 4285–4296. doi:10.1109/TPWRS.2020.2988069

CrossRef Full Text | Google Scholar

Wei, W., Wang, D., Jia, H., Wang, R., Guo, B., Qu, B., et al. (2016). A hierarchical and distributed control strategy of thermostatically controlled appliances for city park based on load model prediction. Proc. Chin. Soc. Electr. Eng. 36, 2049–2056. doi:10.13334/j.0258-8013.pcsee.2016.08.003

CrossRef Full Text | Google Scholar

Yang, Y., Wu, W., Wang, B., and Li, M. (2021). Chance-constrained economic dispatch considering curtailment strategy of renewable energy. IEEE Trans. Power Syst. 36, 5792–5802. doi:10.1109/TPWRS.2021.3081120

CrossRef Full Text | Google Scholar

Yao, Y., and Zhang, P. (2018). Coordinated control method for ancillary services of power system with participation of large-scale inverter air-conditioner. Automation Electr. Power Syst. 42, 127–134. doi:10.7500/AEPS20180131001

CrossRef Full Text | Google Scholar

Yao, Y., and Zhang, P. (2020). Unified control strategy of heterogeneous thermostatically controlled loads with market-based mechanism. J. Mod. Power Syst. Clean Energy 8, 1231–1239. doi:10.35833/MPCE.2018.000828

CrossRef Full Text | Google Scholar

Ye, L., Li, Z., Sun, B., Tang, Y., Lan, H., Wu, L., et al. (2018). Optimal dispatch of system integrated wind farm clusters based on stochastic model predictive control considering temporal correlation of wind power. Proc. Chin. Soc. Electr. Eng. 38, 3172–3183. doi:10.13334/j.0258-8013.pcsee.171162

CrossRef Full Text | Google Scholar

Zeng, B., and Zhao, L. (2013). Solving two-stage robust optimization problems using a column-and-constraint generation method. Operations Res. Lett. 41, 457–461. doi:10.1016/j.orl.2013.05.003

CrossRef Full Text | Google Scholar

Zhang, Z., Chen, Y., Liu, X., and Wang, W. (2019). Two-stage robust security-constrained unit commitment model considering time autocorrelation of wind/load prediction error and outage contingency probability of units. IEEE Access 7, 25398–25408. doi:10.1109/ACCESS.2019.2900254

CrossRef Full Text | Google Scholar

Zhou, A., Yang, M., Zhai, H., Zhao, Y., Zhang, C., and Wei, S. (2018). Distributionally robust real-time dispatch considering moment uncertainty of wind generation. Proc. Chin. Soc. Electr. Eng. 38, 5937–5946. doi:10.13334/j.0258-8013.pcsee.180735

CrossRef Full Text | Google Scholar