Microgrid power distribution system reliability evaluation method based on Bayesian network and sequential Monte Carlo method

Abstract

The invention discloses a reliability evaluation method for a micro-grid power distribution system based on a Bayesian network and a sequential Monte Carlo method, and the method comprises the steps: building a micro power supply multi-state combination output model according to the operation state of a power distribution element; carrying out load partitioning on the load in the microgrid through a multi-target particle swarm optimization algorithm; calculating a load output time sequence problem, and enabling the micro-grid distributed power supply to be equivalent to a piecewise function; carrying out analogue simulation on the wind speed and the illumination intensity by adopting an improved ARIMA (Autoregressive Integrated Moving Average), and establishing a time sequence composite model; calculating the ratio of the weekly load to the annual maximum load, the ratio of the daily load to the weekly maximum load, and the ratio of the load per hour to the daily maximum load; establishing a logic relation graph of each node through a Bayesian network, and establishing a time sequence source-load change model; in combination with the time sequence composite model and a time sequence source-load change model, sampling each node of a system element, and judging an operation state; and performing simulation analysis on the reliability index of the power distribution system to obtain a final result. The method can be used for improving the power supply performance by improving the reliability of the fixed-point weak link.

Description

Technical Field

[0001] The present invention relates to the technical field of power supply reliability assessment of a microgrid distribution system, and in particular to a reliability assessment method of a microgrid distribution system based on a Bayesian network and a sequential Monte Carlo method. Background Art

[0002] The reliability of the distribution system is an important indicator in power grid planning and is directly related to the user's power consumption experience. Traditional reliability assessment methods are mostly used for radial structures and are no longer applicable to new microgrids. Common methods for traditional power grid reliability assessment include the minimum path method, enumeration method, Markov, Monte Carlo method, and artificial intelligence method. These algorithms have gradually been improved and applied to the calculation of microgrid reliability, but it is difficult to calculate the key factors affecting reliability, which is not conducive to improving reliability in power grid planning. In addition, for complex multi-microgrid interconnected distribution network structures, the use of traditional methods is extremely complex and cumbersome.

[0003] There are some challenges of uncertainty in the reliability assessment of microgrid distribution systems. Distributed power sources in microgrids, especially renewable energy sources such as wind power and photovoltaics, have significant uncertainty and volatility. In addition, the current microgrid lacks a comprehensive evaluation model. Existing evaluation methods often focus on a single factor, such as equipment reliability or a single type of distributed power source, and lack an evaluation model that can comprehensively consider multiple factors in the microgrid. When the microgrid is disconnected from the main grid and can operate independently, that is, in the island operation mode, the reliability assessment in this operation mode is more complicated, and it is necessary to consider the dynamic matching of multiple factors such as power supply, load, and energy storage within the microgrid. At the same time, the uncertainty on the demand side or the implementation of demand response measures will increase the complexity of microgrid reliability assessment because it is necessary to consider the uncertainty of user behavior and demand changes.

[0004] In the prior art, a quantitative reliability assessment method based on uncertainty has been proposed, which takes into account the uncertainty on both the supply and demand sides and proposes a new reliability index to measure the risk of power shortage under uncertainty. In addition, an optimization-based method is proposed to calculate the standard network-related reliability indicators widely used in distribution network optimization models under reliability constraints. Distribution network reliability assessment methods based on optimization models, these methods attempt to improve the accuracy and efficiency of the assessment through mathematical optimization techniques. In addition, some studies have proposed reliability assessment methods that consider the island operation characteristics of microgrids. These methods capture the dynamic characteristics of microgrids by establishing more accurate mathematical models. With the development of technology, more innovative methods will be proposed in the future to address the shortcomings of existing technologies and provide more reliable guidance for the planning and operation of microgrids. Summary of the invention

[0005] In order to solve the above technical problems, the present invention provides a reliability assessment method for a microgrid distribution system based on a Bayesian network and a sequential Monte Carlo method. The method establishes a logical relationship between multiple microgrid network nodes based on a Bayesian network, judges the microgrid distribution operation mode based on the Monte Carlo sampling upper node state, and calculates and accumulates the load operation state in combination with the microgrid timing state; then calculates reliability indicators such as system outage time and outage frequency to evaluate the power supply reliability of the multi-microgrid system. Through case analysis, compared with the traditional power grid, the annual average outage rate and the duration of the outage of the microgrid load are significantly improved, and the system power supply reliability is significantly improved. According to the conditional probability of the node state of the event occurring by reverse reasoning based on the logical relationship, the weak link that affects the normal power consumption of the load is obtained, which can be used to improve the reliability of the fixed-point weak link to improve its power supply performance.

[0006] The technical solution adopted by the present invention is:

[0007] The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method includes the following steps:

[0008] Step 1: According to the operating status of the distribution components, a multi-state combined output model of the micro power source is established;

[0009] Step 2: Use the multi-objective particle swarm optimization algorithm to divide the loads in the microgrid into blocks according to the importance, electrical distance and load size as the objective function;

[0010] Step 3: Based on the load block calculation of step 2, the output state of the microgrid distributed power source within one hour can be equivalent to a piecewise function;

[0011] Step 4: Considering the influence of wind speed, light intensity and temperature, the improved ARIMA is used to simulate the wind speed and light intensity to establish a time series composite model;

[0012] Step 5: Based on the piecewise function of step 3 and the time series composite model of step 4, calculate the ratio of weekly load to annual maximum load, the ratio of daily load to its weekly maximum load, and the ratio of hourly load to its daily maximum load;

[0013] Step 6: Establish a logical relationship diagram of each node through the Bayesian network, infer the influence of the upper-level nodes on the system through reverse reasoning, and establish a model of time-series source-load changes to determine the weak links of the microgrid;

[0014] Step 7: Combining the time series composite model established in step 4 with the time series source-load variation model established in step 6, sampling each node of the system component to determine the operating status;

[0015] Step 8: Perform simulation analysis on the reliability index of the power distribution system to obtain the final result.

[0016] In step 1, the possible operating states of the distribution components are divided into three types, including normal operating state, maintenance shutdown state, and fault shutdown state; the failure probability of the distribution component in a year is λ f , when the fault causes downtime, the fault repair time μ f Normal operation is restored. The probability of a planned maintenance of a distribution component within one year is λ r , when the distribution components undergo planned maintenance for a period of μ r After that, the distribution components resume normal operation.

[0017] Since external factors will affect the probability of failure of distribution components, the expression of timing failure rate is obtained as follows:

[0018]

[0019] In formula (1), k0 is the influencing factor; t is the service life of the distribution component; t B is the running-in period of the component; c is the failure rate constant; λ is the average number of failures under normal conditions; ω is the proportional constant.

[0020] The step 2 comprises the following steps:

[0021] Step 2.1: Calculate the weights of different loads, use the same weights to calculate the same load, and express the difference in load cluster weights as:

[0022]

[0023] In formula (2), ω i is the weight when the load is i; is the expected weight of the load in the load block; s is the standard deviation value. When there is only one load in the load block, s = 0; n represents the sample size, that is, the number of observations in the data set.

[0024] Step 2.2: Calculate the electrical distance between different loads as follows:

[0025] When multiple loads in microgrids belong to the same microgrid, the smaller the load electrical distance is, the higher the possibility of being in the same cluster is. The electrical distance between loads is expressed as:

[0026]

[0027] In formula (3), Indicates WG i The distance between loads i and j in the region; S c Indicates the location of the microgrid connected to the grid, dist|Ldi ,S c | indicates load point Ld i Distance to the power grid; Ld i Indicates WG i Load point i in the region; Ld j Indicates WG i Load point j in the region; dist|Ld i ,Ld j | indicates WG i The distance between loads i and j in the region.

