David Da Silva

M Brancaccio

Bruno Duplessis

J Adnot

Room air conditioner load control under summer comfort constraint

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Increasing interest in demand side management strategies appear more and more as a solution to electrical network problems. These problems include the need for future network investments, an increasing peak power consumption, transmission and distribution congestion problems which can be solved or reduced with the implementation of DSM strategies and smartgrid equipments.

Strategies such as load control can have a large impact in the resolution of these problems, especially in system stress situations. However there is a lack of methodologies to quantify the benefits and costs of DSM [START_REF] Strbac | Demand side management: Benefits and challenges[END_REF]. DSM strategies can be applied to several electrical appliances, like HVAC systems, hot water heaters, washing machines and others.

In France, there is an important summer electrical gradient of 500 MW/°C [START_REF] Rte | Bilan prévisionnel de l'équilibre offre-demande d'électricité en France[END_REF]. Regarding the load shift potential, and due to its coincidence with the summer electrical peak (Figure 1), around 13h (1p.m), heat pumps control load can have big impact in stressed electric networks. 30000 35000 40000 45000 50000 55000 60000 0 0 : 3 0 0 2 : 0 0 0 3 : 3 0 0 5 : 0 0 0 6 : 3 0 0 8 : 0 0 0 9 : 3 0 1 1 : 0 0 1 2 : 3 0 1 4 : 0 0 1 5 : 3 0 1 7 : 0 0 1 8 : 3 0 2 0 : 0 0 2 1 : 3 0 2 3 : 0 0 (MW) Regarding flexibility, systems that can store electric energy in another type of energy (for example thermal energy) are among the most interesting for DSM [START_REF] Stadler | Power grid balancing of energy systems with high renewable energy penetration by demand response[END_REF] because they allow longer curtailment times and can be combined with heat pumps. In this article we focused on heat pumps ("reversible split systems"), without any type of heat pump storage component (like cold water tanks), considering only the building thermal inertia as a flexibility component.

Reddy (Reddy, T.A et al., 1991) developed a methodology to evaluate load curtailment potential of an air conditioner as a function of building thermal inertia. Nevertheless his analysis was based in a simple equivalent thermal parameter (ETP) model (one resistance and one capacity) that does not reflects the real temperature fluctuation in a building [START_REF] Orphelin | Méthodes pour la reconstituition de courbes de charge agrégées des usages thermiques de l'électricité[END_REF] Srinivas [START_REF] Srinivas | Evaluation of Residential HVAC Control Strategies for Demand Response programs[END_REF]) evaluated different control strategies applied to HVAC systems. However the impact of different external conditions, during a control, regarding peak power and the subsequent comfort impact are not quantified.

Advanced control strategies where studied by [START_REF] Kyoung-Ho | Model-based demand-limiting control of building thermal mass[END_REF], where the heat pump control takes into account the building thermal mass and occupants to create an optimized pre-cooling. Although this strategy seems to be very attractive, advanced control is needed to correctly apply this strategy. Since this article is focused on small tertiary and residential buildings this kind of strategies will not be evaluated.

In this paper a methodology is developed to quantify the value of different load control strategies for room air conditioners (RAC), in small tertiary and residential buildings, taking into account occupant 's comfort, defined in EN 15251 (EN 15251, 2007), and the electrical grid needs for load shift.

SIMULATION METHODS AND EQUIPMENT SPECIFICATION

Simulation characteristics

To be able to couple Multizone building thermal evolutions with the RAC operation, a set of simulations was done with TRNSYS software [START_REF] Trnsys | TRNSYS 16 with Isibat[END_REF], where the energy effects (positives and negative) and comfort impacts of heat pump control load shift could be determined. These simulations make possible a random characterization of the control impacts due to weather evolutions. A small office building with two floors and a total surface area of 1000 m², representing 25% of the French stock, was chosen [START_REF] Andre | Analysis of building heating and cooling demands in the purpose of assesing the reversibility and heat recovery potentials[END_REF]. The building plan is the same for the two floors and it is shown hereafter. The building is simulated for the "Nice -France" climate using Meteonorm data, since this region suffers congestion and peak consumption problems during summer.

In order to simulate the heat pump (room air conditioner) operation, full-load and part-load, a model created by Marchio and Filfil (Marchio & Filfli, 2004) was implemented in TRNSYS and coupled to the building. Standard manufacturer parameters were used in this study.

For each building room the heat pump power was dimensioned for the maximum sensible and latent needs. Since almost all split heat pumps available in the market do not control humidity in the room, we considered a surplus of 20% of capacity to take into account the latent loads.

Simulations were made with a 1 minute time step to be able to correctly evaluate the comfort impact of each load control strategy.

Heat pump control and comfort criterion

To correctly consider the occupants comfort, the recommended interior operative temperatures for energy calculation, category II, from the EN 15251 were considered (EN 15251, 2007). There is a difference between the heat pump control temperature (air temperature) and the operative temperature, which is used to define comfort conditions within EN 15251 (EN 15251, 2007). Caillet [START_REF] Caillet | Simulation of control options for HVAC managements of typical office building[END_REF] has shown that a constant air temperature set point is not able to guaranty comfortable operative temperatures, taking into consideration the EN 15251 temperature ranges.

To simulate occupants' behavior, it was supposed that the occupants can interact with each RAC air set point at the end of each hour to achieve an operative temperature of 24.5 °C, during occupation hours, achieving this way a comfortable operative temperature. As it can be seen in Figure 3, to achieve an operative temperature around 24.5 °C, at 1 p.m, the heat pump set point needs to be set to 22.3 °C, e.g. to maintain a constant operative temperature the set point must vary with the building cooling loads. However during a control strategy, it is assumed that the occupants can stand an operative temperature of 27°C, passing from the equivalent 10 % PPD to 15 % PPD1 , or assuming that a ventilator is available in the room, increasing the temperature tolerance at the maximum of one degree (according to EN 15251).

The thermal comfort is evaluated in the room with the larger cooling needs ("Office 2" south orientation -Figure 2). This way the maximum discomfort is quantified and it is ensured that other rooms will have a lower comfort impact.

Control strategies for load curtailment

For each strategy applied there can be more or fewer gains in terms of power consumption reduction, however their comfort impacts can be very different. To evaluate the potential of load shift and load reduction strategies during the French summer electrical network peak, three load control strategies were selected: i. Total Load curtailment (TLC) -one hour operation suspension ii. Multiple load curtailments (MLC) -two times 30 minutes operation interdiction with a time interval of 30 minutes between them (combined interdiction: 30 + 30 = 1 hour) iii. Capacity reduction (CR) -Reduction of the heat pump capacity by 50 % for two hours Each one of the selected strategies is evaluated, in terms of energy impacts, for the entire building. Our study is focused in the French electric network, where the electrical network peak occurs normally around 1 p.m every summer day working. So the described strategies were applied each day of the cooling season at 1 p.m and compared with the reference case (normal operation without control strategy). This way the energy and comfort value of each strategy can be evaluated for the entire cooling season.

During heat waves the electric network becomes less reliable due to congestion problems and peak consumption (electric thermal gradient). Therefore emergency load controls can play an important role, avoiding severe problems like black-outs. To evaluate the role of the selected strategies during a heat wave, simulations were made for during ten consecutive days, where the daily average temperature was 24°C.

Analysis indicators definition

The heat pump operation is a function of the building thermal state. It is expected that, if a load curtailment control is applied to a heat pump, the building thermal will change. Imagining load curtailment on a hot day, the building temperature will raise, meaning that when the heat pump is turned on again, it will have to compensate the energy that was absorbed by the building (Figure 4 overconsumption). In order to analyze the energy effects and the discomfort caused by a heat pump load control, we define, hereafter, the following indicators:

• Total consumption [Wh / m²]: Total heat pump consumption during the whole day per square meter (total building surface)

• Energy savings during control [Wh / m²]: energy consumption per square meter (total building surface) difference between the reference case (normal operation) and the case with load control, for the period where the control is applied.

• Overconsumption [Wh / m²]: energy consumption per square meter (total building surface) difference between the reference case (normal operation) and the case with load control, for the period after the control.

• Recovery (%):

: indicates the quantity of energy that is consumed after the control compared to the energy saved during the control. • Degree-hours of discomfort:

DH = (1)
, where Tindoor and t is time

• Non comfortable days: percentage of the days, during the whole cooling season where there room is considered as uncomfortable (when the operative temperature is higher than 27°C)

RESULTS

Hereunder we present the average results for the whole cooling season.

Energy and power saved by load curtailment

In table 2 we present the average results for each control strategy regarding their energy effects. The TLC strategy will generate the highest demand relief with an average of 16.04 Wh/m² saved. However this strategy has the largest overconsumption after the load control (11.25 Wh/m²). On the other hand the CR strategy will have a lower demand relief but it will have also a lower "payback effect". Despite the differences between strategies the recovery phenomena seems to be almost the same for the three strategies with value around 75 %, meaning that 75 % of the saved energy during the control will be consumed once the RAC system comes on again.

Although energy effects give us an idea of the quantity of energy consumed, it does not give the load curve effects. Next, figure 5 shows the average consumption for each strategy. In figure 5 it can be seen that the MLC strategy creates two shorter and distinct demands reliefs that the TLC strategy. These peak consumptions are due to the fact that the building got warmer, making the RAC to operate at full load when it restarts. As we get further away from the restart time the over consumption will decrease gradually, so around 17-18h (5-6 p.m) (Figure 4) the consumption curves, for all strategies have almost the same value that for the reference case.

On the Figure 6, we show the avoided power (positive values), during a control, which reaches 16 W/m² for TLC and MLC. After the control period a surplus of power is needed (negative values), and a negative peak of 7 W/m² is created. Despite the CR strategy gives the lowest load reduction (4.5 W/m²), it has also the lowest impact after the control (4.5 W/m²).

Figure 6 -Average avoided power (difference between strategies and the reference case)

Comfort in load curtailment condition

In this paper the occupants comfort is analyzed for the room with the highest cooling loads. The following figure gives an example of operative temperature evolution when each strategy is applied and for the reference case. It can be seen that with the TLC strategy the temperature reaches the higher temperature peak at almost 30 °C. On the other hand the CR strategy exceeds only by little 27°C, and its time over 27°C is short.

In table 3 we show the comfort impacts of each control strategy for the whole cooling season. As expected the CR strategy has the lowest impact on occupant's thermal comfort with the lowest number of minutes per day out of comfort range (12.89 min). Although the TLC and the MLC strategies seem to have the same energetic effects, in terms of comfort impact the MLC strategy has a smaller impact in occupants comfort. The minutes per day out of comfort range for the MLC strategy is only 5 min/day less than for the TLC strategy, but when we compare the degree hours of discomfort, it can be seen that there is a reduction of 75 % when comparing MLC to TLC.

The EN 15251 specifies that a building can be non-comfortable for 5% of its occupation hours per day. In our case the non-comfortable time-period represents 24 minutes, meaning that only the CR strategy offers comfortable figures for this comfort criterion. The MLC strategy despite having an average number of minutes out of comfort range below the comfort criterion limit (24 minutes), it remains very close to the limit and the MLC standard deviation (±8.34 minutes) shows that the comfort criterion will be exceeded very easily.

Heat Wave

Next we present the results for the energy impacts when load controls are applied during a heat wave, which represent the ten consecutive hottest days for the "Nice-France" climate. The results in table 4 show that due to higher consumption of the RAC in the reference case (no control), we can achieve larger savings during a control. For TLC and MLC savings around 22.6 Wh/m² can be achieved for this building type. The CR strategy has a similar energy saving when compared with the TLC and the MLC strategies during a heat wave.

The recovery effect has an average value of 55 % for all strategies in a heat wave condition, which is lower than the annual average. The reduction of the recovery effect is explained by the fact that the RAC is almost operating at full-load for the reference case (Saturated). When the RAC is restarted only a little surplus of power can be consumed. This feature will harm even more the comfort since the room temperature will decrease slowly, resting this way more time outside of the comfort range (Table 5). The comfort impacts increases significantly when compared with the averages for the all cooling season. In the event of a heat wave no strategy can be considered as comfortable in agreement with EN 15 251.

DISCUSSION

The results give a value for the flexibility that can be provided by RAC systems in small French office buildings. Nevertheless this methodology can also be applied to residential buildings, by modifying building characteristics. We consider that the ongoing study will not change strongly the conclusions since the phenomena of air change and conduction through building shell will not be very different according to building types. However special attention should be taken, regarding special building types because they can have different internal gains.

Total load curtailment (TLC) strategy provides the largest demand relief. However large peak consumption is created immediately after the load control. The Multiple Load Curtailments (MLC) strategy presents the same problem, though instead of one single peak consumption, two peaks are created with the same magnitude. For the electric network this feature can create a problem because the peak consumption is only avoided during a limited time period and when the RAC restarts another peak consumption is created.

On the other hand the CR strategy has a lower peak consumption reduction and a lower peak consumption problem after the load control. With this strategy high current problems caused by the start of the electric motor can be avoided.

The comfort results raise the problem of acceptance by occupants of a load control, which can be a barrier to the deployment of air conditioners load controls. It was seen that only for TLC strategy the time out of comfort range was exceeded, considering our comfort range hypotheses, for the whole cooling season. However MLC strategy has an average time out of comfort very close to the defined comfort limits, its standard deviation shows that the comfort limits can be exceeded very easily.

For load control during a heat wave, it was seen that none of the studied control strategies complies with the comfort ranges. Since the DSM strategies applied during a heat have are normally emergency strategies, some kind of monetary reward could be suitable. Grignon-Massé (Grignon-Masse, 2010) has developed a methodology to evaluate costs of the discomfort in ambiance for office buildings. In his study, he analyzed the working performance decrease as a function of the interior temperature, providing this way a cost to the discomfort. However sociological studies are necessary to define the value of this reward in the residential sector.

Regarding the energy effects during a heat wave, the demand relief will be higher since the RAC systems will operate near the full-load. The energy savings will be more or less identical, around 20 Wh/m², for the three strategies. However the reconstitution is lower than for the average cooling season due to the RAC systems saturation, i.e. during the heat wave the RAC will operate at full-load for the reference case, so when the RAC is restarted only a little surplus of power can be consumed. Because of this feature, over dimensioned systems can increase the recovery problem but decrease also the comfort impact during a control.

CONCLUSIONS

In this article a methodology was developed to analyze the energy potential of RAC systems taking into account occupant's comfort. This methodology was applied to a small office buildings representative of the French stock during the whole cooling season.

The results show that a large demand relief can be achieved as high as 16 W/m², during one hour, for TLC and MLC strategies. However important peak consumption and recovery effect will appear after the load control period. This can create a network problem, when several buildings are controlled at the same time, where strategies to avoid cold load pick up2 should be taken into account.

CR strategy will have a lower demand relief and less peak consumption, however the reconstitution effect will be almost the same as for the other strategies, remaining around 70 %.

