Emmanuel Ramasso Femto-St
email: emmanuel.ramasso@femto-st.fr

Investigating Computational Geometry for Failure Prognostics in Presence of Imprecise Health Indicator: Results and Comparisons on C-MAPSS Datasets

Prognostics and Health Management (PHM) is a multidisciplinary field aiming at maintaining physical systems in their optimal functioning conditions. The system under study is assumed to be monitored by sensors from which are obtained measurements reflecting the system's health state. A health index (HI) is estimated to feed a data-driven PHM solution developed to predict the remaining useful life (RUL). In this paper, the values taken by an HI are assumed imprecise (IHI). An IHI is interpreted as a planar figure called polygon and a case-based reasoning (CBR) approach is adapted to estimate the RUL. This adaptation makes use of computational geometry tools in order to estimate the nearest cases to a given testing instance. The proposed algorithm called RULCLIP-PER is assessed and compared on datasets generated by the NASA's turbofan simulator (C-MAPSS) including the four turbofan testing datasets and the two testing datasets of the PHM'08 data challenge. These datasets represent 1360 testing instances and cover different realistic and difficult cases considering operating conditions and fault modes with unknown characteristics. The problem of feature selection, health index estimation, RUL fusion and ensembles are also tackled. The proposed algorithm is shown to be efficient with few parameter tuning on all datasets.

INTRODUCTION

Knowledge-based systems and Case-Based Reasoning approaches (CBR) have appeared as suitable tools for data-driven Prognostics and Health Management (PHM) [START_REF] Saxena | Integrated diagnosis and prognosis architecture for fleet vehicles using dynamic case-based reasoning[END_REF]; T. [START_REF] Wang | A similaritybased prognostics approach for remaining useful life estimation of engineered systems[END_REF][START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF][START_REF] Ramasso | Joint prediction of observations and states in time-series based on belief functions[END_REF]. In Emmanuel Ramasso et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 United States License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

CBR, historical instances of the system -with condition data and known failure time -are used to create a library of degradation models or health indices. Then, for a test instance, the similarity with the degradation models is evaluated generating a set of Remaining Useful Life (RUL) estimates which are finally aggregated.

The required assumptions for CBR implementation are limited, the main issues consisting in, on the one hand, the choice of an appropriate similarity measure and, on the other hand, the selection of the relevant training instances. CBR approaches are also flexible since it is simple to incorporate quantitative and qualitative pieces of knowledge such as measurements and expertise.

We consider applications for which the noise due to various sources, such as operational conditions or fault modes, can not be well characterised and where filtering may change the meaning of the health index. We assume that the health index can not be well defined by a single real value but only by Imprecise Health Index (IHI). To fix ideas, an illustration taken from the turbofan engine dataset [START_REF] Saxena | Damage propagation modeling for aircraft engine runto-failure simulation[END_REF]) (used and detailed in experiments) is given in Figure 1. The figure pictorially represents the IHI taken from the fourth dataset (made of two fault modes and six operating conditions) for the 8th training data (P 1 ), the 100th training data (P 2 ) and the 1st testing data (P 3 ) of this dataset. As depicted, fault modes may generate

• sudden changes in wear (e.g in P 1 , t ∈ [225,275]) that may increase the lifetime of the unit. It may be due to both fault modes and operating conditions, for example a drastic decrease of speed to cope with mechanical incidents or meteorological phenomenons.

• Unexpected changes in the trend, such as increasing instead of decreasing (e.g. P 2 , t > 125) that may disturb the algorithm. It may be due to component failures such as sensors or electronics. • Sudden bursts characterised by low signal-to-noise ratio (SNR) on a possibly short duration which deeply affect the HI (e.g. on P 3 with t ∈ [10,75]) that may affect the lifetime accordingly to the fault type which is generally unknown. The representation and the propagation of imprecision (or uncertainty) is of paramount importance in engineering analyses [START_REF] Vachtsevanos | Intelligent fault diagnosis and prognosis for engineering systems[END_REF][START_REF] Orchard | Advances in uncertainty representation and management for particle filtering applied to prognostics[END_REF]Beer, Ferson, & Kreinovich, 2013). Several mathematical theories [START_REF] Klir | Uncertainty-based information. elements of generalized information theory[END_REF] have been used in prognostics such as probability theory (including Bayesian approaches) [START_REF] Peng | Integrated fatigue damage diagnosis and prognosis under uncertainties[END_REF], set-membership approaches (including fuzzy-sets) [START_REF] Chen | Machine condition prediction based on adaptive neuro-fuzzy and high-order particle filtering[END_REF][START_REF] El-Koujok | Reducing arbitrary choices in model building for prognostics: An approach by applying parsimony principle on an evolving neuro-fuzzy system[END_REF] and Dempster-Shafer's theory of belief functions [START_REF] Serir | E2GKpro: An evidential evolving multi-modeling approach for system behavior prediction with applications[END_REF]Ramasso et al., 2013). Facing imprecision in HIs for prognostics is thus not new but the way to handle it can be considered differently.

We assume 1-D health index to be available but obtained from noisy measurements. The data points are supposed to represent vertices of a simple planar polygon. The IHI is thus a polygon-shaped health index represented by a planar figure. Three polygons are depicted in Figure 1. Using computational geometry tools, a prognostics method is proposed that handles IHI without knowing nor estimating the noise properties. The method is based on CBR for which a similarity measure adapted to IHI and polygon is developed. The set of cases is made of training instances represented by polygons and the similarity with a testing instance recorded on the in-service system is made dependent on the degree of intersection between both training and testing polygon instances.

The prognostics algorithm introduced is called "RULCLIP-PER" (Remaining Useful Life estimation based on impreCise heaLth Index modeled by Planar Polygons and similarity-basEd Reasoning").

The next Section is dedicated to the presentation of a methodology to build imprecise health index and perform prognostics. The methodology is then applied on C-MAPSS datasets.

PROGNOSTICS BASED ON IMPRECISE HEALTH INDEX: A CBR APPROACH

A health index (HI) takes the form of a 1-dimensional realvalued signal H = [x 1 x 2 . . . x j . . . x T ] T , x j ∈ R obtained at some instants t 1 , t 2 . . . t T .

Polygon-shaped representation of IHI

An IHI is defined as a polygon where each vertex is represented by a data point estimated from the original HI. A polygon (representing an IHI) is thus defined as a set of pairs (x j , t j ) made of HI values x j at time index t j .

The upper envelope of H denoted S is defined by

S = {(x j , t j )|x j ≥ xj } ∪ {(x j-1 , t j )|x j < xj } , (1)
meaning that, for a given data point j, if the HI value x j at time t j is greater than the filtered value xj then the upper envelope is equal to the HI value, otherwise it takes its previous value. The lower envelope I is defined similarly by

I = {(x j , t j )|x j < xj } ∪ {(x j-1 , t j )|x j ≥ xj } . (2)
The ordered pairs of vertices listed counterclockwise represents a bounding closed polygonal chain that separates the plane into two regions. The word "polygon" refers to a plane figure bounded by the closed path defined as:

P = (x 1 , t 1 ) S , (x 2 , t 2 ) S . . . (x j , t j ) S . . . (x T , t T ) S , (x T , t T ) I , (x T -1 , t T -1 ) I . . . (x 1 , t 1 ) I , (x 1 , t 1 ) S
(3) with (x j , t j ) S ∈ S and (x j , t j ) I ∈ I. To close the polygon, the first and last vertices are the same. The pairs of vertices define a finite sequence of straight line segments representing the polygon.

More specifically, a polygon is a region of the plane enclosed by a simple cycle of straight line segments where nonadjacent segments do not intersect and two adjacent segments intersect only at their common endpoint [START_REF] Rosen | Handbook of discrete and computational geometry[END_REF]. However, the second part of the definition of the bounds may generate some segment intersections. These inconsistencies can be corrected easily by exchanging the corresponding values of the lower and upper bounds when an intersection is detected. When consistent bounds are obtained, the polygon is made of nonintersecting line segments which characterise a Jordan's simple closed curve also called simple polygon [START_REF] Filippov | An elementary proof of Jordan's theorem[END_REF]. This category of polygon enables one to apply some standard algorithms from Computational Geometry [START_REF] Rigaux | Spatial databases with application to gis[END_REF][START_REF] Rosen | Handbook of discrete and computational geometry[END_REF][START_REF] Longley | Geospatial analysis : A comprehensive guide to principles, techniques and software tools[END_REF]. Note that some of the most efficient algorithms for operations on polygons can manage self-intersections [START_REF] Vatti | A generic solution to polygon clipping[END_REF][START_REF] Greiner | Efficient clipping of arbitrary polygons[END_REF]) but these inconsistencies generally increase time-consumption.

CBR approach for prognostics based on IHI

Training dataset

We assume the training dataset to be composed of N training instances:

L = {P i , K i } N i=1 ( 4
)
where P i is the ith polygon attached to the ith imprecise health index H i and K i = [y 1 y 2 . . . y j . . . y T ] T , y j ∈ N represents a discrete-valued signal reflecting a system's state. The component K i may be useful in some applications where the system can be described by means of latent variables [START_REF] Ramasso | Making use of partial knowledge about hidden states in hidden Markov models: an approach based on belief functions[END_REF][START_REF] Javed | Novel failure prognostics approach with dynamic thresholds for machine degradation[END_REF]. In that case, K i may represent a partial knowledge about the state. For example, in (Ramasso et al., 2013), partial knowledge was encoded by belief functions to express imprecision and uncertainty about the states.

