Vincent Le Guen
email: vincent.le-guen@edf.fr

Nicolas Thome
email: nicolas.thome@cnam.fr

Shape and Time Distortion Loss for Training Deep Time Series Forecasting Models

This paper addresses the problem of time series forecasting for non-stationary signals and multiple future steps prediction. To handle this challenging task, we introduce the Shape and Time Distortion Loss (STDL), a new objective function dedicated to training deep neural networks. STDL aims at accurately predicting sudden changes, and explicitly incorporates two terms supporting precise shape and temporal change detection. We introduce a differentiable loss function suitable for training deep neural nets, and provide a custom back-prop implementation for speeding up optimization. We also introduce a variant of STDL, which provides a smooth generalization of temporally-constrained Dynamic Time Warping (DTW). Experiments carried out on various non-stationary datasets reveal the very good behaviour of STDL compared to models trained with the standard Mean Squared Error (MSE) loss function, and also to DTW and variants. STDL is also agnostic to the choice of the model, and we highlight its benefit for training fully connected networks as well as specialized recurrent architectures, showing its capacity to improve over state-of-the-art trajectory forecasting approaches.

Introduction

Time series forecasting [START_REF] George Ep Box | Time series analysis: forecasting and control[END_REF] consists in analyzing the dynamics and correlations between historical data for predicting future behavior. In one-step prediction problems [START_REF] Yao Qin | A dual-stage attention-based recurrent neural network for time series prediction[END_REF][START_REF] Lai | Modeling long-and shortterm temporal patterns with deep neural networks[END_REF], future prediction reduces to a single scalar value. This is in sharp contrast with multi-step time series prediction [START_REF] Ben Taieb | A review and comparison of strategies for multi-step ahead time series forecasting based on the NN5 forecasting competition[END_REF][START_REF] Hoang | Comparison of strategies for multi-step-ahead prediction of time series using neural network[END_REF][START_REF] Ben | A bias and variance analysis for multistep-ahead time series forecasting[END_REF], which consists in predicting a complete trajectory of future data at a rather long temporal extent. Multi-step forecasting thus requires to accurately describe time series evolution. This work focuses on multi-step forecasting problems for non-stationary signals, i.e. when future data cannot only be inferred from the past periodicity, and when abrupt changes of regime can occur. This includes important and diverse application fields, e.g. regulating electricity consumption [START_REF] Zheng | Electric load forecasting in smart grids using long-short-term-memory based recurrent neural network[END_REF][START_REF] Masum | Multi-step time series forecasting of electric load using machine learning models[END_REF], predicting sharp discontinuities in renewable energy production [START_REF] Ghaderi | Deep forecast: Deep learning-based spatio-temporal forecasting[END_REF] or in traffic flow [START_REF] Lv | Traffic flow prediction with big data: a deep learning approach[END_REF][START_REF] Li | Diffusion convolutional recurrent neural network: Data-driven traffic forecasting[END_REF], electrocardiogram (ECG) analysis [START_REF] Chauhan | Anomaly detection in ECG time signals via deep long short-term memory networks[END_REF], stock markets prediction [START_REF] Ding | Deep learning for event-driven stock prediction[END_REF], etc.

Deep learning is an appealing solution for this multi-step and non-stationary prediction problem, due to the ability of deep neural networks to model complex nonlinear time dependencies. Many approaches have recently been proposed, mostly relying on the design of specific one-step ahead 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada. architectures recursively applied for multi-step [START_REF] Girard | Multiple-step ahead prediction for non linear dynamic systems -a gaussian process treatment with propagation of the uncertainty[END_REF][START_REF] Hussein | Multi-step-ahead chaotic time series prediction using coevolutionary recurrent neural networks[END_REF][START_REF] Chandra | Co-evolutionary multi-task learning with predictive recurrence for multi-step chaotic time series prediction[END_REF][START_REF] Borovykh | Conditional time series forecasting with convolutional neural networks[END_REF], on direct multi-step models [START_REF] Bao | Multi-step-ahead time series prediction using multipleoutput support vector regression[END_REF] such as Sequence To Sequence [START_REF] Li | Diffusion convolutional recurrent neural network: Data-driven traffic forecasting[END_REF][START_REF] Yu | Long-term forecasting using tensor-train RNNs[END_REF][START_REF] Wen | A multihorizon quantile recurrent forecaster[END_REF][START_REF] Yu | Learning chaotic dynamics using tensor recurrent neural networks[END_REF] or State Space Models for probabilistic forecasts [START_REF] Salinas | DeepAR: Probabilistic forecasting with autoregressive recurrent networks[END_REF][START_REF] Sundar Rangapuram | Deep state space models for time series forecasting[END_REF].

Regarding training, the huge majority of methods use the Mean Squared Error (MSE) or its variants (MAE, etc) as loss functions. However, relying on MSE may arguably be inadequate in our context, as illustrated in To train deep neural nets with STDL, we derive a differentiable loss function for both shape and temporal terms (section 3.1), and an efficient and custom back-prop implementation for speeding up optimization (section 3.2). We also introduce a variant of STDL, which provides a smooth generalization of temporally-constrained Dynamic Time Warping (DTW) metrics [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF][START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF]. Experiments carried out on several synthetic and real non-stationary datasets reveal that models trained with STDL significantly outperform models trained with the MSE loss function when evaluated with shape and temporal distortion metrics, while SDTL maintains very good performance when evaluated with MSE. Finally, we show that STDL can be used with various network architectures and can outperform on shape and time metrics state-of-the-art models specifically designed for multi-step and non-stationary forecasting.

Related work

Time series forecasting Traditional methods for time series forecasting include linear autoregressive models, such as the ARIMA model [START_REF] George Ep Box | Time series analysis: forecasting and control[END_REF], and Exponential Smoothing [START_REF] Hyndman | Forecasting with exponential smoothing: the state space approach[END_REF], which both fall into the broad category of linear State Space Models (SSMs) [START_REF] Durbin | Time series analysis by state space methods[END_REF]. These methods handle linear dynamics and stationary time series (or made stationary by temporal differences). However the stationarity assumption is not satisfied for many real world time series that can present abrupt changes of distribution. Since, Recurrent Neural Networks (RNNs) and variants such as Long Short Term Memory Networks (LSTMs) [START_REF] Hochreiter | Long short-term memory[END_REF] have become popular due to their automatic feature extraction abilities, complex patterns and long term dependencies modeling. In the era of deep learning, much effort has been recently devoted to tackle multivariate time series forecasting with a huge number of input series [START_REF] Laptev | Time-series extreme event forecasting with neural networks at Uber[END_REF], by leveraging attention mechanisms [START_REF] Lai | Modeling long-and shortterm temporal patterns with deep neural networks[END_REF][START_REF] Yao Qin | A dual-stage attention-based recurrent neural network for time series prediction[END_REF][START_REF] Tao | Hierarchical attentionbased recurrent highway networks for time series prediction[END_REF][START_REF] Choi | RETAIN: An interpretable predictive model for healthcare using reverse time attention mechanism[END_REF] or tensor factorizations [START_REF] Yu | Long-term forecasting using tensor-train RNNs[END_REF][START_REF] Yu | Temporal regularized matrix factorization for high-dimensional time series prediction[END_REF][START_REF] Sen | Think globally, act locally: a deep neural network approach to high dimensional time series forecasting[END_REF] for capturing shared information between series. Another current trend is to combine deep learning and State Space Models for modeling uncertainty [START_REF] Seeger | Bayesian intermittent demand forecasting for large inventories[END_REF][START_REF] Salinas | DeepAR: Probabilistic forecasting with autoregressive recurrent networks[END_REF][START_REF] Sundar Rangapuram | Deep state space models for time series forecasting[END_REF][START_REF] Wang | Deep factors for forecasting[END_REF]. In this paper we focus on deterministic multi-step forecasting. To this end, the most common approach is to apply recursively a one-step ahead trained model. Although mono-step learned models can be adapted and improved for the multi-step setting [START_REF] Venkatraman | Improving multi-step prediction of learned time series models[END_REF], a thorough comparison of the different multi-step strategies [START_REF] Ben | A bias and variance analysis for multistep-ahead time series forecasting[END_REF] has recommended the direct multi-horizon strategy. Of particular interest in this category are Sequence To Sequence (Seq2Seq) RNNs models1 [START_REF] Salinas | DeepAR: Probabilistic forecasting with autoregressive recurrent networks[END_REF][START_REF] Laptev | Time-series extreme event forecasting with neural networks at Uber[END_REF][START_REF] Yu | Long-term forecasting using tensor-train RNNs[END_REF][START_REF] Wen | A multihorizon quantile recurrent forecaster[END_REF][START_REF] Fox | Deep multi-output forecasting: Learning to accurately predict blood glucose trajectories[END_REF] which achieved great success in machine translation. Theoretical generalization bounds for Seq2Seq forecasting were derived with an additional discrepancy term quantifying the non-stationarity of time series [START_REF] Kuznetsov | Foundations of sequence-to-sequence modeling for time series[END_REF]. Following the success of WaveNet for audio generation [START_REF] Van Den | WaveNet: A generative model for raw audio[END_REF], Convolutional Neural Networks with dilation have become a popular alternative for time series forecasting [START_REF] Borovykh | Conditional time series forecasting with convolutional neural networks[END_REF]. The self-attention Transformer architecture [START_REF] Vaswani | Attention is all you need[END_REF] was also lately investigated for accessing long-range context regardless of distance [START_REF] Li | Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting[END_REF]. We highlight that our proposed loss function can be used for training any direct multi-step deep architecture.