[0028] Step 2.3: Calculate the capacity allocation of different load blocks. The total amount of load blocks determines the amount of load reduction. The more detailed the load blocks are divided, the smaller the power shortage of the system. The objective function of the load block capacity is as follows:

[0029]

[0030] In formula (4), C L Represents a set of load blocks; L D Indicates the load block capacity.

[0031] Step 2.4: According to the multi-objective particle swarm algorithm, calculate the block situation of the load node under the constraint conditions.

[0032] Among them, the particle velocity and position are expressed as:

[0033]

[0034] In formula (5), x i represents the set of particles in the particle swarm; x1,x2,...,x n represents different particles in the particle group; v i represents the set of particle velocities in a particle group; v1,v2,...,v n Represents the speed of different particles in the particle group.

[0035] Taking into account the inertia weight, individual learning factor, group learning factor and the basis of the best historical position and speed, the position and speed of the next particle are determined, and the following expression is obtained:

[0036]

[0037] In formula (6), and are the particle velocity values ​​for the nth and n+1th iterations respectively; and are the position values ​​of the particles at the nth and n+1th iterations, respectively; ω is the inertia weight; γ1 and γ2 are the difference adjustment inertia weights; c1 and c2 are the individual learning factor and the group learning factor, respectively; and are the historical best positions of individual particles and group particles respectively.

[0038] In order to improve the calculation accuracy and convergence of MOPSO, the appropriate flight path of the particles is set, and the appropriate Gbest is used for optimization calculation; Gbest optimization first determines the similarity distance between the particles in the population and the particles in the external document, obtains a value less than the average similarity distance, and finds the corresponding non-dominated solution in the external document. The calculation process is as follows:

[0039]

[0040] In formula (7), d(x i ,y j ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document; SDi, are the similarity distance set and the average similarity distance respectively; x i,k represents the value of the kth dimension of the ith particle; y j.k represents the value of the kth dimension of the jth particle; n represents the sample size; d(x i ,y1) represents the distance between the i-th particle in the population and the first non-dominated solution of the external document; d(x i ,y2) represents the distance between the i-th particle in the population and the second non-dominated solution of the external document; d(x i ,y m ) represents the distance between the ith particle in the population and the mth non-dominated solution in the external document; m means the mth non-dominated solution; d(x i ,y i ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document.

[0041] Updated SD i :

[0042]

[0043] In formula (8), It means existence, indicating that d(x i ,y i ) exists in the value range of SDi.

[0044] The method based on crowding distance selects the global optimal solution Gbest. Its principle is to determine the optimal solution according to the sparseness and density of the distribution of the non-dominated solution set. iThe solution set is brought into the objective function to calculate the crowding distance of each solution, and the factor with the largest distance is selected as the global optimal solution Gbest.

[0045] Step 2.5: Based on the calculation results of the multi-objective particle swarm algorithm in step 2.4, the loads in the microgrid are clustered and divided into blocks, and the priorities of different loads are calculated and analyzed to ensure the stability of the microgrid during island operation.

[0046] In step 3, the output state of the microgrid within one hour is segmented into n stages, and the segmentation function is:

[0047]

[0048] In formula (9), x0 is the starting node of the first stage; x1 is the ending node of the first stage, that is, the starting node of the second stage; x2 is the ending node of the second stage, that is, the starting node of the third stage; x n Indicates the starting node of the nth stage; x n+1 Indicates the end node of the nth stage; x indicates the stage node between partitions. m represents the node at time m; f(x m ) represents the output state of the microgrid at time m. i represents the starting node at time i, x i+1 Indicates the end node at time i and the start node at the next time.

[0049] Since load, wind speed, and light collection are discrete functions, m represents the interval [x i ,x i+1 ) The number of samples, Indicates the average value.

[0050] In step 4, the time series composite model is an ARIMA (p, d, q) model, and the d-order difference processing obtains a stationary time series ARIMA (p, q), which satisfies the white noise sequence condition, wherein the difference expression is:

[0051]

[0052] In formula (10), d represents the difference order; Y t represents the value of the time series at time point t after difference processing; W t Represents the value of the original time series at time point t; W t-1 Represents the value of the original time series in the time period before time point t; W t-2 Represents the value of the original time series two time periods before time point t.

[0053] The expressions of ARIMA(p,d,q) model and ARIMA(p,q) model are:

[0054]

[0055]

[0056] Wherein, formula (11) is the ARIMA (p, d, q) model, and formula (12) is the ARIMA (p, q) model. Represents the d-order difference of the original time series, where is the difference operator, and d is the order of the difference. t is a random sequence W t The stationary sequence after d-order difference processing, where the criterion of the stationary sequence satisfies the p<0.05 corresponding to the ADF test value;

[0057] p and q are the autoregressive order and sliding order of the model. When q = 0, the ARIMA (p, q) model is converted to an AR model. When p = 0, the ARIMA (p, q) model is converted to an MA model.

[0058] are the parameters in the autoregressive model, where is a constant term, is the autoregression coefficient; Y t-1 represents the value at time t-1 in the past, Y t-p represents the value at the past time tp; θ1 represents the parameter of the moving average model, which represents the coefficient related to the past forecast error; a t-1 ,a t-q They represent the prediction errors of the time series at time point t-1 and time series tp respectively.

[0059] a t is an independent residual term,

[0060]

[0061] In formula (13), h t is the variance of the time series t; ε t is the standardized residual at time t, which follows a standard normal distribution.

[0062] Combining the parameters, the ARIMA (p, q) model prediction formula is:

[0063]

[0064] In formula (14), are the parameters of the autoregressive model, where i is the index of the autoregressive term; θ jis the average model parameter, where j is the index value of the moving average term; ε t-j Represents the prediction error term at time point tj, that is, the difference between the actual value at time point tj and the model predicted value.

[0065] According to the original data ARIMA (p, d, q) model, the prediction formula depends on the difference order d, so the prediction formula is as follows:

[0066]

[0067] In formula (15), W i represents the value of the original time series at time point t; Y represents the value of the time series after difference processing; W t-i Represents the value of the original time series at time point t-1; a i represents the difference operator at time i; a t-q represents the difference operator at time tq; W t-i-1 Represents the value of the original time series at time point ti-1.

[0068] In step 5, the load at a certain time t is:

[0069] L(t)=L y ·P w ·P d ·P h (t) (16);

[0070] In formula (16), L(t) represents the power load at time t; L y Indicates the original load value; P w It indicates the ratio of weekly load to annual maximum load; P d It indicates the ratio of daily load to the maximum weekly load; P h (t) represents the ratio of hourly load to the maximum load of the day.

[0071] 1) Based on historical load data, the ratio of weekly load to annual maximum load is calculated as:

[0072]

[0073] In formula (17), max(L y ) represents the annual peak load; L w (i) represents the load of the ith week; P w (i) represents the ratio of the load in the ith week to the annual peak value.

[0074] 2) Calculate the ratio of daily load to the weekly maximum load based on the weekly load changes:

[0075]

[0076] In formula (18), max(L w (i)) represents the weekly peak load of the ith week; L d (j) represents the load on the jth day of the ith week; P d (j) represents the ratio of daily load to the peak load of the week in which it is located.

[0077] 3) According to the change of daily load, calculate the ratio of hourly load to the maximum load of the day:

[0078]

[0079] In formula (19), L h (k) represents the load at the kth hour; max(L d (j)) represents the peak load on the jth day; P h (k) represents the ratio of the load in hour k to the peak load on that day.

[0080] In step 6, the relationship between each node is established according to the Bayesian network, and the weight of the upper node on a certain state is reversely inferred through the relationship between the upper node and the microgrid system to obtain the influence of the upper node on the system;

