Hosting Capacity Assessment in Distribution Networks Considering Wind–Photovoltaic–Load Temporal Characteristics

Under the background of clean and low-carbon energy transformation, renewable distributed generation is connected to the distribution system on a large scale. This study proposes a probabilistic assessment method of hosting capacity considering wind–photovoltaic–load temporal characteristics in distribution networks. First, based on time series of wind, photovoltaic, and load demands, a discretization–aggregation technique is introduced to generate and filter extreme combinations. The method can effectively reduce the scenarios that need to be evaluated. Then a holomorphic embedding method considering generation and load scaling directions is proposed. The holomorphic function of voltage about an embedding variable is established, and it is analytically expanded in the form of series. The hosting capacity restrained by the voltage violation problem is calculated quickly and accurately. Finally, the proposed stochastic framework is implemented to evaluate hosting capacity involving renewable energy types, penetration levels, and locations. The hosting capacity of single energy and hybrid wind–solar renewable energy systems is evaluated from the perspective of probability analysis. The results verify the outstanding performance of the hybrid wind–solar energy system in improving the hosting capacity.


INTRODUCTION
In response to climatic deterioration and energy shortage, all countries are accelerating the process of new energy. Distributed renewable energy sources have become the mainstay to promote the development of new energy with the advantages of being clean, green, flexible, and efficient (IEA, 2019). Wind energy and solar energy are the most promising renewable energy sources. However, their access to the distribution network also brings uncertainty and intermittence. The booming development of distributed generation (DG) may lead to the violation of system operation constraints such as overvoltage (Ismael et al., 2019;Zhu et al., 2020), overloading of transformers and feeders , conductor thermal capacity (Zhang and Luo, 2018), and protection failure (Singh, 2017;Zobaa et al., 2020). In order to overcome the challenges of renewable energy source integration, it is of significant importance to evaluate the number of DGs that can be integrated into a given distribution network without violating the operating standards.
The concept of hosting capacity (HC) was first proposed by André Even in the context of distributed generation and improved by Bollen and Hassan (2011). The hosting capacity is defined as the maximum capacity of DGs that can be integrated into the distribution system, above which the performance of the system becomes unacceptable. Recently, many scholars have studied hosting capacity assessment in distribution networks. There are four main methods: the deterministic method, the stochastic method (Yang et al., 2019), the optimization-based method (Shen et al., 2017;Injeti and Thunuguntla, 2020), and the time series method (Abideen et al., 2020;Mulenga et al., 2020).
In earlier studies, the analysis methods were often used to calculate the HC at a specific DG access location by deriving the performance index of the system. In the study of Fan et al. (2017) and Li et al. (2021), the formula of voltage difference values at continuous buses of three-phase feeders is derived, and the maximum number of DGs at a specific bus is calculated. Ampofo et al. (2017) studied the impact of voltage rise and thermal loading on HC considering DG access to the end of feeders or the load center. HC is calculated in different scenarios by iteratively increasing the number/capacity of wind generation units and continuously calculating the power flow (PF) until one of the performance standards is violated (Papaioannou and Purvins, 2014;Gonzaga et al., 2019). It is evident that the deterministic method cannot consider the uncertainty of modern power systems, and its application range is limited.
When DGs with high uncertainty characteristics are connected to the distribution network, there are many unknown variables in the calculation of HC. Thus, the randomness of these variables needs to be considered; Monte Carlo simulation (MCS) is often used to generate different scenarios. Zio et al. (2015) proposed a probabilistic power flow method and simulated the variability of customer demand based on MCS but did not consider the variability of DG. The randomness of both DG access locations and load demand are considered (Kolenc et al., 2015;. In Al-Saffar et al. (2019), the probabilistic power flow is implemented under the scenarios with different photovoltaic (PV) penetration levels, and the HC of three real regions is determined, respectively. Mulenga et al. (2021) classify two types of uncertainties, namely, aleatory uncertainties and epistemic uncertainties. The HC is estimated by applying the transfer impedance matrix and the superposition principle to determine the voltage rise due to PV. In addition, using spatial and temporal uncertainties associated with PV, a new spatiotemporal probabilistic voltage sensitivity is proposed. It can calculate the probability distribution of voltage change at a specific bus, due to random change of PV power in the random position of the network (Munikoti et al., 2022).
The optimization-based method is also a common approach to determine the HC. The objective is to maximize the DG injection while constraints are met. In the study of Zou et al. (2016), Alturki et al. (2018), and Shen et al. (2020), based on deterministic optimization algorithms, the best access location is regarded as the main solution. But in fact, the inherent uncertainty of DG needs to be considered. Therefore, the trend is combining the stochastic method and the optimization-based method . A stochastic multi-objective optimization model was proposed in the study by Rabiee et al. (2017), which aims to maximize the HC for wind power and minimize the energy procurement costs, and then it is solved with the NSGA-II algorithm. Otherwise, the chance-constrained method was adopted, and the probabilistic power flow method was used to deal with the randomness problem (Sun et al., 2018;Wu et al., 2019). However, the optimization-based model is generally highly complex and non-linear, and for actual networks, the existing methods may not produce global optimal solutions.
Besides, in some studies (Khoshkbar-Sadigh et al., 2015;Fan et al., 2016;, the historical data of both demand and renewable production are used as the input, and it can provide a more realistic distribution network. Chen et al. (2018) considered temporal characteristics of wind power, PV, and load; the joint probability distribution method and the scene reduction technique were used to solve the DG capacity. Mulenga et al. (2021) studied the influence of time of day on the HC calculation results. However, the time series method considering time-varying renewables and demands are highly dependent on data, and a large amount of data enlarge the computing scale, which tends to be laborious or intractable. Some scholars study the security-constrained unit commitment (SCUC) (Yang et al., 2018;Liu et al., 2020;Yang et al., 2021). Yang et al., 2021 is a pioneer study for SCUC problems that proposes an expanded sequence-to-sequence (E-Seq2Seq)-based data-driven SCUC expert system. It can accommodate the mapping samples of SCUC and consider the various input factors that affect SCUC decision-making, possessing strong generality, high solution accuracy, and efficiency over traditional methods. To mitigate the excessive computational burden, Ochoa et al. (2010) proposed a processing technique for long-term time data, namely, the discretization-aggregation method. It can generate and screen out the reasonable combinations of renewables and demand to simplify data.
When excessive DG penetrates in the distribution network, the radial distribution system with the single power becomes a complex system with multiple power supplies. Then there are reverse power flows, which may lead to voltage rise (Mulenga et al., 2020;Shen et al., 2020;Wang et al., 2021). The studies have shown that the voltage rise is the main restriction considered in the research of HC (Torquato et al., 2018;Dong et al., 2019).
In this study, a stochastic framework of hosting capacity assessment is proposed considering the uncertainty of DG penetration levels, locations, and types, and extreme combinations are introduced to process wind-photovoltaic-load time series data. This effectively reduces the number of scenarios to be evaluated. Moreover, traditional methods of hosting capacity assessment are scenario-based and complex as they rely on the iterative PF algorithm. To avoid a large number of PF calculations, a novel holomorphic embedding method (HEM) based on the recursive algorithm is used to obtain the equivalent analytical formula of voltage (Trias. 2012). The HC corresponding to voltage violation can be directly solved without checking a large number of scenarios, which further significantly reduces the computational burden. In the simulation analysis, hosting capacity assessments of both single resource and hybrid cases are performed, and the results provide planners with a better understanding of the energy integration. The rest of this article is organized as follows: Processing of Renewables and Demand Data discusses the processing technology of time series data of renewables and demand. Holomorphic Embedding Method introduces the holomorphic embedding method. Then the stochastic framework of the hosting capacity assessment is illustrated in Framework for Hosting Capacity Assessment. Numerical Results presents the results and discussions to evaluate the hosting capacity of single and hybrid cases on the IEEE 33-bus system. Finally, Conclusion summarizes the main conclusions.

Discretization-Aggregation Method
Due to the uncertainty and volatility of renewable generation and load, the temporal characteristics of both generation and load demand need to be considered in the hosting capacity assessment. However, long time series will bring a significant number of calculations. Therefore, the discretization-aggregation method is introduced to reduce the computational burden. The technology was first proposed by Ochoa et al. (2010), which only considers wind and load. Furthermore, if we consider the correlation between wind power, PV power, load demand, and time, each data point needs a multidimensional representation. The method has the potential to deal with the problems of multidimensionality.
The method mainly includes two steps: 1) in the discretization process, the historical data of renewables and demand are allocated into a series of bins covering the range between zero and the peak value; and 2) in the aggregation process, the bins of renewables and load demand are grouped into multiple combinations. To illustrate the approach, Figure 1 presents the discretization-aggregation process with only two dimensions in the following example. Figure 1A shows a 5day historical data sample of wind power and load with an interval of 15 min, and their values are normalized against respective peak values. Figure 1B shows the discrete time series. When the width of bins is set to 0.1 p.u, six load demand ranges (e.g., [0.4, 0.5], (0.5, 0.6], . . .) and nine generation ranges (e.g., [0, 0.1], (0.1, 0.2], . . .) are used. Then the time-varying data are allocated to a series of bins. Figure 1C presents the distribution with combinations of wind power and load. The combines of "similar" characteristics are aggregated into the same bin. For instance, the yellow block indicates the data where demand is 0.6 and wind is 0.2.