Evaluation and training metrics

The largely dominant loss function to train and evaluate deep models is the MAE, MSE and its variants (SMAPE, etc). Metrics reflecting shape and temporal localization exist: Dynamic Time Warping [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] for shape ; timing errors can be casted as a detection problem by computing Precision and Recall scores after segmenting series by Change Point Detection [START_REF] Chang | Kernel change-point detection with auxiliary deep generative models[END_REF][START_REF] Li | M-statistic for kernel change-point detection[END_REF], or by computing the Hausdorff distance between two sets of change points [START_REF] Garreau | Consistent change-point detection with kernels[END_REF][START_REF] Truong | Supervised kernel change point detection with partial annotations[END_REF]. For assessing the detection of ramps in wind and solar energy forecasting, specific algorithms were designed: for shape, the ramp score [START_REF] Florita | Identifying wind and solar ramping events[END_REF][START_REF] Vallance | Towards a standardized procedure to assess solar forecast accuracy: A new ramp and time alignment metric[END_REF] based on a piecewise linear approximation of the derivatives of time series; for temporal error estimation, the Temporal Distortion Index (TDI) [START_REF] Frías-Paredes | Assessing energy forecasting inaccuracy by simultaneously considering temporal and absolute errors[END_REF][START_REF] Vallance | Towards a standardized procedure to assess solar forecast accuracy: A new ramp and time alignment metric[END_REF]. However, these evaluation metrics are not differentiable, making them unusable as loss functions for training deep neural networks. The impossibility to directly optimize the appropriate (often non-differentiable) evaluation metric for a given task has bolstered efforts to design good surrogate losses in various domains, for example in ranking [START_REF] Durand | Mantra: Minimum maximum latent structural svm for image classification and ranking[END_REF][START_REF] Yue | A support vector method for optimizing average precision[END_REF] or computer vision [START_REF] Nowozin | Optimal decisions from probabilistic models: the intersection-over-union case[END_REF][START_REF] Yu | The lovász hinge: A novel convex surrogate for submodular losses[END_REF].

Recently, some attempts have been made to train deep neural networks based on alternatives to MSE, especially based on a smooth approximation of the Dynamic time warping (DTW) [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF][START_REF] Mensch | Differentiable dynamic programming for structured prediction and attention[END_REF][START_REF] Abid | Learning a warping distance from unlabeled time series using sequence autoencoders[END_REF].

Training DNNs with a DTW loss enables to focus on the shape error between two signals. However, since DTW is by design invariant to elastic distortions, it completely ignores the temporal localization of the change. In our context of sharp change detection, both shape and temporal distortions are crucial to provide an adequate forecast. A differentiable timing error loss function based on DTW on the event (binary) space was proposed in [START_REF] Rivest | A new timing error cost function for binary time series prediction[END_REF] ; however it is only applicable for predicting binary time series. This paper specifically focuses on designing a loss function able to disentangle shape and temporal delay terms for training deep neural networks on real world time series.

Training Deep Neural Networks with the Shape and Time Distortion Loss

Our proposed framework for multi-step forecasting is depicted in Figure 2. During training, we consider a set of N input time series A = {x i } i∈{1:N } . For each input example of length n, i.e. x i = (x 1 i , ..., x n i ) ∈ R p×n , a forecasting model such as a neural network predicts the future k-step ahead trajectory ŷi = (ŷ 1 i , ..., ŷk i ) ∈ R d×k . Our Shape and Time Distortion Loss (STDL), which compares this prediction ŷi with the actual ground truth future trajectory *

y i = ( * y i 1 , ..., * y i k ) of length k, is composed of two terms balanced by the hyperparameter α ∈ [0, 1]: L ST DL (ŷ i , * y i ) = α L shape (ŷ i , * y i ) + (1 -α) L temporal (ŷ i , * y i ) (1)
Notations and definitions Both our shape L shape (ŷ i , * y i ) and temporal L temporal (ŷ i , * y i ) distortions terms are based on the alignment between predicted ŷi ∈ R d×k and ground truth * y i ∈ R d×k time series. We define a warping path as a binary matrix A ⊂ {0, 1}

k×k with A h,j = 1 if ŷh i is associated to * y i j , and 0 otherwise. The set of all valid warping paths connecting the endpoints (1, 1) to (k, k)

with the authorized moves →, ↓, be the pairwise cost matrix, where δ is a given dissimilarity between ŷh i and * y i j , e.g. the euclidean distance.

