Flexible region aggregation of adjustable loads via an adaptive convex hull strategy☆

0 Introduction

The increasing penetration of variable and uncertain renewable energy sources such as wind and solar power has resulted in a growing need for flexibility in power systems[1].The utilization of additional flexibility sources has become an urgent requirement for reliable and economical scheduling of future power systems [2].The regulation potential of diverse types of loads has attracted considerable attention because the capacity for flexible loads is likely to attain a large scale [3,4].

In this new context, system-level power scheduling should incorporate the power consumed by flexible loads as decision variables.Quantifying the flexible region of load aggregation (i.e., the feasible region of the aggregate load power) is a fundamental but crucial problem.

In many early studies concerning demand response(DR), flexible regions were characterized by a category of response types, such as shiftable and sheddable loads[5].This is because in the preliminary development stage,the majority of DR participants are large industrial users,and flexible regions can be customized individually.

However, the sources of flexible loads would become significantly more diverse and smaller in the future, e.g.,thermostatically controlled loads(TCLs)and electric vehicles (EVs) on the commercial and residential sides [6,7].Different types of flexible sources are aggregated by load aggregators to attain the capacity limit of the electricity market.Thus, the power system scheduling of the system operator relies on the aggregated flexible region provided by the load aggregators.For each point within the aggregated flexible region,in terms of active power,there should be at least one control decision for each flexible load such that all the operational constraints are satisfied.Hence,previous studies that focused on the aggregated flexible region modeling of a certain category of flexible loads[6,8] are unsuitable for general practical applications.

A typical method to assess the aggregate flexibility of heterogeneous loads is to use the up-down method.It involves constructing an external characteristic model and identifying parameters using historical data [9,10].However, for many types of loads such as TCLs and EVs,the current adjustable potential is related to their previous states or cumulative internal states.Aggregation cannot be performed by increasing the external characteristic indices such as the adjustable capacity.Otherwise,the estimated flexibility may exceed the attainable limits[11,12].In addition, numerous DR samples are required to enable the data-mining process.This may be difficult to achieve at present [13].

Recently,new bottom-up modeling techniques based on computational geometry have been introduced.The operational domains of the flexible loads are modeled using geometric feasible sets.The sets are then aggregated via the Minkowski sum of the feasible sets.Thus, successive internal-state dependencies can be captured.However,calculating the Minkowski sum of two polytopic sets specified by facets is an NP-hard problem [15].Hence, the approximation of the Minkowski sum is currently the only feasible solution for practical applications.

Different types of specific polytopes have been evaluated to approximate actual feasible sets of individual loads.This is because their Minkowski sums can be calculated accurately in polynomial time.Boxes[14],zonotopes[15], and homothets [16] were popular options in previous studies.However, shape features should be specified beforehand according to the parameters of individual polytopes [17-19].This results in sensitivity to the types of flexible loads because we cannot anticipate simultaneous high-accuracy performance for heterogeneous loads.

Inspired by Ref.[20],which used the union of homothets to approximate the true aggregated flexible region,this study proposes an approximation method based on the convex hull of multilayer maximum volume boxes and integrates the result into a system-level day-ahead scheduling problem.

The difference between this work and [20] is that a greedy strategy is adopted.Herein, maximum volume boxes rather than presupposed homothets are used, and a convex hull generation algorithm for high-dimensional cases is proposed.In addition, the aggregation result is applied to the power system scheduling optimization problem, and the practical advantages are measured.

The highlights of this study include the following:

1) A heterogeneity-adaptive and high-accuracy method is proposed for aggregate flexibility evaluation.

2) The flexibility of heterogeneous loads at the same node can be utilized fully in day-ahead scheduling with high reliability and economic performance.

The remainder of this paper is organized as follows.In Section 1, the general framework is presented.Section 2 describes the aggregation algorithm in detail.Section 3 introduces the day-ahead optimal scheduling model that considers the aggregate flexibility from heterogeneous loads.Section 4 presents a disaggregation model for the load aggregator.Section 5 presents case studies.Finally,the conclusions are presented in Section 6.

1 General framework