Extreme Combinations
The discretization-aggregation method allocates time series data into a finite number of combinations, which reduces the number of combinations to be evaluated, and retains the relevance between renewables and demand. Importantly, it does retain extreme characteristics. The "coincidence" between maximum generation and minimum demand is normally regarded as the extreme combination for voltage violation and the main constraint of the hosting capacity assessment. The discretization-aggregation process goes through each possible combination and sums the occurrence periods, which captures the full range of generation and load. Figure 2 presents all combinations of data above and their occurrence periods. The combinations labeled with red represent the extreme conditions of maximum power generation and minimum load demand. If the voltage constraint is not violated in these combinations, it is unlikely to be violated in other combinations.
It is obvious that different combinations will be obtained by selecting different widths. The smaller the width, the more detailed the bins to be evaluated. Figure 3 shows the combinations of wind and load when the width is set to 0.05 p.u. The total number of combinations increases significantly, and the number of extreme combinations has only increased by one compared with the results in Figure 2. Therefore, the selection of bin width may affect the scale and accuracy of the  hosting capacity assessment. Moreover, when an additional PV power is added, the discretization-aggregation process remains unchanged, but the dimension is increased.

HOLOMORPHIC EMBEDDING METHOD
The holomorphic embedding method was applied to the PF problem by Dr. Antonio Trias for the first time (Trias, 2012), to avoid the non-convergence problem of the traditional iterative PF methods. In this study, the conventional HEM is improved considering the direction of generation and load change. The equivalent analytical formulations of voltage can be obtained by only one PF calculation. It can establish the dependence between the embedding parameter and the actual operation level.

