Juan Sun
email: juan.sun@ulb.be

Pierre-Etienne Labeau
email: pierre.etienne.labeau@ulb.be

Probabilistic Assessment of the Connection Capability of Distributed Units in a Power System under Active Network Management Estimation Probabiliste de la Capacité de Connexion d'Unités Distribuées à un Réseau Électrique en Gestion Active

Keywords: active network management, distributed generation, congestion risk, curtailment, efficient Monte Carlo sampling gestion active du réseau, production distribuée, risque de congestion, curtailment, techniques d'échantillonnage Monte Carlo efficaces

The variability of power supply from distributed generation (DG) units results in various issues when incorporating more DG units to the existing distribution networks, particularly congestion risks. Active Network Management (ANM) could provide (almost) real-time control, by possibly curtailing their production in case of grid congestion so as to allow more DG units integration, while deferring costly and time-consuming network upgrades. This paper provides a methodology for the fast assessment of the connection capability of DG units to a grid, based on efficient Monte Carlo sampling. Besides, based on correlated sampling, it is possible to simultaneously estimate the congestion risk with and without connecting a new DG unit of variable capacity. This significantly reduces the computation burden in assessing the connection capability of a grid. The effectiveness of the proposed method is demonstrated on a test power grid.

I. INTRODUCTION

In the ongoing energy transition, the amount of decentralized production connected to the Medium-Voltage (MV) power system is increasing. As a consequence, when the power produced is not consumed locally, reverse power flows are injected into the High-Voltage (HV) grid. More frequent line congestions and voltage problems are thus likely to occur. To overcome this issue, the network could be reinforced, but it is an expensive and time-consuming task that Transmission System Operators (TSOs) are willing to plan optimally or to postpone. Additionally, this reinforcement would be useful, in a first stage at least, only in low-probable situations, associated with a large amount of generation and low loads. Alternatively, in order to optimize both the use of the present grid infrastructure and the number of connected Distributed Generation (DG) units, an Active Network Management (ANM) scheme can be envisioned [START_REF] Järventaustac | Smart grid power system control in distributed generation environment[END_REF]. It aims to control the injection of energy produced by DG units to the grid, in almost real time, by possibly curtailing their production in case of grid congestion. These smart connections are also known as non-firm or interruptible connections offered to DG customers. In such a way, the grid is always retained secure, while it can provide flexible connection to facilitate the integration of more DG units in the grid.

As for the assessment of the congestion risk entailed by new DG unit connections, the main drawback of the previous deterministic "fit-and-forget" approach is that conventional (or firm) connections assume that the generators are able to output full capacity. The capability of MV buses to accept power injection from DGs, which are in form of Firm Generation (FG), is largely determined by the local minimum load and the export circuit capacity during the contingency of an outage of the highest-rated distribution circuit (i.e. classical N-1 criterion). This approach to network operation means that much of the distribution network exists without monitoring and control and has spare capacity during normal operation [START_REF] Currie | Active power-flow management utilising operating margins for the increased connection of distributed generation[END_REF]. In contrast, flexible connections (i.e. Non-Firm Generators, NFGs), combined with a set of rules about the order they are dispatched or curtailed under an ANM scheme, provide a solution preferred in particular for transmission network zones to make an efficient use of renewable generation whenever possible and allow connection of larger generators.

In fact, the specific curtailment rules are the Principle of Access (PoA) rules, chosen by system operators following different approaches in achieving desired objectives, such as cost-effectiveness and social optimality [START_REF] Anaya | Experience with Smarter Commercial Agreements for Distributed Wind Generation[END_REF]. The PoA rules of Last-In-First-Out (LIFO) and Pro-Rata schemes are the main two rules applied in practical settings to guide the generation curtailment. The LIFO rule means that the first NFG to sign a contract in the ANM scheme is always the last one to be curtailed. This rule is easily implemented and does not affect existing generators, but it might discourage investment in future DG units in the same area. On the other hand, in the Pro-Rata rule, the required curtailment is shared between all NFG units, proportionally to the rated capacity or actual output of the generators. It also signifies that revenue losses are shared equally amongst all curtailed generators. This plays a key role in maximizing the renewable generation capacity built, which will lead to a higher network infrastructure utilization.

With an ANM arrangement deploying different PoA schemes among the connected units, more DG units are likely to be incorporated into the distribution power system. However, the impact of ANM must be assessed in order to guarantee, for each DG unit's owner, a sufficient profitability. The latter can be measured by the Utilization Factor (UF) of each unit. The UF accounts for the expected effect of curtailment on the more classical Capacity Factor (CF), i.e. the ratio of the actual production over the ideal production the DG unit would generate if it was operating all the time at its installed power. Due to the variability of DG units, a probabilistic methodology should be used, in order to propagate uncertainties on the loads and generations, on the grid model, while using an Optimal Power Flow (OPF) to estimate the most economical curtailment resolving a case of congestion when observed. The congestion risk will be reflected by our proposed probabilistic Risk Indices (RI), associated with power curtailment and the probability of congestion occurrence. Therefore, the impact of a new DG unit integration to the existing power grid with an ANM scheme is evaluated by the evolution of these indicators: RIs and UFs (to be compared with the CFs) for all DG units. If the UF value of each generator still lies within a pre-defined guaranteed range, the operator will accept this smart connection strategy for a new DG unit incorporation; otherwise, the network reinforcement will become unavoidable, what always leads the TSO in practice to reject the connection request of the new DG unit.

Consider an existing power grid operated with an ANM scheme, where there are q DG units connected to n substations. The following assumptions on the current power grid are introduced: 1) the bulk power system is assumed to be fully available and capable of supplying at any time what is not produced by DG units; 2) all DG units are connected to the MV side; 3) no lack of energy supply is expected (as in reliability assessment), only curtailment costs of solving congestions due to too high production with respect to load ; 4) DG units are assumed 100% available (no scheduled maintenance or outages for the DG units are considered, as the assessment is done from the TSO standpoint); only N-1 contingencies as a result of an outage of transmission elements, e.g. the lines and the transformers in the substations or between different substations. Technical reliability issues are therefore not considered. From the viewpoint of TSOs, they want to get a quick evaluation on how much installed capacity of this new DG unit is most suitable to a given node in their present grid, while not affecting too much the economic benefit of Distribution System Operators (DSOs) and DG unit owners. Along this line of consideration, there will give rise to some potential questions of interest:

• What is the impact of a (q+1)th DG unit, of a given type and of given installed capacity, to be connected to a specific node of the grid, on the previous value of the risk indices of the grid with q DG units?

• What is the connection capacity at a given node of the grid of a specific type of DG unit? Such questions should be treated while accounting for each season (summer, winter, inter-season) associated to different ratings of the grid elements, and for each grid topology considered (N and N-1 conditions on the grid elements).

