Feasible operation region of an electricity distribution network ☆

the analytical boundaries are quantified. An 11 kV radial distribution network from the United Kingdom Generic Distribution System (UKGDS) was used for the case study. The simulation results show that the quadratic analytical boundaries well approximated the real boundaries of FOR. The maximum errors for thermal and voltage boundaries would maximally cause an overcurrent up to 116 % and an undervoltage down to 0.96p.u., which are able to satisfy the requirements of engineering practice. Compared to the existing linear approximation (termed as hyperplane expressions) of FOR boundaries, the proposed quadratic expressions were proved to have higher accuracy.

• Feasible operation regions (FORs) are identified for power distribution networks.
• Novel quadratic analytical boundaries are proposed for FORs.
• A general high-dimensional error analysis approach is provided for validation.
• The proposed quadratic boundaries outperform the traditional hyperplane boundaries.

Keywords: Feasible operation region Hosting capacity Voltage boundaries Thermal boundaries High-dimensional error analysis Electricity distribution networks A B S T R A C T
In the course of net zero carbon transition, electricity distribution networks are faced with great challenges brought by increasing total and peak demand due to the electrification of transport and heat, as well as the significant uncertainties from renewable power generation and customer behaviors (such as electric vehicle travelling behaviors).Therefore, distribution network operators (DNOs) need effective tools to assess the capability of distribution networks in integrating generation and demand to conduct active management and efficient expansion of networks.However, conventional scenario-based assessment methods only provide conservative and limited information on the network capability, also with exponentially increased computational burden.By contrast, another stream of methods with different philosophy, named as operation region-based methods, describe the overall picture of the network capability analytically with little computational power needed.However, the existing linearized analytical regions, which are mainly applied to electricity transmission networks, cannot describe the capability of distribution networks accurately enough.To solve these problems, a novel feasible operation region (FOR) method with quadratic analytical expressions was proposed to characterize the range of the operating states of distribution networks, where the thermal and voltage constraints will not be violated.The FOR is a geometry in a high-dimensional space, and thus a high-dimensional error analysis approach was further developed for validating the proposed method.The boundary errors are described by multiple distance functions and operational indices, and the conservativeness of the analytical boundaries are quantified.An 11 kV radial distribution network from the United Kingdom Generic Distribution System (UKGDS) was used for the case study.The simulation results show that the quadratic analytical boundaries well approximated the real boundaries of FOR.The maximum errors for thermal and voltage boundaries would maximally cause an overcurrent up to 116 % and an undervoltage down to 0.96p.u., which are able to satisfy the requirements of engineering practice.Compared to the existing linear approximation (termed as hyperplane expressions) of FOR boundaries, the proposed quadratic expressions were proved to have higher accuracy.