Determining the nearest case

A testing instance takes the form of a health index H * from which the envelopes can be estimated as explained in the previous paragraph, leading to the polygon representation P * . As in usual CBR approaches for prognostics (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]Ramasso et al., 2013), the goal is to sort the training instances with respect to their similarity to the testing instance. However, since the training instances are made of polygons, the usual Euclidean distance is not adapted. We propose the following similarity measure.

Getting inspired from the Computer Vision community [START_REF] Powers | Evaluation: From precision, recall and F-factor to ROC, informedness, markedness & correlation[END_REF], recall, precision and F β indices are used to quantify the relevance of a training instance compared to the testing one. Precision represents the fraction of the retrieved instance that is relevant, while recall is the fraction of the relevant instance that is retrieved. The F β is an harmonic mean which gives equal weight to recall and precision when β = 1. Note that the three indices are normalised into [0, 1].

More precisely, for the ith training instance:

1. Estimate the area of the intersection between the polygon P i and P * :

A ∩ = Area (P i ∩ P * ) (5)
2. Compute the "recall":

R ec = A ∩ A i (6)
3. Compute the "precision":

P rec = A ∩ A * (7)
4. Compute the "F β,i ", in particular for β = 1, characterizing the similarity with the ith training instance:

F 1,i = 2 R ec • P rec R ec + P rec ( 8
)
where A i , A * , A ∩ represent the area of polygons P i , P * and of their intersection respectively.

Example 1 An illustration of intersection is given in Figure 2 where the darkest polygon represents a training instance and the two other polygons are testing ones. The whitest polygon is within the testing instance meaning that the precision is high, but the recall is pretty low since it covers only a small part of the testing instance. On the opposite, the third polygon covers entirely the testing instance leading to a high recall but its spread decreases the precision. Practically, intersection construction is the main difficulty and was tackled quite recently in computational geometry for arbitrary planar polygons. It consists in determining the region of geometric intersection which can be performed in three phases [START_REF] Rosen | Handbook of discrete and computational geometry[END_REF]) (Chap. 38):

1. Compute the intersection between the boundaries of the objects using the linearithmic plane sweep algorithm (Bentley & Ottmann, 1979); 2. If the boundaries do not intersect then determine whether one object is nested within the other; 3. If the boundaries do intersect then classify the resulting boundary fragments gathered to create the final intersec-tion region [START_REF] Margalit | An algorithm for computing the union, intersection or difference of two polygons[END_REF]Chazelle & Edelsbrunner, 1992), which can be performed in linearithmic time. Regularized Boolean operations ensure the closure of the interior of the set-theoretic intersection.

In this paper, the Vatti's algorithm [START_REF] Vatti | A generic solution to polygon clipping[END_REF] has been used because it is generic and can manage most of pratical cases. Several implementations of this algorithm have been proposed, especially in [START_REF] Greiner | Efficient clipping of arbitrary polygons[END_REF] which was shown to be particularly efficient.

Estimating the Remaining Useful Life (RUL)

The F 1 measure is used to sort the N training polygon instances in descending order:

P (1) > P (2) • • • > P (j) • • • > P (N ) so that P (1)
is the closest instance to the testing one and P (N ) the farthest one. The index (i) in P (i) represents a reordering and the symbol > in P (i) > P (j) means that the ith polygon is more similar to the testing instance that the jth.

CBR assumes that a limited number of instances, say K, are enough to approximate the testing instance. The K closest training instances can then be combined to estimate the RUL.

The length of a training instance minus the length of a testing instance provides an estimation of the RUL (Figure 3). Given the definition of a polygon (Section 2.1) and of the training dataset (Eq. 4), the length of both the training and testing polygon instances is given by T i and T * respectively. Therefore, the estimated RUL is given by

RUL = T i -T * . (9)
Example 2 Two polygons are illustrated in Figure 3, one coming from the training dataset #1 (the tenth instance) and one from the testing dataset #1 (the first instance). Since T 1 = 222 and T * = 31, the estimated RUL is 191 time-units. Each closest training instance P (i) can be accompanied by a state sequence K (i) so that K estimations of the RUL, denoted RUL K , can be obtained from the state sequences in addition to the ones obtained with P (i) and denoted RUL P .

Using K (i) , the last transition in the sequence is supposed to represent a jump of the system to a faulty state. This assumption relies on the fact that the last part of a training instance represents the system's end-of-life (Ramasso et al., 2013;[START_REF] Ramasso | RUL estimation by classification of predictions: an approach based on a neuro-fuzzy system and theory of belief functions[END_REF][START_REF] Javed | Novel failure prognostics approach with dynamic thresholds for machine degradation[END_REF].

The 2K estimations of the RUL can then be pooled in one set: RUL PK = { RUL P , RUL K } and an information fusion process can then be performed to combine these partial RUL estimates. According to the application, the fusion rule can be adapted [START_REF] Kuncheva | Combining pattern classifiers: Methods and algorithms[END_REF].

A plot chart of RULCLIPPER algorithm is depicted in Figure 4. Some of the elements will be illustrated in the next section dedicated to experiments.

EXPERIMENTS: METHOD

RULCLIPPER is tested on the datasets obtained from the turbofan engine degradation simulator (Saxena, Goebel, et al., 2008). Before presenting results, several details about the datasets have to be presented, in particular how to select the features and how to compute the health index.

Turbofan engine degradation simulator

The simulation model (Saxena, Goebel, et al., 2008) was built on the Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) developed at NASA Army Research Lab., able to simulate the operation of an engine model of the 90.000 lb thrust class. A total of 21 output variables were recorded to simulate sensor measurements to the system. Another 3 variables representing the engine operating conditions were recorded, namely altitude (kilo feet), Mach number (speed) and Throttle Resolver Angle (TRA) value which is the angular deflection of the pilot's power lever having a range from 20% to 100%. In the sequel, references to variables are made by using their column position in the data files as provided on the data repository of the Prognostics Center of Excellence website: it begins by number 6 and finishes to 26 (see (Saxena, Goebel, et al., 2008) for details).

Datasets

Six datasets generated from independent simulation experiments were provided, each with some specificities (Saxena, Goebel, et al., 2008).

Datasets #1 and #2 include only one fault modes (HPC degradation) while datasets #3 and #4 include two (HPC degradation and fan degradation). Datasets #1 and #3 include a single operational condition against six for datasets #2 and #4.

Dataset #4 represents the most complex case study. Datasets #5 T (semi-final testing dataset) and #5 V (final validation dataset) were generated for the 2008's PHM data challenge with one fault mode and six operating conditions. The two last datasets have common training instances. A summary of the six datasets are shown in Table 1 according to information taken from (Saxena, Goebel, et al., 2008).

Each dataset is divided into training and testing subsets. The training set includes instances with complete run-to-failure data (to develop life prediction models), and the actual failure mode for training instances in #3 and #4 is not labeled.

The testing datasets include instances with data up to a certain cycle and are used for RUL estimation and algorithm performance evaluation.

The testing instances are also simulated run-to-failure and only an earlier portion of the history is provided. The actual life (RUL) of the testing instances are known only for datasets #1, #2, #3 and #4 and can only be used for testing algorithm. For datasets #5 T and #5 V , results have to be uploaded to the data repository: uploading is allowed only once a day for #5 T whereas only a single try is possible for dataset #5 V .

The validation can be performed by many performance measures [START_REF] Saxena | Metrics for evaluating performance of prognostic techniques[END_REF] among which accuracybased measures such as the timeliness, also called scoring function in the sequel since it has been used in the data challenge to sort participants algorithm. The review of papers using the C-MAPSS datasets show that the timeliness was the most used performance measure (about 30% of papers). Note that, for datasets #5 T and #5 V , this performance measure is returned for each submission by the data challenge chairs.

For comparison purpose, the scoring function is also used in this paper with the same parameters as in the challenge:

S = N n=1 S n (10a)
S n = e -dn/13 -1, d n ≤ 0 e dn/10 -1, d n > 0 , n = 1 . . . N (10b)

d n = estimated RUL -true RUL (10c)
This function, that assigns higher penalty to late predictions, has to be minimised. In addition to the scoring function (computed for all datasets), a second performance measure was used (on datasets #1 to #4 for which we know the RUL) called accuracy measure A that evaluates the percentage of testing instances for which the RUL estimate is within the interval [-13, +10] around the true RUL [START_REF] Saxena | Metrics for evaluating performance of prognostic techniques[END_REF]. These values are the same as the scoring function and was used in several papers such as (Ramasso et al., 2013) for dataset #1.

Related results on C-MAPSS

For comparison purpose, results of predictions from other researchers (as exhaustive as possible) on these datasets are summarised below for each dataset. Note that some authors also used the simulator to create their own datasets [START_REF] Sarkar | Data-driven fault detection in aircraft engines with noisy sensor measurements[END_REF][START_REF] Zein-Sabatto | Statistical approach to online prognostics of turbine engine components[END_REF]Al-Salah, Zein-Sabatto, & Bodruzzaman, 2012). References have been put on the NASA PCOE website.

To our knowledge, the full testing dataset of #1 was only used in two papers: In [START_REF] Liu | A data-level fu-sion model for developing composite health indices for degradation modeling and prognostic analysis[END_REF] where the authors reported results by using an average error between true RUL and prediction; and to evaluate the EVIPRO algorithm in (Ramasso et al., 2013) where the performance was assessed by using the accuracy measure which was equal to 53% on the testing dataset #1. The full testing datasets of #2, #3, #4 were not used in the past (only a few instances were considered in a few papers).