This study is based on the premise that different types of loads can actively participate in the power balance as flexible resources through aggregation.Load aggregators are responsible for managing diverse loads at the load nodes.The states of the individual loads are visible to the load aggregators but not to the system operator.The load aggregators should announce the feasible regions of the aggregate power (i.e., the aggregated flexible region) to the system operator.Thereby, the system operator can schedule from the system level and dispatch the adjustable resources of the power system.This process is illustrated in Fig.1.

To achieve this goal, a load aggregator should aggregate the individual flexible regions of heterogeneous loads into an aggregated flexible region.The individual flexible region of each load is first characterized accurately according to its physical operation requirements.As an illustration, flexible regions considering the power consumption of the first and second time steps (P1 and P2, respectively)are shown.The heterogeneity of the loads is reflected in the differences in the shapes of the individual flexible regions.The aggregated flexible region is then calculated approximately using the convex hull of the multilayer maximum volume boxes based on the Minkowski sum.

The aggregated flexible region affects the feasible region of the scheduling optimization problem of the system operator.This subsequently affects the quality of the scheduling solution.

After the power dispatch results from the system operator are obtained, the load aggregators perform disaggregation.That is, these determine the power-consuming modes of all their loads.

1.1 Polyhedron representation of individual flexible regions

The power characteristics of the predominant types of flexible loads can be described linearly.For the constraints in the differential form,such as the constraints on the ramp rate, different operations can be performed by selecting a suitable step size.The constraints in the integral form,such as the constraints on cumulative energy,can be reformulated via discretization.For sources such as battery storage systems with nonzero (dis)charging efficiency, an approximate linear model [21] should be applied to adapt the following approach.Their flexible regions can be expressed using polytopes, as shown in (1):

where D denotes the dimensions of power vector p.Because it physically represents the number of time steps,we refer to it as the time dimension in this study.

To clarify this, we infer the polyhedral representations of EVs and TCLs based on their widely adopted physical models.The polyhedral representations of other types of loads can be derived similarly.

(1) EVs

The physical constraints of EVs include the upper and lower bounds of charging power (as shown in (2)) and charged energy (as shown in (3)):

where the superscript EV is used to specify the EVs; min and max are applied to denote the lower and upper bounds, respectively; the subscript i is used to distinguish the EVs; and t represents the serial number of the time steps.p represents the active power,E is the storage energy of the battery,and is the initial energy of the battery.k is the energy limit coefficient, and Δt is the time duration between the time steps.

Combining (2) and (3), the polyhedron representation of the EVs PEV can be derived as shown in (4)-(6)

whereand 1D are D-dimensional column vectors with all elements being one, ID is a Ddimensional identity matrix, and LD is a D-dimensional unit lower-triangular matrix.

For visualization purposes and to address the concerns of the system operator, a 2-D feasible space is shown in Fig.2.The two sequential active power variables [P(T)and P(T+1)]of the individual EVare on the horizontal and vertical axes, respectively.

Fig.2(a)shows the adjustable ranges of EV power during adjacent periods [T, T + 1] with charged energy constraints.For example, as shown in Fig.2(b), if the discharging power P(T) has previously reduced the SOC of the EV battery, the adjustable range of P(T + 1) in the next step may not expand tobecause of the lower bound of the charged energy.We can observe a truncation in the bottom-left corner of the rectangle, and vice versa.In Fig.2(c),the upper-right corner of the feasible rectangle is truncated owing to the upper bound of the charged energy if the charging power P(T) increases the SOC of the EV battery.

(2) TCLs

The physical constraints of the TCLs include the upper and lower bounds of power consumption(as shown in(7))and ambient temperature (as shown in (8)):

where the superscript TCL specifies the TCLs,θ represents the ambient temperature, and q is the temperature limit parameter.

The relationship between the ambient temperature and power consumption can be derived based on the first-order dynamic equation of ambient temperature.This is shown in (9):

where CTCL denotes the thermal capacity,RTCL denotes the thermal resistance, θTCL,out denotes the outdoor temperature vector, and ηTCL denotes the efficiency vector.

The full feasible space of the TCL considering the temperature variables θ and power variables p is highdimensional.Because only the adjustable active power variables are considered by the system operators,the variables θ are eliminated, and the polyhedron representation of the TCLs PTCL can be derived as shown in (10)-(18).The elimination of variables corresponds to the projection of the initial feasible region onto the power-variable space in a lower-dimensional form.

