Extended Flow-Based Security Assessment for Real-Sized Transmission Network Planning

The evolution of electric power systems involves several aspects, dealing with policy and economics as well as security issues. Moreover, due to the high variability of operating conditions, evolution scenarios have to be carefully defined. The aim of this paper is to propose a flow-based procedure for the preliminary security analysis of yearly network evolution scenarios at the real scale level. This procedure is based on hourly load and generation conditions given by market solutions, and exploits Power Transfer Distribution Factors and Line Outage Distribution Factors to determine N and N−1 conditions, properly accounting for possible islanding in the latter case. The analysis of overloads is carried out by dealing with big data analysis through statistic indicators, based on power system background, to draw out critical operating conditions and outages. The procedure is applied to a provisional model of a European high voltage network.


Introduction
The provisional evolution process of power transmission network is driven by different factors. The studies are based on the definition of generation and demand scenarios [1][2][3], and of transmission expansion plans, involving different kinds of decisions and approaches [4,5]. Possible strategies for generator operation should account for present and future market platforms [6][7][8]. Moreover, the impact of operation patterns on the network has to be assessed, usually through steady-state network analysis accounting for a range of different load/generation outlines [9,10]. Network operation security should be evaluated, highlighting possible violations of operating ranges (voltage values, power ratings) with the whole network in service (N condition) as well as in the presence of single independent outages of network elements (N−1 condition) [11][12][13].
The handling of a huge number of operating conditions and a remarkable number of outages could imply a long simulation time to draw exact results. Sensitivity factors can be employed to determine quick and affordable indications of line flows due to changes in generation/load patterns and in network configurations [14]. In particular, power transfer distribution factors (PTDFs) represent the sensitivity of a line flow to a variation of power injection at a couple of nodes, whereas line outage distribution factors (LODFs) represent the variation of a line flow due to the outage of another line. Methods for the calculation of PTDF and LODF matrices from network data, even with generalized matrix formulation in the presence of islanding, are presented in [15,16], whereas a technique is synthesized in [17] to reduce the computational burden of the matrix formulation, according to topological analysis. - The treatment of big data in realistic transmission network planning is faced by means of threshold-based sampling and parameter correlation. -Proper indices and metrics are defined, to individuate in a readable and synthetic way the relation of N and N−1 overloads with specific factors, and to individuate critical cases and conditions. The proposed methodology is described in Section 2. The case study and test results are reported in Section 3. Section 4 concludes the paper.

Determination of Distribution Factors
Network sensitivity factors are nonlinear, due to the interdependence of active and reactive power levels on nodal voltage amplitude and phase angles. However, in high voltage networks with normal loading, linear DC power flow equations can be exploited involving only phase angles, and sensitivity factors can be obtained by matrix formulation [14].
In a transmission network composed of N i nodes and N b branches, the N b × 1 vector of active power flow on branches F in a defined operating condition, represented by time interval t, can be determined as follows: where P(t) is the vector of N i − 1 nodal power injections in time interval t, excepting the slack bus. Moreover, H represents the N b × (N i − 1) branch-node matrix of PTDFs, and it is determined from B dc , the (N i − 1) × (N i − 1) reduced nodal admittance matrix, and from B b , the N b × (N i − 1) branch-node admittance matrix.
In the case of the outage of the branch k, from node i to node j, with initial flow f (k, t), the postoutage flow on the branch bf (b, k, t) can be determined through PTDFs [14]: where f (b, t) is the pre-outage flow in branch b, h(b, i) and h(k, i) are the elements of PTDF matrix for branches b and k, respectively, at node i, whereas l(b, k) is the LODF giving the flow variation in branch b after the outage of branch k.

