Nadir Farhi
email: nadir.farhi@ifsttar.fr

Habib Haj-Salem
email: habib.haj-salem@ifsttar.frffijean-patrick.lebacque@ifsttar.fr

Jean-Patrick Lebacque

Multi-Anticipative Piecewise Linear Car-Following Model

Keywords:

We propose in this article an extension of the piecewise linear car-following model to multianticipative driving. As in the one-car-anticipative model, the stability and the stationary regimes are characterized thanks to a variational formulation of the car-dynamics. We study the homogeneous driving case. We show that in term of the stationary regime, the multianticipative model guarantees the same macroscopic behavior as for the one-car-anticipative one. Nevertheless, in the transient traffic, the variance in car-velocities and accelerations is mitigated by the multi-anticipative driving, and the car-trajectories are smoothed. A parameter identification of the model is made basing on NGSIM data and using a piecewise linear regression approach.

INTRODUCTION

We present a multi-anticipative car-following traffic model, where drivers control their velocities by taking into account the positions and the velocities of many cars ahead. In basic car-following models [START_REF] Chandler | Traffic dynamics: Studies in car following[END_REF][START_REF] Gazis | Nonlinear follow-the-leader models of traffic flow[END_REF][START_REF] Denos | Car-Following Theory of Steady-State Traffic Flow[END_REF][START_REF] Herman | Traffic dynamics: Analysis of stability in car following[END_REF], the car dynamics are described by stimulus-response equations that express the control process of drivers. Each driver accelerates or decelerates depending on his speed, and on the relative speed and the inter-vehicular distance with respect to the driver of the car ahead.

Multi-anticipative car-following models are often extensions of existing one car anticipative models. Since 1968, Bexelius [START_REF] Bexelius | An extended model for car-following[END_REF] extended the model of Chandler et al. [START_REF] Chandler | Traffic dynamics: Studies in car following[END_REF] to the multianticipative case. Lenz et al. [START_REF] Lenz | Multi-anticipative car-following model[END_REF] extended the model of Bando et al. [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF]. More recently, Hoogendoorn et al. [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF] have extended the model of Helly [START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF][START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF].

The model we present here is an extension of the piecewise linear car following model [START_REF] Farhi | Modelisation Minplus et Commande du Trafic de Villes Regulieres[END_REF][START_REF] Farhi | Piecewise linear car-following modeling[END_REF]. It is a first-order discrete-time model where the car-velocities are given in function of the intervehicular distances. We show here that the variational formulation made in [START_REF] Farhi | Modelisation Minplus et Commande du Trafic de Villes Regulieres[END_REF][START_REF] Farhi | Piecewise linear car-following modeling[END_REF] holds also for the multi-anticipative extension we propose in this article. That is to say that the car-dynamics are again interpreted as dynamic programming equations associated to stochastic optimal control problems of Markov chains, as in [START_REF] Farhi | Modelisation Minplus et Commande du Trafic de Villes Regulieres[END_REF][START_REF] Farhi | Piecewise linear car-following modeling[END_REF]. Thanks to that formulation, we are able to characterize the stability of the car-dynamics and to calculate the stationary regimes.

In the transient traffic, some qualitative results are obtained from our model. First, we simulated the car dynamics on a one lane-road of about 10,000 meters, where we imposed the trajectories of a given number of leaders, and simulated the trajectories of all the followers, basing on our model, and varying the number of leaders taken into account in anticipation. The effect of our anticipation modeling on the transient traffic is thus shown. As expected, we observed that, as the number of leaders taken into account in anticipation, increases, the car trajectories are smoothed, the distant followers are slowed down in unstable traffic phases, but they retrieve their non-anticipative trajectories once the traffic is stabilized.

Second, we proposed a parameter identification approach for our multi-anticipative model, and showed, in a basic example, how this method is applied. We rely on NGSIM data of vehicle trajectories on a segment of U.S. highway 101. The purpose of this part is rather to present the process of parameter identification, since the data used here are not exhaustive to draw conclusions, and the model is likely to be improved. However, some interesting observations are made on the data sample considered. For example, we observed that the scatter plot for space headway and average velocities seems to be easy to approximate by a piecewise linear curve in the case where the space headway is computed with respect to two leaders, compared to the case where the space headway is computed with respect to only one leader.

In the remainder of this introduction, we give a short review on the car-following modeling, we fix notations, and we introduce the next sections. We use the notations 𝑡 for time (discrete or continuous), 𝑥 for distance (car positions), and 𝑛 for the number of cars. The cars are numbered such that the first car (car number 1) is the leader. We consider the variables  𝑥(𝑛, 𝑡): the cumulative traveled distance of car 𝑛, from time zero to time 𝑡.  𝑦(𝑛, 𝑡): the inter-vehicular distance 𝑥 𝑛 -1, 𝑡 -𝑥 𝑛, 𝑡 .  𝑣(𝑛, 𝑡): the velocity of car 𝑛 at time 𝑡.

Car-following models are often based on a behavioral law 𝑉 𝑒 (equilibrium speed spacing function), that gives, at the equilibrium traffic, the velocity 𝑣 of a car 𝑛 as a function of the intervehicular distance 𝑦 between cars 𝑛 and 𝑛 -1. It is then assumed that the law 𝑉 𝑒 also holds on the transient traffic. A kind of general form of first-order car-following models can then be derived as done for the macroscopic first-order modeling (LWR models [START_REF] Lighthill | On kinematic waves ii: A theory of traffic flow on long, crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF]).

𝑣 𝑛, 𝑡 = 𝑉 𝑒

′ 𝑦 𝑛, 𝑡 ∆𝑣 𝑛, 𝑡 ,

where ∆𝑣 𝑛, 𝑡 = 𝑣 𝑛 -1, 𝑡 -𝑣 𝑛, 𝑡 .

The simplest form for the equilibrium speed-spacing function 𝑉 𝑒 𝑦 is the linear one 𝑉 𝑒 𝑦 = 𝛼𝑦 + 𝛽, where 𝛼 and 𝛽 are parameters. In this is case, the linear car-following model [START_REF] Chandler | Traffic dynamics: Studies in car following[END_REF][START_REF] Herman | Traffic dynamics: Analysis of stability in car following[END_REF] is obtained:

𝑣 𝑛, 𝑡 = 𝛼𝑦 𝑛, 𝑡 + 𝛽, (2)
The model ( 2) is, of course, not satisfactory. Another well known car-following model is due to Gazis, Herman, and Rothery [START_REF] Gazis | Nonlinear follow-the-leader models of traffic flow[END_REF].

𝑣 𝑛, 𝑡 + 𝑇 = 𝑎 𝑣(𝑛, 𝑡 + 𝑇) 𝑝 𝑣 𝑛 -

where a reaction time 𝑇 is considered, and where 𝑝 and 𝑙 are parameters. For example, if 𝑝 = 1 and 𝑙 = 2, the model can simply be obtained by taking 𝑉 𝑒 𝑦 = 𝑏 exp -𝑎/𝑦 in (1), satisfying

𝑉 𝑒 ′ 𝑦 = 𝑎 𝑉 𝑒 (𝑦)/𝑦 2 .
Other car-following models that do not necessarily match the form (1) exist, such as the models of Bando et al. [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF] (optimal velocity model), Helly [START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF][START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF], Treiber et al. [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF] (intelligent driver model), etc. Bando et al. [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF] proposed the optimal velocity model 𝑣 𝑛, 𝑡 + 𝑇 = 𝜆(𝑉 𝑒 𝑦 𝑛, 𝑡 -𝑣 𝑛, 𝑡 + 𝑇 )

This model has been much studied recently for being easily analyzed with mathematical tools.

Helly [START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF][START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF] considered the linear model 𝑣 𝑛, 𝑡 + 𝑇 = 𝛼 𝛥𝑣 𝑛, 𝑡 + 𝛽 𝑦 𝑛, 𝑡 -𝑆 𝑛 [START_REF] Farhi | Piecewise linear car-following modeling[END_REF] where 𝛼 and 𝛽 are parameters, and where 𝑆 𝑛 is the desired distance which can be linear 𝑆 𝑛 = 𝑆 0 + 𝑇𝑣 𝑛 , with 𝑆 0 the minimum gross distance between two cars. We base here on the piecewise linear car-following model proposed in [START_REF] Farhi | Piecewise linear car-following modeling[END_REF], where the behavioral law 𝑉 𝑒 is approximated with a (min-max)-piecewise linear curve

𝑉 𝑒 𝑦 = min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 𝛼 𝑢𝑤 𝑦 + 𝛽 𝑢𝑤 (6)
and where a one-car-anticipative discrete-time car-dynamics has been obtained

where 𝑀 𝑢𝑤 and 𝑐 𝑢𝑤 , for (𝑢, 𝑤) ∈ 𝑈 × 𝑊, are matrices and column-vectors respectively.

Two cases have been distinguished in [START_REF] Farhi | Piecewise linear car-following modeling[END_REF].

 The 𝜈 cars move on a ring road. In this case, 𝑀 𝑢𝑤 and 𝑐 𝑢𝑤 are given by

𝑀 𝑢𝑤 = 1 -𝛼 𝑢𝑤 0 ⋯ 𝛼 𝑢𝑤 𝛼 𝑢𝑤 1 -𝛼 𝑢𝑤 0 ⋮ ⋱ ⋱ ⋮ 0 0 𝛼 𝑢𝑤 1 -
where 𝑣 denotes the asymptotic car-velocity (the same for all cars), and 𝑦 denotes the average inter-vehicular distance in the ring road ( 𝑦 = 1/𝑑).  The 𝜈 cars move on an "open" road, where the velocity 𝑣 1 (𝑡) of the first car (the leader one) varies over time but is stationary. In this case, 𝑀 𝑢𝑤 and 𝑐 𝑢𝑤 are given by

𝑀 𝑢𝑤 = 1 0 ⋯ 0 𝛼 𝑢𝑤 1 -𝛼 𝑢𝑤 0 ⋮ ⋱ ⋱ ⋮ 0 0 𝛼 𝑢𝑤 1 -𝛼 𝑢𝑤
and 𝑐 𝑢𝑤 (𝑡) = 𝑣 1 (𝑡) , 𝛽 𝑢𝑤 , … , 𝛽 𝑢𝑤 𝑡 Again, the dynamics (8) is stable under the condition 𝛼 𝑢𝑤 ∊ [0,1] and the inverse behavior law at the stationary regime is obtained as follows.

𝑦 = max 𝑢 ∊𝑈 min 𝑤 ∊𝑊 𝑣 1 -𝛽 𝑢𝑤 𝛼 𝑢𝑤 (10)
where 𝑣 1 denotes the asymptotic velocity of the first car.

The dynamics [START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF] have been interpreted in [START_REF] Farhi | Piecewise linear car-following modeling[END_REF], under the assumption 𝛼 𝑢𝑤 ∈ 0,1 , ∀ (𝑢, 𝑤) ∈ 𝑈 × 𝑊, as a dynamic programming equation associated to a stochastic game on a controlled Markov chain; see [START_REF] Farhi | Piecewise linear car-following modeling[END_REF] for more details.

ANTICIPATION MODELLING

We present in this section an extension of the model [START_REF] Denos | Car-Following Theory of Steady-State Traffic Flow[END_REF] to multi-anticipative traffic, where each car chooses its velocity depending on the inter-vehicular distance with respect to a given number m of cars ahead the considered car (multi-leaders). In order to situate our model with respect to the existing multi-anticipative models, and to explain the extension we do, let us first give a short review on multi-anticipative car-following models.