Consequently, a methodology was proposed [START_REF] Labeau | Mathematics of PRA applied to Distributed Generation curtailment in saturated grids[END_REF][START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF] for the fast assessment of the impact of connecting a given type of DG unit with given capacity to the existing grid without any network reinforcement or upgrade, under an ANM scheme. A typical grid loop is considered, to which wind farms (WFs) and CHP (Combined Heat and Power) units are connected. Probability density functions (pdf's) for their respective generation as well as for the loads at the different nodes of the grid are first elicited. As a direct Monte Carlo Sampling (MCS) of these pdf's turns out to be prohibitively time-consuming, an alternative approach is used: possible congestions have to be analyzed, not on the basis of the detailed generations and loads, but on the Net Balance (NB) at each node of the grid. This net balance is the algebraic sum of all productions minus the total load in a substation. A domain of the NB space that is free from any congestion risk can be identified. In order to do so, the net balance space is discretized, and safe and unsafe cells resulting from this discretization are then obtained, by checking if an acceptable solution of the Load Flow equations exists at each corner of those cells. The optimal curtailment of DG units must on the contrary be calculated using the detailed variant of generations and loads associated with any unsafe net balance variant. This is achieved based on a targeted (systematic) importance sampling of only those detailed variants (all individual productions and loads) likely to lead to congestion, i.e. those corresponding to NBs belonging to unsafe cells [START_REF] Labeau | Mathematics of PRA applied to Distributed Generation curtailment in saturated grids[END_REF][START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF]. This work provides foundation for our research. Hence, risk indices can be calculated efficiently for a given situation of the grid with q connected DG units.

However, this method provides only an assessment of the ANM scheme performance for an existing situation. If a new DG unit is likely to be connected to the grid, its impact on the performance indicators mentioned above is much harder to estimate. Let us denote the case with q+1 DG as reference case, while the initial case without any new DG integration (i.e. with q DG units) will be the perturbed scenario. In order to estimate e.g. the evolution of the UF of each unit caused by the addition of one specific DG unit at a given node, it would be necessary to independently compute estimations of the indicators for both the reference and perturbed scenarios with a high accuracy, as the difference between those indicators might be small. This leads to a significant increase in computing load. But resorting to a Correlated Sampling (CS) method, a unique computation will simultaneously bring the risk estimations for the reference case and for the perturbed case (q+1 and q DG units, resp.), with an acceptable accuracy. This idea can be generalized: by considering a fictitious new DG unit to be connected to a given node with a variable installed power, the connection capacity can be derived from the same type of computations in a limited amount of time.

The paper is organized as follows. Section II defines the risk indices and the mathematical theory supporting risk estimation for different scenarios. Section III presents the application of this methodology to a grid test case with realistic data in order to illustrate the approach. Results and analyses are given in Section IV and conclusion is provided in Section V.

II. DEFINITION AND CALCULATION OF RISK INDICES

A. Definition of Risk Indices

We adopt the classical definition of risk as the product of the probability of occurrence of the undesired event (i.e. contingency) and the related consequence (i.e. severity) [START_REF] Rocchetta | Risk assessment and risk-cost optimization of distributed power generation systems considering extreme weather conditions[END_REF]. To take into account all undesired events, the definition is extended by summing all contributions as: ( ) ( )

ii i RI P E f E =  (1)
Where P(Ei) is the probability of occurrence of the undesired event Ei and f(Ei) is the severity of the related consequence. In probabilistic risk assessment, contingency frequencies are used as probabilities and severity functions as consequences. In the context of power systems, a contingency is defined as the unexpected loss of one or more electric elements (e.g. line, cable, or transformer) in the power system [START_REF] Gooding | Probability analysis of distributed generation for island scenarios utilizing Carolina's data[END_REF]. Overload, influenced by seasons, related with the feeders' thermal limits, and bus voltage magnitude, related with frequency and system balance, are both indicators of power system stress and and are used to represent the consequences for the risk calculation [START_REF] Rocchetta | Risk assessment and risk-cost optimization of distributed power generation systems considering extreme weather conditions[END_REF]. Thus, the total risk index (i.e. sum of all seasons) associated with one contingency can be expressed as follows for the whole power network:

1 ( ) ( ) ( , ) s n k l l k l k l RI C q P C f C    = =    ( 2
)
where χ is the set of all seasons conditions (e.g. summer, inter-season, winter). Ck is the kth contingency (e.g. N/N-1 configuration). ql is the probability of the lth season. fl(Ck,χ) is the performance function under congestion (e.g. the amount of power to be curtailed) for each generator in the conditions of the lth season. Pl(Ck ̸ χ) is the probability of occurrence of congestion in the lth season. ns is the total number of seasonal conditions, nc is the total number of configurations. Let pk be the probability of configuration k. The total risk indices due to all contingencies are then obtained as:

11 ( ) ( , ) cs nn k l kl k kl k kl RI p q P C f C  == =     (3)
An estimation of RI can be done by using Monte Carlo sampling to propagate uncertainties on loads and productions on the OPF calculation module in all possible cases, as depicted in Fig. 1. Yet a direct Monte Carlo algorithm turns out to be inefficient given the dimensions of the problem and the unlikely situations where curtailment has to be carried out.

B. Application to a Given Grid with Given Connected Units

If our power grid comprises n substations, to which WFs and/or CHP units are connected, the total load in substation i is denoted Li. ni is the total number of DG units connected to node i. The production of the jth DG unit connected to node i is denoted Pij. The dimension of these production and consumption vectors is related to the total number of DG units and substations. The more DG units connected, the higher the dimensions of the problem to be analyzed, what will increase the computational burden. We use the concept of Net Balance (NB), defined as the algebraic sum of all generator productions minus the total load in a substation, and calculated with (4) at node i. These NBs correspond to the power possibly exported to the main grid at each substation. They correspond to the quantities likely to cause congestions in the HV grid, independently of the details of the productions and load behind. Identifying possible congestion situations based on the NBs will therefore dramatically decrease the dimension of the problem analysis, and the dimensionality will be only related to the total number of substations. These NB variants will naturally form the NB space, where a boundary between uncongested and congested regions appear. In order to estimate where this boundary lies, the NB space will be discretized into cells, approximating both safe (i.e. no congestion occurrence) and unsafe regions (i.e. congestion region). The congestion region in the NB space is identified by load flow calculations at each corner in those cells, and a cell is deemed unsafe as soon as the NB variant of one of its corners leads to congestion. In a second stage, only those detailed variants (all individual productions and loads) likely to lead to congestion, i.e. those corresponding to NBs belonging to unsafe cells, will be sampled as inputs of OPF calculations providing the total curtailment cost minimization. The curtailed power of each DG unit under specific configuration will be deduced from the OPF for each sampled variant. These variant samples are generated through a targeted (systematic) importance sampling approach. At last, the congestion risk in a given configuration will be estimated by averaging the curtailment power weighted by the probability of congestion occurrence.

Accordingly, the performance function fkl will be the curtailment power of each DG unit in a specific season and a specific configuration (e.g. N/N-1 configuration), and evaluated by OPF computation, which are associated with different PoA arrangements and the curtailment cost of each DG unit, illustrated in Fig. 1. Pkl is the probability of congestion happening in a specific season and a specific configuration (e.g. N/N-1 configuration), which are related with the joint pdf of WF generation, CHP production and load consumption.

Thus, the total risk under all contingencies under consideration in a grid can be expressed as [START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF].

where 𝑃 ̅ and 𝐿 ̅ are the vectors of the productions of the DG units and of the loads connected to the grid nodes, respectively. These two vectors compose the so-called variants of the problem, i.e. the various cases of production and consumption in the grid. 𝜑 𝑘𝑙 (𝑃 ̅ , 𝐿 ̅ )is the joint pdf of the variants in configuration k and season l. 𝑓 𝑘𝑙 (𝑃 ̅ , 𝐿 ̅ ) represents the severity function associated to the curtailment solving the possible congestion associated with this variant. For a given situation kl, (5) can be simplified as:

where 𝜃(𝑁𝐵 ̅̅̅̅ ) is as an indicator function: if 𝑁𝐵 ̅̅̅̅ causes a grid congestion, it is equal to 1; otherwise it is 0. This function divides the NB space into two regions, a safe region where there is no congestion and an Unsafe Region (UR) where congestion exists. As mentioned above, the UR is delimited based on a discretization of the accessible region in the NB space. Hence the UR consists of ncel Unsafe Cells (UC). Replacing the integral on the NB space by the summation (8) on all UCs, equation ( 7) can be written as ( 9) and ( 11).

,min , , ,max

UC UC i w i CHP i i i NB P P L NB  + - (10)
where Pw,i is the total production of WFs connected to substation i; PCHP,i is the total generation of CHP units located in node i; Li is the total load at substation i. 𝑁𝐵 𝑖,𝑚𝑖𝑛 𝑈𝐶 and 𝑁𝐵 𝑖,𝑚𝑎𝑥 𝑈𝐶 are the extreme values of net balance for unsafe cell UC for the dimension corresponding to node i.

The loads at the different nodes are correlated, and so are the wind productions, but correlations between wind productions and loads are negligible. CHP productions are independent of the other variables. Therefore, the total joint pdf of the full variant, i.e. 𝜑(𝑃 ̅ 𝑤 , 𝑃 ̅ 𝐶𝐻𝑃 , 𝐿 ̅ ), consisting of the product of the joint pdf 𝜑(𝑃 ̅ 𝑤 )of all wind productions, the individual pdf's 𝜑(𝑃 𝐶𝐻𝑃,𝑖 ) of the CHP units connected to substation i, and the joint pdf 𝜑(𝐿 ̅ ) of all loads, is represented by (12). Then our risk index is calculated as follows, , 1

( , , )= ( ) ( ) ( ) n w CHP w w L CHP i i P P L P L P     =   (12)
(P , P , L) (P , P , L) P P L

w CHP w CHP w CHP UC RI f d d d  =   (13)
An estimate of RI is obtained when loads and wind generations are sampled from their respective joint pdf, i.e. 𝜑(𝑃 ̅ 𝑤 ),𝜑(𝐿 ̅ ), while the sampling of the CHP generations is forced to take place from the truncated intervals inside a given UC. The CHP production samples are generated from each UC reachable from the partial variant made of wind generation and load. However, it is still likely that many samples of 𝑃 ̅ 𝑤 and 𝐿 ̅ will not be compatible with any output of the existing CHP productions, for specific UCs [START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF]. Yet a set of CHP productions could be compatible, for the sampled values of 𝑃 ̅ 𝑤 and 𝐿 ̅ , with other UCs. Therefore, after sampling wind production and load, CHP contributions in various UCs can be obtained. The risk index estimate will write as:

, , , , 1 1 1 [
(P , P , L ) ( , )] n N w j CHP j j CU i wi j ij j UC i RI f P P L N = = =    ()
where the sampling size is N; 𝑃 ̅ 𝐶𝐻𝑃,𝑗 is the CHP variant forced to fall into this UC, given 𝑃 ̅ 𝑤,𝑗 and 𝐿 ̅ 𝑗 ; PCU,i means the probability of the existing CHP production falling into a given UC at the ith substation. In this way, each variant 𝑃 ̅ 𝑤,𝑗 and 𝐿 ̅ 𝑗 might lead to contributions in all UCs: either the CHP variant is incompatible with the UC, or it leads to a congestion case that contributes to the estimate with the output of the OPF evaluation, i.e. 𝑓(𝑃 ̅ 𝑤,𝑗 , 𝑃 ̅ 𝐶𝐻𝑃,𝑗 , 𝐿 ̅ 𝑗 ) , weighted by ∏ 𝑃 𝐶𝑈,𝑖 (𝑃 𝑤𝑖,𝑗 , 𝐿 𝑖𝑗 ) 𝑛 𝑖=1

.

C. Impact of a New Given Connection to a Given Node

A trend towards a growing integration of DG units will urge the operators to consider when a new DG unit is likely to be connected to a substation, and its impact on the previously assessed risk indices for the existing power grid should be estimated. Suppose the existing power grid is managed with an ANM scheme, where there is a total number of q DG units connected to n substations. If a new DG unit is connected to the power grid, this case is treated as reference case (i.e. ref case) in our study, with a total number of DG units equal to q+1. The initial case without any new DG integration will be seen as perturbed scenario (i.e. per case). The reason for this counter-intuitive qualifications will appear below. The UR related to (q+1) DG units is an extension of that obtained with q units, in the dimension corresponding to the node where this new unit could be connected. Of course, the joint pdf's of the uncertain variants, i.e. 𝜑(𝑃 ̅ 𝑤 , 𝑃 ̅ 𝐶𝐻𝑃 , 𝐿 ̅ ), are affected by the addition of the new DG unit.

We will resort to Correlated Sampling (CS) in order to more accurately assess the supposedly small difference in the risk indices between two situations (ref and per) in one single simulation. CS consists in applying the same sequence of random numbers for the two simultaneous estimations (ref and per) and to introduce compensations in the results of the per case with weighting factors, in order to keep its estimate unbiased [START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF]. CS will not only ensure the estimation accuracy (rather than summing the statistical uncertainties of two independent estimates on the estimation of their difference), but it also avoids repeating the whole estimation process and reduces the computation burden. So with the help of CS, we can simultaneously obtain the results in the ref and per scenarios.

It means that all runs of the OPF module, i.e.𝑓(𝑃 ̅ 𝑤 , 𝑃 ̅ 𝐶𝐻𝑃 , 𝐿 ̅ ), can be just calculated once, both for the ref and the per cases. The risk indices for the ref case will thus be estimated according to (13) and ( 14), as described in ( 15) and (17), if the sampling size is N; while the indices associated to the "q situation" (per case) can be derived from the ref case based on CS, changed as in ( 16) and ( 18

   =    () , , , , 1 1 1 [
(P , P , L ) ( , )] n N ref ref w j CHP j j CU i wi j ij j UC i RI f P P L N = = =   () ,, ,, ,, 1 ,, 1 1 ( , , ) [ (P
  = = =      ()
Indeed, the risk indices assessment for the per situation are computed from the ref case with the correction factor φ per (P ̅ w,j ,P ̅ CHP,j ,L ̅ j ) φ ref (P ̅ w,j ,P ̅ CHP,j ,L ̅ j )

. The correction factor is the ratio between the probability of the jth sampled variant in the per case over that in the ref case. In the detailed calculation process of the correction factors, the correlated DG production sampling vector in the ref case with q+1 units corresponds to that with q DG units in the per scenario.

For the sake of completeness, the previous developments and the estimation procedure of the risk index implicitly assume that there is a load and at least a WF and a CHP unit connected to each node of the grid. The sampling procedure and the calculation of the correction factors can however be adapted when these assumptions are not met.

For instance, assuming there are WFs, CHP units, and load consumption installed in this given connection node, the correction factor calculation is specially defined for two scenarios:

• the added DG unit is CHP, the corresponding correction factor will be represented by .

Once a new CHP unit is connected to a given substation without any previous CHP generator, the corrector factor between the per and ref cases will be equal to 1. As a consequence, the risk indices for the per scenario, i.e. the present power grid before new CHP unit integration, will be estimated in accordance with (19).

- = = =     ()

D. Assessment of the Connection Capacity of a New DG unit

In order to select the optimal installed capacity of a generator incorporated to the current power grid, we take a fictitious DG unit of a given type, to be connected to a given node in the existing power grid and consider various values of its installed capacity. Resorting to the CS algorithm, we can estimate the evolution of these indicators (i.e. the RI, UF of each DG unit) in presence of this new varying installed capacity. Ultimately, at some point, we decide, based on the evolution of these indicators, which installed power is acceptable in the present power grid with the ANM scheme.

For example, could a 8 MW of WF be accepted to a given node of the existing power grid without any network reinforcement? From the TSO standpoint, we are still trying to look at the connection capability, i.e. whether a new project of 5, 10 MW of a given type could be connected without too much additional impact to the existing power grid and DG units, so that they have an idea of how much capacity could still be connected to the grid at a node if a request is sent to them. Accordingly, the case with a new 10MW WF integration is taken as ref case, evaluated through (15) and (17), while the projects of 8, 5 MW are considered as per scenarios, i.e. the impact of a smaller nominal capacity DG unit, connected to the existing grid, could be derived from the ref case by ( 16) and (18). So far, the best optimal capacity generator will be selected under the flexible strategy of ANM considering a specific PoA arrangement.

III. APPLICATION OF THE METHODOLOGY

The structure of the proposed methodology mainly consists of three parts: preprocessing, sample generation and risk estimation [START_REF] Sun | Probabilistic Impact of the Connection of a New DG unit in a Power System Subject to Active Network Management[END_REF], as shown in Fig. 2. The accuracy criteria of the results will be also provided. The preprocessing is actually for the network constraint event analysis in the ref case in order to identify the UR based on the discretization of the NB space; then the corresponding systematic power production and load consumption samples from these congestion areas (or URs) in the NB space are generated through a systematic sampling approach, for curtailment calculation by OPF, considering the seasons, N/N-1 configurations, the PoA options, and the curtailment cost of each DG unit. In parallel with the OPF evaluation, the weighting factors in each UC (i.e. the probability of CHP samples falling into each UC), i.e. the correction factors calculation in each situation, will be computed; finally we can obtain the risk indices in the ref case in different conditions, as well as in per scenarios to determine a suitable capacity for a newly installed DG unit to the existing power grid. Each part is explained in details in the following paragraphs.

A. Preprocess

As shown in Fig. 2, the preprocessing stage, structuring the variants sampling, aims to identify the UR based on the discretization of the NB space, and to partition the NB space into safe and UCs. In this way, the UR in the different configurations may be constructed through these corresponding UCs. At this stage, three steps need to be handled.

• Data collection. A SCADA system performs the network monitoring, which collects every 15 minutes the actual generations and loads. In this way, the NB domain can be determined, based on ( 4) and [START_REF] Shayne | Handbook of toxicologic pathology[END_REF]. Fig. 3(a) shows the NB space with 2 substations. The dots are NB variants (either safe or unsafe situation).

• Identification of possible congestions. A mesh is defined on the NB domain, see Fig. 3(a), and the space is then divided into a set of hyper-cubes (cells). The situations associated to the corners of the generated cells are analyzed by an AC -load flow so as to identify possible congested cases.

• Determination of the UCs. If all the corners of a cell are safe (i.e. if they cause no congestion on the grid), then this cell is considered as a safe cell; otherwise it is considered as a UC. The UR in the NB space consists of the corresponding UCs, as shown in Fig. 3(b). In the meantime, we obtain the discretized approximation of the security border by checking the corners of the cells in the mesh. As shown in Fig. 3(a), a security boundary divides the NB space into the corresponding safe and unsafe regions. Adding a DG unit capacity to node 1 corresponds to extending the NB domain, without affecting the security boundary.

B. Sample Generation Process

Once the congested region has been identified and discretized in the NB space, how to make sure, in the estimation of risk indices, that enough sampled variants go into these unsafe We present a validated systematic sampling scheme for estimating the probabilities of 23ème Congrès Lambda Mu de l'IMdR au 13 octobre 2022, EDF Lab Paris Saclay violation of operating constraints due to congestions in grid elements.

As illustrated in Fig. 4, from the DG production measurement data, we find the stochastic dependence structure of wind units through correlation coefficient matrix, and we model Weibull marginal distribution probability density function (pdf/PDF) for wind speed. Based on these marginal distributions, their joint pdf is constructed thanks to a Gaussian Copula. Since Gaussian copula is easy to implement for high-dimensional variables, it seems to be a good choice for modeling multivariate wind distributions [START_REF] Leuthod | When the wind blows over Europe -a simulation analysis and the impact of grid extensions[END_REF]. Then, a Monte Carlo sampling approach is applied to generate a large number of correlated wind speed samples, and sampled wind speeds are converted to electric power using the power curve. As for the load model, it is obtained by normalizing the historical load data of each year to the maximum consumption of the corresponding substation in that year. After sampling wind and load from the corresponding joint pdf's, significant values of CHP productions correspond to values likely to make the NB fall in a UC [START_REF] Labeau | Mathematics of PRA applied to Distributed Generation curtailment in saturated grids[END_REF][START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF][START_REF] Sun | Probabilistic Impact of the Connection of a New DG unit in a Power System Subject to Active Network Management[END_REF]. So we check the possibility of CHP samples that can contribute to each UC. For each UC, shown in Fig. 5, a compatible value of the CHP production at each node i is sampled from the corresponding truncated pdf, if possible, denoted by (20). Finally, all related samples under each configuration (i.e. seasons, N/ N-1 contingencies) will be produced. represents the joint pdf of the CHP generation at node i. P CU,i is probability of having the CHP production at node i in the right range. At this step, all related production and consumption samples under each configuration have been generated.

C. Risk Estimation

In parallel with the sample generation, the weighting factors in each UC, i.e.𝑃 𝐶𝑈,𝑖 , will be computed. It actually means the probability of CHP samples taken from its truncated pdf inside this specific cell. And it is inevitable to take the weighting factor calculation of each situation (e.g. season condition, N/N-1 contingency) into account.

An OPF tool, whose objective is the curtailment cost minimization, is applied to the generated variants in order to calculate the (possible) curtailments of all DG units. The curtailment calculation of each DG unit considers seasons (for line ratings), N-0/N-1 contingencies, the curtailment cost of each specific DG type as well as its type of contract (PoA).

As for the correction factor calculation, it is expressed in these 4 steps, adopted to decrease the curtailment computation time for the per case.

• Identify the scenarios of perturbation (perturbed pdf's); the type of the added DG unit (i.e. WF or CHP) and the type of connected substation (i.e. no WF, no CHP unit or no load consumption) must be specified in order to generate the correlated sample in the per case;

• Sample variants from the joint pdf of the reference case;

• Calculate the corresponding correction factors for each ref variable;

• Estimate the risk indices of the per case in the different situations (i.e. seasons and N-0/N-1 configuration).

Finally, the risk indices in the ref case can be estimated by averaging the corresponding curtailments weighted by the probability of congestion occurrence. The risk indices for the per case are then expressed from (18) as: , , ,

(P , P , L ) ( , )] CHP i n per CHP ij per per N wj j w i L n ref ref w j j ref j w L CHP ij i n UC UC W i i CU i wj i ji UC i P P L RI N P L P f P P L       = = = = =         () 1 , 1 , 1 , , , 1 ( ) ( ) () 1 [ ( ) ( ) ( )

D. Stopping Criteria

A Student's t-test is most commonly applied when the test statistics would follow a normal distribution if the value of a scaling term in the test statistics were known. Every t-value has a p-value to go with it. A p-value is the probability that the results from your sample data occur by chance, from 0% to 100%. In most cases, a p-value of 5% is accepted to mean the data is valid. Seen from Fig. 6, the confidence level is 95%, i.e. 1-α=95%, that means the t-value under some interval (i.e. 𝑡 ∈ [-𝑐, 𝑐]) corresponds to a p-value of 95%, denoting c as the 97.5th percentile of this distribution, which are obtained by the t test table. Fig. 6. 95% confidence interval by Student's t-test [START_REF] Shayne | Handbook of toxicologic pathology[END_REF] In power system reliability assessment, the coefficient of variation an index is often used a convergence criterion or stopping rule for MCS methods [START_REF] Jonnavithula | Composite System Reliability Evaluation Using Sequencial Monte Carlo simulation[END_REF]. The simulation break occurs if the coefficient of variation reaches the pre-specified tolerance value. We used ε as the coefficient of variation and it is the ratio between the variation from the average RI and the RI mean. Once the coefficient of variation will be less than 10%, the fluctuation on the mean of the risk will decrease and can stop this simulation, that will ensure the accuracy of risk indices results. Thus, the definition of ε, and its values in different scenarios are estimated from:

() c D RI N RI  =  (22)
Where there is sample size N; D(RI) represents the variance of RI; the RI sample mean is RI. To test the validation of this application method, a threenode system with 5 wind farms (2 non-firm wind generation, WF2A and WF2B), 2 CHP units and 1 PV unit is used as a case study. The simplified topology of this reference case (i.e. Ref. Case) is shown in Fig. 7, that corresponds to connecting a new 20MW wind generation (i.e. WF2C) to the NODE 2, while consider two PoA arrangements: LIFO and Pro Rata.

IV. RESULTS AND ANALYSIS

A. Three-node System

The detailed installed capacity of DG units and the range of gross of load consumption in each substation are listed in TABLE І. The previous situation of the power grid (i.e. Per 1. case) and another scenario of a new 10MW installed capacity wind unit accepted at NODE2 (i.e. Per 2. case), both are per cases, in which could be assessed the congestion risk from the ref case based on CS. In addition, the solar generation here is assumed to be constant. Actually, after the analysis of preprocessing, the overload constraint happens in Node 2 due to DG outputs exceeding the total ratings of the T8 and T9 transformers in normal situation. If the additional installed wind unit (i.e. WF2C) is NFG, it will abide by the priority order of the generation curtailments under LIFO scheme, where WF2C will be firstly curtailed and followed by WF2A, and finally WF2B, depending on the total amount of curtailment needed. By contrast, with Pro Rata arrangement, all flexible wind units will be at the same time curtailed based on the size of their own output. TABLE II lists the conditional curtailment risk values of the ref case in summer under different N-1 configurations with the LIFO PoA. In addition, if the system loses one transformer connected to NODE2, i.e. (N-1)8 or (N-1)9 contingency, it will need more curtailment than for the other N-1 situations. Analyzed from Fig. 8, the last connection NFG of WF2C is the first and most curtailed under LIFO arrangement, be it in the ref case or in the Per 2. Case. But with the Pro-Rata rule, three NFGs are allocated fairly to relieve the constraint. By the way, the results of total curtailment risk are multiplied by the probability of different seasons (summer, inter-season and winter) and different N/N-1 contingencies on the foundation results of conditional curtailment risk. It seems that the 10MW wind unit (i.e. Per 2. case) is more suitable to connect to the existing grid, either in the LIFO or Pro-Rata schemes.

TABLE III provides the decrease of the CF of each unit, resulting in their UF, under the LIFO and Pro-Rata rules in the Per 2. Case. The highest difference, obtained for WF2C, is 0.0692% in the LIFO scheme. And under the Pro-Rata rule, WF2B and WF2C display the largest curtailment. Their reduction are both 0.02%. Since the TSOs would like less DG units to participate in the curtailment control, they would in this case prefer continuing to apply LIFO contracts for the extra 10 MW wind unit given the lower decrease in CF, which is between the predefined range of allowable values for which they do not need to give compensation for the power curtailment. ) is about 12%. λ is the variation from the average RI. If the initial reference sample size is increased, the fluctuation on the mean of the risk index can decrease and the convergence will be improved. For example, if the sample size of the ref case is set to 5000 and 15000, resp., then in the Per 2. case, the sample size is 3139 and 9394, resp., and the coefficient of variation of the RI mean is decreased to 7.2% and 4.5%, resp. This satisfies the criterion for the accuracy and the convergence of the risk indices calculation. The Belgian local TSO launched a programme in the Eastern part of Belgium, where there is a high potential for DG deployment. A typical loop network of such an area is represented in Fig. 9. The objective is to enable new DG units to be incorporated in the base network by using feasible ANM solutions without any network reinforcement/upgrade, whilst preventively addressing congestion problems caused by connecting some new power injections into the initial transmission system designed to meet lower rural demand.

B. Typical Loop Network

The current situation without accounting for new DG unit incorporation is defined as Per 1. case, where there is a total of 9 WFs, 12 CHP units and 9 other generators (e.g. solar photovoltaic generation and hydro generation plant) installed at 7 substations. Given the limited share of these other production types, the latter are treated as a constant value in each substation. Via the load flow calculation analysis in each PoA Factors configuration (i.e. seasons and N/N-1 configurations), the locations of active constraints in the Per 1. case are identified and marked in Fig. 9 by ①, ②, ③, ④, ⑤, ⑥, ⑦. Actually, there is only one constraint occurrence under normal/N-0 conditions and in some seasons (i.e. summer and inter-season), that is at location ① marked with dark yellow. The remaining constraints signaled with light yellow ② ~ ⑦ can only occur in some N-1 configurations in specific seasons. Particularly, the constraints located in ① ~ ④, where these branches are connected to the mainland bus, all NFGs, may participate in the curtailment estimation (or OPF calculation process) within the Pro-Rata rule.

Let us select a new 10MW WF unit (that is WF10) added to bus 13 as Ref. case. The constraint analysis identifies a new possible congestion at place ⑧under some N-1 conditions in a given season, in addition to the existing constraints at ① ~ ⑦. TABLE V exhibits the priority order of these flexible generators for power curtailment under a LIFO arrangement in case of a NFG WF10 incorporation. TABLE VI presents the expected values of the power curtailment risk of WF10 in different seasons under each N/N-1 configuration, and these values imply the average power to be curtailed from the DG units in each situation/condition. Due to high electricity demand from users and customers in winter, the risk of power curtailment in winter is the least. Fig. 10 shows the total energy curtailment of these generators undergoing possible curtailment. In the LIFO scheme, WF10 will encounter almost 16 MWh/year energy curtailment, while all the others curtailed DG units will share a 7.5 MWh/year cut. Nevertheless, for the Pro Rata rule, each of the WF4, WF7, WF10 and CHP16 units will undergo approximately 2 MWh/year of energy curtailed. The remaining 2 MWh /year is shared by the other units, based on the ratio of their power output to the total power curtailment. With the help of the CS approach, we could estimate the congestion risk in Per.1 case from that in the Ref. Case. As for the computation time, if set 30000 wind-load pattern samples under the LIFO scheme, it will take 5 days to complete the whole process for both the ref and per scenarios in one computation. If the network operator switches to the Pro-Rata rule, it will be another 2.5 days to run the remaining processes after the step of sample generation. Thanks to the importance sampling scheme with CS, results for both the ref and per scenarios are simultaneously obtained in one computation, resulting in a significant reduction of the computational burden. The key contribution of this research is to propose an efficient methodology for the probabilistic assessment of the risk of connecting a new DG unit to an existing grid, using an ANM scheme. Furthermore, the use of correlated sampling allows significantly reducing the computation time so that one computation is able to simultaneously provide the risk estimation for both the reference case and other perturbed scenarios. For instance, if the TSO needs a quick reply to a possible request of connecting a new DG unit of a given type and installed capacity to a given substation in the current grid, we could provide a whole feedback analysis and flexible strategy related to a possible candidate DG unit integration. Besides, the congestion level due to a new DG unit connection into the current grid is reflected by risk indices and the revenue of each generator could be reflected by its UF, which contains the economic benefits. Last but not least, if the operator needs to implement other PoA arrangements, they could easily insert additional constraints or logical commands to this general methodology.

From the perspective of investors and grid operators, this methodology does not only satisfy the economic expectation of the investment assessment, but it also gives practical information about the export and transmission of the distributed generation in excess. And it also analyzes the impact of the possible acceptance of new DG units, what provides foundation for the future grid planning and its reinforcements.

Fig. 1 .

 1 Fig.1. Schematic diagram for risk index estimation[START_REF] Faghihi | Probabilistic curtailment analysis for transmission grid planning using Active Network Management[END_REF]

  the loads with the NB variants, the risk index formula in a given configuration kl now writes: = ( ( , ), ) ( ) ( ( , ), ) RI f L NB P P NB L NB P P d NBd P  

  f L NB P P L NB P P d P  

  ), ) ( ( , ), ) RI d NB f L NB P P L NB P P d P   

  φ per (P ̅ CHP,j ) φ ref (P ̅ CHP,j ) ;• the added DG unit is WF and the joint pdf of WFs, φ w (P ̅ w ), has been changed, the correction factor is simplified by φ per (P ̅ w,j ) φ ref (P ̅ w,j )

Fig. 2 .

 2 Fig. 2. Flowchart of the methodology

Fig. 4 .

 4 Fig. 4. An example of wind speed marginal distribution

2 Fig. 7 .

 27 Fig. 7. Test grid topology

Fig. 8 .

 8 Fig. 8. Total energy curtailment under different PoA rules

Fig. 9 .

 9 Fig. 9. East Loop simplified topology

  Fig.10shows the total energy curtailment of these generators undergoing possible curtailment. In the LIFO scheme, WF10 will encounter almost 16 MWh/year energy curtailment, while all the others curtailed DG units will share a 7.5 MWh/year cut. Nevertheless, for the Pro Rata rule, each of the WF4, WF7, WF10 and CHP16 units will undergo approximately 2 MWh/year of energy curtailed. The remaining 2 MWh /year is shared by the other units, based on the ratio of their power output to the total power curtailment. With the help of the CS approach, we could estimate the congestion risk in Per.1 case from that in the Ref. Case. Fig.11compares the expected values of the total energy curtailment in the Ref. and Per 1. cases.

Fig. 10 .

 10 Fig. 10. Total energy curtailment under the two rules

  ).

		RI	ref	=	UC 	f	(P , P , L) w CHP			ref	( , w P P L dP d , ) CHP w	P CHP	dL	()
	RI	per		UC	f	(P ,P ,L) w CHP	per ref	( , , ) ( , , ) w CHP w CHP P P L P P L	ref	( , , ) w CHP P P L dP d w	P CHP	dL

TABLE I .

 I CAPACITY AND CONSUMPTION AT EACH NODE

	Substation	Load (MW)	Production (MW) WF CHP Solar
	Node 1	Maximum Minimum	2.5 13.5	10	5	╳
	Node 2	Maximum Minimum	1.6 11.4	52	10	1
	Node 3	Maximum Minimum	0.9 12.2	15	╳	╳

TABLE II .

 II CONDITIONAL RISK OF REF. CASE UNDER LIFO (IN KW)

		DG Units
			WF2A	WF2B	WF2C
		(N-1)0	0	0	9.14
		(N-1)1	6.48	0	407.44
		(N-1)2	6.48	0	404.99
		(N-1)3	0	0	214.94
		(N-1)4	0	0	8.85
	Summer	(N-1)5 (N-1)6	0 0	0 0	8.85 11.35
		(N-1)7	0	0	11.35
		(N-1)8	183.99	1.88	2420.87
		(N-1)9	183.99	1.88	2420.87
		(N-1)10	0	0	9.79
		(N-1)11	0	0	9.79

TABLE III .

 III CAPACITY AND UTILIZATION FACTORS IN PER 2.

CASE LIFO Pro-Rata WF2A WF2B WF2C WF2A WF2B WF2C

	CF(%)	23.793 24.692 24.562 23.793	24.692	24.562
	UF(%)	23.791 24.692 24.545 23.789	24.687	24.557
	Down	0.0084	0	0.0692 0.0168	0.02	0.02

TABLE IV .

 IV CONDITIONAL RISK OF PER 2. CASE UNDER LIFO (IN KW)

	Sample Size	Expected Value	λ	Upper Bound	Lower Bound	λ / RI (%)	per
	1284	675.87	78.76	754.63	597.1	12
	3139	651.14	46.92	698.06	604.22	7.2
	9394	644.38	28.85	673.23	615.52	4.5

  TABLE IV shows the WF2C curtailment risk with a 95% confidence interval in the Per 2. scenario (i.e. a 10 MW wind unit added to the grid under the LIFO rule). When generating 2000 initial samples in the ref case, the sample size of the Per 2. case variants, where more than 10 MW wind power should be filtered, will be decreased to 1284. Then the coefficient of variation of RI mean in the Per 2. case (i.e. ε=λ / RI per

TABLE V .

 V NFG CAPACITY AND INSTALLATION SEQUENCE

	Priority Order	Generator Name	Installed Capacity
	1	CHP12	1
	2	CHP15	0.2
	3	WF6	12
	4	WF3	10
	5	WF8	10
	6	CHP18	2.8
	7	WF9	0.7
	8	CHP22	4.5
	9	CHP19	4.6
	10	WF4	12
	11	CHP16	25
	12	WF7	10
	13	WF10	10

TABLE VI .

 VI WF10 POWER CURTAILMENT IN LIFO RULE (IN KW)

Seasons/(N-1)i


Mounir Bendali-Braham

Jonathan Weber

Germain Forestier

Lhassane Idoumghar

Pierre-Alain Muller

Ensemble classification of video-recorded crowd movements

Keywords: Ensemble classification, Deep learning, Video analysis, Crowd behavior analysis

Ensemble learning methods often improve results in problems addressed by single Machine Learning models. In this work, we apply Ensemble Learning on video-recorded crowd movements. First, we build Ensembles of homogeneous Convolutional Neural Networks (CNN) to compare their performance on the Crowd-11 dataset and show the gain of performance demonstrated by Ensembles compared to single CNN models. Secondly, we evaluate all the possible combinations of these homogeneous Ensembles to build a global Ensemble of heterogeneous models, and we analyze the combination of Ensembles that achieves the best results. Our experiments reveal that Ensemble classification often obtains better results than single models and combining different Ensembles can make the predictions accuracy even better.

I. INTRODUCTION

More and more cities are subject to massive crowd movements due to cultural or political events. During a crowd movement, law enforcement services rely on video surveillance and can count on recent approaches allowing an optimal deployment of cameras [START_REF] Kritter | On the realworld applicability of state-of-the-art algorithms for the optimal camera placement problem[END_REF]. However, although the video data collection is increasingly widespread, the automatic analysis of these videos is not always performed in real time which can delay the intervention of law enforcement services. One of the reasons for this is the lack of generic statistical models that can be used in real time to detect any type of anomalies that may arise from crowd events [START_REF] Thida | A literature review on video analytics of crowded scenes[END_REF]. Furthermore, the development of such models is a tedious task. This can be explained by the scarcity of annotated datasets [START_REF] Carreira | Quo vadis, action recognition? a new model and the kinetics dataset[END_REF]. However, recent years have seen the emergence of datasets illustrating massive and various crowd movements such as Crowd-11 [START_REF] Dupont | Crowd-11: A dataset for fine grained crowd behaviour analysis[END_REF].

The authors of Crowd-11 trained models to classify crowd clips. The model that obtains the best classification results in their article derives from the C3D architecture [START_REF] Tran | Learning spatiotemporal features with 3d convolutional networks[END_REF]. In previous works, we obtained better results [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF], by using a model derived from the TwoStream Inflated 3D architecture (2S-I3D) which already outperforms the C3D models on action recognition datasets [START_REF] Carreira | Quo vadis, action recognition? a new model and the kinetics dataset[END_REF].

In this article, we aim to improve the classification results on the Crowd-11 dataset by applying Ensemble Learning methods. We first create Ensembles of homogeneous models deriving from different architectures and benefiting from different training pre-conditions. The theoretical aspect of this comparison is developed in Section IV-A. Secondly, an evaluation of all the possible combinations of these models allows us to elect a global Ensemble of heterogeneous models that gathers several homogeneous Ensembles. The theoretical aspect of this combination is developed in Section IV-C.

This article is organized as follows: in Section II, we present Ensemble methods applied to video classification and crowd analysis. In Section III, we present our approaches. We discuss our experiments in Section IV.

II. REVIEW

Ensemble methods perform very well in several machine learning tasks [START_REF] Fawaz | Deep neural network ensembles for time series classification[END_REF]. Zhou [START_REF] Zhou | Ensemble learning[END_REF] divides the Ensemble methods into three major categories:

• Boosting, illustrated by its most famous algorithm Ad-aBoost [START_REF] Freund | A desicion-theoretic generalization of on-line learning and an application to boosting[END_REF], which consists in learning T models by associating each time different weights to the training examples. At the beginning, of similar values, these weights change at the t th iteration of the AdaBoost algorithm by taking into account the error obtained from the model trained on t th -1 iteration. In the end, a weighted majority vote is used to combine the decisions of the T models. • Bagging, which is a contraction of Bootstrap Aggregating, where statistical methods are trained on samples created by Bootstrap sampling [START_REF] Efron | Bootstrap methods: another look at the jackknife[END_REF]. Subsequently, these methods are combined into an Ensemble by a majority vote. • Stacking, where different statistical methods are trained on a dataset. Subsequently, a second statistical method, called a meta-classifier, learns to combine the trained models. Furthermore, Zhou considers that some Ensemble methods do not fall into any of these three major categories.

Here we explore some recent Ensemble approaches applied to image processing, video analysis, and crowd analysis.

For image processing and video analysis, Lia et al. [START_REF] Liu | An ensemble deep learning method for vehicle type classification on visual traffic surveillance sensors[END_REF] apply an Ensemble method to provide a solution to classes that lack examples for a vehicle image classification problem. They apply balanced sampling and data augmentation. Their Ensemble method consists of a combination of multiple ResNet models [START_REF] He | Deep residual learning for image recognition[END_REF] and the decision is inferred using the majority vote.

Pouyanfar et al. [START_REF] Pouyanfar | Automatic video event detection for imbalance data using enhanced ensemble deep learning[END_REF] propose EDL (Ensemble Deep Learning) which they use for the classification of videos on the Trecvid [START_REF] Smeaton | Evaluation campaigns and trecvid[END_REF] and Disaster [START_REF] Pouyanfar | Semantic concept detection using weighted discretization multiple correspondence analysis for disaster information management[END_REF] datasets. EDL is a suite of deep feature extractor models from images, which are Convolutional Neural Networks (CNNs) [START_REF] Lecun | Convolutional networks for images, speech, and time series[END_REF] pre-trained on ImageNet, and each extractor is followed by a Support Vector Machine (SVM) [START_REF] Cortes | Support-vector networks[END_REF] which serves as a weak learner in the Ensemble. The learned features are extracted from the last Fully Convolutional Network (FCN) layer of each model. The used architectures are: AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], CaffeNet [START_REF] Jia | Caffe: Convolutional architecture for fast feature embedding[END_REF], Region based CNN [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF], GoogleNet [START_REF] Szegedy | Going deeper with convolutions[END_REF], and ResNet. The decision is inferred following a weighted vote.

Inspired by Liu et al. [START_REF] Liu | Exploratory undersampling for class-imbalance learning[END_REF], Chen et al. [START_REF] Chen | Action recognition using ensemble weighted multi-instance learning[END_REF] propose an Ensemble approach named Ensemble Weighted Multi-Instance Learning. They start by sampling several subsets of the majority class, and by combining each time a subset of the majority class with a minority class, they train a model using AdaBoost [START_REF] Freund | A desicion-theoretic generalization of on-line learning and an application to boosting[END_REF]. The trained models are combined for the final decision.

In the context of crowd analysis, Walach et al. [START_REF] Walach | Learning to count with cnn boosting[END_REF] apply gradient boosting and selective sampling on a simple CNN architecture to perform objects counting in images. Their approach is applied to microscopic bacterial cells datasets, and crowd counting datasets.

Wu et al. [START_REF] Wu | Ensemble learning for crowd flows prediction on campus[END_REF] stack several models whose outputs are used as new features which will be sent as inputs of a new model.

Due to the paucity of annotated data, Gong et al. [START_REF] Gong | Learning human pose in crowd[END_REF] learn in a semi-supervised manner an Ensemble of pose-sensitive DPM (Deformable Part-based Model) mixtures [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF] for pedestrian detection. Several postures are taken into account such as: front, rear, left, right. Each DPM mixture specialized to a specific pedestrian posture is trained using a Latent-SVM [START_REF] Yu | Learning structural svms with latent variables[END_REF].

Contrary to previously mentioned approaches [START_REF] Liu | An ensemble deep learning method for vehicle type classification on visual traffic surveillance sensors[END_REF], [START_REF] Liu | Exploratory undersampling for class-imbalance learning[END_REF], we do not aim to solve the problem of unbalanced data. We do not do Boosting, as in Walach et al. [START_REF] Walach | Learning to count with cnn boosting[END_REF], because in Boosting the training sessions are repeated several times to change the weights of the training examples. The same cannot be reproduced for video data without requiring a huge calculation time. We opt for a compromise between a form of Stacking [START_REF] Wu | Ensemble learning for crowd flows prediction on campus[END_REF], without a meta-classifier because we combine the models at the evaluation phase, and a form of Bagging, because we perform an aggregation of models without applying Bootstrap sampling. Here the samples are the folds already obtained following the cross validation that we did for our previous work [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF]. Our split is stratified, which means that each fold maintains the classes distribution of the original dataset. We do not do semi-supervised learning combined with the use of Ensemble methods, as in Gong et al. [START_REF] Gong | Learning human pose in crowd[END_REF], because our dataset does not suffer from missing annotations. All the video clips are weakly annotated here.

Zhou et al. [START_REF] Zhou | Ensembling neural networks: many could be better than all[END_REF] argue that it is not useful to put a very large number of models in an Ensemble. Choosing a small number of models that are already yielding good results is enough to yield better results when they end up gathering into an Ensemble. Therefore, we decided to split the dataset into 5 folds, as it was already done in Bendali-Braham et al. [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF], which allows each Ensemble to be equipped with 4 single models that extract different knowledge from the Crowd-11 dataset.

III. ENSEMBLE CLASSIFICATION

In the light of the approaches mentioned in Section II, and based on the definition proposed by Zhou [START_REF] Zhou | Ensemble learning[END_REF], Ensemble learning consists in training, or evaluating, a set of statistical models, whether they are of similar nature or not, trained in similar conditions or not.

First, we set up Ensemble approaches made up of homogeneous models that possess the same pre-training conditions. Afterwards, we propose to create global Ensemble approaches which mix heterogeneous models of different architectures and that have various pre-training conditions.

A. Creation of homogeneous models ensembles

In a previous work on the Crowd-11 dataset, we showed that models from the 2S-I3D network perform better than the C3D network and the Inflated 3D Nets (I3D) [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF]. However, the 2S-I3D results peak approximately at 68% accuracy. These results were confirmed by a 5-fold cross-validation.

In this work, we split the dataset into 5 folds, and we train, for each possible combination of these folds, a model from one of the following architectures:

• The 2S-I3D architecture [START_REF] Carreira | Quo vadis, action recognition? a new model and the kinetics dataset[END_REF],

• The I3D architecture [START_REF] Carreira | Quo vadis, action recognition? a new model and the kinetics dataset[END_REF].

• The C3D architecture [START_REF] Tran | Learning spatiotemporal features with 3d convolutional networks[END_REF].

• The Resnet 3D (R3D) architecture [START_REF] Hara | Learning spatio-temporal features with 3d residual networks for action recognition[END_REF]. In each combination, 3 folds are used for training, 1 fold for validation, and the last fold for testing. By selecting a test fold each time, we can produce 4 models by combining the remaining folds. During the evaluation on a test fold, the decisions of the 4 models are summed to deduce the decision.

With this procedure, we can make up to 20 different combinations of the training, validation, and test sets. We mention in the following paragraphs the composition of each of the used architectures.

The I3D architecture is composed of a basis of two 3D convolution layers. Each of these layers is supported by a Batch Normalization and is followed by a 3D MaxPooling operation. These 2 layers are followed by 9 Inception modules whose internal composition changes slightly from one module to another. The last Inception module is connected to a 3D AveragePooling whose outputs are sent to a SoftMax function for the classification task.

The 2S-I3D architecture is made up of two streams. Each branch reproduces the architecture of the I3D network. One of the two streams extracts features from an RGB video clip, and the other one extracts features from an optical flow version of the video clip. The outputs of these two streams are connected to the Softmax classification function.

The C3D architecture, proposed by Tran et al. [START_REF] Tran | Learning spatiotemporal features with 3d convolutional networks[END_REF], consists of 5 layers of 3D convolutions, followed by two layers of FCN whose outputs are sent to a Softmax classification function.

The R3D architecture, proposed by Hara et al. [START_REF] Hara | Learning spatio-temporal features with 3d residual networks for action recognition[END_REF], is made up of several residual blocks. Each residual block is composed of two 3D convolution layers. We choose the version with 34 hidden layers because of its good performance which is demonstrated by Hara et al. [START_REF] Hara | Learning spatio-temporal features with 3d residual networks for action recognition[END_REF] on the Sports-1m action recognition dataset. This version of the R3D architecture consists of a first layer of 3D convolutions followed by 16 residual blocks, and ends with an FCN layer before the Softmax classification function.

1) Constitution of training, validation and test sets: The version of the Crowd-11 dataset available to us corresponds to the version that we worked on in our previous project [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF]. This version of Crowd-11 consists of 1641 scenes that can be considered as contextual boundaries between video clips. These scenes are split into 5769 video clips which are labeled into 11 classes by Dupont et al. [START_REF] Dupont | Crowd-11: A dataset for fine grained crowd behaviour analysis[END_REF]. These 11 classes are: Gas Free, Gas Jammed, Laminar Flow, Turbulent Flow, Crossing Flows, Merging Flows, Diverging Flow, Static Calm, Static Agitated, Interacting Crowd, No Crowd.