Analysis of Islanding Conditions
The matrix calculation of sensitivity factors can be applied as far as B dc is invertible to obtain the PTDFs from Equation (1), and each PTDF value is not equal to 1, to obtain LODFs from Equation (2). These conditions are not satisfied in the presence of network islands in N and N−1 conditions, respectively. In order to deal with the N−1 islanding cases, a specific method is carried out. If electric islands are present in the N condition (e.g., non-synchronous areas), a specific determination of PTDFs of the isolated portion of the network is required, by choosing a further slack bus within it.
The set Ω is of branches whose outage causes an islanding is determined by carrying out a network connection analysis based on a branch-node incidence matrix, and for each branch k in the set Ω is , the same analysis is exploited to determine the set of nodes in the island Ω k [37,38].
In order to determine the post-outage steady-state flow in branch b, within the non-islanded part of the network Ω p , a set of methods can be considered, involving the contribution of primary frequency response [39,40] or individuating the necessary increasing/decreasing generation reserve [41,42]. These methods necessarily involve the knowledge of system operation strategies and of generation plant characteristics, which are usually not known in the planning stage. Hence, the adopted approach aims at simulating the outage, i.e., nullifying the power flow on the branch k, by subtracting the contribution of the net power injection at the nodes of the island Ω k . This implies the exploitation of a slack generator, usually giving low contribution in order to have room for compensating losses in AC formulation. Therefore, the power flow in branch b after branch k outage is determined by means of PTDFs: Hence, for post-outage analysis, a modified LODF matrix L and an extended PTDF matrix H (dimensions N b × N b × (N i − 1)) are defined. Their elements are determined as: With these assumptions, and defining O as the N b × 1 vector of ones, the N b × N b matrixF of post-outage flows in all branches after the disconnection of any other branch in a single operating condition is formulated as: where diag(F(t)) has N b × N b dimensions.

Security Analysis in Varying Operating Conditions
For long-term development planning scopes, a set of different conditions should be studied in order to evaluate system performances with varying inputs, such as load demand and generator availability, analyzing yearly horizon instead of a single point in time [24,43]. Even if some studies, aimed at overloading estimation, define loading samples to reduce the computational burden [44], in this paper the yearly horizon is divided into N t consecutive time intervals to take into account the dependence on the period. Moreover, the perspective of this paper is to study network behavior under a defined network configuration resulting from an evolution plan which is valid for the full yearly horizon under consideration; therefore, the sensitivity factor matrices do not vary for each time interval. The outcomes of preliminary security analysis would be useful to develop specific network management strategies for operational planning, varying over the time horizon and involving further technical details, to analyze the risk of failure and to reduce the detected overloads.
From this perspective, let P and F be the (N i − 1) × N t matrix of yearly nodal power injections and the N b × N t matrix of yearly branch flows, respectively. The PTDF matrix can be used in a single stage to relate the two quantities: Moreover, by extending Equation (6) in the dimension of time intervals, the post-outage flow matrixF, having N b × N b × N t dimensions, is determined as follows: where the elements of the N h × N b × N h matrix Γ and of the N b × N b × N h matrix Φ are defined as follows (τ is a further time index for the sake of notation): Energies 2020, 13, 3363 From pre-and post-contingency flows, the loading indices λ(b, h) andλ(b, k, h) for the b-th branch in each time interval t, in N and N−1 conditions, respectively, are obtained by the following matrix formulation: where the symbol | | stands for the matrix where each element is the absolute value of the corresponding element of the matrix in the argument, and the N b × 1 vector M includes the reciprocal of the steady-state rating of the N b branches r(b).
Once the N b × N h matrix Λ and the N b × N b × N h matrixΛ are evaluated, proper big data management techniques are employed to analyze network behavior and to point out specific features, with particular care given to overloads. Therefore, data processing based on threshold-based sampling, averaging, peaking and parameter correlation are exploited according to power system operation. In N conditions, the following indices can be determined; a relevant, detailed formulation is provided in Appendix A for sake of readability: The global network loading analysis is carried out through two hourly indices u R (t) and u + R (t) relating power flows and overloads to line rating r(b), respectively, whose determination is detailed in Equations (A3) and (A4) of Appendix A.
Moreover, network analysis can aim at highlighting different behaviors according to voltage levels, since the overload of a branch at a higher nominal voltage, linked to higher transfer capacities, can imply more concerns with respect to shorter branches at lower levels. Therefore, defining for each voltage range v the relevant set Ω(v) of branches, the number of overloaded lines a(v) can be obtained, along with average number of hours in overload n(v) and average overload amount λ(v), defined as in Equations (A5) and (A6) of Appendix A.
On the other hand, the treatment ofΛ in order to highlight different implications of the study is synthetically depicted in Figure 1.
where the symbol | | stands for the matrix where each element is the absolute value of the corresponding element of the matrix in the argument, and the × 1 vector includes the reciprocal of the steady-state rating of the branches . Once the × matrix and the × × matrix are evaluated, proper big data management techniques are employed to analyze network behavior and to point out specific features, with particular care given to overloads. Therefore, data processing based on threshold-based sampling, averaging, peaking and parameter correlation are exploited according to power system operation. In N conditions, the following indices can be determined; a relevant, detailed formulation is provided in Appendix A for sake of readability: The global network loading analysis is carried out through two hourly indices and relating power flows and overloads to line rating , respectively, whose determination is detailed in Equations (A3) and (A4) of Appendix A.
Moreover, network analysis can aim at highlighting different behaviors according to voltage levels, since the overload of a branch at a higher nominal voltage, linked to higher transfer capacities, can imply more concerns with respect to shorter branches at lower levels. Therefore, defining for each voltage range the relevant set Ω of branches, the number of overloaded lines can be obtained, along with average number of hours in overload and average overload amount , defined as in Equations (A5) and (A6) of Appendix A.
On the other hand, the treatment of in order to highlight different implications of the study is synthetically depicted in Figure 1. The analysis can be focused at the influence of outages, as described by the green arrow in the left part of Figure 1, calculating the cumulative indices on the green plane:  The analysis can be focused at the influence of outages, as described by the green arrow in the left part of Figure 1, calculating the cumulative indices on the green plane: The number of time intervals where the branch b is interested by at least an N−1 overloadn j (b) can be deduced fromn s (b, t), whereas the average N−1 overload of branch b is obtained as detailed in Similarly to N conditions, further indices can be evaluated: For the v-th voltage range, once the average numberâ(v) of intervals with an N−1 overload is determined, indices representing the number of N−1 overloaded branchesn(v) and average N−1 overloadλ(v) can be defined, as detailed in Equations (A15) and (A16) in Appendix A.
The analysis of lines in service and of outage implications can also be carried out for a single operating condition, with a bi-dimensional analysis, following the blue and green directions depicted in Figure 1 over a slice of the matrixΛ representing a single time interval t (analogous to the orange plane in Figure 1). Thus, the following indices can be obtained: - The number of overloads after branch k outage,n c (k, t); - The average overload due to the outage of branch k,λ c (k, t); - The number of N−1 overloads of branch b,n s (b, t); - The average N−1 overload of branch b,λ s (b, t).
In this way, a generalized procedure for network analysis can be carried out by comparing the obtained indices with standard thresholds. Further case-specific classes and indices can be defined according to the network peculiarities.