Introduction
The net zero agenda worldwide requires and facilitates the development of low carbon technologies in electricity generation and supply.Large-scale low carbon distributed generation (DG) such as wind turbines and photovoltaic panels, grid-connected energy storage and new electricity demand like electric vehicles and heat pumps are being increasingly integrated into electricity distribution networks.Conventionally, "worst-case" analysis (e.g., cases with minimum and maximum loading conditions) is used for assessing the "hosting capacity" of distribution networks (i.e., the capability of distribution networks to accommodate increasing generation and/or demand) [1,2].However, the hosting capacity assessment under "worst-case" analysis method is very conservative, which can result in inefficient use of the headroom of the electricity distribution network.Moreover, the uncertainties in the potential locations and real-time power generation/consumption of DG and new demand significantly obstruct the identification of the "worst cases".Without concerning the many potential locations of DG, pernode hosting capacity of distribution networks becomes a compromised method used by some utilities [3].An alternative way is to consider the maximum total hosting capacity of all candidate locations [2,4-6].In these methods, Monte Carlo simulation methods are usually used [7][8][9].However, neither enough representative DG deployment, nor enough cases of power generation/loading conditions can be generated (even if using Monte Carlo simulation) because of the curse of dimensionality, which makes the above scenario-based methods compromise to conservative or limited results for network capability.
To overcome the deficiency of the scenario-based methods, the operation region-based methodology has been proposed to characterize the allowable range of power injections at nodes of a distribution network.The operation region-based methodology was first developed for assessing the security of electricity transmission networks [10,11] but has started to be applied to electricity distribution networks [12][13][14] and integrated energy systems [15][16][17][18][19] in recent years.It should be noted that we used the term "operation region" in this paper rather than "security region" used in those studies to avoid the confusion from "security".Although being different from the research of this paper, some basic ideas of those studies on "security region", such as the methods to depict a region in the power injections space, has enlightened us to develop the work of this paper.Regarding an electricity distribution network, the boundaries of its operation region, limited by various network constraints such as bus voltage and line flow constraints, contain the whole information of the capability of the network to accommodate more generation and demand.Within the boundaries of the operation region are all feasible operation states (i.e., allowable power injections), while outside of the boundaries some network constraints will be violated.In this regard, identifying and analyzing the region boundaries can be a more accurate way to globally evaluate the hosting capacity of a distribution network compared to scenario-based methods.The methodology for obtaining the region boundaries can be categorized in two categories, that is, simulation approaches and analytical approaches.
Simulation approaches obtain accurate boundary points through simulation.In [20], the boundary points were searched in 2-dimensional or 3-dimensional cross-sections of the operation region by iteratively changing one power injection with equal intervals while fixing other power injections until reaching the region boundaries.Here x-dimensional cross-section is the subspace of the operation region where x power injections (active, reactive or apparent power injections) are variable, while other power injections are fixed.To accelerate the computation, the operation point with the maximum apparent power loading for all nodes of the distribution network was suggested as the initial operation point in the searching process [20].In [17], Li et al. proposed to use an optimization method to search for the boundary points on a specified region boundary.Referring to the proposed optimization method in [17], different boundary points can be obtained along different searching directions with searching angles varying uniformly in the 2-dimensional cross-sections of the operation region.To improve the computational efficiency, the authors further proposed to use the previous obtained boundary point as the initial point for the optimization model in [18,19] to explore the next boundary point, following an orbiting route around the surface of the operation region.
The boundary points obtained in [17][18][19] rely on changing the searching direction angles, which cannot be easily determined uniformly in the high-dimensional operation region, and thus the methods normally apply to 2-dimensional or 3-dimensional cross-sections of the operation region.In our previous paper [21], we approximated the boundary points by using the outermost feasible operation points after examining random operation points through power flow equations and the network constraints.However, the computational complexity of this method grows exponentially with the increase of nodes in the electricity distribution network, preventing the method from searching boundary points of the high-dimensional operation region of a network with many nodes.In summary, the existing methods for exploring the boundary points normally compromise on only identifying 2-dimensional or 3-dimensional cross-sections of an operation region.
Another category of approaches is to use analytical expressions to approximate the operation region boundaries.There have been many studies on developing analytical boundaries for transmission networks [22], but they cannot be directly applied to distribution networks because distribution networks are with high R/X ratio (i.e., the ratio of network resistance to network reactance) and the power flow equations cannot be simplified as decoupled active/reactive power flow equations.For distribution networks, linear approximation (termed as hyperplane expressions) for voltage and thermal boundaries of an operation region were derived in [13] and [14], respectively.Combined with simulation approaches, analytical boundaries can also be obtained through curve fitting methods based on the boundary points obtained by simulation.Specific curve fitting methods applied in approximating the operation region boundaries include the least square method [20,23], the piecewise approximation method [18] and the convex hull method [24].
Compared to simulation approaches, analytical approaches can easily characterize the feasible operation region (of a distribution network) by its boundaries' analytical expressions to provide the allowable range of power injections in the network.Moreover, the analytical approaches establish correspondence between region boundaries and network constraints.As a result, the analytical boundaries can be used to replace network constraints in optimal power flow models to accelerate the calculation [14].Considering the advantages of the analytical approaches, this paper attempts to propose higher-order analytical expressions (i.e., quadratic expressions) of the boundaries of the feasible operation region (FOR) of a distribution network, in contrast to the existing hyperplane expressions, to accurately express the network capability and efficiently exploit the network headroom for integration of more DG/demand.This paper also provides a highdimensional error analysis approach for validating any forms of analytical boundaries in high-dimensional space, which is not available in existing literature.The main contributions of this paper are summarized as follows: (1) Compared to the conventional scenario-based assessment methods, the operation region-based methodology is used for accurately evaluating the capability of distribution networks to accommodate increasing generation and demand.(2) Quadratic expressions are proposed for describing both voltage and thermal boundaries of FORs, considering both active and reactive power injection at each node.Compared with the existing hyperplane expressions, the proposed quadratic expressions are more accurate especially for thermal boundaries.(3) An effective high-dimensional error analysis approach is provided for validating any forms of analytical boundaries of FORs, in contrast to the existing methods mainly applicable to lowdimensional cross-sections of the operation region.
The rest of this paper is organized as follows.Section 2 introduces the concept of feasible operation regions.Quadratic analytical expressions of FOR boundaries are deduced in Section 3. In Section 4, the four-step validation method is provided for high-dimensional validation.Section X. Jiang et al. 5 presents the validation results for the quadratic expressions of FOR boundaries compared to hyperplane expressions.Finally, conclusions are given and the remaining problems are discussed in Section 6.

Definition of feasible operation regions
A feasible operation region (FOR) is defined as the set of feasible operation states of a distribution network, where the network constraints are not violated.Considering the operation states are defined as power injections at different nodes of the distribution network, a feasible operation region can be described as follows: Ω ∈ R 2n is the feasible operation region in the complex power injection space, where n is the number of nodes (excluding the slack bus) in the distribution network.x is the complex power injection vector in a distribution network, where P=(P 1 , …, P n ) T and Q=(Q 1 , …, Q n ) T are active power vector and reactive power vector.Since each active/ reactive power injection at any node of the distribution network corresponds to one dimension of FOR, the FOR of a n-node distribution nework is a 2n-dimensional power injection space.f (V0 ,θ0) (V, θ) = x represents the power flow equations with V 0 and θ 0 as predefined voltage magnitude and phase angle of the slack bus.g(I,V, θ) = 0 represents the relationship between branch currents and node voltages expressed by Ohm's law.V and I are the node voltage vector and branch current vector, which satisfy the voltage constraints C V and thermal constraintsC T , respectively: V i is the voltage magnitude at node i, constrained by its lower limit ⃒ is the magnitude of the branch current constrained by its limit I M ij .N and B denote the set of the nodes and branches in the network respectively, the number of which are n and nb respectively.
According to the definition of the FOR, within the boundaries of the FOR are all feasible operation states, while outside of the boundaries any operation states are infeasible (see Fig. 1).In this regard, the boundaries of FOR represent all the limits to the power injections that can be hosted by a distribution network, which can provide the information of the network capability.FOR is only associated with the network topology and component parameters.As long as the operation states (i.e., net nodal power injections) are inside FOR, however the generation/consumption vary, the node voltages/branch currents will be confined within their upper/lower limits.