In order to split the dataset into folds to apply crossvalidation, we divided the dataset starting from the scenes in our previous work [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF]. To ensure that the division of the dataset is stratified, we apply the algorithms, presented below:

1) The first algorithm splits the scenes of the dataset into multiple folds. In our case, we split it into 5 folds.

2) The second algorithm is by the first algorithm to update the score of each fold w.r.t its distribution.

These two algorithms are used to constitute the folds that will be used for the creation of the training, validation, and test sets, as illustrated in Figure 1 2) Global ensembles of heterogeneous models: We create global Ensembles of models having either different architectures, e.g. fine-tuned 2S-I3D Ensembles coupled with C3D Ensembles trained from scratch, or different training conditions, e.g. fine-tuned I3D Ensembles and I3D Ensembles trained from scratch, or various Ensembles accumulating the two differences, for example C3D Ensembles trained from scratch, fine-tuned I3D Ensembles, and 2S-I3D Ensembles trained from scratch. We compose global Ensembles from the Ensembles of homogeneous models compared in Section IV-A. Next, we evaluate on the test set all the possible combinations from these Ensembles of homogeneous models. Equation (1) computes the number of combinations without repetition that can give rise to global Ensembles.

nb combinations = K i=2 C(K, i) (1)
Where K represents the maximum size of a combination constituting a global Ensemble. C(K, i) represents the function that computes a combination without repetition where the number of choices is set to i. i is the length of the tuple which represents the number of homogeneous Ensembles combined into a global Ensemble of heterogeneous models. As we already evaluate the Ensembles of homogeneous models in Section IV-A, i starts from 2 which is considered as the minimum size of a combination1 .