Test System and Results
The analysis is carried out on a provisional model of ENTSO-E country perimeter at year 2030, coupled with a provisional model of the transmission network focused on Italian border regions. In particular, generation spreading per technology and per country is taken from evolution scenarios, whereas network development is determined according to cross border interest projects and internal projects with a forecast entry in service before 2030 [45].
A synthetic scheme of the test system is reported in Figure 2, involving 6 countries and 11 market zones, with 16 internal exchanges (light blue) and 24 exchanges with external zones (dark blue). The network includes 9045 nodes, 8817 lines-detailed in ranges of voltage level in Table 1   The method is implemented in a MSExcel ® -MatLab ® environment, exploiting Matpower [46] for network analysis. The operating condition of the test system is determined according to a nodal shifting of electricity market results, in order to obtain power production levels for each generator according to technology and zone, and to obtain load level at each node, finally determining nodal power injections [47]. In particular, HVDC links are represented by coupled active power injectionwithdrawal at extreme nodes, whereas phase shift transformers (PSTs) are modelled by equivalent power injections at extreme nodes, depending on shift angle [48].
The whole procedure, running on a 64-bit workstation with 3.50 GHz processor and 16 GB RAM, and exploiting virtual parallel calculus, takes roughly 50 h from yearly market input data to outputs of N and N−1 network analysis. To carry out the preliminary security analysis, more than 1.5·10 12 cases are analyzed (13,123 × 13,123 × 8760). In the following, results refer to line loading, neglecting the representation of transformer loading levels for sake of readability; therefore, 681·10 9 cases (8817 × 8817 × 8760) are handled.  The method is implemented in a MSExcel ® -MatLab ® environment, exploiting Matpower [46] for network analysis. The operating condition of the test system is determined according to a nodal shifting of electricity market results, in order to obtain power production levels for each generator according to technology and zone, and to obtain load level at each node, finally determining nodal power injections [47]. In particular, HVDC links are represented by coupled active power injection-withdrawal at extreme nodes, whereas phase shift transformers (PSTs) are modelled by equivalent power injections at extreme nodes, depending on shift angle [48].