Boundaries of a feasible operation region
FOR of a distribution network is enclosed by several highdimensional surfaces, which are determined by voltage and thermal constraints in (2) and (3) [22].These surfaces confine the network to its normal operation without violating the network constraints.In this study, these high-dimensional surfaces are defined as the boundaries of the FOR.Considering the types of constraints in (2) and (3), the boundaries of the FOR can be further categorized into voltage boundaries and thermal boundaries.In the next section, we propose a new form of analytical expressions of FOR boundaries in a quadratic form of power injections.

Analytical expressions of thermal boundaries of a FOR
The schematic diagram of a general radial distribution network is shown in Fig. 2. The direction of the arrows at each node of the network represents the net power loading.It is noteworthy that the negative value of the net power loading at one node, in other words the net power injection, indicates that the power generation is larger than the power loading at this node.In this section (and also in Section 3.2), we use net nodal power loading in the deduction process but will replace it with net nodal power injection (i.e., the positive direction of active and reactive power at each node is defined as the injection direction) in the final expressions.
For conciseness, we define P j,eq and Q j,eq as the equivalent power loading at the receiving end node of any branch ij in the distribution network in Fig. 2. The expressions of P j,eq and Q j,eq are as follows: P j and Q j are the active and reactive power loading at the receiving end node of the branch ij, while P jk and Q jk are the active and reactive power flow at the branch jk.A j denotes the set of the adjoining downstream nodes of node j (e.g., nodes k, l instead of node i in Fig. 2) in the distribution network.
The voltage drop on the branch ij can be expressed as the voltage phaser diagram in Fig. 2. Compared to the previous study on hyperplane expressions [14] which ignores δV ij , this study considers the impact of the difference of voltage phase angles (i.e., θ ij ) and δV ij is retained.Then the relationship between the voltage drop and the power flow on the branch in Fig. 2 can be expressed as: Vi and Vj are the voltages at the sending end and the receiving end of the branch (V i and V j denote the magnitude of them).R j and X j are the resistance and reactance of the branch.
Based on Ohm's law, the current of the branch ij can be obtained by: From ( 6)-( 9), the magnitude of the branch current can be expressed as: (10) is a symmetric equation in terms of P j,eq and Q j,eq , which indicates that P j,eq and Q j,eq have the same degree of influence on the branch current.In a distribution network with high R/X ratio, the active power losses are normally larger than reactive power losses.Considering the percentage of active power losses over the power loading is small for a distribution network [25], the power losses at the downstream branches in P j.eq and Q j,eq are ignored.In addition, since the allowable variation of node voltages is small (normally within ± 3 % or ± 5 %), the voltage magnitude V j is assumed to be V 0 [14].Following these two assumptions, (10) can be simplified as: D j denotes the set of the downstream nodes of node j (including node j for clarity) in the distribution network.Setting the branch current at its upper limit, i.e., ij , the analytical expressions for the thermal boundaries of FOR can be derived as: From ( 12), the thermal boundaries are in the quadratic form of the net power loading at different nodes of the distribution network, which are different from the expressions of hyperplanes such as those in [14].When defining P k and Q k in (12) as the net power injection at node k in the deduction process, the form of the analytical expressions of thermal boundaries stays the same as (12).

