Probability risk assessment of island operation event for large scale photovoltaic plant

ABSTRACT Original methodology for the calculation of island operation event occurrence possibility is proposed for the photovoltaic power plants connected to medium voltage grid. Within the methodology description, event fault tree with number of pairs of basic initial events is discussed, where probability for the basic events are defined and flow chart for the proposed calculation is given. The novelty in the approach for probability calculation includes also proposal for development of frequency distribution functions for concerned power generation and consumption variables which have to be modelled for the probability calculations, for which the triangular distributions curves have been chosen. Being boundary conditions for definition of non-detected island event, borders of non-detection zones are defined for the common island detection protection functions: over frequency and under frequency, over voltage and under voltage and rate of change of frequency functions. Comparison of the island event probability is presented on case study example, with conclusion that standard voltage and frequency protection alone are not satisfying, but ROCOF protection could be suitable, subject to particular local grid characteristic and calculative verification of probability.


Introduction
Various islanding detection schemes have been proposed, such as method based on exponentially damped sinusoidal signal estimation [1], two level islanding detection method [2], deep learning hybrid method [3], and islanding detection based on probabilistic principal component analysis [4]. Many of them are based on passive methods for recognizing islanding condition. The main problem for passive methods is non-detection zone (NDZ) for small changes in power, voltage and frequency values and angles. Some authors analysed change of power during islanding [5,6] and general performance assessment of inverterbased photo-voltaic distributed systems [7] along with efficient energy management for PV systems [8] and maximum power point tracking in PV systems [9]. Regarding islanding operation analysis and recommendations for distribution networks, UK scenario is recommended for studying [10].
Common practice for network operators worldwide at the moment is not to allow the possibility of island operation of distributed generation plants which do not have regulation features although in recent times approaches for islanded operation have been proposed such as [11] for example, with a lot of further research in this direction [12,13]. But classical grid operator approach requires dedicated protection to detect the island operation as soon as possible and disconnect the plant [14,15]. Few main reasons for which such operation is unwanted are the following: • possible danger for maintenance and restoration personnel, when not aware that part of the network, although disconnected, is still powered, • voltage and frequency of islanded network part can slip out of defined range, power quality worsens, which adversely affect the loads and can result in damage of both customer and utility side equipment, • automatic re-closing of utility breaker may create a condition of asynchronous closure which can lead to serious damage of the equipment, • the very purpose of automatic re-closing of utility breaker, which is to clear the transient faults on the feeders, is not fulfilled, as fault arc will not extinguish if fed from distributed generation plant in island operation.
This paper deals with the assessing of probability risks for the occurrence of such island operation event in particular cases with following characteristics: • large scale photovoltaic plants, • connected to medium voltage grids, • with standard passive protection methods applied as island operation protection.
The original method for calculation of islanding event probability for the distributed generation plant is proposed by this article. The important factor in it is the characteristic power output profile of the plant. Previous research has been conducted in this direction [19], for output profile of inverter-based generation. Large multi-inverter photovoltaic plants are selected for their specific characteristics which influence stability of distributed active anti-islanding functions [20], but the method is applicable also to other types of generation plants. However, specifics of islanding detection with PV plants are considered [21,22].
The novelty of the approach proposed in this paper is in the terms of recognizing the importance and influence of consumer characteristics and particularly of the local electrical grid characteristics onto the probability of occurrence of the distributed source island event. It is proposed how to relatively simply assess, quantify and model this influence within the probability calculation.