Direction-of-Change Scaling Holomorphic Embedding Model
Consider an N-bus system, the power balance equation (PBE) can be expressed as follows: where Y ik is the (i, k) element of the bus admittance matrix, _ S i and _ V i are the complex power injection and voltage at bus i, respectively.
The non-holomorphic PBE is converted into holomorphic functions by embedding a complex parameter _ s. Considering different types of buses, the improved holomorphic embedding formulas are given, where Eq. 2 represents the voltage magnitude constraint for slack bus, Eq. 3 represents the PBE for the PQ buses, Eq. 4 represents the PBE for the PV buses, and Eq. 5 represents the voltage magnitude constraint for the PV buses. The formulas allow the load at all buses and the real power generation at the PV buses to be scaled directionally.
where P i0 , ΔSi _ , ΔP i , ΔQ i are given as follows: where V sp i is the reference voltage amplitude; P Gi0 and Q Gi0 represent the active injection power and active load of bus i under the initial loading level, respectively; Q Li0 is the reactive load of bus i under the initial loading level; and k Gi and k Li are generation growth coefficient and load growth coefficient, respectively, which can represent the change direction of generation and load. N slack , N PQ , and N PV represent the sets of slack bus, PQ buses, and PV buses, respectively.
Since _ V i ( _ s) and Q Gi ( _ s) are holomorphic functions of the parameter _ s, they can be expanded in the following Maclaurin series form: where the voltage sequence coefficients _ V i [n] are complex numbers, and the power sequence coefficient Q gi [n] are real numbers.
The Maclaurin series for _ V p i ( _ s p ) is given as follows: Additionally, let W(s) represent the inverse of the voltage function V(s), defined as follows: The relationship between _ W i ( _ s) and _ V i ( _ s) is shown as follows: The relationship between _ W i [n] and _ V i [n] is obtained as given in Eq. 15 by equating the coefficients of the same order of s on both sides of Eq. 14.
where Eq. 15b can also be formulated as follows: By substituting Eqs 10-13 into Eqs 2-5, we obtain the following: Thus, we establish the recursion relationship from the aforementioned holomorphic embedding formulas to obtain the equations. Then _ W i [n] and _ V i [n] are divided into real parts and imaginary parts for calculation, respectively, and the voltage sequence is solved. Finally, the equivalent analytical expression of voltage can be obtained.

