Sebastian Basterrech
email: sebastian.basterrech@vsb.cz

Gerardo Rubino
email: gerardo.rubino@inria.fr

Evolutionary Echo State Network: evolving reservoirs in the Fourier space

Keywords: Recurrent Neural Networks, Echo State Networks, Genetic Algorithms, Reservoir Computing, Fourier Transform

The Echo State Network (ESN) is a class of Recurrent Neural Network with a large number of hidden-hidden weights (in the so-called reservoir). Canonical ESN and its variations have recently received significant attention due to their remarkable success in the modeling of non-linear dynamical systems. The reservoir is randomly connected with fixed weights that don't change in the learning process. Only the weights from reservoir to output are trained. Since the reservoir is fixed during the training procedure, we may wonder if the computational power of the recurrent structure is fully harnessed.

In this article, we propose a new computational model of the ESN type, that represents the reservoir weights in the Fourier space and performs a fine-tuning of these weights applying genetic algorithms in the frequency domain. The main interest is that this procedure will work in a much smaller space compared to the classical ESN, thus providing a dimensionality reduction transformation of the initial method. The proposed technique allows us to exploit the benefits of the large recurrent structure avoiding the training problems of gradient-based method. We provide a detailed experimental study that demonstrates the good performances of our approach with well-known chaotic systems and real-world data.

I. INTRODUCTION

A Neural Network (NN) with a recurrent topology (RNN) is a powerful tool for modeling time series and solving machine learning problems with sequential data. Recurrences (circuits) are relevant for the computational power of neural architectures; thanks to them, hidden states are able to store more information from input streams. Numerous works explain the distinct advantages presented by RNNs [START_REF] Martens | Learning Recurrent Neural Networks with Hessian-Free Optimization[END_REF]- [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF]: good theoretical properties, biologically plausibility, success in several real-world applications. However, it is still hard to properly define optimal network architectures in RNNs [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF]. The difficulties are even higher when the network is large and the input signal contains long-term data dependencies [START_REF] Martens | Learning Recurrent Neural Networks with Hessian-Free Optimization[END_REF]. In summary, the principal well-identified problems are: not guaranteed convergence, long training times, and difficulties in storing information in internal states for long time periods (long-range memory) [START_REF] Martens | Learning Recurrent Neural Networks with Hessian-Free Optimization[END_REF], [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF], [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF].

A class of single-layer RNNs called Reservoir Computing (RC) has attracted considerable attention during the last years due to its remarkable success in the modeling of complex non-linear dynamical systems [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF]. An RC model relies on the assumption that supervised adaptation of all hidden-hidden weights is not necessary [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF]. A single-layer recurrent network (so-called reservoir) is driven by an input signal, and projects the input information into a high-dimensional space. Reservoir weights are randomly This work was supported by the GACR-Czech Science Foundation project no. 21-33574K "Lifelong Machine Learning on Data Streams". initialized and scaled for satisfying technical conditions of dynamic stability [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF]. A distinctive characteristic of an RC model is that the training algorithm only adjusts the readout structure, in general recurrence-free weights that map the reservoir projection to the output space. Nowadays, learning algorithms for adapting the recurrence-free parameters are fast and robust.

The research community has been actively studying methods for creating efficient reservoir models [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF]. Specific cyclic topologies were studied in [START_REF] Rodan | Simple deterministically constructed cycle reservoirs with regular jumps[END_REF], [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF], and pre-training techniques with self-organized strategies were analyzed in [START_REF] Lukosevicius | Reservoir Computing and Self-Organized Neural Hierarchies[END_REF], [START_REF] Basterrech | Self-organizing Maps and Scale-invariant Maps in Echo State Networks[END_REF]. In [START_REF] Schrauwen | Improving Reservoirs using Intrinsic Plasticity[END_REF], [START_REF] Steil | Online reservoir adaptation by intrinsic plasticity for backpropagation-decorrelation and echo state learning[END_REF], reservoir units with parametric activation functions were trained for maximizing the entropy using the Intrinsic Plasticity (IP) rule. Other approaches based on Evolutionary Algorithms (EAs) for optimizing hyper-parameters of the reservoir structure, such as number of reservoir units, scaling parameters, connectivity properties (density and spectral radius of the reservoir weights matrix), have also been introduced in the literature [START_REF] Ferreira | Comparing evolutionary methods for reservoir computing pre-training[END_REF], [START_REF] Ferreira | An Approach to Reservoir Computing Design and Training[END_REF]. EAs have also been used for optimizing interconnected reservoirs [START_REF] Dale | Neuroevolution of hierarchical reservoir computers[END_REF], [START_REF] Ma | DeePr-ESN: A deep projection-encoding Echo State Network[END_REF], where GAs are applied for searching hyper-parameters of hierarchical reservoirs. Besides, Bayesian optimization for finding hyper-parameters of particular reservoirs was also analyzed in [START_REF] Maat | Efficient Optimization of Echo State Networks for Time Series Datasets[END_REF], [START_REF] Trierweiler Ribeiro | Bayesian optimized echo state network applied to short-term load forecasting[END_REF]. Furthermore, a relatively new prominent research in the RC area consists of models with interconnected reservoirs. A tree structure where the tree-nodes satisfies the reservoir characteristics was developed in [START_REF] Gallicchio | Tree Echo State Networks[END_REF]. The weight connections among the reservoirs are set in order of creating a global contractive mapping. Furthermore, layered/hierarchical architectures composed of multi-layered reservoirs have also been studied in [START_REF] Dale | Neuroevolution of hierarchical reservoir computers[END_REF], [START_REF] Ma | DeePr-ESN: A deep projection-encoding Echo State Network[END_REF], [START_REF]Deep echo state network (deepesn): A brief survey[END_REF]- [START_REF] Gallicchio | Deep reservoir computing: A critical experimental analysis[END_REF].

Motivations. Defining good hyper-parameters of ESNs is still often made by experience using a grid search strategy, partially brute force, and in general requires numerous trails [START_REF] Rodan | Simple deterministically constructed cycle reservoirs with regular jumps[END_REF]. Another drawback is that different reservoir matrices with the same hyper-parameters may produce substantially different results [START_REF] Schrauwen | Improving Reservoirs using Intrinsic Plasticity[END_REF]. In addition, since the reservoir is fixed during the training procedure, the computational power of the model is not fully harnessed [START_REF] Martens | Learning Recurrent Neural Networks with Hessian-Free Optimization[END_REF], [START_REF] Rodan | Simple deterministically constructed cycle reservoirs with regular jumps[END_REF].

Contributions.

In this article, we introduce the EVOlutionary Echo State Network (EvoESN), a novel computational model rooted in RC concepts. We combine ideas regarding dimensionality reduction using the Fourier transform for searching weights of RNNs (fully connected architectures) [START_REF] Koutnik | Evolving neural networks in compressed weight space[END_REF], optimization of network weights using an evolutionary search (such as in EVOLINO [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF]), and the potential of ESNs for fast and robust temporal learning [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF].

Given a reservoir with a fixed pattern of connectivity (as in standard ESNs), we apply an indirect encoding procedure for representing the reservoir weights into the Fourier space (specifically, by means of the Discrete Cosine Transform, DCT) [START_REF] Koutnik | Evolving neural networks in compressed weight space[END_REF], [START_REF] Koutník | Evolving large-scale neural networks for vision-based reinforcement learning[END_REF]. Then, we perform a fine-tuning of the reservoir weights applying genetic algorithms over the DCT coefficients. It is important to note that we are not looking for a good reservoir directly working with the weights, that becomes changing the RC paradigm and leaving the fast learning advantages of an RC model. Instead, we solve an optimization problem in the frequency domain. Using this Fourier representation drastically reduces the size of the searching space, and makes it possible to efficiently look for good reservoir weights. Our framework is a strong integration between learning in the RC paradigm and using EAs to improve the reservoir (to find better weights) following a global iterative process moving at each iteration between the weights space and the frequency space using the DCT transform and its inverse.

We demonstrate the utility of EvoESN by experiments on two well-known chaotic systems (the Mackey-Glass system and the Lorenz systems [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]) and a real-world problem (Sunspot series [START_REF] Silso World | The international sunspot number[END_REF]). Preliminary results of part of the content of this paper were briefly presented in the 2-page poster [START_REF] Basterrech | Evolving Reservoir Weights in the Frequency Domain[END_REF].

The rest of the article is organized as follows. In the next section we present the ESN framework. Our main contributions are the object of Section III, where we discuss the reservoir encoding and the evolutionary search approach. Experimental results are provided in Section IV. The article ends in Section V with a discussion about the presented work, its limitations and possible future research avenues.

II. PRELIMINARIES A. The Echo State Network framework