Fault tree
The island operation of a power plant with the local grid part is an event which may occur during the plant operation. As such a situation is not allowed, this event can also be considered as fault. As a first step in determining the fault probability, a flowchart that leads to the occurrence of the fault is required. Such a diagram is called a fault tree. There are two basic (initial) events that need to occur simultaneously in order to cause the occurrence of the plant island operation.
One of the events is the state of balanced power production of the plant and power consumption in the local grid area, and the second event is the disconnection event of the switching device which will disconnect from the main grid exactly the part of the local grid with the balanced conditions.
At that point in time, in such local part of the grid, practically all locally generated energy is consumed also locally and there is no or almost no power flow through the switching device on which the disconnection occurs. Given that the protection methods have a certain non-detection zone, due to the protection settings, there may be a certain small exchange of energy with a main grid. This exchange must be such small that the disturbance caused by this excess or lack of energy in the properties (measured parameters) of just formed island system remains within the limits of the detection of the protection method, i.e. within its non-detection zone. If this is the case, an island operation event has occurred.
However, it is possible to have a large number of such pairs of initial events, depending on the configuration of the local grid, and the simultaneous occurrence of any two paired events will lead to the island operation event. The number of pairs depends directly on the number of switching devices that may cause the local grid to disconnect from the main grid. The fault tree is therefore not complicated and is shown in Figure 1.
In the first step, therefore, it is necessary to define all possible pairs of initial events in the observed part of the network to which the distributed generation plant is connected, i.e. to define how many possible "local grid subsystems" is there (indicated by n on Figure 1).

Number of local subsystems
As a first step, number of possible subsystems, local grids, on which island operation is possible, has to be determined. The criteria to determine this number (number n from Figure 1) is the following: P min _grid_1 < . . . < P min _grid_i < . . . < P min _grid_n < P max _plant < P min _grid_n+1 (1) where, P min_grid_1 is minimum load active power for the local grid behind first switching device, P min_grid_i is minimum load active power for the local grid behind switching device i, P min_grid_n is minimum load active power for the local grid behind switching device n, P max_plant is maximum generation active power for the connected power plant, P min_grid_n+1 is minimum load active power for the local grid behind switching device n + 1. The next step is to determine the probabilities of individual initial events. It is necessary to determine the probability of the event of balanced powers of generation and load in every possible part of local grid on which occurrence of island operation is possible. It is the probability of the first event in each of the defined pairs of events. Also it is necessary to define the probability of the switching-of event for every individual switching device in local grid on which the disconnection of the local grid from the main grid (electrical system) can be done. It is the probability of the second event in each of the defined pairs of events.

Probability of the balanced power condition
It can be assumed that the powers (both generation and load power, and also active and reactive power) are random variables. The values of these variables are in a certain range of possible values. The interval is from minimum to maximum power, the values of which are known or can be determined.
Therefore it is possible, on this defined interval, to represent the probability of power values with the random variable continuous probability density function.
Probability density function f of random variable x has following attributes: The probability for such random variable to have value within interval [a,b] is: As a next step, it can also be assumed that the generation power variables and load power variables are mutually independent variables. Without question, there can be some indirect dependency, but it would be impossible to quantify it (e.g. in case of solar power plants, when there is high sun the generation power is higher, but at the same time it is more hot and it affects the rise of the load power influenced by more air condition devices working). Therefore, it is disregarded, and for the calculations the assumption is that these are independent variables. Let us assume that variable x is the power plant generation power which can be in a range from x min to x max , and variable y is the load power which can be in a range from y min to y max , with respective probability density functions f X (x) and f Y (y).
For individual value of X from noted interval x min to x max , the probability for y to be in interval from Xx − to X+ x + is, according to (4): (5) where, P() is probability of given event; x − is the border of protection function non-detection zone for power of value X, in negative direction; x + is the border of protection function non-detection zone for power of value X, in positive direction.
The protection function non-detection zone is the zone in which the given protection function can not detect the island operation as the change of measured values is less then the protection setting threshold. It is equivalent definition to balanced conditions of generation and load.
The above probability given by (5) is related only to exact value of X, and is not of use for calculation of total probability, i.e. it is probability for individual fixed value of X that the value of Y is within the protection function non-detection zone, which does not mean a lot.
In order to calculate the probability of initial event, it is necessary to notice that the initial event of balance of generation and load power is actually a twodimensional variable, i.e. a random two-dimensional continuous vector composed of two random variables X and Y.
If X and Y are independent random variables, which is starting assumption, and if both have their probability density function f X (x) and f Y (y), then for the two-dimensional random vector composed of these two variables it is valid that the joint probability density function is [23]: where, f X (x) is probability density function of random variable X; f Y (y) is probability density function of random variable Y; f XY (x,y) is probability density function of two-dimensional random vector composed of two independent variables X and Y. The above formula is for the probability density function of the event. The probability itself can be calculated if the formula (4) is adapted: Now the known border conditions have to be used as an input to above formula: • X is variable of plant generation power, which can be in a range from x min to x max , • the range of non-protection zone for the protection function, for variable Y, is from X-x − to X+ x +.
When these border conditions are entered into (6) the probability of the first initial event (from the pair of initial events), i.e. the probability of the balanced generation and load power in defined part of local grid is: Where, P i1 is probability of the first initial event from the i pair of initial events, the probability of the balanced generation and load power in defined part i of local grid; x min is minimum possible generation power; x max is maximum possible generation power.

