Kevin Pasini
email: kevin.pasini@irt-systemx.fr

Mostepha Khouadjia

Allou Same

Fabrice Ganansia

Allou Samé

Latifa Oukhellou

LSTM encoder-predictor for short

Keywords: Machine learning, Time series forecasting, Representation learning, Mobility data, Transport, Train load

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

In recent years, the population growth in metropolitan areas has led to overcrowding on trains. Transport operators are working on enriching real-time passenger information systems by providing passengers with train loads in addition to train schedules. This information can allow passengers to better plan their daily trips, which can improve overall comfort and avoid overcrowding on trains.

Moreover, such forecasting can be used by public transport authorities and transport operators either to enrich public transport route planning or to improve the synchronization of train traffic with passenger flows. Transport operators will increasingly need to assess and predict network passenger load to improve train regulation processes and service quality levels.

The development of smart technologies and the rapid growth in data storage abilities have increased the availability of massive transport data, such as passenger affluence, train load, real-time train schedules and so forth. This increase in data contributes to leveraging the developments in data mining and machine learning approaches for processing such spatio-temporal data to extract valuable information with the aim of providing better services to passengers or to match the transport supply with the demand. This paper addresses the forecasting of train load at a railway station considering a historical dataset that includes two data sources: train load data and automatic vehicle location. The latter source contains all information related to the train operation (delay, time of arrival/departure of vehicles and so on). Most of the prediction problems in this domain address the prediction of passenger affluence at an aggregated level (per 15 minutes or 30 minutes time horizon) [START_REF] Toqué | Forecasting dynamic public transport origin-destination matrices with long-short term memory recurrent neural networks[END_REF][2] [START_REF] Ziat | Spatio-temporal neural networks for space-time series forecasting and relations discovery[END_REF]. In contrast to these studies, we focus here on the prediction at the non-aggregated level, taking into account real-time train schedules. This induces variability in the time step of the time series that we should predict. Furthermore, the prediction model has to take contextual factors impacting train load into account, such as calendar information (day, time, holiday and so forth) and train operation.

We address this prediction task as a multi-step short-term forecasting problem on irregularly structured times series influenced by several contextual factors. We work at the station level for each passage of train, which involves a temporal variability, making difficult the application of techniques that usually exploit structure regularity of time series. To handle these specificities, we rely on the abstraction capabilities of neural networks linked to the concept of representation learning [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]. The underlying idea is to build a mobility representation of our known influential factors. The model takes the form of an encoder-predictor neural network architecture associated with representation learning on contextual factors. It aims to predict the train load of the next trains of a station from the values of the last trains and all the contextual features characterizing these trains.

This paper is organized as follows. Section 2 is devoted to related work conducted on prediction models dedicated to public transport demand forecasting and on some main works on deep neural networks. In Section 3, we detail the case study that we consider and our dataset. Section 4 presents the proposed methodology formulated in a general framework. The evaluation of the proposed model is then conducted in Section 5 through different experiments, aiming to compare the performance of the proposed model to those provided by four baseline models on the one hand and to illustrate the learning representation space of the network on the other hand. Section 6 concludes the paper.

Related work

Numerous studies have addressed the mining of large-scale mobility data for exploration, clustering or prediction purposes. Depending on the available data, on the scale of analysis and on the targeted goal, different types of methodologies can be distinguished. For example, the authors in [START_REF] Ceapa | Avoiding the crowds: Understanding tube station congestion patterns from trip data[END_REF] propose predicting crowding levels from automated fare collection data by using simple techniques based on historic aggregates. Models proposed in [START_REF] Tsai | Neural network based temporal feature models for short-term railway passenger demand forecasting[END_REF] take the form of an neural network architecture using feature engineering capturing daily and monthly trends to perform daily passenger demand forecasting. The study conducted in [START_REF] Lathia | Mining mobility data to minimise travellers' spending on public transport[END_REF] explores the viability of building a fare-recommendation system for public transport to avoid incorrect fares. For this purpose, the authors propose a two-step approach: predicting future travel habits and then matching travel habits to fares.

Extensive research on predictive models in the transport domain has been conducted. Most of the works address the prediction of passenger flows. In [START_REF] Ke | Short-term forecasting of passenger demand under on-demand ride services: A spatio-temporal deep learning approach[END_REF], the authors propose a deep neural network model to forecast passenger demand related to an on-demand ride service station. The combination of convolutional layers with LSTM layers allows taking spatial, temporal and exogenous dependencies into account. An LSTM recurrent neural network was proposed in [START_REF] Toqué | Forecasting dynamic public transport origin-destination matrices with long-short term memory recurrent neural networks[END_REF] to address the short-term forecasting of passenger flows in a transport network. The authors in [START_REF] Ding | Predicting short-term subway ridership and prioritizing its influential factors using gradient boosting decision trees[END_REF] work on short-term subway ridership prediction by means of a gradient boosting decision tree model. On the basis of smart card data, the authors build a prediction model using both temporal features (time and calendar) and historical data related not only to subway activities but also to bus transfer activities. Following the same line of research, the authors in [START_REF] Roos | Short-term urban rail passenger flow forecasting: A dynamic bayesian network approach[END_REF] examine causal relationships between the adjacent flows on a public transport station with transport service features. The proposed methodology, which is based on a dynamic Bayesian network, allows highlighting such causalities and performing the prediction. Recently, the authors in [START_REF] Heydenrijk-Ottens | Supervised learning: Predicting passenger load in public transport[END_REF] formalized the problem of tram load passenger prediction as a classification task, where the passenger load is labelled into different classes depending on the percentage of occupied seats. Once this labelling is performed, the authors build classical machine learning classifiers, such as k-nearest-neighbours, multi-layer perceptron, grading boosting decision tree and random forest, to predict the level of crowding in the transport. Here, temporal and historical data were used as model inputs.

Advanced machine learning models, particularly deep learning models, have also recently been proposed to address forecasting problems. Recurrent neural networks [START_REF] Gers | Learning to forget: Continual prediction with lstm[END_REF] are potential tools capable of capturing the dynamics of the time series. These models have even been extended to handle regular spatio-temporal data by [START_REF] Xingjian | Convolutional lstm network: A machine learning approach for precipitation nowcasting[END_REF] with a convolutional LSTM for weather prediction. In parallel, early research work has been accomplished by [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF]. The authors propose a recurrent neural network encoder-decoder for a statistical machine translation system. This model is capable of capturing both semantic and syntactic structures of phrases. Considering the spatio-temporal dependencies in the forecasting, [START_REF] Ziat | Spatio-temporal neural networks for space-time series forecasting and relations discovery[END_REF] proposes a dynamical spatio-temporal neural network for forecasting the time series of spatial processes. The idea investigated in this model is to learn both temporal and spatial dependencies between the series to be predicted through a combined use of a latent embedding structured by the temporal dynamics of the series and a decoder mechanism to make the prediction. In [START_REF] Zhang | Deep spatio-temporal residual networks for citywide crowd flows prediction[END_REF], the authors work on forecasting the flow of crowds in all regions of a city by means of a deeplearning-based model. Historical trajectory data, weather and events are used to build the model. Residual neural networks and a parametric fusion mechanism are employed to design the forecasting model of crowd traffic.

Focusing on the transport domain, most of the aforementioned research works propose achieving the load/affluence prediction using classical classification or regression machine learning models. These models do not allow fully exploiting the sequential structure of the time series to be predicted. Moreover, they do not consider variability in time steps of the data, as it is the case here. The flexible transportation schedule, the variability in transport demand, and the contextual factors lead to complex dynamics of the series that the model should capture.

To address the structural variability in the passenger load series and influence factors, we rely on the abstraction capabilities of deep neural network models linked to the concept of representation learning [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]. The underlying idea is to learn a meaningful representation of mobility flows taking contextual factors into account. The proposed model takes the form of an RNN encoder-decoder neural network [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF] associated with the representation learning of contextual factors. It aims to predict the passenger load of the next trains at a station from measures of the last trains and all the contextual features characterizing all of these trains at the same station. Note that the experiments will be performed on a real dataset collected over a long period (one and a half years), which leads to a robust evaluation of the models. Before detailing the proposed model, the next section will describe the real dataset used in the experiments.

Data description

Our study focuses on a dataset collected from a French railway line that serves approximately fifty stations located in the northern area of suburban Paris. The railway line carries approximately 250,000 passengers daily. The dataset covers a period of 18 months from January 2015 to June 2016 on 40 stations for daytime exploitation from 5 am to 2 am of the next day. It includes both timetable information and count data of passengers boarding and alighting at each station collected by radar sensors on trains (2000000 records covering 86% of trains). These heterogeneous sources of data that have been enriched with calendar information enable us to reconstruct the passenger load on each train departing from a station.

The main goal of this study concerns the forecasting of univariate train load series for each station. To have an idea of the time series to be predicted, Figure 1 shows an example of weekday and weekend daily train passenger load profiles collected from two stations. The suburb station accounts for approximately 22000 train stops with a particularity low train frequency and few train routes that serve this station. Conversely, the inner-city station accounts for approximately 84000 train stops with an important train frequency and multiple train routes that serve the station. Figure 1 provides insights on the forecasting problem to be solved and highlights the particularities of our dataset, namely:

-A variable sampling period due to the train timetables and railway operation.

Each station has its own train frequency evolution. -A specific temporal behaviour of each time series, which was found to be linked to the spatial location of public transport stations and geographical aspects of the city (population & employment densities, leisure and so forth). -Train load series are impacted by calendar factors such as the type of day (weekday or weekend), holiday, public holiday and so on. -Train load series are also impacted by train characteristics that are closely linked to their services (multi-destination line, various train services).

In addition to these contextual factors, public transport demand and therefore train load passengers can also be impacted by events (social, cultural, sport and so forth). The time series forecasting model has to face all the temporal, spatial and exogenous factors listed above. In this paper, a prediction model will be built for each railway station. The next section details the methodology developed to achieve this prediction task.

LSTM Encoder-Predictor

In this section we formalize the application of the recurrent encoder-predictor architecture to our particular structural constraints: a train sequence with variable time steps and heterogeneous attributes. Let (y 1 ,..,y t ) denote the sequence to be predicted, where y t corresponds to a passenger train load in the same station, t referring to the trains arrival order. It is assumed that each realization y t is associated with an observation S t which includes contextual features e t and past measures m t . We have tackled the issue of the variability of the time between consecutive trains by encoding it as a contextual feature associated with specific coefficients in the model. We also use the notation y I , S I , e I and m I to designate the subsequences (y t ) t∈I , (S t ) t∈I , (e t ) t∈I , (m t ) t∈I with I ⊂ [1; T ]. Given a time window