Reference State Calculation
To solve the system of equations above, the reference state at s 0 is given as follows: Notice that the meaning of the reference solution of the improved HEM is different from that of the conventional HEM (Rao et al., 2016). The reference solution of the conventional HEM represents the power system with no load and no generator, while the reference solution of the improved HEM represents the voltage and reactive power injections at the buses for the power system under the initial loading level. The solution process is as follows: and Q Gi [0] are expressed as the holomorphic function, so a complex parameter _ σ is embedded in Eqs 22-25, and we obtain: where Y sh i corresponds to the shunt part of the admittance matrix, Y tr ik corresponds to the "non-shunt-branch" part of the admittance matrix, and Y tr To establish the recursive relationship of variables, a new variable δ ni is defined as follows: For the slack bus, the power series coefficient expression is written as follows: The voltage power series coefficients of PQ bus and PV bus are solved, we can obtain the following equations: where can be solved. Obviously, the process of solving the reference solution of the improved HEM is the same as that of the conventional HEM when the embedded variable is 1, which is actually the voltage solution of the traditional power flow equation (Eq. (1)). Therefore, the Newton-Raphson method can also be used to solve the PF at the initial loading level to obtain the reference solution.

Calculation Process
The process of using improved HEM to solve a PF problem and voltage violation is as follows: 1) The PBEs are embedded with the parameter _ s, the voltage and the active power injections become the holomorphic function of _ s, and the holomorphic embedding models are established. 2) Taking the direction of generation and load change into account, the growth coefficients kGi and kLi are defined, and then ΔṠ i , ΔP i , and ΔQ i are calculated.

FRAMEWORK FOR HOSTING CAPACITY ASSESSMENT
Considering the temporal characteristics of renewable energy generation and load demand, this study proposes a framework for evaluating the hosing capacity in distribution networks with DGs. The framework consists of three modules, as shown in Figure 4.

1) Module 1: Deployment schemes of DGs.
This module generates multiple potential DG deployment schemes. The variables include DG location penetration, the locations of DGs, and the types of DGs. The steps are as follows: Step 1: Identify the location penetration Loc i . DG location penetration is defined as the ratio of the number of selected DG locations to the number of all potential locations. The location penetration is increased by a fixed step (e.g., 10%) from 0 to 100%. Let Loc i ∈ 10%, 20%, . . . , i × 10%, 100%} { , (i 1, 2, . . ., 10).  Step 2: Generate DG locations. For each location penetration level Loc i , MCS is performed to generate k DG deployment schemes. Then the deployment scheme is represented as D ij (j 1, . . ., k).
Step 3: Set the types and shares of DGs. For example, 50% wind and 50% PV.
Step 4: Determine a base DG capacity (e.g., 1 MW). For each deployment scheme, the initial rated power of DGs is allocated based on the corresponding peak load.
2) Module 2: Calculation method of the hosting capacity. This module studies the calculation method of the hosting capacity based on the improved HEM, considering temporal characteristics of DGs and load demand. The steps are as follows: Step 1: Process the historical data of wind, PV, and load. By the discretization-aggregation method, select the proper bin width and obtain m extreme combinations Sn (n 1, . . ., m).
Step 2: For a specific deployment scheme Dij, perform HEM calculation on different extreme combinations, respectively. Then obtain the equivalent analytical functionV i ( _ s) of all buses on each combination.
Step 3: If none of the bus voltages exceed 1.05 p.u, let | V i ( _ s) | 1.05, and calculate the value of _ s.
Step 4: Compare the minimum value of _ s on each extreme combination S n , that is, _ s min corresponds to the maximum hosting capacity under deployment scheme D ij .