where 0 is an all-zero column vector.

The two-dimensional feasible space of the TCL is shown in Fig.3(a).Without temperature constraints, the feasible region space formed by the adjustable poweron [P(T), P(T + 1)] should be a square with upper and lower boundaries of and,respectively.The square is truncated owing to temperature control constraints.This is shown in Fig.3(b) and (c).

Intuitively, the individual flexible regions can vary for different types of loads and operating states of the same load category.Moreover, because of the existence of temporal coupling relationships,the feasible regions vary from boxes to complex shapes.

1.2 Mathematical formulation of the aggregated flexible region based on the Minkowski sum

The target of aggregation is to obtain a feasible region of the aggregate load power at the same node.

Considering a population of flexible loads indexed by j ∈J:= {1,...,J} (where the individual flexible region of each load is denoted as Pj),the purpose of the aggregation is to obtain an aggregated flexible-region Pagg that satisfies(19):

where pagg is the column vector of the aggregate load power.

Eq.(19)exactly fits the definition of the Minkowski sum[22].That is, mathematically, we need to calculate (21).

where ⊕is the Minkowski sum operator.

The Minkowski sum of the polytopes remains a polytope.That is, Pagg can also be expressed as shown in (22):

The Minkowski sum of any two polytopes can be derived by traversing all feasible facet combinations and calculating the Minkowski sum[23].The accurate calculation of the Minkowski sum based on facets is usually an NP-hard problem.Several special cases are exceptions,e.g., when the polytopes are all boxes, zonotopes, or homothets.However, special cases are typically ineffective for practical applications.Therefore, approximating the results of the Minkowski sum is an affordable engineering approach.The focus is on improving the accuracy of the approximation.

2 Aggregation using the convex hull of multilayer boxes via load aggregators

2.1 Overall flowchart

The proposed method is characterized by the union of multilayer maximum volume boxes to internally approximate individual flexible regions.It enables a convex polytopic representation of the aggregated flexible region over Minkowski-summed boxes.The aggregation method comprises three steps.This is graphically depicted in Fig.4.

Step 1: Decomposition of individual flexible regions.The goal is to use a series of boxes to approach the true flexible region.Fig.4 shows the individual flexible regions for loads i and j.Both these are decomposed into a union of two-layer boxes as an example.

Step 2: Minkowski sum of boxes.This process calculates the Minkowski sum over boxes along the same axial direction.It outputs a union of boxes as an approximate result of an aggregated flexible region.

Step 3:Facet-growing algorithm(FGA)for convex-hull approximation.The objective is to obtain a convex representation of the aggregated flexible region because it is the basis for optimization applications.

Using multilayer maximum volume boxes, highaccuracy representations can be achieved for different types of flexible loads.Based on the properties of the Minkowski sum of boxes,the Minkowski sum of a large collection of loads can be computed straightforwardly with linear complexity.Its output is the union of a set of boxes.The FGA for convex hull computation further enables the convex polyhedron representation of the aggregated flexible region.Thus,it can be integrated into scheduling optimization problems.The detailed formulations are as follows.

2.2 Step 1: Decomposition of individual flexible regions

Unlike previous methods in which only one polytope is used to internally approximate an individual flexible polytope, this method applies multiple layers of maximum volume boxes.This is because true flexibility polytopes can be approached with higher accuracy when these are not in a box form.

During the stacking process, the layers of boxes may overlap with each other.However, the overlap does not affect the calculation of the Minkowski sum.

The box used in this study was formulated as shown in(22) [24].

where B represents the box,p is the consumption power vector,c is the geometric center vector of the box,andis the half-edge length vector.

A box can be rewritten as shown in(25).That is,a box describes the upper and lower power limits of each time step.

Decomposition determines a series of boxes to jointly approximate an individual flexible region infinitely.This process involves the concept of layers.