A straightforward multi-leader extension of the model of Chandler et al. [START_REF] Chandler | Traffic dynamics: Studies in car following[END_REF] is the Bexelius model [START_REF] Bexelius | An extended model for car-following[END_REF] 𝑣 𝑛, 𝑡 + 𝑇 = 𝛼 𝑗 ∆𝑣 𝑗 𝑛, 𝑡 𝑚 𝑗 =1 [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF] where 𝑣 denotes the acceleration, 𝛼 𝑗 , 𝑗 = 1,2, … , 𝑚 are sensitivity parameters with respect to the 𝑗th car ahead, and where ∆𝑣 𝑗 𝑛, 𝑡 = 𝑣 𝑛 -𝑗, 𝑡 -𝑣(𝑛, 𝑡). This model is very simple but permits some mathematical analysis. Hoogendoorn et al. [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF] have noted the non convenience of the additive form of Bexelius model [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF], and proposed the modification

𝑣 𝑛, 𝑡 + 𝑇 = min 1≤𝑗 ≤𝑚 𝛼 𝑗 ∆𝑣 𝑗 𝑛, 𝑡 (12)
Hoogendoorn et al. [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF] have also proposed a multi-anticipative generalization for the Helly model ( 5)

𝑣 𝑛, 𝑡 + 𝑇 = 𝛼 𝑗 ∆𝑣 𝑗 𝑛, 𝑡 + 𝛽 𝑗 [∆𝑥 𝑗 𝑛, 𝑡 -𝑆 𝑗 (𝑛)] 𝑚 2 𝑗 =1 𝑚 1 𝑗 =1 (13)
where ∆𝑥 𝑗 𝑛, 𝑡 = 𝑥 𝑛 -𝑗, 𝑡 -𝑥(𝑛, 𝑡). Lenz et al. [START_REF] Lenz | Multi-anticipative car-following model[END_REF] have generalized the Bando model ( 4) as follows.

𝑣 𝑛, 𝑡 = 𝜅 𝑗 𝑉 𝑒 ∆𝑥 𝑗 𝑛, 𝑡 𝑗 -𝑣 𝑛, 𝑡 𝑚 𝑗 =1 (14)
where 𝜅 𝑗 expresses the sensitivity with respect to the 𝑗th leader.

We propose here a multi-leader extension for the piecewise linear car-following model [START_REF] Denos | Car-Following Theory of Steady-State Traffic Flow[END_REF]. We use a minimum form as in [START_REF] Lenz | Multi-anticipative car-following model[END_REF] (rather than an additive form as in [START_REF] Serge P Hoogendoorn | Empirics of multianticipative carfollowing behavior[END_REF]). Moreover, we use a uniform form for the sensitivity with respect to the inter-vehicular distance as in ( 14) (the intervehicular distance with respect to the 𝑗th leader is divided by 𝑗). We consider the dynamics:

𝑥 𝑛 𝑡 + 1 = 𝑥 𝑛 𝑡 + min 1≤𝑗 ≤𝑚 (1 + 𝜆) 𝑗 -1 min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 𝛼 𝑢𝑤 𝑥 𝑛-𝑗 𝑡 -𝑥 𝑛 𝑡 𝑗 + 𝛽 𝑢𝑤 ( 15
)
where 𝑚 is the number of leaders taken into account in anticipation, and 𝜆 ≥ 0 is a discount parameter with respect to the leader index. The dynamics (15) can be written simply

𝑥 𝑛 𝑡 + 1 = 𝑥 𝑛 𝑡 + min 1≤𝑗 ≤𝑚 min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 𝛼 𝑗𝑢𝑤 𝑥 𝑛-𝑗 𝑡 -𝑥 𝑛 𝑡 𝑗 + 𝛽 𝑗𝑢𝑤 ( 16
)
where ∀ 𝑢, 𝑤 ∈ 𝑈 × 𝑊, (𝛼 𝑗𝑢𝑤 ) 𝑗 , 1 ≤ 𝑗 ≤ 𝑚 are increasing non negative sequences, and (𝛽 𝑗𝑢𝑤 ) 𝑗 , 1 ≤ 𝑗 ≤ 𝑚 are increasing sequences. The interpretation of the minimum operator with respect to the 𝑗th leader in ( 16) is that a car 𝑛 maximizes its velocity under the constraints

𝑥 𝑛 𝑡 + 1 -𝑥 𝑛 𝑡 ≤ min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 𝛼 𝑗𝑢𝑤 𝑥 𝑛-𝑗 𝑡 -𝑥 𝑛 𝑡 𝑗 + 𝛽 𝑗𝑢𝑤 , 1 ≤ 𝑗 ≤ 𝑚.
One consequence of anticipation in driving is that the information that a car 𝑖, for 𝑖 = 𝑛 -1, 𝑛 -2, … , max 1, 𝑛 -𝑚 , decelerates at time 𝑡, is immediately transmitted to the car 𝑛 that reacts at time 𝑡 + 1, instead of 𝑡 + 𝑛 -𝑖. The discounting with respect to the leader indices, made by introducing the multiplicative term (1 + 𝜆) 𝑗 -1 , permits to favor closer leaders over distant ones. If 𝜆 = 0, then the cars respond equally to the stimulus of all the leaders 𝑗, with 𝑗 = 1,2, … , 𝑚.

In the following section (section 2), we study the stability of the car dynamics [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF], and characterize the existence of stationary regimes. Two cases are distinguished: traffic on a ring road, and traffic on an open road. In both cases, we give the asymptotic car positions when stationary regimes exist. The transient traffic for the car dynamics ( 15) is treated in section 3.

STABILITY ANALYSIS AND STATIONARY REGIMES

As in [START_REF] Farhi | Modelisation Minplus et Commande du Trafic de Villes Regulieres[END_REF][START_REF] Farhi | Piecewise linear car-following modeling[END_REF], we consider 𝜈 cars moving on a 1-lane road without passing. We first study the case where the cars move on a ring road, and then explore the "open" road case.

Traffic on a ring road

The cars being moving on a ring road, the indices 𝑛 -𝑗, in the dynamics [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF], are cyclic in the set {1,2, … , 𝜈}. The idea here is that the two minimum operators in ( 15) can be summarized in only one minimum operator, and then retrieve the one car-anticipative form for the dynamics. Let us denote by 𝑍 the set of all pairs of indices (𝑗, 𝑢), with 1 ≤ 𝑗 ≤ 𝑚 and 𝑢 ∈ 𝑈 𝑍 = 𝑧 = 𝑗, 𝑢 , 1 ≤ 𝑗 ≤ 𝑚, 𝑢 ∈ 𝑈 . The dynamics ( 15) is then written

𝑥 𝑛 𝑡 + 1 = min 𝑧∈𝑍 max 𝑤 ∊𝑊 { 𝑀 𝑧𝑤 𝑥 𝑡 ] 𝑛 + 𝑐 𝑛 𝑧𝑤 , 1 ≤ 𝑛 ≤ 𝜈 (17)
where the matrices 𝑀 𝑧𝑤 = 𝑀 𝑗𝑢𝑤 and the column vectors 𝑐 𝑧𝑤 = 𝑐 𝑗𝑢𝑤 are given as follows.

𝑀 𝑗𝑢𝑤 = The dynamics (18) have the same form as [START_REF] Helly | Simulation of bottlenecks in single lane traffic flow[END_REF]. It is then interpreted as a dynamic programming equation associated to a stochastic game on a controlled Markov chain. The stability is guaranteed under the condition 𝛼 𝑗𝑢𝑤 ∈ 0,1 , ∀ 𝑗, 𝑢, 𝑤 ∈ {1,2, … , 𝑚} × 𝑈 × 𝑊; see [START_REF] Farhi | Modelisation Minplus et Commande du Trafic de Villes Regulieres[END_REF][START_REF] Farhi | Piecewise linear car-following modeling[END_REF] for more details. The stationary regime is characterized by the additive eigenvalue problem

𝑣 + 𝑥 𝑛 = min 1≤𝑗 ≤𝑚 min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 { 𝑀 𝑗𝑢𝑤 𝑥] 𝑛 + 𝑐 𝑛 𝑗𝑢𝑤 , 1 ≤ 𝑛 ≤ 𝜈 ( 18
)
where 𝑣 is the asymptotic car-velocity, the same for all cars, the vector 𝑥 is the asymptotic carpositions, given up to an additive constant. The following result gives a solution for the system (19). ■ A particular case is important here, which is the non discounting case where 𝜆 = 0 in [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF]. In this case, the matrices 𝑀 𝑗𝑢𝑤 and the vectors 𝑐 𝑗𝑢𝑤 still depend on 𝑗, whilst the parameters 𝛼 𝑗𝑢𝑤 and 𝛽 𝑗𝑢𝑤 are independent of 𝑗 for all 𝑗 ∈ {1,2, … , 𝑚}. Thus the average car speed 𝑣 coincides with the average car-speed obtained in the one car anticipative model [START_REF] Denos | Car-Following Theory of Steady-State Traffic Flow[END_REF]:

𝑣 = min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 𝛼 𝑢𝑤 𝑦 + 𝛽 𝑢𝑤 ,
where 𝛼 𝑢𝑤 = 𝛼 𝑗𝑢𝑤 and 𝛽 𝑢𝑤 = 𝛽 𝑗𝑢𝑤 , ∀ 𝑗 ∈ {1,2, … , 𝑚}.

Traffic on an open road

We suppose here given the speed 𝑣 1 (𝑡) of the first car over time, since the cars move on an open road. In addition, in order to analyze the stability and the stationary regime of the car dynamics, we assume that the velocity of the first car approaches a constant value 𝑣 1 . That is to say that lim 𝑡→+∞ 𝑣 1 (𝑡) = 𝑣 1 . Moreover, the number of anticipation cars for the first 𝑚 cars cannot be 𝑚 (it is less than 𝑚). More precisely, the number of anticipation cars for a car numbered 𝑛 is min 𝑛 -1,𝑚.

The dynamics [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF] is written here

𝑥 𝑛 (𝑡 + 1) = min 1≤𝑗 ≤𝑚 min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 { 𝑀 𝑗𝑢𝑤 𝑥 𝑡 ] 𝑛 + 𝑐 𝑛 𝑗𝑢𝑤 , 1 ≤ 𝑛 ≤ 𝜈 (19)
where the matrices 𝑀 𝑗𝑢𝑤 and the vectors 𝑐 𝑗𝑢𝑤 are given by 𝑀 𝑗𝑢𝑤 :

1 0 0 ⋯ ⋯ 0 0 1 0 0 ⋯ 0 ⋮ ⋱ ⋮ 0 0 ⋯ 1 0 0 ⋯ 0 𝛼 𝑗𝑢𝑤 /𝑗 0 ⋯ 0 1 -𝛼 𝑗𝑢𝑤 /𝑗 ⋯ 0 0 𝛼 𝑗𝑢𝑤 /𝑗 0 ⋯ 0 1 -𝛼 𝑗𝑢𝑤 /𝑗 ⋮ ⋮ 0 ⋱ 0 ⋯ ⋱ 0 0 ⋯ 0 𝛼 𝑗𝑢𝑤 /𝑗 0 ⋯ 0 1 -𝛼 𝑗𝑢𝑤 /𝑗 1 2 ⋮ 𝑗 𝑗 + 1 𝑗 + 2 ⋮ 𝜈 𝑐 𝑗𝑢𝑤 = 𝑣 1 (𝑡) +∞ ⋮ +∞ 𝛽 𝑗𝑢𝑤 𝛽 𝑗𝑢𝑤 ⋮ 𝛽 𝑗𝑢𝑤 1 2 ⋮ 𝑗 𝑗 + 1 𝑗 + 2 ⋮ 𝜈
The entries (𝑀 𝑗𝑢𝑤 ) 𝑖𝑘 and (𝑐 𝑗𝑢𝑤 ) 𝑖 for 𝑖, 𝑘 ≤ 𝑗 do not play any role in the car dynamics since (𝑐 𝑗𝑢𝑤 ) 𝑖 = +∞, ∀𝑖 ≤ 𝑗. Those entries correspond to anticipation of a car 𝑖 with respect to its 𝑗 th leader that does not exist since 𝑖 ≤ 𝑗.