[0081] According to the total probability formula, the Bayesian formula is:

[0082]

[0083] In the formula, P(B i |A) means that under the condition that event A occurs, event B i The probability of occurrence is also called the posterior probability.

[0084] P(B i )>0, i=1,2,…,n represents the prior probability; P(A|B i ), i = 1, 2, 3, ..., n represents the posterior probability, n represents event B i The total number of categories, that is, all possible B i The number of events. Applied to the Bayesian network structure of the distribution system, event A can be represented as a load fault, event B can be represented as an upper-level component node, and the influence of the component state probability on the load node state, that is, the causal relationship.

[0085] Through this relationship, the known result state is reversed to the conditional probability caused by various reasons. The larger the value, the greater the impact on the result, so as to identify the weak part of the load. According to the node division in the distribution system, the corresponding set is reversed to the component node in turn, the component node state under the set is accumulated, and the relationship between the corresponding system operation state and the component state is established.

[0086] The model of time-series source-load change is established by using the wind speed and light intensity forecast data provided by the ARIMA model in step 4. Based on historical load data and predicted load demand, different time periods are predicted to meet the demand and obtain the load time-series model. Combined with the DG output power forecast, the model of time-series source-load change is obtained to analyze the matching between the output of power generation resources and the demand at different time points.

[0087] The output power prediction model of the solar photovoltaic system in the DG output power prediction model is as follows:

[0088]

[0089] Where: P PV (t) is the predicted output power at time t; η is the efficiency of the photovoltaic panel; G(t) is the solar irradiance at time t; A is the total area of ​​the photovoltaic panel; A sh is the area covered by the shadow; R soiling is the percentage reduction in performance due to fouling.

[0090] The output power prediction model of the solar wind system in the DG output power prediction model is as follows:

[0091]

[0092] Where: P wind (t) is the predicted output power at time t; ρ is the air density; A rotor is the swept area of ​​the wind turbine rotor; C p is the power coefficient, which depends on the design of the turbine and the wind speed; V(t) is the wind speed at time t; V ref is the reference wind speed, usually corresponding to the rated wind speed of the turbine.

[0093] In step 7, the distribution network system is initialized, and the logical relationship of each node of the Bayesian network is input; the sequential Monte Carlo sampling method is used to simulate the operation of the microgrid, and the model of the time series source-load change established by the time series model and the Bayesian network is used to sample the system components at time t. The state at a certain time is expressed as:

[0094] Xt i ={x1,x2,…,xk},k=1,2,…,m(23);

[0095] In formula (21), X ti t i The state set of the component nodes at the moment; x1, x2, …, x k Indicates that at time point t i The state of the system components is k, where k is the index of the component and m is the number of original components.

[0096] The node state of the sampled component is expressed as: k ∈[-1,0,1],x k =-1 indicates maintenance status, x k =0 means normal state, x k =1 indicates fault status.

[0097] Calculate the reliability-related parameters of each node of the system components to determine the operating status of the microgrid, including:

[0098] (1) According to t k The expression of load outage rate within a time period is:

[0099]

[0100] Where: In the time period t k The failure rate or outage rate of the i-th load point within . t k The number of outages of load i in the M samplings at time instant; In the time period t k The failure rate or downtime of the i-th load point within m a 、m b are the number of type a and type b failures respectively; r a 、r b They are the repair time for type A and type B faults respectively.

[0101] (2) The expression of the annual average reliability index of the load is obtained as follows:

[0102] Annual load failure outage rate λ of microgrid distribution system (unit: times / year):

[0103]

[0104] In formula (24): represents the annual load failure outage rate of the i-th load point; It represents the failure outage rate of the i-th load point in the k-th time period.

[0105] Microgrid distribution system load annual fault outage duration u (unit: h / year):

[0106]

[0107] In formula (25), u Lpi represents the annual fault outage duration of the i-th load point, that is, the total number of hours that the load point is out of service due to faults in one year; It represents the average failure downtime of the i-th load point in the k-th time period.

[0108] Average power outage duration r of microgrid distribution system load (unit: h / time):

[0109]

[0110] In formula (26), r Lpi It represents the average outage duration of the ith load point, that is, the average number of hours of outage for each fault.

[0111] The step 8 comprises the following steps:

[0112] Step 8.1: First, establish the following distribution system reliability indicators:

[0113] The average power outage frequency SAIFI is obtained from the number of power outages in the system in a year, unit: times / (household / year):

[0114]

[0115] In formula (27), N Lpi represents the number of users at the i-th load point.

[0116] The average outage duration SAIDI is obtained from the system's power outage duration in one year, unit: h / (household·year):

[0117]

[0118] The average power outage duration CAIDI is obtained from the power outage duration of each household in one year, unit: h / (household·year):

[0119]

[0120] From this, we can get the average power supply availability rate ASAI and the average power supply unavailability rate ASUI indicators:

[0121]

[0122] Expected power shortage index EENS, unit: MW / (household / year):

[0123]

[0124] In the above formula, N Lpi 、M i are the number of users at load point i and the number of users forced to shut down; P Lpi (t k ) is the real-time load of load point i at time k; R is the set of load points in the distribution system.

[0125] Step 8.2: Perform specific simulation analysis in the test system.

[0126] The BUS6-F4 of RBTS is used as the test system. Three independent microgrids including wind power, photovoltaic power, energy storage device, diesel generator and other micro power sources are added. The operation status of each node in the whole year is extracted, the load changes throughout the year are accumulated, and the power supply reliability index is obtained.

[0127] Step 8.3: Analyze the simulation system through the model obtained in step 8.2, and calculate the improved results in combination with the reliability index in step 8.1.

[0128] The present invention provides a reliability assessment method for a microgrid distribution system based on a Bayesian network and a sequential Monte Carlo method, and the technical effects are as follows:

[0129] 1) Step 1 of the present invention can comprehensively consider the characteristics and operating states of different power supplies, provide accurate input data for subsequent optimization and analysis, and the multi-state combination includes multiple states such as power supply startup, shutdown, and failure, providing more detailed power supply management;

[0130] 2) Step 2 of the present invention can consider multiple objective functions at the same time through multi-objective optimization to achieve more reasonable load distribution. The particle swarm optimization algorithm provides a new search strategy to improve the search efficiency and solution quality.

[0131] 3) Step 3 of the present invention equates the output state to a piecewise function, which can simplify the calculation process while maintaining sufficient accuracy. This method provides a new output state modeling method, which helps to better predict output changes.

[0132] 4) Step 4 of the present invention takes into account key environmental factors that affect the output of distributed power sources, improves the accuracy of the model, and the improved simulation includes a more accurate physical model, thereby improving the reliability of the simulation.

[0133] 5) Step 5 of the present invention provides a new load analysis method, which helps to optimize load management and prediction, and provides a new load analysis method, which helps to optimize load management and prediction.

[0134] 6) The Bayesian network in step 6 of the present invention can handle uncertainty and probabilistic problems, which helps to identify and analyze the weak links of the system. The use of reverse reasoning provides a new method for logical relationship analysis, which helps to deeply understand the dynamic behavior of the system.

[0135] 7) Step 7 of the present invention combines the timing model and the Bayesian network to more comprehensively analyze the operating status of the system, providing a new system status analysis framework that is helpful for real-time monitoring and fault diagnosis.

[0136] 8) Step 8 of the present invention can evaluate the reliability of the system without actually interrupting the service through simulation analysis. This method includes new reliability indicators or simulation technologies and provides a more comprehensive system evaluation. BRIEF DESCRIPTION OF THE DRAWINGS