The large recognition of ESNs comes from their record-breaking success in predicting observations through the Mackey-Glass attractor [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. The most common application of ESNs is time series prediction, where the operation is divided into two phases. An initial training phase where the model is externally driven by an input signal and the parameters are trained used supervised learning, followed by an autonomous phase (so-called free-run prediction) where the input signal is provided by the model itself using previous predicted values [START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF].

Let us consider an ESN with K input neurons, N > K hidden neurons and L output neurons. The dynamics of the model is described in terms of a discrete-time system with an input signal u(t) ∈ R K and a hidden state x(t) ∈ R N . The recurrent state is computed according to the following recurrence [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF]:

x(t) = f W in u(t) + W r x(t -1) + W fb y(t -1) + ε(t) , ( 1
)
where t is the discrete time index, W in is a N×K matrix with input-to-reservoir weights, the N×N matrix W r has the reservoir internal weights, W fb is a weight N×L matrix with the feedback connections going from the output neurons to the reservoir, and ε(t) is a noise vector with dimension N . Function f (•) is the activation function, usually an hyperbolic function, applied element-wise. Matrix W fb and vector ε are optional [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. The model's output y(t) ∈ R L is computed by a readout layer defined as:

y(t) = g(W out [u(t); x(t)]), (2)
where W out is the N×(K+L) reservoir-to-output weight matrix, g(•) is a coordinate-wise function (usually an hyperbolic tangent or, in some cases, the identity) and [•; •] denotes vector concatenation. The readout weights are trained using a classical tool such as Tikhonov regularization [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF], [START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF]. In the case of an ESN with feedback connections (non-null W fb ), during the training procedure the feedback values are teacher-forced: instead of feeding back the reservoir with the predicted value y(t), the real target is used as the feedback value [START_REF] Reinhart | Regularization and stability in reservoir networks with output feedback[END_REF]. In exploitation mode, the model's output y(t) is fed back to the reservoir. Since the reservoir weights are frozen during learning, the initial reservoir design may have an impact on the model's performance. The most relevant hyper-parameters to be chosen are: reservoir size, input scaling factor, activation functions and reservoir spectral radius. Most research effort has gone into understanding the impact of those parameters. For more information, especially about the fundamental property of the reservoir dynamics (so-called Echo State Property (ESP)), we suggest [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF], [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF], [START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF], [START_REF] Lukoševičius | A Practical Guide to Applying Echo State Networks[END_REF], [START_REF] Basterrech | Empirical analysis of the necessary and sufficient conditions of the echo state property[END_REF].

B. Evolutionary Algorithms applied on RC models

There has been a rising interest in large-scale neural architecture engineering using evolutionary algorithms [START_REF] Elsken | Neural Architecture Search: A Survey[END_REF]. Several EAs have been developed to optimize NNs' architectures [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Rawal | Evolving Deep LSTM-based Memory networks using an Information Maximization Objective[END_REF], [START_REF] Stanley | Designing neural networks through neuroevolution[END_REF]. The RC field has also taken advantage of the recent advances in evolutionary computation. EAs have been employed for finding the hyper-parameters of an RC model, such as: sparsity, input scaling factor, injected noise, spectral radius, leaky rate and number of reservoir units. Particle Swarm Optimization (PSO) and Genetic Algorithms (GAs) have been used as metaheuristics for selecting optimal hyperparameters in [START_REF] Ferreira | Comparing evolutionary methods for reservoir computing pre-training[END_REF], [START_REF] Ferreira | An Approach to Reservoir Computing Design and Training[END_REF], [START_REF] Sergio | PSO for Reservoir Computing Optimization[END_REF]. Global parameters of an extended ESN (DRESN) were optimized using GAs for solving medical signal prognosis [START_REF] Zhong | Genetic algorithm optimized double-reservoir echo state network for multi-regime time series prediction[END_REF]. EAs have also been applied for evolving the topology of some RC models. An evolutionary pre-training procedure based on the extreme learning approach has been studied in [START_REF] Basterrech | An Experimental Analysis of the Echo State Network Initialization Using the Particle Swarm Optimization[END_REF], where the authors applied PSO for adjusting a subset of the reservoir weights. Recently, NEAT and HyperNEAT algorithms were also applied for evolving the pattern of connectivity of reservoirs in [START_REF] Chatzidimitriou | A neat way for evolving echo state networks[END_REF], [START_REF] Matzner | Neuroevolution on the edge of chaos[END_REF]. In [START_REF] Matzner | Neuroevolution on the edge of chaos[END_REF] the evolution of the reservoir topology was made considering a fitness function that combines the memory capacities of the model and its accuracy. HyperNEAT algorithm over a full recurrent topology has a significant computational cost [START_REF] Matzner | Neuroevolution on the edge of chaos[END_REF]. GAs have also been used for optimizing the hyper-parameters of hierarchical reservoirs and for finding connection parameters between reservoir structures [START_REF] Dale | Neuroevolution of hierarchical reservoir computers[END_REF], [START_REF] Ma | DeePr-ESN: A deep projection-encoding Echo State Network[END_REF]. Furthermore, a microbial GA was used for optimizing multilayered architectures of interconnected reservoirs [START_REF] Dale | Neuroevolution of hierarchical reservoir computers[END_REF]. It was also employed for finding good architectures of single ESNs with leaky integrator units. Recently, a hybrid approach using swarm optimization and local search was applied for evolving hierarchical ESNs [START_REF] Long | Evolving Deep Echo State Networks for Intelligent Fault Diagnosis[END_REF].

Differences between EvoESN and previous works.

There are also other works that tune the reservoir weights using evolutionary techniques. We identify two relevant differences between our approach and the proposed in [START_REF] Matzner | Neuroevolution on the edge of chaos[END_REF]. EvoESN performs evolution over a fixed reservoir topology and the fine-tuning is done only over the weights, while in [START_REF] Matzner | Neuroevolution on the edge of chaos[END_REF] the evolution is applied in both structure and connection weights of the reservoirs using HyperNEAT. In addition, we perform the optimization of the weights in the Fourier space applying a considerable reduction of the searching space. As already mentioned above, there are also other works where was applied EA over the reservoir weights, e.g. microbial GA was used on a single ESN [START_REF] Dale | Neuroevolution of hierarchical reservoir computers[END_REF] and PSO in [START_REF] Basterrech | An Experimental Analysis of the Echo State Network Initialization Using the Particle Swarm Optimization[END_REF]. The difference between EvoESN and these previous works is the encoding technique, in the mentioned works the evolution was applied directly over the weight connections.

III. EVOESN: EVOLUTIONARY ECHO STATE NETWORK A. Dimensionality reduction of reservoir weights applying a Fourier-type transform

Genetic encoding is a fundamental aspect of EAs, with special interest when algorithms are applied for optimizing large complex structures. A direct encoding (that is, a one-to-one mapping where every weight is represented in the gene) is clearly unsuitable for large-scale networks [START_REF] Koutník | Evolving large-scale neural networks for vision-based reinforcement learning[END_REF]. Furthermore, it has been shown that naive crossover operations present deficiencies when used for combining recurrent networks [START_REF] Lehman | Safe mutations for deep and recurrent neural networks through output gradients[END_REF]. In [START_REF] Koutnik | Evolving neural networks in compressed weight space[END_REF], an indirect encoding procedure for NNs was introduced, where the network weights are represented in the frequency domain applying a Fourier transform, more precisely, the type-II DCT of the weight matrix. The evolutionary search is conducted in a reduced space composed by the coefficients of the low frequencies only [START_REF] Koutník | Evolving large-scale neural networks for vision-based reinforcement learning[END_REF]. This leads to a significant dimensionality reduction. The approach has been successfully applied over rather small fully connected RNNs, on three benchmark problems of the area of control theory: pole balancing, ball throwing and octopus arm [START_REF] Koutnik | Evolving neural networks in compressed weight space[END_REF]. The authors have also applied this procedure on computer vision problems and in reinforcement learning [START_REF] Koutník | Evolving large-scale neural networks for vision-based reinforcement learning[END_REF]. In [START_REF] Butcher | Reservoir Computing and Extreme Learning Machines for Non-linear Time-series Data Analysis[END_REF], the authors studied the trade-off between non-linear mappings and memory capacity of ESNs. The authors applied FFTs to estimate the deviation from linearity in the reservoir activations. Furthermore, the state transition structure of the reservoir was analyzed in the frequency domain in [START_REF] Aceituno | Tailoring Echo State Networks for Optimal Learning[END_REF].

In our case, we also exploit the well-known benefits of Fourier transform for compressing data information. This efficient indirect encoding allows us to represent large reservoir matrices in a reduced space using only the low frequency coefficients of the transform. As a consequence, instead of having to tune hundreds of thousands or millions of reservoir weights, we will look for hundreds or thousands of Fourier coefficients, which implies a huge gain in complexity.

Discrete Cosine Transforms (DCTs)

. Let s = (s 0 , s 1 , . . . , s J-1 ) be a sequence of real numbers having size J. In the sequel, when we say DCT we mean the normalized DCT transform of type-II. The DCT S of s (we write S = DCT(s)) is also a sequence of J reals S = (S 0 , S 1 , . . . , S J-1 ), defined as follows:

S 0 = 1 √ J J-1 j=0 s j ,
and for = 1, 2, . . . , J -1,

S = 2 J J-1 j=0 s j cos π (2j + 1) 2J .
The inverse transform, going from S to s, denoted IDCT, is given by the following expression: for = 0, 1, . . . , J -1,

s = S 0 √ J + 2 J J-1 j=1 S j cos πj(2 + 1) 2J .
This IDCT transform is also called normalized DCT of type-III. So, S = DCT(s) and s = IDCT(S). DCT satisfies some properties that make the method suitable for encoding large recurrent structures and performing GAs. Observe that both mappings are linear, so they are continuous. Also, observe that the fact that one is inverse of the other means that they define one-to-one mappings from the space of vectors having size J to itself.

The DCT transform has a nice "energy compaction" property, informally meaning that if we keep in S = DCT(s) just a first few terms (called "low frequency" terms) and replace the rest by 0s, obtaining a new vector S , then, reconstructing the initial vector from S using the IDCT has a small error. That is, IDCT(S ) ≈ s. This characteristic of the DCT is strongly used in signal compression, and is a key element of our approach in this paper.

B. Evolving reservoirs in the Fourier space

Recall that in the canonical ESN, we select a subset of weights in the reservoir to be randomly initialized and that the rest is set to 0. We denoted by N the number of neurons in the reservoir; let M be the number of weights randomly chosen, M < N 2 (for instance, M ≈ 0.2N 2 [START_REF] Lukoševičius | Reservoir Computing Approaches to Recurrent Neural Network Training[END_REF], [START_REF] Lukoševičius | A Practical Guide to Applying Echo State Networks[END_REF]).

Let us denote by p the positions of the sampled weights in the reservoir matrix W r , a vector having size M . In our approach, we will keep fixed the vector p but we will modify the corresponding weights during the procedure. For that reason, we call these weights the "unfrozen" ones, and we call the remaining N 2 -M weights the frozen-to-zero ones. After choosing an arbitrary order on those unfrozen weights, the kth coordinate in p corresponds to some element (i, j) of W r , that is, p k = (i, j). Define now vector ν having dimension M , containing the value of those weights, that is, ν k = W r p k = W r i,j . The proposed method moves the vector ν to the frequency domain using the DCT, and there, a GA will produce new vectors, always in the frequency space. Call P this population of candidates. Each vector in P (denoted by α) will be transported to the network weights using the IDCT, occupying the fixed positions given in p. The ESN training procedure (e.g. ridge regression over the readout structure) will then be run, and the resulting network evaluated using a validation data set. The training can be done using the techniques already mentioned in Section II-A. This will allow us to select the best individuals of P. If its quality is not good enough, the best elements in P are used to start a new generation of candidates through GAs. The selection procedure used to find the best element in P must move each C-dimensional vector α to the weights space, in order to evaluate the corresponding ESN. For this purpose, we first extend α with zeros, to obtain a vector with size M . This procedure is called padding, which is necessary if C < M because the input and output of an inverse DCT have equal dimensions. The dimension C of the frequency space is a control parameter. The "energy compaction" property of DCTs makes that using some C value not too small, those coefficients will have enough information about the corresponding weights vectors, in spite of having much less components.

Figure 1 visualizes the global EvoESN procedure. Summarizing, the padding operation is performed if necessary (if C < M ), then sent to the weight space as a new feasible solution through the inverse DCT. The ESN is evaluated, after training, and its error is used in the evolutionary operations. Aiming to increase the readability, we describe the proposed algorithm into two pseudo-codes. The evolutionary search is presented in pseudo-code 1, and the computation of the fitness value of each feasible solution by the evolutionary algorithm is given in pseudo-code 2.

IV. EXPERIMENTS

To measure the computational power of the presented approach, we have selected three well-known benchmark problems: Mackey-Glass system, Lorenz system and monthly sunspot index data. Due to their interesting properties both chaotic systems have been largely studied in the RC community [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF], [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF], [START_REF] Butcher | Reservoir Computing and Extreme Learning Machines for Non-linear Time-series Data Analysis[END_REF], [START_REF] Jaeger | Controlling recurrent neural networks by conceptors[END_REF]- [START_REF] Gallicchio | Architectural and Markovian factors of Echo State Networks[END_REF]. In addition, monthly Sunspot data has information of the magnetic field in the Sun's outer regions. It also has been analyzed using RC models in at least [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF], [START_REF] Ma | DeePr-ESN: A deep projection-encoding Echo State Network[END_REF], [START_REF] Sun | TE-ESN: Time Encoding Echo State Network for Prediction Based on Irregularly Sampled Time Series Data[END_REF]. An open source monthly sunspot series is available at [START_REF] Silso World | The international sunspot number[END_REF].

A. Benchmark problems

The Mackey-Glass system (MGS). As far as we know, ESN is still the reference RNN for modeling the MGS chaotic time series [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF], defined by the equation:

∂y ∂t = αy(t -τ ) 1 + y(t -τ ) β -γy(t),
where the parameters are usually α = 0.2, β = 10, γ = 0.1. For these parameters' values, the system has a chaotic attractor if the delay τ satisfies τ > 16.8 [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. We study the system with a delay τ = 17, the most commonly used value. For this problem we replicate the experimental setting used in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. The first 1000 time steps were discarded to washout the initial transients. The training length was 3000 time steps. The remaining 2084 time steps were used for testing. From those, the first 2000 steps were used in a teacher-forced mode, then the network was in exploitation mode.

Lorenz attractor. The sequence data set is based on the following equations:

∂x ∂t = σ(y -x), ∂y ∂t = rx -y -xz, ∂z ∂t = xy -bz.
The parameters used are r = 28, b = 8/3, σ = 10, and the step size is 0.01. We replicated the preprocessing procedure for sequence x and the experimental setting described in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. The sequence was rescaled by a factor of 0.01, and the results are presented with data scaled back. There was a washout of 1000 time steps, and the training data size had 6000 time steps. For testing, we used 1000 points in a teacher-forced mode and 600 time steps in exploitation mode.

Monthly sunspot series. We studied the sunspot index from Jan. 1749 to Dec. 2021, a total of 3276 data points, available at [START_REF] Silso World | The international sunspot number[END_REF]. In this problem, we replicated the experimental setting investigated in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF]. With a minor difference that our sequence is slightly larger than the series studied in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF]. Initial washout of 100 samples, 1600 points for training, 500 points of validation and rest of the sequence for testing. // Initialization of I individuals // in GA population 1 Random initialization of P = {α (1) , . . . , α (I) }; 2 Compute fitness for each α (i) using Algorithm 2;