3) Module 3: Analysis of hosting capacity results. Repeat steps in
Module 2 to obtain the hosting capacity results for each deployment scheme D ij , and obtain the hosting capacity results HC {HC 1 , . . . , HC k }. Then perform statistical analysis for the obtained hosting capacity results. The steps are as follows: 1) Histogram of HC is obtained based on the results HC {HC 1 , . . . , HC k }, and the probability density function (PDF) based on Kernel density estimation is helpful to understand the probabilistic HC at a specific penetration level.  2) Cumulative distribution function (CDF) curve of HC is helpful for planners to estimate the probability that the HC does not exceed a specific value. 3) Histogram of total delivered generation. The total delivered generation of DG is a valuable quantitative indicator of the use of renewable energy (Bawazir and Cetin, 2020). Therefore, some statistical data of total delivered generation in a year can provide the energy utilization of different types of renewable energy.

NUMERICAL RESULTS
The simulations are carried out on the IEEE 33-bus distribution system. The detailed parameters of the test system are available in Baran and Wu. (1989). Bus 1 is set as the slack bus and the voltage is set to 1.0 p.u. Other nodes are PQ buses. The reference voltage is 12.66 kV and the reference capacity is 10 MV A. The upper voltage limit of each bus is set to 1.05 p.u. Nodes 2-33 are candidate buses accessible to DG. In this section, initially, the historical data of wind, PV, and load demand are processed by the discretization-aggregation method, and the extreme combinations are filtered. Then the proposed HEM is used to solve the voltage violation problem on extreme combinations. For a specific deployment scheme, the hosting capacity results with different bin widths are discussed. Finally, detailed hosting capacity assessments of both single resource and hybrid cases are performed.

Renewables and Demand Data
The simulations use the historical data of wind speed, solar irradiation, and load demand from a typical distribution system. The data of one year have a total of 35,040 data points with 15-min temporal resolution. The levels of load demand, wind, and PV output are normalized against peak values, as shown in Figure 5. It can be seen that load demand and PV have both obvious seasonal characteristics. The load in summer is relatively lower than that in winter, and the PV generation in summer is significantly higher than that in winter, but for wind power, the feature is not so obvious. It should be noted that different buses are close geographically in the test system, and the potential of renewable energy power generation is similar to a certain extent. So it is assumed that DGs follow the same time series curves.
The combinations of wind, PV, and load demand are obtained by the discretization-aggregation technology in Introduction. Figure 6 shows the combinations and occurrence periods of renewable generation and load when the bin width is 0.05 p.u.
For the wind-load case, there are a total of 400 combinations, but only 175 contain non-zero occurrence periods. Similarly, for  the PV-load case, the occurrence periods of only 133 combinations are non-zero. For renewable energy, wind power and PV are negatively correlated. PV generation mainly depends on solar radiation, and solar energy at night can be ignored. In contrast, the wind power during the day is usually less than that at night. Therefore, there is a certain complementarity between wind power and PV (Miglietta et al., 2017;Guozden et al., 2020). The extreme cases of maximum generation and minimum demand are critical on the hosting capacity assessment restrained by voltage rise. The results of extreme combinations are given in Table 1. For the wind-load case and the PV-load case, only three and four extreme combinations need to be considered, respectively. For the wind-PV-load case, there are 2,567 combinations. It is difficult to show them in visual graphics, but 394 extreme combinations can be screened by the discretization-aggregation technique. Compared with the use of original historical data, the introduction of extreme combinations can significantly shrink the calculation scale of multidimensional problems.   Hosting Capacity Assessment