We assume that the velocity 𝑣 1 (𝑡) of the first car reaches a fixed value 𝑣 1 at the stationary regime. The stationary regime is thus characterized as follows.

𝑣 + 𝑥

𝑛 = min 1≤𝑗 ≤𝑚 min 𝑢 ∊𝑈 max 𝑤 ∊𝑊 { 𝑀 𝑗𝑢𝑤 𝑥] 𝑛 + 𝑐 𝑛 𝑗𝑢𝑤 , 1 ≤ 𝑛 ≤ 𝜈 (20)
The following result gives a solution for the system (21).

Theorem 2. For all 𝑦 ∈ 𝑅 satisfying min 𝑢∈𝑈 max 𝑤∈𝑊 𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 = 𝑣 1 , the couple (𝑣, 𝑥) is a solution for the system (21), where 𝑣 = 𝑣 1 and 𝑥 is given up to an additive constant by

𝑥 = 𝜈 -1 𝑦, 𝜈 -2 𝑦, … , 𝑦, 0 . 𝑡 (21)
Proof. The proof is similar to that of Theorem 3 of [START_REF] Farhi | Piecewise linear car-following modeling[END_REF]. ■ Theorem 2 gives the car-velocity at the stationary regime with one stationary configuration of cars (uniform configuration) and gives the optimal strategy for drivers at that regime. An important remark here is that the car-velocity obtained is the same as the one obtained for the one car anticipative model [START_REF] Denos | Car-Following Theory of Steady-State Traffic Flow[END_REF]. Moreover, the "optimal strategy" of driving at the stationary regime in the case of multi-anticipative model is to drive by taking into account only one leader: 𝑗, 𝑢, 𝑤 = (1, 𝑢, 𝑤). That is, at the stationary regime, once the traffic is stabilized, it is not necessary for drivers to take into consideration more than one leader.

Moreover, the stationary configurations of cars in the two cases of one car anticipation and multi-anticipation models may coincide. This case is interpreted as follows. Even though the cars reduce their approach in the multi-anticipative dynamics (due to the minimum operator over the leader indices), comparing to their movement under the one car anticipative dynamics; as long as the cars approach the stationary regime, where the traffic is stabilized, they retrieve what they have lost on the transient regime. Therefore, by introducing the minimum on the multianticipative dynamics, the traffic becomes smoother, without decreasing the stationary car speed.

TRANSIENT TRAFFIC

We take here the same example presented in [START_REF] Farhi | Piecewise linear car-following modeling[END_REF], which we adapt to multi-anticipative case, in order to make a comparison. We simulate the car-dynamics [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF]. We take as the time unit half a second (1/2 s), and as the distance unit 1 meter (m). The parameters of the model are the same as those of Example 1 of [START_REF] Farhi | Piecewise linear car-following modeling[END_REF] (the parameters have been determined by approximating a given behavior law). More precisely, The behavior law considered here is approximated by the following piecewise linear curve of six segments.

𝑉 𝑦 = max 𝛼 1 𝑦 + 𝛽 1 , min 𝛼 2 𝑦 + 𝛽 2 , 𝛼 3 𝑦 + 𝛽 3 , 𝛼 4 𝑦 + 𝛽 4 , 𝛼 5 𝑦 + 𝛽 5 , 𝛼 6 𝑦 + 𝛽 6 , where the parameters 𝛼 𝑖 and 𝛽 𝑖 for 𝑖 = 1,2, … , 6 are given by We simulate the car dynamics on a one lane road of about 10,000 meter. We vary the number of leaders and the velocity of the first car (over time) in order to show the effect of multi-anticipation on the transient traffic. We give the results on Table 1. We see easily on that table that the trajectories are smoothed by anticipation. We notice here that the trajectory of the first car is the same for all views of Table 2. Although the number of leaders taken into account by drivers cannot exceed 5 in practice, we simulated here the car dynamics with anticipation with up to 100 leaders. We did it for curiosity, but it can be interesting in the case, for example, where one likes to study the traffic of communicating cars or automatic ones, etc. We denote by 𝑦 (𝑚 ,𝜆) the variable 𝑦 𝑚 ,𝜆 (n, t) = min 1≤j≤m 1 + λ j-1 x n -j, t -x n, t j .

For fixed values of 𝑚 and 𝜆 and with measured velocities 𝑣(𝑛, 𝑡) and inter-vehicular distances 𝑥 𝑛 -𝑗, 𝑡 -𝑥(𝑛, 𝑡) for every car 𝑛, we approximate the scatter plot 𝑉(𝑦 ) by a piecewise-linear curve

𝑉 𝑦 = min 𝑢 ∈𝑈 max 𝑤 ∈𝑊 𝛼 𝑢𝑤 𝑦 + 𝛽 𝑢𝑤 . (22)
The min-max piecewise linear approximation ( 23) is based on a piecewise linear regression approach where the number of segments as well as the points where they intersect are determined optimally by deterministic dynamic programming. The optimization with respect to the parameters 𝑚 and 𝜆 is done numerically by varying the two parameters in convenient intervals, and then determine the optimal parameters.

For fixed values of parameters 𝑚 and 𝜆, the scatter plot 𝑉(𝑦 (𝑚 ,𝜆) ) is approximated by linear regression on separated intervals. The intervals are determined by a dynamic programming approach. More precisely, we divide the axis of inter-vehicular distances 𝑦 into unity intervals (𝑦 𝑖 , 𝑦 𝑖+1 ). We start from the first interval, make a linear regression on that interval, and then decide for the second interval whether we make only one linear regression for the two intervals, or approximate the scatter plot on the second interval with another linear segment. Then we do the same for the third unity interval, and so on. We solve here the optimal control problem associated to the decision process. Note that the inter-vehicular distance plays the role of time in this decision process.

The decision at the unity interval 𝑦 𝑖 , denoted here by 𝑟(𝑦 𝑖 ), is in 0,1 :

 𝑟 𝑦 𝑖 = 0, if we decide to make one linear regression for the interval 𝑦 𝑖 together with the intervals before it.

 𝑟 𝑦 𝑖 = 1 if we decide to start a new linear regression from the unity interval 𝑦 𝑖 .

We define costs 𝑘(𝑦, 𝑁, 𝑟) in order to minimize regression errors and penalize large segmentation (limit the number of segments used in the approximation).

 𝑘(𝑦 𝑖 , 𝑁, 𝑟) the error of regression at the unit interval (𝑦 𝑖 , 𝑦 𝑖+1 ), when the interval (𝑦 𝑖-𝑁 , 𝑦 𝑖 ) is wholy approximated by one segment, and when the decision 𝑟 is taken at the stage 𝑦 𝑖 .

Then the costs 𝑘(𝑦, 𝑁, 𝑟) are given by  𝑘(𝑦, 𝑁, 𝑟) = Linear regression error in (𝑦 𝑖-𝑁 , 𝑦 𝑖+1 ) if 𝑟 = 0,  𝑘(𝑦, 𝑁, 𝑟) = Linear regression error in 𝑦 𝑖 , 𝑦 𝑖+1 + 𝜑 if 𝑟 = 1, where 𝜑 is a penalization of starting a new linear regression. Then the following optimal control problem is solved min 𝛾 ∈𝛤 𝑘(𝑦, 𝑁, 𝑟)