Datasets C-MAPSS DATASETS TURBFOFAN CHALLENGE #1 #2 #3 #4 #5 T #5 V Nb. of
Testing datasets #5 T (corresponding to a "semi-final" testing dataset) and #5 V (corresponding to the "final" validation dataset) represent datasets for which the true RULs is not known. These datasets were used in many papers summarised in Table 2 (for published work after 2008) and in Table 3 (for results of challengers during the competition in 2008). The complete review of scores on these datasets were found on the web or obtained by request to the conference chairs. In Table 3, methods (1), ( 2) and (3) were published in (T. [START_REF] Wang | A similaritybased prognostics approach for remaining useful life estimation of engineered systems[END_REF], [START_REF] Heimes | Recurrent neural networks for remaining useful life estimation[END_REF] and [START_REF] Peel | Data driven prognostics using a Kalman filter ensemble of neural network models[END_REF] respectively.

It can be observed that no score has been mentioned in the literature on the final validation dataset #5 V since 2008, whereas the semi-final testing dataset #5 T was used in several papers. The final dataset is complex and the performances obtained by the challengers are high. According to our knowledge, good performances (in terms of scoring) can be obtained on the final dataset only if the algorithm is robust. Indeed, a few important mistakes (too late or too early predictions) can lead to bad scores. This was also observed with RULCLIPPER on the other datasets. Robustness can be evaluated by computing several PHM metrics [START_REF] Saxena | Metrics for evaluating performance of prognostic techniques[END_REF] as proposed in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF].

Therefore, the generalisation capability of the algorithm should be ensured before trying the final dataset. This is illustrated in Tables 2-3 and Figure 6 which depict the scores obtained on the semi-final dataset #5 T and on the final dataset #5 V . Some algorithms exhibited very low score on #5 T (made of 218 instances), whereas a relatively poor score was obtained on the final dataset. The winner obtained 737 on #5 T (according to the conference chairs) which is not the best score, but only 5636 on the final dataset #5 V .

Algo. (pseudo.) #5 T #5 V RULCLIPPER 752 11572 SBL (P. Wang, Youn, & Hu, 2012) 1139 n.a. DW [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 1334 n.a. OW [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 1349 n.a. MLP [START_REF] Riad | Evaluation of neural networks in the subject of prognostics as compared to linear regression model[END_REF] 1540 n.a. AW [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 1863 n.a. SVM-SBI [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 2047 n.a. RVM-SBI [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 2230 n.a. EXP-SBI [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 2282 n.a. GPM3 [START_REF] Coble | Merging data sources to predict remaining useful life -an automated method to identify prognostic parameters[END_REF] 2500 n.a. RNN [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 4390 n.a. REG2 [START_REF] Riad | Evaluation of neural networks in the subject of prognostics as compared to linear regression model[END_REF] 6877 n.a. GPM2B [START_REF] Coble | Merging data sources to predict remaining useful life -an automated method to identify prognostic parameters[END_REF] 19200 n.a. GPM2 [START_REF] Coble | Merging data sources to predict remaining useful life -an automated method to identify prognostic parameters[END_REF] 20600 n.a. GPM1 [START_REF] Coble | Merging data sources to predict remaining useful life -an automated method to identify prognostic parameters[END_REF] 22500 n.a. QUAD [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF] 53846 n.a.

Table 2. Performance of the state-of-the-art approaches on #5 T (semi-final dataset) and #5 V (final dataset) after 2008 (published work).

Note that some papers using the datasets of the data challenge are not mentioned in the table because error measures (accuracy-based) were given and that is not possible by using the original testing datasets for which the true RULs are not known: testing errors are not possible on testing datasets #5 T and 5 V , but only on the training dataset. This rule (defined from 2008 to 2014) may change in the near future so that other metrics (in addition to the scoring function) could be obtained on demand to the data challenge chairs.

Priors about the datasets

Some rules were used to improve prognostics on these datasets, some have been proposed in previous papers:

R1: The first rule is related to the fact that, according to (Saxena, Goebel, et al., 2008), the maximum RUL in testing instances for #5 T was greater than 10 and lower than 150 time-units, while being greater than 6 and greater than 190 in testing instances for #5 V . Moreover, most of previous approaches agreed on limiting the RUL estimates around 135 (depending on papers (T. where RUL was the estimated RUL. R2: The difference between 1 and the average of the first 5% of an instance was used as an offset to compel the health index (HI) to begin around 1. Even though the health index function (Eq. 11) already compels it, there are some cases, in particular for #2, #3 and #4, for which the health index was strongly disturbed by fault modes and operating R3: To limit the impact of fault modes (datasets #3 and #4), a detection of the monotonicity [START_REF] Coble | Merging data sources to predict remaining useful life -an automated method to identify prognostic parameters[END_REF] is performed. When the testing instance is less than the half of the training instance and if more than 25 consecutive samples are above the training instance, then the training instance is not taken into account. This simple rule was applied on all datasets considered (even without fault modes or without operating conditions). R4: To decrease the risk of overpredictions, the sequence of states K were made of two states, the second state being activated only 15 samples before the end-of-life. This setting similar to [START_REF] Ramasso | RUL estimation by classification of predictions: an approach based on a neuro-fuzzy system and theory of belief functions[END_REF], was the same for all tests and all datasets.

Local/global health index estimation

To reflect a real-world and practical cases, the health indices (HI) for both training and testing datasets were not given by the organisers (Saxena, Goebel, et al., 2008). An adaptation of the approach proposed in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF] is presented below to estimate the HI for each instance. These HIs (highly corrupted by noise) are the basis of the proposed methodology described in previous sections (Fig. 4).

The set of features for the ith unit is

X i = [x 1 x 2 . . . x t . . . x Ti ] T
where x t = [x t,1 x t,2 . . . x t,m . . . x t,q ] is the q-dimensional feature vector at t (composed of sensor measurements), and u t is the vector of operating conditions at t. The operational conditions variables can be clustered into a finite number of operating regimes (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]. Crisp outputs are obtained such that the current regime at time t, C t , is precisely known. Then, for samples (u t , x t ) collected at early age of the system, e.g. t < σ 1 , the health index attached to the ith training unit is HI(x t , θ θ θ p ) = 1, where the set of parameters θ θ θ p depends on the model used to link regimes and sensor measurements.

At late age of the system, e.g. t > σ 2 , the corresponding output is HI(x t , θ θ θ p ) = 0. In (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], the author used only the data at t > σ 2 and t < σ 1 in addition to 6 models (one for each operating mode) built on all data. In comparison, we propose to make use of samples between σ 1 and σ 2 while building a local model for each operating mode in each training instance. Moreover, we have used one HI for each training instance while in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]) a global HI model was estimated using all instances.

The corresponding output of the index is set to

ĤI i (x t , θ θ θ p ) ≡ 1 -exp log(0.05) 0.95 • T i • t , t ∈ [σ 1 , σ 2 ]. (11)
This function allows to compel the health index to be globally decreasing, from 1 (healthy) to 0 (faulty). As proposed in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], σ 1 = T i • 5% and σ 2 = T i • 95% where T i is the length of the ith training instance. We used local linear models for multi-regime health assessment so that θ θ θ p i = [θ p i,0 θ p i,1 . . . θ p i,q ] represents the parameters of a linear model defined conditionnally to the pth regime. The health index at time t given the pth regime can be estimated as

HI i (x t , θ θ θ p i ) = θ p i,0 + q n=1 θ p i,n • x t,n (12)
where θ θ θ p can be estimated by standard least-squares algorithms. In experiments, in case the estimation of HI is performed by considering the three operating conditions, then it will be called a local approach (Fig. 4) and global otherwise. HIs are then transformed into IHIs as proposed in previous sections (Fig. 4).

Information fusion for improved RUL estimation

The first family of rules is a combination of minimum and maximum RUL estimates suggested in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]:

αmM (R) = α • min R + (1 -α) • max R ( 13
)
where R is a set of RUL estimates and αmM (R) the combination result. For example, in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], α = 13/23. In this paper, we considered α ∈ {0.1, 0.2, 0.3, . . . 0.9, 13/23}. The authors in (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]) also added two outlier re-moval (OR) rules to keep RULs within the interquartile range:

OR : {a ∈ R : a ∈ [q 25 , q 75 ]} (14)
and W L : {a ∈ R : q 50 -3•(q 50 -q 25 ) < a < q 50 +2•(q 75 -q 50 )} (15) The set of RUL estimates provided by the algorithm, considering either discrete (K) or continuous predictions (P), is denoted R ≡ RUL

[OR|W L],[th],M P[K] (16)
Only the M first RULs estimates were taken into account (sorted according to the F 1 measure) with M ∈ {11, 15} in this study. OR|W L means that one of the outlier removal operators was applied. The optional parameter [th] means that only training instances with F 1 measure greater than 0.5 were kept.

Weighted average is the second family of rules:

mw [e],[OR] L = L i=1 ω [e],[OR] i • R (i) ( 17
)
where the weights are made dependent on the similarity F 1,i (Eq. 8) between the testing instance and the ith training instance; R (i) is the ith RUL estimate in set of RULs R sorted in descending order with respect to the similarity (F 1,i ) ; L ∈ {3, 5, 7, 9, 11, 15} is the number RULs kept to compute the average while applying or not the outlier removal rule OR.

The weights are given by the following equations:

ω i = F 1,i / L k=1 F 1,k , (18)
with softmax normalisation:

ω e i = exp(F 1,i )/ L k=1 exp(F 1,k ) , (19)
using outlier removal (OR):

ω OR i = OR(F 1,i )/ L k=1 OR(F 1,k ) , (20)
and combining OR and softmax:

ω e,OR i = exp(OR(F 1,i ))/ L k=1 exp(OR(F 1,k )) . ( 21
)
The third kind of rules is a combination of the previous ones:

RUL = 0.5 • αmM (R) + 0.5 • mw [e],[OR|W L] L (22)
Considering several combinations of parameters, about 3168 rules were considered.

Selecting the subset of sensors

As shown by the literature review presented beforehand, many combinations of features can be used (among 21 variables), and a subset was particularly used made of features {7, 8, 12, 16, 17, 20} (involving key sensors for the turbofan degradation [START_REF] Sarkar | Data-driven fault detection in aircraft engines with noisy sensor measurements[END_REF]). To this preselection, a subset of sensors was added from every possible subsets with cardinality equal to 1, 2, 3 and 4 in {∅, 9, 10, 11, 13, 14, 18, 19, 22, 25, 26} as well as subsets of cardinality 5 comprising sensor 9 leading to a total of 511 cases. For each combination (511 cases for each dataset), we applied the prognostics algorithm RULCLIPPER introduced previously and the best subset was selected by minimising the scoring function.

Testing datasets

Given the training instances of a given dataset, the first task is to create a testing dataset in order to estimate 1) the input features and 2) the fusion RUL of RULCLIPPER. The training instances were truncated at a time instant randomly selected from a uniform distribution between 10% and 80% of the half-remaining life. This procedure allowed to obtain instances with small enough RULs to allow the occurrence of substantial degradation, and also large enough RULs to test the robustness of algorithms [START_REF] Hu | Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life[END_REF]. The obtained testing datasets were used in RULCLIPPER with all subsets of features (511 subsets, 3168 fusion rules) and with two subsets of features (511×511 combinations for each fusion rule).

RESULTS AND DISCUSSIONS

Results are presented and compared to past work for all datasets (turbofan and data challenge). More details are given for datasets #1, #3 and #5 T and #5 V .

Performances on datasets #1 and #3

The results can be represented in the penalty -accuracy plane for all combinations of features. A point in that plane with coordinates (P 1 (S 1 , A 1 )) is obtained by considering the accuracy (A 1 ) for the lowest penalty score (S 1 ) given a subset of features. In order to represent the imprecision concerning the performances, a second point P 2 (S 2 , A 2 ) is taken and defined by the lowest score plus 25% (S 2 = S 1 + 25%S 1 ) with accuracy (A 2 ): these two points define a rectangle in the penalty -accuracy plane.

Figure 5 represents the accumulation of these rectangles for all combinations of features in the testing datasets #1 and #3. The whitest part corresponds to the area where most of rectangles are located and thus to the likeliest performances given several subsets of features. If the white area is large then it means that the subset of features should be carefully selected. If the area is concentrated then several subsets of features provided similar performances: it is an image of the robustness with respect to the choice of the subsets. The scores have been divided by the number of testing instance for comparison purpose.

For dataset #1, the performance's centroid is located around (60%; 4.0) (or (60%; 400)). One can draw any subset of features (among the 511 combinations considered) and can expect a score between S = 310 and S = 440 with an accuracy between A = 56 and A = 64. A few "optimal" subsets led to better performances (reported in Table 4 for the testing datasets). The effect of the fault mode on the performances is important. The scores are more spread and a clear global decrease of the accuracy is observed (translation of the cluster of performances to the left hand-side). The level of the colorbar indicates that the choice of the features becomes more and more crucial as the difficulty of the dataset increases: it is simpler to find a subset of features for dataset #1 leading to low penalty and high accuracy because the level is quite similar on a large area (with a value around 8). However, it is more sensitive for dataset #3 for which the level is around 12 on a very local area. A similar (and magnified) observation was obtained on the other datasets. Based on these results obtained on the testing datasets, the fusion rule and the subset of features were selected for final evaluation of the testing datasets by minimising the scoring function (as done for the PHM data challenge) and maximising the accuracy. The results obtained on the testing datasets #1 and #3 are summarised in Table 4 (note that the features indicated in the table have to be assembled with features 7, 8, 12, 16, 17, and 20). For each dataset, the combinations of features are given with respect to the two best scores ("Best S") and the two best accuracies ("Best A"). For example, the first line of Table 4 concerns dataset #1 for which the best score is S = 261 (with A = 63%) when using features 9, 10, 14, 25 and 26, and the RUL fusion "0.9mM ( RUL th,11 P

) ⊕mw 7 " which corresponds to the combination of two elements: 1) the output of the min/max operator (Eq. 13) with parameter α = 0.9 applied on the 11 first RUL estimates and keeping only estimates with a similarity greater than 0.5, and 2) the weighted average (Eq. 17) of the L = 7 first RUL estimates after outlier removal. The high value of α (0.9) implies more weight to the minimum (early) estimate. An accuracy of 70% on #1 was obtained with the same subset of features while keeping a low score at S = 301. This accuracy obtained by the RULCLIPPER algorithm is significantly higher (+16%) than the previous known results given by the EVIPRO-KNN algorithm (Ramasso et al., 2013) which yielded 53%. Other metrics were computed (see Table 6) for performance comparison with previous approaches: An exponential-based regression model with health index estimation proposed in [START_REF] Liu | A data-level fu-sion model for developing composite health indices for degradation modeling and prognostic analysis[END_REF] that provided MAPE = 9% on #1 and an Echo State Network with Kalman filter and submodels of instances presented in [START_REF] Peng | Integrated fatigue damage diagnosis and prognosis under uncertainties[END_REF] with1 MSE = 3969.