For an individual flexible region, we can determine the maximum volume box inside.This is indicated by the blue rectangle in Fig.4.This box is defined as the first-layer box.It divides the D-dimensional individual flexible regions into 2D subregions.We can also determine the respective maximum volume boxes within these subregions.This is illustrated by the green rectangles in Fig.4.These are denoted as the second-layer boxes.In a similar fashion,there are(2D)m-1 boxes for the m-th layer.We denote the n-th box of the m-th layer as and use （·）j to distinguish the j-th load.

Eq.(26) should be solved times in the worst cases.Therefore,it is important to reduce the computation time.The following skills are applied to achieve this: (a)The matrix A in the linear constraint is sparse.Therefore,we can employ sparse algorithms.(b) To apply this method applied for day-ahead scheduling, it is necessary to limit the number of time dimensions D for controllable loads to at most four.

2.3 Step 2: Minkowski sum of boxes

The Minkowski sum of boxes is still a box.It can be calculated rapidly based on the addition operation.

Consider a population of boxes indexed by j ∈J:= {1,...,J}.Its geometric centers are denoted as cj.The half-edge length is .The Minkowski sum results of these parameters were computed as shown in(25)and(26):

According to the distributive law of Minkowski sum[25], the Minkowski sum of the union of boxes equals the union of the Minkowski-summed boxes.Hence, we first calculated the Minkowski sum and outputted the union result for the following convex hull generation.

The Minkowski sum of boxes should be performed on the Cartesian product of the decomposed boxes of the individual flexible regions.As the layer number m increases, the number of elements in the Cartesian products increases at an exceptionally high rate.Considering that the subsequent steps generate a convex hull over the union of the Minkowski-summed boxes, discarding boxes with internal redundancies has a negligible effect on the final result.

Therefore,Step 2 performs only the Minkowski sum for boxes with different flexible loads in the same layer and along the same direction.This is shown in (27):

where is the Minkowski sum of boxes in the same layer and in the same direction for all flexible loads.It is still a box.The serial number n is assumed to be equal for the same directional box of the same layer m under different flexible loads.

Considering Fig.4 as an example, for the twodimensional case, only the two boxes along the +y, +x,-y,and-x axes and that of the first layer are Minkowski summed.

The final output of Step 2 is.

2.4 Step 3: Facet-growing algorithm (FGA) for convex hull approximation

This step was performed to obtain the convex hull of the union output in Step 2.The basic concept of FGA is to use the vertices of the boxes farthest along the coordinate axes to generate an initial approximate convex hull and then expand it outward step-by-step based on the facets.Pk denotes the results after facet growth was conducted k times.Apparently, represents the growth from Pk to Pk+1, as shown in Fig.4.

A flowchart of the FGA is presented in Fig.5.The FGA contains the following processes: First, the farthest vertices outside the facets of Pk are determined.Second,the facets that need to be replaced are determined.Third,is generated using these farthest vertices and the replaced facets.Finally, Pk+1 is obtained by combining Pk and .The details are as follows:

1) Identify the farthest vertices outside the facets of Pk.Assuming that the equation of a facet of Pk is λTx=τ, identifying the farthest vertex outside the facet is equivalent to the optimization problem shown in (28).

where V is the vertex set of all the boxes derived from Step 2.

2) The facets that need to be replaced are to be determined.Assuming that the farthest vertex of a facet=υi is x*, all the facets of Pk are examined to determine if these should be replaced.For a given facet with the function =υj, it should be replaced if cj =-υj ＞0 is true.

3) is generated using the farthest vertices and their replaced facets.The facets that are indicated to be replaced are used to generate an incremental convex hull in conjunction with the farthest vertices.

4) Pk+1 is obtained by combining Pk and.

The procedure stops after the given generating time limit K is attained.

The initialization of the convex hull and subsequent incremental convex hull generation at each step were implemented by calling the Qhull toolkit [26].Because Qhull is based on the fast convex hull algorithm,the complexity of solving the convex hull of all the vertices obtained in Step 2 is O （|V|D）.Here, |V| represents the number of vertices, and D is the time dimension.When partial vertices V small are considered to generate a convex hull, the algorithm complexity becomes .

2.5 Theoretical analyses of algorithm performance

1) Conservatism of the algorithm

Conservatism implies that the obtained flexible region is a stringent inner approximation of the actual aggregate flexibility polytope.