𝑦 𝑚𝑎𝑥 𝑦 =0 , ( 23
)
where 𝛤 is the set of boolean strategies on {0,

The parameter identification method is thus summarized as follows. For each couple of parameters (𝑚, 𝜆), we calculate 𝑦 (𝑚 ,𝜆) . Then the curves 𝑉(𝑦 (𝑚,𝜆) ) are approximated by using the piecewise linear regression approach explained above (that is by solving the dynamic programming equation ( 26)). We obtain a total regression error for each approximation (𝑚, 𝜆). Finally, we determine the best couple (𝑚, 𝜆) that gives the minimal total error.

We show below, a first application of the identification method proposed above. We base here on NGSIM data of vehicle trajectories on a segment of U.S. Highway 101 (Hollywood Freeway) in Los Angeles, California. The data are collected between 7:50 a.m. and 8:05 a.m. on June 15, 2005. A preliminary analysis of the trajectories showed that, according to the multi-anticipative model presented here, the 𝑛th leaders' positions for 𝑛 > 3 are redundant in the data considered here, even with a null discount parameter (𝜆 = 0). That is to say that ∀ 𝜆 ≥ 0, ∀ 𝑛, 𝑡, 𝑦 𝑚 ,𝜆 𝑛, 𝑡 ≥ 𝑦 3,𝜆 𝑛, 𝑡 .

By consequent, we do not consider more than 3 leaders here (𝑚 ∈ {1,2,3}). Indeed, we can see from the dynamics [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF], that, for 𝜆 = 0, taking into account more than one leader in anticipation does not change anything in the case where the traffic is accelerating, because the spacing to the 𝑗th leader is bigger than 𝑗 times the spacing to the first leader. That is 𝑥 𝑛 -𝑗 𝑡 -𝑥 𝑛 𝑡 ≥ 𝑗. 𝑥 𝑛-1 𝑡 -𝑥 𝑛 𝑡 , 𝑗 = 1,2, … Nevertheless, this is realistic because, in case of accelerating traffic, anticipation is not significant. First, we show in Table 3 the scatter plot for 𝑦 (𝑚 ,0) with the average car-velocities (the average over all cars), for 𝑚 ∈ 1,2,3 . It seems that 𝑉(𝑦 (2,0) ) and 𝑉(𝑦 (3,0) ) (where we denote also by 𝑉 the average car-velocity over all cars) can easily be approximated with piecewise linear curves, comparing to 𝑉(𝑦 (1,0) ). We do not have justification for that presently.

TABLE 3 Average car-velocity function of the inter-vehicular distance. 𝝀 = 𝟎 1 leader 2 leaders 3 leaders

In Figure 1, we give the results of the parameter identification. The figure gives the total errors from the piecewise linear regression, obtained for varied values of parameters 𝑚 ∈ {1,2,3} and 𝜆 ∈ [0, 5], increased by 0.1. The optimal parameters obtained here are 𝑚 = 2 and 𝜆 = 1.5. As mentioned above, it is not worthy to consider values of 𝑚 that exceed 3 (for the data considered here). Also, for large values of 𝜆, we retrieve the same results as if we take 𝑚 = 1. Therefore, the optimization we did here is significant.

FIGURE 1 Identification of parameters 𝒎 and 𝝀

The scatter plot for 𝑦 (2,1.5) with 𝑉, is approximated by the following curve of four segments.

𝑉 𝑦 = max 0, min 0.38𝑦 -1.90, 0.11𝑦 + 2.95, 10 .

(

) 26
The approximation is shown in Figure 2. As shown in Figure 2, it is not easy to identify one behavior for a big number of drivers. In fact the parameter identification should be made for each driver. The identification of the behavior law for a randomly selected driver has given the following result. The optimal parameters (𝑚, 𝜆) are 𝑚 = 2 (anticipation with two leaders) and 𝜆 = 0 (no discounting). The curve 𝑉(𝑦 (2,0) ) is approximated as follows.

𝑉 𝑦 2,0 ≈ max 0, min 0.33𝑦 2,0 -1.71, 8 .

The approximation is shown in Table 4, where we have shown also the approximation of the curve 𝑉 𝑦 1,0 .

𝑉 𝑦 1,0 ≈ max 0, min 0.26𝑦 1,0 -0.9, 8 .

As obtained by the parameter identification approach, we can see on Table 4 that the curve 𝑉 𝑦 2,0 is well fitted by a piecewise-linear curve, comparing to the curve 𝑉 𝑦 1,0 . Note that the analytical results on the stability of the dynamics and the stationary regimes, presented in section 2 do not necessarily hold for the case of several driver behaviors. We will consider in a next step, the driver heterogeneity.

TABLE 4 Parameter identification of one driver behavior

Driver behavior 𝑉(𝑦 (1,0) ) -one leader anticipation Driver behavior 𝑉(𝑦 (2,0) ) -two leaders anticipation

CONCLUSION

We presented in this article an extension of the piecewise linear car following model to multianticipative driving. The minimum form used for taking into account more than one leader, with the discounting parameter used to favor closest leaders over distant ones, seem to be convenient to recapture the main characteristics of the anticipative traffic. We have, in particular, shown that the trajectories are smoothed by the anticipation on the transient traffic, without affecting the stationary regimes. That is, the anticipation do not slow down the stationary traffic. The minimum form used for anticipation allows us to make the same variational formulation as done in the onecar anticipative traffic, and by this, it allows us to characterize the stability of the car-dynamics and calculate the stationary regimes. The identification test we made here is only a first step of analyzing the proposed model. In-depth analyses on exhaustive data should be done in the future, in order to improve the modeling approach presented in this article. In particular, we shall try to extend the model in a way that it takes into consideration heterogeneity in driving.

TABLE 2

 2 Traffic on a 1-lane road. On the 𝒙-axis: time. On the 𝒚-axis: car-position. The number of cars taken into account in anticipation are 1, 5, 10, 20, 50 and 100. The length of the road is 10,000 meter. The total simulation time is 500 seconds

FIGURE 2

 2 FIGURE 2 Approximation of the law 𝑽(𝒚 (𝟐,𝟏.𝟓) ) with a min-max piecewise linear curve. The number of leaders taken into account is 𝒎 = 𝟐. The discounting parameter is 𝝀 = 𝟏. 𝟓

Theorem 1 .

 1 If ∀ 𝑗, 𝑢, 𝑤 ∈ 1,2, … , 𝑚 × 𝑈 × 𝑊, 𝛼 𝑗𝑢𝑤 ∈ [0,1], then the system (19) admits a solution (𝑣, 𝑥) given by: 𝑑 is the average inter-vehicular distance in the ring road. 𝛼 𝑗𝑢𝑤 ) 𝑗 and (𝛽 𝑗𝑢𝑤 ) 𝑗 are increasing sequences (with respect to 𝑗) and 𝑦 ≥ 0, we have

					𝑣 = min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 ,
					𝑥 =	𝜈 -1 𝑦, 𝜈 -2 𝑦, … , 𝑦, 0 .
	Proof. Following the same approach as in [4, 5], we obtain the following solution for the system
	(19).
	𝑣 = min 𝑧 ∊𝑍	max 𝑤 ∊𝑊	𝛼 𝑧𝑤 𝑦 + 𝛽 𝑧𝑤 = min 1≤𝑗 ≤𝑚	min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	{𝛼 𝑗𝑢𝑤 𝑦 + 𝛽 𝑗𝑢𝑤 } ,
	𝑥 =	𝜈 -1 𝑦, 𝜈 -2 𝑦, … , 𝑦, 0 .
	min 1≤𝑗 ≤𝑚	min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	{𝛼 𝑗𝑢𝑤 𝑦 + 𝛽 𝑗𝑢𝑤 } = min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 .

𝑡

where 𝑦 = 1/𝑡 Then since (

  Let 𝑦 ∈ 𝑅 satisfying

			min 𝑢∈𝑈	max 𝑤 ∈𝑊	𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 = 𝑣 1 .
	min 1≤𝑗 ≤𝑚	min 𝑢∈𝑈	max 𝑤 ∈𝑊 [𝑀 𝑗𝑢𝑤 𝑥] 𝑛 + 𝑐 𝑛 𝑗𝑢𝑤	= min 1≤𝑗 ≤𝑚	min 𝑢∈𝑈	max 𝑤 ∈𝑊	𝛼 𝑗𝑢𝑤 𝑦 + 𝛽 𝑗𝑢𝑤 + 𝑥 𝑛
						= min 𝑢 ∈𝑈 = 𝑣 1 + 𝑥 𝑛 . max 𝑤 ∈𝑊	𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 + 𝑥 𝑛
								1𝑢 𝑤 𝑦 +
	𝛽 1𝑢 𝑤 = 𝑣 1 . Indeed
		[𝑀 1𝑢 𝑤 𝑥] 𝑛 + 𝑐 𝑛 1𝑢 𝑤	= 𝛼 1𝑢 𝑤 𝑦 + 𝛽 1𝑢 𝑤 + 𝑥 𝑛
				= min 𝑢∈𝑈	max 𝑤∈𝑊	𝛼 1𝑢𝑤 𝑦 + 𝛽 1𝑢𝑤 + 𝑥 𝑛
				= min 1≤𝑗≤𝑚	min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	𝛼 𝑗𝑢𝑤 𝑦 + 𝛽 𝑗𝑢𝑤 + 𝑥 𝑛
				= min 1≤𝑗≤𝑚	min 𝑢 ∊𝑈	max 𝑤 ∊𝑊	[𝑀 𝑗𝑢𝑤 𝑥] 𝑛 + 𝑐 𝑛 𝑗𝑢𝑤 .

Let 𝑥 be given by (22). Then ∀𝑛 ∈ {1,2, … , 𝜈} we have Moreover, the optimal strategy at the stationary regime is 𝑗, 𝑢, 𝑤 = (1, 𝑢, 𝑤) such that 𝛼

TABLE 1 Approximation of the behavior law 𝑽 (𝒚) with a piecewise-linear curve

 1

	Segments	1	2	3	4	5	6
	𝛼 𝑖	0	0.54	0.32	0.13	0.34	0
	𝛽 𝑖	0	-8.1	-1.47	6.11	10.6	14

  𝑦 𝑖 , 𝑁) = min 𝑘 𝑦 𝑖 , 𝑁, 0 + 𝐺 𝑦 𝑖+1 , 𝑁 + 1 , 𝑘 𝑦 𝑖 , 𝑁, 1 + 𝐺 𝑦 𝑖+1 , 1 .

	𝐺(𝑦 𝑚𝑎𝑥 ) = 0
	𝐺(
	1, … , 𝑦 𝑚𝑎𝑥 }. That is 𝛾 ∶ 0,1, … , 𝑦 𝑚𝑎𝑥 ∋ 𝑦 ↦ 𝑟 ∈
	{0,1}. The value function associated to (24) is
	𝑦 𝑚𝑎𝑥
	𝐺 𝑦, 𝑁 =	𝑘 𝑧, 𝑁, 𝑟	(24)
	𝑧=𝑦
	𝐺 satisfies the dynamic programming equation

PARAMETER IDENTIFICATIONWe propose here a parameter identification method based on (piecewise) linear regression. Given measured data of the car positions and velocities on a given section, the method permits to determine the optimal parameters that match our model with the measured data. Since the dynamics of the model is simply the car-velocities given as functions of the inter-vehicular distances, then the optimal parameters would be the ones that approximate the scatter plot of instantaneous inter-vehicular distances and velocities.


Y B Én Édic

P Wira

J Merckl
email: jean.merckle@uha.fr

A NEW METHOD FOR THE RE-IMPLEMENTATION OF THRESHOLD LOGIC FUNCTIONS WITH CELLULAR NEURAL NETWORKS

A new strategy is presented for the implementation of threshold logic functions with binary-output Cellular Neural Networks (CNNs). The objective is to optimize the CNNs weights to develop a robust implementation. Hence, the concept of generative set is introduced as a convenient representation of any linearly separable Boolean function. Our analysis of threshold logic functions leads to a complete algorithm that automatically provides an optimized generative set. New weights are deduced and a more robust CNN template assuming the same function can thus be implemented. The strategy is illustrated by a detailed example.

Introduction

Cellular Neural Networks (CNNs) represent an important paradigm of Artificial Neural Networks (ANNs). Indeed, all cells of a CNN are identical and operate simultaneously in addition to every other friendly properties of ANNs. Through their particular architecture, CNNs are well suited for hardware implementation and their parallel computing properties can thus be fully exploited. As a consequence, the theoretical computational speed of CNNs can be very fast compared to other types of ANNs, specially in real-time applications.

Nevertheless, their correct operation is sensitive to noise. This happens mainly because of the electrical components that compose the cell and which are not ideal ones. Indeed, electrical components are affected by noise, imperfections in the fabrication process, or post-manufacturing disturbance like temperature. A CNN implementation is an analog processor, informations are therefore represented by electrical signals and operations are achieved by electrical circuits. Some quantification errors can therefore appear and be easily propagated resulting in erroneous behavior of cells for some tasks. This is referenced to as the CNN internal noise. Binary-output CNNs, which share this property of being realizable in hardware, achieve final activations of ±1 with certain restrictions over the original CNN model. Furthermore, with inputs and initial internal states restricted to Boolean values, their overall behavior is then the one of a Threshold Logic Function (TLF).

Hence, this article addresses the problem of implementing arbitrary threshold logic functions with noise-robust CNN templates. Therefore, we introduce a new method based on an original mathematical framework for linearly separable Boolean functions referred to as generative set. Synthesizing a CNN template with generative sets allows to optimize its weights and to increase its robustness to the internal noise. As an advantage, the pro-posed method uses any correct template as an a priori knowledge to decrease the time needed for its optimization.

The paper is structured as follows. In section 2, the CNN concept is presented and different ways to compute weights are reviewed. The paradigm of the generative set is then introduced in section 3 as an alternative representation of TLFs. ¿From an analysis of a CNN template, we deduce an algorithm able to compute its generative set. Section 4 proposes an algorithm that optimizes any known generative set. This section also shows how to use it in order to resynthesize the original template into a more robust one. Each step of this procedure is illustrated with an example. Finally, concluding remarks are offered in the last section.

Cellular Neural Networks

The CNN concept

CNNs where first introduced by Chua and Yang in [START_REF] Chua | Cellular neural networks: Theory[END_REF] . They are large arrays of identical nonlinear dynamic systems called cells which satisfy two properties: interactions are local within a finite radius, and all state variables are continuous-valued signals. CNNs thus consist of cells connected only to the cells in their neighborhood. They are arranged in an internal layer and are connected to both the input layer and the output layer through neighboring cells, as shown by Fig. 1. Fig. 1. Schematics of a two-dimensional CNN with a 3 × 3-wide neighborhood.

The main feature of CNNs, which are considered as a special class of recurrent neural networks, is therefore that information is directly exchanged just between neighboring cells. As a result, the connectivity of CNNs is greatly reduced and makes them particularly suitable for on-chip implementations, so allowing high speed parallel processing for real-time applications in many fields 2 : image processing, artificial vision, partial differential equation solution, nonlinear phenomena modeling, biological systems modeling, robot control, etc.