Analytical expressions of voltage boundaries of a FOR
In Fig. 3, we draw a circle with radius of V i (i.e., the voltage magnitude at the sending end node of the branch in Fig. 2) to observe the difference of the voltage magnitude between the sending end node and the receiving end node.From Fig. 3, (6) can be replaced by the following scalar equation: For conciseness, XY denotes the line segment between point X and point Y. Due to the vertical relation between OC andAB, BC in can be expressed as: Extend AB to intersect the circle at A ′ .From the fact that the points A and A ′ on the circle are symmetrical about theOB, which is perpendicular toAA ′ , we have Due to the small difference of the phase angles at the two endpoits of a branch (i.e., θ ij is normally small), the numerator ( 2 , a second-order element ofθ ij .We retain the first-order small element yet ignore the second-order one and Eq. ( 15) is simplified as: Substituting ( 14)-( 16) in (13), the scalar equation of the voltage drop can be expressed as: ΔV ij and δV ij refer to (7) and ( 8) respectively.Then consider the same assumptions in Section 3.1 (i.e., the ignorance of power losses at the downstream branches and the voltage magnitude to be V 0 [14]), (17) can be simplified as: Applying (18) to all branches, we have U j denotes the set of the upstream nodes of node j (including node j for clarity) in the distribution network.By substituting (4)-( 5) in (19) and considering the effect of the power losses are small (compared with the effect of the power loading), the voltage magnitude for each node can be expressed as: Fig. 3. Relations between the voltage magnitudes for voltage drop through a branch.
From ( 20), the voltage magnitude at any node includes a first-order term and a second-order term of power loading.In line with the definition of FOR in Section 2, here we define the positive direction of active and reactive power at the node as the injection direction.Then by reversing the signs for all the load terms P l and Q l (in both the linear term and the quadratic term of (20)), the relationship between the node voltage and power injections can be replaced by (21).It should be mentioned that (20) and (21) are the same except that the positive directions for net power injections are defined differently.
Regarding the upper/lower limit of the voltage magnitude at node j, i.e., V j : = V M j /V j : = V m j , from (21) we have the analytical expressions of voltage boundaries:

High-dimensional error analysis approach
As stated in Section 3, the analytical expressions of FOR boundaries in this study are the quadratic expressions of power injections, formulating the surfaces of the FOR.The analytical expressions of FOR boundaries are validated in the high-dimensional power injection space in this section.It should be mentioned that the high-dimensional error analysis is provided to validate or compare the analytical FOR boundaries, but the network operators can directly use the analytical FOR boundaries for practical application.
A general validation method is illustrated in Fig. 4. Just as a special case for demonstration purpose, the analytical boundaries in Fig. 4 are drawn on the outside of FOR.However, caused by the errors, the analytical boundaries might be at both sides of the real boundaries in practice.Supposing a FOR in the high-dimensional space such as that in Fig. 4, the validation for the analytical boundaries of the FOR follows four steps: Step 1: Generate power-increasing directions (d 1 , d 2 , …, d k , …) from a point inside the FOR in the power injection space (P, Q) T ∈ R 2n .
Step 2: Obtain real boundary points of the FOR along the powerincreasing directions.
Step 3: Obtain the operation points on the analytical FOR boundaries along the power-increasing directions.
Step 4: Analyze the error of the analytical FOR boundaries through comparing the set of real boundary points (on each thermal/voltage boundary of the FOR) obtained in Step 2 and the set of operation points (on each corresponding analytical thermal/voltage boundary of the FOR) obtained in Step 3.
The four steps are specified in the subsection 4.1-4.4respectively.

Generation of power-increasing directions inside a FOR
As with Fig. 4, power-increasing directions (d 1 , d 2 , …, d k , …) can be represented by power injection vectors in the high-dimensional power injection space of the FOR.Considering power injections at any nodes can be bidirectional (either positive or negative), we select the origin of the power injection space (i.e., (P, Q) T = (0, 0) T ∈ R 2n ) as the starting point for all the power injection vectors.For the end points of these power injection vectors, they can be uniformly distributed in the highdimensional power injection space from the statistical view by exploiting the Marsaglia's algorithm [26].
Following Marsaglia's algorithm, we firstly generate an 2n-dimensional vector (from the origin): where each component of the power injection vector x k (i.e., p k i or q k i ) follows normal distribution N (0,1).Then we have one power-increasing direction: Provided enough power-increasing directions are generated following this step, the real boundary points on all thermal boundaries and voltage boundaries enclosing the FOR can be obtained along these directions.The operation points on different analytical FOR boundaries can also be obtained.

Obtaining real boundary points of a FOR
As shown in Fig. 4, along each power-increasing direction in Step 1, one real boundary point can be obtained by extending the initial power injection vector d k to its maximum such that the power injection vector is still inside the FOR but a small increase will make it exceed the FOR.A positive constant γ k is introduced as the coefficient to extend d k .The real boundary point along any power-increasing direction is calculated through the optimization model as follows: The objective of the optimization model is to derive the maximum γ k towards the predefined power-increasing direction, while the constraints confine the extended power injection vector γ k d k inside the FOR.With the maximum γ k , the power injection γ k d k is on the real boundary of FOR.
Moreover, the critical constraint at the optimum in the model provides the specific FOR boundary to which the obtained real boundary point belongs.By determining which inequality constraint for node voltages/branch currents in is "active" at the optimal solution, the critical constraint can be found.Here "active" refers to the fact that the optimal solution causes the inequality to be an equality.For example, if i at the optimal solution, the obtained real boundary point (i.e., the maximum power injections along the specified power-increasing direction) causes V i to reach its upper limit.This indicates that the obtained real boundary point belongs to the voltage upper boundary that is determined by It is worth noting that the optimization solvers (e.g., solvers in Matlab) for the constrained optimization problem normally determine the optimality of the solution in the iteration process by using the Lagrangian function of the optimization problem.Therefore, a convenient method to decide which constraint is active at the optimal solution is to observe which Lagrangian multiplier for the constraints is not zero.Repeatedly solving the optimization model considering different power-increasing directions, the set of real boundary points on each thermal/voltage boundary can be obtained.

Obtaining operation points on the analytical FOR boundaries
In the similar way to that in Section 4.2, the operation points on the analytical FOR boundaries can be obtained by extending the power injection vector d k such that it intersects the analytical FOR boundaries (see Fig. 4).The process can be expressed as: The positive constant ω k in ( 27) denotes the extension coefficient; g (P, Q) = 0 expresses a specified analytical FOR boundary (see (12) for each thermal boundary and ( 22)-( 23) for each voltage upper/lower boundary respectively).It is noteworthy that this process can be conducted for any forms of analytical expressions of FOR boundaries.
By solving ω k , one operation point ω k d k on this analytical FOR boundary can be obtained.Selecting different power-increasing directions which are identified towards one particular FOR boundary in Section 4.2, the set of operation points on each corresponding analytical FOR boundary can be obtained.

Boundary errors
The boundary error for an analytical FOR boundary is defined as the error between the corresponding real boundary and the analytical boundary itself.For clarity, define the set of points on one real FOR boundary and the set of operation points on its corresponding analytical FOR boundary as set A :={ a i ∈ R 2n } and set B :={ b i ∈ R 2n } respectively.The boundary error for this analytical boundary can be represented by the distance between set A and set B , which can be further expressed as the set of the distances between each real boundary point a i and set B : D(set A , set B ) denotes the distances between set A and set B .It is noteworthy that each FOR boundary can be analyzed independently when set A and set B corresponds to one FOR boundary.D(a i , set B ) denotes the distance between the ith real boundary point a i and set B , which can be represented by the distance between a i and its nearest operation point in set B as shown in (29): where : D (a i , b j ) is the distance between the ith real boundary point a i and the jth operation point b j on the analytical boundary.In this study, we use three distance functions in (29) for the measurement of the distance between two points a i and b j .a i (k) and b j (k) are the components in the kth dimension of a i and b j respectively.Based on the absolute error between a i (k) and b j (k), Chebyshev distance function is used for measuring the largest error among the errors in all the dimensions, Euclidean distance function represents the length of the line segment from a i to b j in the high-dimensional power injection space, and Manhattan distance function is calculated by the total errors in all dimensions.
To facilitate the comparison of the errors between distribution networks with different scales, D (a i , set B ) can be normalized by dividing S DN which represents the size of the FOR.In this study, S DN is estimated by the average distance (under different distance functions) between the origin and the FOR boundaries in (30): N R total denotes the total number of real boundary points of set A .Through the normalization in (30), the percentage errors between the real FOR boundaries and the analytical FOR boundaries are obtained and can be analyzed statistically.In this study, the mean, minimum, and maximum boundary errors between the real boundary points and the analytical FOR boundaries in (31)-(32) are considered for analysis.The calculations in (31)-(32) can apply to any of the three different distance functions:

Conservativeness of analytical FOR boundaries
Since the error indices in Section 4.4.1 cannot provide the information on whether the analytical FOR boundaries are inside or outside the real FOR (i.e., whether the analytical FOR boundaries are conservative or not), another index, "conservative proportion" (CP) is further proposed in (34) for measuring the conservativeness of the analytical boundaries, expressed by the proportion of operation points which are inside the real FOR over the total operation points on the analytical FOR boundaries.
N A in is the number of operation points on the analytical FOR boundaries and inside the real FOR.N A toal is the total number of the operation points on the analytical FOR boundaries.
The operation points on the analytical FOR boundaries and inside FOR can be identified by comparing γ k in Section 4.2 and ω k in Section 4.3.Towards the same power-increasing direction d k , if ω k is smaller than γ k , then the operation point ω k d k on the analytical boundary is inside the real FOR.Otherwise, the operation point is outside the real FOR.
If the whole analytical FOR enclosed by the analytical boundaries is inside the real FOR, then we regard the analytical FOR as a conservative approximation of the real FOR, in which all operation points will not violate the constraints of the real distribution network.In most cases, however, approximating the high-dimensional nonlinear real FOR boundaries results in the fact that the analytical FOR boundaries might lie on the safe side or the unsafe side or even both sides of the real FOR.Therefore, the network operators are suggested to conduct offline analysis of the operation errors to learn the worst operating risks.

Operation errors
Operation errors are defined as the physical consequences (including the overvoltage, undervoltage and overcurrent) of boundary errors.According to the definition of the FOR, all the operation points inside a FOR should be feasible and satisfy the network constraints including the thermal and voltage constraints.However, due to the boundary errors in the analytical FOR boundaries, the operation points confined by the analytical FOR boundaries might violate the network constraints.Therefore, it is important to learn how large this kind of violation will be.
The boundary points of a FOR actually correspond to the critical operation states (i.e., the worst operating conditions) of a distribution network.Therefore, we can observe the maximum operation errors by analyzing the power flow results of all the analytical boundary points.Specifically, power flow calculation is repeatedly conducted for all the analytical boundary points, so that the overcurrent, overvoltage and undervoltage problems can be observed.

Summary of the error analysis process
The overall process of high-dimensional error analysis is summarized in Fig. 5.Despite the four steps as illustrated in Section 4.1-4.4,the flow chart also provides the statistical stopping criteria [27] for sampling power-increasing directions.Once the standard deviation of the distances between the origin and the real boundary points along the sampled power-increasing directions satisfies -, the sampling process can stop: where s k denotes the distance between the origin and the real boundary point γ k d k , while { s k } represents the set of distances with k samples.z α/2 is suggested to be set at 1.96 in accordance with the confidence level of 95 %; ε is the acceptable error of the mean of distances to the predefined confidence level and it is set as 0.001MVA in this study.

Case study
In this section, a 5-node feeder and a 27-node feeder of an 11 kV high-voltage underground (HV UG) network from the United Kingdom Generic Distribution System (UKGDS) [28] were used for the case study.The feasible operation regions of these two feeders were studied, and the analytical expressions of the FOR boundaries were validated.
The computation of the case study was performed in Matlab R2019b on a PC with an Intel(R) Core(TM) i7-9700 CPU @ 3.00 GHz processor and 16 GB RAM.The simulation method in was implemented by fmicon solver in Matlab and the analytical method in used the algebraic solution in Matlab.

UKGDS 5-node feeder
The 5-node feeder studied is shown in Fig. 6.For simplification, the node 301 (slack bus) and nodes 1100-1103 in [28] were numbered as node 1 and nodes 2 to 5 in Fig. 6.The branch thermal limits and the node voltage limits of the feeder were listed in Table 1.

Results of analytical FOR boundaries
Since there are 4 power injection nodes (excluding node 1 as the slack bus) in the feeder, the dimension of its FOR in the complex power injection space is eight.In the 8-dimensional operation region, there are twelve boundaries in total, including four thermal boundaries and eight voltage boundaries, which are determined by the thermal limits and voltage limits in Table 1, enclosing the FOR of the feeder.
Applying the proposed analytical expressions of thermal boundaries in (12) and voltage boundaries in ( 22) and ( 23), these twelve boundaries can be approximated.Considering the difficulties of visualization in the high-dimensional complex power injection space, Fig. 7 illustrates the twelve analytical FOR boundaries in a 2-dimensional P 4 -P 5 crosssection.Fig. 7 also shows how these analytical boundaries proposed in Fig. 6.The schematic diagram of the 5-node feeder selected from the 11 kV UKGDS distribution network.

Table 1
Thermal and voltage limits for the 5-node feeder.this paper characterize the FOR (i.e., the allowable range of power injections in yellow) of the feeder.We further depict the real FOR of the feeder and the FOR obtained by the hyperplane analytical expressions in 2-dimensional cross-sections (including P 4 -P 5 , P 4 -Q 5 , P 5 -Q 4 and Q 4 -Q 5 cross-sections) for comparison.The results are shown in Fig. 8.The results in these four crosssections show that the quadratic analytical boundaries are more accurate than hyperplane ones, especially considering the impact of reactive power injections on the branch currents.However, we cannot make these conclusions in the whole power injection space since the errors between the analytical FOR boundaries and the real FOR boundaries differ in different cross-sections.This requires the high-dimensional error analysis of the analytical FOR boundaries in the whole power injection space.

Demonstration of high-dimensional error analysis
In this subsection, the high-dimensional error analysis method developed in this study was demonstrated in the 2-dimensional P 4 -P 5 cross-section as shown in Fig. 9.In Fig. 9, 50 power-increasing directions were generated by Marsaglia's algorithm [26] for Step 1 of the method (see Section 4), along which the real boundary points and the operation points on the anlaytical FOR boundaries can be obtained following Step 2 and Step 3 respectively.The FOR and the innermost analytical boundaries in P 4 -P 5 cross-section are also displayed in Fig. 9 for reference.It is noteworthy that the developed high-dimensional error anlaysis can apply to any forms of anlaytical FOR boundaries besides the quadratic analytical FOR boundaries presented in Fig. 9. Following Step 4 in Section 4, the errors of the analytical FOR boundaries can be analyzed through comparing the set of the real boundary points and the set of the operation points on the anlaytical FOR boundaries.The error analysis results in the whole power injection space were presented in the Fig. 10. of the boundary errors between quadratic expressions and hyperplane expressions of (a) thermal boundaries and (b) voltage boundaries in the 5-node feeder measured by multiple distance functions.following section.

Error analysis results
100,000 power-increasing directions were generated in the whole power injection space for the high-dimensional error analysis.Through obtaining and analyzing the real boundary points, eight boundaries (three thermal boundaries and five voltage boundaries that are determined by thermal limits of branches 2-3, 3-4 and 4-5, voltage upper limits of nodes 2-5 and voltage lower limit of node 5) were found to constitute the real FOR for the 5-node feeder.The composition of FOR boundaries indicates that the allowable rage of the power injections of the case network is confined by these thermal/voltage limits.It is noteworthy that the lower voltage boundaries that are determined by the lower voltage limits for nodes 2-4 are beyond the FOR and are thus not FOR boundaries.The reason for this is that the increase of power loading will first increase the branch currents to their limits before trigering low voltages at nodes 2-4.
The error analysis for the quadractic expressions of the eight FOR boundaries was presented as follows.The quadratic expressions proposed in this study were also compared with the hyperplane expressions used.From a statistical point of view, the analytical expressions of FOR boundaries in Fig. 10 shows similar error results with different distance functions (i.e., Chebyshev, Euclidean and Manhattan).Whichever distance function is used, the mean boundary errors for quadratic thermal boundaries and voltage boundaries are less than 3.0 % and 1.6 %, respectively.The maximum boundary errors of quadratic thermal boundaries can be confined within 12.8 %, while the maximum boundary errors for quadratic expression of voltage boundaries are 14.9 %.
In comparison, the results of hyperplane expressions show larger boundary errors in both thermal and voltage boundaries.In particular, hyperplane expressions are not able to describe thermal boundaries accurately, whereas the boundary errors of the voltage boundaries are close to those of quadratic expressions.

2) Conservativeness of analytical FOR boundaries
The analytical boundary points which are inside FOR can be identified by determining whether the obtained ω k in is smaller than γ k in along the same power-increasing directions.By counting the number of  the analytical boundary points which are inside the FOR, the conservativeness of the analytical FOR boundaries can be assessed with detailed results listed in Table A.1.For the 5-node feeder, 61.5 % of the quadratic analtyical FOR boundaries is conservative, while 56.3 % of the hyperplane anlaytical boudnaries is conservative.Since the overall convervative proportion of the analtyical FOR boundaries is not 100 % (for both quadratic and hyperplane expressions), there must be some critical operation points, which are within analytical FOR boundaries but will violate the network constraints.Therefore, analyzing the opertaion errors is important to decide whether the errors of the anlaytical FOR boundaries are acceptable in practice.