Proof.Assume that there are J flexible loads.These are distinguished by the subscript j.

According to Step 1,

From Step 2,

After Step 3,

This implies that the optimization solution based on the approximated flexible region is always feasible.

2) Approximate accuracy assurance for the algorithm

The approximate accuracy assurance is that the volume proportion of the true aggregated flexible region occupied by the approximate result rv continues to increase provided that rv ＜1 has an increasing number of layers m.When m →∞, rv →1.

The proof can be obtained from the literature [20].

3) Computational efficiency analysis

For Step 1,when calculating the maximum boxes of the m th layer, model (24) should be solved （2D ）m-1 times in the worst-case scenario.

For Step 2, the Minkowski sum should be calculatedtimes.

For Step 3, there would bevertices in the worst case, and a D-dimensional initial convex hull should be generated using D2D vertices.Because the vertices of the boxes farthest along the coordinate axes are used for initialization, it is independent of m.

An increase in the time dimension caused an abrupt increase in the computational load.Rolling calculations and parallel computation strategies can be applied to analyze the demand for longer durations.

3 Day-ahead optimal scheduling by system operators

The aggregation method can be applied in many cases.In this study,we used it in a day-ahead optimal scheduling problem.The aim was to fully utilize the flexibility of small-capacity loads.

3.1 Objective

The objective of the system operator is to minimize the overall operational cost co.This includes the thermal generation cost ch, load-loss penalty punishment closs, and renewable energy abandonment cost ca.The quadratic formulation of ch is adopted from [28].

where Ωh,Ωb,and Ωw are the sets of thermal plants,buses,and wind farms, respectively; pht is the thermal power at time step t; , and are the cost coefficients for the thermal generation cost;is the load loss power at time step t; aloss is the load loss; is the abandoned power of renewable energy at time step t;and aa is the abandonment cost coefficient.

3.2 Constraints

The constraints of the model include the following:

1) Constraints of the thermal power plants

where(36)describes the bounds of the power generation of thermal plant i,represents the minimum technical output of thermal plant i, represents the maximum attainable power output of thermal plant i at time step t,(37) represents the upper bound of the maximum attainable power output, represents the rated power of thermal plant i, (38) represents the upward ramp constraint,represents the maximum upward ramp power of thermal plant i, and (39) represents the maximum downward amp power of thermal plant i.

2) Constraints of the wind power farm

where(40)describes the bounds of wind power generation, is the maximum wind power generation of the wind farm j at time step t,and(41)is the method for calculating the abandoned wind power.

3) Constraints of the aggregate flexible loads for each load bus

Here, (42) describes the constraint of the aggregate loads of bus k.Ωagg is the set that records buses with flexible loads.The symbol ‘‘～” is used to represent the approximation results.

4) Constraints of the power transmission lines

where and represent the lower and upper bounds of transmission line l, respectively; L denotes the power transfer distribution factor (PTDF); Ωl is the line set; ΩfL is the flexible load group set;and ΩrL is the set that records the buses of rigid loads.

5) Constraints of the power balance

6) Reserve capacity constraints

where St is the reserve capacity demand at time step t.

A rolling strategy was adopted to render the calculations more tractable.The 24 h period was divided into eight time steps.The nodal flexible loads were updated at the beginning of each time step.

4 Disaggregation model for the load aggregator k

The disaggregation model of the load aggregator k determines the optimal power allocation scheme for its managed loads.

4.1 Objective

The objective is to maximize the energy consumption of the TCLs, as shown in (48).This is because the total load power curve is stipulated by the system operator, and a preferential electricity supply to the TCLs can enhance the user experience.

whereis the TCL set of aggregator k and pTCL is the load power of the TCLs.

4.2 Constraints

The constraints of the model include the following:

1) Operational constraints of individual flexible loads

where is the EV set of aggregator k.The operational constraints(47)and(48)can be derived as described in Section 1.1.

2) Constraint of the total load power

5 Case studies

The proposed model was tested on a modified six-bus system and a future-generation mix from the North China Grid 2050 (NCG 2050).

1) This study verified the proposed aggregation method on a six-bus system with flexible residential loads including 20 EV aggregation groups and 28 TCL aggregation groups.The approximate and computational efficiencies were compared.