A mathematical formal description of the internal dynamics of a cell is given by the following equation with t the time variable:

dx(t) dt = -x(t)+A * N(y(t))+B * N(u(t))+z, (1)
where the scalars x(t), u(t) and y(t) represent respectively the internal state, the input and the output of the cell. The matrix A holds the weights from the feedback connections (output) and likewise the matrix B holds the weights from the feedforward connections (input). z is a bias term. Since all the neurons of a CNN are identical, the values A, B, and z are the same for every cells and are therefore called a template of the CNN. N(.) is a neighborhood function and the symbol " * " stands for a linear convolution product. The output y(t) is computed from x(t) thanks to the following piecewise-linear activation function:

y(t) = 1 2 (|x(t) + 1| -|x(t) -1|). (2)
Binary-output CNNs are derived from the latter model by giving the self-feedback weight a0 in A a value greater than 1 1 . Restricting the values of the inputs and the initial internal states to Boolean values, then leads to the conclusion that a binary-output CNN behaves like a threshold logic gate [START_REF] Monnin | Boolean design of binary initialized and coupled CNN image processing operators[END_REF] :

y(t → ∞) = 1 ⇔ A * N(y(t)) + B * N(u(t)) -x(t) > -z ∀ t ≥ 0 (3)
The implementation of a TLF is achieved by setting A, B and z.

Weights computation

Using conventional learning algorithms to find a template results in real-valued weights not necessarily optimal. The problem then remains that of finding appropriate weights. Appropriate weights are weights that yield robust CNNs, i.e., more tolerant against internal noise. Once determined, these weights are kept constant during the application of the CNN template to a specific problem.

Since CNNs are able to implement any arbitrary linearly separable Boolean function [START_REF] Chen | Realization of Boolean functions via CNN: Mathematical theory, LSBF and template design[END_REF] , they can thus be analyzed with similar mathematical theories. This has led to many approaches, which have been developed over these last years.

Of these approaches, the most widely used is the look-up tables of Chow's parameters [START_REF] Chow | On the characterization of threshold functions[END_REF] . Indeed, a set of parameters is shown to fully characterize TLFs of up to eight variables. Thus, it is possible to recover the TLF implementation from a look-up table which would store the implementation of every possible TLF, along with their respective Chow's parameters. In a practical sense, such tables exist for TLF with less than six variables 6 , but methods have been developed to break higher dimensional TLF into smaller ones [START_REF] Picton | Deriving weights for single threshold logic gates using decomposition[END_REF] . Other techniques have been proposed, either inspired from Karnaugh map minimization procedures [START_REF] Celinski | A mapping technique for the synthesis of linear threshold gate networks used to implement Boolean functions[END_REF] , or based on a formal analysis [START_REF] Chen | Realization of Boolean functions via CNN: Mathematical theory, LSBF and template design[END_REF] , on an analytical description [START_REF] Hänggi | Cellular neural networks: Analysis, design and optimization[END_REF] or on an algebraic representation [START_REF] Monnin | Boolean design of binary initialized and coupled CNN image processing operators[END_REF] .

Despite the strength of these methods, none of them is able to optimize the robustness of a given template without starting from scratch, thus leading to time consuming procedures. As opposed to them, the method developed thereafter is able to alter an existing CNN design in order to enhance its robustness to noise. It is based on specific properties of TLFs from which it derives a new arithmetic expression of a TLF: its generative set [START_REF] Bénédic | Fast analysis method for both coupled and uncoupled binary-output cloning templates[END_REF] . A simple algorithm is then used to tweak the template into a more robust one, according to the properties revealed by its generative set 11 .

Generative Sets

Introduction

A generative set is an alternative representation of a threshold logic function, making good use of the latter's algebraic properties such as symmetry and monotony. It can be directly computed from a cloning template thanks to an analatical process detailed in 10 , and which is briefly recalled in this section.

A threshold logic function F is formally defined by the following thresholded weighted sum:

F (s) ≡ (w * s > T h) ≡ n i=1 wisi > T h , (4)
where T h is the so-called threshold and the weights wi and inputs si are the elements of w and s respectively.

One can notice the similar look of Eq. ( 3) and Eq. ( 4). In the latter, the n elements of A and B were simply gathered into a single weighting vector w, whereas the inputs and outputs from N(y(t)) and N(u(t)) were collected into s. The expression of T h is known since 1 :

T h = (1 -a0)x(0) -z, (5)
where a0 > 1 is the self-feedback weight value of A.

As a result, a binary-output cloning template ends up being fully characterized by the threshold logic function's attributes: w, s and T h. Without loss of generality, we will assume from now on that the weights wi are positive, and that they are sorted such that wi ≥ wj if i > j. Before dealing with the generative set extraction algorithm, the following subsection introduces the definitions of symmetry and monotony in the Boolean algebra.

Definitions

The previous general expression of a TLF implies two important properties: symmetry and monotony. The smart usage of these properties makes the generative set representation to be more efficient than the truth table by being less redundant.

Definition 1 (The symmetry of a TLF)

A TLF is called symmetric if and only if it is invariant under any permutation of its inputs. A TLF is only partially symmetric in a set of inputs if and only if the permutation of these inputs leaves the function unchanged [START_REF] Kohavi | Switching and finite automata theory[END_REF] .

Definition 2 (The monotony of a TLF)

A TLF is 12 :

• constant in one input si if and only if its output does not depend on si;

• increasing/decreasing in si if and only if its output increases/decreases with the increase of si.

Thanks to the commutativity property of the addition in Eq. ( 4), the output of any arbitrary threshold logic function F is invariant under any permutation of equally weighted inputs and therfore symmetric in these inputs. As a result the value of F can only be affected by the number of true inputs among a set of symmetrical ones rather than by the way these true values are mapped. Now let c be the number of sets of equally weighted variables and ni the number of variables gathered in the set labeled i. Notice that a set might contain just one variable and that c i=1 ni = n. Any input vector s can then be replaced by a vector k, called reduced input vector and built as follows.

Definition 3 (Reduced input vector) A reduced input vector is a vector k = (kc, . . . , k1) made of reduced variables ki which are the number of true inputs within the i th set of symmetrical inputs of s. Opposed to the values of the inputs si which are restricted to ±1, the values of ki are defined by ki = 0, 1, . . . ni.

As a consequence to the Def. 3, two input vectors s and s ′ , derived from each others by one or more permutations of equally weighted inputs, can both be replaced by the same reduced input vector k. This leads to a loss of redundancy in the description of a TLF made by the reduced input vectors compared to the one made with the usual input vectors. As a result, F (s) can be re-written as a function of k: F (s) ≡ F (k) and the Eq. ( 4) becomes:

F (k) ≡ c i=1 wi(2ki -ni) > T h , (6)
where the equal weight values have been factorized and the respective set of symmetrical inputs have been replaced by their reduced input counterpart [START_REF] Bénédic | Fast analysis method for both coupled and uncoupled binary-output cloning templates[END_REF] .

The symmetries are thus used to define a reduced input space in which the description of any TLF is much less redundant than in the usual input space. The following uses the monotony property to bring another improvement to the description of TLF.

Let k + be a reduced vector such that F (k + ) ≡ true. Such a vector is called a minterm of F . Since all the weights are considered positive, increasing any element k + i of k + clearly increases the value of the weighted sum of the Eq. ( 6). As a result, F is increasing in each one of its c reduced variables. Therefore, increasing either element of a given reduced vector k will only affect F (k) in one out of the following three ways:

• the weighted sum is not risen enough and F (k) remains f alse;

• the weighted sum reaches the threshold and F (k) switches from f alse to true;

• the weighted sum remains over the threshold and F (k) remains true.

This shows that once a minterm k + of F is known, many more can be derived from it, simply by increasing the value of at least one of its elements. These derived minterms are thus useless to the representation of a TLF. Hence, there must be a quite small set of minterms that can generate every single minterm of F . This set is called the generative set of the function.

Definition 4 (The generative set of a TLF)

The set of minterms that generate every single minterm of F is referred to as the generative set of F and is written {F}. It gathers generative vectors, i.e., any minterm k of F that cannot be derived from another minterm. As a convention, the generative vectors are indexed so as to be sorted in an anti-lexical ordering.

According to the previous definition, the generative set {F} of a TLF obviously describes F completely. Furthermore, the number of generative vectors included in {F} can be considered as a way to measure the complexity of the function F . The following subsection describes a way to extract the generative set of the TLF performed by a CNN template.

Analyzing algorithm

The analyzing algorithm discussed in this subsection aims at extracting the generative set of a TLF implemented by a template w and T h. It is assumed that k = (kc, . . . , k1) and w = (wc, . . . , w1), with wi > wj if and only if i > j.

For a given value κc of kc, Eq. ( 6) writes:

c-1 i=1 wi(2ki -ni) > T h -wc(2κc -nc), (7)
which is the definition of a TLF of c -1 reduced variables. As a result, the generative set of a TLF of c variables is the concatenation of admissible values of kc with the generative sets obtained from their respective TLF of c -1 variables. This remark leads to the idea that the analyzing algorithm should be recursive. Such a recursion would stop when a TLF of one variable and a threshold value T h ′ are reached:

w1(2k1 -n1) > T h ′ . (8)
In this case, the generative set consists in a single mono-dimensional reduced vector, which component is the lowest integer value of k1 that satisfies the latter inequality. It formally writes:

k1 = 1 + 1 2 • w1 T h ′ + w1 • n1 . (9)
An adequate threshold value has to be forwarded from one recursive level to another in order to be evaluated and eventually stop the algorithm. This new threshold is defined at each recursive level by:

T hc = T h + c i=1 wini. (10)
Thanks to this, Eq. ( 7) becomes:

F (k) ≡ c-1 i=1 2 • wi • ki > T hc + 2wcκc . ( 11
)
Finally, the threshold value which must be forwarded to a recursive launch for a given value κc of kc is:

T hc,κ c = T hc -2 • κc • wc. (12)
One can see that this algorithm requires exactly Π c i=1 (1 + ni) recursive launches, which is usually lower than 2 n , assuming that there exists i such that ni = 1. This number can be significantly reduced. It is the case when the value of a reduced variable is high enough to ensure that the underlying TLF is always true, which means that the resulting generative set is the null vector and does not need to be computed. Speaking of kc, this limit value is:

κc,max = 1 + T hc 2 • wc . (13)
Likewise, a reduced input value might be too small to allow the underlying TLF to be true, thus leading to an empty generative set. The smallest value of kc for which this does not happen is given by:

κc,min = 1 + 1 2 • wc T hc -2 c-1 i=1 wi • ni .
(14) The way the expressions of both Eq. ( 13) and ( 14) are obtained is detailed in [START_REF] Bénédic | Fast analysis method for both coupled and uncoupled binary-output cloning templates[END_REF] .

These two limits define a set of allowed integer values of kc: [min(nc, κc,max -1); max(0, κc,min)] for which recursive launches must be made. The analyzing algorithm finally writes:

Illustrative example

This example considers an arbitrary CNN template to illustrate, through real values, how the analyzing method works. The weights can result from any of the synthesis methods previously discussed in section 2.2. For example: The inputs values are ±1 and the initial state is x(0) = 0. Notice that the weights of B are sorted and are all positive. In this example, no weights are equal, thus ni = 1 for every i and c = 4.

A =   0 0 0 0 a0 0 0 0 0   =   0 0 0 0 3 0 0 0 0   , B =   w1 w2 w3 w4 0 0 0 0 0   =   1
The algorithm is initialized with T h given by ( 5) and with T h4 given by (10). In the present case T h = -z = -2.5 and T h4 = 11.5.