// Apply evolutionary search 3 while convergence criterion is not satisfied do

B. Hyper-parameter selection and experimental setup

The selection of the most important hyper-parameters of the proposed EvoESN model was made following a grid search strategy. In the case of the ESNs we also used suggestions from the literature that are described below according to each specific benchmark problem. For all the benchmark problems we applied a grid search in the interval [0.1, 1.4] for finding a good spectral radius value. Note that we consider spectral radius also larger than 1. As a consequence, the ESP is not guaranteed and the reservoir can exhibit unstable dynamics. This strategy was also applied in [START_REF] Schrauwen | Improving Reservoirs using Intrinsic Plasticity[END_REF], [START_REF] Ferreira | An Approach to Reservoir Computing Design and Training[END_REF], [START_REF] Butcher | Reservoir Computing and Extreme Learning Machines for Non-linear Time-series Data Analysis[END_REF]. In the MGS and Lorenz problems some parameters, such as injected noise, input bias, activation functions, and scaling feedback connections, were directly taken from the seminal article [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF].

MGS problem. We adopted the network configuration employed in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. We designed a reservoir with 1000 units and the hyperbolic tangent activation function. The three weight matrices W in , W r and W fb were sampled from the uniform distribution on (-1, 1). According to our results, a good candidate for the spectral radius of the reservoir matrix is 0.8, a result coherent with [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. Left figure in Fig. 2 shows the relationship between spectral radius and reservoir sparsity. The curves were produced with the average of testing errors computed over 30 independent simulations. Last, the following parameters were taken from [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]: injected uniform noise was scaled by 10 -10 , and input bias had value 0.2.

Lorenz problem. For this task we also followed the hyper-parameters selection studied in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. Input scaling and data normalization was done following [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. We used a reservoir with 600 units, with the hyperbolic tangent activation function. Right figure in Fig. 2 shows the testing error for different values of reservoir spectral radius and reservoir sparsity. Curves were created averaging the results of 30 independent simulations. We observed that the best spectral radius candidates were close to 1. The reservoir matrix had 20% of non-zero weights. Then, we set up a spectral radius to 0.97 (the same value was selected in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]). The input and reservoir matrices were initialized with the uniform distribution over (-1, 1). Uniform injected noise had values ±10 -7 [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. The feedback connections were sampled following an uniform distribution in (-4, 4) (this particular range was suggested in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]).

Sunspot series problem. We adopted the setting used in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF], and we also analyzed some parameters using grid search. An important difference in this problem with the previous task, is that the goal is to predict only one single time-step ahead (as it was studied in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF]). In previous tasks we applied free run prediction over large testing time windows. For this task we studied only networks without feedback connections. Left figure in Fig. 3 shows the relationship between spectral radius and sparsity of the reservoir. Besides, for this problem we also analyzed leaky units in the range [0.01 : 0.05 : 0.5]. For the regularization parameter we evaluated the values 10 -i with i ∈ {3, 4, 5, 6, 7, 8}, and selected 10 -5 . The selected input scaling was 0.1 after evaluation over the grid {0.01, 0.02, 0.03, 0.05, 0.07, 0.09, 0.1, 0.15, 0.2, 0.3, 0.5}.

Computational aspects. The IDCT was chosen with the orthogonal norm using the fftpack Python package. We conducted GA as EA, and it was implemented using the DEAP package [START_REF] Fortin | DEAP: Evolutionary algorithms made easy[END_REF]. We analyzed the impact of C (number of coefficients in the frequency-domain vector) in a size range between 25 and 750. Besides, different numbers of generations were evaluated: 70, 150, 300, 600. After carrying out multiple runs of the GA, we selected the hyper-parameters in order to have relatively fast convergence and to diminish the probability of being trapped in a local minima. Here we present the results for the following GA parameters. We used a Gaussian mutation with parameter mutpb=0.15, two-point crossover, probability of mating two individuals of 0.5, number of individuals participating in each tournament = 3, and population size = 20. The readout was computed using a ridge regression with a regularization factor of 10 -9 (MGS task) and 10 -6 (Lorenz task) [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF], and 10 -5 (Sunspot task).