Yearly Operation-N Conditions
The whole procedure, running on a 64-bit workstation with 3.50 GHz processor and 16 GB RAM, and exploiting virtual parallel calculus, takes roughly 50 h from yearly market input data to outputs of N and N−1 network analysis. To carry out the preliminary security analysis, more than 1.5 × 10 12 cases are analyzed (13,123 × 13,123 × 8760). In the following, results refer to line loading, neglecting the representation of transformer loading levels for sake of readability; therefore, 681 × 10 9 cases (8817 × 8817 × 8760) are handled.

Yearly Operation-N Conditions
The graphical representation of the N b × N t matrix Λ of yearly loading conditions in the absence of contingencies, obtained by means of PDTFs, and with 77.2 × 10 6 cases (8817 lines × 8760 h), is reported in Figure 3, where overload cases are represented by markers whose different colors represent overload intensity levels. It can be deduced that the network is interested by slight N overloads affecting a limited set of lines, although no overload is detected in 68 h only.
Energies 2020, 13, x FOR PEER REVIEW 8 of 19 The graphical representation of the × matrix of yearly loading conditions in the absence of contingencies, obtained by means of PDTFs, and with 77.2·10 6 cases (8817 lines × 8760 h), is reported in Figure 3, where overload cases are represented by markers whose different colors represent overload intensity levels. It can be deduced that the network is interested by slight N overloads affecting a limited set of lines, although no overload is detected in 68 h only.
A synthesis of performance indices is reported in Table 2. It can be seen that the total number of overloads corresponds to the 0.25% of the total combinations line-hour, affecting only 1.5% of lines. A particular focus should be given to extreme conditions, registered at hours 849-852 for maximum and for lines 2992-2993, pertaining to Voltage Range C, for maximum . Maximum is observed at hour 5773 for 23 lines.  Values of the indices by voltage range are reported in Table 3, where the number of branches in overload for less than 15% of the year and the number of branches with average overload higher than 1.2 p.u.
are reported as well. It can be observed that overloaded lines in Range A are 2.2% of the total number, but they show sporadic overloads, since only 5 of them are in overload for more than 15% of the year. Lines in Range A have low overloads as well; considering a 1.2 p.u. overload as acceptable in planning standards [49,50], only 4 lines are troubling. On the other hand, the 9 focused lines in Range C are subject to overloads for more than 40% of hours on average, and with a mean loading value higher than 2. Lines in Range B show an intermediate behavior in each aspect. A synthesis of performance indices is reported in Table 2. It can be seen that the total number of overloads corresponds to the 0.25% of the total combinations line-hour, affecting only 1.5% of lines. A particular focus should be given to extreme conditions, registered at hours 849-852 for maximum n q (t) and for lines 2992-2993, pertaining to Voltage Range C, for maximum n t (b). Maximum λ q (t) is observed at hour 5773 for 23 lines. Values of the indices by voltage range are reported in Table 3, where the number of branches in overload for less than 15% of the year a 1 (v) and the number of branches with average overload higher than 1.2 p.u. a 2 (v) are reported as well. It can be observed that overloaded lines in Range A are 2.2% of the total number, but they show sporadic overloads, since only 5 of them are in overload for more than 15% of the year. Lines in Range A have low overloads as well; considering a 1.2 p.u. overload as Energies 2020, 13, 3363 9 of 19 acceptable in planning standards [49,50], only 4 lines are troubling. On the other hand, the 9 focused lines in Range C are subject to overloads for more than 40% of hours on average, and with a mean loading value higher than 2. Lines in Range B show an intermediate behavior in each aspect. For global network analysis, given that the sum of the ratings of all lines b r(b) is equal to 26.2 GW, trends of indices are reported in Figure 4. It can be noted that u R (t) ranges between 3.1% and 9.6% (maximum at hour 849), with an average of 5.4%, resulting in a weakly charged network. This conforms to the assumptions of the method, based on DC power flow equations. On the other hand, u + R (t) reaches a maximum of 0.097% at hour 848, although the mean value is 0.0190%, the median value is 0.0133% and the 90-percentile is 0.0491%, highlighting good performances.  For global network analysis, given that the sum of the ratings of all lines ∑ is equal to 26.2 GW, trends of indices are reported in Figure 4. It can be noted that ranges between 3.1% and 9.6% (maximum at hour 849), with an average of 5.4%, resulting in a weakly charged network. This conforms to the assumptions of the method, based on DC power flow equations. On the other hand, reaches a maximum of 0.097% at hour 848, although the mean value is 0.0190%, the median value is 0.0133% and the 90-percentile is 0.0491%, highlighting good performances.