2) Further simulations were conducted on the planning system to demonstrate its practicality.In addition to the 20 groups of residential loads, flexible resources such as batteries and power-to-hydrogen (P2H) are involved in the day-ahead schedule study.

For the six-bus case, the time dimension D was set to two-four, the number of layers M was set to two-five,and the number of FGA iterations was set to two-five for comparison.For NCG 2050, the time dimension D was set to three, the number of layers M was set to four,and the number of FGA iterations was four.

Yalmip(version R20210331)[30]with CPLEX(version 12.9)[31]was used to solve the optimization problems on a laptop with a 2.90 GHz CPU and 32.0 GB RAM.

5.1 Six-bus system

The modified six-bus system is shown in Fig.6.In this system, three thermal power units are connected to buses 1, 2, and 6.The total installed capacity of the thermal power units is 420 MW.A wind farm with an installed capacity of 200 MW is connected to bus 5.The details are available in [27].

The loads are integrated from buses 3 to 5,and all these include rigid and flexible parts.Flexible loads contain TCLs and EVs and are updated every 3 h.

Considering that each TCL and EV is too small for a large-scale power system,the TCLs and EVs are first combined into groups.The loads in one group are assumed to work in conjunction and have identical internal states.However, the initial states and parameters of the different groups are generated randomly, conforming to a uniform distribution.It is assumed that 18 TCL groups are integrated into bus 3, 10 TCL groups into bus 5, and 15 EV groups are connected to buses 4 and 5.The flexible load mix is presented in Table 1.In Table 1, the power and energy ratios represent the power and energy consumption proportions without DR, respectively, compared with the total electricity load demand.

The parameters for each TCL and EV are generated based on a uniform distribution assumption.The settings are shown in Tables 2 and 3.The load curves for each load bus are presented in Fig.7.

The following case studies are divided into two parts:The first part compares the accuracy and efficiency performances of different aggregation methods.The second part shows the advantages of improved accuracy in aggregation for day-ahead scheduling.

In this section, previous methods including box-based[14], zonotope-based [15], and homothet-based [16] methods and the union of homothets [20] are compared with the proposed method in terms of the aggregation accuracy and efficiency.

1) Accuracy of the proposed methods

The accuracy index can be defined based on the volume ratio of the flexible polytopes [11,15], as shown in (50).A power of 1/D was added as a normalization process for different dimensions to highlight the impact of the algorithm itself.Larger results indicate better performance.A ratio of 1 indicates that the approximation is more accurate.

where αv denotes the accuracy index.The volume of ～Pagg is derived via the Monte Carlo method, whereas that of Pagg is calculated via the MPT toolbox [29].

The statistical mean and standard deviation of the accuracy index αv calculated for different load types are shown in Fig.8.For the results of the proposed method,the number of box layers is set to four,and the FGA is applied iteratively four times.Out of consistency, four-layer results are shown for the union-of-homothets method.Because the volume of Pagg is difficult to obtain for a large collection of heterogeneous loads when the dimension is larger than four, the dimension is constrained in the range of two-four, and the number of flexible loads set to three only in this case is.Load types of only TCLs, only EVs,and TCLs and EV mixtures are computed.The parameters of the TCLs and EVs are generated randomly and conform to a uniform distribution [11].

From Fig.8(a), the mean values of the accuracy index of the proposed method are the best for the different time dimensions among the different methods.Moreover, the improvement in accuracy is more evident when the time dimension increases.It gains from the adoption of a greedy strategy using maximum volume boxes as the decomposition base and applying convex hull generation.

Further details regarding the accuracy are listed in Table 4.Because all the Minkowski sum approximation methods are inner expansion algorithms, a high accuracy implies that the optimal solution of the power system’s optimal operation approaches the real value.The proposed method can obtain the result closest to the real range.The more the layers of boxes and larger the number of executions of the FGA, the better is the approximation effect.

Furthermore,based on Fig.8(b),the standard deviation of the accuracy indices of the proposed method was the smallest among the different methods.This indicates that the proposed method is insensitive to the heterogeneity of flexible loads.

2) Day-ahead optimal scheduling results