[0137] Figure 1 Schematic diagram of the component state transition model.

[0138] Figure 2 This is the flow chart of the MOPSO algorithm.

[0139] Figure 3 This is the optimized load connection status diagram.

[0140] Figure 4 It is a schematic diagram of piecewise function.

[0141] Figure 5 This is the ARIMA model flow chart.

[0142] Figure 6 Build a relationship graph of nodes for the Bayesian network.

[0143] Figure 7 Flowchart for reliability assessment of multi-microgrid distribution system.

[0144] Figure 8 To improve the feeder network structure diagram.

[0145] Fig. 9 This is the Bayesian network diagram corresponding to the non-microgrid system.

[0146] Fig.10 This is the Bayesian network diagram corresponding to the non-microgrid system.

[0147] Fig.11 It is a bar chart of the average annual failure outage rate of load points. DETAILED DESCRIPTION

[0148] A reliability assessment method for microgrid distribution system based on Bayesian network and sequential Monte Carlo method. This method establishes the logical relationship between multiple microgrid network nodes through Bayesian network, judges the microgrid distribution operation mode based on the upper node state of Monte Carlo sampling, calculates and accumulates the load operation state in combination with the microgrid timing state, and then calculates reliability indicators such as system outage time and outage frequency to evaluate the power supply reliability of the multi-microgrid system. Through case analysis, compared with the traditional power grid, the annual average outage rate and the duration of the outage of the microgrid load are significantly improved, and the system power supply reliability is significantly improved. According to the logical relationship, the conditional probability of the node state of the reverse reasoning event is obtained to obtain the weak link that affects the normal power consumption of the load, which can be used to improve the reliability of the fixed-point weak link to improve its power supply performance. Specifically, it includes the following steps:

[0149] Step 1: According to the operating status of the distribution components, a multi-state combined output model of the micro power source is established;

[0150] Step 2: Use the multi-objective particle swarm optimization algorithm to divide the loads in the microgrid into blocks according to the importance, electrical distance and load size as the objective function;

[0151] Step 3: Based on the load block calculation of step 2, the output state of the microgrid distributed power source within one hour can be equivalent to a piecewise function;

[0152] Step 4: Considering the influence of wind speed, light intensity and temperature, the improved ARIMA is used to simulate the wind speed and light intensity to establish a time series composite model;

[0153] Step 5: Based on the piecewise function of step 3 and the time series composite model of step 4, calculate the ratio of weekly load to annual maximum load, the ratio of daily load to its weekly maximum load, and the ratio of hourly load to its daily maximum load;

[0154] Step 6: Establish a logical relationship diagram of each node through the Bayesian network, infer the influence of the upper-level nodes on the system through reverse reasoning, and establish a model of time-series source-load changes to determine the weak links of the microgrid;

[0155] Step 7: Combining the time series composite model established in step 4 with the time series source-load variation model established in step 6, sampling each node of the system component to determine the operating status;

[0156] Step 8: Perform simulation analysis on the reliability index of the power distribution system to obtain the final result.

[0157] In step 1, the possible operating states of the distribution components are divided into three types, including normal operating state, maintenance shutdown state, and fault shutdown state; Figure 1 As shown:

[0158] Where: f is the annual failure rate, that is, how many failures occur in a year; μ f is the annual fault repair time, that is, how many hours are spent repairing faults in a year. r is the annual maintenance rate, i.e. the number of maintenance times in a year; μ r It is the annual maintenance time, that is, how many hours a year the equipment is maintained.

[0159] like Figure 1 As shown, the failure probability of the distribution component in one year is λ f , when the fault causes downtime, the fault repair time μ f Normal operation is restored. The probability of a planned maintenance of a distribution component within one year is λr , when the distribution components undergo planned maintenance for a period of μ r After that, the distribution components resume normal operation.

[0160] Since external factors will affect the probability of failure of distribution components, the expression of timing failure rate can be obtained as follows:

[0161]

[0162] In formula (1), k0 is the influencing factor; t is the service life of the distribution component; t B is the running-in period of the component; c is the failure rate constant; λ is the average number of failures under normal conditions; ω is the proportional constant, which is affected by the type of failure.

[0163] The step 2 comprises the following steps:

[0164] Step 2.1: Calculate the weights of different loads, use the same weights to calculate the same load, and express the difference in load cluster weights as:

[0165]

[0166] In formula (2), ω i is the weight when the load is i; is the expected weight of the load in the load block; s is the standard deviation value. When there is only one load in the load block, s = 0; n represents the sample size, that is, the number of observations in the data set.

[0167] Step 2.2: Calculate the electrical distance between different loads as follows:

[0168] When multiple loads in microgrids belong to the same microgrid, the smaller the load electrical distance is, the higher the possibility of being in the same cluster is. The electrical distance between loads can be expressed as:

[0169]

[0170] In formula (3), Indicates WG i The distance between loads i and j in the region; S c Indicates the location of the microgrid connected to the grid, dist|Ld i ,S c | indicates load point Ld i Distance to the power grid; Ld i Indicates WG i Load point i in the region; Ld j Indicates WG i Load point j in the region; dist|Ld i,Ld j | indicates WG i The distance between loads i and j in the region.

[0171] Step 2.3: Calculate the capacity allocation of different load blocks. The load blocks are isolated from the distribution network system through intelligent switches. The total amount of load blocks determines the amount of load reduction. The more detailed the load blocks are divided, the smaller the system power shortage is. The objective function of the load block capacity is as follows:

[0172]

[0173] In formula (4), C L Represents a set of load blocks; L D Indicates the load block capacity.

[0174] Step 2.4: According to the multi-objective particle swarm algorithm, calculate the block situation of the load node under the constraint conditions.

[0175] Among them, the particle velocity and position can be expressed as:

[0176]

[0177] In formula (5), x i represents the set of particles in the particle swarm; x1,x2,...,x n represents different particles in the particle group; v i represents the set of particle velocities in a particle group; v1,v2,...,v n Represents the speed of different particles in the particle group.

[0178] In order to speed up the convergence of the algorithm and avoid the algorithm result being a local optimum, the inertia weight, individual learning factor, group learning factor and the basis of the historical best position and speed are comprehensively considered to determine the position and speed of the next particle, and the following expression is obtained:

[0179]

[0180] In formula (6), and are the particle velocity values ​​for the nth and n+1th iterations respectively; and are the position values ​​of the particles at the nth and n+1th iterations, respectively; ω is the inertia weight; γ1 and γ2 are the difference adjustment inertia weights; c1 and c2 are the individual learning factor and the group learning factor, respectively; and are the historical best positions of individual particles and group particles respectively.

[0181] In order to avoid the result being a local optimum, the convergence of the algorithm is accelerated in many aspects. In order to improve the calculation accuracy and convergence of MOPSO, the appropriate flight path of the particles is set, and the appropriate Gbest is used for optimization calculation; Gbest optimization first determines the similarity distance between the particles in the population and the particles in the external document, obtains a value less than the average similarity distance, and finds the corresponding non-dominated solution in the external document. The calculation process is as follows:

[0182]