(step condition) is noted A k,k . Let ∆(ŷ i , * y i ) := [δ(ŷ h i ,

Shape and temporal terms

Shape term Our shape loss function is based on the Dynamic Time Warping (DTW) [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF], which corresponds to the following optimization problem: DT W (ŷ i , *

y i ) = min A∈A k,k A, ∆(ŷ i , * y i ) . A * = arg min A∈A k,k A, ∆(ŷ i , * y i )
is the optimal association (path) between ŷi and * y i . By temporally aligning the predicted ŷi and ground truth * y i time series, the DTW loss focuses on the structural shape dissimilarity between signals. The DTW, however, is known to be non-differentiable. We use the smooth min operator min γ (a 1 , ..., a n ) = -γ log( n i exp(-a i /γ)) with γ > 0 proposed in [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF] to define our differentiable shape term L shape :

L shape (ŷ i , * y i ) = DT W γ (ŷ i , * y i ) := -γ log   A∈A k,k exp   - A, ∆(ŷ i , * y i ) γ     (2)
Temporal term Our second term L temporal in Eq (1) aims at penalizing temporal distortions between ŷi and * y i . Our analysis is based on the optimal DTW path A * between ŷi and * y i . A * is used to register both time series when computing DTW and provide a time-distortion invariant loss. Here, we analyze the form of A * to compute the temporal distortions between ŷi and * y i . More precisely, our loss function is inspired from computing the Time Distortion Index (TDI) for temporal misalignment estimation [START_REF] Frías-Paredes | Assessing energy forecasting inaccuracy by simultaneously considering temporal and absolute errors[END_REF][START_REF] Vallance | Towards a standardized procedure to assess solar forecast accuracy: A new ramp and time alignment metric[END_REF], which basically consists in computing the deviation between the optimal DTW path A * and the first diagonal. We first rewrite a generalized TDI loss function with our notations:

T DI(ŷ i , * y i ) = A * , Ω = arg min A∈A k,k A, ∆(ŷ i , * y i ) , Ω (3)
where Ω is a square matrix of size k × k penalizing each element ŷh i being associated to an * y j i , for h = j. In our experiments we choose a squared penalization, e.g. Ω(h, j) = 1 k 2 (h -j) 2 , but other variants could be used. Note that prior knowledge can also be incorporated in the Ω matrix structure, e.g. to penalize more heavily late than early predictions (and vice versa).

The TDI loss function in Eq (3) is still non-differentiable. Here, we cannot directly use the same smoothing technique that for defining DTW γ in Eq (2), since the minimization involves two different quantities Ω and ∆. Since the optimal path A * is itself non-differentiable, we use the fact that A * = ∇ ∆ DT W (ŷ i , * y i ) to define a smooth approximation A * γ of the arg min operator, i.e. :

A * γ = ∇∆DT Wγ(ŷi, * y i ) = 1/Z A∈A k,k A exp - A,∆(ŷ i , * y i ) γ , with Z = A∈A k,k exp - A,∆(ŷ i , * y i ) γ
being the partition function. Based on A * γ , we obtain our smoothed temporal loss from Eq (3):

L temporal (ŷ i , * y i ) := A * γ , Ω = 1 Z A∈A k,k A, Ω exp - A,∆(ŷ i , * y i ) γ (4)

STDL Efficient Forward and Backward Implementation

The direct computation of our shape and temporal losses in Eq (2) and Eq (4) is intractable, due to the cardinal of A k,k , which exponentially grows with k. We provide a careful implementation of the forward and backward passes in order to make learning efficient.

Shape loss Regarding L shape , we rely on [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF] to efficiently compute the forward pass with a variant of the Bellmann dynamic programming approach [START_REF] Bellman | On the theory of dynamic programming[END_REF]. For the backward pass, we implement the recursion proposed in [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF] in a custom Pytorch loss. This implementation is much more efficient than relying on vanilla auto-differentiation, since it reuses intermediate results from the forward pass.

Temporal loss For L temporal , note that the bottleneck for the forward pass in Eq (4) is to compute A * γ = ∇ ∆ DT W γ (ŷ i , * y i ), which we implement as explained for the L shape backward pass.

Regarding L temporal backward pass, we need to compute the Hessian ∇ 2 DT W γ (ŷ i , * y i ). We use the method proposed in [START_REF] Mensch | Differentiable dynamic programming for structured prediction and attention[END_REF], based on a dynamic programming implementation that we embed in a custom Pytorch loss. Again, our back-prop implementation allows a significant speed-up compared to standard auto-differentiation (see section 4.4).

The resulting time complexity of both shape and temporal losses for forward and backward is O(k 2 ).

Discussion A variant of our approach to combine shape and temporal penalization would be to incorporate a temporal term inside our smooth L shape function in Eq (2), i.e. :

L ST DL t (ŷ i , * y i ) := -γ log   A∈A k,k exp   - A, α∆(ŷ i , * y i ) + (1 -α)Ω γ     (5)
We can notice that Eq (5) reduces to minimizing A, α∆(ŷ i , * y i ) + (1 -α)Ω when γ → 0 + . In this case, L ST DL t can recover DTW variants studied in the literature to bias the computation based on penalizing sequence misalignment, by designing specific Ω matrices: Sakoe-Chiba DTW hard band constraint [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] Ω(h, j) = +∞ if |h -j| > T , 0 otherwise

Weighted DTW [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF] Ω(h, j) = f (|i -j|), f increasing function L ST DL t in Eq (5) enables to train deep neural networks with a smooth loss combining shape and temporal criteria. However, L ST DL t presents limited capacities for disentangling the shape and temporal errors, since the optimal path is computed from both shape and temporal terms. In contrast, our L ST DL loss in Eq (1) separates the loss into two shape and temporal misalignment components, the temporal penalization being applied to the optimal unconstrained DTW path. We verify experimentally that our L ST DL loss outperforms its "tangled" version L ST DL t (section 4.3).

Experiments 4.1 Experimental setup

Datasets: To illustrate the relevance of our STDL approach, we carry out experiments on 3 nonstationary time series datasets from different domains (see examples in Fig 4). The multi-step evaluation consists in forecasting the future trajectory on k future time steps. Synthetic (k = 20) dataset consists in predicting sudden changes (step functions) based on an input signal composed of two peaks. This controlled setup was designed to measure precisely the shape and time errors of predictions. We generate 500 times series for train, 500 for validation and 500 for test, with 40 time steps: the first 20 are the inputs, the last 20 are the targets to forecast. In each series, the input range is composed of 2 peaks of random temporal position i 1 and i 2 and random amplitude j 1 and j 2 between 0 and 1, and the target range is composed of a step of amplitude j 2 -j 1 and stochastic position i 2 + (i 2 -i 1 ) + randint(-3, 3). All time series are corrupted by an additive gaussian white noise of variance 0.01.

ECG5000 (k = 56) dataset comes from the UCR Time Series Classification Archive [START_REF] Chen | The UCR time series classification archive[END_REF], and is composed of 5000 electrocardiograms (ECG) (500 for training, 4500 for testing) of length 140. We take the first 84 time steps (60 %) as input and predict the last 56 steps (40 %) of each time series (same setup as in [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF]).

Traffic (k = 24) dataset corresponds to road occupancy rates (between 0 and 1) from the California Department of Transportation (48 months from 2015-2016) measured every 1h. We work on the first univariate series of length 17544 (with the same 60/20/20 train/valid/test split as in [START_REF] Lai | Modeling long-and shortterm temporal patterns with deep neural networks[END_REF]), and we train models to predict the 24 future points given the past 168 points (past week).

Network architectures and training:

We perform multi-step forecasting with two kinds of neural network architectures: a fully connected network (1 layer of 128 neurons), which does not make any assumption on data structure, and a more specialized Seq2Seq model [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] with Gated Recurrent Units (GRU) [START_REF] Cho | Learning phrase representations using RNN encoderdecoder for statistical machine translation[END_REF] with 1 layer of 128 units. Each model is trained with PyTorch for a max number of 1000 epochs with Early Stopping with the ADAM optimizer. The smoothing parameter γ of DTW and TDI is set to 10 -2 . The hyperparameter α balancing L shape and L temporal is determined on a validation set to get comparable DTW shape performance than the DT W γ trained model: α = 0.5 for Synthetic and ECG5000, and 0.8 for Traffic. Our code implementing the STDL is available on line from https://github.com/vincent-leguen/STDL.

STDL forecasting performances

We evaluate the performances of our STDL loss function, and compare it against two strong baselines: the widely used Euclidean (MSE) loss, and the smooth DTW introduced in [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF][START_REF] Mensch | Differentiable dynamic programming for structured prediction and attention[END_REF]. For each experiment, we use the same neural network architecture (section 4.1), in order to isolate the impact of the training loss and to enable fair comparisons. The results are evaluated using three metrics: MSE, DTW (shape) and TDI (temporal). We perform a Student t-test with significance level 0.05 to highlight the best(s) method(s) in each experiment (averaged over 10 runs).

Overall results are presented in Table 1. MSE comparison: STDL outperforms MSE when evaluated on shape (DTW) in all experiments, with significant differences on 5/6 experiments. When evaluated on time (TDI), STDL also performs better in all experiments (significant differences on 3/6 tests). Finally, STDL is equivalent to MSE when evaluated on MSE on 3/6 experiments.

DTW γ [START_REF] Cuturi | Soft-dtw: a differentiable loss function for time-series[END_REF][START_REF] Mensch | Differentiable dynamic programming for structured prediction and attention[END_REF] comparison: When evaluated on shape (DTW), SDTL performs similarly to DTW γ (2 significant improvements, 1 significant drop and 3 equivalent performances). Regarding time (TDI) and MSE evaluations, STDL is significantly better than DTW γ in all experiments, as expected.

We display a few qualitative examples for Synthetic, ECG5000 and Traffic datasets on Fig 4 (other examples are provided in supplementary 2). We see that MSE training leads to predictions that are non-sharp, making them inadequate in presence of drops or sharp spikes. DTW γ leads to very sharp predictions in shape, but with a possibly large temporal misalignment. In contrast, our STDL predicts series that have both a correct shape and precise temporal localization. Evaluation with external metrics To consolidate the good behaviour of our loss function seen in Table 1, we extend the comparison using two additional (non differentiable) metrics for assessing shape and time. For shape, we compute the ramp score [START_REF] Vallance | Towards a standardized procedure to assess solar forecast accuracy: A new ramp and time alignment metric[END_REF]. For time, we perform change point detection on both series and compute the Hausdorff measure between the sets of detected change points T * (in the target signal) and T (in the predicted signal):

Hausdorff(T * , T ) := max(max t∈ T min t * ∈T * | t -t * |, max t * ∈T * min t∈ T | t -t * |) (6)
We provide more details about these external metrics in supplementary 1.1.

In Table 2

Comparison to temporally constrained versions of DTW

In Table 3, we compare the Seq2Seq STDL to its tangled variants Weighted DTW (STDL t -W) [START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF] and Band Constraint (STDL t -BC) [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] on the Synthetic dataset. We observe that STDL performances are similar in shape for both the DTW and ramp metrics and better in time than both variants. This shows that our STDL leads a finer disentanglement of shape and time components. Results for ECG5000 and Traffic are consistent and given in supplementary 3. We also analyze the gradient of STDL vs STDL t -W in supplementary 3, showing that STDL t -W gradients are smaller at low temporal shifts, certainly explaining the superiority of our approach when evaluated with temporal metrics. Qualitative predictions are also provided in supplementary 3.

Eval loss STDL (ours) STDL

STDL Analysis

Custom backward implementation speedup: We compare in Fig 5(a) the computational time between the standard Pytorch auto-differentiation mechanism and our custom backward pass implementation (section 3.2). We plot the speedup of our implementation with respect to the prediction length k (averaged over 10 random target/prediction tuples). We notice the increasing speedup with respect to k: speedup of × 20 for 20 steps ahead and up to × 35 for 100 steps ahead predictions.

Impact of α (Fig 5(b))

: When α = 1, L ST DL reduces to DTW γ , with a good shape but large temporal error. When α -→ 0, we only minimize L temporal without any shape constraint. Both MSE and shape errors explode in this case, illustrating the fact that L temporal is only meaningful in conjunction with L shape . 4 for the traffic dataset reveal the superiority of TT-RNN over LSTNet-rec, which shows that dedicated multi-step prediction approaches are better suited for this task. More importantly, we can observe that our Seq2Seq STDL outperforms TT-RNN in all shape and time metrics, although it is inferior on MSE. This highlights the relevance of our STDL loss function, which enables to reach better performances with simpler architectures.

Eval loss

LSTNet-rec [START_REF] Lai | Modeling long-and shortterm temporal patterns with deep neural networks[END_REF] TT-RNN [START_REF] Yu | Long-term forecasting using tensor-train RNNs[END_REF][START_REF] Yu | Learning chaotic dynamics using tensor recurrent neural networks[END_REF]

Conclusion and future work

In this paper, we have introduced the Shape and Time Distortion Loss (STDL), a new differentiable loss function for training deep multi-step time series forecasting models. The STDL combines two terms for precise shape and temporal localization of non-stationary signals with sudden changes. We showed that the STDL is comparable to the standard MSE loss when evaluated on MSE, and far better when evaluated on several shape and timing metrics. STDL compares favourably on shape and timing to state-of-the-art forecasting algorithms trained with the MSE.

For future work we intend to explore the extension of these ideas to probabilistic forecasting, for example by using bayesian deep learning [START_REF] Gal | Dropout as a bayesian approximation: Representing model uncertainty in deep learning[END_REF] to compute the predictive distribution of trajectories, or by embedding the STDL loss in a deep State Space Model architecture suited for probabilistic forecasting. Another interesting direction is to adapt our training scheme to relaxed supervision contexts, e.g. semi-supervised [START_REF] Robert | Hybridnet: Classification and reconstruction cooperation for semi-supervised learning[END_REF] or weakly supervised [START_REF] Durand | Exploiting negative evidence for deep latent structured models[END_REF], in order to perform full trajectory forecasting using only categorical labels at training time (e.g. presence or absence of change points).

Aknowledgements We would like to thank Stéphanie Dubost, Bruno Charbonnier, Christophe Chaussin, Loïc Vallance, Lorenzo Audibert, Nicolas Paul and our anonymous reviewers for their useful feedback and discussions.

Figure 1 :

 1 Figure 1: Limitation of the euclidean (MSE) loss: when predicting a sudden change (target blue step function), the 3 predictions (a), (b) and (c) have similar MSE but very different forecasting skills. In contrast, the STDL loss proposed in this work, which disentangles shape and temporal decay terms, supports predictions (b) and (c) over prediction (a) that does not capture the sharp change of regime.

Fig 1 .

 1 Here, the target ground truth prediction is a step function (in blue), and we present three predictions, shown in Fig 1(a), (b), and (c), which have a similar MSE loss compared to the target, but very different forecasting skills. Prediction (a) is not adequate for regulation purposes since it doesn't capture the sharp drop to come. Predictions (b) and (c) much better reflect the change of regime since the sharp drop is indeed anticipated, although with a slight delay (b) or with a slight inaccurate amplitude (c). This paper introduces the Shape and Time Distortion Loss (STDL), a new objective function for training deep neural networks in the context of multi-step and non-stationary time series forecasting. STDL explicitly disentangles into two terms the penalization related to the shape and the temporal localization errors of change detection (section 3). The behaviour of STDL is shown in Fig 1: whereas the values of our proposed shape and temporal losses are large in Fig 1(a), the shape (resp. temporal) term is small in Fig 1(b) (resp. Fig 1(c)). STDL combines shape and temporal terms, and is consequently able to output a much smaller loss for predictions (b) and (c) than for (a), as expected.

Figure 2 :

 2 Figure 2: Our proposed framework for training deep forecasting models.

Figure 3 :

 3 Figure 3: Shape and Time Distortion Loss computation for separating the shape and temporal errors.

Figure 4 :

 4 Figure 4: Qualitative forecasting results.

Figure 5 (

 5 Figure 5(a): Speedup of STDL

Table 1 :

 1 Forecasting results evaluated with MSE, Shape and Time metrics, averaged over 10 runs (mean ± standard deviation). For each experiment, best method(s) (Student t-test) in bold.

			Fully connected network (MLP)	Recurrent neural network (Seq2Seq)
	Dataset Eval	MSE	DTWγ [13]	STDL (ours)	MSE	DTWγ [13]	STDL (ours)
		MSE	1.65 ± 0.14	4.82 ± 0.40	1.67± 0.184	1.10 ± 0.17	2.31 ± 0.45	1.21 ± 0.13
	Synth	DTW	38.6 ± 1.28	27.3 ± 1.37	32.1 ± 5.33	24.6 ± 1.20	22.7 ± 3.55	23.1 ± 2.44
		TDI	15.3 ± 1.39	26.9 ± 4.16	13.8 ± 0.712	17.2 ± 1.22	20.0 ± 3.72	14.8 ± 1.29
		MSE	31.5 ± 1.39	70.9 ± 37.2	37.2 ± 3.59	21.2 ± 2.24	75.1 ± 6.30	30.3 ± 4.10
	ECG	DTW	19.5 ± 0.159	18.4 ± 0.749	17.7 ± 0.427	17.8 ± 1.62	17.1 ± 0.650	16.1 ± 0.156
		TDI	7.58 ± 0.192	38.9 ± 8.76	7.21 ± 0.886	8.27 ± 1.03)	27.2 ± 11.1	6.59 ± 0.786
		MSE	0.620 ± 0.010	2.52 ± 0.230	1.93 ± 0.080	0.890 ± 0.11	2.22 ± 0.26	1.00 ± 0.260
	Traffic	DTW	24.6 ± 0.180	23.4 ± 5.40	23.1 ± 0.41	24.6 ± 1.85	22.6 ± 1.34	23.0 ± 1.62
		TDI	16.8 ± 0.799	27.4 ± 5.01	16.7 ± 0.508	15.4 ± 2.25	22.3 ± 3.66	14.4± 1.58

Table 2 :

 2 , we report the comparison between Seq2Seq models trained with STDL, DTW γ and MSE. We see that STDL is always better than MSE in shape (Ramp score) and equivalent to DTW γ in 2/3 experiments. In time (Hausdorff metric), STDL is always better or equivalent compared to MSE (and always better than DTW γ , as expected). Forecasting results of Seq2Seq evaluated with Hausdorff and Ramp Score, averaged over 10 runs (mean ± standard deviation). For each experiment, best method(s) (Student t-test) in bold.

			MSE	DT W γ [13]	STDL (ours)
		Hausdorff	2.87 ± 0.127	3.45 ± 0.318	2.70 ± 0.166
	Synthetic Ramp score (x10)	5.80 ± 0.104	4.27 ± 0.800	4.99 ± 0.460
		Hausdorff	4.32 ± 0.505	6.16 ± 0.854	4.23 ± 0.414
	ECG5000 Ramp score	4.84 ± 0.240	4.79 ± 0.365	4.80 ± 0.249
		Hausdorff	2.16 ± 0.378	2.29 ± 0.329	2.13 ± 0.514
	Traffic	Ramp score (x10) 6.29 ± 0.319	5.78 ± 0.404	5.93 ± 0.235

Table 3 :

 3 Comparison to the tangled variants of STDL for the Seq2Seq model on the Synthetic dataset, averaged over 10 runs (mean ± standard deviation).

	43]

t -W [28] STDL t -BC [

Table 4 :

 4 Comparison with state-of-the-art forecasting architectures trained with MSE on Traffic, averaged over 10 runs (mean ± standard deviation).

	Seq2Seq STDL

A Seq2Seq architecture was the winner of a

Kaggle competition on multi-step time series forecasting (https://www.kaggle.com/c/web-traffic-time-series-forecasting)

We use the available Github code for both methods.


Xu-Guang Li
email: li@lss.supelec.fr

Arben C ¸ela
email: a.cela@esiee.fr

Silviu-Iulian Niculescu
email: silviu.niculescu@lss.supelec.fr

Abdellatif Reama
email: a.reama@esiee.fr

A Switched Control Method for Networked Control Systems

Keywords: Networked control systems (NCSs), switched systems, optimal control

The performance of networked control systems based on a switched control law is studied. The control law is switched between sampled-data feedback and zero control, the switch instant of which is determined by a switching parameter. The objective of the paper is to analyze 1) under what condition the switched control law can lead to better performance than the standard one and 2) how to select the switching parameter if the condition is satisfied. Numerical examples show the effectiveness of the proposed results.

INTRODUCTION

Employing a control network will bring various advantages, such as increased system reliability, ease of diagnosis and maintenance. The Networked Control System (NCS) has been an active research topic during the past several decades. A detailed introduction to NCSs can be found in [START_REF] Hespanha | A survey of recent results in networked control systems[END_REF], Hristu-Varsakelis and Levin (2005), and [START_REF] Zhang | Stability of networked control systems[END_REF].

In NCSs, zero-order-hold (ZOH) devices are widely used since the sampling rate can not be infinitely large. As a result, between two consecutive sampling instants, kT ≤ t < (k + 1)T , the controller can only use the sampled data x(kT ). Thus, a linear feedback control law u(t) = Kx(kT ) for kT ≤ t < (k + 1)T . This kind of control law, which will be called the standard control law in this paper, can be seen almost everywhere in the studies of NCSs, see e.g. [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF], [START_REF] Mirkin | Some remarks on the use of timevarying delay to model sample-and-hold circuits[END_REF], [START_REF] Fujioka | Stability analysis for a class of networked/embedded control systems: A discrete-time approach[END_REF].

In this paper, we will employ an easily-implemented switched control law instead of the standard one. The basic idea is that when kT ≤ t < kT + a, let u(t) = Kx(kT ); while when kT + a ≤ t < (k + 1)T , u(t) = 0, where a (0 ≤ a ≤ T ) is the switching parameter to be designed.

When using the standard control law, the closed-loop system is equivalent to an input-delayed system with delay τ (t) = t-kT , as analyzed in [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF]. Due to the effect of delay, we believe that in a certain time interval the open-loop dynamic (i.e. u(t) = 0) is more favorable for the system to "converge". This is the reason we employ the new control law. Here, we give a motivating example. Consider the system with matrices below equation ( 6) and with sampling T = 1.7. The system is asymptotically stable and has a Lyapunov function V (t) = x T (t)P x(t) with P = [4.03 5.09; 5.09 13.49]. The trajectory of V (t) is shown in Fig. 1. It is seen that at sampling instants, V (t) monotonously decreases. However, between sampling instants, V (t) first decreases and then increases. For this system, as time increases between two sampling instant, the standard control law gradually loses its stabilizing effect. We then employ a switched control with switching parameter a = 0.94 and make a simulation for V (t), see the bold curve in Fig. 1. It is clear that the switched control law lead better performance than the standard one.

The switched control law has been used in [START_REF] Sun | A novel method of stability analysis for networked control systems[END_REF], where it is called even-time-driven control. Therein, a Lyapunov-Krasovskii method is used to estimating the decay rate of the closed-loop system and diverge rate of the open-loop system. The stability is guaranteed if the system is on average convergent over the whole sampling period. The method proposed in this paper is different with that used in [START_REF] Sun | A novel method of stability analysis for networked control systems[END_REF].

In this paper, we will address the following problems: Under what condition, the standard control law is better and under what condition, the new control law is better? (Through the motivating examples given later, we will see that the new control law does not always lead to better performance) If the new control law can lead to better performance, how should we choose a switching time? Optimal control for networked control systems has been widely considered, see e.g., Ben Gaid, C ¸ela, and Hamam (2006) and [START_REF] Chen | H 2 -optimal sampleddata control[END_REF]. In this paper, we will use a different method, with the existing ones, to design a switching parameter. x T (iT )Qx(iT ), we can have dV /da = x T (0)Θ (a)x(0), as will be shown later in the text. Thus, in an interval where Θ (a) is negative-definite (positive-definite), the system performance monotonously decreases (increases) with respect to the switching parameter a. If in an interval where Θ (a) is indefinite (neither negative-semidefinite or positivesemidefinite), there does not exist an optimal switching parameter. In this case, the optimal switching parameter varies with respect to different initial condition. According to the property (negative-definite, positive-definite or indefinite) of Θ (a), we can tell if the switched control law can lead to better performance than the standard one and how to choose the switching parameter. Numerical examples show the effectiveness of the proposed method.

The remaining paper is organized as follows. In Section 2, the optimization problem of a networked control system with the switched control is formulated and some motivating examples are given. In Section 3, main results of the paper are proposed. Numerical examples are presented in Section 4, and some concluding remarks end the paper.

PROBLEM FORMULATION AND MOTIVATING EXAMPLES

Consider a networked control system with a switched control law

ẋ(t) = Ax(t) + Bu(t) (1) u(t) = Kx(kT ), kT ≤ t < kT + a 0, kT + a ≤ t < (k + 1)T (2)
where A, B, K are n × n, n × m, m × n real matrices, T is the sampling period, k are nonnegative integers and a is the switching parameter satisfying 0 ≤ a ≤ T .

It is assumed that matrix A is not Hurwitz and K is designed such that A+BK is Hurwitz, without considering the effect of network.

The standard control law is given by Under the standard control law (3), which is implemented by the zero-order-hold (ZOH) devices, the control value is fixed as Kx(kT ) until the new sampled data x((k + 1)T ) arrives.

u(t) = Kx(kT ), kT ≤ t < (k + 1)T (3)
The controller (2) used in this paper has two switched control status: between kT ≤ t < kT + a, it uses the value Kx(kT ), as the standard controller does; however between kT + a ≤ t < (k + 1)T , it uses the value 0. System (1) under controller (2) thus can be expressed in the following switched system:

ẋ(t) = Ax(t) + BKx(kT ), kT ≤ t < kT + a Ax(t), kT + a ≤ t < (k + 1)T (4)
We choose control law (2) instead of the standard one (3) based on the following analysis:

When employing the standard control law (3), the networked control system can be equivalently written as

ẋ(t) = Ax(t) + BKx(t -τ (t)), kT ≤ t < (k + 1)T (5) where τ (t) = t -kT, kT ≤ t < (k + 1)T .
Between two sampling instants, the delay τ (t) increases with respect to time. Therefore, it is possible that in some time interval the open-loop dynamic ẋ(t) = Ax(t) is more favorable than the closed-loop dynamic ẋ(t) = Ax(t) + BKx(kT ) for the system performance, due to the effect of delay. If this is true, after an appropriate switching instant kT + a, we set the control values as 0 (which is easy to implement in practice) instead of the sampleddata feedback such that an improved performance can be obtained.

In this paper, we choose the following performance index:

V = N i=1 x T (iT )Qx(iT ) (6)
And the system matrices are given by -3.75 -11.5 ) , Q = 4.03 5.09 5.09 13.49 .