This subsection illustrates the impact of the approximation accuracy of the feasible region on the optimization capability of the day-ahead optimal operations.Thecost-minimization problem was solved using a branchand-bound algorithm.The total operational costs of the different methods are listed in Table 5.

A comparison of the results of the one-by-one modeling(No.6)and those of the flexible loads without DR(No.7)revealed that the total operating cost was reduced effectively using the DR potential.The maximum valueattained 11.79 %.This reduction is a reference result assuming that these loads are fully visible and controllable by the system operator.However, this assumption is invalid for practical scenarios.Therefore,one-by-one modeling is not applicable in practice.The higher the accuracy of the Minkowski sum approximation method, the closer the reduction was to 11.79 %.

In the case studies, the cost reduction result of the proposed method with five-layer boxes and five-time FGA was the closest to that of one-by-one modeling and showed higher cost reduction advantages.Fig.9 visually demonstrates how the feasible region of controllable loads at node 3 gradually expands in all the methods in Table 5.

If the timing-coupling constraints are not considered(i.e.,only the upper and lower bounds of power are considered),an infeasible scheme is produced.From Fig.9,when temporal coupling constraints are overlooked,the approximated bounds are outside the reference bounds.This indicates that the response requirements cannot be satisfied in practice and that the scheduling scheme is unreliable.

3) Efficiency performance

The previous subsections demonstrate how the improved aggregation accuracy enhanced the system’s optimization capability.However, tradeoffs exist.When the Minkowski sum is used to approximate highdimensional feasible spaces, an inherent tradeoffexists between computational time and accuracy.

Although the simple system described in this subsection includes only six nodes, it incorporates 43 controllable loads (28 TCLs and 20 EVs) with different internal states.This significantly exceeds the number of variables in existing mixed-integer programming cases involving controllable loads.Table 6 demonstrates the additional computational cost associated with improving the accuracy.

Previous Minkowski sum-based methods (Nos.1-3 in Table 6) specify the topology and number of layers of the polytope in a uniform manner.Although this reduces the computational cost, it does not ensure the desired accuracy.

To achieve sufficient accuracy, a certain form of optimization should be applied to expand the volume of the boxes and increase the number of layers.The union-ofhomothets-based method (No.5) can achieve an accuracy similar to proposed method(No.6)by increasing the number of layers.However,as derived in Section 2.5,the computational complexity increases super-exponentially.This results in a two-order increase in computational cost.

5.2 The case of North China Grid 2050 (NCG 2050)

This test system was based on the future case of the North China Grid 2050.According to the plan, NCG 2050 should be a zero-carbon system with significantly high renewable penetration.The total installed power capacity is 173.78 GW.It consists of 43.2 GW of conventional thermal power units (101 units), 23.48 GW of wind power, and 107.1 GW of PV.The maximum net load (net load=load-renewables)is 50.97 GW,and the minimum net load is -6.5 GW.This results in high PV penetration.The parameters of the flexible load used in this case are the amplification of those in the six-bus system according to the capacity of the flexible loads.The wind and solar power generation data were obtained from actual operational data for a specific historical year.The details are presented in Table 7.

Studies[6-8,27]have provided optimal operation methods for high-penetration renewable energy power systems considering controllable loads.These studies focused primarily on the contribution of controllable loads in enhancing system flexibility, whereas this case emphasizes the contribution of controllable loads in improving the system’s power balance probability.In high-penetration renewable energy power systems,a small amount of power imbalance is permitted to control the overall costs.However,the estimation of the imbalance rate should be based on an accurate calculation of feasible regions.Therefore,this study aimed to improve the estimation accuracy of the feasible regions of controllable loads.

Calculations and analyses were conducted for four subcases (denoted as Cases 1-4) based on the inclusion of all adjustable resources such as sources,loads,and storage,aswell as on whether temporal coupling constraints were considered.

Case 1:All the flexible resources are included,considering temporal coupling.

Case 2: All the flexible resources are included, without considering temporal coupling.

Case 3: Only the thermal power regulation capacity is included, considering temporal coupling.

Case 4: Only the thermal power regulation capacity is included, without considering temporal coupling.