Yearly Operation-N−1 Conditions
The N−1 analysis is therefore carried out by means of the LODF matrix, and the × × loading matrix is obtained. It is worthwhile to mention that 1583 cases of N−1 islanding are detected, i.e., 18% of total outages, and they are efficiently dealt with through the developed procedure.
According to the definition of indices in Section 2.3, for the sake of representation, the values of , ℎ and , ℎ are depicted in Figures 5 and 6, respectively, representing the values assumed by the indices in 77.2·10 6 cases. It can be derived that at least one N−1 overload in each hour is registered for 9 lines (rows full of markers in Figures 5 and 6).

Yearly Operation-N−1 Conditions
The N−1 analysis is therefore carried out by means of the LODF matrix, and the N b × N b × N t loading matrixΛ is obtained. It is worthwhile to mention that 1583 cases of N−1 islanding are detected, i.e., 18% of total outages, and they are efficiently dealt with through the developed procedure.
According to the definition of indices in Section 2.3, for the sake of representation, the values of n s (b, h) andλ s (b, h) are depicted in Figures 5 and 6, respectively, representing the values assumed by the indices in 77.2 × 10 6 cases. It can be derived that at least one N−1 overload in each hour is registered for 9 lines (rows full of markers in Figures 5 and 6). Energies 2020, 13, x FOR PEER REVIEW 10 of 19  The summary of indices is reported in Table 4. It can be noted that N−1 overloads are seen in 0.8% of branch-hour combinations (the total number of markers reported in Figures 5 and 6), and affect 3.8% of branches. Moreover, the cases with overload for the maximum number of outagesi.e., when , reaches the maximum value of 8817 (red markers in Figure 5)-are detected in 5.9% of total overload cases, affecting 12 lines, mostly in Voltage Range B. The maximum level of is registered at hour 851, and the maximum level of is observed for lines 2992-2993. Regarding overload intensity, , ℎ ≥ 2 is registered for 7.3% of total overloads, affecting 31 lines, whereas the maximum levels of have analogous features of maximum .  The summary of indices is reported in Table 4. It can be noted that N−1 overloads are seen in 0.8% of branch-hour combinations (the total number of markers reported in Figures 5 and 6), and affect 3.8% of branches. Moreover, the cases with overload for the maximum number of outagesi.e., when , reaches the maximum value of 8817 (red markers in Figure 5)-are detected in 5.9% of total overload cases, affecting 12 lines, mostly in Voltage Range B. The maximum level of is registered at hour 851, and the maximum level of is observed for lines 2992-2993. Regarding overload intensity, , ℎ ≥ 2 is registered for 7.3% of total overloads, affecting 31 lines, whereas the maximum levels of have analogous features of maximum . The summary of indices is reported in Table 4. It can be noted that N−1 overloads are seen in 0.8% of branch-hour combinations (the total number of markers reported in Figures 5 and 6), and affect 3.8% of branches. Moreover, the cases with overload for the maximum number of outages-i.e., when n s (b, t) reaches the maximum value of 8817 (red markers in Figure 5)-are detected in 5.9% of total overload cases, affecting 12 lines, mostly in Voltage Range B. The maximum level ofn q (t) is registered at hour 851, and the maximum level ofλ t (b) is observed for lines 2992-2993. Regarding overload intensity,λ s (b, h) ≥ 2 is registered for 7.3% of total overloads, affecting 31 lines, whereas the maximum levels ofλ j (b) have analogous features of maximum λ j (b). Regarding time distribution, from Figure 6 it stems that, for lines with identifier 2500-3000 and 6500-7000, in Ranges B and C, higher overloads are registered in December, January and February, whereas the critical period for lines with identifier 7500-7800, in Range B, is between May and June.
In Table 5, the analysis of the 340 lines in N−1 overload by voltage range is reported. Further to indices analogous to the ones proposed in the N condition (see Table 3), the average number of N−1 events causing and overloadê 1 (v) is estimated, whereasê 2 (v) represents the number of lines of the vth voltage range with averagen s (b, t) lower than 10% of total N−1 events. Table 5. Characterization of overloaded lines in the N−1 condition. It can be seen that the increase of overloaded lines is most evident in Range A, where more than 6% of lines are affected. However, it is confirmed that higher voltage levels experience less frequent and less intense overloads, seldom reaching remarkable peaks, as compared to Range C.