A = 0 1 0 -0.1 , B = 0 0.1 , K = (
Two motivating examples are given before the end of this section.

Motivating example 1: When T = 0.9 and N = 2. We choose four switching parameters 0.3, 0.5, 0.7, 0.9 and three initial state (1, 1) T , (1, -1) T , (1, 2) T . We compute the performance under different switching parameters and initial conditions, see Table .1 for the results.

From the table, we can see that the optimal performance is always obtained when a = T = 0.9. That is to say, the control law (2) can not lead to better performance than the standard control law (3), for this example.

Motivating example 2: When T = 1.2 and N = 100. We choose four switching parameters 0.6, 0.9, 1.0, 1.2 and three initial state (1, 1) T , (1, -1) T , (1, 2) T . The computed performance are listed in Table .2.

From the table, we can see that if switched at some but not arbitrary instant, the control law (2) can lead to better performance than the standard control law (3). Now, one may naturally have two questions on the switched control law (2): 1) Under what condition, we can obtain better performance by using the switched control law (2) than using the standard one (3)? How should we choose an appropriate switching parameter to obtain an improved performance? The objective of this paper is to answer these two questions. Two emphasize the effect of switching parameter, we do not take the transmission delay into account. We only consider a simple form of performance index here. The method of this paper can be further extended to the cases with transmission delay and with other performance indices.