3) Operation errors
The network operating states of the analytical boundary points were calculated to observe the maximum operation errors.The results for quadratic expressions and hyperplane expressions are compared in Fig. 11.The maximum operation errors of quadratic expressions of FOR boundaries for the overcurrent are up to 105.3 % and for the undervoltage are down to 0.968p.u..This level of thermal and voltage violation is acceptable in practice.
The results in Fig. 11 show that there are some aggressive operation points on hyperplane analytical boundaries, which cause the network constraint violation up to an unacceptable level.Most aggressive operation points come from the significant boundary errors in the hyperplane thermal boundaries.

UKGDS 27-node feeder
The longest feeder from the UKGDS HV UG network (i.e., the 27node feeder) [28] is shown in Fig. 12, where nodes 1-27 denote the node 301 (slack bus), node 1100 and nodes 1151-1175 respectively in [28].The branch thermal limits and the node voltage limits of the feeder refer to [28].Since the 27-node feeder is much longer than the 5-node feeder in Section 5.1, the total power losses and the voltage difference between the power injection nodes (especially the end node) and the slack bus can be larger.Hence, this feeder is used to further verify the accuracy of the quadratic expressions of FOR boundaries in this study.
In the 27-node feeder there are 26 power injection nodes and 26 branches.As a result, 26 thermal boundaries and 52 voltage boundaries in the 52-dimensional complex power injection space can be identified to characterize the FOR of the feeder.The quadratic analytical expressions for these boundaries can be written as in (12), (22) and (23).
The validation for these analytical boundaries can be conducted by the high-dimensional error analysis method in Section 4. 100,000 power-increasing directions were used for searching the real boundary points of the FOR and the operation points on the analytical FOR boundaries.By comparing the set of real boundary points and the set of operation points on analytical boundaries, the boundary errors were obtained in Fig. 13.Comparing to the results of the 5-node feeder, the mean boundary errors increase to 4.6 % and 10.3 % for quadratic expressions of thermal and voltage boundaries respectively, while the maximum boundary errors to 22.7 % and 123.3 %.However, the quadratic expressions proposed in this study show greater advantages over the hyperplane expressions in characterizing the FOR of the longer 27-node feeder.
The conservative proportion of the quadratic analytical boundaries is 77.4 %, while the conservative proportion of the hyperplane analytical boundaries is 51.4 %.Hence the operation errors were further analyzed and the results were shown in Fig. 14.The maximum operation errors of the quadratic analytical boundaries for the overcurrent are up to 116.0 % and for the undervoltage are down to 0.960p.u., which is also acceptable in practice.The reason why the maximum boundary error reaches 123.3 % but the operation errors are small is that, the maximum boundary error happens at the lower voltage boundary of node 20, where the conservativeness of the corresponding analytical boundary is 100 %.For hyperplane analytical boundaries, most aggressive operation points come from the significant boundary errors in the hyperplane thermal boundaries.