For the sake of clarity, the operations involved by the analyzing algorithm are detailed line by line:

l. 1 since c = 4, the algorithm jumps to line 10. l. 10 κ4,max = 2 according to (13).

l. 11 κ4,min = 0 from (14).

l. 13 since κ4,max > n4, the vector (2, 0, 0, 0) is not added to the generative set.

l. 17 the iteration will start at k4 = 1 down to k4 = 0.

l. 19 recursive launch for k4 = 1:

l. 1 since c = 3, the algorithm jumps to line 10. l. 10 κ3,max = 1 according to (13). l. 11 κ3,min = 0 from (14). l. 13 since κ3,max ≤ n3, the vector (1, 0, 0) is added to the generative set. l. 17 the iteration will start at k3 = 0 down to k3 = 0. l. 19 recursive launch for k3 = 0: l. 1 since c = 2, the algorithm jumps to line 10. l. 10 κ2,max = 1 according to (13). l. 11 κ2,min = 0 from (14). l. 13 since κ2,max ≤ n2, the vector (1, 0) is added to the generative set. l. 17 the iteration will start at k2 = 0 down to k2 = 0. l. 19 recursive launch for k2 = 0: l. 1 since c = 1, the algorithm jumps to line 5. l. 5 κ1,max = 1 according to (13). l. 6,7 since κ1,max ≤ n1, the vector (1) is added to the generative set. l. 25 the generative set is returned. l. 22 the vector (0, 1) is added to the generative set. l. 23 no more iterations. l. 25 the generative set is returned. l. 22 the vectors (0, 1, 0) and (0, 0, 1) are added to the generative set. l. 23 no more iterations. l. 25 the generative set is returned. l. 22 the vectors (1, 1, 0, 0), (1, 0, 1, 0) and (1, 0, 0, 1) are added to the generative set.

l. 23 recursive launch for k4 = 0:

l. 1 since c = 3, the algorithm jumps to line 10. l. 10 κ3,max = 2 according to (13). l. 11 κ3,min = 1 from (14). l. 13 since κ3,max > n3, the vector (1, 0, 0) is not added to the generative set. l. 17 the iteration will start at k3 = 1 down to k3 = 1. l. 19 recursive launch for k3 = 1:

l. 1 since c = 2, the algorithm jumps to line 10. l. 10 κ2,max = 1 according to (13). l. 11 κ2,min = 0 from (14). l. 13 since κ2,max ≤ n2, the vector (1, 0) is added to the generative set. l. 17 the iteration will start at k2 = 0 down to k2 = 0. l. 19 recursive launch for k2 = 0: l. 1 since c = 1, the algorithm jumps to line 5. l. 5 κ1,max = 1 according to (13). l. 6,7 since κ1,max ≤ n1, the vector (1) is added to the generative set. l. 25 the generative set is returned. l. 22 the vector (0, 1) is added to the generative set. l. 23 no more iterations. l. 25 the generative set is returned. l. 22 the vectors (1, 1, 0) and (1, 0, 1) are added to the generative set. l. 23 no more iterations. l. 25 the generative set is returned.

l. 22 the vectors (0, 1, 1, 0) and (0, 1, 0, 1) are added to the generative set.

l. 23 no more iterations. l. 25 the generative set is returned.

Finally, putting together the five generative vectors which were returned by the first recursive level of the algorithm gives the following generative set:

{F} =            (1, 1, 0, 0) (1, 0, 1, 0) (1, 0, 0, 1) (0, 1, 1, 0) (0, 1, 0, 1)            . (15)
The generative set (15) completely describes the operator studied in this example. The low redundancy of this representation is visible. Indeed, five vectors are only necessary to describe the operator compared to the sixteen rows usually involved in a truth table of four inputs. It is also noticable that the algorithm took care of the anti-lexical ordering mentioned in Def. 4.

Like this example, any binary-output cloning template can be analyzed into a generative set with the proposed algorithm.

Optimization of CNN templates

The behavior of a binary-output cloning template, i.e. of a TLF, is described as the linear separation of the 2 n vertices of a n-dimensional hypercube into the two classes labeled -1 and 1. As previously stated by Eq. ( 4), the separation is achieved by an hyperplane H defined by w and T h. The task implemented by a given CNN is considered non-optimal, or non-robust to internal noise, if the hyperplane is poorly positioned in the hypercube, i.e., being unnecessarily close to one of its vertices. This property is crucial since CNNs generally deal with analog signals, meaning that even binary values (namely, the inputs and the initial internal states) will slightly vary about their mean value. We now propose an algorithm which partially solves this optimization problem. It is based on an analysis of the generative set of the TLF achieved by a CNN. By looking for non-implemented symmetries, this algorithm tweaks the original hyperplane to make it fit them. It results in an improved implementation of the TLF.

Generative sets reduction

The previously described algorithm is an efficient tool to know which TLF a given binary-output CNN performs. As seen in section 3, the returned Boolean expression is given in a reduced input space of c dimensions, as opposed to the input space which is the usual n-dimensional hypercube that was discussed so far. This reduction relies on the values n = (nc, . . . , n1) introduced in the previous section, and which are representative of the way symmetries were implemented in the original hyperplane H.

If it turned out from the analysis of the generative set implemented by H that a symmetry is not implemented, say between the reduced inputs ki and ki+1, then it would mean that n can be even more reduced into n ′ = (nc, . . . , ni + ni+1, . . . , n1). This extra reduction can then be applied to the generative vectors: ∀k ∈ {F}, k → k ′ = (kc, . . . , k i + k i+1 , . . . , k1).

(16) This extra reduction affects the weights in the following manner:

w → w ′ = wc, . . . , niwi + ni+1wi+1 ni + ni+1 , . . . , w1 , (17)
whereas the threshold value is maintained:

T h ′ = T h.
The resulting hyperplane H ′ = {w ′ , Th ′ } was proven to implement the same TLF as H in [START_REF] Bénédic | Optimisation of binary-output CNNs: First step of an analytical design process[END_REF] . The following subsection explains why H ′ is more robust to internal noise than H.

Toward an optimal hyperplane

Given an implementation of a given TLF, the aim of any analytical optimization process is to find a better (if not the best) location T h and orientation w of the separating hyperplane H. An adequate representation of the TLF which implementation is supposed to be optimized must therefore be chosen. The most straightforward one probably is the truth-table. The flaw behind this approach is that it is not threshold-logic-specific: every piece of informations about the equation of the sought hyperplane remain spread out, and yet, difficult to collect. On the contrary, the generative set representation is very close to the hyperplane's coordinates (see Fig. 2).

true false

(0, 0) (2, 4) {F} = (0, 3) (2, 2)
Fig. 2. A TLF represented in its two-dimensional reduced space. The vertices from its generative set are pointed by the arrows. Amongst the many hyperplanes which implement a given TLF, the optimal ones are defined as those which distances to their respective closest vertices of the hypercube, are the furthest. More formally:

Definition 5 (Optimal hyperplanes) Consider a hyperplane H defined by {w, T h}. Let s [START_REF] Chua | Cellular neural networks: Theory[END_REF] be the regular input vector assigned to the closest vertice of the n-dimensional hypercube, s [START_REF] Chua | Cellular neural networks: Applications[END_REF] the second closest one, and so on. H is said optimal, for the TLF it implements, if and only if 11 :

1. the distance D (1) = T h -n i=1 s (1)
i wi , between H and s [START_REF] Chua | Cellular neural networks: Theory[END_REF] , is greater than or equal to the distance between any other hyperplane implementing the same TLF, and its respective closest vertice;

amongst the hyperplanes satisfying item 1,

H is the hyperplane which distance D [START_REF] Chua | Cellular neural networks: Applications[END_REF] to its second closest vertice s [START_REF] Chua | Cellular neural networks: Applications[END_REF] is again greater than or equal to the other corresponding distances;

3. and so on with every remaining vertice.

The following theorem establishes the fact that an optimal hyperplane necessarily fits the symmetries of the TLF it implements. This means that it gives the same value to the respective weights.

Theorem 1 Let H (o) be an optimal hyperplane defined by {w (o) , T h (o) } for a given TLF. If this TLF is partially symmetrical in two distinct inputs si and sj (with i = j), then w

(o) i = w (o) j .
The proof of this theorem is given in [START_REF] Bénédic | Optimisation of binary-output CNNs: First step of an analytical design process[END_REF] . It relies on the fact that if an optimal hyperplane, which does not satisfy Th. 1, is found, then it can systematically be turned into a more robust hyperplane under the meaning of Def. 5 thanks to the weight alteration of the Eq. ( 17).

The consequence of Th. 1 thus is that the closest an hyperplane H is from being optimal, the closest its descriptor n is to the real symmetries of the implemented TLF, and reciprocally. The optimization process hence consists in the two following stages:

1. identify any non-implemented symmetry between two or more reduced inputs;

2. reduce them using Eq. ( 17).

This procedure is called: generative set reduction.

Optimization algorithm

The first stage of the generative set reduction is to check the symmetry of every pair of reduced inputs ki, kj, i = j, using the generative set of the TLF which implementation is to be optimized. Testing every combination leads to c(c -1)/2 tests. The real number of tests is actually much lower than this. Indeed, having the weights sorted implies that if ki and ki+j were symmetrical, then so would ki, ki+1, ki+2, . . . and ki+j be. As a result, one only needs to check kc with kc-1, kc-1 with kc-2 and so on until k2 is finally checked with k1.

In this manner, the number of symmetry tests is reduced to c -1 which is linear in the number of reduced inputs. The reduction of a generative set is achieved through the Th. 2, using the following notations and terms:

{F} =         
(k1,c, . . . , k1,2, k1,1) (k2,c, . . . , k2,2, k2,1)

. . .

(kq,c, . . . , kq,2, kq,1)

         =          k1 k2 . . . kq          , ( 18
) where q is number of generative vectors included in the set, and c is the dimension of the reduced input space.

In Eq. (18), we call the j th i-order root of {F} the head (kj,c . . . kj,c-i+1) obtained from the generative vectors kj with 1 ≤ i ≤ c and 1 ≤ j ≤ q.

To each i-order root value is associated a iorder sub-generative set extracted from {F}. It is made of the generative vectors which have a same iorder root value, and with their root part removed.

For example, the 2-order roots and subgenerative sets of the generative set (15) are highlighted by Fig. 3. processing time depends on the complexity of the threshold logic function, i.e., on its number of generative vectors as well as on the number of symmetries.

We also applied the proposed technique to the templates listed in the CNN bank presented in [START_REF] Chen | Realization and bifurcation of Boolean functions via cellular neural networks[END_REF] . A wide range of templates were improved. Indeed, over the 20,097 templates which were published, 17,799 were successfully optimized in an overall time of 2 minutes and 22 seconds.

Conclusion

This paper introduces a method for the reimplementation of threshold logic functions applied to CNNs. Based on a set of mathematical theories, the paradigm of generative sets is developed as a new representation of linearly separable Boolean functions. Thanks to it, the symmetries of any TLF can be detected more efficiently than with a truth table approach and their CNN implementation can thus be optimized against internal noise through symmetry reduction. To this effect, this paper provides the two core algorithms involved in the overall optimization process: the first one extracts the generative set from a correctly operating CNN template while the second detects symmetries in the generative set structure and reduces them to lead to a more robust CNN template.