The stability analysis, not the purpose of this paper, for NCSs based on the switched control (2) is reported in [START_REF] Li | Stability analysis of networked control systems based on a switched control[END_REF]. It is shown that larger stability interval can be obtained by employing the switched control (2). For the system considered in this paper, it the sampling period is constant, the stability interval is T ∈ [0, 1.724) by the standard control (3). By the switched control (2), the stability interval is T ∈ [a, +∞), if a ∈ (0, 0.534). If the sampling periods are time-varying, the stability interval is T k ∈ [0, 1.72] by the standard control (3). By the switched control (2), the stability interval is T k ∈ [0, 10.25].

MAIN RESULTS

The system described by ( 1) and ( 2), with the initial state x(0) has the following response:

x(t) = Φ(t -kT, t -kT )Φ k (a, T )x(0), kT ≤ t < kT + a Φ(a, t -kT )Φ k (a, T )x(0), kT + a ≤ t < (k + 1)T where Φ(a, T ) = e A(T -a) (e Aa + a 0 e Aθ dθBK)
And the discrete-time expression of the system is

x((k + 1)T ) = Φ(a, T )x(kT )
Remark 1. It is easy to see that when a = T , the switched control law (2) is equivalent to the standard one (3).

Given a sampling period T , the performance index ( 6) can be rewritten as

V (a) = x T (0)Θ(a)x(0) (7)
where

Θ(a) = N i=1 (Φ i (a, T )) T QΦ i (a, T )
It can be seen that for given x(0), T and N , the system performance depends on the switching parameter a. Some results are given below to show the relationship between the performance V (a) and the switching parameter a. Throughout this paper, A > 0 (A < 0) denotes that matrix A is positive-definite (negative-definite). Theorem 1. If there exists an interval [a -, a -] such that

Θ (a) < 0 for a ∈ [a -, a -]
where

Θ (a) = N i=1 (d Φ i (a, T ) T da QΦ i (a, T ) + Φ i (a, T ) T Qd Φ i (a, T ) da ), d Φ(a, T ) da = (e AT -Ae A(T -a)
a 0 e Aθ dθ)BK, it follows that V (a 1 ) > V (a 2 ) for any a 1 , a 2 with a -≤ a 1 < a 2 ≤ a -and any initial state x(0).

If there exists an interval [a + , a + ] such that Θ (a) > 0 for a ∈ [a + , a + ], it follows that V (a 1 ) < V (a 2 ) for any a 1 , a 2 with a + ≤ a 1 < a 2 ≤ a + and any initial state x(0).

Proof. One can see that d V (a) da = x T (0)Θ (a)x(0). For any initial state x(0), we have V (a) decreases (increases) as a increases if Θ (a) < 0(Θ (a) > 0). The proof is complete. Corollary 1. If Θ (a) < 0 for a ∈ [0, T ], the control law (2) can not lead to better performance than the control law (3).

Proof. If Θ (a) < 0 for a ∈ [0, T ], we always have the minimal performance when a = T , which is exactly the standard control law (3). Corollary 2. If there exists an interval [T 1 , T ](0 ≤ T 1 < T ) such that Θ (a) > 0 for a ∈ [T 1 , T ], the control law (2) can lead to better performance than the control law (3) under any initial state.

Proof. One can see that there exist a ∈ [T 1 , T ) such that V (a) < V (T ) under any initial state, if Θ (a) > 0 for a ∈ [T 1 , T ].
Theorem 2. If there exists an interval [a 0 , a 0 ] such that Θ (a) is indefinite (neither negative-semidefinite nor positivesemidefinite) for a ∈ [a 0 , a 0 ], there does not exist an optimal switching parameter, say a * , in this interval such that V (a * ) < V (a) for a = a * ∈ [a 0 , a 0 ] under any initial condition.