Occurrence of the disconnection event
The probability of the second initial event in every pair of events is the probability of the disconnection event of individual switching device on which local grid part together with the power plant can be disconnected from the rest of the grid. Statistical data about switching off of the switching devices can be obtained from available databases of the network operator or in some other way. It can be a case-by-case approach and assigning characteristics to a single switching device based on operator/dispatcher experience, or general approach based on general statistical data for the type of switching device and the number of its disconnections over a given period of time.
Although the statistical data for the switching device, and the derived characteristic (number of disconnections in a given time period) of an individual switching device is obtainable, the likelihood of the disconnection event (during that time period) is not easy to calculate, i.e. for that variable the expected value is known, but the probability density function is not known. It will be shown below that this will not be necessary either. Now when the characteristic for each observed switching device, the number of expected switching offs in a given time period, is known, this number for the switching device i can be indicated with p SUi , therefore: p SUi is number of expected switching off events of switching device i during observed period.
This number should include all disconnections of the switch, manipulation caused by intent, unintentional accidental switching off, switching off by protection, etc.
The probability of occurrences of island operation event has to be calculated for a certain period of time. For the most common period of observation, the period of one year will be taken, for many reasons most suitable.

Total probability
There are many pairs of initial events simultaneous occurrence of which leads to the occurrence of island operation event of the concerned power plant with part of local power system (shown in the fault tree on Figure 1). Initial events in pairs are independent, i.e. there is no dependence between the balance of powers in the local area of the grid and the disconnection of a particular switching device.
The probability of simultaneous events of a particular event pair (in the observed period) is therefore equal to the product of probability of one and probability of another event. The probability of this concurrent event will be marked with P i .
For the first event, the balance of generation and load power, the probability is calculated according to the formula (7). Since the probability of a second event is knot known, switching off of the switching device, but the expected number of disconnections (yearly) is known, the likelihood of the total event will be formulated in a slightly different way.
Let us take a look at the probability that there will not be a simultaneous occurrence of both initial events from one pair of events, in the observed period. This probability is equal to the probability that during the observed period (one year) neither one of the p SUi switchings of the switching device does not occur in the time when the balanced power conditions exist. The event "non existance of balanced power conditions" is complementary to the "event of balanced conditions" and therefore its probability is: Probability of the event that there will not be a simultaneous occurrence of both initial events from the pair i of events, not once in the observed period, is: where, P i ' is probability of the event that there will not be a simultaneous occurrence of both initial events from the pair i of events; P i is probability of the simultaneous occurrence of both initial events from the pair i of events; P i1 ' is probability of non-occurrence of the first initial event from the pair i of events; P i1 is probability of occurrence of the first initial event from the pair i of events, i.e. probability of balanced power conditions in the part i of local grid. The next step towards relation for total probability follows the identical logic. The probability that there will not be a simultaneous occurrence of both initial events from any of n pairs in observed period equals to the product of all probabilities calculated as per formula (9) for each of the pairs from 1 to n. Complementary to this probability is the value that is looked for, the probability for simultaneous occurrence of initial events in at least one of the event pairs. Respective relation is: where, P io is total probability of occurrence of power plant island operation event with part of local grid; n is number of local grid parts on which there exist a possibility for occurrence of island operation; P i1 , p SUi are as defined above.
It should be noted that the number n is determined by analysing the local grid configuration, statistical data p SUi is acquired from operational statistics, and each individual value of P i1 should be calculated using relation (7).
Description given in previous sections of this chapter can be formed into a process applicable for calculation of island operation probability for the power plant and can be presented in a form of flow chart. Such flow chart is shown on Figure 2