The part entitled (#1, #3)/S indicates the best scores for the same subset of features tested with the same fusion method on both datasets. Considering simultaneously #1 and #3 is equivalent to a situation where the engine is degrading while developing a fault. As the score is low and the accuracy high on both datasets using the same subset of features and the same method, it means that this parameterisation makes the prognostics robust to the introduction of the fault mode.

RULCLIPPERs ensemble to manage sensors faults

Two RULCLIPPERs were considered, each with one particular subset of features. All couples of subsets of features were studied (about 130000 combinations) on each testing dataset. The best couples are given in Table 5.

Beyond the important improvement of scores and accuracies compared to the previous results (Table 4), it shows it is not enough to take the two subsets leading to the two best results and expecting an improvement of the performances. Indeed, in most cases, performances for the single feature subsets selected are not the best ones, but their combination yielded to significative improvement of the performances compared to Table 4. For example, for dataset #1, combining RUL estimates provided by subset of features {10, 11, 14, 22} (in addition to 7, 8, 12, 16, 17, 20) (in addition to 7, 8, 12, 16, 17, 20) leading to the best performances in terms of scores (and the corresponding accuracies) for each dataset using the fusion of two RULCLIPPERs. formances obtained in Table 4 with subset {13, 14, 18, 25}, and more when considering the performances of single subsets (S 1 = 301 and P 1 = 64%, or S 2 = 325 with P 2 = 62%). Similar observations can be made on #3.

4.3.

Results on the PHM data challenge (#5 T , #5 V )

Based on the 218 training instances provided, RULCLIPPER was run on the testing dataset using the 511 combinations of features with the 3168 fusion rules. The results were sorted by ascending order with respect to the scoring function. The first five best subsets of features were then selected: {9, 11, 26}, {9, 18, 22, 25}, {9}, {9, 10, 13, 25}, {9, 10, 18, 25, 26} (in addition to 7, 8, 12, 16, 17, 20 for each subset). These combinations of features were considered and evaluated on the dataset #5 T . The best score was given by averaging three configurations of RULCLIPPERs, each with ensembles based on three subsets of features:

• RULCLIPPER 1 with inputs {9, 11, 26}, {9, 18, 22, 25} and {9};

• RULCLIPPER 2 with inputs {9, 11, 26}, {9, 18, 22, 25} and {9, 10, 13, 25};

• RULCLIPPER 3 with inputs {9, 11, 26}, {9, 18, 22, 25} and {9, 10, 18, 25, 26}.

The RUL limit was set to 135 as described in Section 3.4 and the fusion rule was the same for all individual RULCLIP-PER, namely 0.9mM ( RUL 11 PK )⊕mw OR 15 . The score obtained on dataset #5 T (on the NASA's website) was equal to 752, which is the 3rd score compared to published works. An alternative was considered by increasing the RUL limit from 135 to 175. The fusion rule was the same as previously and the score obtained was 934 which is quite low relatively to the high risk taken by setting the RUL limit to 175.

The average of the three configurations given above provided a set of RULs parameterised by both a RUL limit (135, 175) and a fusion method. Three parameterisations were considered and combined:

Ω 1 = (135, 0.8mM ( RUL 11 PK )⊕mw OR 5 ), Ω 2 = (175, 0.9mM ( RUL 11 PK ) ⊕ mw OR 9
), and Ω 3 = (175, 0.9mM ( RUL 11 PK ) ⊕ mw 15 ). The motivation of this configuration was to make long-term predictions possible while minimising the risk of making late predictions. The RULs obtained by Ω 2 and Ω 3 were averaged and the resulting combined by a weighted average with with Ω 1 . The weights were set by a sigmoid (with shape parameter: 0.3 and position: 120) to increase the importance of RULCLIPPERs Ω 2 and Ω 3 when the estimation is greater than 120 while giving more importance to Ω 1 otherwise. This methodology was then applied with the final testing dataset (#5 V ) yielding 11672. The comparison with approaches can be quantified on Figure 6. The generalisation of RULCLIP-PER parameterised as proposed in this paper is lower than the first five algorithms (see square markers on the left-hand side). Indeed, some of these algorithms provided higher scores on #5 T but lower on the final dataset #5 V . One explanation accounting for the lack of generalisation capability compared to the first five algorithms may hold in the "rules" integrated in RULCLIPPER (section 3.4). These rules have been tuned according to observations on the five other datasets but may be not relevant for dataset #5 V if the statistics governing the generation of instances have been modified (Saxena, Goebel, et al., 2008). In order to show the applicability of RULCLIP-PER algorithm with as less parameterisation as possible, the author intentionally kept the same settings for all datasets without distinction in particular concerning the number of fault modes or thresholds on RUL limits. However, the generalisation is better than the 23 remaining algorithms, for which lower scores on #5 T have been obtained with higher ones on #5 V (see square markers on the right-hand side). RULCLIPPER provided a relatively low score on both datasets using the same parameters (816 on #5 T and 11672 on #5 V ). The authors remarked on the previous datasets (#1 to #4) that a few instances can disturb the algorithm (in particular to test the robustness), generating very late or very early predictions, degrading drastically the scores. The important difference on scores between #5 T and #5 V can be due to this particularity.

A summary of results of RULCLIPPER on C-MAPSS datasets is given in Table 6. The best performances were selected according to the scoring function (better accuracies can be obtained as shown in previous tables but with lower scores). Metrics are defined in [START_REF] Saxena | Metrics for evaluating performance of prognostic techniques[END_REF].

CONCLUSION

The RULCLIPPER algorithm is proposed to estimate the remaining lifetime of systems in which noisy health indices are assumed imprecise. RULCLIPPER is made of elements inspired from both the computer vision and computational geometry communities and relies on the adaptation of casebased reasoning to manage the imprecise training and testing instances. The combination of these elements makes it an original and efficient approach for RUL estimation as shown in experiments.

RULCLIPPER was validated with the six datasets coming from the turbofan engine simulator (C-MAPSS), including the so-called turbofan datasets (four datasets) and the data challenge (two datasets), and compared to past work. These datasets are considered as complex due to the presence of fault modes and operating conditions. In addition to RUL-CLIPPER, a method was proposed to estimate the health indi-cator (in presence of faults and operating conditions) and the problem of the selection of the most relevant sensors was also tackled. Information fusion rules were finally studied to combine RUL estimates as well as ensemble of RULCLIPPERs. The review of past work, the presentation of the datasets, as well as the results on sensor selection, health index estimation, information fusion rules and RULCLIPPER ensembles are expected to give a hand to other researchers interested in testing their algorithms on these datasets.

The use of the same algorithm (RULCLIPPER) for all datasets lets suppose that, more generally, computational geometry seems promising for PHM in presence of noisy HIs. However, as for all similarity-based matching algorithm (T. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], the computational cost associated to sort instances is the most important limitation of RULCLIPPER. Two solutions can be considered. Firstly, since operations on convex polygons are simpler (and faster), a procedure can be used to approximate the bounds and to decrease the number of segments.

The second solution concerns implementations, particularly of the intersection algorithm. Computational geometry has become a very active field in particular to improve memory and time requirements, with applications in multimedia (computer graphics such as games). CUDA implementations on processor arrays (using graphic cards) can be pointed out as a promising solution. With such implementations, real-time and anytime prognostics can be considered. The extension of RULCLIPPER to multiple health indices is also under study, in particular by using polytopes.

ACKNOWLEDGMENT

This work has been carried out in the Laboratory of Excellence ACTION through the program "Investments for the future" managed by the National Agency for Research (references ANR-11-LABX-01-01). The author also would like to thank A. Saxena and K. Goebel to produce and make available the C-MAPSS datasets. Comput., C28, 643-647. Chazelle, B., & Edelsbrunner, H. (1992). An optimal algorithm for intersecting line segments in the plane. J.

RULCLIPPER performance

Figure 1 .

 1 Figure 1. Effect of fault modes and operating conditions on health indices estimation. HIs (here obtained from training instances) are described with planar figures called polygons.

Figure 2 .

 2 Figure 2. Illustration of recall and precision.

Figure 3 .

 3 Figure 3. Polygon instances: training (P 1 ) and testing (P * ).

Figure 4 .

 4 Figure 4. The sequence of operations involved in the proposed approach.

Figure 5 .

 5 Figure 5. Performances for #1 (top) and #3 (bottom).

  Al-Salah, T.,Zein-Sabatto, S., & Bodruzzaman, M. (2012). Decision fusion software system for turbine engine fault diagnostics. In Southeastcon, 2012 proceedings of ieee (p. 1-6). Beer, M., Ferson, S., & Kreinovich, V. (2013). Imprecise probabilities in engineering analyses. Mechanical Systems and Signal Processing. Bentley, J., & Ottmann, T.(1979). Algorithms for reporting and counting geometric intersections. IEEE Trans.

Figure 6 .

 6 Figure 6. Comparison of RULCLIPPER with other state-of-the-art approaches. Some scores on the testing dataset #5 T are missing. Scores have to be minimised.

  The set of vertices is obtained by first rearranging the data points to define an ordered sequence that is made possible by extracting the upper and lower envelopes of the noisy HI. For that,

let's define H = [x 1 x2 . . . xj . . . xT ] T a smooth HI obtained by applying a filter over H such that the extraction of both envelopes of H is made easier. The filter used in this paper was a 15-point moving average.

Table 1 .

 1 Datasets characteristics according to the organisers. In this paper, results for all datasets are provided in the experiments, but more details are given specifically for datasets #1 and #3. Note that the datasets called "data challenge" have a common training datasets made of 218 instances. The "semi-final" testing dataset (#5 T ) is made of 218 instances and the "final" validation dataset (#5 V ) is made of 435 instances.

	faults	1	1	2	2	1	1
	Nb. of operating conditions	1	6	1	6	6	6
	Nb. training instances	100 260 100 249	218
	Nb. testing instances	100 259 100 248	218	435
	Minimum RUL	7	6	6	6	10	6
	Maximum RUL	145 194 145 195	150	190

  [START_REF] Wang | A similaritybased prognostics approach for remaining useful life estimation of engineered systems[END_REF]

	Algo. (pseudo.) / Data	#5 T	#5 V
	heracles (1)	737 (3rd)	5636 (1st)
	FOH (2)	512 (2nd)	6691 (2nd)
	LP (3)	n.a.	25921
	sunbea	436.8 (1st) 54437 (22nd)
	bobosir	1263	8637
	L6	1051	9530
	GoNavy	1075	10571
	beck1903	1049	14275
	Sentient	809	19148
	A	975	20471
	mjhutk	2430	30861
	RelRes	1966	35863
	phmnrc	2399	35953
	SuperSiegel	1139	154999
	Table 3. Pseudonyms and scores (known on both #5 T and
	#5 V ) during the 2008's PHM data challenge. Methods (1),
	(2) and (3) were published
	T. Wang, 2010; Heimes, 2008; Riad et al., 2010)) be-
	cause too large and late estimates are greatly penalized
	by the scoring function. So, for most of tests presented
	below, the RUL was given by max(6, min( RUL, 135))

Table 4 .

 4 with {13, 18, 19, 22} led to S = 216 and P = 67%. It represents 27% of improvement on the score and +4% on accuracy compared to the best per-Subset of features (in addition to7, 8, 12, 16, 17, 20) leading to the best performances in terms of scores (and the corresponding accuracies) for each dataset using a single RULCLIPPER.

			TYPE	DATA FEATURES	FUSION			S	A
			Best /S Best /A Best /S Best /A Best /A (2)	#1 #1 #3 #3 #3	9, 10, 14, 25, 26 0.9mM ( RUL 9, 10, 14, 25, 26 0.9mM ( RUL 9, 13, 14, 22, 26 0.9mM ( RUL 9, 19, 25 0.8mM ( RUL 18, 25, 26 mw e,OR e,5 ⊕ mw 5 th,11 P th,11 P W L,11 ) ⊕ mw 7 ) ⊕ mw 13 PK ) ⊕ mw OR 3 W L,11 ) PK	272 68 301 70 353 57 632 63 580 63
			Best /A (3)	#3	9, 10, 14, 18, 25 0.8mM ( RUL	W L,11 PK	) ⊕ mw 3	476 60
			(#1, #3) /S	#1 #3	9, 11, 14, 25, 26 0.9mM ( RUL --	th,11 P	) ⊕ mw e,OR 9	294 64 480 55
			(#1, #3) /S(2)	#1 #3	9, 10, 14, 25, 26 0.9mM ( RUL --	th,11 P	) ⊕ mw OR 13	299 63 480 56
			(#1, #3) /S(2)	#1 #3	9, 10, 14, 25, 26 0.9mM ( RUL --	th,11 P	) ⊕ mw OR 5	315 66 435 54
	Data	Type	Subset 1	Features S 1 P 1 Subset 2	S 2 P 2			Fusion	S	P
	#1	Best /S 10, 11, 14, 22 301 64 Best /A 10, 11, 14, 22 301 64	13, 18, 19, 22 13, 18, 19, 22	325 62 325 62	0.9mM ( RUL 0.9mM ( RUL	th,11 P th,11 P	) ⊕ mw OR 5 ) ⊕ mw OR 7	216 67 224 70
	#3	Best /S 13, 19, 25, 26 525 61 9, 13, 14, 22, 26 353 58 Best /A 14, 18, 19, 26 499 61 9, 10, 13, 14, 26 399 58 0.9mM ( RUL 0.9mM ( RUL W L,15 th,11 P P ) ⊕ mw OR ) 11	317 59 332 63

Table 5 .

 5 Combination of subsets of features

Table 6 .

 6 Summary of results.

	Dataset	#1	#2	#3	#4	#5 T	#5 V
	Score	216 2796 317 3132 752 11672
	Accuracy (%) 67	46	59	45	n.a.	n.a.
	FPR (%)	56	51	66	49	n.a.	n.a.
	FNR (%)	44	49	34	51	n.a.	n.a.
	MAPE (%)	20	32	23	34	n.a.	n.a.
	MAE	10	17	12	18	n.a.	n.a.
	MSE	176 524 256 592	n.a.	n.a.

Authors in(Peng et al.,

2012) actually provided the best RMSE obtained equal to 63, so MSE was computed as 3969 = 63 2 .


Keywords: Irreversible port Hamiltonian systems, Passivity based control, Irreversible thermodynamics, Entropy creation, CSTR

Irreversible port Hamiltonian systems are a class of pseudo Hamiltonian systems that expresses both the conservation of energy and the irreversible entropy production as a structural property. These systems encompass a large class of irreversible themordynamic systems, such as heat exchangers and chemical reactors, and also multi-energy systems such as coupled mechanicthermodynamic systems. In recent work the irreversible port-Hamiltonian formulation has been used to derive a closed-loop stability condition using an energy based availability function, generated by the internal energy, as Lyapunov function. This paper presents an important extension of the previous results: the system theoretic interpretation of the stability condition in terms of conjugated inputs and outputs and the formulation of the control as an interconnection and damping assignment -passivity based control problem. A constructive method to derive the stabilizing control law is proposed and the formalism is illustrated on a general CSTR example.

Interconnection and Damping Assignment -Passivity Based Control of Irreversible Port Hamiltonian Systems

Hector Ramirez * Yann Le Gorrec * Bernhard Maschke * *

INTRODUCTION

The attempts to model irreversible processes as PHS have lead to a class of system called quasi PHS [START_REF] Hangos | Hamiltonian view on process systems[END_REF][START_REF] Otero-Muras | Local dissipative Hamiltonian description of reversible reaction networks[END_REF][START_REF] Ramirez | On the control of non-linear processes: An IDA-PBC approach[END_REF][START_REF] Dörfler | An introduction to interconnection and damping assignment passivity-based control in process engineering[END_REF][START_REF] Hoang | The port Hamiltonian approach to modelling and control of continuous stirred tank reactors[END_REF]. These systems retain much of the port Hamiltonian structure, but differ by their structure matrices and input vector fields which depend explicitly on the gradient of the Hamiltonian. An important remark is that, although the forms of PHS and quasi PHS are very similar and both embed, by skew-symmetry of the structure matrix, the conservation of the Hamiltonian, in the latter the drift dynamic is a nonlinear function in the gradient ∂U ∂x (x). In this sense the symplectic structure of the PHS, given by the Poisson tensor associated with the structure matrix J(x), is destroyed. From a control perspective it is usually more complicated to impose a desired closed-loop dynamic on quasi PHS since passivity based techniques can not be easily applied due to the non-linearity with respect to the gradient of the Hamiltonian.

In [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF] a class of quasi PHS, denoted Irreversible Port Hamiltonian Systems (IPHS), has been proposed to model a large class of thermodynamic systems and to embed the first (conservation of the Hamiltonian) and the second principle (irreversible increase of an entropy like function). IPHS have two important properties: they provide an irreversible non-linear system a physically meaningful structure (just as PHS for mechanical systems); and they are defined with respect to the total energy, which makes it possible to interconnect them with reversible non-linear systems (mechanical systems) (Ramirez et al., 2013c), which is not the case for formulations defined with respect to the total entropy for instance.

In Ramirez et al. (2013a) the framework of IPHS has been combined with the framework of thermodynamic availability function to derive an asymptotic stability condition for irreversible thermodynamic systems. Indeed, since the closed-loop equilibrium differs from the openloop equilibrium it is necessary to perform the closed-loop stability analysis with respect to a modified Hamiltonian, which in our case is derived departing from the total energy (open-loop Hamiltonian). In this paper we specialize the previous result for input-affine IPHS (Section 2) and interpret the stability condition in terms of interconnection and damping assignment -passivity based control [START_REF] Ortega | Putting energy back in control[END_REF][START_REF] Ortega | Interconnection and damping assignment passivity based control of port-controlled Hamiltonian systems[END_REF] (Section 3). It is interesting to see that the closed-loop system can again be interpreted as an irreversible thermodynamic systems with respect to a modified thermodynamic potential (Section 3). A constructive method to derive the stabilizing control law is proposed and the formalism is illustrated on a general CSTR example (Section 4). Finally (Section 5), some closing remarks and comments on future work are given.

IRREVERSIBLE PORT-HAMILTONIAN SYSTEMS

Irreversible Port Hamiltonian Systems (IPHS) have been defined in [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF] as a modification of Port Hamiltonian systems (PHS) and allow to represent simultaneously the conservation of energy and the entropy creation due to an irreversible processes. An IPHS uses two thermodynamic potentials in order to define the dynamics: the total energy and the total entropy of a thermodynamic system and is defined as follows. Definition 1. An input affine IPHS is defined by the dynamic equation

ẋ = R x, ∂U ∂x , ∂S ∂x J ∂U ∂x (x) + g x, ∂U ∂x u, y = g x, ∂U ∂x ∂U ∂x (x) (1)
where x ∈ R n is the state vector, U (x) :

C ∞ (R n ) → R
is the Hamiltonian function (for thermodynamic systems, the internal energy) and S(x) : C ∞ (R n ) → R is a real function (for thermodynamic systems the total entropy function), J ∈ R n×n is the constant skew-symmetric structure matrix of the Poisson bracket acting on any two C ∞ (R n ) functions Z and G as:

{Z, G} J = ∂Z ∂x (x)J ∂G ∂x (x). (2)
The drift term of the system (1) differs from a Hamiltonian system by the multiplication by the real function R = R x, ∂U ∂x , ∂S ∂x . This function is composed of the product of a positive definite function γ and the Poisson bracket of the entropy function S (x) and the energy function U (x): R x, ∂U ∂x , ∂S ∂x = γ x, ∂U ∂x {S, U } J , with γ(x, ∂U ∂x ) = γ(x) : C ∞ (R n ) → R, γ ≥ 0, a nonlinear positive function. The input map is defined by g(x, ∂U ∂x ) ∈ R n×m with the input u(t) ∈ R m a time dependent function.

The balance equations of the entropy function S (x) and the energy function U (x) of the IPHS express the first and second principles of the irreversible thermodynamics: the conservation of energy and the irreversible entropy creation due to the irreversible processes. Indeed consider the balance equation of the internal energy, one computes dU dt = y u.

where u(t) and y(z) are the input and output port of the system. This is precisely the balance equation corresponding to lossless dissipative systems with the (energy) supply rate y u. The balance equation of the entropy function is given by

dS dt = R x, ∂U ∂x , ∂S ∂x ∂S ∂x J(x) ∂U ∂x + ∂S ∂x g x, ∂U ∂x u = γ x, ∂U ∂x {S, U } 2 J + g x, ∂U ∂x ∂S ∂x u
From the assumptions defining the IPHS (Definition 1) the first term is positive: γ x, ∂U ∂x {S, U } 2 J = σ x, ∂U ∂x ≥ 0 which represents for models of irreversible thermodynamic systems, the internal entropy production. The second term corresponds to the definition of an entropy supply rate. For further details on IPHS and its thermodynamic interpretation we refer the reader to [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF].

IDA-PBC FOR IPHS

In Ramirez et al. (2013a) the framework of IPHS and available storage [START_REF] Willems | Dissipative dynamical systems part I: General theory[END_REF] has been used to propose a class of PBC for irreversible thermodynamic systems. Using the convexity of the internal energy, the Hamiltonian function of an IPHS has been used as convex extension to define an energy based availability function which serves as Lyapunov function for the closed-loop system. This has been done using the same thermodynamic arguments as for instance in [START_REF] Alonso | Process systems, passivity and the second law of thermodynamics[END_REF]; [START_REF] Ydstie | Process systems and passivity via the clausius-planck inequality[END_REF]; [START_REF] Alonso | Stabilization of distributed systems using irreversible thermodynamics[END_REF]; [START_REF] Ydstie | Passivity based control via the second law[END_REF]; [START_REF] Hoang | The port Hamiltonian approach to modelling and control of continuous stirred tank reactors[END_REF][START_REF] Hoang | Lyapunov-based control of non isothermal continuous stirred tank reactors using irreversible thermodynamics[END_REF], being the difference that in our case the closed-loop Lyapunov function is generated by the energy of the system and not the entropy. This permits to state the control problem similarly as for mechanical systems, with the difference that the irreversible nature of IPHS introduces an "interconnection obstacle" (Ramirez et al., 2013c), also refereed as "dissipation obstacle" in the case of dissipative PHS [START_REF] Ortega | Putting energy back in control[END_REF]. In this section we shall elaborate on the IDA-PBC method for deriving a stabilizing control for which the closed-loop Hamiltonian function is the energy-based availability function, i.e., we will investigate the existence of controls that renders the closed-loop system into the form

ẋ = M(x) ∂A ∂x (x)
where A is the energy-based availability function and M(x) the closed-loop structure matrix, possibly defining an (dissipative) IPHS.

The Energy based availability function

Let us briefly recall some fundamental properties of irreversible homogeneous thermodynamic systems and more precisely the notion of the thermodynamic availability function associated with the internal energy (Ramirez et al., 2013c). The variation of the internal energy of an homogeneous system is defined by Gibbs' fundamental equation on a N -dimensional state space:

dU = T dS -P dV + N -2 i=1 µ i dn i (3)
where the extensive variables are the internal energy U , the entropy S, the volume V and the mole numbers n i , and the intensive variables are the temperature T , the pressure P and the chemical potentials µ i . By defining the vectors of extensive and intensive variables, z = [S, V, n 1 , . . . , n N -2 ] and w

(z) = [T (z), -P (z), µ 1 (z), . . . , µ N -2 (z)],
(3) may be written as dU = w(z) dz.

(4) The internal energy U is a homogeneous function of degree 1, so from Euler's Theorem we obtain, U = w(z) z, (5) which implies that w(z) is a homogeneous function of degree 0 of z and w(z) = ∂U ∂z (z). For homogeneous systems, as a consequence of the second law of Thermodynamics [START_REF] Callen | Thermodynamics and an introduction to thermostatistics[END_REF][START_REF] Alonso | Stabilization of distributed systems using irreversible thermodynamics[END_REF], the internal energy corresponds to a convex function. Furthermore the internal energy function is strictly convex with respect to z and independent of the dynamic behaviour of the system if one of the extensive variables, i.e. one of the elements in z, is fixed [START_REF] Jillson | Process networks with decentralized inventory and flow control[END_REF]. For thermodynamic systems, the total energy function (and hence also the total entropy function) only possesses a strict minimum at the thermodynamic equilibrium, which does not coincide with the desired dynamic equilibrium. Indeed, at thermodynamic equilibrium the system does not undergo any irreversible transformation, hence it is not of practical interest.

Since the Hamiltonian of IPHS is the total energy it is natural to derive a closed-loop Lyapunov function using the convexity of the total energy, rather than the concavity of the total entropy. Recalling Definition 1, and taking x as the state vector of the IPHS of we give the definition of the energy based availability function and the closed-loop stability condition. Definition 2. The energy based availability function is given by

A(x, x * ) = U (x) -U (x * ) + ∂U ∂x (x * ) (x -x * ) ≥ 0 (6)
where x * is a reference and possibly a desired equilibrium. For a thermodynamic system A becomes a strictly convex with minimum value A = 0 if one of the extensive variables is fixed [START_REF] Jillson | Process networks with decentralized inventory and flow control[END_REF]

, i.e., if z ∈ R N -1 = R n
and the state space is chosen as x = z. Furthermore, the only possible value for which

A = 0 is x = x * .
The reader is refereed to [START_REF] Alonso | Process systems, passivity and the second law of thermodynamics[END_REF]; Ydstie and Alonso (1997); [START_REF] Alonso | Stabilization of distributed systems using irreversible thermodynamics[END_REF]; [START_REF] Ydstie | Passivity based control via the second law[END_REF] for further details on the construction and thermodynamic properties of the availability function. It is clear from Definition 2 that the energy based availability function qualifies as a Lyapunov function candidate for controlled IPHS. Proposition 1. The closed-loop equilibrium x * of a controlled IPHS is asymptotically stable if A is strictly convex and the control u (x) satisfies:

γ {S, U } J {A, U } J + ∂U ∂x (x) - ∂U ∂x (x * ) gu ≤ 0, (7)
with strict equality only at (x * , u * ), the steady state value at the desired equilibrium.

Proof. The time derivative of A is given by

dA dt = ∂U ∂x (x) - ∂U ∂x (x * ) dx dt . (8)
Using Definition 1 we obtain the balance equation of the availability function

dA dt = -R ∂U ∂x (x * ) J ∂U ∂x (x) + ∂U ∂x (x) - ∂U ∂x (x * ) gu, = γ {S, U } J {A, U } J + ∂U ∂x (x) - ∂U ∂x (x * ) gu, (9)
By positivity of γ and since the availability A qualifies as a Lyapunov function for the system it is stable if (7) holds. If x * is an isolated minimum, which is the case if A is strictly convex, then asymptotic stability follows invoking La Salle's invariance principle on a region around x * .

The IDA-PBC problem

Proposition 1 provides an asymptotic stability condition for the closed-loop IPHS, but provides no constructive information on how to choose the control input u. Moreover, it is interesting to investigate what (7) represents in the sense of passivity based control. Using the port-conjugated output, we may rewrite (7) as

γ {S, U } J {A, U } J + y u - ∂U ∂x (x * ) gu ≤ 0. ( 10
)
Let us define the availability-conjugated output defining an availability supply rate in the balance equation ( 9) of the availability function

ỹ = g ∂U ∂x (x) -∂U ∂x (x * ) = g ∂A ∂x (x), (11)
which actually only becomes zero at the desire equilibrium.

Then, the stability condition (10) becomes

γ {S, U } J {A, U } J + ỹ u ≤ 0. ( 12
)
Now, define the control input as

u = γ {S, U } J ũ, with ũ a new control input, to obtain γ {S, U } J {A, U } J + γ {S, U } J ỹ ũ ≤ 0 γ {S, U } J {A, U } J + ỹ ũ ≤ 0 (13)
Notice that this condition may be compared to the stabilization condition for reversible Hamiltonian systems. In this case one considers 1 instead of R = γ {S, U } J and the condition {A, U } J + ỹ ũ ≤ 0 expresses the decrease of the closed-loop Hamiltonian function, chosen as candidate Lyapunov function. In that case, it suffices to select the control input as the static output feedback ũ = -K ỹ, with K = K > 0, to stabilize the closed-loop system, since {U, U } J (x) = ỹ(x * ) = 0. For Irreversible Port Hamiltonian systems, from (13) we see however that this condition is modulated by the nonlinear term R = γ {S, U } J , which is not signed defined since the thermodynamic driving force {S, U } J may take positive or negative values. A particular solution to ( 13) is

ỹ ũ = -{A, U } J -{S, U } J [A, A] M + {A, A} J (14)
where J (x) is any skew symmetric matrix defining a Poisson bracket, and M (x) is a symmetric and positive semi-definite matrix defining a pseudo-Riemannian metric expressed by the Ginzburg-Landau dissipative bracket, defined in some local coordinates for two C ∞ (R n ) functions Z and G as,

[Z, G] M = ∂Z ∂x (x)M (x) ∂G ∂x (x). ( 15
)
Notice that {A, A} J = 0 and does not contribute in (13). This bracket is included however since it shows that the closed-loop vector field may be shaped without injecting more damping into the system. This is referred to as the interconnection assignment. Now, using the definition of ỹ (Eq. ( 11)) we have that the unique solution to (14) is 16) which exists if and only if

β(x) = g † (x) -{S, U } J M (x) + J (x) ∂U ∂x (x) - ∂U ∂x (x * ) -g † (x)J ∂U ∂x (x), (
-g ⊥ (x){S, U } J M (x) ∂U ∂x (x) - ∂U ∂x (x * ) + g ⊥ (x) J -J ∂U ∂x (x) -J ∂U ∂x (x * ) (17) with g † (x) = [g (x)g(x)] -1 g (x)
the Moore-Penrose pseudo inverse and g ⊥ (x) a left full rank annihilator of g(x), i.e., g ⊥ (x)g(x) = 0. Equations (16 and ( 17) express the well known IDA-PBC problem, but for irreversible control system.

Let us compare more in details this solution with the solution for a reversible port-Hamiltonian system ẋ = J ∂U ∂x + gũ (18)

keeping the same target system corresponding to the equivalent IDA-PBC problem defined by equations ( 16) and ( 17). The target system is then given by ẋ

= -{S, U } J M (x) + J (x) ∂A ∂x , (19)
where the matrix {S, U } J M (x) is not necessary positive semi-definite, since {S, U } J is of arbitrary sign. This "contradiction" comes of course from the definition of the control input, u = γ {S, U } J ũ, but it also expresses the irreversible nature of the control system. The system is already dissipative, hence damping injecting is not the issue. The issue is to make the closed-loop system dissipate in such a way that it is stable. This contradiction has been reported in the literature as the "dissipation obstacle" for reversible PHS [START_REF] Ortega | Putting energy back in control[END_REF], and more recently, as the "interconnection obstacle" for IPHS (Ramirez et al., 2013c). The previous results can be summarized as follows.

Proposition 2. Consider the input affine IPHS (1). Assume A is strictly convex and that there exist matrices M (x) ≥ 0 and J (x) = -J (x), and a full-rank left annihilator g ⊥ (x) of g(x) that verify (17). Then the control u = γ {S, U } J β(x), with β(x) defined as in ( 16), asymptotically stabilizes the closed-loop equilibrium x * . Furthermore, the closed-loop system is given by the pseudo-PHS

ẋ = -σM (x) + RJ (x) ∂A ∂x . ( 20
)
Proof. The proof follows directly from the previous developments and ( 20) is obtained by replacing the control law

u = γ {S, U } J β(x) in (1).
In some cases it is possible to use the physical properties of a system to solve (17). For instance, if we don't want to change the open-loop structure/interconnection matrix, and only express the natural dissipation of the system with respect to A, then the simplest choice is some constant M > 0 and J = J. Corollary 3. Assume J = J, then the closed-loop equilibrium x * is asymptotically stable if

g ⊥ (x)J ∂U ∂x (x * ) = 0, and (21)
g ⊥ (x)M (x) = 0. ( 22
)
The second condition of Corollary 3 is in general easy to fulfil since the only restriction on M is that it is positive semi-definite and symmetric. The first condition is however very restrictive, since J, U (x) and g(x) (and hence also g ⊥ (x)) are defined by the physical system. Nevertheless, since J expresses the conserved quantities of the thermodynamic process we may find this condition automatically fulfilled for systems whose inputs are related with the physical invariants. This point will be illustrated on the example of the CSTR.

EXAMPLE: THE CSTR

Assume a chemical reaction in a CSTR with the following reversible reaction scheme

ν 1 B 1 + . . . + ν l B l ν l+1 B l+1 + . . . + ν m B m , m > l ≥ 1.
where where B i denotes the ith chemical species, νi is the signed stoichiometric coefficient: νi = -ν i if it appears on the left hand side of the reaction scheme, νi = ν i in the other case. We shall consider that the CSTR satisfies the following standard operation assumptions [START_REF] Aris | Elementary chemical reactor analysis[END_REF][START_REF] Favache | Thermodynamics and chemical systems stability: The CSTR case study revisited[END_REF] Assumption 4.

• The volume, denoted by V , in the reactor is constant as well as the pressure, • The initial number of moles of a specie in the reactor is equal to the number of moles of the inlet of the same specie, • For a given steady state temperature and steady state input there is only one possible steady state for the mass (numbers of moles) balance. • The reactor operates at gas phase,

• The controlled inputs are the the dilution rate u 1 = F/V and the heat flux from the cooling jacket u 2 = Q.

We have assumed a reaction in gas phase, but the developments may be applied identically to a reactor with a reaction in liquid phase. The assumptions of constant volume and pressure impose a constraint over the total outlet flow [START_REF] Couenne | Bond graph modelling for chemical reactors[END_REF][START_REF] Couenne | Structured modeling for processes: A thermodynamical network theory[END_REF], making the outlet flows state dependent. Furthermore this assumption guarantee that the energy based availability function A is strictly convex.

The mass balance equations

The time variation of the species in the reactor are given by [START_REF] Aris | Elementary chemical reactor analysis[END_REF])

ṅi = F V (n ei -n i ) + r i V i = 1, . . . , m (23)
where n i and n ei are, respectively, the number of moles of the species i in the reactor and at the inlet (and n = (n 1 , . . . , n m ) and n ei = (n e1 , . . . , n em ) ), r i = νi r where r (n, T ) is the reaction rate which is the difference of the forward reaction rate r f and the backward reaction rate r b : r = (r f -r b ) and depends on the temperature and on the reactant mole number and F is the volumetric flow rate.

The energy and entropy balance equations

Under Assumption 4 the internal energy of the CSTR is given by [START_REF] Couenne | Bond graph modelling for chemical reactors[END_REF][START_REF] Favache | Thermodynamics and chemical systems stability: The CSTR case study revisited[END_REF])

U = m i=1 n i [c pi (T -T 0 ) + u 0i ], (24)
where c pi , u 0i , T 0 are respectively the heat capacity at constant pressure, reference molar energy and reference temperature. The entropy balance equation may be deduced from Gibbs' relation and is given by

Ṡ = F V m i=1 (n ei s ei -n i s i ) + Q T e + σ, (25)
where s ei and s i are respectively the inlet molar entropy and the molar entropy of species i, σ is the irreversible entropy creation due to mass transfer, heat transfer and chemical reactions:

σ = F V m i=1 n ei T (h ei -T s ei -µ i ) + Q T - Q T e - m i=1 µ i ν i r T ,
where Q = λ(T e -T ) is the heat flux from the jacket with λ the heat conduction coefficient, T e the temperature of the jacket, h ei the inlet specific molar enthalpy.

IPHS model [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF] The derivation of the IPHS model of the CSTR has been presented in details in [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF]. The dynamical equation of the CSTR may be expressed as the IPHS

ẋ = RJ ∂U ∂x (x) + g x, ∂U ∂x u (26)
= [u 1 , u 2 ] ,
where u 1 is the dilution rate F V and u 2 the heat flux from the cooling jacket Q. The first row of the input map is the an input and output mole number vector and second row a thermal interaction vector,

g = ñ 0 φ(x) 1 T with ñ = n e -n, and φ(x) = m i=1 (n ei s ei -n i s i )+ nei T (h ei - T s ei -µ i ).

IDA-PBC of the CSTR

In this subsection we present how to use the IPHS model of the CSTR (open-loop model) to derive a stabilizing control law which imposes a desired IPHS structure on the closed-loop system. We will use Corollary 3 to show that the state feedback 16 asymptotically stabilizes the CSTR at a desired equilibrium point x * . For this purpose we select J = J and should find a matrix M such that ( 22) is satisfied. The first step is to find a full rank left annihilator for the input map. One possible choice is the (m -1) × n matrix

g ⊥ (x) =     ñ2 -ñ 1 0 . . . 0 0 0 0 ñ3 -ñ 2 . . . 0 0 0 0 0 . . . . . . . . . 0 0 0 0 . . . 0 ñm -ñ m-1 0     (27)
The most restrictive condition in Corollary 3 is ( 21), from which we obtain, when using the previously chosen annihilator,

g ⊥ (x)J =     0 . . . 0 ν1 ñ2 -ν2 ñ1 0 . . . 0 ν2 ñ3 -ν3 ñ2 . . . . . . . . . . . . 0 . . . 0 νm-1 ñm -νm ñm-1     .
It suffices to verify that g ⊥ (x)J = 0, which is true if

ñ1 ν1 = ñ2 ν2 = • • • = ñm-1 νm-1 = ñm νm . (28)
From Assumption 4 we have that the initial numbers of moles of each specie equals the numbers of moles at the inlet, i.e., n(t = 0) = n 0 = n e . Hence, ( 28) is actually the expression of De Donder's extent of reaction [START_REF] Callen | Thermodynamics and an introduction to thermostatistics[END_REF])

n 0i -n i νi = ξ, (29)
which is always true and hence (21) also. It is interesting that ( 21) is automatically fulfilled for this example. This is due to that J expresses the stoichiometry of the reaction and g(x) the mole (mass) balance. Since the reactor operates at constant volume, the total mass becomes an invariant for the reaction, and g ⊥ (x)J = 0 is simply the mathematical expression of this invariant.

The second condition, given by ( 22), can be solved by noting that the last column of g ⊥ (x) is zero. This implies that any matrix M (x) = M (x) > 0 for which the first m rows and columns form a null submatrix is solution to (21). A (simple) possible choice is M = diag(0, . . . , 0, 1), which corresponds to a matrix with all elements equal to zero except the last element of its diagonal which is 1. The closed-loop system takes then the form

ẋ = -σM + RJ ∂A ∂x (x, x * )
and the time derivative of A along the trajectories of the CSTR is

dA dt = -σ(T -T * ) 2 .
Hence, the "dissipation" of the energy based availability function is related to the irreversible energy production, which is positive definite and only zero at the thermodynamic equilibrium (which is not included in the region of operation of the CSTR). In order to complete the stability proof, we need to apply La Salle's invariance Theorem in a sufficiently small region around T = T * . From Assumption 4 we have that there is only one equilibrium for each temperature, hence we can conclude that the closed-loop systems is asymptotically stable.

CONCLUSION

The stabilization of irreversible thermodynamic systems has been studied using the framework of Irreversible port-Hamiltonian Systems (IPHS) and the thermodynamic availability function. The availability function is derived from the intrinsic properties of homogeneous thermodynamic systems, and under certain operation conditions can be used to analyse the dynamic properties, such as the stability, of a thermodynamic system, i.e., it may be used as a Lyapunov function for thermodynamic systems. In a first step, an energy based availability function has been defined and used to derive a general asymptotic stability condition for IPHS. It is indeed natural to use the energy, rather than the entropy, to derive the closed-loop availability function since the Hamiltonian of IPHS corresponds to the total energy. Then, in a second step, the stability condition has been developed an interpreted in terms of the matching of vector fields, leading to an Interconnection and Damping Assignment -Passivity Based Control (IDA-PBC) problem. The matching condition and the definition of the control law associated to the IDA-PB controller clearly expresses the irreversible nature of the open-loop system, and it is shown that the closed-loop target system is an dissipative IPHS with the energy based availability as Hamiltonian function. An IDA-PBC synthesis method for IPHS has been proposed and the control design approach has been illustrated on a general model of a chemical CSTR.

  with state vector x = [n , S] , the internal energy U (x) as Hamiltonian function, -symmetric matrix whose elements are the stoichiometric coefficient of the chemical reaction mapping the network structure of the reaction, andR = γ x, ∂U ∂x {S, U } J = rV T A Awith γ = rV T A and {S, U } J = A, where A = -m i=1 νi µ i is the chemical affinity of the reaction and corresponds to the thermodynamic driving force of the chemical reaction. The input port of the IPHS is given by gu, with u

This work was supported by French and French-Chilean sponsored projects HAMECMOPSYS ANR-11-BS03-0002, Labex ACTION ANR-11-LABX-01-01 and ECOS-CONICYT C12E08.