IV. EXPERIMENTS

In this section, we detail the different types of experiments that we have conducted:

• We compared the performance of Ensembles made up of fine-tuned models to other Ensembles made up of models which did not benefit from pre-training; • We compared Ensembles that were trained on augmented data from the Crowd-11 dataset and compared them to Ensembles that did not benefit from augmented data. In this work, the chosen hyperparameters for the training process correspond to the hyperparameters that we used in Bendali-Braham et al. [START_REF] Bendali-Braham | Transfer learning for the classification of video-recorded crowd movements[END_REF]. The video clips of the Crowd-11 dataset last approximately for ≈ 5 seconds. For I3D and 2S-I3D architectures, 20 frames are selected from a video clip. These frames are found at regular intervals all along the clip. The size of each image is fixed to 224 × 224 pixels. For R3D and C3D architectures, 16 frames are selected from a video clip and the size of each frame is fixed to 112 × 112 pixels. For the 2S-I3D models, the optical flow version of each clip is obtained via the TV-L1 algorithm [START_REF] Zach | A duality based approach for realtime tv-l 1 optical flow[END_REF]. During the comparison of the Ensembles that benefited from data augmentation or not, we will substitute the Farnebäck optical flow extraction algorithm for the TVL1 algorithm because the former is quicker than the latter. However, we also undertake experiments to check if the use of the Farnebäck algorithm does not significantly reduce the performance of the 2S-I3D models. The Farnebäck algorithm is theoretically known to be less accurate than the TVL1 algorithm.