C. Results analysis

Accuracy evaluation. It is common to see the NRMSE metric on experiments conducted with MGS and Lorenz data [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF], [START_REF] Butcher | Reservoir Computing and Extreme Learning Machines for Non-linear Time-series Data Analysis[END_REF]. Therefore, for comparison purposes we also used NRMSE for these two tasks, and for the sunspot series problem we use NMSE as it was applied in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF]. We denote by NMRSE(H) the error computed in the testing phase using the residuals over the whole horizon H (by residual we mean the difference between target and predicted values). We also present another metric denoted by NRMSE H following the notation in [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF], [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. NRMSE H is computed in the testing phase after H autonomous steps. The error is computed only using the residual at the H time step. NRMSE 84 has become the standard metric for evaluating the MGS prediction [START_REF] Schmidhuber | Training Recurrent Networks by Evolino[END_REF]. I presents a non-exhaustive overview of the reported state-of-the-art results. Previous studies have analyzed the problems using different error metrics. Then, the last column in the table shows the information about the error measure. Furthermore, we also implemented the canonical ESN, the ESN with leaky integrator units [START_REF] Jaeger | Optimization and applications of Echo State Networks with leaky-integrator neurons[END_REF], and one variation of ESN with feedback connections introduced in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. In [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF], the authors describe a learning approach for correcting the bias introduced by the teacher-forcing schema in the case of feedback connections. Jaeger et al. obtained the best performance for MGS prediction using this particular learning schema introduced as refined learning method. Let ESN with FB (I) denote the ESN with feedback connections using the teacher-forcing schema, and ESN with FB (II) denote the ESN with feedback connections using the refined learning method.

Accuracy comparison. Table

Results. We partially investigated the impact of the parameter C. Some results are shown in the right figure in Fig. 3 where the prediction error, in a logarithmic scale, is plotted against time (generation number), for different values of C. The fitness function was computed using 300 points of generalization on 12 independent prediction tasks. The curves show the evolution of the average of NRMSE across the 12 independent experiments. They also show that the accuracy is increasing with C and with time. See that if C ≥ 200, for instance, then after generation 50, the performance of the method becomes roughly independent of C. The results are consistent, the curves behave as expected, and globally, the procedure shows its efficiency once the frequency-domain vectors are large enough. Figure 4 shows results over Lorenz data with C = 150. There is a group of 6 figures. The top-left figure has an example of Lorenz testing data. Next graphics show EvoESN predictions (of a randomly selected individual) for different generations of the GA. In a few iterations, the model is able to reproduce with significant accuracy a large segment of the Lorenz attractor. Figure 5 shows the absolute error for different C values on 600 time units in the Lorenz system task. For the three values of C the results outperform the ones reported in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF] using different variations of ESNs. Note that the red curve (C = 300) has a better accuracy than the other two curves (C = 50 and C = 150) until the time horizon t ≈ 300. For longer time horizons than t = 300 both the red and blue curves seem to have similar behavior. The explanation probably lies in that the fitness function used in the GA algorithm was computed using 200 points of generalization, i.e. the reservoir weights are optimized for predicting the next 200 points. After that threshold limit the models behaves similarly (independently of C). Further possible improvements will be looked for by studying different fitness functions in the evolutionary search.

Table II summarizes the results of experiments over MGS. In the case of EvoESN, we present the results for 500 coefficients, and 70 and 300 generations in the GA algorithm. A total of 50 randomly created ESNs were evaluated. We show the average of the error produced by independent runs. EvoESN reaches much lower accuracy than the obtained with ESNs using our implementation. EvoESN also reaches the results reported in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. However, even though the variance is very low, we cannot say that the difference is significant with the error reported in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. We expect to do further studies using statistical analysis and a full evaluation of the GA parameters, in order to check if the difference is statistically significant.