Hosting Capacity Calculation of a DG-Specific Deployment Scheme
Taking a specific DG deployment scheme as an example, the influence of bin width is discussed. And the rapidity and effectiveness of the hosting capacity calculation method based on HEM are verified. Wind power is connected at Bus 2, 7, 24, and 33. Four widths are considered, namely, 0.1 p.u, 0.05 p.u, 0.01 p.u, and 0.001 p.u. The hosting capacity is calculated according to Module 2 in Framework for Hosting Capacity Assessment. For comparison, the Newton-Raphson power flow method is used to calculate the hosting capacity on each extreme combination, by increasing the total DG capacity until the upper voltage is violated. The hosting capacity result is 10.872 MW from all historical data, which is regarded as the accurate value. Using the calculation method based on HEM in this study, the results of extreme combinations and the hosting capacity with different bin widths are shown in Table 2.
The results in Table 2 show that the method proposed in this study can greatly shorten the calculation time and improve the calculation efficiency. When the bin width is set to 0.1 p.u, the result is 12.97% smaller than the accurate value. Therefore, larger width may underestimate the hosting capacity and lead to conservative results. On the contrary, when the width is less than or equal to 0.01 p.u, the relative error is less than 1%. We can conclude that when the appropriate bin width is selected, the calculation scale can be simplified by using the proposed extreme combinations. More importantly, the hosting capacity calculation method based on HEM can further shorten the calculation time. Figure 7 shows that the voltage violation occurs when the width is 0.01 p.u. Under the deployment scheme, when the extreme scenarios are {(0.81, 0.82], (0.59, 0.60]}, Bus 24 first exceeds the upper voltage limit.

Probability Assessment of the Hosting Capacity
In the test system, the number of DG candidate buses is 32. When the location penetration level is 50%, the total possible number of DG deployment schemes is more than 6 × 10 8 . A large number of potential DG deployment schemes bring a significant computational burden. Therefore, MCS is used to simulate relatively few scenarios to obtain approximate results. In this context, we adopt the variance coefficient β (Prusty and Jena, 2017;Shen et al., 2020). β ≤ 0.5% is set as the stopping criterion of MCS, and the results of MCS are considered to be accurate.
We performed the stochastic framework to access the hosting capacity of a single DG case in Framework for Hosting Capacity Assessment. The bin width is set as 0.01 p.u, and the location penetration level is increased from 10 to 100% by a fixed step of 10%. Table 3 shows the statistical results of the wind system and PV system at location penetration levels of 30, 50, 70, and 90%.
The mean values of the hosting capacity increase with higher location penetration levels, while the standard deviations show the opposite trend. The reasons are as follows: first, with the increase of the location penetration level, DG locations increase and the capacity allocated to each location decreases, so the total hosting capacity of the system increases. Second, the higher penetration level reduces the uncertainty of DG locations, so the results of different deployment schemes are more concentrated. Figure 8 shows the probability density distribution of the hosting capacity at 50% location penetration level. For the hosting capacity of the wind system, the minimum and maximum values are 3.995 and 14.224 MW, respectively. For the hosting capacity of the PV system, the results of the hosting capacity are 4.236 and 14.472 MW, respectively. The difference between extreme values is due to the location distribution of DGs. Figure 9 shows the cumulative probability distribution of the hosting capacity, which can provide planners with the probability that the hosting capacity in distribution networks is lower than a specific value. For example, the probability of the hosting capacity, which is no more than 10 MW in the wind system, is 0.9359, while in the PV case, the result is 0.9097.
The total delivered generation can realistically reflect the energy utilization. Therefore, when planning the distributed system, not only the maximum hosting capacity of the system but also the total delivered generation is needed to be considered. The comparison of the mean values of total delivered generation of the wind system and the PV system is presented in Table 4. Although the hosting capacity of the PV system is slightly higher than that of the wind system, the total delivered generation of PV  is less than 50% of that of the wind system. The main reason for the difference is the higher correlation and capacity coefficient between wind resource and load demand. There is complementarity between wind energy and photovoltaic energy, and their joint action may affect the hosting capacity of the system. However, the existing literature rarely discusses the hosting capacity of hybrid wind-PV energy systems. The proposed stochastic framework is used to analyze whether hybrid renewable energy helps improve the available hosting capacity.
A variety of deployment schemes are generated by MCS. Figure 10 shows the probability distribution of HC results obtained with different shares of wind and PV at the 50% penetration level. It can be seen that the PDF curve of the single energy power generation hosting capacity is on the left side of the PDF curve of the hybrid energy hosting capacity. With the increase of wind power share, the PDF curve of the hosting capacity moves to the right until the share of wind power generation reaches 75%, and then the PDF curve of the hosting capacity moves to the left. It indicates that hybrid power generation has advantages in improving the level of the hosting capacity.
According to the results, when the wind power share reaches 75%, the hosting capacity of the hybrid wind-PV system reaches the maximum. Table 5 shows the mean values of hosting capacity and delivered generation. The mean value of the total hosting capacity is 9.396 MW, while the mean values of wind and PV hosting capacity are 7.040 and 2.356 MW, respectively. Compared with the results of the single wind case and PV case, the hosting capacity increases by about 1.22 times and 1.18 times, respectively. Meanwhile, the mean value of total annual delivered generation is 17,999.262 MWh, of which wind power generation accounts for about 87% and PV accounts for about 13%. Compared with the single wind case and PV case, the total delivered generation increases by about 1.04 times and 2.17 times, respectively.
Wind power plays a leading role in hybrid wind-PV systems. PV accounts for a relatively small proportion, but as a supplement to energy, it is also essential to increase the total hosting capacity and total delivered generation. The results show that the complementarity between wind power and PV is conducive to distribution networks to accommodate more distributed renewable resources. It can leverage more renewable generation capacity to be utilized, thereby promoting higher energy export.