The daily operating costs for the entire year (365 days)and the power balance indicators including balanced power and unbalanced risk were calculated.These are listed in Table 8.

The balance probability reflects the probability that power balance can be achieved while considering the uncertainty of the load and renewable power generation.It is calculated using the integral of the combined probability density function of the uncertain sources over the aggregate feasible region of all regulation sources.This is shown in (53) and (54).The unbalanced risk indicates the risk of energy imbalance.It can be computed as the expectation of unbalanced energy (as shown in Eqs.(55)and (56)).

where BP and UR represent the balance probability and unbalanced risk, respectively.xsum is the value of the summed random variable Xsum of the random variables of the nn uncertainty sources, and f （·） is its combined probability density function.∁is a complementary set.

Comparing the results of Cases 1 and 2,as well as those of Cases 3 and 4 in Table 8,the estimated system operating costs are lower when the temporal coupling constraints are omitted.This difference is even more pronounced when all types of flexible resources are considered.However, the power balance probability of the system in Case 2 is significantly lower than the required 90 %.This indicates a higher risk of imbalance.This, in turn, implies that the risks of load shedding and curtailment, in conjunction with their corresponding economic costs, are not anticipated effectively during the actual operation.

However, a comparison of the results of Cases 1 and 3 in Table 8 reveals that if all flexible resources can be utilized, the operating costs of a high-proportion system can be reduced effectively while ensuring that the power balance performance satisfies the requirements.In this example, the daily operating cost decreases from RMB 1.0439 billion to RMB 0.7496 billion.This amounts to a reduction of 28.2 %.

From Fig.10,the operational records of that year indicate that the regional power grid experienced extreme weather conditions from the sixth to the eighth day.This resulted in the shutdown of many wind turbines and prediction errors of up to 70 %.This section focuses on the analysis of these three days to demonstrate the effectiveness of a flexible load.

Table 8 lists the computation time application of the proposed method in a large system.The full model case(Case 1) costs 6.5 h for the entire year (365 days).These observations are reasonable for long-term studies.Comparatively, in the classical optimization operation model[6-8,27],the controllable loads and energy storage are considered as a single controllable load or energy storage.However, the case incorporates 10 + 10 = 20 controllable loads with different internal states.This significantly exceeds the number of variables in existing mixed-integer programming cases involving controllable loads.

6 Conclusion

To improve the approximation accuracy of the aggregated flexible region for day-ahead scheduling, this paper proposes a novel aggregation method using a convex hull of multilayer maximum-volume boxes.

Compared with existing aggregation algorithms based on the Minkowski sum of special polytopes, the proposed method shows high adaptability to the heterogeneity of flexible resources.Compared with the method using the union of homothets, a smaller number of box layers is required under the same requirement of approximation accuracy,and at the same time,the output result is in convex form.Compared to the external performance of conventional methods, this algorithm can strictly guarantee the conservative property of the approximate aggregated flexible regions.

When applied in the day-ahead scheduling problem,the scheduling schemes can always benefit from the conservative property, which means that the dispatch instructions can be followed reliably.The adjustable ranges of load aggregations can be enlarged effectively; in the case study,the average value can be increased from 34.64 % to 92.76 %, and the performance can be improved further by changing the algorithm parameters.The full use of load flexibility can help decrease the operation cost at the system level.In the 6-bus case study, an 11.72 % cost can be saved compared with the nonresponse case, and the proportion is effectively increased compared with the 9.11 % value of the conventional box-based method.The NCG 2050 case shows that flexible loads improve the power balance in a zero-carbon system,especially on some extreme days.

CRediT authorship contribution statement

Yisha Lin: Writing - review & editing, Writing - original draft,Visualization,Validation,Methodology,Formal analysis,Data curation,Conceptualization.Zongxiang Lu:Writing-review&editing,Supervision,Resources,Investigation, Funding acquisition, Conceptualization.Ying Qiao:Writing-review&editing,Validation,Supervision,Methodology, Investigation, Conceptualization.Ruijie Chen: Writing - review & editing, Visualization, Methodology, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by State Grid science and technology projects ‘‘Research on energy and power supply and demand interactive simulation technology for new power system (5100-202257028A-1-1-ZN)”.