Table III shows the results for the Lorenz problem. The NRMSE 84 error has the same order of magnitude as obtained by Jaeger et al. using a feedback connection and a refined learning method [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. It also shows the absolute error over a 600 time horizon, that was studied in [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]. Table IV presents the test errors for the Sunspot series problem. EvoESN outperforms the results reported in [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF], and also the results obtained by our implementations of ESN with leaky-integrators. We ran 50 independent simulations with the same ESN architecture, and the standard deviation of error across simulation has been really low (≈ 1 × 10 -6 independently of the leaky rate. Standard deviation values for EvoESN experiments are shown in the table. Regarding complexity, EvoESN has a higher cost than standard ESNs. Note that the GA performs, at each generation, training of the readouts using linear regression, computation of reservoir spectral radius and of inverse DCT. Other GA operations such as crossover, selection and mutation are less costly than the mentioned above. However, it seems that a few iterations are enough for achieving much lower errors with respect to ESNs. The optimization is performed in a frequency space whose dimension is much smaller than that of the nonnull reservoir weights. This leads to a significant dimensionality reduction, in turn producing a significant computer time reduction, while keeping the quality of the optimization procedure. These are consequences of the properties of Fourier-like transforms. Since we use an indirect encoding the dimension of the searching space is basically independent of the reservoir size, an additional advantage of the proposed technique. Another advantage is that the fine-tuning of the reservoir weights is made using the fitness function of the evolutionary search, which can be different than the one used for training the readouts. In addition, by applying GA we avoid the difficulties associated with the training of RNNs. The quality of these first experiments clearly show the potential of the approach, and the first work to be done in the close future is to push the procedure towards better performances by a careful analysis of both weight space encoding and evolutionary search procedures.

Figure 1 :

 1 Figure 1: High-level visualization of the EvoESN model.

Algorithm 1 :

 1 Evolutionary search for a fine-tuning the reservoir weights in the Fourier space.

4 Selection7 2 :

 42 Compute fitness for each α (i) using Algorithm 2; 8 end while Algorithm Evaluation of a frequency-domain C-vector α. 1 if C<M then 2 Padding operation for α; 3 end if // Mapping from frequency domain // to weight domain 4 ν = IDCT(α); 5 Assign ν values to W r according to positions in p; 6 Train readout weights; 7 Return ESN evaluation error of the given individual;

Figure 2 :

 2 Figure 2: Impact of spectral radius versus reservoir sparsity.

Figure 3 :

 3 Figure 3: Impact of hyper-parameters. Left figure shows the relationship between spectral radius and sparsity in the Sunspot series problem. Right figure shows the impact of the number of coefficients C in the MGS problem.

Figure 4 :

 4 Figure 4: Prediction of the Lorenz attractor: convergence of the proposed method according to the GA generations.

Table I :

 I Accuracy comparison with other RC methods.

	ProblemReference		Model	Error	Metric
		H. Jaeger et al. [25] ESN-FB	≈ 10 -4.2	NRMSE 84
		Schmidhüber	et	Evolino	≈ 1.9 ×	NRMSE 84
	MGS	al. [23]			10 -3
		J. Steil [51]		BPDC	0.318	NMSE
		A.	Ferreira	et	RCDESIGN	0.00051	NRMSE
		al. [12]
		C. Gallicchio et	RDESN	1.112	×	MSE
		al. [47]			10 -8
	Lorenz Jaeger et al. [25]		ESN-FB	≈ -4.27	log 10 (NRMSE 84 )
	Sunspot	Rodan et al. [6] Rodan et al. [6]		ESN DLR, DLRB, SCR ≈ 0.1039 ≈ 0.1042	NMSE NMSE

 [START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF]Table II


Soumeya Kaada
email: soumeya.kaada@nokia.com

Marie Line

Alberi Morel

Gerardo Rubino
email: gerardo.rubino@inria.fr

Sofiene Jelassi
email: sofiene.jelassi@irisa.fr

Resilience analysis and quantification method for 5G-Radio Access Networks

Keywords: Resilience, Coverage, CTMCs, Quantification, Outage prediction, 5G RAN, RSRP

A 5G Radio Access Network (5G-RAN) can be disturbed or shutdown due to a variety of failures, in spite of advanced optimization techniques and self-healing methods. Recently, operators started to take an interest in improving the resilience of communication networks with adaptive compensation techniques to mitigate outageinduced performance degradation. However, for the sake of effective and efficient resilience management, it is vital to be able to measure current and prospective resiliency levels of a given 5G-RAN using relevant and explicit metrics. Thus, the characterisation of resilience goes over a thorough analysis of 5G-RAN performance indicators followed by a rigorous quantification of current and future levels of resilience. In this work, we perform an analysis and a quantification of 5G-RAN resilience using a coverage indicator. It is known as a main performance indicator for network planners and operators, coverage is a necessary prerequisite to ensure a certain level of Quality of Service.

For that, we model the network coverage using Continuous Time Markov Chains (CTMCs) where coverage status is characterized with multiple states defined with Reference Signal Received Power (RSRP) signal. The proposed Markov model is analytically studied allowing to perform quantitative analysis, predict coverage outage and provide resilience quantification. Using our model, we conduct numerical analysis of several usage scenarios and propose a resilience framework to show the usability of our proposed approach.

I. INTRODUCTION

5G and beyond networks are expected to provide higher quality of services and to ensure more reliable communications especially in the context of Ultra-Reliable Low-Latency Communication (URLLC) systems or Ultra Dense Networks (UDNs). Given the complexity of the radio environment, multiple intelligent management techniques are deployed in 5G-RAN to ensure the radio performance requirements, like Physical layer Radio Resource Management (RRM) schemes, interference and mobility management algorithms, Coordinated Multipoint (COMP), etc, as well as multiple 5G enabling technologies like Massive Multiple-Input Multiple-output (MIMO) [START_REF] Adedoyin | Combination of ultra-dense networks and other 5g enabling technologies: A survey[END_REF]. Moreover, Self-Organizing Networks (SON) concept has been considered to offer 4G and 5G-RAN automatic processes of network management (optimization, configuration, diagnosing and healing). In turn, these Self-Management techniques help engineers to reduce operationally expensive manual oversight and to repair anomalies, degradation or failures like e.g. Coverage and Capacity Optimization (CCO) function. These management mechanisms play a pivotal role to provide an adequate Quality of Service/Experience (QoS, QoE) for end users. However, multiple risk factors can still disturb the radio environment (or RAN) resulting in failure occurrences and poor QoS such as insufficient capacity at BS, fading effects, large variation in antenna performance, energy consumption overload, etc. When the outage occurs, the main ways for engineers to deal with service blackouts are either performing manual repairs or activating self-healing methods [START_REF] Firdaus | Sleeping cell analysis in lte network with self-healing approach[END_REF], which are known to be time consuming. In fact, according to Regulatory Authority for Electronic Communications, Posts and Press Distribution (ARCEP), an average repair time of 46 hours was recorded in 2021 for 85% of quickly detected anomalies [START_REF]Autorité de régulation des communications électroniques et des postes (ARCEP)[END_REF]. Thus, current efforts in building reliable and optimised Base Stations (BSs), have to go further to protect RAN against disruptive events and service interruption and enhance the network resilience.

Generally speaking, Network Resiliency is thought of as the ability of a network to absorb the adverse impacts, adapt to and quickly bounce back from disruptions to continue keeping an acceptable level of service from users' perspective. Traditional resilience solutions used by operators, consisted on deploying costly redundant solutions to address temporally cell outage, like replacing under-performing network equipment with redundant ones. Recently, researchers showed a major interest in developing adaptive mechanisms that enhance resilience. Based on the virtualized infrastructures and the dynamic reconfiguration, prior works proposed to activate Cell Outage Compensation (COC) techniques, as a function of SON, to mitigate outage induced performance in case of failure detection [START_REF] Amirijoo | Cell outage compensation in lte networks: Algorithms and performance assessment[END_REF]. However, these few developed mechanisms are limited when the outage is hard or impossible to detect. Furthermore, network failures face a lack of management and of anticipation strategy [START_REF] Mottahedi | Resilience analysis: A formulation to model risk factors on complex system resilience[END_REF]. Besides, most prior works in RAN resilience enhancement propose resilient mechanisms without quantifying the improvements observed with resilience metrics [START_REF] Mottahedi | Resilience analysis: A formulation to model risk factors on complex system resilience[END_REF]. This makes the majority of resilience solutions in the RAN domain purely qualitative. Thus one of the major challenges addressed in this paper is to provide an analysis and quantification method of resilience in the radio domain to help current or future RAN in risk management. We assume that in the future, service or network providers or operators will need to define the network resilience requirements in Service Level Agreement (SLA) contracts based on target attributes. This implies a clear vision of 5G network resilience and how to manage it quantitatively [START_REF] Mottahedi | Resilience analysis: A formulation to model risk factors on complex system resilience[END_REF]. The analysis is proposed in the perspective to be part of a resilient management system which its main goal is to supervise the infrastructure and act according to observations. This concept is called proactive resilience [START_REF] Jackson | Proactive and Reactive Resilience: A Comparison of Perspectives[END_REF]. This approach is based on collecting and performing real-time data analytics of system abnormal events to anticipate network blackouts before specifying how to handle and plan for, or how to reduce effects of future disturbances.

Analysing resilience in RAN is related to the analysis of quality of services over the network. The quality of services is characterized by multiple Key Performance Indicators (KPIs) like: good coverage and capacity, low latency, high throughput, availability, etc [START_REF] Stoica | Automating diagnosis of cellular radio access network problems[END_REF]. Therefore, to obtain an accurate measurement of global service resilience, it is important to integrate all main KPIs in the definition of resilience. To limit complexity and cost of resilience investigation, we adopt a sub-optimal approach that consists of decomposing the initial complex problem into multiple simpler and easier sub-problems. Thus, as a Fig. 1: Radio Coverage scheme and architecture of 5G cellular site primary approach, we analyze resilience under the scope of coverage metrics that is the main network performance indicator. Coverage is one of the key factors that an operator considers when commercializing cellular communication networks due to its direct impact on service [START_REF]measuring mobile network performance: Coverage, quality of service and maps[END_REF]. Besides, coverage failures are hard to troubleshoot and diagnose in practice since logs cannot be retrieved when there is absolutely no coverage [START_REF] Stoica | Automating diagnosis of cellular radio access network problems[END_REF]. Therefore, an anticipation of coverage outage coupled with quantitative resilience metrics is required.

In this work, we provide a quantitative resilience analysis for 5G-RAN to predict coverage outage and provide resilience quantification. To do this, we are interested in Markov processes to model network coverage states. The multiple states are defined with Reference Signal Received Power (RSRP) signal of the Channel State Information-Reference Signal (CSI-RS). RSRP is a 3GPP standard radio signal [START_REF]Evolved universal terrestrial radio access (E-UTRAN); requirements for support of radio resource management[END_REF] and is equivalent to 4G-LTE RSRP. It is a good indicator of the signal strength received by the user and thus of the coverage status. We define and discuss several candidate metrics of resilience using the coverage performance indicator. A numerical analysis of several types of networks is conducted in order to study and discuss the obtained results coming from a variety of RAN use-cases. Finally, an overview of a RAN resilience framework is discussed to leverage the proposed modeling and resilience analysis within a 5G RAN system architecture.

The paper is structured as follows. Section II discusses related work about resilience and proactive analysis in the RAN domain. Section III contains a brief resilience overview in the RAN domain and presents the 5G radio coverage use case. In Section IV we present resilience analysis using Markov processes for modelling network coverage, and we define relevant resilience metrics for our quantitative analysis. In Section V, we discuss numerical results in different scenarios and in Section VI we detail the RAN resilience framework. Conclusions are given in Section VII.

II. RELATED WORK

The first and maybe crucial step of resilience management is known as resilience analysis [START_REF] Mottahedi | Resilience analysis: A formulation to model risk factors on complex system resilience[END_REF]. The analysis allows engineers to characterize the system behaviour via modeling techniques, and to quantify the resilience using suitable metrics [START_REF] Ramirez-Marquez | Generic metrics and quantitative approaches for system resilience as a function of time[END_REF]. In the literature, we noticed that most works seeking to enhance resilience propose reactive mechanisms without quantifying the improvements. We refer the reader to works on compensation methods of SON for coverage issues [START_REF] Amirijoo | Cell outage compensation in lte networks: Algorithms and performance assessment[END_REF][START_REF] De La Bandera | Adaptivz cell outage compensation in self-organizing networks[END_REF]. For example, [START_REF] Amirijoo | Cell outage compensation in lte networks: Algorithms and performance assessment[END_REF] proposes a COC method that mitigates outage impacts and restore coverage via controlling different radio parameters such as antenna tilt and the up-link target received power level P0. However, in the absence of an analysis and anticipation of disruptive events, these mechanisms are limited when facing failures that are difficult to detect, such as sleeping cells [START_REF] Firdaus | Sleeping cell analysis in lte network with self-healing approach[END_REF], and remain purely qualitative without resilience quantification.

To analyse resilience, it is interesting to get inspired from dependability studies1 that are well elaborated for the analysis of repairable systems than resilience discipline2 . For instance, a previous work on reliability analysis was proposed in [START_REF] Imran | Continous time markov chain based reliability analysis for future cellular networks[END_REF]. The authors used Continuous Time Markov Chains (CTMCs) to model the status of a specific BS into 3 states (optimal, sub-optimal and outage), and extract reliability metrics from the model. Through a fault prediction framework, they analyse the transient and steady state probability distributions of 3 study cases to predict network status. Based on metric values the framework enhance reliability via activation of root cause analysis to avoid outage. However, the reported model captures only high-level features of many categories of failure events affecting a BS, and neglects the complex cellular network environments. Besides, only reliability metrics are computed from the model. It has to be noted, that quantifying system resilience is more complex than what implies simple dependability attributes like reliability known metrics such as Mean-time-to-failure (MTTF) [START_REF] Avizienis | Basic concepts and taxonomy of dependable and secure computing[END_REF]. It is not a simple Boolean function or a binary view of the system, it is always a matter of degree, and for sure system resilience is typically not measurable on a single ordinal scale [START_REF] Avizienis | Basic concepts and taxonomy of dependable and secure computing[END_REF].

In RAN, there is no common and well-defined resilience concept which explains why we find many variations of the idea adapting it to each case [START_REF] Craig | Resilient control systems: Practical metrics basics for defining mission impact[END_REF][START_REF] Daniel | Structure-basedresilience metrics for serviceoriented networks[END_REF]. In [START_REF] Malanchini | Reliability in future radio access networks: from linguistic to quantitative definitions[END_REF], the authors proposed a definition of availability, resilience and reliability in the RAN domain. Resilience has been expressed as an integral of availability of hardware and radio resources for service, over a period of time. Others propose to compute survivability 3 (a sub-field of resilience) in wired networks, in terms of the loss probability or the delay distribution of non-lost packets in wired networks [START_REF] Heegaard | Network survivability modeling[END_REF].

Unlike previous work, we propose modelling using CTMCs for coverage states instead of considering at once all failures on BS [START_REF] Imran | Continous time markov chain based reliability analysis for future cellular networks[END_REF]. We consider radio data when modelling the network coverage, to reflect as much as possible the reality of the radio environment in the proactive analysis. We then adapt existing measures and define new ones to quantify resilience in terms of coverage. We think of resilience metrics as an attribute uncorrelated with other reliability or survivability metrics. To the best of our knowledge, this is the first work that addresses resilience quantification using coverage indicators extracted from CTMCs in the RAN domain. Besides, our analysis method can be useful to combine with reactive mechanisms like COC methods using our framework.

III. RAN RESILIENCE CONTEXT

The management of proactive resilience in RAN may be separated into: a first part that defines methods for proactive system analysis, and a second part that specifies the resilient mechanisms to fulfill the requirements in QoS. In this section, we give an overview of the general resilience management scheme in a 5G architecture and focus on the resilience analysis part for coverage issues. In section VI a detailed framework of the resilient management scheme is discussed.

A. 5G cellular network architecture

Consider a 5G network comprising multiple macro cells (MCs) and a large number of small BSs. Such network deploys an architecture with large small cells and handles multiple heterogeneous devices. UEs are connected to Remote Radio Heads (RRH) or Radio Units (RU), which sends information to 5G network via the Distributed Unit (DU) and then to the Control Unit (CU) (see figure 1). The management plan of the radio system is executed by the Network Management System (NMS). This module can integrate a resilience management framework. Depicted in figure 2, the framework gives an overview of how proactive resilience can be supported in 5G-RAN or in future cellular network architecture. It can be applied on demand to enhance resilience of a group of cells, a geographical area or a slice. It comprises two principle modules. A first module responsible for analysis of resilience. Using network performance indicators the analysis module allows predicting outage and quantifying the impacts. The second estimator module uses the results of the resilience analysis to perform decision making on what type of adaptive resilience mechanisms to activate. The framework works using data collected from BSs and stores all the relevant information needed for the process into a knowledge base. The framework will be further discussed in section VI to demonstrate its usability to enhance resilience of coverage aspect. The purpose of this work is to provide a method for the analysis module while the estimator module is left as a future work.

B. 5G cellular network use case: radio coverage

In this section, we discuss about radio coverage of a given group of serving cells deployed in a given geographical area. The area is considered as covered when the User Equipment (UE) is able to connect to their mobile network, establish and maintain a call for a minimum period of time, achieve a specific data transmission speed, and access different services. An UE can connect to network if it receives a signal strength from BS that is above a minimum acceptable coverage threshold. In urban areas, UE can be covered simultaneously from multiple gNodeBs to reduce coverage related-failures. But, rural areas have sparse deployments where no or rare overlapping coverage is available so coverage problems are more frequent. In both cases, coverage degradation or loss can still be observed within a geographical area. For 5G, the CSI -radio received signal RSRP is a valuable layer 1 measurement to provide information about the signal strength received by users within a service area. It is defined as a linear average over the power contributions (Watts) of a single Reference Signal (RS) resource element (in dBm) [START_REF]Evolved universal terrestrial radio access (E-UTRAN); requirements for support of radio resource management[END_REF]. The measures lie between -140 dBm and -44 dBm and represent a good metric of the radio coverage experimented by users.

In our work, we assume a coverage area in which a set of phone users are distributed uniformly. The number of users within the area is considered static. In real applications it may correspond to a constant traffic during a fixed period of time where the number of users vary slowly. Under the hypothesis of static traffic and uniform distribution of users, we split users into 3 categories according to their levels of RSRP. Figure 1 shows in the 5G cellular architecture the positions of each category users in terms of RSRP level that can be observed over a coverage area.

Let us denote by x = x(t), y = y(t) and z = z(t) the percentage of users that has Good, Fair, Poor levels of RSRP respectively at t. The traffic is considered static and the sum of all percentages is equal to 100%. We investigate 2 configurations of RSRP levels (C1 & C2). Table I shows the radio RF conditions of RSRP values in both cases 1 & 2. The first type corresponds to a network with minimal requirements in terms of coverage, whereas the second one corresponds to networks with higher requirements in terms of coverage, and thus higher RSRP values boundaries.

IV. COVERAGE BASED RESILIENCE ANALYSIS

In this section, we present the analysis of the radio system resilience using CTMCs, a powerful tool widely used in repairable systems study. It is based on stochastic processes allowing to predict network status and characterizing system behaviour.

Methodology for resilience analysis and quantification (A) Definition of a Markov model (states and transitions) to characterize the behaviour of the coverage status within a service area. (B) Prediction of outage occurrence while computing transient and steady state probability distributions. (C) Quantification and analysis of system resilience using resilient metrics based on coverage performance indicator.

A. Model Development

We propose a CTMC model to represent the different coverage states of L cells of 5G in a given geographical area. The model will serve to analyse the network resilience. It is formed by states and transitions between these states. The CTMC model relies on RSRP measurements that are collected within RAN to track network coverage evolution.

G F A P O λ d λ d λ d λ d λ d λ d λ c λ c λ d µ r µ o µ o µ o µ h Good Fine Acceptable Poor Outage
Fig. 3: Markov model with resilience solution Let X(t) be a CTMC with 5 states (see Figure 3) with state space S = {G = "Good", F = "Fine" , A = "Acceptable", P = "Poor", O = "Outage"}. These states are defined to reflect the behaviour of the radio coverage observed by engineers with a fine grained degrees. Each state is defined in terms of the tuple (x, y, z) previously described within the coverage area of a gNodeB. The model considers the radio system to be resilient when the majority of users within a service area satisfies a minimum acceptable coverage threshold specified with RSRP levels given in Table I. We set out to 5% the admissible percentage of users, which receive RSRP below the threshold associated with poor RSRP . This implies that the operator network guarantees a minimal radio coverage of 95% to be resilient. It is directly inspired from assumptions of mobile network radio planning for which the term of location probability is introduced. It describes the probability of the receiver being able to capture the signal with a power higher than the receiver sensitivity. On the contrary, a location probability in outdoor area of below to 95% is qualified as a non admissible coverage [START_REF] Mishra | Advanced Cellular Network Planning and Optimisation: 2G/2.5G/3G[END_REF].

The coverage model of L cells is made of 5 states (3 main states: G, A, O and 2 intermediate states: F, P) as follows:

1) X(t) = G for "Good" -Cells carrying traffic with best performances. A coverage area is Good if and only if the percentage of users with poor RSRP is less than 5% and the percentage of users with good RSRP is the highest one, i.e, iff x > y, x ≥ z, x > y + z and z < 5%. 2) X(t) = F for "Fine" -Cells carrying traffic with fine performances. This intermediate state is quite similar to G state where percentage of users with good RSRP is the highest one i.e, iff x > y, x ≥ z and z < 5% but x ≤ y + z. 3) X(t) = A for "Acceptable" -Cells still carrying some traffic, but with slightly lower performances.