[0183] In formula (7), d(x i ,y j ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document; SDi, are the similarity distance set and the average similarity distance respectively; x i,k represents the value of the kth dimension of the ith particle; y j.k represents the value of the kth dimension of the jth particle; n represents the sample size; d(x i ,y1) represents the distance between the i-th particle in the population and the first non-dominated solution of the external document; d(x i ,y2) represents the distance between the i-th particle in the population and the second non-dominated solution of the external document; d(x i ,y m ) represents the distance between the ith particle in the population and the mth non-dominated solution in the external document; m means the mth non-dominated solution; d(x i ,y i ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document.

[0184] Updated SD i :

[0185]

[0186] In formula (8), It means existence, indicating that d(x i ,y i ) exists in the value range of SDi.

[0187] The method based on crowding distance selects the global optimal solution Gbest. Its principle is to determine the optimal solution according to the sparseness and density of the distribution of the non-dominated solution set. i The solution set is brought into the objective function to calculate the crowding distance of each solution, and the factor with the largest distance is selected as the global optimal solution Gbest. The algorithm flow is as follows Figure 2 shown.

[0188] Initialize the particle population Gbest, input the particle Pbest, speed and position, set the distribution network structure, number of iterations, termination conditions and other related parameters.

[0189] Step 2.5: Based on the calculation results of the multi-objective particle swarm algorithm in step 2.4, cluster the loads in the microgrid, calculate and analyze the priorities of different loads, and ensure the stability of the microgrid during island operation. Taking the IEEE-RBTSBUS6 feeder network as an example, the loads in the microgrid area are divided into blocks. Figure 3 As shown. Figure 3 It can be seen that equipping both ends of the distribution line with section switches to ensure that the fault is removed when a line fault occurs can reduce the non-fault power outage area, and the optimized load block increases the switch division load, improving the power supply reliability of the distribution network.

[0190] The load shedding priority calculated in combination with the electrical distance between the micro power source can be expressed as shown in Table 1:

[0191] Table 1 Load shedding priority

[0192] Load point Priority LP11,LP14,LP19 1 LP12,LP15,LP22 2 LP13,LP16,LP23 3 LP17-LP18,LP20-LP21 4

[0193] In step 3, the output state of the microgrid within one hour is segmented into n stages, and the segmentation function is:

[0194]

[0195] In formula (9), x0 is the starting node of the first stage; x1 is the ending node of the first stage, that is, the starting node of the second stage; x2 is the ending node of the second stage, that is, the starting node of the third stage; x n Indicates the starting node of the nth stage; x n+1 Indicates the end node of the nth stage; x indicates the stage node between partitions. m represents the node at time m; f(x m ) represents the output state of the microgrid at time m. i represents the starting node at time i, x i+1 Indicates the end node at time i and the start node at the next time.

[0196] Since load, wind speed, and light collection are discrete functions, m represents the interval [x i ,x i+1 ) The number of samples, represents the average value. Its piecewise function diagram is as follows Figure 4 shown.

[0197] In step 4, the time series composite model is an ARIMA (p, d, q) model, and the d-order difference processing obtains a stationary time series ARIMA (p, q), which satisfies the white noise sequence condition, wherein the difference expression is:

[0198]

[0199] In formula (10), d represents the difference order; Y t represents the value of the time series at time point t after difference processing; W t Represents the value of the original time series at time point t; W t-1 Represents the value of the original time series in the time period before time point t; W t-2 Represents the value of the original time series two time periods before time point t.

[0200] The expressions of ARIMA(p,d,q) model and ARIMA(p,q) model are:

[0201]

[0202]

[0203] Wherein, formula (11) is the ARIMA (p, d, q) model, and formula (12) is the ARIMA (p, q) model. Represents the d-order difference of the original time series, where is the difference operator, and d is the order of the difference. t is a random sequence W t The stationary sequence after d-order difference processing, where the criterion of the stationary sequence satisfies the p<0.05 corresponding to the ADF test value;

[0204] p and q are the autoregressive order and sliding order of the model. When q = 0, the ARIMA (p, q) model is converted to an AR model. When p = 0, the ARIMA (p, q) model is converted to an MA model.

[0205] and θ are unknown coefficients that are not zero; are the parameters in the autoregressive model, where is a constant term, is the autoregression coefficient; Y t-1 represents the value at time t-1 in the past, Y t-p Represents the value at time tp in the past;

[0206] θ1 represents the parameter of the moving average model, which represents the coefficient related to the past forecast error; a t-1 ,a t-qThey represent the prediction errors of the time series at time point t-1 and time series tp respectively.

[0207] a t is an independent residual term,

[0208]

[0209] In formula (13), h t is the variance of the time series t; ε t is the standardized residual at time t, which follows a standard normal distribution.

[0210] Combining the parameters, the ARIMA (p, q) model prediction formula is:

[0211]

[0212] In formula (14), are the parameters of the autoregressive model, where i is the index of the autoregressive term; θ j is the average model parameter, where j is the index value of the moving average term; ε t-j Represents the prediction error term at time point tj, that is, the difference between the actual value at time point tj and the model predicted value.

[0213] According to the original data ARIMA (p, d, q) model, the prediction formula depends on the difference order d, so the prediction formula is as follows:

[0214]

[0215] In formula (15), W i represents the value of the original time series at time point t; Y represents the value of the time series after difference processing; W t-i Represents the value of the original time series at time point t-1; a i represents the difference operator at time i; a t-q represents the difference operator at time tq; W t-i-1 Represents the value of the original time series at time point ti-1.

[0216] The ARIMA model building process is as follows Figure 5 shown.

[0217] In step 5, the load at a certain time t is:

[0218] L(t)=L y ·P w ·P d ·P h (t) (16);

[0219] In formula (16), L(t) represents the power load at time t; L y Indicates the original load value; P w Indicates the ratio of weekly load to annual maximum load; P d It indicates the ratio of daily load to the maximum weekly load; P h (t) represents the ratio of hourly load to the maximum load of the day.

[0220] 1) Based on historical load data, the ratio of weekly load to annual maximum load is calculated as:

[0221]

[0222] In formula (17), max(L y ) represents the annual peak load; L w (i) represents the load of the ith week; P w (i) represents the ratio of the load in the ith week to the annual peak value.

[0223] 2) Calculate the ratio of daily load to the weekly maximum load based on the weekly load changes:

[0224]

[0225] In formula (18), max(L w (i)) represents the weekly peak load of the ith week; L d (j) represents the load on the jth day of the ith week; P d (j) represents the ratio of daily load to the peak load of the week in which it is located.

[0226] 3) According to the change of daily load, calculate the ratio of hourly load to the maximum load of the day:

[0227]

[0228] In formula (19), L h (k) represents the load at the kth hour; max(L d (j)) represents the peak load on the jth day; P h (k) represents the ratio of the load in hour k to the peak load on that day.

[0229] In step 6, a relationship diagram of each node is established according to the Bayesian network, such as Figure 6 As shown in the figure, through the relationship between the upper node and the microgrid system, the weight of the upper node on a certain state can be reversely inferred to obtain the influence of the upper node on the system;

[0230] According to the total probability formula, the Bayesian formula is:

[0231]

[0232] In the formula, P(B i |A) means that under the condition that event A occurs, event B i The probability of occurrence is also called the posterior probability.

[0233] P(B i )>0, i=1,2,…,n represents the prior probability; P(A|B i ), i = 1, 2, 3, ..., n represents the posterior probability, n represents event B i The total number of categories, that is, all possible B i The number of events. Applied to the Bayesian network structure of the distribution system, event A can be represented as a load fault, event B can be represented as an upper-level component node, and the influence of the component state probability on the load node state, that is, the causal relationship.

[0234] Through this relationship, the known result state is reversed to the conditional probability caused by various reasons. The larger the value, the greater the impact on the result, so as to identify the weak part of the load. According to the previous introduction, the nodes in the distribution system are divided into corresponding sets, which are reversed to the component nodes in turn, and the component node states under the accumulated sets are established to establish the relationship between the corresponding system operation state and the component state.