Load modelling
Usually in the statistics when the actual probability density function for the variable is not known, normal distribution is assumed. Even probability distribution of load power is often assumed to follow normal distribution [24,25]. But also, some studies have shown that generally the statistical distribution of the load power variables does not follow any of the common distribution functions used to mathematically describe the natural phenomena [26]. On the other hand, research on modelling focused on specific parts of the daily load diagram shows the possibility of applying normal and log-normal distribution for probabilistic modelling [27]. According to various researches described in [27], the load power by types of consumers/parts of the grid has a certain probability distribution and the probability of occurrence is certainly not the same for all power values within the range. To determine the characteristics of these distributions (mean and standard deviation in normal distribution) or the characteristics of some other distribution, a lot of data is needed. As outlined previously, what can be expected as easily available is the most commonly recorded medium, minimum and maximum values achieved. These three values are not much, but from them three characteristic points are obtained. In such cases, triangular distribution is most often used in probabilistic research [28]. Triangular distribution is the first approximation of all the functions stated above, but is simpler and easier for use in calculations. Probability density function (PDF) of the triangular distribution is described by relation: where a is lower limit of possible values of variable; b is upper limit of possible values of variable; c is expected value of variable (mode). In case that data is available only for the point on the beginning of the feeder, for subordinate parts of the local grid function it can be extrapolated. This is applicable to both active and reactive load power.

Generation modelling (MV connected PV plant)
The methodology described in chapter 3 is generic and can be applicable for any type of the power plant. In this case research is done only for the photovoltaic plants due to specific generation profile.
As the purpose of this assessment is to determine probability for the occurrence of subject islanding events in long term (installation life time) periods, suitable periods of observation are (repetitive) annual periods, and these are therefore of primary interest here. Also, in such way seasonal oscillations are eliminated by using the annual average daily diagram which gives representative enough result for this purpose. Of course, variations of generation output in different seasons of the year occur and the modelling procedure could instead be applied on e.g. monthly periods, but that would require much more input data (possibly not available) and effort for modelling without the according improvement in the result.
The average annual daily diagram of a PV plant shows generation from 6 am to 6 pm, and the rest of the daily diagram is the nighttime in which there is no generation. One half of the observed time (year base) is therefore a power plant without generation. In addition, it is known that the power ranges from minimum to maximum, whereby this data can easily be obtained, whether in the preparation and planning phase of the plant or by measurements during the exploitation of the power plant. Likewise, it is always possible to estimate the data of average annual generation, that is, of the average annual power of the power plant. This data (and the expected annual generation diagram) is one of the basic data when estimating the investment in power plant construction and is therefore always available from the planning phase.
References to very few studies conducted in the specific direction of the statistical (probabilistic) description of the photovoltaic generation power can be found in the available literature [29], but the possibility of modelling the quantities of sun-radiant radiation on the Earth's surface for different needs has been investigated by several researchers [30][31][32]. One research suggests that the distribution of the solar radiation probability density (directly related to photovoltaic power generation) is closest to the beta distribution [33]. For further research it would probably be possible to corelate such model of beta distribution function with input parameters that depend on other variables, e.g. geographical latitude or available data on the quantity of solar radiation at micro/macro locations. The approach in this work is to give one possibility for approximating the probability density function, and certainly this function can also be approximated in other ways, which does not significantly affect the methodology (calculating the probability of the power plant island operation) proposed herewith.
Also, in the case of power generation of photovoltaic plants, as well as in modelling the probability density of the load power presented in previous section, the goal is to easily model a function that will reflect the probability density of the output power variables, based on a limited input data. As first consideration in this research, normal distribution was used for modelling the probability density function [24], which is not exactly appropriate as the realistic distribution is not symmetrical.
If the average power of the power plant exceeds one third of the rated power, then there is possibility to present the curve by a classical triangular distribution, where from the average power data the expected power value c is calculated. Afterwards, it is proceeded in a similar way as in the modelling of the load probability distribution to get a function that describes the density of the distribution. The triangular distribution can be applied, according to the criterion of the average value of the variable, if the average value ranges between one and two thirds of the minimum and maximum values difference. With photovoltaic power plants, the average power in the daytime period does not exceed two thirds of rated value nowhere in the world [34]. As an example, on the upper diagram in Figure 3 it can be seen that the triangular distribution could roughly approximate the probability density distribution at Abuja, Nigeria, for which the average power is certainly greater than onethird of the rated power [33]. In the Figure 3 triangular distribution is compared to other suggested functions.
If the average power of the power plant is less than one third of the rated power, then even with the marginal triangular distribution (linear) there is error. In such cases, the distribution is more inversely proportional, i.e. the probability of the expected minimum power value is greater than for the marginal triangular distribution, while the probability of higher power is smaller. A triangular distribution can be applied here, but proportionally divided by the ratio of the average power to one third of the rated power.
The rest of the distribution is assigned to the minimum value, zero, i.e. to the time without power generation. This gives fairly defined expected mean value of the variable in this case. As an example, on the lower diagram in Figure 3 it can be seen that the triangular distribution could roughly approximate the probability density distribution in Belfast, UK, for which the average power is less than one third of the rated power.
In any case, the probability density function of the power is modified triangular distribution, as the outcome of the integration with formula (3), is one half. Other half is assigned to value zero (night mode) and does not influence the calculation.
It is important to note that the minimum possible reactive power (maximum inductive power of the power plant) can occur only at the same time as the maximum active power of the power plant. Also, the maximum reactive power (maximum capacitive power of the power plant) can only occur at the same time as the minimum active power plant power. However, for the maximum and minimum active power of the power plant there is a full range of potential reactive power values, i.e. it is not fixed for given active power.
As for reactive power, the photovoltaic plants of large size which are connected to MV grid have reactive power component even if inverter power factor is set to one [35]. Reactive power behaviour itself is not directly correlated with the active power. What is intuitive is, however, that the night mode of operation (with maximum capacitance) prevails because it consumes at least 50% of the time a year. Thus, it is known that the expected value of the reactive power variables is the maximum capacitive power. The probability distribution of other variable values is not intuitive, and it is only possible to obtain them from the long-term measurements of the reactive power of the plant. However, in the absence of such data, in this case, also a triangular (boundary triangular, linear) distribution for probability density function can be used, which, as stated, is predominantly used in situations of insufficient knowledge of the variable behaviour.