A cell coverage is acceptable if and only if the percentage of users with poor RSRP is between 0 ≤ z < 4% and the percentage of users with fair RSRP is the highest one, i.e y > x, y ≥ z 4) X(t) = P for "Poor" -Cells carrying on some traffic, but with poor performances. This intermediate state is observed when percentage of users with poor RSRP is between 4 ≤ z < 5% and the percentage of users with fair RSRP is the highest one, i.e y > x, y ≥ z.

5) X(t) = O for "Outage" -A severely degraded coverage or with no traffic carried. A coverage outage is observed if and only if the percentage of users having poor RSRP is greater than 5%, i.e, z ≥ 5%. The classification of coverage status in 5 states allows having a fine grained observation and participates to a better anticipation of outages. The reader can verify that all states are fulfilled by the tuple (x, y, z), where the sum of of all percentages is 100%, and the coverage status of the service area can only have one state at a time without correlation. In the proposed Markov model, we classify failures as follows:

• Degradation: perturbations (e.g. interference or inaccurate antenna downtilt) that degrade good or fine coverage states to sub-optimal states A or F without causing failures. We denote the rate of their occurrence by λ d . • Severe degradation: perturbations (e.g. sleeping cell due to failures) that cause likely the cell outage, and are severe than simple degradation in terms of the resulting states (P or O). It may occur at sub-optimal states F,A,P. The corresponding transition rate is λ d . • Critical failures: major and unpredictable cell outages occurring less frequently. It might correspond to disasters, due to environmental conditions (hurricanes, ..) or human made (cable/power cut, ...) etc. The associated transition rate is λ c . We note that similar transitions rates appear between several pairs of states like λ d , because the same types of disruptive events are assumed to be more likely to happen. We categorize recovery transitions in the model into 3 types:

• Optimization mechanisms: are the currently applied coverage enhancements techniques in RAN like Access Traffic Steering, switching and Splitting (ATSSS) used in overlapping coverage with carrier aggregation or any other technique proposed in Rel-15/17 of 5G. They are responsible for increasing coverage quality from state (F, A, P) to G due to degradation. These techniques are practically not designed to face major and severe failures. The corresponding transition rate is µ o . We assume that the same types of mechanisms are applied in F, A, P, so in average, it takes the same speed µ o to switch to the state G. The absence of µ o between intermediate states is explained with the fact that optimizations are theoretically responsible to restore network status to an optimal and not sub-optimal or intermediate states.