CONCLUSION
Due to the inherent uncertainty of renewable energy resources, the hosting capacity in distributed networks is not an immutable value. Meanwhile, it is important to consider the uncertainties of penetration levels, locations of DGs, and shares of wind and PV. Thus, HC needs to be expressed in some statistical ways. In this study, based on the wind-photovoltaic-load temporal characteristic, a stochastic framework for the hosting capacity is proposed. The main conclusions are as follows: 1) The uncertainty of wind power, PV, and load demand is considered through time series data. The discretization-aggregation method is introduced to process time series data and generate extreme combinations. It reduces the number of scenarios to be evaluated and significantly mitigates computational complexity.
2) The holomorphic embedding model is proposed considering the direction of generation and load change. The equivalent analytical formula of voltage establishes the corresponding relationship between the actual operation level and the embedding parameter. The improved HEM can improve the efficiency of the hosting capacity assessment. 3) Hosting capacity of the wind system, PV system, and hybrid wind-PV system is studied from a probabilistic view. Compared with the single resource case, the hybrid case has the advantage in power generation. Due to the complementarity between wind power and PV, the hybrid wind-PV system has the potential to increase the hosting capacity and energy production in distributed networks. The performance in promoting energy integration and improving utilization varies according to different shares of wind and PV.
From the development trend of the low-carbon goal, a large amount of distributed renewable energy will inevitably lead to more significant changes in distribution networks. This study proposes a method to quantify the hosting capacity in distribution networks with DGs based on the holomorphic embedding method. It offers assistance in understanding the level of renewable energy generation, making better use of the available renewable resources, and promoting the application of distributed hybrid power generation in the power grid. Furthermore, we will combine the optimization algorithm to plan the optimal access scheme of distributed generation with the method proposed in this study.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
ND was responsible for methodology, formal analysis, validation, and editing of the manuscript. FT and QL were responsible for review and supervision. CW contributed to the conception and design of the study. XG and JX wrote sections of the manuscript. JZ and RL were responsible for part of data curation and project administration. All authors contributed, read, and approved the submitted version.