A. Comparison of Ensembles of models with homogeneous architectures

We want to verify whether an Ensemble of fine-tuned models perform better than an Ensemble trained from scratch. Fine-tuned 2S-I3D and I3D models were pre-trained on the ImageNet [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF] and the Kinetics [START_REF] Kay | The kinetics human action video dataset[END_REF] datasets. The fine-tuned C3D models were pre-trained on the Sports-1m dataset [START_REF] Karpathy | Large-scale video classification with convolutional neural networks[END_REF]. Furthermore, we train from scratch models from the 2S-I3D, I3D, C3D, and R3D architectures.

In all these situations, we create 5 major training/validation and testing contexts, where we fix a test sample upstream, and we vary the validation samples (that must be different from the test sample) downstream. Under these circumstances, for each Ensemble, three samples are selected for training. These latter are different from the test sample and the selected validation sample. At the end, for the 5 test folds, 20 individual models are trained from scratch or fine-tuned. Each group of 4 singles models, whether trained from scratch or fine-tuned, constitutes an Ensemble of models at the evaluation phase.

The prediction results of these models are evaluated in terms of accuracy in Table I. Overall results demonstrate that the constitution of Ensemble models increases the average accuracy for the classification task.

B. Data augmentation

By augmenting the video data of Crowd-11, we aim to evaluate the impact of data augmentation on the classification of video-recorded crowd movements. Data augmentation was only applied on the training sets of the Crowd-11 dataset. To augment data, we applied the following operations: random crop, salt and pepper effect, video flip. In our experiments, we compare two data augmentation strategies: II, we find that using the Farnebäck extraction algorithm slightly reduces the performance of 2S-I3D models. The average score goes from 69.02% to 68.41%, which is an insignificant reduction that can be neglected for the rest of our experiments. On-the-fly data augmentation does not improve the results. Conversely, this type of data augmentation worsens the performance of the 2S-I3D models. The pre-computed data augmentation greatly improves the results of fine-tuned 2S-I3D models whose second branch feeds on flow video clips obtained via the Farnebäck algorithm. These Ensembles even progress, on average, by 1 point of accuracy from 68.41% to 69.81% thus beating the Ensembles of fine-tuned 2S-I3D models whose second streams were fed by the pre-computed flow clips extracted using the TVL1 algorithm. This good performance is offset by the 5 times greater training time required by models benefiting from the precomputed data augmentation.