• Healing mechanisms: restore coverage from state O to G through total healing/repairs or self-healing mechanisms. It may also involve adaptive resilient mechanisms that mitigate failures in an optimal way to not restrict resilience in its basic definition. The transition rate is µ h . • Resilience mechanisms: restore coverage from state O to A. In fact, the main and by default goal of resilience is to maintain an acceptable level of service under failures. These mechanisms are essentially implemented by adaptive resilient mechanisms (i.e. alternative systems, redundancy, rollback/rollforward functions, compensation or mitigation functions, failure isolation, ...). We denote the associated transition rate by µ r . The infinitesimal generator matrix Q of X(t) is:

Q =     -[2λ d +λc] λ d λ d 0 λc µo -[µo+λ d +λ d +λc] λ d λ d λc µo 0 -[2λ d +µo] λ d λ d µo 0 0 -[λ d +µo] λ d µ h 0 µr 0 -[µ h +µr]     (1)
To investigate the case of an infrastructure with no resilience mechanisms, the same model can be used except for µ r which is set to µ r = 0. In the following, the transient and the steady state distribution equations are deduced form the model in order to compute state probabilities, to predict the network coverage outage and status and to define resilience metrics.

B. Transient and steady state analysis

Transient state analysis: Using the generator matrix Q, the dynamic behaviour of the resilient CTMC model can be described by the Kolmogorov differential equation P (t) = P (t)Q where P (t) = (P G (t), P F (t), P A (t), P and i∈S P i (t) = 1. The transient state probability vector P (t) of the resilient model can be obtained by solving the differential equations system defined as follows:

                           P G (t) = -(2λ d + λ c )P G (t) + µ o P F (t) + µ o P A (t) + µ o P P (t) + µ h P O (t), P F (t) = λ d P G (t) -(λ d + λ d + λ c + µ o )P F (t), P A (t) = λ d P G (t) + λ d P F (t) -(2λ d + µ o )P A (t) + µ r P O (t), P P (t) = λ d P F (t) + λ d P A (t) -(λ d + µ o )P P (t), P O (t) = λ c P G (t) + λ c P F (t) + λ d P A (t) + λ d P P (t) -(µ h + µ r )P O (t). (2)
Steady state analysis: The steady state distribution reflects the behaviour of the system in the long run. It requires to solve the linear equation system πQ = 0 where π = (π G , π F , π A , π P , π O ), plus the normalization condition i∈S π i = 1. The linear system is:

             -(2λ d + λ c )π G + µ o π F + µ o π A + µ o π P + µ h π O = 0, λ d π G -(λ d + λ d + λ c + µ o )π F = 0, λ d π G + λ d π F -(2λ d + µ o )π A + µ r π O = 0, λ d π F + λ d π A -(λ d + µ o )π P = 0 λ c π G + λ c π F + λ d π A + λ d π P -(µ h + µ r )π O = 0. (3)
To obtain the transient and steady state probability vectors of a network that doesn't integrate resilient mechanisms, the differential equations system (2) and linear equations system (3) should be solved after setting µ r = 0.

C. Resilience metrics using coverage performance indicator

The problem of quantifying system resilience has been elusive [START_REF] Ramirez-Marquez | Generic metrics and quantitative approaches for system resilience as a function of time[END_REF]. There is a lack of standardization in resilience quantification especially in the RAN domain, and many researchers are coming up with their own definitions of resilience metrics. For a proactive resilience analysis, we investigate metrics that enable an estimation of future changes in the network (see part [START_REF] Adedoyin | Combination of ultra-dense networks and other 5g enabling technologies: A survey[END_REF]). In fact, we propose two stochastic metrics that quantify the resilience level. They allow to measure the distance between the estimated coverage status and the targeted minimal coverage requirements. In part [START_REF] Firdaus | Sleeping cell analysis in lte network with self-healing approach[END_REF], we propose a metric that assesses, how the network can "bounce back" to a normal state via various Resilience Mechanisms (RMs). We refer it as a score of resilience. So, It is a measurement of the improvement in term of resilience brought by RMs.

(1) Estimating the evolution of resilience in terms of coverage within a geographical zone over time dependable metrics:

-We define resilience metric, referred as M(t), as the global stochastic mean value of RSRP medians at time t: This metric considers the probability of being in each state P i (t) with i ∈ {G, F, A, P, O}, which is obtained from the analysis of the Markov chain. Also, med θ RSRP (t) where θ ∈ {Good, F air, P oor} represent the RSRP median of state Good(x), Fair(y) and Poor(z) described in Table I. According to state definition, we consider only the population of users having in majority good or fair RSRP for computing the RSRP median in G, F and A states. Instead, we consider the population of users having in majority fair or poor RSRP for computing the RSRP median in states P and O states -Based on [START_REF] Malanchini | Reliability in future radio access networks: from linguistic to quantitative definitions[END_REF], we propose another resilience metric that measures the availability of an acceptable coverage over a period T c . It's written as an integral over a duration T c around time t of RSRP above a minimal acceptable threshold denoted as γ as follows:

PR Tc (t) = 1 T c t+Tc/2 t-Tc/2 S R P(RSRP > γ) dt, (5)
T c is a sliding period considered to track coverage availability, and S R is the set of states {G, F, A, P} where most RSRP values are greater than a given threshold γ to have fair or good RSRP level.

(2) Quantifying the gain that might be brought by a network integrating a RM: -We denote the resilience value by RM Res(t)) when the network integrates RM and by N RM Res(t)) when there is no RM. We define the metric of resilient score as the ratio between the resilience level between RM and NRM cases. The score inspired from a previous work in [START_REF] Ramirez-Marquez | Generic metrics and quantitative approaches for system resilience as a function of time[END_REF] in equation [START_REF] Jackson | Proactive and Reactive Resilience: A Comparison of Perspectives[END_REF]. It is expressed as follows:

(6) Res score(t) = RM Res(t)) N RM Res(t)
This metric is general but in this paper, RM Res(t) and N RM Res(t) will be computed using equations ( 4) or ( 5) and expressed in dB.

Remark about methodology and resilience metrics 1) Once the transition rates in the model are available, the analyst can answer "what if" questions and evaluate any kind of metric, even those that will be defined "tomorrow". For instance, consider, as a performance criteria, the time spent at the set of main states {G, A}. Call τ that random variable, assuming the model starts at G. Standard Markov properties give us the queue of the distribution of τ (a kind of reliability at t metric), or its average. As an illustration, the expectation Exp(τ ) is:

Exp(τ ) = λ d + 2λ d + µ o (2λ d + λ c )(2λ d + µ o ) -λ d µ o .
The same can be done if we look at the subset of states {G, F, A}, or if we are interested in other more complex metrics. We don't pursuit this line here, for lack of room.

2) For sure, other resilience metrics can also be defined to enrich the analysis. Our work fits a first step of the whole resilience vision. We are currently in the steps of investigating resilience metrics considering only a coverage performance indicator. In future work, another category of metrics enriched with multiple QoS indicators as data transmission speed or latency will be considered. Their inclusion is more appropriate to ensure an accurate resilience measurement and to enhance risk management during decision making of RM solutions.

V. NUMERICAL RESULTS AND DISCUSSIONS

For numerical analysis, we apply our methodology on four practical study cases that depict different configurations of transient rates. Firstly, we discuss about the proposed rate values for each study case given in Table II. Then, the rate values are used to analyse transient and steady state probabilities given in equations ( 2) and (3). They help to predict outage occurrence and quantify the resilience. They are drawn in figures 4 and 5 in the case of basic and critical networks, respectively. Then, the evolution of resilience over time, given in equations ( 4) and [START_REF] Mottahedi | Resilience analysis: A formulation to model risk factors on complex system resilience[END_REF], is studied quantitatively. It is shown in figures 6 and 7(a), respectively. In figure 7(b), we analyse the score of resilience (eq. 6) that quantifies the improvements brought by RMs in figure 7.

(a) Study cases configurations: We consider 4 study cases that correspond to various real network types. We distinguish 2 types of network groups that are basic network and critical network. The first group corresponding to cases 1 and 2 represent common networks that carry basic network services with minimal levels of QoS requirements. The second group corresponding to cases 3 and 4 represent distinguished networks that carry critical services with higher requirements in terms of QoS. Critical networks are currently well optimized against degradation and apply frequently preventive dependability methods to reduce failure events. To evaluate the impact of resilience, we assume that each type has two different resilience configurations. Cases 1 and 3 don't integrate resilience schemes while cases 2 and 4 have RM enabled. Table II contains theoretical values in hours of the failure transient rates λ i and recovery values µ i for all study cases. They are set based on analysis of real RAN performance problems reported in works [START_REF]Autorité de régulation des communications électroniques et des postes (ARCEP)[END_REF][START_REF] Imran | Continous time markov chain based reliability analysis for future cellular networks[END_REF]. In real implementation, these rates are computed statistically from failure log data. II shows that the rate λ c of major failures is observed in average once every 2 months for basic network services, and around every 3 months for critical network services. In case 1, where no resilience is enabled, we set the degradation rate λ d to 1 every 10 hours and severe degradation λ d , that can likely cause the outage, every 24h. To limit degradation, optimization functions operate in background. We assume a recovery transition happens with rate 1/4h, in other words, those optimization functions take, on the average, 4h to restore coverage (and µ o = 1/5h for case 1). In case of outage occurrence, healing and repairs take roughly 2 days including the detection time and the root cause analysis either automatically or manually (that is, µ h = 1/48h). In case 2 where the network uses RMs, we should expect observing less failure occurrences and faster recovery processes than in a non resilient network. In fact, we assume it takes 2h on average to deploy an adaptive resilience solution to restore the coverage to an acceptable level (µ r = 1/2h). This would increase occupancy time at state A where optimisation is applied efficiently within 4h. It has to be noted that RM activation impacts the network through decreasing slightly rates of degradation or severe degradation. Here, we choose to consider it by decreasing λ d to 1/20h to illustrate the impact of RMs on degradation transitions. Rate µ h is also increased compared to case 1 because it integrates total repairs plus optimal resilience solutions as explained in Section IV. For the second type of networks, we assume that we will observe less failure events in critical infrastructures compared with the first type since they already include rigorous optimizations and more efficient healing functions, but still not sufficient to cope with all failures encountered, and sometimes, it is possible that they don't have adaptive resilient solutions to face quickly the failures (case 3). Therefore, we increase λ d to once every 3 days, the arrival time of degradation λ d to 20h and we assume that the optimizations take in average 4h. These values are identical for both cases 3 and 4. However, in case 4 we focus on studying the impact of RM activation on the resilient transition µ r only set to 1h and the optimal resilient solution or total healing transition µ h that takes in average 24h.