W i = [i -k, i + k ] composed of a past horizon P i = [i -k, i[ and a future horizon F i = [i, i + k ],
the goal of our multi-step forecasting approach is to infer a realization on the horizon y Fi from information available on S Wi as shown in Figure 2. Table 1 summarizes the notation used in this article. This forecasting is particularly challenging because it asks to understand the laws behind realizations (y t ) taking into account the multiple influencing factors. The model must be able to dissociate the influencing factors on a structurally irregular sequence by exploiting the contextual attributes. Following the line of research on the RNN encoder-decoder proposed by [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF], we propose a neural network LSTM encoder-predictor for short-term multi-step prediction including a representation learning of the contextual factors.

u t = (u Pi , u Fi ) (2)
Enc : A many-to-one LSTM 'encoder' is dedicated to capture a past latent dynamic (h i ) from the past measures m Pi and the past contextual representations (u Pi ).

Enc(m Pi , u Pi ) = h i (3)
Dec : A many-to-many LSTM 'decoder' decodes recurrently latent reconstruction states r Pi of past observations from latent dynamics of the past horizon (h i ). Each latent reconstruction state is then interpreted by 'Reconst', a linear reconstruction layers that infer past observation realization. From 'Reconst' outputs we get ŷPi . These outputs are used as an intermediate objective during the training phase to facilitate capture of past latent dynamics.

Dec(h i ) = r Pi (4) Reconst(r Pi ) = t∈Pi Reconst(r t ) = ŷPi (5)
Enc and Dec form an encoder-decoder structure that synthesizing the dynamics of past observations from theirs contextual and measurement features. P red : A many-to-many LSTM 'predictor' infers latent prediction states (z Fi ) of future observations from their contextual representations (u Fi ) taking into account the latent dynamics of the past horizon (h i ). Each latent prediction state is then interpreted by 'Predict', a linear prediction layers that infer future observation realization. From 'Predict' outputs we get ŷFi which corresponds to the multi-step prediction aim. Note that since the model is designed to address variability in the time step, this makes it straightforward to remove observations of the dataset due to missing data. Moreover, once the LSTM encoder-predictor is trained, predictions can be performed on missing data in the future horizon if we are able to reconstruct contextual information.

P red(h i , u Fi ) = z Fi (6) P redict(z Fi ) = t∈Fi P redict(z t ) = ŷFi (7)

Optimization

A deep neural network is trained through end-to-end gradient back-propagation by minimizing the following loss function:

L(θ) = α p * t∈Pi ||y t -ŷt || 2 + α f * t∈Fi ||y t -ŷt || 2 , W ith θ = (θ F act , θ Enc , θ Dec , θ P red , θ Reconst , θ P redict ). ( 8
)
The first term measures the ability of the model to reconstruct the past observations from the latent past dynamics. It is an intermediate objective that facilitates the learning of the past dynamics. The second term measures the prediction ability of the model. Hyper-Parameters α p and α f are the weights of the reconstruction and prediction objectives.

For the learning phase, we realize mini-batch optimization thanks to a Nadam optimizer [START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF]. Two gradients (prediction and reconstruction) are propagated from their output layers (Predict and Reconst) to the upstream layers towards the context factory through LSTM layers. The encoder-predictor is implemented based on the TensorFlow [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF] environment and Keras [START_REF] Chollet | [END_REF] as a library and highlevel neural network API.

The parameters have been chosen empirically after several experiments based on model performance and learning convergence. Fact is composed of 3 dense layers of size [50, 100, 200] with a sigmoid activation function for a total of 27000 parameters. Enc, Dec, and Pred are 3 LSTM layers of size 200 with a sigmoid activation function for a total of 880000 parameters. Reconst and Predict are composed of 2 dense layers of size [100, 1] with a linear activation function for a total of 40000 parameters. The whole neural network has approximately 900,000 parameters.

Training is empirically realized on a batch of size 128 on several thousand iterations, which takes few hours on one standard GPU card depending on the dataset and time depth. Further work on the choice of parameters is required to improve the convergence.

Experimental results and discussion

For the experimental part, we evaluate two models for train load forecasting, each one on a given dataset respectively related to a station situated in the suburb and in the inner city. The datasets concern the period from January 2015 to June 2016, and they are split into training (year 2015 ≈ 66% ) and test sets (year 2016 ≈ 33%). Both models use several contextual information e t as long-term (calendar) features:

-Day of the year (8-dimensional): Position of the day in the year (365 possible values) encoded by cosine and sine of (2 ×4) frequencies. -Day type (8-dimensional): Position of the day in the week with an additional attribute if the day corresponds to a holiday or not.

-Minutes (8-dimensional): Minute of the day (1440 possible values) encoded by cosine and sine of (2×4) frequencies. -Train services (8-dimensional): Feature related to the train routes that serve the considered station.

Moreover, we also consider short-term features m t by considering a lag window that ranges between 1 to 6 past observations:

-Delay (1-dimensional): Difference in minutes between the real and the theoretical schedule of the train at the station. -Load (6-dimensional): Number of passenger on the train for each of the last 6 passages at the considered station. -Board (6-dimensional): Numbers of boarded passengers for of the last 6 train passages at the station. -Alight (6-dimensional): Numbers of alighted passengers for each of the last 6 train passages at the station.

Compared Forecasting Models

We compare on both suburb and inner-city stations different classical and machine learning models that serve as baselines for our LSTM EP:

-Last Value (LV): It is the simplest forecasting that consists of forwarding the last observed load on the train to the next one at the same station. -Contextual Average (CA): It consists of using of the average load of trains that are committed on the same day type and time slice. -Gradient Boosting (XGB) [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF]: A regressor model that produces a prediction model in the form of an ensemble of weak decision trees. Two models are proposed depending on the input features:

• XGB LT: Model that is trained based on long-term features e t .

• XGB ST: Model trained on both long e t and short-term m t features. -LSTM: Basic recurrent network using both types of features.

-LSTM EP: The proposed RNN trained on context representation using both types of features.

The best parameters of XGB have been selected with a grid search procedure in conjunction with cross validation for K = 3 as the number of folds.

We evaluate the performance of models on each time step by root mean square error (RMSE) and weighted absolute percentage error (WAPE) measures:

W AP E score : t ||y t -ŷt || ȳ (9)
WAPE is a derivative of the MAE that can be interpreted as the percentage of the overall error compared to the average value of the actual observation.

Forecasting Results

The evaluation of the forecasting models is conducted by making comparisons based on performance metrics. These metrics are expressed in terms of RMSE and WAPE and are given for both the training and the test phases for the suburb and inner-city stations. The five models defined in section 4.1, namely, the basic last value (LV), the contextual average (CA), the gradient boosting (XGB) short term (ST), long term (LT) and the LSTM-EP are compared. The errors obtained for both the training and test sets are given in Table 2.

The results show that advanced models (XGB, LSTM) outperform the LV and CA models. This performance improvement could mainly be explained by the fact that the XGB and LSTM models have better generalisation abilities and are able to fit more complex models than LV and CA, which simply predict by forwarding the last observed value or averaging historical data. The basic LSTM performs less in the inner-city station compared to the suburb station. This can be explained by its difficulties to deal with the heterogeneity in term of train services with that kind of station. Overall, LSTM EP leads to the best results since it is able to capture in better way the underlying dynamics of the temporal irregularity related to the heterogeneity of train services by means of its encoder-decoder component. Looking at the prediction error of the LSTM EP according to the load to predict (see Figure 5), we observe that errors are increasing according to the load. The model tends to slightly overestimate weakly loaded trains and greatly underestimate the highly loaded trains. Heavily loaded trains are rare and present contextual information similar to many less loaded trains, which explains the difficulty of the model in predicting large loads. To remedy this problem, provisioning features to distinguish these trains appears to be necessary. One could imagine indicators related to the disturbance of the network and the known presence of events near the station. As shown in Figure 6, the model makes errors of the same order of magnitude for weekdays and weekends with different difficulties. The variance of the error over a time slot is correlated to the importance of the load. On weekdays, we observe larger errors in the morning and afternoon peak hour linked to the strong variance and high load. The model makes average errors in the middle of the day and low errors in the morning and evening. On weekends, except in the morning, a relatively similar error variance is observed with a maximum at noon and in the middle of the afternoon. We also observe that the model has more trouble predicting weekend evenings than weekday evenings. When we examine the performance of the models on multi-step temporal prediction, the LSTM EP outperforms XGB and basic LSTM for the next 6 time steps (Table 3 andTable 4). These time steps correspond to the train passages at the station and range between 14 to 182 minutes for the suburb station, whereas it ranges between 2 to 61 minutes when considering the inner-city station since the train passages are more frequent. Note that these performances are obtained with a single model that simultaneously predicts the load at all the time steps for both LSTM models. The XGB LT is time-step invariant since it only considers long term features. For XGB ST, we have as many models as the considered time steps. The performance of short-term models are slightly degraded when we move forward in time, excepted for basic LSTM in the case of inner city station, where we can notice a strong degradation of its performance over time steps due to the heterogeneity of train services.The LSTM EP shows very competitive results and better robustness compared to other models for both the suburb and inner-city stations for all steps considered. This can be explained by a better understanding of contextual factors through a latent representation that helps to capture the underlying dynamics of train service at the station.

Representation learning exploration

In this section, we propose exploring the latent spaces provided by our neural network. These latent spaces correspond to the projection of contextual features and predictive state representations learned during the training phase of our train load target.

For this purpose, we first project the learned representation of the penultimate layer of our LSTM EP obtained during the training of the train load at the suburb station. Figure 7 shows the scatter plot of the dimensional reduction of the latent representation. For each sub-plot, the considered representation is projected depending on a given feature (day, hour, month, train load). Each colour block represents a common category of the handled feature. Dimensionality reduction is performed by preserving large pairwise distances between the points with the help of principal component analysis (PCA) to reduce 200-dimensional data into 50 dimensions, followed by a T-distributed stochastic neighbour embedding (tSNE) to reduce the dimensionality from 50 to 2 or 3 dimensions. As shown in Figure 7, the contextual representation (u t ) is strongly structured according to different levels depending on each contextual feature. Each point of the obtained structure reflects a train passage, and depending on the gradient of the underlying structure, we can highlight the profile of train service during weekday and weekend (see Figure 7.a), the frequency of service with a distinction between the rush and off-peak periods (see Figure 7.b), the seasonality (see Figure 7.c) and the load range (see Figure 7.d). When we reduce the predictive latent space (z t ) to 3 dimensions, we can observe the structure of the manifold as a spiralling band representing the daily evolution of the train load (Figure 8).

Conclusion

In this paper, we investigated train load forecasting on both single-and multistep time horizons. The single step aims to predict the load for the next train at the station, whereas the goal of multi-step is to predict ahead of time the load for further train passages. The forecasting problem is particularly challenging since we have variability in train services and several contextual factors that could have a major impact on the time series to be predicted. For this purpose, we proposed a model called LSTM EP based on an encoder-decoder combined with a representation learning of contextual factors. This network has the particularity of being able to learn a contextual representation from contextual features, to capture the latent past dynamics through the underlying sub-structure of the encoder-decoder, and then forecast the forthcoming dynamics with the help of the predictive layer. The representation learning is a key value for the proposed architecture; it contributes to understanding the features and representation of the underlying dataset, which is essential for selecting the best neural network architecture for this prediction task.

The obtained results have demonstrated the potential of the LSTM encoderpredictor to address the short-term prediction of train load at stations compared to the Gradient Boosting model and basic LSTM. We evaluated the performance of the proposed model on two real datasets related to suburb and inner-city stations for single and multi-step forecasting horizons of the train load. On both configurations, the LSTM EP outperforms the LSTM, XGB and baseline models by maintaining robustness in the quality of the forecasts throughout the time horizon.

Future research should investigate the exploration of the learned representation. In particular, it would be interesting to investigate the ability of the predictive latent space to characterize abnormal situations, such as disturbances and traffic anomalies.

Fig. 1 .

 1 Fig. 1. Train loads in year 2015 per hour on suburb and inner-city stations

Fig. 2 .

 2 Fig. 2. Illustration of notations used in the forecasting model Table 1. Notation and variable Notation t A time step t ∈ [1, T ] y1, ..., yT (yt) Realization series S1, ..., St (St) Observation sequences e1, ..., eT (et) Sequence of feature contextual vectors m1, ..., mT (mt) Sequence of feature measure vectors Windows Wi [i -k, i + k ] : Window associated to the ith observation Pi [i -k, i[: Past horizon of window Wi Fi [i, i + k ] : Prediction horizon of window Wi Xi (mP i , eP i , eF i ) Input features from the window Wi Latent space (see subsection 3.1) u1, ..., uT (ut) Contextual representation h1, ..., hT (ht) Latent past dynamic r1, ..., rT (rt) Latent reconstruction state z1, ..., zT (zt) Latent prediction state Model and sub-part LST M EP Neural network model F act MLP factoring contextual features Enc Recurrent encoder of past observation Dec Recurrent decoder of past observation. P red Recurrent predictor of future observation Reconst MLP to reconstruct past realizations P redict MLP to predict future realizations

3. 1 1 )Fig. 3 .

 113 Fig. 3. General architecture of the LSTM encoder-predictor network

Fig. 4 .

 4 Fig. 4. Details on the layout of the LSTMs

Fig. 5 .

 5 Fig. 5. Prediction errors depend on load class for suburb station

Fig. 6 .

 6 Fig. 6. Prediction errors depend on hour class for suburb station

Fig. 7 .Fig. 8 .

 78 Fig. 7. Latent representation (ut) according to contextual features obtained after dimensionality reduction (PCA + tSNE)

Table 2 .

 2 Model performance on both studied stations

	Model	Suburb	Inner city
		WAPE RMSE	WAPE RMSE
	LV	17.9	35.8 Train score	41.9	186.7
	CA	13.7	28.7	14.2	73.1
	XGB LT	8.4	17.2	8.3	44.75
	XGB ST	7.5	15.1	8.2	43.5
	LSTM	11.5	24.3	8.9	51.5
	LSTM EP	10.7	22.1	10.9	57.7
	LV	24.1	47.2 Test score	46.9	205.0
	CA	19.0	40.0	18.5	96.5
	XGB LT	18.8	38.9	13.4	76.0
	XGB ST	16.8	35.7	12.7	73.0
	LSTM	16.2	34.0	13.7	75.3
	LSTM EP	16.0	33.8	12.9	72.4

Table 3 .

 3 RMSE test score of the suburb station for the multi-step forecasting models

	Model	t+1	t+2	t+3	t+4	t+5	t+6
	Time interval*	14-32 29-62 44-92 59-122 75-152 90-182
	XGB LT	38.9	38.9	38.9	38.9	38.9	38.9
	XGB ST	35.7	36.6	36.7	36.7	37.6	38.1
	LSTM	34.0	34.4	34.8	35.5	36.3	36.9
	LSTM EP	33.8 34.0 34.1 34.4	34.7	34.9

Table 4 .

 4 RMSE test score on the inner-city station for the multi-step forecasting

	Model	t+1 t+2 t+3 t+4	t+5	t+6
	Time interval*	2-13 5-23 9-31 12-43 15-53 18-61
	XGB LT	76.0 76.0 76.0 76.0	76.0	76.0
	XGB ST	73.0 72.8 73.3 73.8	73.4	73.5
	LSTM	75.3 75.4 80.2 83.9	90.5	92.9
	LSTM EP	72.4 72.1 72.1 72.2 72.6 72.8

*The 5th and 90th percentiles of the time interval in the passage of trains at the time T and T+n


Romain Mathonat
email: romain.mathonat@insa-lyon.fr

Diana Nurbakova
email: diana.nurbakova@insa-lyon.fr

Jean-Francois Boulicaut
email: jfboulicaut@gmail.com

Mehdi Kaytoue
email: mehdi.kaytoue@gmail.com

SeqScout: Using a Bandit Model to Discover Interesting Subgroups in Labeled Sequences

Keywords: Pattern Mining, Subgroup Discovery, Sequences, Upper Confidence Bound

It is extremely useful to exploit labeled datasets not only to learn models but also to improve our understanding of a domain and its available targeted classes. The so-called subgroup discovery task has been considered for a long time. It concerns the discovery of patterns or descriptions, the set of supporting objects of which have interesting properties, e.g., they characterize or discriminate a given target class. Though many subgroup discovery algorithms have been proposed for transactional data, discovering subgroups within labeled sequential data and thus searching for descriptions as sequential patterns has been much less studied. In that context, exhaustive exploration strategies can not be used for real-life applications and we have to look for heuristic approaches. We propose the algorithm SeqScout to discover interesting subgroups (w.r.t. a chosen quality measure) from labeled sequences of itemsets. This is a new sampling algorithm that mines discriminant sequential patterns using a multi-armed bandit model. It is an anytime algorithm that, for a given budget, finds a collection of local optima in the search space of descriptions and thus subgroups. It requires a light configuration and it is independent from the quality measure used for pattern scoring. Furthermore, it is fairly simple to implement. We provide qualitative and quantitative experiments on several datasets to illustrate its added-value.

I. INTRODUCTION

In many data science projects we have to process labeled data and it is often valuable to discover patterns that discriminate the class values. This can be used to support various machine learning techniques the goals of which are to predict class values for unseen objects. It is also interesting per se since the language for descriptions is, by design, readable and interpretable by analysts and data owners. Among others, it can support the construction of new relevant features. The search for such patterns has been given different names like subgroup discovery, emerging pattern or contrast set mining [START_REF] Novak | Supervised descriptive rule discovery: A unifying survey of contrast set, emerging pattern and subgroup mining[END_REF]. We use hereafter the terminology of the subgroup discovery framework [START_REF] Wrobel | An algorithm for multi-relational discovery of subgroups[END_REF].

Labeled sequential data are ubiquitous and subgroup discovery can be used throughout many application domains like, for instance, text or video data analysis, industrial process supervision, DNA sequence analysis, web usage mining, video game analytics, etc. Let us consider an example scenario about the maintenance of a cloud environment. Data generated from such a system are sequences of events. Applying classification techniques helps answering to "will an event occur" (See, e.g., [START_REF] Xing | A brief survey on sequence classification[END_REF]), while applying sequential event prediction helps determining "what event will occur" (See, e.g., [START_REF] Letham | Sequential event prediction[END_REF]). Nevertheless, another need is to explain, or at least, to provide hypotheses on the why. Addressing such a descriptive analytics issue is the focus of our work. Given data sequences, labeled with classes, we aim at automatically finding discriminative patterns for these classes. Given our example, sequences of events could be labeled by breakdown presence or absence: the goal would be to compute patterns "obviously occurring with breakdowns". Such patterns provide valuable hypotheses for a better understanding of the targeted system and it can support its maintenance.

Let us now discuss research issues when considering subgroup discovery on labeled data. Given a dataset of object descriptions (e.g., objects described by discrete sequences), a sequential pattern is a generalization of a description that covers a set of objects. An unusual class distribution among the covered objects makes a pattern interesting. For example, when a dataset has a 99%-1% distribution of normalabnormal objects, a pattern covering objects among which 50% are abnormal is highly relevant (w.r.t. a quality measure like the Weighted Relative Accuracy measure WRAcc [START_REF] Lavrac | Rule evaluation measures: A unifying view[END_REF]). However, such patterns cover generally a small number of objects. They are difficult to identify with most of the algorithms that exploit an exhaustive exploration of the search space with the help of a minimal frequency constraint (See, e.g., SD-MAP [START_REF] Atzmller | SD-Map -a fast algorithm for exhaustive subgroup discovery[END_REF]). Therefore, heuristic approaches that are often based on beam search or sampling techniques are used to find only subsets of the interesting patterns. An open problem concerns the computation of non-redundant patterns, avoiding to return thousands of variations for Boolean and/or numerical patterns [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF]. Heuristic methods for sequential data have not yet attracted much attention. [START_REF] Diop | Sequential Pattern Sampling with Norm Constraints[END_REF] introduced a sampling approach that is dedicated to the frequency measure only. A promising recent proposal is [START_REF] Egho | A user parameter-free approach for mining robust sequential classification rules[END_REF]: the sampling method misère can be used for any quality measure when exploiting sequences of events. Their key idea, is to draw patterns as strict generalizations of object descriptions while a time budget enables it, and to keep a pool of the best non redundant patterns found so far. Such a generic approach has been the starting point of our research though we were looking for further quality assessment of discovered subgroups.

Our contribution provides a search space exploration for labeled sequences of itemsets (and not just items) called SeqScout. It is based on sampling guided by a multiarmed bandit model followed by a generalization step, and a phase of local optimization. We show that it gives better results than a simple adaptation of misère with the same budget when considering huge search spaces. Our algorithm has several advantages: it gives results anytime and it benefits from random search to limit redundancy of results and to increase subgroup diversity. It is rather easy to implement and can be parallelized. It is agnostic w.r.t. the used quality measure and it finds local optima only (w.r.t. the considered quality measure).

Section II discusses related work. We formally define the problem in Section III. Section IV describes SeqScout. Section V presents an empirical study on several datasets, including applications on Starcraft II data (game analytics) and abstracts of articles published in the Journal of Machine Learning and Research.

II. RELATED WORK

Sequential pattern mining [START_REF] Agrawal | Mining sequential patterns[END_REF] is a classical data mining task. It remains a challenging task due to the size of the search space. For instance, Raïssi and Pei [START_REF] Raïssi | Towards bounding sequential patterns[END_REF] have shown that the number of sequences of length k is

w k = k-1 i=0 w i |I| k-i ,
with |I| being the number of possible items. As an example, if we consider the dataset promoters [START_REF] Dua | UCI machine learning repository[END_REF], |I| = 4 and k = 57 (see Table II), the size of the search space is approximately 10 41 . Various methods have been proposed to mine interesting sequential patterns. We review them briefly and we discuss their relevancy when considering our need for discriminative patterns.

A. Enumeration-based methods

Many enumeration techniques enable to mine patterns from Boolean, sequential and graph data [START_REF] Giacometti | 20 years of pattern mining: a bibliometric survey[END_REF]. They can be adapted for the case of discriminative pattern mining. For instance, the SPADE algorithm [START_REF] Zaki | Spade: An efficient algorithm for mining frequent sequences[END_REF] has been adapted for sequence classification based on frequent patterns [START_REF] Zhou | Pattern based sequence classification[END_REF]. The main idea of such methods is to visit each candidate pattern only once while pruning large parts of the search space. Indeed, we know how to exploit formal properties (e.g., monotonicity) of many user-defined constraints (and not only the minimal frequency constraint). The quality measure is thus computed either during or after the discovery of all frequent patterns [START_REF] Bosc | A Pattern Mining Approach to Study Strategy Balance in RTS Games[END_REF]. This is inefficient for the discovery of the best discriminative patterns only. To overcome this limitation and to help pruning the search space, upper bounds on the quality measure can be used. However, they remain generally too optimistic and are specific to a particular measure (See, e.g., [START_REF] Nowozin | Discriminative subsequence mining for action classification[END_REF]). Moreover, enumeration techniques coupled to upper bounds aim at solving the problem of finding the best pattern in the search space, and not the best pattern set.

B. Heuristic methods

An interesting current trend to support pattern discovery is to avoid exhaustive search and to provide anytime sets of high quality patterns, ideally with some guarantees on their quality. Examples of such guarantees are the distance to the best solution pattern [START_REF] Belfodil | Anytime subgroup discovery in numerical domains with guarantees[END_REF] or the guarantee that the best solution can be found given a sufficient budget [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF]. Let us discuss some of the heuristic approaches that have been proposed so far.

Beam search is a widely used heuristic algorithm. It traverses the search space level-wise from the most general to the most specific patterns and it restricts each level to a subset of non redundant patterns of high quality [START_REF] Duivesteijn | Exceptional model mining[END_REF]. The greedy nature of beam search is its major drawback. Boley et al. proposed a two-step sampling approach giving the guarantee to sample patterns proportionally to different measures: frequency, squared frequency, area, or discriminativity [START_REF] Boley | Direct local pattern sampling by efficient two-step random procedures[END_REF]. However this method only works on those measures. Diop et al. proposed an approach which guarantees that the probability of sampling a sequential pattern is proportional to its frequency [START_REF] Diop | Sequential Pattern Sampling with Norm Constraints[END_REF]. It focuses on the frequency measure only. Egho et al. have proposed the measure agnostic method misère [START_REF] Egho | A user parameter-free approach for mining robust sequential classification rules[END_REF]. Given a time budget, their idea is to generate random sequential patterns covering at least one object while keeping a pool of the best patterns obtained so far. It provides a result anytime, empirically improving over time.

Sequential pattern mining can be modeled as a multiarmed bandit problem enabling an exploitation/exploration trade-off [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. Each candidate sequence is an "arm" of a bandit. For instance, Bosc et al. have developed such a game theory framework using Monte Carlo Tree Search to support subgroup discovery from labeled itemsets [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF]. They proposed an approach based on sampling, where each draw improves the knowledge about the search space. Such a drawn object guides the search to achieve an exploitation/exploration trade-off.

To the best of our knowledge, the problem of mining discriminative sequences of itemsets with sampling approaches has not been addressed yet in the literature. We introduce our SeqScout method that computes top-k non-redundant discriminative patterns, i.e., tackling a NP-hard problem [START_REF] Qin | Diversifying top-k results[END_REF]. For that purpose, we want hereafter to maximize the well known quality measure called W RAcc [START_REF] Lavrac | Rule evaluation measures: A unifying view[END_REF]. Even though we focus on it, our method is generic enough for using any quality measure, without requiring specific properties. SeqScout does not require parameter tuning, unlike Beam Search. In contrast to misère, we have the guarantee to find local optima thanks to a hill-climbing search.

III. NON-REDUNDANT SUBGROUP DISCOVERY:

PRELIMINARIES AND PROBLEM STATEMENT Let I be a set of items. Each subset

X ⊆ I is called an itemset. A sequence s = X 1 ...X n is an ordered list of n > 0 itemsets. The size of a sequence s is denoted n, and l = n i=1 |X i | is its length. A database D is a set of |D| sequences (See an example in Table I). Given a set of classes C, we denote by D c ⊆ D the set of sequences in D that are labeled by c ∈ C. Definition 1 (Subsequence). A sequence s = X 1 ...X ns is a subsequence of a sequence s = X 1 ...X n s , denoted s s , iff there exists 1 ≤ j 1 < ... < j ns ≤ n s such that X 1 ⊆ X j1 ...X ns ⊆ X jn s . In Table I, {a}{b, c} is a subsequence of s 1 and s 2 .
Definition 2 (Extent, support and frequency). The extent of a sequence s is ext

(s) = {s ∈ D | s s }. The support of a sequence s is supp(s) = |ext(s)|. Its frequency is f req(s) = supp(s)/|D|. In Table I, ext( {a}{b, c} ) = {s 1 , s 2 }. Definition 3 (set-extension). A sequence s b is a set- extension of s a = X 1 X 2 ...X n by x if ∃i, 1 ≤ i ≤ n + 1 such that s b = X 1 ...{x} i ...X n+1 .
In other words, we have inserted an itemset X i = {x} in the i th position of s a . {a}{c}{b} is a set-extension of {a}{b} .

Definition 4 (item-extension). A sequence s

b is an item- extension of s a = X 1 X 2 ...X n by x if ∃i, 1 ≤ i ≤ n such that s b = X 1 ...X i ∪ {x}, ..., X n+1 . {a, b}{b} is an item-extension of {a}{b} . Definition 5 (Reduction). A sequence s b is a reduction of s a if s a is an set-extension or item-extension of s b .
Definition 6 (Quality measure). Let S be the set of all subsequences. A quality measure ϕ is an application ϕ : S → R that maps every sequence s ∈ S with a real number to reflect Table I

: An example database D. id s ∈ D c s 1 {a}{a, b, c}{a, c}{d}{c, f } + s 2 {a, d}{c}{b, c}{a, e} + s 3 {e, f }{a, b}{d, f }{c}{b} - s 4
{e}{g}{a, b, f }{c}{c} -its interestingness (quality score). For instance, Precision, defined by P(s → c) = supp(s,Dc) supp(s,D) , is a quality measure. Definition 7 (Local Optimum). Let N (s) be the neighborhood of s, i.e., the set of all item-extensions, set-extensions and reductions of s. r is a local optimum of S w.r.t. the quality measure ϕ iff ∀r ∈ N (r ), ϕ(r ) ≥ ϕ(r).

Definition 8 (Non θ-redundant patterns). A set of patterns, (also called subgroups or subsequences here), S p ⊆ S is non θ-redundant if given θ ∈ [0; 1] and ∀s1, s2 ∈ S p , where s 1 = s 2 , we have: sim(s 1 , s 2 ) ≤ θ. In the following, we use the Jaccard index as a similarity measure:

sim(s 1 , s 2 ) = |ext(s 1 ) ∩ ext(s 2 )| |ext(s 1 ) ∪ ext(s 2 )| .
Problem Statement (Non-redundant subgroup set discovery). For a database D, an integer k, a real number θ, a similarity measure sim, and a target class c ∈ C, the non redundant subgroup discovery task consists in computing the set S p of the best non θ-redundant patterns of size |S p | ≤ k, mined w.r.t the quality measure ϕ.

IV. METHOD DESCRIPTION

SeqScout is a sampling approach that exploits generalizations of database sequences, and search for local optima w.r.t. the chosen quality measure. Fig. 1 provides a global illustration for the method 1 .

A. General overview

The main idea of the SeqScout approach is to consider each sequence of the labeled data as an arm of a multi-armed bandit when selecting the sequences for further generalization using the Upper Confidence Bound (UCB) principle (See Algorithm 1). The idea of the UCB is to give a score to each sequence that quantifies an exploitation/exploration trade-off, and to choose the sequence having the best one.

First (Lines 2-4), priority queues, π and scores, are created. π stores encountered patterns with their quality, and scores keeps in memory the list of U CB scores of each sequence of the dataset, computed by using Equation 1(See Section IV-B). data + contains the list of all sequences of the data set labeled with the target class. Indeed, taking sequences having the target class will lead to generalizations having at least one positive element. Then, the main procedure is launched as long as some computational budget is available. The best sequence w.r.t U CB is chosen (Line 9). This sequence is 'played' (Line 10), meaning that it is generalized (See Section IV-C) and its quality is computed (See Section IV-G). The created pattern is added to π (Line 11). Finally, the U CB score is updated (Line 12). As post processing steps, the top-k best non-redundant patterns are extracted from scores using the filtering step (See Section IV-D). Finally, these patterns are processed by a local optimization procedure (See Section IV-E). Moreover, SeqScout needs other modules that concern the selection of the quality measure (See Section IV-F) and the quality score computation (See Section IV-G).

B. Sequence selection

We propose to model each sequence of the dataset as an arm of a multi-armed bandit slot machine. The action of playing an arm corresponds to generalizing this sequence to obtain a pattern, and the reward then corresponds to the quality of this pattern. Following an exploitation/exploration trade-off, sequences leading to bad quality patterns will be avoided, leading to the discovery of better ones.

The problem of the multi-armed bandit is well known in the Game Theory literature [START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF]. We consider a multiarmed bandit slot machine with k arms, each arm having its own reward distribution. Our problem is then formulated as follows. Having a number N of plays, what is the best strategy to maximize the reward? The more an arm is played, the more information about its reward distribution we get. However, to what extent is it needed to exploit a promising arm, instead of trying others that could be more interesting in the long term (explore)? Auer et al. proposed a strategy called UCB1 [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. The idea is to give each arm a score, and to choose the one that maximizes it:

U CB1(i) = xi + 2ln(N ) N i , ( 1
)
where xi is the empirical mean of the i th arm, N i is the number of plays of the i th arm and N is the total number of plays. The first term encourages the exploitation π ← P riorityQueue()

3:
scores ← P riorityQueue()

4:
data + ← F ilterData()

5:
for all sequence in data + do 6:

scores ucb .add(sequence, ∞)

7:
end for 8:

while budget do scores.update(seq,

Ni * qual+qualp Ni+1

, N i + 1)

13:

end while 14:

π.add(Optimize(π))

15:

return π.topKN onRedundant() 16: end function return pattern, qual 26: end function of arms with good reward, while the second encourages the exploration of less played arms by giving less credit to the ones that have been frequently played.

Performing an exploitation/exploration trade-off for pattern mining has already been applied successfully to itemsets and numerical vectors by means of Monte Carlo Tree Search [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF]. When dealing with a huge search space, using sampling guided by such a trade-off can give good results. However, contrary to [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF], we consider the search space of extents, not intents. Exploring the extent search space guarantees to find patterns with non null support while exploring the intent search space leads towards multiple patterns with a null support. This is a crucial issue when dealing with sequences.

C. Pattern generalization

After the best sequence w.r.t. UCB1 is chosen, it is generalized, meaning that a new more general pattern is built. It enables to build a pattern with at least one positive element. Indeed, most of the patterns in the search space have a null support [START_REF] Raïssi | Towards bounding sequential patterns[END_REF]. SeqScout generalizes a sequence s in the following way. It iterates through each item within each itemset X i ∈ s, and it removes it randomly according to the following rule:

remain, if z < 0.5 remove, if z ≥ 0.5
, where z ∼ U(0, 1).

D. Filtering Step

To limit the redundancy of found patterns, a filtering process is needed. We adopt a well-described set covering principle from the literature (See, e.g., [START_REF] Bosc | Anytime discovery of a diverse set of patterns with monte carlo tree search[END_REF], [START_REF] Van Leeuwen | Diverse subgroup set discovery[END_REF]), summarized as follows. First, we take the best element, and then we remove those that are similar within our priority queue π. Then, we take the second best, and continue this procedure until the k best non-redundant elements are extracted.

E. Local optimum search

Finally, a local optimum search is launched w.r.t. Definition 7. Various strategies can be used. The first possible strategy is the Steepest Ascend Hill Climbing. It computes the neighborhood of the generalized pattern, i.e., all its itemextensions, set-extensions and reductions. Then, it selects the pattern among those of the neighborhood maximizing the quality measure. This is repeated until there is no more patterns in the neighborhood having a better quality measure. Another possible strategy is the Stochastic Hill Climbing: a neighbor is selected at random if its difference with the current one is large "enough". Notice however that it introduces a new parameter. Depending on the dataset, the branching factor can be very important. Indeed, for m items and n itemsets in the sequence, there are m(2n + 1) patterns in its neighborhood (See Theorem 1). To tackle this issue, we use First-Choice Hill Climbing [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF]. We compute the neighborhood until a better pattern is created, then we directly select it without enumerating all neighbors.

Theorem 1. For a sequence s, let n be its size, l its length, and m the number of possible items, the number of neighbors of s, denoted |N (s)|, is m(2n + 1).

Proof: The number of item-extensions is given by:

|Iext| = n i=1 |I| -|X i | = nm - n i=1 |X i | = nm -l.
We have now to sum the number of reductions, setextensions and item-extensions:

|N (s)| = l + m(n + 1) + |Iext| = m(2n + 1).

F. Quality measure selection

The choice of the quality measure ϕ is application dependent. Our approach can deal with any known measures that support class characterization, such as, e.g., the F 1 score, informedness or the Weighted Relative Accuracy [START_REF] Lavrac | Rule evaluation measures: A unifying view[END_REF]. The later, the W RAcc, is commonly used for discriminant pattern mining and subgroup discovery. It compares the proportion of positive elements (i.e., sequences labeled with the target class) to the proportion of positive elements in the whole database. Let c ∈ C be a class value and s be a sequence:

W RAcc(s, c) = f req(s) × supp(s, D c ) supp(s, D) - |D c | |D| .
It is a weighted difference between the precisions P(s → c) and P( → c). The weight is defined as f req(s) to avoid the extraction of infrequent subgroups. Indeed, finding very specific subgroups covering one positive element would result in a perfect quality value but a useless pattern. W RAcc value ranges in [-0.25, 0.25] in the case of a perfect balanced data set, i.e., containing 50% of positive elements.

We consider objective quality measures that are based on pattern support in databases (whole dataset, or restricted to a class). Using random draws makes it particularly difficult as each draw is independent: we cannot benefit from the data structures that are so efficient for classical exhaustive pattern mining algorithms do (See, e.g., [START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF]).

G. Efficient computation of quality scores

To improve the time efficiency of support computing, bitset representations have been proposed. For instance, SPAM [START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF] uses a bitset representation of a pattern when computing an item-or set-extension at the end of a sequence. In our case, we consider that an element can be inserted anywhere. Therefore, we propose a bitset representation that is independent from the insertion position. Its main idea lies in keeping all bitset representations of encountered itemsets in a hash table (memoization), and then combining them to create the representation of the desired sequence. The main idea is given in Fig. 2. Assume we are looking for the bitset representation of {ab}, {c} . Let {c} be an already encountered pattern (i.e., its representation is known) while {ab} was not. This can not be handled by the SPAM technique as a new element has to be added before a known sequence. Our algorithm will first try to find the bitset representation of {ab} . As it does not exist yet, it will be generated and added to the memoization structure. Then, position options for the next itemset are computed (Line 2 in Fig. 2). The latter is then combined with a bitset representation of {c} using bitwise AND (Line 4). The support of the generated sequence can then be computed.

V. EXPERIMENTS

We report an extensive empirical evaluation of SeqScout using multiple datasets. All the experiments were performed on a machine equipped with Intel Core i7-8750H CPU with 16GB RAM. Algorithms mentioned hereafter are implemented in Python 3.6.

Datasets: We perform the evaluation using some benchmark datasets that are well known in the pattern mining community, namely aslbu [START_REF] Papapetrou | Discovering frequent arrangements of temporal intervals[END_REF], promoters [START_REF] Dua | UCI machine learning repository[END_REF], blocks [START_REF] Dua | UCI machine learning repository[END_REF], context [START_REF] Dua | UCI machine learning repository[END_REF], splice [START_REF] Dua | UCI machine learning repository[END_REF], and skating [START_REF] Mörchen | Efficient mining of understandable patterns from multivariate interval time series[END_REF]. We also apply our algorithm on the real life dataset sc2 that has been used in [START_REF] Bosc | A Pattern Mining Approach to Study Strategy Balance in RTS Games[END_REF]. It was extracted from Starcraft IIa famous Real Time Strategy (RTS) game. The goal of each Starcraft match is to destroy units and buildings of the opponent. Three different factions exist, each with its own features (i.e., combat units, buildings, strategies). Roughly speaking, there is a duality between economy and military forces. Investing in military forces is important to defend or attack an opponent, while building economy is important to have more resources to invest into military forces. Sequences in sc2 correspond to the buildings constructed during a game, and the class corresponds to the winner's faction. Once the target class (i.e., the winning faction) has been chosen, our algorithm will look for patterns of construction that characterize the victory of this faction.

Moreover, we use jmlr, a dataset consisting of abstracts of articles published in the Journal of Machine Learning Research [START_REF] Tatti | The long and the short of it: Summarising event sequences with serial episodes[END_REF]. Table II summarizes the statistics of datasets.

Baselines: We are not aware of available algorithms that address the problem of discriminative pattern mining in sequences of itemsets but several exist for processing sequences of events. Therefore, we compare SeqScout to two algorithms that we have modified to tackle sequences of itemsets, namely misère [START_REF] Egho | A user parameter-free approach for mining robust sequential classification rules[END_REF] and BeamSearch [START_REF] Van Leeuwen | Diverse subgroup set discovery[END_REF]. First, we implemented a simple extension of misère [START_REF] Egho | A user parameter-free approach for mining robust sequential classification rules[END_REF]. Second, we implemented a sequence-oriented version of a beam-search algorithm, denoted BeamSearch. To deal with sequences of itemsets, we consider item-extensions and set-extensions at each given depth. Moreover, for the sake of non-redundancy in the returned patterns, we modify its best-first search nature so that the expanded nodes get diverse, as defined in [START_REF] Van Leeuwen | Diverse subgroup set discovery[END_REF]. Indeed, without this modification, BeamSearch selects directly the best patterns at each level, that are often redundant. It may lead to a situation where almost all patterns are removed in the non-redundancy postprocessing step, giving rise to solutions with less than k elements.

Settings: If not stated otherwise, we use the following settings. Each algorithm has been launched 5 times, and the reported results are averaged over these runs. For BeamSearch, we empirically set the parameter width = 50. For all algorithms, we set θ = 0.5, time budget = ∞, iteration num = 10, 000, and top k = 5. Note that instead of giving a fixed time budget for running an algorithm on each dataset, we chose to limit the number of iterations iteration num, one iteration corresponding to a single computation of the quality measure. Indeed, this computation is the most time consuming one: it needs for computing the extent w.r.t. the whole dataset. Therefore, using the same time budget on different datasets would not provide a fair comparison: having 50,000 iterations on a small dataset versus 50 on a complex one with the same time budget is not interesting.

A. Performance evaluation using W RAcc

To assess the performance of the algorithms, let us first use the mean of the W RAcc of the top-k non redundant patterns given by algorithms misère, BeamSearch and SeqScout. Fig. 3 provides the comparative results. Our approach gives better results on the majority of datasets. Moreover, it remains quite stable, contrarily to BeamSearch which can be sometimes fairly inefficient (See on splice).

B. Quality with respect to the number of iterations

We discuss the result quality in terms of W RAcc over the number of iterations. Fig. 4, 5, 6 depict the results for the top-5 non-redundant patterns on datasets promoters, skating and sc2. Note that for the same data, the results may vary from run to run, due to the random component of misère and SeqScout. It explains some fluctuations of the quality. Nevertheless, for each iteration num setting, SeqScout has shown better results.

C. Using other quality measures

To empirically illustrate the measure agnostic characteristic of SeqScout, we have used other quality measures: F 1score and Inf ormedness. The results are shown in Table III. Our algorithm gives generally better results.

D. Performance study under varying θ

We evaluate the performance of the algorithms when varying the similarity threshold θ. Fig. 7 shows the performance on the dataset promoters. We did not included the results for θ = 0 because it would mean finding patterns with totally disjoint extents. It results in finding a number of patterns lesser than k for all algorithms, such that the mean would be misleading. We can see from the plot that SeqScout outperforms other algorithms for all θ values.

E. Performance study under varying k (top-k)

We investigate the performance of the search for top-k patterns when changing the k parameter. Fig. 8 shows the results when considering the context dataset. SeqScout gives better results. Note that the mean W RAcc decreases for all algorithms, as increasing k leads to the selection of lower quality patterns.

F. Quality versus search space size

We evaluate the average W RAcc of the top-5 patterns w.r.t. the maximum length of sequences (See Fig. 9). To do so, we have truncated the dataset to control the maximum lengths of sequences. We demonstrate it on the skating dataset. The plot shows that SeqScout gives better average W RAcc values whatever the search space size is. We also note a general decrease of quality when increasing the search space size. Indeed, some patterns that are good for a smaller maximum length can appear in negative elements for larger maximum lengths, resulting in a decreasing quality of patterns. Note that the opposite phenomenon can also appear.

G. Comparison to ground truth

We now compute the ratio of the W RAcc result quality of SeqScout to the one of the ground truth. Therefore, we have launched an exhaustive algorithm (Breadth-First-Search) on the dataset sc2. However, the search space size makes the problem intractable in a reasonable time. To tackle this issue, we truncate our dataset to have sequences with a length of 10. Note that even with this severe simplification, our exhaustive algorithm took 35 hours to conclude, and found more than 40 billions of patterns. As we can see in Fig. 10, SeqScout improves over time, and it reaches more than 90% of the best pattern set in less than 2,000 iterations (35 seconds on our machine).

H. Sequence length

The pattern lengths on all datasets are reported in Fig. 11. Let us consider the splice dataset: BeamSearch gives short patterns (max 8), which is significantly less than SeqScout and misère. This may explain why the BeamSearch result quality is bad on this dataset (See Fig. 3). One hypothesis could be that it does not have enough iterations to go deep enough in the search space. Another hypothesis is that BeamSearch cannot find good patterns having bad parents w.r.t. W RAcc: its good patterns are short ones.

I. Bitset vs. Integer set representation

We investigate the usefulness of bitset representation by comparing it against an integer set representation where each item is represented by an integer. Therefore, we have compared the number of iterations SeqScout made for a fixed time budget on each dataset. We set time budget = 10 seconds. The results are summarized in Table IV. We can see that the bitset representation gives a performance gain for datasets with smaller search space size upper bound. Indeed, context, skating, and jmlr have sequences of large lengths, leading to large bitset representations. Those bitsets are then split into different parts to be processed by the CPU. If the number of splits is too large, the bitset representation becomes inefficient, and using a classical integer set representation is a better option. During the local optima search step, an additional number of iterations is made to improve the qualities of the best patterns. In Fig. 12, we plot the ratio of the additional iterations necessary for local optima search w.r.t. the number of iterations given in the main search. The more iterations we have, the more negligible the ratio is. Based on that analysis, we may infer that a rather fair number of iterations for our experiments is 10,000. However, note that we did not plot the additional cost of jmlr. Indeed, in the particular case of text data, the number of possible items is quite large, leading to a fairly long local optima search. Consequently, we note that SeqScout may not be the relevant choice with this kind of dataset.

K. Quality evaluation on sc2

Given a faction, what are the best strategies to win? SeqScout will look for patterns of construction which are characteristic to the victory of this faction. Let us consider the "Terran" faction, as an example. For this faction, one of the best pattern found is: {Hatchery}, {Supply}, {F actory}, {Supply} . We use the colors as indicators of fractions: blue for "Terran" and purple for "Zerg". We can see that the "Terran" player is investing in military units ({Supply} and {F actory}). The "Zerg" player chooses to invest in its economy by creating a {Hatchery}: she fosters the so-called late game, sacrificing its military forces in the so-called early game. In such a scenario, the "Terran" strikes in early game, knowing its military advantage, and she tends to win the match. This example shows that our algorithm can provide relevant patterns that may be useful, e.g., to identify unbalanced strategies [START_REF] Bosc | A Pattern Mining Approach to Study Strategy Balance in RTS Games[END_REF].

L. Quality evaluation on jmlr

Let us now use a dataset consisting of abstracts of articles from JMLR journal [START_REF] Tatti | The long and the short of it: Summarising event sequences with serial episodes[END_REF]. Sentences have been tokenized and stemmed. Each sequence is then a sequence of items, where an item corresponds to a stemmed word. To label sequences, we chose to add a class '+' to sequences containing a target word, and a class '-' to others. We also removed the target word from the sequences, as it would lead to the discovery of patterns having the target word, which is not interesting. The goal is then to find sequences of items (words), characteristic for the target word. The results for two target words (svm and classif ) are presented in Table V. As we can see, these patterns are indeed characteristics of the target words.

VI. CONCLUSION

We presented the algorithm SeqScout to discover relevant subgroups in sequences of itemsets. Though we are not aware of available algorithms to solve the same problem, we have implemented adaptations of two other algorithms when considering sequences of itemsets, namely misère and BeamSearch. Our experiments showed that SeqScout gave better results, without the need for additional parameter tuning as in the case of Beam Search. However we can note that due to the nature of the bandit approach, this algorithm perfoms better on dataset of reasonnable size. Parallelizing the algorithm with an efficient programming language could be interesting to tackle this issue in future works. We plan to further improve our method by incorporating the Monte Carlo Tree Search, a logical evolution of bandit based method to sample the search space. Indeed, instead of repeatedly generalizing promising dataset sequences, it could be more interesting to generalize promising patterns that occur in such promising sequences, identifying the best parts of the search space as the number of iterations increases.

Figure 1 :

 1 Figure 1: Illustration of SeqScout.

9 :

 9 seq, qual, N i ← scores.bestU CB() 10: seq p , qual p ← P layArm(seq) 11: π.add(seq p , qual p ) 12:

17 : 18 :

 1718 function OPTIMIZE(π) 19: topK ← π.topKN onRedundant() 20: for all pattern in topK do 21: while pattern is not a local optima do 22: pattern, qual ← BestN eighbor()

Figure 2 :

 2 Figure 2: Bitset representation and support computing.

Figure 3 :

 3 Figure 3: W RAcc of top-5 best patterns (10K iterations)

Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 456789101112 Figure 4: W RAcc for top-5 patterns w.r.t. iterations (skating)

Table II :

 II Datasets

	Dataset	|D|	|I|	lmax	Search Space Size
	aslbu [27]	441	124	27	4.45 * 10 60
	promoters [12]	106	4	57	1.58 * 10 41
	blocks [12]	210	8	12	2.74 * 10 12
	context [12]	240	47	123	5.24 * 10 224
	splice [12]	3,190	8	60	3.28 * 10 62
	sc2 [16]	5,000	30	30	6.48 * 10 48
	skating [28]	530	41	120	1.15 * 10 212
	jmlr [29]	788	3,836	228	1.84 * 10 853

Table III :

 III Mean values of measures for top-5 patterns.

	Dataset		Informedness		F1
		misere BeamS	SeqScout	misere BeamS	SeqScout
	aslbu	0.195	0.198	0.198	0.505	0.505	0.505
	promoters	0.039	0.089	0.113	0.512	0.545	0.580
	blocks	0.399	0.394	0.393	0.381	0.370	0.382
	context	0.436	0.437	0.469	0.528	0.532	0.584
	splice	0.342	0.006	0.356	0.394	0.082	0.471
	sc2	0.004	0.004	0.010	0.524	0.519	0.526
	skating	0.396	0.385	0.415	0.361	0.380	0.389
	jmlr	0.224	0.403	0.355	0.143	0.236	0.243

Table IV :

 IV Number of found local optima -10s budget

	Dataset	Integer Set	Bitset	Variation(%)
	aslbu	16,453	17,330	5
	promoters	7,185	8,858	24
	blocks	30,725	45,289	51
	context	4,651	2,667	-47
	splice	289	254	-12
	sc2	943	605	-36
	skating	4,001	1,283	-68
	jmlr	704	31	-95

Table V :

 V Top-10 non-singleton patterns given by SeqScout on jmlr for target words: "svm" and "classif "

	svm	classif
	effici method	classifi set
	requir support	classifi value
	propos paper	gener real
	formul classif	problem condit
	individu show	learn error
	scale vector	support gener
	base number support	paper propos
	support machin present show problem classification set
	machin techniqu	problem appli
	solut number	data gener

In the context of sequential pattern mining, the search space is a priori infinite. However, we can define the border of the search space (the bottom border in Fig.1) by excluding patterns having a null support. We can easily prove that each element of this border is a sequence within the database. Therefore, the search space shape depends on the data.

Acknowledgment This research has been partially funded by the French National project FUI DUF 4.0 2017-2021.


Marwa Turki

Nicolas Langlois

Adnan Yassine

An Analytical Tuning of MPC Control Horizon Using the Hessian Condition Number

Keywords: Predictive control, control horizon tuning, analytical approach, effective rank, condition number, stability guarantee, SISO system

HAL is

I. INTRODUCTION

MPC has proven to be an excellent candidate for controlling complex systems and is now widely implemented in industry for many years [START_REF] David | Application of generalized predictive control to industrial processes[END_REF], [START_REF] Qin | A survey of industrial model predictive control technology[END_REF]. Its parameters N c , N p and λ influence significantly the closed-loop behavior, stability and robustness in a complex manner [START_REF] Qi | Robust stability of model predictive control[END_REF], [START_REF] Ankush | Effect of tuning parameters on performance of first-order plus dead-time processes using generalized predictive control[END_REF]. MPC parameters tuning is, therefore, a challenging issue. Over the last decades, many research efforts have led to the development of MPC tuning approaches, [START_REF] David W Clarke | Generalized predictive control-part i. the basic algorithm[END_REF], [START_REF] Jorge | Model predictive control tuning methods: A review[END_REF], [START_REF] Rc Gutiérrez-Urquídez | Systematic selection of tuning parameters for efficient predictive controllers using a multiobjective evolutionary algorithm[END_REF]. We can classify these methods on three categories:

• The first category includes the analytical methods (more or less numerical-based). Only few work has been published due to the complexity of the problem [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF], [START_REF] Shridhar | A novel tuning strategy for multivariable model predictive control[END_REF], [START_REF] Bagheri | Analytical approach to tuning of model predictive control for firstorder plus dead time models[END_REF], [START_REF] Bagheri | An analytical tuning approach to multivariable model predictive controllers[END_REF], [START_REF] Gholaminejad | Adaptive tuning of model predictive control based on analytical results[END_REF].

• The second category includes the heuristic methods. It aims to find an approximation of optimal values for the MPC. Different approaches have been extensively published in the literature. Some of them are based on fuzzy logic [START_REF] Mamboundou | Commande prédictive à base d'état adaptive indirecte supervisée par la logique floue[END_REF], [START_REF] Mamboundou | Indirect adaptive model predictive control supervised by fuzzy logic[END_REF], [START_REF] Ali | Automatic tuning of model predictive controllers based on fuzzy logic[END_REF] on genetic algorithms [START_REF] Haber | Optimal choice of horizons for predictive control by using genetic algorithms[END_REF], [START_REF] Ben Aicha | A multivariable multiobjective predictive controller[END_REF] and on neural networks [START_REF] Tohidi | Self-tuning adaptive multiple model predictive control with application to ph control process[END_REF], [START_REF] Han | Nonlinear model predictive control based on a selforganizing recurrent neural network[END_REF]. Other {adnan.yassine@univ-lehavre.fr} choose to use the Analysis of Variance (ANOVA) and nonlinear regression [START_REF] Ebrahimi | Tuning of generalized predictive controllers for first order plus dead time models based on anova[END_REF].

• The third category includes the empirical approaches where the MPC parameters are determined regarding the designer's experience [START_REF] Mohammad B Shadmand | Auto-tuning the cost function weight factors in a model predictive controller for a matrix converter var compensator[END_REF], [START_REF] Valencia | Programmable logic controller implementation of an auto-tuned predictive control based on minimal plant information[END_REF], [START_REF] Chang | Online parameter tuning technique for predictive current-mode control operating in boundary conduction mode[END_REF], [START_REF] Boucher | La commande prédictive[END_REF].

The majority of these methods have a trial-and-error nature, which does not always permit to identify explicitly the robustness area of the control system. None of the published methods deals with a general case of a process whatever its order and none of these methods offers a setting of the control horizon while taking into account the constraints. In this paper, we intend to overcome these limits thanks to our original approach by computing the control horizon while enhancing the numerical condition of the controlled process regarding an analytical method. Important abilities of our approach are highlighted when linear inequality constraints have to be considered. Indeed, the tuning strategy easily integrates the constraints without increasing the calculation effort or undoing the proposed optimization algorithm. Its advantage lies on the fact that it is applicable to any controllable and observable SISO linear system. This paper is outlined as follows: section II reminds the theoretical background on MPC based on state-space representation. Section III highlights the proposed analytical tuning approach for the control horizon. In section IV, a comparative study of the obtained performances on a simulated pH neutralization process is carried out. Simulations results are shown to emphasize the effectiveness of our approach.

II. REMIND ON MODEL PREDICTIVE CONTROL DESIGN

A. Augmented state-space model

Let consider a SISO system represented by the following discrete-time state-space model:

x m (k + 1) = A m x m (k) + B m u(k) y(k) = C m x m (k) (1)
Where y ∈ R is the output system, u ∈ R the manipulated variable and the row matrix x m is the state variable vector of size n Am . k is the sampling instant (positive value). In (1), A m ∈ n Am ×n Am , B m ∈ n Am ×1 and C m ∈ 1×n Am are the state-space matrices.

In order to design predictive controller, let adopt the formulation of augmentated-state-model with embedded integrators whose advantages have been already discussed [START_REF] Wang | Model predictive control system design and implementation using MATLAB R[END_REF]. Such a formulation has been inspired by the integral functionality [START_REF] David W Clarke | Generalized predictive control-part i. the basic algorithm[END_REF].

The difference of the state and control variables are:

∆x m (k) = x m (k) -x m (k -1) ∆u(k) = u(k) -u(k -1) (2)
1

By integrating the increments of the variables x m (k) et u(k). Equation (2) becomes:

∆x m (k + 1) = A m (x m (k) -x m (k -1)) + B m (u(k) -u(k -1)) = A m ∆x m (k) + B m ∆u(k).
Now, let consider a new state space vector as:

x(k) = ∆x m (k) y(k)
Finally, the augmented state-space model is as following:

x(k+1)

∆xm(k + 1) y(k + 1) = A Am 0 t m CmAm 1 x(k) + B Bm CmBm ∆u(k) y(k) = C 0m 1 ∆xm(k) y(k) (3)
Where

0 m = [0 0 • • • 0] is a row matrix of size n Am . The matrices A, B and C of size respectively (n A × n A ), (n A × 1)
and (1 × n A ) (with n A = n Am + 1), constitute the discrete time augmented state-space representation.

B. MPC formulation

As an hypothesis, the system is supposed to be observable and controllable. The incremental control signal vector ∆U of size (1 × N c ) is defined as:

∆U = ∆u(k) ∆u(k + 1) • • • ∆u(k + N c -1) t (4)
While superscript t denotes matrix transpose. The desired output Y des of size (N p × 1) is:

Y des = y des (k + 1) y des (k + 2) • • • y des (k + N p ) t
Assuming the predicted output vector Ŷ defined by [START_REF] Wang | Model predictive control system design and implementation using MATLAB R[END_REF]

Ŷ = F x(k i ) + Φ∆U (5)
where

Ŷ = ŷ(k + 1 | k) ŷ(k + 2 | k) • • • ŷ(k + N p | k) t , F = CA CA 2 CA 3 • • • CA Np t , (6)
and

Φ =        CB 0 0 ... 0 CAB CB 0 ... 0 CA 2 B CAB CB ... 0 . . . CA Np-1 B CA Np-2 B CA Np-3 B ... CA Np-Nc B        (7)
Let the cost function J to be minimized as follows:

J = (Y des -Ŷ ) t (Y des -Ŷ ) + ∆U t R∆U ( 8
)
where R is a (N c × N c ) matrix defined by: R = λI (Nc×Nc) [START_REF] Shridhar | A novel tuning strategy for multivariable model predictive control[END_REF] and I denotes the identity matrix.

In order to calculate the optimum value of the manipulated variable, let the partial derivation of J with respect to ∆U as:

∂J ∂∆U = -2Φ t (Y des -F x(k)) + 2(Φ t Φ + R)∆U (10)
The optimal control sequence ∆U is found by solving [START_REF] Bagheri | Analytical approach to tuning of model predictive control for firstorder plus dead time models[END_REF] equals to zero, which leads to:

∆U = (Φ t Φ + R) -1 Φ t (Y des -F x(k)) (11)
According to the receding horizon principle, the control signal applied to the process will be the first element of ∆U as:

u(k) = u(k -1) + ∆u(k) u(k) = u(k -1) + I (1×Nc) (Φ t Φ + R) -1 Φ t (Y des -F x(k))
Note that the estimation of the optimal control requires the state vector x(k). Here the control signal is subject to linear inequality constraints as proposed by [START_REF] Wang | Model predictive control system design and implementation using MATLAB R[END_REF]:

U min ≤ U ≤ U max . (12)

III. ANALYTICAL TUNING OF THE CONTROL HORIZON

This section is devoted to explain the analytical method we propose to tune N c .

Most often in the literature, N c is taken equal to unity as proposed in the guidelines provided by [START_REF] David W Clarke | Generalized predictive control-part i. the basic algorithm[END_REF] and lower than the predicted horizon N p . In fact, considering N c = 1 gives acceptable control performance. A higher N c is recommended when systems to control have unstable poles [START_REF] David W Clarke | Generalized predictive control-part i. the basic algorithm[END_REF].

A. The Hessian condition number

In order to compute analytically the optimal value of N c , the concept of numerical stability is considered. Indeed, the numerical stability concerns mainly the condition number of a square nonsingular matrix [START_REF] Michael | Scientific computing: an introductory survey[END_REF]. Generally speaking, this specific matrix is the Hessian matrix (or just the Hessian) of an algorithm. Dealing with MPC, the Hessian matrix is present in the formulation of the optimal control sequence [START_REF] Keyser | Model based predictive control for linear systems[END_REF]. With the assumption that (Φ t Φ + R) -1 exists, let consider the Hessian matrix H of size (N c × N c) defined as:

H = (Φ t Φ + R) -1 (13)
To evaluate the conditioning of the Hessian matrix H, one calculates its condition number defined by:

cond(H(k)) = H(k) 2 . H(k) -1 2 = σ max (k) σ min (k) (14)
Where σ max (k) and σ min (k) are respectively the maximum and the minimum singular values of matrix H. Then, the condition number of a matrix indicates how close a matrix is to be singular: a matrix with a large condition number is nearly singular, wheras a matrix with a condition number close to unity is far from being singular [START_REF] Michael | Scientific computing: an introductory survey[END_REF].

Regarding the literature, different ways exist to improve the condition number of a Hessian matrix: thus [START_REF] Wills | Interior-point methods for linear model predictive control[END_REF] suggests using a numerical stable projection or re-writing the Hessian formulas. [START_REF] Laurí | Data-driven latent-variable model-based predictive control for continuous processes[END_REF] improves the condition number of the Hessian matrix using a latent variable method with MPC. [START_REF] Rojas | An svd based strategy for receding horizon control of input constrained linear systems[END_REF] establishes a strategy combining a singular value decomposition (SVD) and a receding horizon control (RHC) principle to enhance the Hessian conditioning. A sub-optimal control signal is produced by discarding the smallest singular values of the Hessian. De Keyser recommends solving this numerical problem using principal component analysis (PCA) [START_REF] Keyser | Model based predictive control for linear systems[END_REF]. This method looks like the SVD-RHC technique of [START_REF] Rojas | An svd based strategy for receding horizon control of input constrained linear systems[END_REF]. Both approaches cannot systematically leads to satisfactory results in comparison with the method proposed in [START_REF] Wang | Model predictive control system design and implementation using MATLAB R[END_REF], [START_REF] Wang | Use of exponential data weighting in model predictive control design[END_REF] and [START_REF] Wang | Discrete model predictive control using laguerre functions: numerical sensitivity analysis[END_REF]. This one uses an exponential data weighting in the cost function to enhance the Hessian conditioning. This additive weighting technique yields a more straightforward and practical method for engineers and researchers. Here in this paper, we intend to overcome these limits and enhance the condition number of the Hessian matrix by computing analytically the control prediction.

B. Improving the condition number: a dimension reduction problem

As shown in [START_REF] Scattolini | On the choice of the horizon in long-range predictive control-some simple criteria[END_REF] and [START_REF] Marco | An infinite horizon model predictive control for stable and integrating processes[END_REF], high MPC horizons guarantee the system closed-loop stability. Ideally, these values tend towards infinity (N c → ∞, N p → ∞). Then, we notice that the issue of improving the condition number of the Hessian matrix is converted into a dimension reduction problem of matrix H. Thus, when the dimension reduction problem is solved, an optimal value of N c is deduced (Figure 1). In order to reduce the dimension of matrix H, the concept of the effective rank (ER) is considered [START_REF] Roy | The effective rank: A measure of effective dimensionality[END_REF] (Appendix 1).

C. Relation between the Effective Rank and N c

In this part, we relate the concept of the effective rank with MPC design in order to evaluate the optimum value of the control horizon. The different steps for calculating the optimal value of N c are as follows:

1) Tending N c → ∞ and N p → ∞ (N c < N p, N c ∈ N * and N p ∈ N * ). 2) Taking A ER = H.
3) Evaluating Q defined as follows: 7) Evaluating finally the optimal value of N c by solving the following minimization problem using Yalmip toolbox [START_REF] Lofberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF]:

Q = min{M ER , N ER } = min{N c, N c} = N c.

4) Decomposing into singular values the matrix H and evaluating the matrix

σ = [σ 1 σ 2 • • • σ ∞ ] t : H = U H D H V H σ 1 ≥ σ 2 ≥ • • • ≥ σ ∞ (15) (16
     N copt = round(e HShannon(p1,p2,••• ,p∞) ) min cond(H(k)) -1 N copt (17)
Under the following constraints:

   N c ∈ N * 1 ≤ N c < N p U min ≤ U U max
As a practical advantage, this strategy does not require an apriority knowledge of the minimum rank and it takes no account of the full singular value spectrum [START_REF] Roy | The effective rank: A measure of effective dimensionality[END_REF]. As a conclusion, the smallest control prediction we take, the better condition number is. Thus, the numerical stability is enhanced. Note 1 Linear inequality constraints considered in the MPC do not influence the tuning approach proposed in this paper. Note 2 [START_REF] Wang | Model predictive control system design and implementation using MATLAB R[END_REF] shows that the minimization of MPC horizons enhance the numerical condition of the controlled process. One way to reduce the control horizon is using the concept of the effective rank.

IV. APPLICATION TO A SIMULATED CHEMICAL PROCESS

As a benchmark, a simulated pH neutralization process is considered in this part [START_REF] Michael | Adaptive nonlinear control of a ph neutralization process[END_REF], [START_REF] Nie | Modeling ph neutralization processes using fuzzy-neural approaches[END_REF].

A. Performances criteria

In order to carry out a comparative study between the conventional guideline proposed by Shridhar [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF] and our approach, the following performance indexes are considered:

1) The stability degree index (SDI) is used to evaluate the system closed-loop stability. It is the difference between the radius of the unit circle and the modulus of the pole most remote from the unit circle centre [START_REF] Mamboundou | Indirect adaptive model predictive control supervised by fuzzy logic[END_REF]:

SDI(k) = 1 -max(|p 1 | , |p 2 | , • • • , |p A N |). ( 18
)
where

{p 1 , p 2 , • • • , p A N } are the eigenvalues of A cl = A -BK mpc . So the closed-loop system is stable when SDI ∈]0, 1[. 2)
The variance of control signal (VARU) makes it possible to observe the mean value of the square deviations of u(k) from its average as follows [START_REF] Ben | Identification et commande numérique des procédés industriels[END_REF]:

V ARU (k) = u(k) 2 -u(k) 2 ( 19
)
where ū is the mean of u.

3) The rise time (RT) is from 10% to 90% [START_REF] Ben | Identification et commande numérique des procédés industriels[END_REF]. 4) The settling time (ST) is within 2%. 5) The overshoot (OV) is the overshoot in the output signal. 6) The static error (SE) 7) The control signal energy (CSE) is:

CSE = ST k=0 u 2 (k) (20)
8) The control effort energy (CEE) is:

CEE = ST k=0 ∆u 2 (k) (21)
9) The computational load (CL) depends heavily on the control horizon [START_REF] Ben | Identification et commande numérique des procédés industriels[END_REF]. 10) the condition number of H (Cond(H) is computed from equation ( 14).

B. Process description

The considered benchmark is a nonlinear multivariable system. For simplification reasons, we deal with a SISO pH neutralization process. The process consists of acid, base and buffer streams mixed in a vessel. A schematic diagram of the studied process is presented in Figure 2. In the SISO case, acid (HN O 3 ) stream represents the measured system disturbance, the base (N aOH) stream is the control signal and the buffer (N aHCO 3 , N aOH) stream represents the unmeasured disturbance of the system.

The main objective of this part is to control the value of the pH of the outlet stream. In addition, it is assumed that the pH of the outlet stream is measured at a distance from the plant, which introduces a measurement time delay θ. The sampling time T s is set to 8 ms as in [START_REF] Bagheri | An anova based analytical dynamic matrix controller tuning procedure for fopdt models[END_REF].

C. Process linearisation

Since the SISO system is highly nonlinear and for control purpose, we adopt the linear model constructed by [START_REF] Bagheri | An anova based analytical dynamic matrix controller tuning procedure for fopdt models[END_REF] around an operation point (pH = 7) whose the augmented-state representation in the discrete time domain is the following:

A = 0.9102 0 0.0397 1 , B = 1 0.0436 , C = [0 1] . (22)

D. Tuning parameters

The prediction horizon and the weighting factor values are chosen according to [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF]. Based on the tuning approach described in Section III, the control horizon is analytically computed as follows:

• Initializing N c = 39, N p = 40 and λ = 0.063.

+ θ Ts + 1 λ f K 2 f K 2 f 0 Nc = 1 Nc 500 ( 3.5τ Ts -Nc-5 2 ) Nc > 1 0 Nc = 1 Nc 500 ( 3.5τ Ts -Nc-5 2 ) Nc > 1
• Decomposing H into singular values, evaluating σ and computing the Shannon entropy value based on ( 15), ( 16) and ( 27) as:

H Shannon = 0.5530

• Computing the optimal value of N c based on [START_REF] Ben Aicha | A multivariable multiobjective predictive controller[END_REF] as:

N copt = e HShannon(p1,p2,••• ,p∞) = 1.7385 2

E. Simulation test

The tuning strategy detailed in [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF] propose guideline for tuning a First Order Plus Dead Time (FOPDT) model whose transfer function is, in the general case as follows:

G F OP DT (s) = Ke -θs τ s + 1 (23)
The proposed tuning equations are shown in Table I. The simulation results are depicted in Figure 3. The parameters II.

The linear inequality constraint of the control signal is taken as: 0 ml/s ≤ u ≤ 25 ml/s. As shown in Figure 3, the fastest response is obtained using Schridhar approach [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF]. However, this latter leads to the highest VARU, CSE and CEE indices. The best SDI is obtained thanks to the guideline of [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF] but this method causes a high overshoot and a considerable static-error. The condition number computed with our approach is about 200 times smaller than the one given by [START_REF] Shridhar | A tuning strategy for unconstrained siso model predictive control[END_REF] which indicate a remarkable improvement in the calculation conditions. In conclusion, only with the computation of the control horizon, few performance criteria have been greatly improved. In the work perspective, an analytical method for calculating the prediction horizon and the weighting factor is expected.

V. CONCLUSION

An analytical approach with an enhancement of numerical condition to tune the control horizon value has been proposed in this paper. This approach is dedicated to nonlinear systems around the operating points. Some advantages of our novel approach are highlighted when constraints have to be considered. Indeed, the tuning strategy easily integrates the linear inequality constraints without increasing the computational effort or undoing the proposed optimization algorithm. Future work will aim to find an analytical approach for computing the prediction horizon and the weighting factor.

VI. APPENDIX 1

Effective Rank concept Let consider a complex-valued non-all-zero matrix A ER of size (M ER × N ER ) whose SVD is given by A ER = U ER D ER V ER . Where U ER and V ER are unitary matrices of size (M ER × M ER ) and (N ER × N ER ), respectively, and D ER is an (M ER × N ER ) diagonal matrix containing the real positive singular values:

σ 1 ≥ σ 2 ≥ • • • ≥ σ Q ≥ 0, with Q = min{M ER , N ER }.
Let define for writing simplification σ = σ 1 σ 2 • • • σ Q t and the singular value distribution:

p k = σ k σ 1 (25)
where . 1 denotes the l 1 norm defined as σ 1 = Q k=1 |σ k | Definition VI.1. Effective Rank The effective rank of the matrix A ER is defined as:

erank(A ER ) = exp H Shannon (p 1 , p 2 , • • • , p ∞ ) (26)
where H Shannon (p 1 , p 2 , • • • , p Q ) is the Shannon entropy:

H Shannon (p 1 , p 2 , • • • , p Q ) = - Q k=1 p k log p k . (27)

Fig. 1 .

 1 Fig. 1. Proposed Nc tuning strategy

) 5 )

 5 Evaluating the singular value distribution pk with k = [1, 2, • • • , ∞]. 6)Computing the Shannon entropy defined in[START_REF] Keyser | Model based predictive control for linear systems[END_REF].

Fig. 2 .

 2 Fig. 2. Sketch of the pH neutralization process

Fig. 3 .

 3 Fig. 3. Output and control signals vs. time

TABLE II .

 II MPC PARAMETERS AND PERFORMANCE COMPARISONS OF MPC TUNING METHODS FOR THE PH NEUTRALIZATION PROCESS

		Shridhar and Cooper 1997 [8] Proposed approach
	Ts(ms)	8	8
	Nc	4	2
	Np	40	40
	λ	0.063	0.063
	SDI	0.2850	0.2789
	VARU	37.3001	12.7142
	RT (s)	7.2537	8.3255
	ST (s)	14.2457	15.2734
	OV (%)	1.7938	0
	SE	0.6159	0.2965
	CSE (e + 03)	6.3796	4.8804
	CEE (e + 03)	4.3775	0.4237
	CL	5.8155	4.6622
	Cond(H)	352.3533	1.7244