C. Evaluation of global sets of models with heterogeneous architectures

We create global Ensembles by combining the models of homogeneous Ensembles whose results are illustrated in Table III. These 8 Ensembles can participate into 247 combinations containing at least 2 heterogeneous Ensembles.

The combination that achieves the best results combines the fine-tuned 2S-I3D Ensembles that have benefited from pre-computed data augmentation, with the fine-tuned 2S-I3D Ensembles that have not benefited from data augmentation, fine-tuned C3D Ensembles, and I3D Ensembles trained from scratch. The results of this combination are shown in Table IV. We find that this combination improves the overall performance by 1.5% in terms of accuracy.

Table IV shows that the Ensembles complete each other for the classification task. In this case, the least performing Ensembles do not curb the performance of the best global Ensemble. The main result of this work is that Ensemble learning improves the performances of the classification task and combining multiple Ensembles in a global Ensemble can make the results even better. Despite these good results, it is worth noting that the computation time for the decision making of the Ensembles is obviously higher than for the individual models.

V. CONCLUSION AND PERSPECTIVES

Ensemble learning improves the accuracy of the video classification task. In this work, we verified this improvement for video-recorded crowd movements illustrated by the Crowd-11 dataset. First, we found that the Ensemble of fine-tuned 2S-I3D models improve the results of their single models which were trained on different training samples. Here, thanks to Ensemble methods, the average accuracy of single models increases from 66.98% to 69.02%.