(c) Outage prediction over a time period: Figures 4 and5 represent the results of transient analysis of states probabilities P (t). The probabilities are computed over 15 days in all cases. The analysis of probabilities reveals that a basic network with no resilience (case 1) encounters a higher failure risk measured by P O (t) than in the other cases with approximately 40% after 5 days, while for critical networks with no resilience (case 3) its value is around 10%, due to the optimisations assumed to occur in this type of network. In addition, introducing resilience reduces P O to 2-5% in case 2 and to 1% in case 4, and enhances performances through increasing P G . We also notice that P A is more stable after 2 days in case 2 and doesn't decrease compared with case 1, meaning that resilience creates more chances to be in state A and increases the occupancy time as highlighted above. CTMCs are powerful tools to track short term network evolution. For instance, the evolution of P O (t) allows the operator to predict the Fig. 4: Transient analysis of cases 1 and 2 outage risk occurrence within a time limit, e.g 5 days. This is helpful to anticipate correctly the disruption according to operator's resilience policy. For example, resilience solution like adaptive RM (e.g. antenna tilt modification) can be prepared for deployment and at a suitable activation time. -80 dBm) predefined within an operator policy, to anticipate and avoid coverage outage. We observe that in the case of basic network with no resilience (case 1) the network resilience is under the threshold fixed by policy once every 4 days. In figure 6(b) we notice different values of resulting global stochastic mean of RSRP medians on same cases of critical networks using different configurations. According to this figure, we verify that the critical network with RM (case 4) is the more resilient. Similarly, we can observe in figure 7(a) the availability over a sliding period of 2 days for all cases. We notice that the availability falls after 48h to around 5%. This metric gives another vision of the resilience evolution in the future than the transient probabilities computed at instant t. Finally, figure 7(b) shows resilience score to quantify the improvements of basic and critical networks with RM. The score is computed using global stochastic mean of RSRP medians. We note an optimal improvement of RSRP with 17dB in basic networks and 7dB in critical networks.

To conclude, proactive resilience analysis can be done through Markov modeling and quantification using various metrics for short term prediction. The usage of average RSRP over the coverage zone allows focusing on the impact of failures instead of their types, and hence all types of failures are considered. The analysis quantifies how much vulnerable the basic networks with non resilient solutions are, and how the RM solutions improve network coverage performance status either in basic network services or strengthen more critical network services. However, CTMC is still limited for the long run and looking on RSRP signal only for providing indication on coverage status is not enough. Thus, other signals and KPIs may be considered (RSRQ, throughput...). Besides, accurate resilience measurement requires to consider other indicators like interference, capacity, etc, to quantify the impact of failures on multiple service aspects. In real application, the challenge is to obtain the correct transient rates from past failure logs of the observed networks. Next section provides an overview of utilisation of this method in RAN.

VI. DATA-BASED RESILIENCE FRAMEWORK

As discussed in Section III, a resilience management scheme can be deployed to enhance RAN resilience within a control loop via a data-based and automated framework located at NMS level in 5G (figure 2). In the case of coverage, figure 8 shows that the framework starts with the analysis module learning the transient rates (λ d , µ r , ..) using ML tools or statistical learning from past failure logs to build the Markov model. These estimated rates will be used to analyse the transient and steady state distribution to predict the network status and quantify resilience. Based on the output of the analysis and the resilience metrics, some scripts containing instructions can be prepared to activate adaptive RMs to restore coverage on the area, such as modifying dynamically antenna tilt or P0 of neighboring cells for Outage compensation [START_REF] Amirijoo | Cell outage compensation in lte networks: Algorithms and performance assessment[END_REF]. The estimator module is driven by the operator resilience policy, to take the right decisions on which RM to activate and at what time. The estimator can use intelligent decision algorithms like regression or reinforcement learning. Both modules run in background and in real time, and can communicate with each other. The historical data, model, transient rates, policies and scripts are stored into a database (Figure 8). Besides, a resilience manager will serve to instantiate the modules, coordinate their communication and control the process. Finally, the resilience process is evaluated with real-time metrics to validate the accuracy of the estimation process. The pipeline is executed within a life cycle.

VII. CONCLUSION

This paper provides a proactive quantitative resilience analysis for 5G-RAN. The analysis uses Markov processes for modelling network coverage using RSRP data. Based on the model, we predict coverage outage occurrence and propose resilience metrics to quantify the level of network resilience. To show the usability of our approach, we present a numerical analysis in 4 different use case scenarios. A databased resilience framework is proposed. According to numerical results, the proactive analysis can help to characterize the behaviour of the system in short or mid-term, and take anticipated decisions to improve the system by activating resilience mechanisms. These resilience mechanisms envision to improve current basic networks and strengthen the already well optimized critical networks. For future work, the proposed predictive model will be tested on real or simulated radio data to compute the transient rates. It will be also extended to consider other key performance indicators (throughput, latency, etc). Finally, machine learning methods will be considered to help decision making in resilience mechanims activation.

Fig. 2 :

 2 Fig. 2: RAN Resilience Framework

(

Fig. 5 :Fig. 6 :

 56 Fig. 5: Transient analysis of cases 3 and 4 (d) Quantitative evolution of resilience: The metrics defined in section IV-C quantify resilience level and gain to help operators or network owners getting useful and familiar resilience measurements insights.We assume that RSRP are uniformally distributed for both RSRP configurations of table I and we plot

Fig. 7 :

 7 Fig. 7: (a) Acceptable coverage availability over a sliding period T c = 2d (b) Basic and critical networks resilience score.

Fig. 8 :

 8 Fig. 8: Proactive resilience analysis and estimation steps the evolution of global stochastic mean of RSRP medians in figure 6. To observe how this metric behaves we use (a) the same RSRP configuration C1 for all study cases and (b) different configurations C1 & C2 on same study cases (we focus on case 3 and 4). The observation gives relevant information for network planners or operators on how RSRP levels decrease over time under potential failures, and what are minimal resilience conditions to ensure. These conditions may serve to be compared with a threshold (eg. (a):-100 dBm, (b):-80 dBm) predefined within an operator policy, to anticipate and avoid coverage outage. We observe that in the case of basic network with no resilience (case 1) the network resilience is under the threshold fixed by policy once every 4 days. In figure6(b) we notice different values of resulting global stochastic mean of RSRP medians on same cases of critical networks using different configurations. According to this figure, we verify that the critical network with RM (case 4) is the more resilient. Similarly, we can observe in figure7(a) the availability over a sliding period of 2 days for all cases. We notice that the availability falls after 48h to around 5%. This metric gives another vision of the resilience evolution in the future than the transient probabilities computed at instant t. Finally, figure7(b)shows resilience score to quantify the improvements of basic and critical networks with RM. The score is computed using global stochastic mean of RSRP medians. We note an optimal improvement of RSRP with 17dB in basic networks and 7dB in critical networks.To conclude, proactive resilience analysis can be done through Markov modeling and quantification using various metrics for short term prediction. The usage of average RSRP over the coverage zone allows focusing on the impact of failures instead of their types, and hence all types of failures are considered. The analysis quantifies how much vulnerable the basic networks with non resilient solutions are, and how the RM solutions improve network coverage performance status either in basic network services or strengthen more critical network services. However, CTMC is still limited for the long run and looking on RSRP signal only for providing indication on coverage status is not enough. Thus, other signals and KPIs may be considered (RSRQ, throughput...). Besides, accurate resilience measurement requires to consider other indicators like interference, capacity, etc, to quantify the impact of failures on multiple service aspects. In real application, the challenge is to obtain the correct transient rates from past failure logs of the observed networks. Next section provides an overview of utilisation of this method in RAN.

TABLE I

 I

: Radio RF condition configurations (RSRP )

  4) M(t) = med Good RSRP (t)P G (t) + med (Good+F air) RSRP (t)P F (t) + med F air RSRP (t)P A (t) + med (F air+P oor) RSRP (t)P P (t) + med P oor RSRP (t)P O (t).

TABLE II :

 II Transient rates of all study cases (b) Discussion on transient values: Table

Dependability and reliability methods are designed to avoid failures and to be mainly active before outage occurrence (e.g predictive maintenance).

Resilience methods try to adapt the network once the outage occurs. They mainly react after failure to quickly restore the service.

Survivability focus on attacks and large scale disasters[START_REF] Sterbenz | Resilience and survivability in communication networks: Strategies, principles, and survey of disciplines[END_REF]