Voltage Range vâ(v)n(v)λ(v)â 1 (v)â 2 (v)ê 1 (v)ê 2 (v)
For N−1 global network indices, no remarkable difference is seen between u R (t) andû R (t). Moreover,û + R (t) has an average of 0.106%, and a maximum of 0.338% at hour 849, its median is 0.0849%, and its 90-percentile is 0.2165%. These values can be acceptable in N−1 planning, considering the probability of events [51]. Distributions of u + R (t) andû + R (t), with respect to total network demand (between 107.7 MW at hour 8107 and 272.4 MW at hour 3846), are shown in Figure 7. As expected, it can be seen that N−1 events make overloads increase, although severe conditions are not linked to maximum demand (values higher than the 90-percentile are also seen for half of the peak load). This can be ascribed to the load distribution among zones. In Table 6, the analysis of indices by outage-hour, , and , , is synthesized.
Overloads  In Table 6, the analysis of indices by outage-hour,n c (k, t) andλ c (k, t), is synthesized. Overloads are seen in the majority of outage cases. The maximum number of N−1 overloadsn c (k, t) is 84, observed at hour 849 for the outage of line 1036, in Voltage Range A, whereas the maximumλ c (k, t) of 2.43 is seen for the outage of line 2924 at hour 2842, and values ofλ c (k, t) ≥ 2 are reached in 2.9% of total overloads. Table 6. Yearly network performances in the N−1 condition (lines in outage-hours).

Index Value
Total overload occurrences 76.6 × 10 6 Maximumn c (k, t) 84 Maximumλ c (k, t) 2.43 Number of cases withλ c (k, t) ≥ 2 2.25 × 10 6 Hourly analysis by means of indicesn h (b, k) andλ h (b, k) is carried out on 77.7 × 10 6 cases (8817 lines × 8817 lines), and reported in Table 7. It can be deduced that N−1 overloads are observed in 1.5% of total cases. The 132 lines in overload in the N condition experience N−1 overload for more than 8600 outages. Moreover, N−1 overload is present in 191 lines for less than 10 outages, in 11 lines for a range of 10 to 100 events and in 6 lines for 100 to 600 outages. As regards the values of indices, 31 lines in service haven h (b, k) > 3000 for almost all outages, whereas other 68 lines are interested byn h (b, k) > 3000 in less than 20 N−1 events. Only in 13 cases, each affecting a single line in service, isn h (b, k) > 8000 is observed. Values ofλ h (b, k) ≥ 2 are present in 3% of overloads, affecting 4 branches in any N−1 event, and another 16 branches in less than 5 outages, but only in 7 cases is â λ h (b, k) ≥ 3 observed. Table 7. Yearly network performances in N−1 condition (lines in service-lines in outage).

Index Value
Total overload occurrences 1.165 × 10 6 Total number of lines affected by overload 340 Number of cases withn h (b, k) ≥ 8000 13 Number of cases withλ h (b, k) ≥ 2 35.3 × 10 3 In order to individuate the most troubling conditions, for each couple line in service-line in outage (i.e., b, k), the hour with the maximum λ + (b, k, t), defined as critical overload, is found. The occurrences of critical overloads in each hourn cr (t) are reported in Figure 8, where the bullet color represents the number of lines affected by critical overloads in that hourâ cr (t). Critical overloads are detected in 144 h (total number of markers in Figure 8), mostly in winter. The highest number is 124.4 × 10 3 for hour 848 (10.7% of the total), whereas the maximum number of 48 lines with critical overload is registered in hour 849. This last hour is the most critical condition for different aspects, deserving a detailed analysis in the next subsection.
occurrences of critical overloads in each hour are reported in Figure 8, where the bullet color represents the number of lines affected by critical overloads in that hour . Critical overloads are detected in 144 h (total number of markers in Figure 8), mostly in winter. The highest number is 124.4·10 3 for hour 848 (10.7% of the total), whereas the maximum number of 48 lines with critical overload is registered in hour 849. This last hour is the most critical condition for different aspects, deserving a detailed analysis in the next subsection.