Proof. Assume that there is an optimal switching parameter a * . For a positive a # , which is sufficiently close to a * , we have Θ(a * ) -Θ(a # ) is indefinite since in [a 0 , a 0 ], Θ (a) is indefinite. Therefore, there exists an initial state x(0) satisfying V (a * ) -V (a # ) > 0. The result is proved by contradiction. Corollary 3. If there exists an interval a ∈ [T 2 , T ](0 ≤ T 2 < T ) where Θ (a) is indefinite, the control law (2) can lead to better performance than the control law (3) under some (not all) initial state. This result follows straightforwardly from Theorem 2. Remark 2. Sometimes, the optimal switching parameter does not exist. However, we can have an indefinite interval such that for any initial state the optimal switching parameter lie in this interval. We call this interval quasioptimal interval in this paper. This will be seen in the numerical example section

The eigenvalues of Θ (a) are the roots of the following characteristic equation

CE = det(λI -Θ (a))
The highest degree of CE is fixed. Thus, as analyzed in Li, C ¸ela, Niculescu, and Reama (2009a), the eigenvalues of Θ (a) are continuous with respect to a and the parameter sweeping technique is effective to check the property (negative-definite, positive-definite or indefinite) of Θ (a). The parameter sweeping technique can be easily realized by standard computing software such as MATLAB. It is worth to note that the parameter sweeping technique has been widely used in the stability analysis, see, e.g., [START_REF] Fazelinia | Stability robustness analysis of multiple time-delayed systems using "Building Block" concept[END_REF] and [START_REF] Gu | On stability of crossing curves for general systems with two delays[END_REF].

NUMERICAL EXAMPLES

In this section, some numerical examples are given to illustrate the proposed method in this paper. We choose the system (1) under the switched control law (2) with matrices A, B, K, Q given in section 2.

Example 1: When T = 0.9 and N = 2. By parametersweeping, we find for a ∈ [0, 0.9], Θ (a) < 0. According to Theorem 1, we know the performance monotonously decreases with respect a and the optimal performance is always obtained at a = T = 0.9. In other words, the control law (2) can not lead to better performance than the standard control law (3), according to Corollary 1. We choose three initial conditions (1, 1) T , (1, -1) T , and (1, 2) T for simulation, see Fig. 2. The simulations coincide with the proposed results in this paper.

Example 2: When T = 1.2 and N = 100.

By parameter-sweeping, we have Θ (a) < 0 for a ∈ [0, 0.810) and Θ (a) > 0 for a ∈ (1.003, 1.2]. Therefore, the quasi-optimal switching parameters should lie in the interval [0.810, 1.003](the indefinite interval). For any initial state, the optimal performance will be obtained in the interval [0.810, 1.003].

For this example, we choose three initial state (1, 1) T , (-2, 1) T , and (3, -1) T . And for each simulation, a closeup graph is added, see Figures 234567. For initial state (1, 1) T , (-2, 1) T , and (3, -1) T , we have the optimal switching parameter 0.99, 0.88, and 0.84, respectively. Through the simulation, we can see the performance monotonously decreases as a increases in the interval [0, 0.810), while monotonously increases as a increases in the interval (1.003, 1.2]. The optimal switching parameter varies depending on the initial state, however it must be inside [0.810, 1.003]. For this example, we can obtain improved performance when employing the switched control law (2) than employing the standard control law (3).

CONCLUSION

The performance of networked control systems based on a switched control law is considered. The relationship between the switching parameter and the system performance is analyzed. The method is based on testing the property (negative, positive or in-definite) of the performance derivative with respect to the switching parameter. As a result, we know if the switched control law can lead to better performance than the standard one and how to select a switching parameter.

Ben Gaid, M.E.M., C ¸ela, A., and Hamam, Y., (2006).

Optimal integrated control and scheduling of networked control systems with communication constraints: appli-

  Fig. 1. A motivating exampleThe analysis in this paper is based on testing the derivative of the performance with respect to the switching parameter. Given a performance index V = N i=1

Fig. 4 .

 4 Fig. 2. Simulation for example 1

  Fig. 5. Simulation for example 2 with x 0 = (-2, 1) T

Fig. 6 .

 6 Fig. 6. Simulation for example 2 with x 0 = (-2, 1) T (closeup)

Fig. 7 .

 7 Fig. 7. Simulation for example 2 with x 0 = (3, -1) T

Fig. 8 .

 8 Fig. 8. Simulation for example 2 with x 0 = (3, -1) T (closeup)

Table 1 .

 1 Performance for motivating example 1

		(1, 1) T	(1, -1) T	(1, 2) T
	0.3	34.5178	1.3156	16.0483
	0.5	16.1623	1.0396	13.4389
	0.7	8.0496	0.8277	11.2502
	0.9	5.2061	0.6687	9.4186

Table 2 .

 2 Performance for motivating example 2

		(1, 1) T	(1, -1) T	(1, 2) T
	0.6	12.8753	2.5601	33.7183
	0.9	4.4923	0.5462	8.3129
	1.0	4.1121	0.5059	7.2336
	1.2	6.0208	0.9400	12.9831


Luca De Cicco
email: ldecicco@poliba.it

Saverio Mascolo
email: mascolo@poliba.it

Silviu-Iulian Niculescu
email: silviu.niculescu@lss.supelec.fr

Robust Stability Analysis of a Class of Smith Predictor-based Congestion Control Algorithms for Computer Networks

Keywords: Time-delay systems, Robust stability, congestion control, Smith predictor

come

INTRODUCTION

Time delays are often present in feedback control systems due to reasons such as the transport of material or information. From the control theoretic point of view it is well-known that an increase of the time delay may lead to instability of the closed loop system and to performance degradation as well.

The Smith principle is a classic approach which is often employed to design controllers for time delay systems [START_REF] Smith | A controller to overcome dead time[END_REF]). It is known that, by assuming exact knowledge of both the plant model and time delay, controllers designed using a Smith predictor are very eective in counteracting the eect of time delays. Robustness of the Smith predictor with respect to uncertainties in the knowledge of the time delay has been extensively studied since 1980 [START_REF] Palmor | Stability properties of Smith deadtime compensator controllers[END_REF], [START_REF] Yamanaka | Eects of mismatched Smith controller on stability in systems with time-delay[END_REF]).

The Internet represents a relevant example of time delay system due to the presence of delays that are caused by the propagation of the information, which is sent in form of data packets, from a source to a destination through a series of communication links and router queues.

A cornerstone component of the Internet protocol stack is the end-to-end congestion control which has been implemented in the TCP by [START_REF] Jacobson | Congestion avoidance and control[END_REF] in order to avoid congestion and preserve network stability. Several uid models have been proposed for the TCP congestion control algorithm in order to analytically study the stability of the network under dierent scenarios. See [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF]; [START_REF] Hollot | Analysis and design of controllers for AQM routers supporting TCPows. Automatic Control[END_REF]; [START_REF] Srikant | The Mathematics of Internet Congestion Control[END_REF] and references therein.

In [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF] a simple model of the plant made of an integrator (modelling the bottleneck queue) plus two time delays (modelling forward and backward delays), has been proposed along with a Smith predictor plus a proportional controller. The paper also shows that the Smith predictor controller with a proportional gain models the congestion control law which is employed in the today running TCP congestion control algorithm. Moreover, the model presented in [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF] has been employed in [START_REF] Grieco | Adaptive rate control for streaming ows over the internet[END_REF] to design and implement a ratebased congestion control algorithm which has been found to produce ows that are TCP-friendly.

A similar plant model is employed in [START_REF] Quet | Rate-based ow controllers for communication networks in the presence of uncertain time-varying multiple time-delays[END_REF] to design a rate-based congestion control algorithm implemented at the router that is robust to uncertain timedelays by employing the H ∞ technique. The controller parameters are quite complex to derive and, at the best of authors knowledge, no real implementations are currently available to assess the eectiveness of the proposed solution.

Measurement of the plant time delay to be used in the Smith predictor can be aected by uncertainties due to the fact that the time delay is made of a constant propagation delay plus time-varying queueing delays. To the purpose, the standard TCP [START_REF] Postel | Transmission Control Protocol[END_REF]) estimates the Round Trip Time (RTT) through time-stamping in order to set the retransmission timeout (RTO) which is needed for detecting heavy congestion episodes in the network. The RTT is dened as the time that elapses from when a segment is sent until the corresponding acknowledgement segment is received by the sender. In the standard TCP implementation, the RTT is measured each RT T seconds, whereas no measurements are taken on retransmitted segments due to the Karn's algorithm in order to avoid spurious timeouts (see [START_REF] Karn | Estimating round-trip times in reliable transport protocols[END_REF]). For these reasons the standard TCP does not provide an accurate measure of RTT. In order to overcome this issue an optional scheme has been proposed and standardized in [START_REF] Jacobson | TCP Extensions for High Performance[END_REF] which makes use of timestamps in an optional eld of the TCP header.