Selected protection functions and non-detection zones
The non-detection zone has a major impact on the occurrence on island operation. Depending on the protection function, non-detection zones differ. For the selected protection functions here are defined the borders of non-detection zones, expressed as the unbalance of the active ( P) and/or reactive ( Q) power, i.e. the range of P and Q in which the island operation will not be recognized by the applied protection functions.
Among the countries that have the longest experience of connecting small power plants to the grid is certainly the UK, whose developed engineering practices serve as a guide to many countries in developing their own rules [36]. According to their recommendations (colloquially known as G59/3) for generation connection to the distribution grid, the recommended protection functions that network operators require are as follows [37]: Recommended LOM function is rate of change of frequency, ROCOF, especially for smaller size plants. All above protections are passive type protections.
Active protections are only considered as additional option.
For under frequency and over frequency protection function boundary values are: (13) where, f max is over frequency protection setting value; f min is under frequency protection setting value; Q − is boundary value of non-detection zone for value of Q load smaller then Q DG ; Q + is boundary value of non-detection zone for value of Q load larger then Q DG ; Q f is local grid quality factor. For under voltage and overvoltage protection function boundary values are: where, V max is over voltage protection setting value; V min is under voltage protection setting value; P − is boundary value of non-detection zone for value of P load smaller P DG ; P + is boundary value of non-detection zone for value of P load larger then P DG .
As it is necessary for islanding protection to be unsensitive for standard disturbances in the grid, such as voltage dips, overvoltages, harmonic distorsions and frequency sliding, protection setting values for voltage and frequency need to allow for certain range of deviation. This is the reason itself for the existence of non-detection zones. It shall be noted that maximum range setting values are used for calculations presented in this paper and these are usually defined by grid normative documents. Here reference is made to EN 50160, according to which these values are: f max = 51 Hz, f min = 49 Hz, V max = 1.1p.u., V min = 0.9p.u.
For rate of change of frequency function boundary values are: where, (df /dt) set is ROCOF protection setting value; SH = i S i H i is total kinetic energy of rotating machines connected to local grid part disconnected from main grid (where S is power base and H respective inertia constant).   Notes: axis scales differ for each graph, load PDFs are relative to load size.