Analysis of Specific Loading Condition
The line loading levels in the N condition at hour 849 show that 74 lines are in overload in the N condition, and 22 of them, in Ranges B and C, have , 849 > 2. The global network performance is described by 849 = 9.63% and 849 = 0.095%.
As regards N−1 analysis, involving 77.7·10 6 cases, synthetic indices are reported in Table 8. Overloads are observed in 0.8% of combinations, affecting 193 lines. Peak overloads, with , ℎ, 849 ≥ 5, are detected only in 6 lines due to outages close to lines already overloaded in the N condition. In 28.6% of the overloads, the amount of N−1 overload is higher than the overload in the N condition, affecting in particular lines 8432 and 8433, in Voltage Range B. Moreover, global N−1 performances are described through the index 849 = 0.338%.

Analysis of Specific Loading Condition
The line loading levels in the N condition at hour 849 show that 74 lines are in overload in the N condition, and 22 of them, in Ranges B and C, have (b, 849) > 2. The global network performance is described by u R (849) = 9.63% and u + R (849) = 0.095%. As regards N−1 analysis, involving 77.7 × 10 6 cases, synthetic indices are reported in Table 8. Overloads are observed in 0.8% of combinations, affecting 193 lines. Peak overloads, withλ(b, h, 849) ≥ 5, are detected only in 6 lines due to outages close to lines already overloaded in the N condition. In 28.6% of the overloads, the amount of N−1 overload is higher than the overload in the N condition, affecting in particular lines 8432 and 8433, in Voltage Range B. Moreover, global N−1 performances are described through the indexû + R (849) = 0.338%. In Table 9, average N−1 overloadsλ s (b, 849) are divided into classes. More than half of overloaded lines haveλ s (b) ≤ 1.2, considered acceptable in N−1 network planning as described in the above, whereas 23 lines, mostly in overload in the N condition, show values higher than 2. The number of lines in overload due to the outage of line k,n c (k, 849), is detailed in Table 10. In 91.7% of outages, the number of overloaded lines is the same as in the N condition. The maximum of 84 overloads is observed for the outage of lines 1035 and 1036, in Voltage Range A, arising as the most  The performed investigation allows an early detection of critical lines, which could represent the target of a specific security analysis during operational planning, along with risk measures depending on device failure rates and further technical details.

Conclusions
A flow-based procedure for the preliminary security analysis of a yearly network evolution scenario has been proposed. Starting from hourly load and generation by market solution, linear PTDFs and LODFs have been exploited to assess N and N−1 power flows and loading levels. In particular, an extended matrix formulation has been carried out to efficiently handle N−1 islanding. Statistic indices have been defined to inspect network behavior in the big data framework. The procedure has been tested on a model of provisional European high voltage network involving thousands of nodes and lines. Simulations have proved the effectiveness of the tool to handle a real-sized network with a plain formulation, and to point out the security effects of varying loading conditions and events, suiting multiple post-elaborations. The proposed method has proved to be a powerful and flexible tool for TSO long-term and secure operation analysis, even in the presence of a remarkable amount of data, leading to point out critical conditions to be further inspected during system operation, and critical branches to evaluate grid developments. The tool could be useful for other entities involved in power systems as well, such as generation companies, traders and regulators, in order to draw, by means of the defined indicators, useful technical-based signals to evaluate new initiatives or strategies in energy and service markets. Future work could aim at regulation strategies for reducing overloads, and at their application to network expansions.

Appendix A. Formulation of Network Operation Indices
Indices for operation in N conditions: Indices for operation in N−1 conditions: λ c (k, t) = 1 n c (k, t)