Afterwards, the comparison between the Ensembles of homogeneous models that have different training pre-conditions and deriving from different architectures shows the superiority of the pre-computed data augmentation strategy for the Ensembles of 2S-I3D models. This augmentation strategy is the appropriate regularization method that allows these Ensembles to generalize well. The pre-computed data augmentation strategy helped the 2S-I3D Ensembles to increase their accuracy from 69.02% to 69.81%.

Finally, the best global Ensemble of heterogeneous models increases the classification accuracy from 69.81% to 71.32%.

Currently, our models do not realize real-time predictions. As a remedy, we intend to focus in our future work transferring the knowledge learned by our Ensembles to lighter and quicker models, by applying knowledge distillation; a concept coined by Hinton et al. [START_REF] Hinton | Distilling the knowledge in a neural network[END_REF].

Fig. 1 .

 1 Fig. 1. Illustration of the constitution of the train, validation and test sets from different combinations of folds of the Crowd-11 dataset

•

  A fixed pre-computed data augmentation before the beginning of the training session. The augmentation methods are chosen randomly, alone or combined with other augmentation methods. As a result of data augmentation, the size of the training set was multiplied by 3. In each epoch of the training session, the model explores the augmented data as well as the non-augmented data. • An on-the-fly data augmentation that is renewed at each epoch of the training session. In this context, at each training epoch, data augmentation has a probability of 75% to occur. The epochs are repeated 4 times to allow the on-the-fly data augmentation to be similar in terms of quantity and variety to that of the pre-computed augmentated data. a) Discussion of the results of data augmentation: As a result of our experiments illustrated in Table