Computer simulations have shown that most of the existing CNN templates can be efficiently rewritten with the proposed method. As a result, any template designed for a specific task can be optimized with our approach. The result is a template more robustness to internal noise.

Fig. 3 . 2 -

 32 Fig. 3. 2-order roots and sub-generative sets of the TLF of Eq. (15).

Theorem 2 (Generative set reduction)

Let {F} denote a generative set with the notations of the Eq. ( 18). The two reduced inputs ki and ki-1 are symmetrical if and only if for all j, the (c -i)-order sub-generative sets {Fj} each satisfy: test 1 every possible integer value of the sum (kj,i + kj,i-1), from the smallest which appears in {Fj} to the greatest one, indeed appear in {Fj};

test 2 for each one of these sum values, every possible combination of values kj,i and kj,i-1, giving the right sum appear in {Fj};

test 3 the (c -i + 2)-order sub-generative sets of {Fj} attached to a same value of (kj,i + kj,i-1) are identical.

A formal proof of this theorem is provided in [START_REF] Bénédic | Optimisation of binary-output CNNs: First step of an analytical design process[END_REF] . It is the translation into the reduced input space of Def. 1.

Once the symmetrical reduced inputs have been detected, the optimized weights can be computed using Eq. ( 17). The reduction algorithm thus is: A New Method for the Re-Implementation of Threshold Logic Functions with Cellular Neural Networks

Illustrative example

In order to illustrate the optimization process, we propose to synthesize a more robust cloning template from the template used in the previous example (see section 3.4). Its analysis resulted in the generative set given in (15) and recalled hereafter:

(copy of Eq. ( 15))

with c = 4 reduced variables, c -1 = 3 symmetry tests are thus needed. Applying the optimization algorithm to this generative set is detailed thereafter.

Symmetry k4 with k3

Sub-generative sets extraction

There is only one 0-order sub-generative set to deal with: the generative set itself. Only sub-generative set Merging k4 with k3 according to Eq. ( 16) yields the potential reduced input k4,3 = k4 +k3. Its value respectively is 2, 1, 1, 1 and 1 for each one of the five generative vectors of Eq. ( 15): test 1 There is no missing integer between the two values max(k4,3) = 2 and min(k4,3) = 1. test 2 The only way to obtain 2 is k4 = 1/k3 = 1 which achieved by the first generative vector of Eq. (15); whereas the two ways to obtain 1 are k4 = 1/k3 = 0 and k4 = 0/k3 = 1 which are achieved by at least one generative vector each. test 3 There is only one 2-order subgenerative set associated with k4 = 1/k3 = 1 and it is therefore equal to itself; the 2-order sub-generative set associated with k4 = 1/k3 = 0 and k4 = 0/k3 = 1 are the same. Reduction of k4 and k3 The two reduced inputs are symmetrical and their weights are reduced to the value w3,4 = 1 2 (w4 + w3) = 5.25. The generative set obtained after the reduction and with the duplicate generative vectors removed is:

Symmetry k3 with k2

Sub-generative sets extraction

There is only one 0-order sub-generative set to deal with: the generative set itself. Only sub-generative set Merging k3 with k2 according to Eq. ( 16) yields the potential reduced input k3,2 = k3 + k2. Its value respectively is 2, 2 and 1 for each one of the three generative vectors of the generative set (19). test 1 There is no missing integer between the two values max(k3,2) = 2 and min(k3,2) = 1. test 2 The two ways to obtain 2 are k3 = 2/k2 = 0 and k3 = 1/k2 = 1 and they are achieved by the first two generative vectors; whereas the two ways to obtain 1 are k3 = 1/k2 = 0 and k3 = 0/k2 = 1 but the second possibility is missing. The symmetry test thus fails.

test 3 Not necessary. Reduction of k3 and k2 The two reduced inputs are not symmetrical and their weights cannot be reduced.

Symmetry k2 with k1

Sub-generative sets extraction The two 1-order sub-generative sets are:

and

First sub-generative set Merging k2 with k1 according to Eq. ( 16) yields the potential reduced input k2,1 = k2 + k1. Its only value is 0. test 1 There is no missing integer between the two values max(k2,1) = 0 and min(k2,1) = 0. test 2 The only way to obtain 0 is k2 = 0/k1 = 0 and it is achieved by the only generative vector. test 3 There is no 2-order sub-generative sets, they are thus equal. Second sub-generative set Merging k2 with k1 according to Eq. ( 16) yields the potential reduced input k2,1 = k2 + k1. Its only value is 1. test 1 There is no missing integer between the two values max(k2,1) = 1 and min(k2,1) = 1. test 2 The only two ways to obtain 1 are k2 = 1/k1 = 0 and k2 = 0/k1 = 1; they are both achieved by at least one generative vector. test 3 There is no 2-order sub-generative sets, they are thus equal. Reduction of k2 and k1 The two reduced inputs are symmetrical in both {F1} ′ and {F2} ′ and therefore, so are they in {F } ′ . Their weights are thus reduced to the value w2,1 = 1 2 (w2 + w1) = 1.75 and the obtained generative set is:

Finally, the cloning template from the subsection 3.4 finally gives the following optimized template:

A and z remaining unchanged.

Discussion

All these operations have been implemented in a C++ application which can be used to evaluate the computational costs of the processing chain analysis/reduction described in this paper. For example, templates with four inputs (like the one used in the previous example) approximately take 13 ms to be optimized whereas templates with 16 inputs need 693 ms. These durations are mean values over a great number of templates and were executed on a Pentium Centrino Duo 1.8 GHz. The exact


Yuri Goegebeur
email: yuri.goegebeur@imada.sdu.dk

Armelle Guillou
email: armelle.guillou@math.unistra.fr

Andréhette Verster
email: verstera@ufs.ac.za

Robust and asymptotically unbiased estimation of extreme quantiles for heavy tailed distributions

A robust and asymptotically unbiased extreme quantile estimator is derived from a second order Pareto-type model and its asymptotic properties are studied under suitable regularity conditions. The finite sample properties of the proposed estimator are investigated with a small simulation experiment.

Introduction

In extreme value statistics, the estimation of extreme quantiles of a distribution function is a central topic. Indeed, many important applications in climatology, finance, actuarial science, hydrology and geology, to name but a few, require extrapolations outside the data range, and extreme value theory provides the only realistic framework for such an exercise. In the present paper we shall address this estimation problem, with special focus on asymptotic unbiasedness and robustness against outliers. F (x) = Cx -1/γ (1 + γ -1 δ(x)), is ultimately nonzero, of constant sign and |δ| ∈ RV τ .

Clearly, condition (R) implies that the tail quantile function U , defined as U (y) := inf{x : F (x) ≥ 1 -1/y}, y > 1, satisfies y -γ U (y) → C γ as y → ∞ and the function a implicitly defined by

U (y) = C γ y γ (1 + a(y)) (1)
satisfies a(y) = δ(C γ y γ )(1 + o(1)) as y → ∞, so |a| ∈ RV ρ , with ρ = γτ .

The second order condition (R) can be used to derive the so-called extended Pareto distribution, EPD [START_REF] Beirlant | Statistics of Extremes -Theory and Applications[END_REF][START_REF] Beirlant | Second-order refined peaks-over-threshold modelling for heavy-tailed distributions[END_REF], with distribution function given by

G(y) =    1 -[y(1 + δ -δy τ )] -1/γ , y > 1, 0, y ≤ 1, (2)
where γ > 0, τ < 0, and δ > max{-1, 1/τ }. As shown in Proposition 2. given that X > u can be approximated by ( 2) with δ = δ(u) up to an error that is uniformly

o(δ(u)) for u → ∞. In Dierckx et al. ( 2013
), a robust and asymptotically unbiased estimator for γ was introduced by fitting the EPD to a sample of relative excesses by the minimum density power divergence (MDPD) criterion [START_REF] Basu | Robust and efficient estimation by minimizing a density power divergence[END_REF]. In particular, let X 1 , . . . , X n be independent and identically distributed (i.i.d.) random variables with a distribution function satisfying (R), and denote by X 1,n ≤ • • • ≤ X n,n the corresponding order statistics. The parameters γ and δ of the EPD are then estimated with the minimum density power divergence criterion applied to the relative excesses over the random threshold u = X n-k,n , namely Y j := X n-k+j,n /X n-k,n , j = 1, . . . , k, i.e. one minimises the empirical divergence

∆ α (γ, δ) := ∞ 1 g 1+α (y)dy -1 + 1 α 1 k k j=1 g α (Y j ),
in case α > 0, and (2009) studied the asymptotically unbiased estimation of small tail probabilities based on the EPD, fitted by the maximum likelihood method. In Gomes and Pestana (2007) an asymptotically unbiased extreme quantile estimator was introduced for heavy-tailed distributions. These approaches are however not robust against outliers. To the best of our knowledge, robust and asymptotically unbiased extreme quantile estimation has not been considered before.

∆ 0 (γ, δ) := - 1 k k j=1 log g(Y j ), in case α = 0,
The remainder of our paper is organised as follows. In the next section we will introduce the robust and asymptotically unbiased estimator for extreme quantiles and study its asymptotic properties under suitable regularity conditions. The finite sample behaviour of the proposed estimator and some alternatives from the literature is illustrated with a small simulation experiment in Section 3.

From the second order condition (R) and using the EPD as approximation to the distribution of X/u n given X > u n we can for F (u n ) → 0 and p n → 0 such that

p n /F (u n ) → c ∈ [0, ∞) introduce U 0 1 p n := u n p n F (u n ) -γ 1 -δ(u n ) 1 - p n F (u n ) -ρ (3)
as approximation for U (1/p n ).

Lemma 1 Assume (R). If F (u n ) → 0 and p n → 0 such that p n /F (u n ) → c ∈ [0, ∞) we have that U 0 (1/p n )/U (1/p n ) → 1.
The proof of this lemma is straightforward and therefore it is for brevity omitted from the paper. Now, let X 1 , . . . , X n be i.i.d. random variables with a distribution function satisfying (R), and

denote by X 1,n ≤ • • • ≤ X n,n the corresponding order statistics. Taking u n = X n-k,n , replacing
F by the empirical distribution function in (3), and using the fact that e -x ∼ 1 -x for x → 0, we can introduce the following extreme quantile estimator

U 1 p n := X n-k,n np n k -γn exp -δ n 1 - np n k -ρn , (4)
where ( γ n , δ n ) is the MDPD estimator for (γ, δ) and ρ n is a consistent estimator sequence for ρ.

In order to study the asymptotic behaviour of U (1/p n ), properly normalised, we need some preliminary results. Firstly, we need the limiting distribution of the MDPD estimator for ( γ n , δ n ).

This was already derived in

δ n := δ(X n-k,n ). Theorem 1 Let X 1 , . . . , X n be a sample of i.i.d. random variables from a distribution function satisfying (R). Then if k, n → ∞ with k/n → 0 and √ ka(n/k) → λ ∈ R, we have that √ k   γn -γ 0 δn -δ n   (Γ, ∆) with (Γ, ∆) ∼ N 2 0, C -1 (ρ 0 )B(ρ 0 )Σ(ρ 0 )B (ρ 0 )C -1 (ρ 0 ) ,
where Σ(ρ 0 ) is a symmetric (3 × 3) matrix with elements