Computation time
The computation time for obtaining the FOR boundary points by the simulation method in (26) and the proposed analytical FOR boundaries in ( 27) is summarized in Table 2.As observed in Table 2, the analytical method spent 0.08 ms to calculate one FOR boundary point for the UKGDS 5-node feeder while the simulation method took 41.97 ms, achieving five-hundredfold speed up.Regarding the UKGDS 27-node feeder, the proposed analytical method is even thousandfold faster than the simulation method.The reason is that compared to the simulation method, the proposed analytical method can directly obtain the boundary point through algebraic computation, avoiding repeatedly iterations in the simulation method.In addition, the proposed method can directly use analytical expressions for complete characterization of the FOR of the network, while the simulation method requires much more work on point-based simulation.

Error analysis of power losses assumption
The overall error anlaysis results in Section 5.1.3have shown that the errors caused from the assumptions of the analytical FOR boundaries are acceptable in practice.This section provides further analysis of the impact of power losses assumption on the results of branch currents and node voltages.
Power losses are assumed to be ignored in the branch current equation (11) and node voltage equation (20).Under this assumption, the equivalent power loading at node j (i.e., P j,eq /Q j,eq ) can be estimated by the summation of the power loading at the downstream nodes of node j (i.e., With the actual peak power loading conditions in the test systems, the errors of the equivalent power loading caused by  the ignorance of power losses were summarized in Table 3.Since the maximum errors should appear at the first PQ node of the system, Table 3 only presents the errors for P j,eq and Q j,eq at the first PQ node in different test systems.
From Table 3, the maximum errors for both equivalent active and reactive power are increased as the feeder length of the test systems increases.In addition, the percentage error for the equivalent reactive power loading seems very large especially in long-feeder system.However, the magnitude of the error for the equivalent reactive power loading is actually smaller than that for active power loading since the reactance is smaller than the resistance in the distribution networks.
We further calculated the errors of all branch currents and the errors of all node voltages, which result from the ignorance of power losses in calculating the equivalent power loading in Table 3.The results of the maximum error of all the branch currents and node voltages were summarized in Table 4. From the obtained results, the errors of the branch currents and the node voltages are acceptable (even if the error of the equivalent reactive power loading in Table 3 can be 11.28 % in the UKGDS 27-node test system).

Conclusions and discussions
This paper lays the foundation of using operation region-based methodology to express the capability of distribution networks.The conclusions are presented as follows: (1) The boundaries of a feasible operation region (FOR) can be used to characterize the allowable range of power injections at all nodes of a distribution network.
(2) The analytical expressions of both thermal and voltage boundaries of a FOR have been formulated in the quadratic form of the power injections at the nodes of a distribution network.A highdimensional error analysis approach has been provided to validate any forms of analytical expressions of FOR boundaries in a high-dimensional operation region.(3) From the results of the UKGDS 5-node and 27-node feeders in the case study, the proposed quadratic expressions of FOR boundaries can well approximate the real boundaries of the FOR.The mean boundary errors are less than 6.1 % and 10.3 % in the 5node feeder and the 27-node feeder respectively, while the maximum boundary errors can be confined within 14.9 % and 123.3 % respectively.Despite the increase of the boundary errors of analytical FOR boundaries in the 27-node feeder compared to those in the 5-node feeder, the worst operation states have the overcurrent up to 116.0 % and the undervoltage down to 0.960p.u., which are acceptable in practice.This is because the maximum boundary error happens at a boundary with 100 % conservativeness.(4) The proposed quadratic FOR boundaries proposed in this paper are more accurate than the existing hyperplane expressions.In particular, the hyperplane expressions of thermal boundaries have large errors and thus should be examined before implemented in practice.
The main challenges with some suggestions are also summarized below for future research: (1) The key to deriving the analytical expressions of FOR boundaries is to obtain the explicit relationship between branch currents (or node voltages) and power injections, which is very difficult without any assumptions.The errors from these assumptions can be quantified by the high-dimensional error-analysis approach proposed in this study.In the future, combined analytical and data-driven methods may further reduce the errors, as an extension of the analytical method and high-dimensional error-analysis approach presented in this study.
(2) In this paper, two UKGDS feeders were used for case studies.
Regarding a large-scale distribution network, users of FOR methodology can consider two strategies to simplify the computation.Firstly, a large distribution network can normally be divided into different feeders or subnetworks, of which the FORs can be analyzed independently and in parallel.In this regard, the FOR of each subnetwork is only determined by the topology and parameters of the subnetwork itself.Moreover, the generation/ consumption at the downstream nodes (or in lower-voltage networks) which are not concerned in detail can be deemed as one aggregated power injection for simplification.Based on these considerations, the test cases in this paper are able to validate the methods sufficiently.
(3) In a distribution network with power electronics, the network operators using the analytical expressions of FOR boundaries in this study should add the power injections from the power electronics into the net nodal power injections before assessing the operation condition of the network.How to characterize the hosting capacity of the network reinforced by power electronic devices is worthy of further studies in the future.

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.

Fig. 1 .
Fig. 1.Relationship between hosting capacity and feasible operation region for a distribution network.

Fig. 4 .
Fig. 4. Schematic diagram of high-dimensional error analysis for validating analytical FOR boundaries.

Fig. 5 .
Fig. 5. Flow chart of high-dimensional error analysis for analytical FOR boundaries.

Fig. 8 .Fig. 9 .
Fig. 8.Comparison of the approximate FOR enclosed by analytical boundaries with the real FOR in 2-dimensional cross-sections.

Fig. 11 .
Fig. 11.Network operating states of analytical FOR boundary points for the 5-node feeder: (a) branch currents and (b) node voltages.Note that power flow calculations through the Newton-Raphson algorithm cannot converge for 3% extreme operation points with very large power injections for the hyperplane thermal boundary points.

Fig. 14 .
Fig. 14.Network operating states of analytical FOR boundary points for the 27-node feeder: (a) branch currents and (b) node voltages.Note that power flow calculations through the Newton-Raphson algorithm cannot converge for 11% extreme operation points with very large power injections on the hyperplane thermal boundaries.

Table 2
Comparison of the time consumption between the analytical method and the simulation method.

Table 3
Errors of the equivalent power loading at the first PQ node in different test systems.

Table 4
Impact of ignorance of power losses on the results of branch currents and node voltages.