TABLE I COMPARISON

 I BETWEEN THE RESULTS OBTAINED BY THE PRE-TRAINED 2S-I3D ENSEMBLES FINE-TUNED ON CROWD-11 AND THEIR SINGLE MODELS

	Test sample involved	0	1	2	3	4	µ	σ
	Validation sample: accuracy per associated individual model	1: 67.86 0: 66.55 0: 66.04 0: 63.09 0: 69.73 2: 69.08 2: 66.55 1: 67.28 1: 66.26 1: 70.60 3: 67.33 3: 66.29 3: 68.52 2: 63.52 2: 69.30 4: 69.86 4: 65.44 4: 66.66 4: 62.15 3: 67.39
	Accuracies standard
	deviations of the models sharing the same test	0.99	0.45	0.91	1.52	1.17
	sample
	Accuracies mean of
	the models sharing	68.53	66.21	67.13	63.76	69.26	66.98 1.92
	the same test sample
	Accuracy per ensemble	70.48	67.57	68.61	66.43	72.00	69.02 1.99

TABLE II PERFORMANCE

 II COMPARISON BETWEEN ENSEMBLES OF 2S-I3D MODELS THAT WERE FINE-TUNED OR WERE NOT FINE-TUNED ON AUGMENTED DATA

	Involved test sample	0	1	2	3	4	µ	σ
	2S-I3D Ensembles FarneBack (Flow) Non Augmented accuracies	70.56 66.55	69.93	64.21	70.78	68.41	2.59
	2S-I3D Ensembles FarneBack (Flow) Augmented Precomputed accuracies 71.00 69.10	71.88	65.58	71.47	69.81	2.31
	2S-I3D FarneBack (Flow) Augmented On The Fly accuracies	68.47 67.23	69.23	64.21	69.30	67.69	1.89

TABLE IV COMPARISON

 IV BETWEEN THE BEST COMBINATION GIVING RISE TO A GLOBAL SET AND THE SET MODELS CONSTITUTING IT

	Test sample	0	1	2	3	4	µ	σ
	(1) C3D pretrained	61.13	60.59	61.00	58.13	61.56 60.48	1.21
	(2) I3D scratch	54.93	55.91	58.53	53.85	58.86 56.42	1.97
	(3) 2S I3D pretrained (TVL1) w/o DA	70.48	67.57	68.61	66.43	72.00 69.02	1.99
	(4) 2S I3D pretrained (Farnebäck) w DA 71.00	69.10	71.88	65.58	71.47 69.81	2.31
	Global ensemble (1) + (2) + (3) + (4)	72.05	70.04	73.20	66.35	74.95 71.32	2.95

The source code of this project is available here: https://github.com/ MounirB/Crowded-scenes-Ensemble-classification

ACKNOWLEDGMENTS

The authors would like to thank NVIDIA Corporation for the GPU Grant and the Mésocentre of Strasbourg for providing access to the GPU cluster. This work was supported by the ANR OPMoPS project (grant ANR-16-SEBM-0004) funded by the French National Research Agency.