However, even if the timestamp option is employed by both peers of the communication, the granularity chosen for TCP timestamps is implementation-dependent. In a recent work [START_REF] Veal | New Methods for Passive Estimation of TCP Round-Trip Times[END_REF] carry out an extensive measurement campaign on RTTs. Authors used 500 servers and found that 76% of the servers had timestamping option enabled, and out of these servers 37% used a 100 ms granularity, 55% a 10 ms granularity and only 7% of them had a granularity of 1 ms.

A preliminary study on robust stability of a proportional Smith predictor used for congestion control in data networks has been carried out by using the Nyquist criterion in [START_REF] Mascolo | Modeling the Internet congestion control as a time delay system: a robust stability analysis[END_REF]. It revealed that in order to guarantee asymptotic stability it is sucient that ∆ < 1/k where ∆ represents the delay uncertainty and k is the gain of the proportional controller.

The goal of this paper is to provide a characterization of the robust stability of system introduced in Mascolo (1999) by applying the geometric approach which has been developed in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF].

The rest of the paper is organized as follows: in Section 2 we briey review the model of the closed loop congestion control in a generic packet switching network presented in [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF]; in Section 3 we apply the geometrical approach developed in [START_REF] Morarescu | On the geometry of stability regions of Smith predictors subject to delay uncertainty[END_REF] in order to nd the stability crossing curves of the system; in Section 4 we present the robust stability analysis; in Section 5 some simulations are presented to support the theoretical results obtained; nally Section 6 concludes the paper.

CONGESTION CONTROL MODEL

A network connection is basically made by a set of communication links and store-and-forward nodes (routers) where packets are enqueued before being routed to the destination . Congestion can arise when packets arrive at a rate r(t) which is above the capacity of the output link so that the router queue builds up until it is full and it starts to drop packets.

In [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF] a model of the Internet ow and congestion control as a time delay system is provided and in [START_REF] Mascolo | Modeling the Internet congestion control using a Smith controller with input shaping[END_REF] it is shown that dierent variants of TCP congestion control algorithms can be modelled in a unied framework by proper input shaping of the proportional Smith predictor controller.

In particular, the model consists of a feedback loop in which two time delays are present as it is shown in Figure 1: τ f w models the propagation time of a packet from source to the bottleneck queue and τ f b models the propagation time from the bottleneck to the destination and then back to the sender. The round trip time of the connection is

τ = τ f w + τ f b . k 1 s e -sτ f b e -sτ f w - - w(t) 1-e -sτ s q(t) - r(t) b(t) Figure 1. Functional block of the congestion control model
The simple integrator 1/s models the bottleneck queue that is lled (or drained) by the rate mismatch r(t) -b(t), where b(t) is the bottleneck available bandwidth.

The controller is a proportional Smith predictor with gain k that computes the rate r(t) to match the available bandwidth b(t) and to produce a stable output. The reason for using a simple proportional controller is that in this way the closed-loop dynamics can be made that of a rst-order system with time constant 1/k delayed by τ f w . Thus, the step response of the system can be made faster by increasing the proportional gain k providing an always stable system without oscillations or overshoots.

This choice provides a controller in which only one design parameter, i.e. the gain k, has to be tuned having a direct inuence on the dynamics of the output. Finally, the input signal w(t) models the congestion window (cwnd) or the advertised window, that is used by the congestion control algorithm to bound in-ight packets [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF]).

Model mismatches are known to aect the closed loop dynamics when a Smith predictor controller is employed.

In this case, it is worth noting that the only source of mismatch between the model and the actual plant is the entity of the delay (see Section 1) whereas the model of the bottleneck queue is an integrator and does not add any uncertainty. In the next sections we will give simple tuning rules for the design parameter k in order to retain asymptotic stability when the measure of time delay τ is uncertain.

Finally, we remark that a Smith predictor controller is recommended when designing a congestion control algorithm for data networks, since using PID controllers would provide an unacceptable sluggish system due to large delays involved in communication networks (see [START_REF] Astrom | PID Controllers: Theory, Design, and Tuning[END_REF] and [START_REF] Mascolo | Congestion control in high-speed communication networks using the Smith principle[END_REF]).

STABILITY CROSSING CURVES IN THE PARAMETERS SPACE

Review of the geometrical approach

We start by briey reviewing the geometrical approach developed in [START_REF] Morarescu | On the geometry of stability regions of Smith predictors subject to delay uncertainty[END_REF] which we will employ to analyze the robust stability of the considered system.

The reader is advised to refer to [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF] for a complete description of the method. We denote with a(s; τ 1 , τ 2 ) the characteristic function of the closed-loop system where τ 1 represents the nominal delay used in the Smith predictor and τ 2 = τ 1 + ∆ represents the actual plant delay aected by a bounded mismatch ∆. It is easy

q(ω) -h ( s ) e
a(s; τ 1 , τ 2 ) = 1 -h(s)e -τ1s + h(s)e -τ2s
(1) where h(s) is the transfer function of the closed loop system when no delays are present in the loop:

h(s) = C(s)G 0 (s) 1 + C(s)G 0 (s)
with G 0 (s) being the delay free plant and with C(s) being the controller transfer function.

In order to analyze the stability of the system we look for the solutions of the characteristic equation:

a(jω; τ 1 , τ 2 ) = 0 (2)
In this way we are able to nd all the conditions under which the system has at least one pole on the imaginary axis. The geometrical approach relies on the observation that the three terms of the characteristic function (1) can be seen as vectors in the complex plane. Therefore, the equality a(s; τ 1 , τ 2 ) = 0 can be represented in the complex plane via an isosceles triangle as it is shown in Figure 2. Thus, equation ( 2) is equivalent to the following three conditions:

(1) The triangular inequality must hold for the triangle shown in Figure 2, which implies that:

|h(jω)| ≥ 1 2 (3)
(2) Equation ( 2) must satisfy the phase rule;

(3) The sum of the internal angles of the isosceles triangle must be equal to π;

The solution of (3), which does not depend on time delays τ 1 or τ 2 , forms the frequency crossing set Ω which is the union of a nite number N of intervals of nite length Ω 1 , Ω 2 , . . . , Ω N . For any ω > 0 which belongs to the frequency crossing set there exists at least a pair (τ 1 , τ 2 ) in the parameters space such that the system has at least one imaginary pole. The conditions 2 and 3 imply that for all ω ∈ Ω all the couples (τ 1 , τ 2 ) ∈ R 2 + satisfying a(jω; τ 1 , τ 2 ) = 0 can be found using the following equations:

τ u± 1 = ∠h(jω) + 2uπ ± q(ω) ω (4) τ v± 2 = ∠h(jω) + (2v -1)π ∓ q(ω) ω (5)
where u and v are integers such that the corresponding τ u± 1 , τ v± 2 are non negative and q(ω) represents the internal angle of the isosceles triangle:

q(ω) = arccos 1 2|h(jω)|
In order to understand the meaning of equations ( 4) and

(5) let us x u = u,v = v and consider the set Ω i ⊆ Ω : if ω varies in Ω i and we evaluate (4) and ( 5) for both positive and negative signs we obtain two curves in the parameter space (τ 1 , τ 2 ) which we denote T i+ u,v and T i- u,v respectively.

It is worth noting that the curves

T i u,v = T i- u,v ∪ T i+ u,v can be
either open curves or closed curves depending on the set Ω i we are considering. In particular, it is easy to show that if the left end of Ω i is 0 then the associated curve is an open curve with both ends approaching ∞ when ω → 0.

On the other hand, if the left end of Ω i is not 0 then T i u,v is a closed curve [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]).

We dene the stability crossing curves T in the τ 1 , τ 2 plane as the union of all the curves T i u,v when i ∈ {1, . . . , N }, and u and v vary in the set of integers.

Finally, it is important to point out that when a stability crossing curve is crossed in the τ 1 , τ 2 plane, two poles cross the imaginary axis on the complex plane [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]).

Stability crossing curves of the computer network congestion control model

In order to characterize the impact of the delay uncertainty on the stability of the considered feedback system we apply the geometric approach we have reviewed in Section 3.1.

It is worth to notice that the delay-free model of the plant is G 0 (s) = 1/s and the controller transfer function is C(s) = k.

We suppose that the system described in Section 2 is aected by a delay uncertainty ∆ which is bounded by δ > 0, i.e. |∆| < δ. By considering the delay uncertainty, the characteristic equation of the system can be rewritten as follows:

1 + k s - k s e -τ1•s (1 -e -∆•s ) = 0 (6)
where τ 1 represents the nominal round trip time (RTT) of the considered connection, which is used in the Smith predictor, and τ 2 = τ 1 + ∆ is the actual plant time delay.

By multiplying by s/(s + k) both sides of ( 6) we obtain:

1 - k s + k e -τ1s + k s + k e -(∆+τ1)s = 0 (7)
so that by considering h(s) = k/(s + k) (7) is in the form of (1). We are interested in characterizing the stability of the system when τ 1 , τ 2 and k vary in R + .

By making the change of variable z = s/k we obtain:

1 - 1 z + 1 e -h1z + 1 z + 1 e -h2z = 0 (8)
where h 1 = kτ 1 and h 2 = kτ 2 , which reduces the free parameters to two. It is worth to notice that the 7) to (8) simply involves a scaling of the closed-loop eigenvalues by 1/k, thus indicating a natural trade-o between gain and delay since when k increases the closed loop poles approach to the imaginary axis [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays. Automatic Control[END_REF]).