σ 11 (ρ 0 ) := α 2 (1 + γ 0 ) 2 [1 + α(1 + γ 0 )] 2 [1 + 2α(1 + γ 0 )] , σ 21 (ρ 0 ) := α(1 + γ 0 )[α(1 + γ 0 ) -ρ 0 ] [1 + α(1 + γ 0 )][1 -ρ 0 + α(1 + γ 0 )][1 -ρ 0 + 2α(1 + γ 0 )] , σ 22 (ρ 0 ) := [α(1 + γ 0 ) -ρ 0 ] 2 [1 -ρ 0 + α(1 + γ 0 )] 2 [1 -2ρ 0 + 2α(1 + γ 0 )] , σ 31 (ρ 0 ) := γ 0 1 [1 + 2α(1 + γ 0 )] 2 - 1 [1 + α(1 + γ 0 )] 3 , σ 32 (ρ 0 ) := γ 0 1 [1 -ρ 0 + 2α(1 + γ 0 )] 2 - 1 [1 + α(1 + γ 0 )] 2 [1 -ρ 0 + α(1 + γ 0 )] , σ 33 (ρ 0 ) := γ 2 0 2 [1 + 2α(1 + γ 0 )] 3 - 1 [1 + α(1 + γ 0 )] 4 , C(ρ 0 ) is a symmetric (2 × 2) matrix with elements c 11 (ρ 0 ) := γ -α-2 0 1 + α 2 (1 + γ 0 ) 2 [1 + α(1 + γ 0 )] 3 , c 12 (ρ 0 ) := γ -α-2 0 ρ 0 (1 -ρ 0 )[1 + α(1 + γ 0 ) + α 2 (1 + γ 0 ) 2 ] + α 3 ρ 0 (1 + γ 0 ) 3 [1 + α(1 + γ 0 )] 2 [1 -ρ 0 + α(1 + γ 0 )] 2 , c 22 (ρ 0 ) := γ -α-2 0 (1 -ρ 0 )ρ 2 0 + αρ 2 0 (1 + γ 0 )[α(1 + γ 0 ) -ρ 0 ] [1 + α(1 + γ 0 )][1 -ρ 0 + α(1 + γ 0 )][1 -2ρ 0 + α(1 + γ 0 )] ,
and

B(ρ 0 ) := γ -α-2 0   γ 0 0 -1 γ 0 -γ 0 (1 -ρ 0 ) 0   .
Secondly, we need the limiting distribution of the intermediate order statistic X n-k,n under (R), properly normalised.

Lemma 2 Let X 1 , . . . , X n be a sample of i.i.d. random variables from a distribution function

satisfying (R). For k, n → ∞ such that k = o(n) and √ ka(n/k) → λ ∈ R we have that √ k X n-k,n U (n/k) -1 X
where X ∼ N (0, γ 2 0 ).

In the next theorem we state the limiting distribution of the extreme quantile estimator (4), when properly normalised.

Theorem 2 Let X 1 , . . . , X n be a sample of i.i.d. random variables from a distribution function

satisfying (R). Then if k → ∞ as n → ∞ with k/n → 0, √ ka(n/k) → λ ∈ R, np n /k → 0 and ln(np n )/ √ k → 0 we have that √ k ln k npn   U 1 pn U 1 pn -1   Γ.
Theorem 2 indicates that the normalised extreme quantile estimator inherits the asymptotic distribution of the MDPD estimator for γ 0 . As shown in [START_REF] Dierckx | An asymptotically unbiased minimum density power divergence estimator for the Pareto-tail index[END_REF], the MDPD estimator for γ 0 based on the EPD is robust against outliers and asymptotically unbiased.

Simulation experiment

In this section we investigate the finite sample properties of U (1/p n ) as given in (4) with different parameter estimators, in particular the MDPD estimator γ n and δ n with α = 0.1, 0.5 and 1, and the maximum likelihood estimator (corresponding to MDPD with α = 0, see also [START_REF] Beirlant | Second-order refined peaks-over-threshold modelling for heavy-tailed distributions[END_REF]. We also consider the Weissman estimator [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF] given by

U W (1/p n ) = X n-k,n np n k -H k,n
, with H k,n being Hill's estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]. For the parameter ρ we use the estimator of Fraga [START_REF] Alves | A new class of semi-parametric estimators of the second order parameter[END_REF]. • Uncontaminated Fréchet distribution (Figure 1): F (x) = exp(-x -β ), x > 0, β > 0, denoted Fréchet(β). For this study β was chosen as 2.

• Contaminated Fréchet distribution: F (x) = (1 -)F (x) + F (x) where F (x) represents the uncontaminated Fréchet(2) and F (x) = 1 -(x/x c ) -β , x > x c where β is chosen as 0.5 and x c = 2 times the 99.99% quantile of the uncontaminated Fréchet(2). We take = 0.01 (Figure 2) and = 0.02 (Figure 3).

• Uncontaminated Burr distribution (Figure 4):

F (x) = 1 -(η/(η + x τ )) λ , x > 0, η, τ, λ > 0,
denoted Burr(η, τ, λ). For this study we have chosen η = 1, τ = 1 and λ = 2.

• Contaminated Burr distribution: F (x) = (1 -)F (x) + F (x) where F (x) represents the uncontaminated Burr(1,1,2) and F (x) = 1 -(x/x c ) -β , x > x c where β = 0.5 and x c = 1.2 times the 99.99% quantile of the uncontaminated Burr(1,1,2). We take = 0.01 (Figure 5) and = 0.02 (Figure 6).

We report only the results for quantile 1 -1/500. The 1 -1/1000 quantile was also considered and resulted in similar outcomes.

In with distribution function H(y) = 1 -1/y, y > 1. Thus, For L 1 , use the well-known fact that √ k(k/nY n-k,n -1) Z where Z ∼ N (0, 1) (see for instance Corollary 2.2.2 in de Haan and Ferreira, 2006) and the delta method to obtain that L 1 X under the conditions of the lemma. For the term L 2 , a straightforward application of Taylor's theorem gives Since a is regularly varying we have that a(tx)/a(t) → x ρ as t → ∞, locally uniformly for x > 0.

√ k ln X n-k,n U (n/k) D = √ k ln U (Y n-k,n ) U (n/k) = γ 0 √ k ln k n Y n-k,n + √ k ln 1 + a(Y n-k,n ) 1 + a(n/k) =: L 1 + L 2 .
L 2 = √ ka(n/k) a(Y n-k,n ) a(n/k) -1 + o(1) + o a(Y n-k,n ) a(n/k) .
Combining this with the fact that k/nY n-k,n → 1 a.s. and the assumption √ ka(n/k) → λ ∈ R we have that L 2 P → 0. Lemma 2 follows then by collecting the terms and another application of the delta method.

Proof of Theorem 2

First we comment on the joint convergence in distribution of the random vector

( √ k( γ n -γ 0 ), √ k( δ n -δ n ), √ k(X n-k,n /U (n/k) -1), ρ n ).
( √ k( γ n -γ 0 ), √ k δ n ) (Γ, ∆), where (Γ, ∆) ∼ N 2 ((0, λ), C -1 (ρ 0 )B(ρ 0 )Σ(ρ 0 )B (ρ 0 )C -1 (ρ 0 )).

  3 of Beirlant et al. (2009), for distribution functions satisfying (R), the distribution function of the relative excess Y := X/u

  where g denotes the density function of G. The parameter ρ is estimated externally, e.g. by using one of the estimators proposed in Fraga Alves et al. (2003) or Goegebeur et al. (2010). Other robust estimators for γ were proposed by e.g. Peng and Welsh (2001), Juárez and Schucany (2004), Vandewalle et al. (2007), Kim and Lee (2008).In the present paper we will consider robust and asymptotically unbiased extreme quantile estimation under model (R), using the MDPD estimator of[START_REF] Dierckx | An asymptotically unbiased minimum density power divergence estimator for the Pareto-tail index[END_REF].Beirlant et al.

  Dierckx et al. (2013), but we repeat the result for completeness here. Let the arrow denote convergence in distribution, and let P → denote convergence in probability. From now on we denote by γ 0 and ρ 0 the true values of the parameters γ and ρ, respectively, and

Figures 1 to 6

 6 Figures 1 to 6 illustrate the results of a small simulation study based on 100 datasets, each of size n = 200, simulated from the distributions given below. The same distributions were considered in Dierckx et al. (2013).

Figures 1 to 6 , 5 - 6 )Figure 1 :Figure 2 :

 65612 Figure 1: Fréchet simulation, quantile 1-1/500. Top: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (dotted), α = 0.5 (solid) and α = 1 (dashed). Bottom: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (solid), MLE based estimator (dotted) and Weissman estimator (dashed).

Figure 3 :

 3 Figure 3: Fréchet simulation, quantile 1-1/500, 2% contamination. Top: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (dotted), α = 0.5 (solid) and α = 1 (dashed). Bottom: median (left) and MSE (right) of the MDPD based estimator with α = 0.5 (solid), MLE based estimator (dotted) and Weissman estimator (dashed).

Figure 4 :

 4 Figure 4: Burr simulation, quantile 1-1/500. Top: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (dotted), α = 0.5 (solid) and α = 1 (dashed). Bottom: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (solid), MLE based estimator (dotted) and Weissman estimator (dashed).

Figure 5 :

 5 Figure 5: Burr simulation, quantile 1-1/500, 1% contamination. Top: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (dotted), α = 0.5 (solid) and α = 1 (dashed). Bottom: median (left) and MSE (right) of the MDPD based estimator with α = 0.5 (solid), MLE based estimator (dotted) and Weissman estimator (dashed).

Figure 6 :

 6 Figure 6: Burr simulation, quantile 1-1/500, 2% contamination. Top: median (left) and MSE (right) of the MDPD based estimator with α = 0.1 (dotted), α = 0.5 (solid) and α = 1 (dashed). Bottom: median (left) and MSE (right) of the MDPD based estimator with α = 0.5 (solid), MLE based estimator (dotted) and Weissman estimator (dashed).

C - 1 → λ and ρ n P → ρ 0 1 + a 1 pn-δ n 1 -d ρn n ( 5 )=: T 1 + T 2 + T 3 -T 4 .

 10151234 From the proof of Lemma 1 and Theorem 2 in[START_REF] Dierckx | An asymptotically unbiased minimum density power divergence estimator for the Pareto-tail index[END_REF] we can deduce that γ n and δ n are independent of X n-k,n , and therefore( √ k( γ n -γ 0 ), √ k δ n , √ k(X n-k,n /U (n/k) -1)) (Γ, ∆, X),where (Γ, ∆, X) ∼ N 3 ((0, λ, 0), Ψ)(ρ 0 )B(ρ 0 )Σ(ρ 0 )B (ρ 0 )C -1 (ρ 0 ) we have also that( √ k( γ n -γ 0 ), √ k( δ n -δ n ), √ k(X n-k,n /U (n/k) -1), ρ n ) (Γ, ∆, X, ρ 0 ). Now, consider ln( U (1/p n )/U (1/p n )). Let d n := k/(np n ). Straightforward calculations give ln n -γ 0 ) ln d n + ln 1 + a n k Clearly T 1 = O P (1/ √ k) byLemma 2 and T 2 = O P By using the regular variation properties of the function a and the fact that d n → ∞ we have that a(1/p n )/a(n/k) → 0, and thus under the conditions of the theorem T 3 = O(a(n/k)). Finally, for T 4 note that d ρn n = o P (1) and δ n = O P (1/ √ k). Collecting all the terms we see thus that the rate of convergence of ln( U (1/p n )/U (1/p n )) is given by ln dn √ k . Multiplying both sides of (5) by √ k/ ln d n , the result of the theorem follows.

Appendix Proof of Lemma 2

Using the inverse probability integral transform we have that

denotes order statistic n -k of a random sample Y 1 , . . . , Y n from the unit Pareto distribution