[0235] When the status of the main components in the power system changes, it will affect the safe operation of the system. For example, when the generator fails, the power supply capacity decreases and the system stability decreases; when the transformer is overloaded, it may cause equipment damage or fire, and when the transformer fails, it may cause power outages; when the circuit breaker fails, the scope of the fault power outage will increase. Through the failure of different components, or different fault manifestations, the relationship between the component status and the distribution system status can be analyzed, and then the impact of each link in the distribution system on the system power supply reliability can be obtained.

[0236] The model of time-series source-load change is established by using the wind speed and light intensity forecast data provided by the ARIMA model in step 4. Based on historical load data and predicted load demand, different time periods are predicted to meet the demand and obtain the load time-series model. Combined with the following DG output power forecast, the model of time-series source-load change is obtained to analyze the matching between the output of power generation resources and the demand at different time points.

[0237] The output power prediction model of the solar photovoltaic system in the DG output power prediction model is as follows:

[0238]

[0239] Where: P PV(t) is the predicted output power at time t; η is the efficiency of the photovoltaic panel; G(t) is the solar irradiance at time t; A is the total area of ​​the photovoltaic panel; A sh is the area covered by the shadow; R soiling is the percentage reduction in performance due to fouling.

[0240] The output power prediction model of the solar wind system in the DG output power prediction model is as follows:

[0241]

[0242] Where: P wind (t) is the predicted output power at time t; ρ is the air density; A rotor is the swept area of ​​the wind turbine rotor; C p is the power coefficient, which depends on the design of the turbine and the wind speed; V(t) is the wind speed at time t; V ref is the reference wind speed, usually corresponding to the rated wind speed of the turbine.

[0243] In step 7, the distribution network system is initialized, and the logical relationship of each node of the Bayesian network is input; the sequential Monte Carlo sampling method is used to simulate the operation of the microgrid, and the model of the time series source-load change established by the time series model and the Bayesian network is used to sample the system components at time t. The state at a certain time can be expressed as:

[0244] Xt i ={x1,x2,…,xk},k=1,2,…,m(23);

[0245] In formula (21), X ti t i The state set of the component nodes at the moment; x1, x2, …, x k Indicates that at time point t i The state of the system components is k, where k is the index of the component and m is the number of original components.

[0246] The node state of the sampled component can be expressed as: k ∈[-1,0,1], x k =-1 indicates maintenance status, x k =0 means normal state, x k =1 indicates fault status.

[0247] The reliability assessment process of the microgrid distribution system can be as follows: Figure 7 shown.

[0248] Calculate the reliability-related parameters of each node of the system components to determine the operating status of the microgrid, including:

[0249] (1) According to t k The expression of load outage rate within a time period is:

[0250]

[0251]

[0252] Where: In the time period t k The failure rate or outage rate of the i-th load point within . t k The number of outages of load i in the M samplings at time; In the time period t k The failure rate or downtime of the i-th load point within m a 、m b are the number of type a and type b failures respectively; r a 、r b They are the repair time for type A and type B faults respectively.

[0253] (2) The expression of the annual average reliability index of the load is obtained as follows:

[0254] Annual load failure outage rate λ of microgrid distribution system (unit: times / year):

[0255]

[0256] In formula (24): Lpi represents the annual load failure outage rate of the i-th load point; It represents the failure outage rate of the i-th load point in the k-th time period.

[0257] Microgrid distribution system load annual fault outage duration u (unit: h / year):

[0258]

[0259] In formula (25), u Lpi represents the annual fault outage duration of the i-th load point, that is, the total number of hours that the load point is out of service due to faults in one year; It represents the average failure downtime of the i-th load point in the k-th time period.

[0260] Average power outage duration r of microgrid distribution system load (unit: h / time):

[0261]

[0262] In formula (26), r Lpi It represents the average outage duration of the ith load point, that is, the average number of hours of outage for each fault.

[0263] The step 8 comprises the following steps:

[0264] Step 8.1: First, establish the following distribution system reliability indicators:

[0265] The average power outage frequency SAIFI is obtained from the number of power outages in the system in a year, unit: times / (household / year):

[0266]

[0267] In formula (27), N Lpi represents the number of users at the i-th load point.

[0268] The average outage duration SAIDI is obtained from the system's power outage duration in one year, unit: h / (household·year):

[0269]

[0270] The average power outage duration CAIDI is obtained from the power outage duration of each household in one year, unit: h / (household·year):

[0271]

[0272] From this, we can get the average power supply availability rate ASAI and the average power supply unavailability rate ASUI indicators:

[0273]

[0274] Expected power shortage index EENS, unit: MW / (household / year):

[0275]

[0276] In the above formula, N Lpi 、M i are the number of users at load point i and the number of users forced to shut down; P Lpi (t k ) is the real-time load of load point i at time k; R is the set of load points in the distribution system.

[0277] Step 8.2: Conduct specific simulation analysis on the established sequential source-load variation model in the test system.

[0278] The BUS6-F4 of RBTS is used as the test system. Three independent microgrids including wind power, photovoltaic power, energy storage device, diesel generator and other micro power sources are added, and the BUS6-F4 feeder is improved as follows: Figure 8 The operating status of each node in the whole year is extracted, the load changes throughout the year are accumulated, and the power supply reliability index is obtained.

[0279] like Fig. 9 As shown, it is the Bayesian network corresponding to the non-microgrid system;

[0280] like Fig.10 As shown, it is the Bayesian network corresponding to the microgrid system;

[0281] Step 8.3: Analyze the simulation system through the model obtained in step 8.2, and calculate the improved results in combination with the reliability index in step 8.1.

[0282] The whole year is divided into 8760 hours, the operation status of the distribution system in each time period is simulated, the probability of failure shutdown and shutdown time in a certain time period are accumulated, and the corresponding power supply reliability index is calculated. As shown in Table 2.

[0283] Table 2 Load point reliability

[0284]

[0285] According to Table 2, a bar chart of the load annual failure outage rate is established as follows: Fig.11 As shown. It can be seen intuitively that the average outage rate of load points in the microgrid is significantly lower than that of the traditional power grid. According to the bar chart, the average annual outage rate distribution of LP19-LP23 loads in the microgrid WG1, LP20 and LP21 average annual outage rates are significantly higher than LP19, LP22 and LP23. The reason is that when the microgrid is islanded, the output power of the micro-source cannot meet the power supply needs of all loads, so some loads are removed according to the load priority to ensure the uninterrupted power supply of important loads. LP20 and LP21 have the highest removal priority, so the average annual outage rate of this load is significantly higher than other loads.

[0286] The calculation results of the power supply reliability index of the distribution system are shown in Table 3 below.

[0287] Table 3 System power supply reliability indicators

[0288]

[0289]

[0290] As shown in Table 3, the power supply reliability index ASAI of the single microgrid system WG1, WG2, and WG3 at the end of the distribution network is higher than that of the traditional power grid structure near the distribution line, and the system reliability indexes SAIFI, SAIDI, EENS, and ASUI are lower than some traditional power grid distribution systems; after the improvement of the entire distribution system of feeder F4, the power supply reliability of the distribution system containing microgrids is higher than that of the traditional distribution system without microgrids. Due to the existence of micro power sources, the microgrid system is like adding a backup power source in the microgrid, and the main grid fault can still be supplied by distributed power sources. The microgrid distribution system uses clean renewable energy to generate electricity without polluting the environment, while improving the system power supply reliability and increasing the user's good experience in the power consumption process.

[0291] According to algorithm analysis, compared with traditional power grids, the annual average outage rate and duration of microgrid load outages have decreased significantly, and the power supply reliability of the system has been significantly improved; reverse reasoning based on the Bayesian network logical relationship can be used to specifically study the weights of factors affecting the power supply reliability of important loads. This method can be used to improve the reliability of weak links and improve their power supply performance.

Claims