Case study example
The principle diagram on Figure 4 shows the operational connection configuration of the grid feeder to which case study plant PVP Kanfanar (in Croatia) is connected [38]. The diagram also shows the possible points of disconnection between the main point of supply, substation HV/MV and the power plant, where the power grid and power plant can separate and thus where the creation of island operated local grid parts is possible. Also shown on the diagram are the rated power values of the transformers for all MV/LV distribution substations connected to the feeders. Figure 5 shows the indication of local grid parts on which island operation is possible. There are three possible local grid parts. Figure 6 shows the representation of probability density functions for the power variables after the modelling, done in accordance with proposed method.

Additional data
Additional data needed and used for calculation is: • The frequency of disconnection events for switching devices in observed local grid for which calculations of island probability are being done, • Quality factor for local grid parts, • Inertia constant for rotating machines in local grid (substitute inertia constant for whole grid part) which is needed for island operation probability calculation in case of applied ROCOF protection df /dt.
In this case the data on disconnection events has been acquired from plant SCADA system data which number of events [39] has been split among the switching devices on the feeder with the majority assigned to circuit breaker on the feeder beginning. The assumption is that feeder breaker is to switch off the feeder much more often then the disconnecting switches on the feeder.
Required power factor for the consumers connected to the grid is from cosϕ = 0.95 inductive to cosϕ = 1, and due to the compensation applied and contribution of cabled grid, expected overall power factor in MV grid is larger than cosϕ = 0.98. The corresponding realistic grid quality factor is Q f = 0.45 which value was taken in calculations.
Regarding the inertia, hereto are noted assumptions and parameters considered in calculations for df /dt NDZ: • Inertia constant for inverters equals to zero, • Base power for individual local grid part equals to installed load power of that grid part, • The protection setting was set at value 0.3 Hz/s and it was assumed that the portion of rotating machine loads in any part of the local grid equals to one third of installed load power, and the assumed inertia constant of such small rotating machines was 0.3 s [40], whilst inertia constant of all other consumers equals to zero.

Results
Using the proposed modelling for the power probability density functions and the proposed methodology for calculation of island operation probability risk, values for the probability for occurrence of island operation event during one year period are calculated, for the case study example. Moreover, the calculation has been done also for whole range of photovoltaic plant sizes theoretically connected to the same point, to compare the eventual differences. Table 1 presents the probability value results. It can be seen that in cases when frequency protections (f < > ) and voltage protections (U < > ) are applied, the island operation probability is extremely high. On the other hand, in case that rate of change of frequency (df /dt) protection function is applied, island operation probabilities fall significantly.  For comparison, further calculations were conducted, one on the model where operating connection of plant was on another feeder (1 and 2 as per Table 2) of the same grid (not presented on figures given in this paper), and also for the classic plant with synchronous generator (SG) instead of multi-inverter PVP, in same connection point. Value used for inertia constant for such machines is as per literature [41]. Results are presented in Table 2.
It can be subject to further discussion, but according to [42] the probabilities of 1-10% for a failure per year are acceptable/tolerable in case of low and moderate consequences of such failure, but not acceptable in case of severe consequences. Standard voltage and frequency protection alone are not satisfying with inverter generation, but ROCOF protection could be suitable, subject to particular local grid characteristic and calculative verification of probability. On the other hand, with synchronous machines, ROCOF in some cases does not provide satisfying level of protection.

Conclusion
The proposed method for assessing the probability risk of island operation for distributed sources, based on island operation fault tree and local grid configuration data, has been applied for the calculation of probability of such operation in case study of large photovoltaic plants connected to medium voltage grid. Probability density function proposed for both load and generation power modelling for active and reactive power is triangular distribution, as input data available is often limited. The slight modification of triangular distribution is proposed in accordance with the expected sun radiation at the location of the plant. Definition of the non-detection zone boundaries for recommended island operation protection functions is presented. Island operation probability calculation results prove that, among the compared protection functions, the most effective protection function in such case is definitely rate of change of frequency.

Disclosure statement
No potential conflict of interest was reported by the authors.