We are now ready to study the stability of the original system in the h 1 , h 2 plane regardless the value of the proportional gain k.

First of all, by applying (3) we nd that the crossing set is made by the single interval Ω = [0, √ 3] which means that the stability crossing curves in the h 1 , h 2 plane are open curves which extend to innity when ω → 0. By using (4) and (5) the stability crossing curves of the considered system result the following:

h u± 1 (ω) = -arctan ω + 2uπ ± arccos √ 1+ω 2 2 ω (9) h v± 2 (ω) = -arctan ω + (2v -1)π ∓ arccos √ 1+ω 2 2 ω ( 10
)
Figure 3 shows the stability crossing curves of the considered system. We start by considering the h 2 axis (h 1 = 0), which means that we are employing a simple proportional controller without the Smith predictor. By starting from the origin and increasing the value of h 2 the rst curve is crossed at h 2 = π/2 which means that the system becomes unstable for h 2 > π/2 as expected. On the other hand, the axis h 1 represents the system in which no delay aects the plant, but the Smith predictor is in the controller. Figure 3 shows that the system is stable for all the delays in the Smith predictor.

Finally, points on the positive bisector represent the case of perfect matching of nominal delay τ 1 with the actual delay τ 2 . Indeed, if we move on this line no curves will be crossed since the Smith predictor in this case provides a stable system regardless the value of the proportional gain k.

ROBUST STABILITY ANALYSIS

In this Section we will develop an analysis of the robust stability of the considered system by using the stability crossing curves we have shown in the previous Section. We already know that the considered system is always asymptotically stable for any delay τ 1 and any proportional gain k as far as the delay uncertainty is zero thanks to the perfect compensation of the time delay τ 1 provided by the Smith predictor. In the h 1 , h 2 plane this condition means that the system is asymptotically stable on all the positive bisector.

In order to characterize the robustness of the system in the face of delay uncertainties we compute the maximum delay mismatch which still preserves stability. Thus, the problem here is to look for the maximum deviation δ with respect to a generic point (τ * 1 , τ * 1 ) with τ * 1 ≥ 0 which lies on the positive bisector such that the system is stable for any (τ 1 , τ 2 ) which satises:

|τ 2 -τ * 1 | < δ
We remark that solving the maximum admissible delay uncertainty problem is equivalent to nd the minimum distance between the stability crossing curves and a generic point on the positive bisector of the h 1 , h 2 plane.

Thus for any τ * 1 > 0 we have to solve:

δ(τ * 1 ) = min u,v min τ v± 2 ∈T |τ v± 2 -τ * 1 | (11)
so that the maximum delay to retain stability is:

δ = min τ * 1 ∈R+ δ(τ * 1 ) (12)
Proposition 1. A necessary and sucient condition for the asymptotic stability of the system regardless the value of the nominal delay τ 1 is:

|∆| < α k (13)
where ∆ is the delay uncertainty, α ∼ = 1.4775 and k is the proportional gain of the controller.

Proof.

We start by considering the stability crossing curves in the parameters space h 1 , h 2 . In order to nd the minimum distance between the stability crossing curves and a generic point of positive bisector of the h 1 , h 2 plane we evaluate the tangent to the crossing curves with direction parallel to the positive bisector:

dh 2 dh 1 = 1 ⇔ dh 2 dω dω dh 1 = 1 ⇔ dh 1 dω = dh 2 dω (14)
To the purpose we look for a subset T of the stability crossing curves T that are the closest curves to the positive bisector. By considering a generic curve T u,v and by applying ( 9) and ( 10) it turns out that for all u and v and for all ω ∈ Ω it holds h v+ 2 -h u+

1 < h v- 2 -h u-
1 so that it is sucient to consider only the curves T + u,v in the region h 2 > h 1 and the curves T - u,v in the region h 2 < h 1 , since they will be the closest ones to the positive bisector. Thus, we can refer without loss of generality to the generic curves of T as T u,u+i for all i and u in the integers. Straightforward computations on ( 9) and ( 10)

give:

h u+i 2 -h u 1 > h u+i-1 2 -h u 1
which means that when i decreases the curves T u,u+i will move downwards in the h 1 , h 2 plane. Figure 4 shows the values of u and v for the curves T u,u and T u,u+1 . It is then easy to show that if we set v = u we obtain the closest curves to the positive bisector in the region h 2 < h 1 whereas the curves with v = u + 1 are those which are closest to the positive bisector in the region h 2 > h 1 . In conclusion we can restrict our search to the set :

T = T - u,u ∪ T + u,u+1
for all u in the integers. Let us consider the region h 2 > h 1 i.e. we consider the subset T + u,u+1 . By considering (14) after straightforward computations we get the following equation:

arccos

√ ω 2 + 1 2 + ω 2 √ ω 2 + 1 √ 3 -ω 2 + +π(v -u - 1 2 ) = 0 (15)
with ω ∈ Ω. By letting v = u + 1 the equation ( 15) has the unique solution ω = 1.3483 rad/s in Ω which is independent of u. If we substitute this value in ( 9) and ( 10) we obtain:

h 1 (ω) = h 1 = 4.6601u -0.2654 h 2 (ω) = h 2 = 4.6601v -3.4480
Thus, all the points belonging to the the curves T + u,u+1

having a tangent which is parallel to the positive bisector, lie on the line:

h 2 = h 1 + 1.4775 (16)
For this reason we can conclude that the maximum uncertainty, in the h 1 , h 2 coordinates is 1.4775. The proof is completed by recalling that h 1 = kτ 1 and h 2 = kτ 2 and that τ 2 = τ 1 + ∆. Thus, we nally obtain:

h 2 -h 1 < 1.4775 ⇒ k∆ < 1.4775 ⇒ ∆ < 1.4775 k (17)
It is worth to notice that the same procedure can be followed in the case v = u which leads to the inequality:

h 1 -h 2 < 1.4775 ⇒ -k∆ > 1.4775 ⇒ ∆ > - 1.4775 k (18)
Thus, by considering both ( 17) and ( 18) we obtain (13).

In order to prove the necessity of the condition (13) let us consider the curves T u,u+1 . The points of the curve T u,u+1 that correspond to the frequency ω = 1.3483 rad/s lie on the line described by ( 16) so that the maximum delay uncertainty admissible for those points is exactly α/k. If we select a larger value for δ the system will become unstable at least on those points. This concludes the proof.

Remark 2. The fact that the maximum uncertainty allowed does not depend on the nominal delay τ 1 is a nice Remark 3. The condition (13) expresses a trade-o between the maximum delay mismatch δ and the proportional gain that can be used to tune the controller gain k.

Remark 4. This result improves the robust stability condition |∆| < 1/k found in [START_REF] Mascolo | Modeling the Internet congestion control as a time delay system: a robust stability analysis[END_REF].

Proposition 5. The system is stable, regardless the value of τ 1 , if the delay uncertainty ∆ satises the following inequality:

-

τ 1 < ∆ < -τ 1 + β k (19)
with β = 1.1188.

Proof.

The proof follows the same arguments of Proposition 1, therefore it is omitted.

Remark 6. The condition (19) implicitly requires the delay uncertainty ∆ to be negative, i.e. the nominal delay τ 1 should be always below the actual delay of the plant τ 2 .

Thus, condition (19) has no particular meaning for the characterization of controller robustness, since the sign of the uncertainty is not known a priori.

SIMULATION RESULTS

In this Section we report simulation results obtained by using a SIMULINK model that implements the system when τ 1 = 1 sec, k = 4 sec -1 , b(t) = 100 packets/sec, w = 150 packets evolution q(t) and the input rate r(t) when the delay uncertainty is either zero, δ/2 ∼ = 0.185 sec or δ ∼ = 0.37 sec.

The gure shows, as expected, that the performance of the closed loop response degrades when the delay uncertainty increases. In particular, oscillations are present when the delay uncertainty is δ/2 still providing an acceptable response, whereas when the delay uncertainty increases to the maximum allowed value persistent oscillations occur.

CONCLUSIONS

In this paper we have analyzed the robust stability of a very important class of congestion control algorithms when delay uncertainties are present. We have shown how the geometrical approach developed in [START_REF] Morarescu | On the geometry of stability regions of Smith predictors subject to delay uncertainty[END_REF] can be easily applied in order to nd stability bounds on the parameter of the controller. Moreover, we found a simple necessary and sucient condition on the gain of the proportional controller k in order to retain asymptotic stability regardless the value of the nominal delay τ 1 . Such a result suggests that congestion control algorithms that employ controllers made by a Smith predictor plus a proportional gain can be easily tuned in order to be robust to a bounded delay uncertainty.

- 1 ReFigure 2 .

 12 Figure 2. Triangle formed by the three vectors when the characteristic equation holds to show that the characteristic function in this case is given by:

Figure 3 .

 3 Figure 3. Stability crossing curves for the considered system transformation from (7) to (8) simply involves a scaling

Figure 4 .

 4 Figure 4. Stability regions for the considered system feature of the Smith predictor based controller. This makes the controller eective even with large delays.

depicted in Figure 1 .Figure 5 .

 15 Figure 5. Queue length (a) and input rate evolution (b)