The comfort impacts have shown that as an average during the cooling season only the TLC strategy does not comply with the defined comfort limits. This means that a potential load control can be made without interfering with the occupants comfort.

The analysis during a heat wave has shown that none of the simulated strategies complies with the comfort limits, at that time monetary rewards could be needed to convince the occupants to allow any load control.

PERSPECTIVES

Future works will be made to evaluate other control strategies, and other air conditioner equipments (Chillers, Variable Refrigerant Volume air conditioning units…) which could allow different control strategies. A study about the control limits should also be made, i.e. which is allowed control time without affecting the thermal comfort?

An analysis of a load aggregation for a large number of buildings should be made in order to estimate the impact of each strategy in the electrical network.

These studies could be preceded by a more detailed analysis of the thermal comfort in a room, not analyzing only the average temperature in a room, but taking into account the different local discomfort and other variables like humidity.

Figure 1 -

 1 Figure 1 -French typical summer load curve ( Réseau de transport d'électricité)

  Figure 2 -Building plan

Figure 3 -

 3 Figure 3 -Effect of occupant's comportment over the set point temperature

Figure 4 -

 4 Figure 4 -Effects of heat pump load control

Figure 5 -

 5 Figure 5 -Average annual consumption for different strategies

Figure 7 -

 7 Figure 7 -Operative temperature evolutions for each strategy

Table 1 -

 1 EN 15251 -Temperature ranges for energy calculations

	Type of building / space Category	Acceptable Top (ºC) range
	Offices and spaces	Heating season	Cooling season
	with	~ 1,0 clo**	~0,5 clo**
	similar activity Sedentary activity ~1,2met*	I 21,0 -23,0 II 20,0 -24,0 III 19,0 -25,0	I 23,5 -25,5 II 23,0 -26,0 III 22,0 -27,0

Table 2 -

 2 Energy effects for the different control strategies

		Total consumption [Wh / m²]	Energy savings during control [Wh / m²]	Overconsumption after the control [Wh / m²]	Recovery [%]
	Total curtailment (TLC)	80.36± 24.26**	16.04 ± 5.58**	11.26 ± 3.10**	75 ± 17**
	Multiple load curtailment (MLC)	80.36± 24.34**	15.67 ± 5.64**	10.89± 3.13**	75 ± 19**
	Capacity reduction (CR)	81.50± 24.64**	9.30 ± 6.84**	5.66 ± 3.80**	71 ± 24**
	Reference case	85.14	-	-	-

**-Standard deviation -The deviation is used to characterize the results distribution

Table 3 -

 3 Comfort impact for each strategy

		Total load curtailment	Capacity reduction	Multiple load curtailments
		(TLC)	(CR)	(MLC)
	Non-comfortable days (Top>27°C)	78 %	8 %	52 %
	Number of minutes per day out of comfort range	27.66 ± 10.86**	12.89 ± 9.90**	22.38 ± 8.34**
	Annual comfort impact	39.96 DH*	2.62 DH*	9.08 DH*

**-Standard deviation -The deviation is used to characterize the results distribution

Table 4 -

 4 Energy effects for the different control strategies during Heat wave

		Energy savings	Overconsumption after the control	Recovery
		[Wh / m²]	[Wh / m²]	[%]
	Total curtailment (TLC)	22.94 ± 0.55**	12.65 ± 1.73**	55 ± 7**
	Multiple load curtailment (MLC)	22.62 ± 0.56**	12.32 ± 1.70**	54 ± 7**
	Capacity reduction (CR)	19 .34 ± 1.73**	10.34± 1.52**	53 ± 4**

**-Standard deviation -The deviation is used to characterize the results distribution

Table 5 -

 5 Comfort impact for each strategy during a heat wave

		Total load curtailment	Capacity reduction	Multiple load curtailments
		(TLC)	(CR)	(MLC)
	Non-comfortable days	10/10	9/10	10/10
	Time out-of-comfort zone per day (min)	35.23 ± 4.31**	26.57 ± 13.45**	30.42 ± 7.38**
	**-Standard deviation -The deviation is used to characterize the results distribution

-PPD -Percentage of Person Dissatisfied

Savings Overconsumption

-Cold load pick-up: After a load disturbance on the electric network, the reconnected load will be, in general, higher in size and will have a different behavior than the electric network before the disturbance.


Zheng Zhang
email: zhangzheng2035@163.com

Romain Tavenard
email: romain.tavenard@univ-rennes2.fr

Adeline Bailly
email: adeline.bailly@uhb.fr

Xiaotong Tang
email: xiaotongtang1994@163.com

Ping Tang
email: tangping@radi.ac.cn

Thomas Corpetti
email: thomas.corpetti@univ-rennes2.fr

Dynamic Time Warping Under Limited Warping Path Length

Keywords: Dynamic Time Warping (DTW), Time series, Distance measure, Warping path, Classification

Dynamic Time Warping (DTW) is probably the most popular distance measure for time series data, because it captures flexible similarities under time distortions. However, DTW has long been suffering from the pathological alignment problem, and most existing solutions, which essentially impose rigid constraints on the warping path, are likely to miss the correct alignments. A crucial observation on pathological alignment is that it always leads to an abnormally large number of links between two sequences. Based on this new observation, we propose a novel variant of DTW called LDTW, which limits the total number of links during the optimization process of DTW. LDTW not only oppresses the pathological alignment effectively, but also allows more flexibilities when measuring similarities. It is a softer constraint because we still let the optimization process of DTW decide how many links to allocate to each data point and where to put these links. In this paper, we introduce the motivation and algorithm of LDTW and we conduct a nearest neighbor classification experiment on UCR time series archive to show its performance.

Introduction

With the explosion of acquisition devices and sensors, time series is becoming more available in a very large number of applications (multi-media, physics, remote sensing, economics, clinical medicine, biodemography, etc) [START_REF] Kartikeyan | Shape description by time series[END_REF][START_REF] Petitjean | Satellite image time series analysis under time warping[END_REF] and it has the advantage of fully capturing the changing process of the observed objects [START_REF] Fu | A review on time series data mining[END_REF]. Therefore not only the information in a specific time t is important but also the temporal evolution of the observed processes is a key information to analyze. A fundamental task in time series analysis is to quantify the similarity (or dissimilarity) between two numerical sequences [START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF]. The performance of such a distance measure is critical to many data mining jobs, for instance, classification [START_REF] Huang | ECG frame classification using dynamic time warping[END_REF][START_REF] Jeong | Weighted dynamic time warping for time series classification[END_REF][START_REF] Petitjean | Dynamic Time Warping averaging of time series allows faster and more accurate classification[END_REF], clustering [START_REF] Liao | Clustering of time series data -a survey[END_REF][START_REF] Petitjean | Satellite image time series analysis under time warping[END_REF], retrieval [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF][START_REF] Bartolini | Warp: Accurate retrieval of shapes using phase of fourier descriptors and time warping distance[END_REF][START_REF] Rakthanmanon | Searching and mining trillions of time series subsequences under dynamic time warping[END_REF], etc. In the context of time series, usual distance measures such as Euclidean distance, are often not suitable [START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF], because time shifts or time distortions exist commonly and unpredictably. Dynamic Time Warping (DTW), introduced three decades ago in the context of sound processing [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF], is a widely accepted distance measure for time series [START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF]. DTW is famous for its ability to manage time distortions by realigning time series when comparing them. However, DTW still suffers from its drawback that it could lead to pathological alignments [START_REF] Keogh | Derivative dynamic time warping[END_REF]. As an example, Fig. 1 illustrates a typical pathological alignment generated by DTW, where we can observe several singularities. A singularity is a data point inside one time series that links to a large subsection of the other time series. The presence of singularities is the main signature of pathological alignment. Obviously, such an alignment is not the "correct" one we could expect and by "correct" we mean intuitively obvious "feature to feature" alignment. The quality of this alignment greatly influences the accuracy of the associated similarity measure. We cannot find accurate similarity based on pathological alignment. To cope with this problem, an intuitive solution is to limit the number of points a point can link to, and this idea has been implemented by a widely used category of constraints on DTW, which is called windowing constraint [START_REF] Rabiner | Considerations in dynamic time warping algorithms for discrete word recognition[END_REF][START_REF] Tappert | Memory and time improvements in a dynamic programming algorithm for matching speech patterns[END_REF][START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF][START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF] (will be discussed in Section 2.3). Although windowing constraint can solve pathological alignment, it also brings another problem that it takes a risk of preventing the correct alignment from being generated.

In fact, a crucial observation on pathological alignment is that if singularities happen, the total number of links (green lines in Fig. 1) between sequences must be abnormally large. Therefore, another intuitive solution, which will be explored in this paper, is to limit the total number of links between two time series. This way, we can also oppress singularities effectively. Moreover, this is a softer constraint than windowing constraint because we still let the optimization process of DTW decide how many links to allocate to each data point and where to put these links, instead of setting a rigid limit. As a result, it allows more flexibilities and avoid the risk of missing the correct alignment. We call this new method LDTW ("L" for "Limited" warping path "Length"), because it essentially limits the length of the warping path (detailed definition of warping path can be found in Section 2) during the DTW process.

LDTW originates from naive DTW. To clearly introduce LDTW and its motivations, we give a deep review of DTW and its variants in Section 2. Then in Section 3, we provide and explain the detailed implementation of the LDTW algorithm based on dynamic programming process. In the experiment section, we use LDTW as a distance measure and we report the one nearest neighbor classification error rate on UCR time series archive, in comparison with Euclidean distance, DTW without constraint and DTW with most popular windowing constraints. The experimental result shows that LDTW can achieve a better classification performance on most datasets. The upper bound of warping path length is the only parameter of LDTW and we also describe how to find the optimal upper bound by LOOCV (leave-one-out cross validation) [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF] on the training set.

Dynamic Time Warping (DTW)

Notation

Before entering into technical details, Table 1 establishes a set of notations that will be used throughout this paper. Note that the relation between current step count s and current path length l is: s = l -1, because the path length takes into account the starting point.

Table 1: List of Notations

Naive DTW

As shown in Fig. 1 and Fig. 2, DTW is alignment-based. The terminology alignment here refers to a collection of all links (or couplings) between two time series. The goal of DTW is to find an optimal alignment between two time series that achieves minimum global cost while ensuring time continuity.

The global cost is the summation of the cost between each pair of points in the alignment. So we should first define the cost between points of time series, which is denoted by the symbol δ. Euclidean distance or squared Euclidean distance is normally used as the base cost δ and we employ squared Euclidean distance as our δ in this work. The detailed definition of δ is shown in Eq. 1.

δ(a r , b c ) = (a r -b c ) 2 (1)
The optimal global cost of DTW, namely, the DTW distance, can be recursively calculated by :

D(A r , B c ) = δ(a r , b c ) + min          D(A r , B c-1 ) D(A r-1 , B c ) D(A r-1 , B c-1 ) (2)
A direct implementation of Eq. 2 leads to an exponential time complexity if we compute each subproblem from scratch. Fortunately, Eq. 2 belongs to a typical dynamic programming problem [START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Felzenszwalb | Dynamic programming and graph algorithms in computer vision[END_REF], because overlapping subproblems widely exhibit. For instance, the computation of D(A r , B c-1 ) and D(A r-1 , B c ) both require D(A r-1 , B c-1 ). The right solution here is to memorize the results of all these subproblems in a cumulative cost matrix, where the cell (r, c) keeps the result of D(A r , B c ). This way the process of DTW is converted into the process of fully filling the cumulative cost matrix in a certain order. Since the dimension of the matrix is R × C, the time complexity of the DTW algorithm is Θ(R × C). Algorithm 1

shows the basic implementation of DTW. Intuitively, we fill the bottom row and the leftmost column first, and then sweep the matrix row by row (from left to right) from bottom to top according to Eq. 2.

The last matrix cell m[R][C] holds the final DTW distance between two entire sequences.

From the perspective of cumulative cost matrix, we can only solve the calculation of DTW distance.

But if we want to gain a deep insight into DTW and the motivations of its variants, we need to review DTW closely from the perspective of warping path. As shown in Fig. 2, a warping path consists of a series of continuous matrix cells passing through the cumulative cost matrix. Each matrix cell, for example (r, c), corresponds to a link (or coupling) between two points A r and B c . In this way, the warping path W defines the alignment between two time series. Initializing the first column and row of m

3: m[1][1] = δ(a 1 , b 1 ) 4:
for r = 2 to R do

5: m[r][1] = m[r -1][1] + δ(a r , b 1 ) 6:
end for 7:

for c = 2 to C do 8: m[1][c] = m[1][c -1] + δ(a 1 , b c ) 9:
end for 10:

Fully filling the cumulative cost matrix 11:

for r = 2 to R do 12:

for c = 2 to C do

13: m[r, c] = δ(a r , b c ) + min          m[ r ][c -1] m[r -1][ c ] m[r -1][c -1] 14:
end for 15:

end for 16: The warping path is generated by certain rules and different rules lead to different distance measures.

return m[R][C]
From this point of view, the Euclidean distance can also be regarded as alignment-based and has its own rule of warping path. As shown in Fig. 2(a-b), Euclidean distance obeys a rigid one-to-one alignment and thus its warping path is only composed of the diagonal cells of the cumulative cost matrix. Implied by Eq. 2, the warping path of DTW, W , must satisfy the following three rules [START_REF] Keogh | Derivative dynamic time warping[END_REF]:

1. Boundary constraint : Fig. 2(c-d) show the alignment and the corresponding warping path generated by DTW. If we follow the idea of LDTW and make the length of warping path not longer than six, then we get Fig. 2(e-f).

w l = (1,
Given the warping path W = w 1 ... w l ... w L , DTW distance can also be defined as :

D(A, B) = min W L l=1 δ(w l ) (3)
where if w l = (r, c), then δ(w l ) = δ(a r , b c ).

This new definition reminds us that the idea of DTW is to find the optimal warping path that leads to the minimum cumulative cost.

Unlike DTW distance, the corresponding warping path does not emerge explicitly during the DTW algorithm. Although the warping path has already been determined, it hides in the cumulative cost matrix and you need another backtracking procedure to reveal it. Starting from the upper-right corner (R, C), we check the three backward neighbors of the current position at each stage. Assume (r, c) is the current position inside the cumulative cost matrix. Its three backward neighbors, (r -1, c), (r, c -1), and (r -1, c -1), will be checked and the one with the smallest value will be added to the warping path. Then we move to the new position and repeat the above steps until the lower-left corner (1, 1) is reached. Note that the warping path comes from backtracking so it is in reverse order. We need to reverse the entire path in the end to get the forward-going path.

Variants of DTW

As mentioned in the Introduction, DTW has the drawback of generating pathological alignments.

Therefore, a large number of variants of DTW have been proposed to tackle with this issue. Generally, we can classify these variants into two major categories.

The first category essentially sets constraints on the warping path. The most widely used one is called windowing constraint [START_REF] Rabiner | Considerations in dynamic time warping algorithms for discrete word recognition[END_REF][START_REF] Tappert | Memory and time improvements in a dynamic programming algorithm for matching speech patterns[END_REF][START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF][START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF]. From the perspective of alignment, it limits the number of points any point can link to, and from the perspective of warping path, it limits the allowable domain of the warping path in a specific area. Fig. 3(a-b) illustrate two most widely used shapes of constraint windows, Sakoe-Chiba band and Itakura parallelogram. In order to introduce more flexibility than a fixed authorized shape, some researchers prefer to learn an adaptive window shape from the data [START_REF] Ratanamahatana | Making time-series classification more accurate using learned constraints[END_REF][START_REF] Yu | Dynamic time warping constraint learning for large margin nearest neighbor classification[END_REF][START_REF] Candan | sDTW: computing DTW distances using locally relevant constraints based on salient feature alignments[END_REF].

For instance, Fig. 3(c) shows a window generated by learning. Other types of constraints include step pattern [START_REF] Myers | Performance tradeoffs in dynamic time warping algorithms for isolated word recognition[END_REF][START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF] and slope weighting [START_REF] Kruskall | The symmetric time warping algorithm: From continuous to discrete[END_REF][START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF].

Step pattern specifies the reachable neighbors for each step, as illustrated in Fig. 3(d). Slope weighting sets a bias toward a certain direction.

The second category replaces the feature DTW considers. Derivative DTW [START_REF] Keogh | Derivative dynamic time warping[END_REF] is a pioneer work in this category. It employs local derivative rather than the raw value, because the authors believe the trend of a data point is more important. For example, two points may have the same value, but they should be considered different if one is part of a ascending trend while the other is part of a descending trend. After Derivative DTW, many other types of features have also been considered, such as local maximum, local minimum [START_REF] Zhang | Local feature based dynamic time warping[END_REF], shape context [START_REF] Zhang | Dynamic time warping under pointwise shape context[END_REF][START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF], and interest-points [START_REF] Bailly | Bag-of-Temporal-SIFT-Words for time series classification[END_REF], etc. It is also possible to combine different distance measures to create new methods, for example, combining Derivative DTW with DTW [START_REF] Górecki | Using derivatives in time series classification[END_REF] or combining DTW with Euclidean distance [START_REF] Chen | Dtw-d: time series semi-supervised learning from a single example[END_REF].

These two categories follow different clues and they are independent with each other. The method we propose in this paper belongs to a new constraint for the warping path, so we concentrate on the methods of the first category in this paper.

LDTW : DTW with limited warping path length

Motivation

Windowing constraints limit the number of points each point can link to. It is an intuitive and straightforward solution to the pathological alignment problem. However, it takes the risk of preventing the correct alignment from being found due to its rigidity. Moreover, it is not easy to find an appropriate size or shape for the constraint window. Other types of constraints such as step patterns and slope weighting also suffer from the similar drawbacks.

In fact, another intuitive type of constraint is still open and seems promising. That is to limit the total number of links between two time series, rather than the number of links each point involves. This way, the pathological alignment will be alleviated because we prevent the tremendous amount of unreasonable links from happening by setting an upper bound to the total number of links. In other words, the number of links is not sufficient enough for singularities to form. In addition, we gain more flexibility because we still let the optimization procedure of DTW decide the number of points each point links to and where these links happen. We call DTW under this new constraint LDTW, because it actually requires us to solve DTW under Limited warping path length.

LDTW Algorithm

The calculation of LDTW is not trivial. As already shown, DTW only concentrates on the minimum cumulative cost at each stage, regardless of the warping path length. However, LDTW puts a constraint on the warping path length. Therefore, we have to consider the current length of the warping path as an additional factor. The path length is decided by the number of steps that have been taken and in order to describe the algorithm more intuitively, we sometimes use the term "step count". Note that step count = path length -1 because path length takes into account the starting point. After considering the factor of step count (denoted by s), the problem of computing DTW under a fixed warping path length can be formulated recursively into Eq. 4. Compared with the formulation of DTW shown in Eq. 2,

an extra dimension dedicated to record the current step count is added to the cumulative cost matrix. To exhaust all possible subproblems at each cell, the primary task is to ascertain the range of path lengths that each matrix cell has to consider. This is equivalent to computing the minimum and maximum steps allowed to reach a certain matrix cell. Given the rule of DTW, the minimum steps to reach a cell, for instance (r, c), is clear. It is just the larger one between r -1 and c -1, as shown in Algorithm 2. The red line in Fig. 4 illustrates one shortest path to reach (4, 3) with three steps, and the blue line illustrates one longest path with four steps. Although a path with five steps is not explicitly forbidden given the rules of DTW, we will never pick such a path because it involves orthogonal steps (shown in Fig. 4 with the dashed line) that will always lead to an extra cost. As we have mentioned, passing a matrix cell will incur a cost, therefore two orthogonal steps always have a higher cost than a diagonal step.

D(A r , B c , s) = δ(a r , b c ) + min          D(A r , B c-1 , s -1) D(A r-1 , B c , s -1)
D(A r-1 , B c-1 , s -1) (4)
The computing of maximum steps to reach a cell is a bit tricky. First let us discuss the natural situation where no limit is set on the warping path length. In this situation, if the cell is located in the bottom row (r = 1) or the left-most column (c = 1), then the maximum steps is also equal to the larger one between r -1 and c -1. Elsewhere, the maximum steps is equal to (c

-2) + 1 + (r -2) = r + c -3,
for example, we first take (c -2) horizontal steps to reach the last but one column, then a diagonal step to reach the last column, and finally (r -2) vertical steps to reach the end. Algorithm 3 summarizes the above process. Then let us discuss the situation where a upper bound L U B on the warping path length is given. In this case, the total quota of steps is L U B -1. After reaching (r, c), we still need at least max(R -r, C -c) steps to get to the end, as shown by the green line in Fig. 4. So we should reserve these steps in advance when we compute the maximum steps to reach (r, c). As a result, the current upper bound of maximum steps should be the total step quota (L U B -1) minus the reserved steps and finally we choose the smaller one between the maximum steps in natural situation and the current upper bound. Algorithm 4 summarizes the complete procedure. candidate 1 = max steps no limit(r, c) After the cumulative cost matrix is filled, the upper-right cell (R, C) contains those results to be compared. Each of them is an optimal cumulative cost when the warping path length is fixed at a different value that is smaller than the upper bound

3: candidate 2 = L U B -1 -max(R -r, C -c)
L U B . For instance, m[R][C][L -1]
is the minimum cumulative cost when the warping path length is fixed at L. We select the smallest value among these results to be the LDTW distance. Algorithm 5 summarizes the complete LDTW procedure.

The algorithms to backtrack the warping path generated by LDTW and DTW are almost identical, but we still need to notice the differences. For LDTW, the cumulative cost matrix has a third dimension to indicate the step count, so the actual step count must be provided as input, which means we should also record the final step count when calculating LDTW. For example, if m[R][C][s] is the smallest value in the final cell (R, C), then s is the final step count. In addition, remember to decrease the current step count by one after each loop.

Experiment

The main purpose of DTW and its variants is to measure the similarity between numerical sequences, so in this section we evaluate the performance of LDTW as a distance measure by nearest neighbor classification experiment 1 .

Datasets

Experiments are conducted on twenty-two datasets from the UCR time series classification archive Let L be the actual length of the warping path

1: function LDTW(A, B, L U B ) 2:
Initializing the first column and row of m

3: m[1][1][0] = δ(a 1 , b 1 ) 4:
for r = 2 to R do

5: m[r][1][r -1] = m[r -1][1][r -2] + δ(a r , b 1 ) 6:
end for 7:

for c = 2 to C do 8: m[1][c][c -1] = m[1][c -1][c -2] + δ(a 1 , b c ) 9:
end for 10:

Exhausting all subproblems at each cell of m 11:

for r = 2 to R do for s = min s to max s do The number of classes ranges from 2 (e.g. Gun-Point, Wafer, etc) to 50 (50Words). The length of time series ranges from 60 (Synthetic Control) to 637 (Lightning-2).

16: m[r, c, s] = δ(a r , b c ) + min          m[ r , c -1, s -1] m[r -1, c , s -1] m[r -1, c -1, s -1]

1NN Accuracy Report

One nearest neighbor classification (1NN for short) is strongly recommended [START_REF] Chen | The UCR Time Series Classification Archive[END_REF][START_REF] Ding | Querying and mining of time series data: experimental comparison of representations and distance measures[END_REF] to demonstrate the performance of a distance measure, because 1NN has no parameters and thus the accuracy solely depends on the distance measure. In Table 3, we report the 1NN error rate of LDTW with the best upper bound of warping path length, which is learned from the training sets by cross validation and the details will be introduced in Section 4.6. We compare LDTW with four most widely used distance measures between time series : Euclidean distance, naive DTW, DTW under best Sakoe-Chiba band (Fig. 3(a)), and DTW under best Itakura parallelogram (Fig. 3(b)). These best parameters are also learned by cross validation on the training sets. From Table 3 we can observe that LDTW get the lowest error rates on 19 out of 22 datasets and the improvements are considerable especially on OSU Leaf (5%), 50Words (6%), Fish (7%), and Car (10%). On several datasets, DTW error rates are close to zero, but LDTW still achieves better accuracies. These results demonstrate the utility of the proposed LDTW.

Fig. 6 shows the pair-wise comparison of 1NN accuracy between LDTW and its four competitors. Each dot represents a dataset, whose x-value is the 1NN accuracy generated by the method in competition and y-value is the 1NN accuracy generated by LDTW. In such a setting, a dot located above the diagonal indicates LDTW outperforms its competitor on that dataset and vice versa. A dot laying on the diagonal means a tie and the farther a dot is above the diagonal, the better LDTW performs. The number of dots above, on, and below the diagonal are depicted on each graph. For instance on Fig. 6(b), LDTW outperforms naive DTW distance on 15 out of 22 datasets, ties on 6 datasets and loses on only 1 dataset.

These graphs visually confirm the high superiority of LDTW on this variety of datasets.

Statistical Significance Test

To prove that the accuracy improvement made by LDTW is statistical significant, we first conducted omnibus Friedman test with Iman and Davenport's correction [START_REF] Friedman | The use of ranks to avoid the assumption of normality implicit in the analysis of variance[END_REF][START_REF] Iman | Approximations of the critical region of the friedman statistic[END_REF] to detect if at least one of the methods performed significantly different from the others. And then we applied pair-wise Bergmann and

Hommel test [START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF][START_REF] Garcia | An extension on"statistical comparisons of classifiers over multiple data sets"for all pairwise comparisons[END_REF] to compare methods. In both tests, we used the accuracy data shown in Table 3 and the R package "scmamp".

We passed the omnibus test with chi -squared = 12.617 and p -value = 4.388e -08 < 0.05. Then Table 4 and Table 5 list the raw and corrected results in the Bergmann and Hommel test, where each p -value represents the possibility that the two methods being compared are not significantly different.

The smaller the p -value is, the more differently two methods perform. If p -value < 0.05, we believe they are significantly different.

Pairwise comparisons of methods show that all observed differences between LDTW and well-established

The Texas Sharpshooter Fallacy

A simple but pervasive logic error in time series classification papers is that they claim a method is useful only because the method has better accuracy on some datasets. However, this is not very useful unless we can tell ahead of time that on which datasets it will be more accurate [START_REF] Batista | A complexity-invariant distance measure for time series[END_REF]. To avoid this logic error, we must prove that we can predict in advance when LDTW will have superior accuracy. We make the prediction by measuring the expected accuracy gain looking only at the training data as:

gain = accuracy LDT W accuracy competitor (5)
On training set, the accuracy is calculated by leave-one-out cross validation. If a parameter is involved, we choose the one with best accuracy. Gain values larger than one indicate that we expect LDTW will outperform the competitor on a given dataset and vice versa. We also calculated the actual accuracy gain using testing data.

To testify our prediction, we compare the expected accuracy gain with the actual accuracy gain. Fig. 7 shows the results of LDTW versus different competitors respectively. Each point represents a dataset.

The comparison results in four cases:

• TP(True Positive): In this region we predicted that LDTW would improve accuracy and we were correct. This is the most beneficial situation for LDTW. Most points fall into this region, which demonstrates the real utility of LDTW.

• TN(True Negative): In this region we correctly predicted that LDTW would decrease accuracy. In this case, we can avoid the loss of accuracy by choosing the other method.

• FN(False Negative): In this region we predicted that LDTW would decrease accuracy, but the accuracy actually improved.

• FP(False Positive): In this region we predicted that LDTW would improve accuracy, but the accuracy actually decreased. This is the only truly bad case. But only few points are in this region.

A Close Look at the Alignment Generated by LDTW

One primary motivation of LDTW is to solve the pathological alignment problem by limiting the total number of links between two time series. To visually demonstrate that LDTW can achieve that goal, let us take a close look at the alignment generated by DTW and LDTW in Fig. 8. We can observe that DTW does lead to several singularities, but they are sufficiently suppressed by LDTW. The limited number of links in LDTW are allocated more evenly while capturing those significant time shifts. space may be too large, and in that case we can set a larger step size instead of testing every candidates.

The step size we use for each dataset is shown within the parentheses in the last column of Table 3, and they are roughly proportional to the length of time series.

During LOOCV, each candidate results in a 1NN classification error count. Fig. 9 shows the LOOCV processes on three datasets and they exhibit two typical situations. Fig. 9(a) belongs to the first situation where we can observe a single optimal candidate that corresponds to the minimum error count, and in this case we choose this candidate directly. Fig. 9(b) belongs to the second situation where we can observe a continuous zone of candidates that all lead to relatively small error counts. The shapes of the error count curves in Fig. 9(a-b) are similar: the error gradually descends until it reaches the minimum point or minimum zone, and then it ascends until the end of the search space. Based on this observation, we decide to choose the median of the optimal zone as our optimal result. Note that if the size of the optimal zone is an even number, there will be two medians, and we choose the smaller one. Fig. 9(c) is a special case of the second situation, where the optimal zone happens at the end of the search space so there is no ascending phase. We still choose the median if the optimal zone happens at the end or the start. The last column of Table 3 lists the optimal value or optimal zone of L for each dataset.

To further demonstrate that selecting the median of optimal zone is fair, we compared the 1NN error rates when selecting minimum, maximum, and median/mean. The step size does not change during the entire LOOCV process and the best candidate must be an integer. In this setting, the mean is always equivalent to the median. For example, in [START_REF] Bailly | Bag-of-Temporal-SIFT-Words for time series classification[END_REF][START_REF] Batista | A complexity-invariant distance measure for time series[END_REF][START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF][START_REF] Candan | sDTW: computing DTW distances using locally relevant constraints based on salient feature alignments[END_REF][START_REF] Chen | The UCR Time Series Classification Archive[END_REF], the median is 5 and the mean is also 5. Or if the number of elements is even, for example, in [START_REF] Bailly | Bag-of-Temporal-SIFT-Words for time series classification[END_REF][START_REF] Bartolini | Warp: Accurate retrieval of shapes using phase of fourier descriptors and time warping distance[END_REF][START_REF] Batista | A complexity-invariant distance measure for time series[END_REF][START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF][START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF][START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF], the smaller median is 3 and the integer part of its mean (f loor(3.5)) is also 3. Table 6 lists the 1NN error rates of selecting minimum, maximum, or median/mean on every dataset involves an optimal zone. After a simple comparison, we can conclude that selecting the median is a better choice.

Contributions of LDTW

Apart from the interesting performances of LDTW in terms of classification, we also would like to point out the theoretical contributions of this technique. To compare running time of different methods, we repeated twenty times the leave-one-out cross validation process on five different datasets, and Fig. 10 shows the average running time. We can see that the running time basically coincide with their own time complexities. LDTW is more time-consuming than existing methods. But even when running time would be the bottleneck, we always have enough tools to deal with it, for example, parallel computing, GPGPU, FPGA, etc. In addition, it should be outlined that the structure of the LDTW algorithm is suitable for parallelization in both bottom and top level. In bottom level, for each matrix cell, we can consider situations with different step counts in parallel. In top level, one can also compute LDTW between different pairs of time series in parallel. The classification of each time series is also independent and can be performed simultaneously.

Conclusion

In

Figure 1 :

 1 Figure 1: A typical example of pathological alignment generated by DTW. The two time series here are the first two in the ECG200 training set.

A

  = a1 ... ar ... aR -a time series of length R B = b1 ... bc ... bC -another time series of length C An -a sub-sequence consists of the first n points of A m[r][c] -a R × C cumulative cost matrix used in DTW, each cell keeps only one value m[r][c][s] -a R × C cumulative cost matrix used in LDTW, each cell keeps multiple values r -row index of the cumulative cost matrix m, also the index of A c -column index of the cumulative cost matrix m, also the index of B W = w1 ... w l ... wL -a warping path of length L w l = (r, c) -a link between ar and bc s -current step count of the warping path, s starts from 0 and s = l -1 LUB -upper bound of warping path length δ() -cost between points of time series

Algorithm 1 DTW

 1 Require: A = a 1 ... a r ... a R Require: B = b 1 ... b c ... b C Let δ() be the cost between two points of different time series Let m[r][c] be the cumulative cost matrix 1: function DTW(A, B) 2:

Figure 2 :

 2 Figure 2: Alignments and their corresponding warping paths for Euclidean distance, DTW, and LDTW.

1 ) 2 .

 12 and w L = (R, C) Monotonicity constraint : Given w l = (r, c), w l+1 = (r , c ) Then r ≥ r and c ≥ c 3. Continuity constraint : Given w l = (r, c), w l+1 = (r , c ) Then r ≤ r + 1 and c ≤ c + 1 The boundary constraint indicates that the warping path should start at the lower-left corner (1, 1) and end at the top-right corner (R, C). The other two constraints define the admissible step pattern, where only the upper, upper-right, and right neighbors of the current position are reachable at each step.

Figure 3 :

 3 Figure 3: Two types of constraints on the warping path: windowing constraint (a-b-c) and step pattern (d).

  To solve Eq. 4, we rely on the similar dynamic programming principle used for DTW and the key point is to memorize the results of all subproblems. However, in DTW each cell inside the cumulative cost matrix only keeps one value, which is the global minimum cumulative cost among all possible path lengths. In contrast, in LDTW each cell has to keep multiple values, each of which is the minimum cumulative cost under a certain path length. For example, in the same cell at position (r, c), m[r][c][s] is the result of subproblem D(A r , B c , s) and m[r][c][s -1] is the result of subproblem D(A r , B c , s -1).Note that Eq. 4 only defines how to calculate DTW under a fixed warping path length. But the goal of LDTW is to find the minimum DTW distance when the warping path length is shorter that a given upper bound L U B . Therefore, a natural solution is to calculate the results for all admissible path lengths, and then choose the smallest one among those results. To achieve this solution, the first job is to ascertain the range of all admissible path lengths. According to the rules of DTW, the strict lower bound of the warping path length is max(R, C) and we already have a given upper bound L U B , so the range of path lengths is (max(R, C), ..., L U B ]. Next, we are supposed to calculate DTW under each

Figure 4 :

 4 Figure 4: Maximum and minimum steps allowed to reach a certain cell.

Algorithm 2 2 : 1 3: end function Algorithm 3 2 :Algorithm 4

 221324 Minimum steps to reach matrix cell (r, c) Require: r row index of the cumulative cost matrix Require: c column index of the cumulative cost matrix 1: function min steps(r, c) return max(r, c) -Maximum steps to reach matrix cell (r, c) with no limit on the warping path length Require: r row index of the cumulative cost matrix Require: c column index of the cumulative cost matrix 1: function max steps no limit(r, c) if r == 1 or c == 1 then end function The next job is to fill the cumulative cost matrix according to Eq. 4. For each matrix cell, we compute Maximum steps to reach matrix cell (r, c) given limited warping path length Require: r row index of the cumulative cost matrix Require: c column index of the cumulative cost matrix Require: L U B the given upper bound of warping path length Require: R length of the first time series Require: C length of the second time series 1: function max steps(r, c, L U B , R, C) 2:

4 :

 4 return min(candidate 1 , candidate 2 ) 5: end function the cumulative costs under different path lengths respectively. Fig.5illustrates how we consider different path lengths respectively. The numbers labeled on each cell show all possible step counts to reach that cell. To reach cell[START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF][START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF], the number of steps ranges from 4 to 6. When computing the 5-step cost of that cell (m[START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF][4][START_REF] Bergmann | Improvements of general multiple test procedures for redundant systems of hypotheses[END_REF]), we only consider the 4-step costs of its left, lower-left, and lower neighbors. In parallel, when computing the 4-step cost, we only consider the 3-step costs of the three neighbor cells. Fig. 5 also shows another situation that is different from DTW. In DTW we always have to consider three adjacent cells at each step, but in LDTW some cells could be inaccessible due to limited step count. As shown in Fig. 5, we cannot get to cell (3, 1) with only one step, but m[3][2][2] only considers the 1-step cost of adjacent cells, so the cell (3, 1) has to be ignored in the computation of m[3][2][2].

Figure 5 :

 5 Figure 5: LDTW considers different path lengths respectively. The numbers labeled on each cell are all possible step counts to reach that cell.

[ 9 ]

 9 , the largest online time series database widely referenced by the community in recent years. One big advantage of the archive is its diversity. It consists of mostly real time series datasets collected from various domains and some synthetic ones. For example, Gun-Point dataset is collected from video trajectory recognition field; ECG200 dataset comes from electrocardiogram signals; Swedish Leaf dataset 1 Related code of this paper is available at http://sites.univ-rennes2.fr/costel/corpetti/LDTW/LDTW CODE.zip Algorithm 5 Complete LDTW procedure Require: A = a 1 ... a r ... a R Require: B = b 1 ... b c ... b C Require: L U B the given upper bound of warping path length Let δ() be the cost between two points of different time series Let m[r][c][s] be the cumulative cost matrix

12 : for c = 2 to C do 13 :

 1213 min s = min steps(r, c) 14: max s = max steps(r, c, L U B , R, C) 15:

20 : 21 :for s = min s f inal to max s f inal do 26 :

 202126 Selecting the best result and the corresponding path length in the last cell (R, C)min s f inal = max(R, C) if m[R][C][s] < ldtw then 27: ldtw = m[R][C][s] end functionis converted from the contours of fifteen different leaf species; while Synthetic Control, CBF, and Two-Patterns datasets are synthetic. Table2summarizes the basic information of all datasets used in the experiments, such as the number of classes, size of training set, size of testing set and time series length.

Figure 6 :

 6 Figure 6: 1NN accuracy of LDTW versus the four competitors respectively.

  vs DTW (Sakoe-Chiba band) (d) vs DTW (Itakura parallelogram)

Figure 7 :

 7 Figure 7: Expected accuracy gain calculated on training data versus actual accuracy gain

4. 6 .Figure 8 :

 68 Figure 8: Comparison between alignments generated by DTW and LDTW. The two time series here are the first two in the CBF training set.

Figure 9 :

 9 Figure 9: Three typical error count curves in LOOCV (leave-one-out cross validation) processes.

First

  , although LDTW and some other existing variants of DTW, for instance, windowing constraint, both set constraints on the warping path, LDTW prefers a global one rather than local ones. As already shown in Section 2.3, given a typical windowing constraint such as Sakoe-Chiba band or Itakura parallelogram, a point can only link to the points temporally close enough to it, so windowing constraint essentially limits the local behavior of the warping path. Other variants of DTW such as step pattern and slope weighting also set constraints locally. In contrast, LDTW limits the total length of the warping path and this obviously belongs to a global constraint. As shown by the experimental results, LDTW demonstrates that global constraints may be more efficient than local constraints. We believe LDTW opens up a new direction of imposing different types of global constraints on DTW. It may take more advantages of the intrinsic optimization process of DTW and be more flexible.Second, LDTW makes us follow a similar pattern of time distortions when comparing different pairs of time series. Excluding noises, the dissimilarities between time series mainly come from two sources: time distortions and intrinsic differences. The ideal distance measure only captures intrinsic dissimilarities while fully eliminates the dissimilarities caused by time distortions. But the definition of DTW implies that DTW tends to attribute all dissimilarities between time series to time distortions. If a set of time series come from the same sensor or are issued from the same underlying mechanism, there is a high possibility that they share a common pattern of time distortions. The extracted warping path reflects the pattern of time distortions among time series. To some degree, limiting or even fixing the length of warping path ensure the distance measure obeys this common pattern of time distortions, and thus make comparisons between sequences mainly affected by intrinsic dissimilarities rather than time distortions.4.8. Time Complexity and ParallelizationIn this subsection, we discuss the time complexity and running time of LDTW in comparison with other well-established methods. Let R and C be the lengths of two time series, and L U B be the parameter (upper bound of warping path length) in LDTW. Euclidean distance can only deal with time series with the same length, and its time complexity is Θ(R or C), R = C. As DTW has to consider all cells in the warping matrix, its time complexity is Θ(R • C). SC-DTW and IP-DTW set constraints on the warping path, so they only need to consider part of the warping matrix and their time complexities are both O(R • C). Similarly to DTW, LDTW considers the entire warping matrix. However, for each matrix cell, LDTW also considers situations under different path lengths. As a result, LDTW has a time complexity of O(L U B • R • C).

Figure 10 :

 10 Figure 10: Running time of the leave-one-out cross validation process on five different datasets.

  this paper, we have proposed a new variant of DTW, entitled LDTW, which limits the length of the warping path in DTW. LDTW sets a global and softer constraint instead of those local and rigid ones used by many existing variants of DTW. LDTW thus can efficiently alleviate the pathological alignment while keeping the flexibility of DTW. To implement LDTW, we have extended the original DTW algorithm by adding a step count dimension to the cumulative cost matrix. All details associated with the LDTW algorithm have been introduced in this paper. To demonstrate the performance of LDTW, we have conducted a series of one nearest neighbor classification experiments on the newest UCR time series archive. Comparing with other most popular time series distance measures, LDTW shows better accuracy on most datasets. We therefore believe this work opens up a new direction of imposing global constraints on DTW instead local ones.

Table 2 :

 2 Basic Information of Datasets

		Number of	Size of	Size of	Time series
	Dataset
		classes	training set testing set	length
	Synthetic Control	6	300	300	60
	Gun-Point	2	50	150	150
	CBF	3	30	900	128
	Face(all)	14	560	1690	131
	OSU Leaf	6	200	242	427
	Swedish Leaf	15	500	625	128
	50Words	50	450	455	270
	Trace	4	100	100	275
	Two Patterns	4	1000	4000	128
	Wafer	2	1000	6164	152
	Face(four)	4	24	88	350
	Lightning-2	2	60	61	637
	Lightning-7	7	70	73	319
	ECG200	2	100	100	96
	Adiac	37	390	391	176
	Yoga	2	300	3000	426
	Fish	7	175	175	463
	Plane	7	105	105	144
	Car	4	60	60	577
	Beef	5	30	30	470
	Coffee	2	28	28	286
	OliveOil	4	30	30	570

Table 4 :

 4 Raw p -value in Bergmann and Hommel test

		Euclidean	DTW	SC-DTW IP-DTW	LDTW
	Euclidean	N/A	0.095205	0.006580	0.000419	5.49e-08
	DTW	0.095205	N/A	0.294266	0.062991 0.000166
	SC-DTW	0.006580	0.294266	N/A	0.417685 0.006580
	IP-DTW	0.000419	0.062991	0.417685	N/A	0.056530
	LDTW	5.49e-08	0.000166	0.006580	0.056530	N/A

Table 5 :

 5 Corrected p -value in Bergmann and Hommel test

		Euclidean	DTW	SC-DTW IP-DTW	LDTW
	Euclidean	N/A	0.190409	0.026321	0.002514	5.49e-07
	DTW	0.190409	N/A	0.294266	0.188974 0.000995
	SC-DTW	0.026321	0.294266	N/A	0.417685 0.026321
	IP-DTW	0.002514	0.188974	0.417685	N/A	0.113061
	LDTW	5.49e-07	0.000995	0.026321	0.113061	N/A

Table 6 :

 6 1NN Error Rates of Selecting Minimum, Maximum, and Median/Mean of Optimal Zone

	Dataset	Optimal Zone Minimum Maximum Median/Mean
	Synthetic Control	71-74	0.0100	0.0067	0.0067
	Gun-Point	156-161	0.0267	0.0200	0.0200
	CBF	142-254	0.0022	0.0033	0.0033
	Trace	286-546	0.0100	0.0000	0.0000
	Two Patterns	139-254	0.0000	0.0000	0.0000
	Wafer	154-157	0.0032	0.0034	0.0028
	Lightning-2	748-818	0.0983	0.0983	0.0983
	Lightning-7	335-380	0.2740	0.2603	0.2055
	ECG	97-99	0.1200	0.1100	0.1200
	Yoga	447-487	0.1423	0.1483	0.1490
	Plane	151-286	0.0000	0.0000	0.0000
	Coffee	287-327	0.0000	0.0357	0.0000
	OliveOil	571-601	0.1333	0.1667	0.1333

Acknowledgments

The authors would like to thank Prof. Eamonn Keogh and all the people who have contributed to the UCR time series classification archive for their selfless work. We also thank the anonymous reviewers for their valuable advice. This work has been supported by Major Project of High Resolution Earth Observation System of China (Grant No.03-Y20A04-9001-15/16), the CNES TOSCA-VEGIDAR Program, and CAS-CNRS Joint Doctoral Promotion Program.

a r is the percentage of time series length.

b d is the half length of the shorter diagonal of the parallelogram.

c this column shows the optimal value or optimal zone of the upper bound of warping path length found by leave-one-out cross validation.