1. A reliability assessment method for microgrid distribution system based on Bayesian network and sequential Monte Carlo method, characterized by The following steps are involved: Step 1: According to the operating status of the distribution components, a multi-state combined output model of the micro power source is established; Step 2: Use the multi-objective particle swarm optimization algorithm to divide the loads in the microgrid into blocks according to the importance, electrical distance and load size as the objective function; Step 3: Based on the load block calculation of step 2, the output state of the microgrid distributed power source within one hour is equivalent to a piecewise function; Step 4: Considering the influence of wind speed, light intensity and temperature, the improved ARIMA is used to simulate the wind speed and light intensity to establish a time series composite model; Step 5: Based on the piecewise function of step 3 and the time series composite model of step 4, calculate the ratio of weekly load to annual maximum load, the ratio of daily load to its weekly maximum load, and the ratio of hourly load to its daily maximum load; Step 6: Establish a logical relationship diagram of each node through the Bayesian network, infer the influence of the upper-level nodes on the system through reverse reasoning, and establish a model of time-series source-load changes to determine the weak links of the microgrid; Step 7: Combining the time series composite model established in step 4 with the time series source-load variation model established in step 6, sampling each node of the system component to determine the operating status; Step 8: Perform simulation analysis on the reliability index of the power distribution system to obtain the final result.

2. The reliability assessment method of a microgrid distribution system based on a Bayesian network and a sequential Monte Carlo method according to claim 1 is characterized in that: In step 1, the possible operating states of the distribution components are divided into three types, including normal operating state, maintenance shutdown state, and fault shutdown state; the failure probability of the distribution component in a year is λ f , when the fault causes downtime, the fault repair time μ f Resume normal operations; The probability of a planned maintenance of a distribution component within one year is λ r , when the distribution components undergo planned maintenance for a period of μ r After that, the power distribution components resume normal operation; Since external factors will affect the probability of failure of distribution components, the expression of timing failure rate is obtained as follows: In formula (1), k0 is the influencing factor; t is the service life of the distribution component; t B is the running-in period of the component; c is the failure rate constant; λ is the average number of failures under normal conditions; ω is the proportional constant.

3. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: The step 2 comprises the following steps: Step 2.1: Calculate the weights of different loads, use the same weights to calculate the same load, and express the difference in load cluster weights as: In formula (2), ω i is the weight when the load is i; is the expected weight of the load in the load block; s is the standard deviation value. When there is only one load in the load block, s = 0; n represents the sample size, that is, the number of observations in the data set; Step 2.2: Calculate the electrical distance between different loads. The calculation method is as follows: When multiple loads in microgrids belong to the same microgrid, the smaller the load electrical distance is, the higher the possibility of being in the same cluster is. The electrical distance between loads is expressed as: In formula (3), Indicates WG i The distance between loads i and j in the region; S c Indicates the location of the microgrid connected to the grid, dist|Ld i ,S c | indicates load point Ld i Distance to the power grid; Ld i Indicates WG i Load point i in the region; Ld j Indicates WG i Load point j in the region; dist|Ld i ,Ld j | indicates WG i The distance between loads i and j in the region; Step 2.3: Calculate the capacity allocation of different load blocks. The total amount of load blocks determines the amount of load reduction. The more detailed the load blocks are divided, the smaller the power shortage of the system. The objective function of the load block capacity is as follows: In formula (4), C L Represents a set of load blocks; L D Indicates the load block capacity; Step 2.4: According to the multi-objective particle swarm algorithm, calculate the block situation of the load node under the constraint conditions. Among them, the particle velocity and position are expressed as: In formula (5), x i represents the set of particles in the particle swarm; x1,x2,...,x n represents different particles in the particle group; v i represents the set of particle velocities in a particle group; v1,v2,...,v n Indicates the speed of different particles in the particle group; Taking into account the inertia weight, individual learning factor, group learning factor and the basis of the best historical position and speed, the position and speed of the next particle are determined, and the following expression is obtained: In formula (6), and are the particle velocity values ​​for the nth and n+1th iterations respectively; and are the position values ​​of the particles at the nth and n+1th iterations, respectively; ω is the inertia weight; γ1 and γ2 are the difference adjustment inertia weights; c1 and c2 are the individual learning factor and the group learning factor, respectively; and are the historical best positions of individual particles and group particles respectively; In order to improve the calculation accuracy and convergence of MOPSO, the appropriate flight path of the particles is set, and the appropriate Gbest is used for optimization calculation; Gbest optimization first determines the similarity distance between the particles in the population and the particles in the external document, obtains a value less than the average similarity distance, and finds the corresponding non-dominated solution in the external document. The calculation process is as follows: In formula (7), d(x i ,y j ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document; SDi, are the similarity distance set and the average similarity distance respectively; x i,k represents the value of the kth dimension of the ith particle; y j.k represents the value of the kth dimension of the jth particle; n represents the sample size; d(x i ,y1) represents the distance between the i-th particle in the population and the first non-dominated solution of the external document; d(x i ,y2) represents the distance between the i-th particle in the population and the second non-dominated solution of the external document; d(x i ,y m ) represents the distance between the ith particle in the population and the mth non-dominated solution in the external document; m means the mth non-dominated solution; d(x i ,y i ) represents the distance between the i-th particle in the population and the j-th non-dominated solution in the external document; Updated SD i : In formula (8), It means existence, indicating that d(x i ,y i ) exists in the range of SD i middle; Select the global optimal solution Gbest, the principle of which is to determine the optimal solution based on the sparseness and density of the non-dominated solution set distribution; update the SD i The solution set is brought into the objective function to calculate the crowding distance of each solution, and the factor with the largest distance is selected as the global optimal solution Gbest; Step 2.5: Based on the calculation results of the multi-objective particle swarm algorithm in step 2.4, the loads in the microgrid are clustered and divided into blocks, and the priorities of different loads are calculated and analyzed to ensure the stability of the microgrid during island operation.

4. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: In step 3, the output state of the microgrid within one hour is segmented into n stages, and the segmentation function is: In formula (9), x0 is the starting node of the first stage; x1 is the ending node of the first stage, that is, the starting node of the second stage; x2 is the ending node of the second stage, that is, the starting node of the third stage; x n Indicates the starting node of the nth stage; x n+1 Indicates the end node of the nth stage; x indicates the phased node between partitions; x m represents the node at time m; f(x m ) represents the output state of the microgrid at time m; x i represents the starting node at time i, x i+1 Indicates the end node at time i and the start node at the next time; Since load, wind speed, and light collection are discrete functions, m represents the interval [x i ,x i+1 ) The number of samples, Indicates the average value.

5. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: In step 4, the time series composite model is an ARIMA (p, d, q) model, and the d-order difference processing obtains a stationary time series ARIMA (p, q), which satisfies the white noise sequence condition, wherein the difference expression is: In formula (10), d represents the difference order; Y t represents the value of the time series at time point t after difference processing; W t Represents the value of the original time series at time point t; W t-1 Represents the value of the original time series in the time period before time point t; W t-2 Represents the value of the original time series in the two time periods before time point t; The expressions of ARIMA(p,d,q) model and ARIMA(p,q) model are: Wherein, formula (11) is the ARIMA (p, d, q) model, and formula (12) is the ARIMA (p, q) model; Represents the d-order difference of the original time series, where is the difference operator, d is the order of the difference; Y t is a random sequence W t The stationary sequence after d-order difference processing, where the criterion of the stationary sequence satisfies the p<0.05 corresponding to the ADF test value; p and q are the autoregressive order and sliding order of the model. When q = 0, the ARIMA (p, q) model is converted to an AR model. When p = 0, the ARIMA (p, q) model is converted to an MA model. are the parameters in the autoregressive model, where is a constant term, is the autoregression coefficient; Y t-1 represents the value at time t-1 in the past, Y t-p represents the value at the past time tp; θ1 represents the parameter of the moving average model, which represents the coefficient related to the past forecast error; a t-1 ,a t-q They represent the prediction errors of the time series at time point t-1 and time series tp respectively; a t is an independent residual term, In formula (13), h t is the variance of the time series t; ε t is the standardized residual at time t, which obeys the standard normal distribution; Combining the parameters, the ARIMA (p, q) model prediction formula is: In formula (14), are the parameters of the autoregressive model, where i is the index of the autoregressive term; θ j is the average model parameter, where j is the index value of the moving average term; ε t-j Represents the prediction error term at time point tj, that is, the difference between the actual value at time point tj and the model predicted value; According to the original data ARIMA (p, d, q) model, the prediction formula depends on the difference order d, so the prediction formula is as follows: In formula (15), W i represents the value of the original time series at time point t; Y represents the value of the time series after difference processing; W t-i Represents the value of the original time series at time point t-1; a i represents the difference operator at time i; a t-q represents the difference operator at time tq; W t-i-1 Represents the value of the original time series at time point ti-1.

6. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: In step 5, the load at a certain time t is: L(t)=L y ·P w ·P d ·P h (t) (16); In formula (16), L(t) represents the power load at time t; L y Indicates the original load value; P w Indicates the ratio of weekly load to annual maximum load; P d It indicates the ratio of daily load to the maximum weekly load; P h (t) represents the ratio of hourly load to the maximum load of the day; 1) Based on historical load data, the ratio of weekly load to annual maximum load is calculated as: In formula (17), max(L y ) represents the annual peak load; L w (i) represents the load of the ith week; P w (i) represents the ratio of the load in the ith week to the annual peak value; 2) Calculate the ratio of daily load to the weekly maximum load based on the weekly load changes: In formula (18), max(L w (i)) represents the weekly peak load of the ith week; L d (j) represents the load on the jth day of the ith week; P d (j) represents the ratio of daily load to the peak load of the week in which it is located; 3) According to the change of daily load, calculate the ratio of hourly load to the maximum load of the day: In formula (19), L h (k) represents the load at the kth hour; max(L d (j)) represents the peak load on the jth day; P h (k) represents the ratio of the load in hour k to the peak load on that day.

7. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: In step 6, the relationship between each node is established according to the Bayesian network, and the weight of the upper node on a certain state is reversely inferred through the relationship between the upper node and the microgrid system to obtain the influence of the upper node on the system; According to the total probability formula, the Bayesian formula is: In the formula, P(B i |A) means that under the condition that event A occurs, event B i The probability of occurrence, also called the posterior probability; P(B i )>0, i=1,2,…,n represents the prior probability; P(A|B i ), i = 1, 2, 3, ..., n represents the posterior probability, n represents event B i The total number of categories, that is, all possible B i The number of events; Applied to the Bayesian network structure of the distribution system, event A can be represented as a load failure, event B can be represented as an upper-level component node, and the impact of the component state probability on the load node state, that is, the cause-effect relationship; Through this relationship, the known result state is reversed to the conditional probability caused by various causes. The larger the value, the greater the impact on the result, thereby identifying the weak part of the load impact; According to the node division in the distribution system, the corresponding set is reversed to the component node in turn, the component node state under the accumulation set is established, and the relationship between the corresponding system operation state and the component state is established.

8. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 5 is characterized by: The model of time-series source-load change is established as follows: Using the wind speed and light intensity forecast data provided by the ARIMA model in step 4, based on historical load data and predicted load demand, forecast the demand in different time periods to obtain the load time series model; combined with the DG output power forecast, obtain the model of time series source-load change, and analyze the matching between the output of power generation resources and the demand at different time points; the output power forecast model of the solar photovoltaic system in the DG output power forecast model is as follows: Where: P PV (t) is the predicted output power at time t; η is the efficiency of the photovoltaic panel; G(t) is the solar irradiance at time t; A is the total area of ​​the photovoltaic panel; A sh is the area covered by the shadow; R soiling is the percentage of performance reduction due to fouling; The output power prediction model of the solar wind system in the DG output power prediction model is as follows: Where: P wind (t) is the predicted output power at time t; ρ is the air density; A rotor is the swept area of ​​the wind turbine rotor; C p is the power coefficient, which depends on the design of the turbine and the wind speed; V(t) is the wind speed at time t; V ref is the reference wind speed.

9. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: In step 7, the distribution network system is initialized, and the logical relationship of each node of the Bayesian network is input; the sequential Monte Carlo sampling method is used to simulate the operation of the microgrid, and the model of the time series source-load change established by the time series model and the Bayesian network is used to sample the system components at time t. The state at a certain time is expressed as: In formula (21), t i The state set of the component nodes at the moment; x1, x2, …, x k Indicates that at time point t i The state of the system components is, where k is the index of the component; m is the number of original components; The node state of the sampled component is expressed as: k ∈[-1,0,1], x k =-1 indicates maintenance status, x k =0 means normal state, x k =1 indicates fault status; Calculate the reliability-related parameters of each node of the system components to determine the operating status of the microgrid, including: (1) According to t k The expression of load outage rate within a time period is: Where: In the time period t k The failure rate or outage rate of the i-th load point; t k The number of outages of load i in the M samplings at time; In the time period t k The failure rate or downtime of the i-th load point within m a 、m b are the number of type a and type b failures respectively; r a 、r b They are the repair time for category A and category B faults respectively; (2) The expression of the annual average reliability index of the load is obtained as follows: Annual load failure outage rate λ of microgrid distribution system (unit: times / year): In formula (24): Lpi represents the annual load failure outage rate of the i-th load point; represents the failure outage rate of the i-th load point in the k-th time period; Microgrid distribution system load annual fault outage duration u (unit: h / year): In formula (25), u Lpi represents the annual fault outage duration of the i-th load point, that is, the total number of hours that the load point is out of service due to faults in one year; It represents the average failure downtime of the i-th load point in the k-th time period; Average power outage duration r of microgrid distribution system load (unit: h / time): In formula (26), r Lpi It represents the average outage duration of the ith load point, that is, the average number of hours of outage for each fault.

10. The reliability assessment method of microgrid distribution system based on Bayesian network and sequential Monte Carlo method according to claim 1 is characterized by: The step 8 comprises the following steps: Step 8.1: First, establish the following distribution system reliability indicators: The average power outage frequency SAIFI is obtained from the number of power outages in the system in a year, unit: times / (household / year): In formula (27), N Lpi represents the number of users at the i-th load point; The average outage duration SAIDI is obtained from the system's power outage duration in one year, unit: h / (household·year): The average power outage duration CAIDI is obtained from the power outage duration of each household in one year, unit: h / (household·year): From this, we can get the average power supply availability rate ASAI and the average power supply unavailability rate ASUI indicators: Expected power shortage index EENS, unit: MW / (household / year): In the above formula, N Lpi 、M i are the number of users at load point i and the number of users forced to shut down; P Lpi (t k ) is the real-time load of load point i at time k; R is the set of load points in the distribution system; Step 8.2: Perform specific simulation analysis in the test system; The BUS6-F4 of RBTS is used as the test system; three independent microgrids including wind power, photovoltaic power, energy storage device, diesel generator and other micro power sources are added, the operating status of each node in the whole year is extracted, the load changes throughout the year are accumulated, and the power supply reliability index is obtained; Step 8.3: Analyze the simulation system through the model obtained in step 8.2, and calculate the improved results in combination with the reliability index in step 8.1.

Images