Guillaume Chapuis

Rayan Chikhi

Dominique Lavenier

Parallel and memory-ecient reads indexing for genome assembly

As genomes, transcriptomes and meta-genomes are being sequenced at a faster pace than ever, there is a pressing need for ecient genome assembly methods. Two practical issues in assembly are heavy memory usage and long execution time during the read indexing phase. In this article, a parallel and memory-ecient method is proposed for reads indexing prior to assembly. Specically, a hash-based structure that stores a reduced amount of read information is designed. Erroneous entries are ltered on the y during index construction. A prototype implementation has been designed and applied to actual Illumina short reads. Benchmark evaluation shows that this indexing method requires signicantly less memory than those from popular assemblers.

Introduction

Until the emergence of next-generation sequencing (NGS) technologies, software for assembling genomes could process up to millions of long (∼ 10 4 bp) reads. Now, a typical genome assembly instance for a vertebrate genome consists of billions of short (100 bp) reads. Despite this technological shift, computational models for assembly are essentially based on constructing and simplifying a genome graph. However, graph-based models have inherent limitations that make them unpractical for assembly of NGS data. They require the construction of either a large string graph containing all the reads, or a de Bruijn graph containing all the k-length substrings (k-mers) of the reads. For human-sized genomes, the de Bruijn graph typically requires hundreds of gigabases of memory [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. Nevertheless, NGS assembly tools rely on optimized implementations of these graph models. For instance, leading assembly programs have implemented ecient heuristics using a de Bruijn graph [START_REF] Zerbino | Velvet: Algorithms for de novo short read assembly using de bruijn graphs[END_REF][START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. For more details concerning these implementations, refer to a recent survey [START_REF] Miller | Assembly algorithms for next-generation sequencing data[END_REF]. In a near future, larger eukaryotic genomes and meta-genomes will be sequenced at a faster pace than computational resources growth. Hence, new assembly models need to be developed to sustain the increasing rate of NGS technologies.

Several theoretical advances have been recently proposed to reduce the memory usage of graph-based assemblers. Simpson et al. implemented compression techniques (FM-index [START_REF] Ferragina | Indexing compressed text[END_REF]) during construction of the string graph [START_REF] Simpson | Ecient construction of an assembly string graph using the FM-index[END_REF] at the expense of running time. Conway et al. used succinct bitmap structures [START_REF] Okanohara | Practical entropy-compressed rank/select dictionary[END_REF][START_REF] Raman | Succinct indexable dictionaries with applications to encoding k-ary trees, prex sums and multisets[END_REF] to construct an immutable de Bruijn graph [START_REF] Conway | Succinct data structures for assembling large genomes[END_REF]. Distributed de Bruijn graph construction using a message passing interface have been implemented in several assemblers [START_REF] Jackson | Parallel short sequence assembly of transcriptomes[END_REF][START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF][START_REF] Kundeti | Ecient parallel and out of core algorithms for constructing large bi-directed de bruijn graphs[END_REF].

Greedy assemblers use a dierent assembly strategy. Instead of constructing a genome graph, they repeatedly perform an extension procedure until branching is detected. Previous implementations of greedy assemblers used a prex tree to store reads [START_REF] Warren | Assembling millions of short DNA sequences using SSAKE[END_REF], which consumes signicantly more memory than a de Bruijn graph. Recent optimized implementations use custom k-mer indexing structures for memory eciency [START_REF] Ariyaratne | PE-Assembler: de novo assembler using short pairedend reads[END_REF][START_REF] Boisvert | Ray: Simultaneous assembly of reads from a mix of High-Throughput sequencing technologies[END_REF][START_REF] Chapman | Meraculous: De novo genome assembly with short Paired-End reads[END_REF]. In particular, these implementation have been applied to complex mammalian genomes, demonstrating that greedy assemblers are not limited to genomes with low repeat content. Unlike de Bruijn graph assemblers, data structures used in greedy assemblers typically contain references to read sequences. Hence, ecient read indexing is necessary to keep memory usage low.

In the next section, we propose a parallel reads indexing procedure designed specically for assembly. Two novel ltering methods are introduced to reduce memory usage: a procedure to remove erroneous k-mers on the y, and a procedure to avoid referencing redundant reads. Finally, a prototype implementation is applied to real Illumina data to validate the method.

Methods

Distributed and multi-threaded indexing

A multi-threaded, multi-node procedure for reads indexing is proposed. A hash table is constructed, where the entries are k-mers, and the values are references to reads. Taking advantage of shared memory between threads, reads sequences are stored separately in memory, without redundancy within a node. Index construction is distributed among N nodes, and each node performs independent computations in parallel. Specically, each node n is running T n threads, each thread t n constructs a separate sub-index I(n, t n ). A binning method adapted from [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF] assigns each k-mer to a unique sub-index. Let h be a k-mer hash value with perfect hashing [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]. The corresponding k-mer belongs to the sub-index I(n, t n ) if:

h mod N = n h mod T n = t n
which ensures that each sub-index contains distinct k-mers. Each thread reads the entirety of the input data to construct its sub-index. When all the subindexes are constructed, an inexpensive merging phase yields the complete index. Hence, the indexing procedure always constructs the same complete index on dierent architectures. In the following, two algorithmic ingredients are described for parallel sub-index construction: k-mers ltering and reads indexing.

On-line parallel k-mers ltering

Memory eciency is crucial when assembling NGS data. In many approaches, including the one proposed here, memory consumption is proportional to the number of indexed k-mers. It is therefore important to lter out erroneous kmers as early in the indexing process as possible. Erroneous k-mers are produced whenever the sequencing process makes a mistake during base calling. The abundance distribution K n t (m) is dened as the number of k-mers seen exactly m times at indexing time t by node n. A key fact is that the hash function used above evenly distributes k-mers among sub-indexes. Hence, each

K n t (m)
is identically distributed as the entire distribution n K n t (m). This observation enables independent, parallel ltering for each sub-index. The superscript n is then omitted in the following.

A typical distribution of K t (m) at nal time t is multimodal. A large number of k-mers occur only a few times: these are mostly sequencing errors. Assuming uniform sequencing coverage, the distribution of correct k-mers is a Gaussian mixture. The most abundant component is centered at the expected coverage of the target genome. Less abundant components are centered at multiples of the coverage, due to repeats in the genome. The proposed method consists in (i) detecting components corresponding to erroneous and correct k-mers as soon as they separate suciently from each other and (ii) nding an appropriate erroneous threshold (cut-o value). Every k-mer that has appeared fewer times than the erroneous threshold so far is then considered as an error and removed. This procedure could be extended to correct errors in reads, but it is outside the scope of the current indexing scheme.

Error detection The following two inequalities must be satised to trigger the ltering procedure. First, erroneous k-mers are identied by their abundance. Theorem 3 from [START_REF] Peng | IDBA -a practical iterative de bruijn graph de novo assembler[END_REF] establishes that, under reasonable sequencing assumptions, an error is signicantly less likely to appear m + 1 times than m times. Thus, the abundance of erroneous k-mers peaks at m = 1 and has a strictly decreasing slope. The low end m low (t) is computed as the largest m that satises K t (m -1) > K t (m) for m ≥ 1. Then, the peak abundance m high (t) of correct k-mers is computed as the parameter at which the maximum value of K t (m) is attained for m > m low (t). Erroneous and correct k-mers are considered to be separated when:

m high (t) -m low (t) > r
where r is a user-dened resolution parameter. Second, to avoid the computational cost of ltering too soon or too often, a constraint is imposed on the amount of erroneous k-mers. Let S min be a minimum amount (user-dened) of erroneous k-mers before the ltering process can be performed:

m low m=1 K t (m) > S min
Calculating the cut-o value During early ltering passes, a small fraction of correct k-mers still contributes to the erroneous component. Hence, removing the entire component at each ltering pass is not a sensible choice. An incrementing value, dened as the cut-o value, is introduced to overcome this problem. All k-mers of abundance lower than the cut-o value are removed by the ltering procedure, others are kept. Formally, a threshold m solid is dened as the number of occurrences below which a k-mer is considered a potential error. All k-mers over this threshold at the end of the indexing phase are solid k-mers. Let tReads and nReads(t) be the total number of reads in the input le and the number of reads processed at time t respectively. The cut-o value F (t) is calculated according to the following formula:

F (t) = m solid • nReads(t) tReads

Reads indexing structure

Each sub-index is populated independently with a ltered set of references to reads, given a ltering function designed for de novo assembly. The extension of a k-mer in a read is dened as the sux immediately following the k-mer (e.g. for a read r = uwv where w is a k-mer and u, v are arbitrary strings, v is the extension of w in r). We introduce a notion of redundancy between extensions. Let (v 1 , v 2 ) be two extensions of the same k-mer, without loss of generality assume that the length |v 1 | is shorter than |v 2 |. Two extensions v 1, v 2 are said to be t-redundant if the Hamming distance between their prexes of length |v 1 | is lower than t. The representative read spectrum with similarity threshold t, noted RRS(k, t), is dened for a set of input reads as follows:

(i) associate a set S w to each solid k-mer w occurring in the reads (ii) S w discards all but one of the reads associated to t-redundant extensions. A read with the longest extension is kept, ties are broken arbitrarily. Figure 1 shows an example of a representative reads spectrum. The reads sequences referenced by the RRS are stored separately. Practically, both a read and its reverse-complement are indexed. References to paired-end reads are explicitly made to the left or the right mate of the read. In essence, this structure records a representative set of reads for each solid k-mer. Note that this indexing does not correct errors in read, but merely ignores errors in reads suxes. Erroneous prexes yield un-solid k-mers, hence these reads are not indexed in the structure. This property is well suited with Illumina reads as sequencing errors are known to mostly occur at read suxes. Provided the sequencing coverage is high, errors in suxes can be corrected at a later stage during a consensus phase. This justies the arbitrary removal of other reads having equally long t-redundant extensions. To maximize the eectiveness of the structure for assembly, sequencing reads should contain solid k-mers corresponding to every position in the genome. Hence, either a high sequencing For instance, reads 4 and 5 are not indexed to the CAT entry because extensions T and TA are 1-redundant with respect to extension GA from read 1. Reverse complements of reads are also indexed, but are omitted in this gure.

coverage or a low error-rate is required. Both criteria are typically met with recent Illumina sequencers. For assembly, it can be veried that basic traversal of a string graph can be performed with this structure. The RRS acts as an incomplete inverted index for the reads. Specically, in the string graph, out-neighbors of a read (i.e., other reads that overlap that read to the right) are retrieved from the RRS by querying each of the read k-mers. In-neighbors (left overlap) are equivalent to out-neighbors of the read reverse complement.

Results

We developed an implementation of the on-line k-mers ltering and the reads indexing algorithms, as part of the Monument assembler [START_REF] Chikhi | Localized genome assembly from reads to scaolds: practical traversal of the paired string graph[END_REF]. The implementation is tested on two actual sequence datasets from R. sphaeroides (SRA reference SRR034530) and N.crassa (all libraries from [START_REF] Shea | ALLPATHS 2: small genomes assembled accurately and with high continuity from short paired reads[END_REF]) sequenced using the Illumina technology. The R. sphaeroides dataset (dataset 1) contains 46 million reads of length 36 bp. The N.crassa dataset (dataset 2) contains 320 million reads of average length 32 bp. Benchmarks were run on a 64-bit 8-cores machine with 66 GB of memory. In this implementation, read sequences are stored in memory on each node as an array of 2-bit encoded sequences. In the case of multi-nodes computation, n 4 bytes are redundantly stored per node, where n is the number of nucleotides in the reads. For the R. sphaeroides reads set, this amounts to 0.462 GB. We rst examined the eect of on-line k-mers ltering on the rst dataset. To this end, only the abundance count is retained for each k-mer. A comparison against a k-mer counting without ltering is made in Figure 2. It is important to note that, when entries corresponding to erroneous k-mers are removed from the hash table, the allocated memory is not freed but is instead made available for new entries. There are 144 M k-mers in the dataset, only 4.5 M (3.1%) of which are correct. On-line ltering enabled to keep the number of k-mers in the hash table under 23 M at any time. We veried that 4,544,973 solid kmers are retrieved without ltering, compared to 4,464,256 (98.2%) solid kmers with ltering (solid threshold m solid = 10). The dierence of 80,717 k-mers corresponds to premature ltering of k-mers that would be solid if given enough time before ltering. Then, we computed the full indexing time for an increasing number of cores (Figure 3). Some constant over-head occurs as reads pre-loading is not parallelized.

We compared memory usage of indexing procedures from other popular ultrashort reads assemblers with our implementation. The Velvet assembler (version 1.1.03) and the SOAPdenovo assembler (version 1.05) are based on de Bruijn graphs and use graph simplication heuristics. SOAPdenovo is specically optimized for memory eciency, it discards reads and pairing information in the initial graph structure. Our implementation uses spectrum parameter t = 4, S min = 10 6 , m solid = 10 and r = 0 for both datasets. All the assemblers are executed with k-mer size of 21. Only the indexing phase of assemblers were run (pregraph for SOAPdenovo, velveth for Velvet). Results are summarized in Table 1. The k-mers ltering step is essential in our method: complete indexing of Dataset 1 without k-mers ltering required 20.1 GB of memory. In terms of wall-clock time, these methods are comparable: for the largest dataset, SOAPdenovo and Monument completed indexing in respectively 41 and 64 minutes using 6 threads. In conclusion, our indexing scheme signicantly reduces the memory bottleneck for assembly, with minor impact on parallel indexing time. We veried that this index permits assembly with comparable quality than existing methods. To this end, we implemented an assembly method which constructs contigs based on extensions recorded in the representative reads spectrum. We computed the following N50 values (contig length such that longer contigs produce 50% of the assembly) with respect to the same assembly size (that of Velvet, 4.30 Mbp). Our assembly of R. sphaeroides yields 4.28 Mbp of contigs with a N50 value of 1.49 kbp. In comparison, we executed Velvet with similar k-mer size, which yields 4.30 Mbp of contigs with a N50 value of 1.41 kbp. In terms of running times, Velvet assembly phase (velvetg) executed in 2 min 46 sec CPU time, whereas our parallel string graph traversal required 1 min 42 sec CPU time per thread using 6 threads.

Discussion

A novel method is proposed for multi-nodes, multi-threaded reads indexing. It introduces two ltering techniques for memory eciency: on-line removal of erroneous k-mers and on-line indexing of only representative reads. To our knowledge, this is the rst read index that provides in-memory O(1) access to overlaps between reads, full read sequences and also pairing information. This novel indexing method allows to design assembly algorithms which overcome several of the performance issues inherent to graph-based assembly. For instance, memory usage is lowered as no graph structure is constructed. Independent indexing of sub-indexes allows for embarrassingly parallel and distributed computation. As implemented in our prototype, the read index is constructed using signicantly less memory than recent, optimized implementations of de Bruijn graphs with comparable indexing time.

These results can also be used to reduce memory usage of single-threaded greedy assemblers. Usually, the same data structure is used to construct the nal reads index and then access it during assembly. However, greedy assemblers typically uses an immutable index. Hence, taking advantage of immutability, one can focus on designing a more compact representation of the reads index once it is fully constructed. For memory eciency, sub-indexes can be simply constructed one at a time. In our index, the hash table can be replaced by a succinct rank/select data structure [START_REF] Okanohara | Practical entropy-compressed rank/select dictionary[END_REF] to represent entries, and a simple array containing xed-length lists of representative reads. The memory overhead of such structure becomes negligible as no pointer is used. A preliminary experiment shows that the entire index of dataset 2 is represented in only 4.2 GB of memory with this method.

Fig. 1 .

 1 Fig.1. The representative reads spectrum for a set of 8 reads with parameters k = 3, t = 1. Entries are solid k-mers from reads. Each k-mer is associated with a list of reads which extend the k-mer to the right. The extensions are ltered for t-redundancy. For instance, reads 4 and 5 are not indexed to the CAT entry because extensions T and TA are 1-redundant with respect to extension GA from read 1. Reverse complements of reads are also indexed, but are omitted in this gure.

Fig. 2 .

 2 Fig.2. Memory usage during the on-line k-mers ltering procedure, compared with un-ltered indexing. Dataset 1 is processed with parameters m solid = 10, r = 10, Smin = 10 7 and using 1 thread. The rst ltering pass is triggered at 11.6% of the dataset. Sporadic jumps in memory consumption correspond to resize operations of the hash table.

Fig. 3 .

 3 Fig.3. Execution time of indexing for our implementation on datasets 1 and 2 using 1 node and 1 to 8 threads.

Table 1 .

 1 Practical memory usage of indexing 46 M reads from R. sphaeroides (dataset 1) and 320 M reads from N.crassa (dataset 2) using Velvet, SOAPdenovo and Monument. Velvet exceeded the memory limit (66 GB) on the second dataset.

		Dataset	Monument	Velvet	SOAPdenovo
	Peak memory (GB)	1	2.7	7.7	3.9
		2	15.3	-	31.4


Boumedyen Boussaid
email: boussaid.boumedyen@cran.uhp-nancy.fr

Hristophe Aubrun
email: christophe.aubrun@cran.uhp-nancy.fr

M Ohamed

Naceur Abdelkrim
email: naceur.abdelkrim@enig.rnu.tn

Mohamed Koni

Ben Gayed
email: mbengayed@yahoo.fr

PERFORMANCE EVALUATION BASED FAULT TOLERANT CONTROL WITH ACTUATOR SATURATION AVOIDANCE

Keywords: fault tolerant control systems, performance degradation, reference management, model predictive control, performance index

In this paper, a new approach regarding a reconfigured system is proposed to improve the performance of an active fault tolerant control system. The system performance is evaluated with an intelligent index of performance. The reconfiguration mechanism is based on a model predictive controller and reference trajectory management techniques. When an actuator fault occurs in the system, a new degraded reference trajectory is generated and the controller calculates new admissible controls. A constraint set and cost function are established to avoid actuator saturation and reduce the control energy spent in closed loop dynamics. The effectiveness of the proposed method is illustrated using a hydrothermal system subject to actuator faults and constraints on actuator dynamic ranges.

Introduction

Industrial systems have become more sophisticated and complex. This complexity leads to an increased request for reliability, reconfigurability and safety of the systems. In order to guarantee the properties quoted above, it is essential to develop methods of supervision such as diagnosis and Fault Tolerant Control (FTC). Recently, the importance of FTC systems becomes increasingly apparent, and significant amount of research has already been done in this area (Patton, 1997;[START_REF] Korbicz | Fault Diagnosis: Models, Artificial Intelligence, Applications[END_REF][START_REF] Guerra | Robust fault detection with unknown input set-membership state estimators and interval models using zonotopes[END_REF][START_REF] Lunze | Reconfigurable fault-tolerant control: A tutorial introduction[END_REF]Zhang and Jiang, 2008;Noura et al., 2009;Puig, 2010).

According to a variety and severity of faults that may affect the system, different levels of performance have to be considered in different fault scenarios. From safety region to danger region, degraded performances are often acceptable. In addition, to ensure that the closed-loop system be able to track a command input or a reference model/trajectory even in the event of faults, a reconfigurable feed-forward controller has to be synthesized to achieve command tracking (Zhang and Jiang, 2003). In the case of performance degradation and actuator saturation avoidance being required, reference management may need to be used to adjust the command input or reference trajectory automatically or provide advisory information to human operators in the event of faults.

To the best of our knowledge, few works have published on the topic of reference management for fault tolerant systems. However, there are still many open issues which have to be addressed. [START_REF] Jiang | Graceful performance degradation in active fault tolerant control systems[END_REF]2006) have proposed a reference management approach based on graceful performance degradation and explicit model-following techniques. When a fault is detected by the Fault Detection and Diagnosis (FDD) algorithm, a reconfigurable controller is automatically designed so that the dynamics of the closed-loop system match those of the reference model known as the performance reduced reference model. Another case of multiple actuator failures is studied by Zhang et al. (2008). The approach proposed by Theilliol et al. (2008) is based on a modified recovery control system in which the reconfigurable reference input, applied when an actuator fault occurs in the system, is determined by considering the error between the reference and the output as an impulse signal. In the work of Theilliol et al. (2009), the modified reference trajectory is computed by minimizing the output-tracking error by a cost function based on the control energy concept.

Other approaches using set-point optimization based on Model Predictive Control (MPC) algorithms are discussed by [START_REF] Marusak | Actuator fault tolerance in control systems with predictive constrained setpoint optimizers[END_REF] as well as Tatjewski (2010), who propose a predictive constrained setpoint optimizer in order to deal with system faults. The capabilities of MPC algorithms to perform a reconfiguration action after a fault occurrence are quite limited by the accuracy of the FDD stage [START_REF] Maciejowski | Predictive Control with Constraints[END_REF][START_REF] Ding | A synthesis approach of on-line constrained robust model predictive control[END_REF][START_REF] Boussaid | FTC approach based on predictive governor[END_REF]. To deal with model uncertainty after a fault occurrence induced by the fault detection mechanism error, robust approaches have been proposed by [START_REF] Bemporad | Robust model predictive control: A survey[END_REF] as well as [START_REF] Cannon | Optimizing prediction dynamics for robust MPC[END_REF].

This paper proposes a method to reconfigure the reference trajectory based on performance evaluation with different levels of achievable performance in the presence of various faults under given potential system performance limitations. In the present work, the capabilities of MPC to handle on-line system constraints and reference trajectory management are exploited to design a fault tolerant controller.

The main contribution of this paper is to introduce a new methodology to generate an explicit degraded reference trajectory and an admissible control set after an actuator fault has occurred. A hydrothermal example is used to illustrate the concepts and the design procedures, and some simulation results are shown. The paper is organized as follows. Some performance backgrounds are considered in Section 2. Problem formulation is detailed in Section 3. Section 4 is dedicated to the design of the fault tolerant controller based on reference trajectory management and model predictive control. Numerical examples with simulation results are presented in Section 5 to illustrate the proposed scheme, followed by conclusions in Section 6.

Performance analysis principle

2.1. System operating modes. A fault tolerant system has three distinct operating modes [START_REF] Blanke | Diagnosis and Fault-Tolerant Control[END_REF]:

• Nominal mode: this mode corresponding to normal operating characteristics. In this mode, all objectives are assumed to be attainable. That means that, if nominal control signals are sent to the actuators, nominal performances will be provided.

• Faulty mode: after a fault occurrence, the system switches to an abnormal operating mode which is called the "Faulty mode". In this condition, the system provides faulty performances for nominal controls. The faulty mode may induce different system behavior depending on fault severity.

• Degraded mode: this mode is a temporary operating mode. In this mode, the system continues to work with acceptable objectives/performances. These performances are considered to be degraded and the control signals are supposed to be admissible.

The objective of fault tolerant control is to establish a strategy of control which has the property to limit or even cancel the effects of a fault on the performances of the plant. This strategy must modify the control structure according to the fault severity. In certain situations, particularly when the degradation of the performances is critical, fault tolerance should be achieved by modifying both the set point of the process and the control parameters in a way that releases constraints on the actuators. Then, the new operating points allow the system to recover performances as close as possible to the nominal ones [START_REF] Aubrun | Design of a supervised control system for sludge dewatering process[END_REF]Zerz, 2008).

Performance consideration.

Usually, the controller is designed for the faultless plant in order to meet given performance specifications in a closed loop. FTC has the ability to react on the existence of the fault in order to satisfy the declared performances. Production systems have some various objectives to reach. These objectives are usually expressed in terms of quality, energy cost, time constraints, etc. On the other hand, objectives should be reached under certain constraints related to

• systems stability,

• the error between the reference signal and the output signal which impacts, e.g., on the product quality,

• robustness against disturbances or parameter variations,

• fault tolerance.

Hence, a fault tolerant system works in a closed loop to maintain the performance specifications expected by the operator (Fig. 1). The system operations can be decomposed into three functioning modes: nominal, faulty and degraded. The 459 principle of performance degradation is illustrated in Fig. 2. Due to the occurrence of a fault, the system deviates from its nominal functioning mode, defined by a pair (π nom ,υ nom ), to a faulty operating point (π f ,υ f ).

The goal of the accommodation procedure is to determine a new control law that takes the degraded system parameters into account and drives the system to a new operating point (π deg ,υ deg ), with respect to the control constraints. Let us consider Π as the set of all possible performances π and υ as the pair (u,o),

Π {π : υ =(u, o) ∈ R m × R n }, (1)
where π, u and o are respectively the performance, control and objective vectors. It is easy to see that

Π=Π nom ∪ Π deg ∪ Π un , (2)
where Π nom , Π deg and Π un are depicted in Fig. 2. In order to explain the difference between each operating mode, we introduce the domain of performance coverage, D C , for a controlled and stable system, which is defined as

D C := {π ∈ Π:(u, o) ∈ (U, O)} , (3)
where U and O are respectively the control and objective sets. Figure 2 (a) shows that this domain is entirely included in the required performance region, Π nom .I n a faulty case (Fig. 2(b)), this domain is moved in space, which leads to cross the degraded and/or the unacceptable performance region depending to the severity of the fault. Let us consider the following sets for nominal and faulty cases:

D nom C {π ∈ Π nom :(u, o) ∈ (U nom ,O nom )} , (4) D f C {π ∈ Π:(u, o) ∈ (U f ,O f )} .
(5)

For a perfect design of the nominal controller, the domain of performance coverage should be completely included in the nominal performance region, D nom

C ∩ Π nom = D nom C .
In the faulty case, it can be decomposed into three partial domains, D nom C , D deg C and D un C ,w h i c h represent the nominal, degraded and unacceptable performance coverage domains, respectively:

D f C = D nom C ∪ D deg C ∪ D un C (6) with ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ D nom C = D f C ∩ Π nom , D deg C = D f C ∩ Π deg , D un C = D f C ∩ Π un .
The following different cases arise:

Case 1: D nom C = ∅.
In this case, the system is completely reconfigurable and the nominal objectives are achieved.

3: D nom C = ∅, D deg C = ∅ and D un C = ∅.
In this case, the fault is severe and the system should be stopped. Then, maintenance operations could be performed.

2.3.

Performance evaluation. The performance is evaluated by means of the following performance index:

ε π = o -π 2 o 2 , (7)
where o and π represent the objective and the measured performance, respectively. In practice, the performance of a system is measured at each sampling time k which induces the computation of the performance index at each sampling time. It can be noticed that, for an abrupt variation in the output, it is meaningful to comment this index. To solve this problem, we propose here to use a moving average of the measured performance π through a receding horizon h:

• simple moving average:

π(k)= 1 h h-1 i=0 π(k -i), (8)
• weighted moving average:

π(k)= h-1 i=0 α i π(k -i), (9)
with

h-1 i=0 α i =1,α i >α i-1 .
Here α i can be chosen as a geometric progression:

α i = q 0 • q i
with 0 <q<1 and

q 0 = 1 -q 1 -q h .
The sampled performance index becomes

ε π (k)= o(k) -π(k) 2 o(k) 2 . ( 10
)
Using the following thresholds: ε nom π , ε deg π and ε d π , the performance sets for each region could be defined as

• Π nom {π :0≤ ε π ≤ ε nom π }, • Π deg {π : ε nom π <ε π ≤ ε deg π }, • Π un {π : ε deg π <ε π }.
These thresholds are fixed by the operator according to the plant requirements such as safety, production quality, energy spent, etc.

Problem statement

Let us consider a nominal and faulty system described by the following state space representations:

S nom : x(k +1)=A n x(k)+B n u(k), y(k)=C n x(k), (11)
S f : x(k +1)=A f x(k)+B f u(k), y(k)=C f x(k), (12)
where x(k) ∈ R n is the state vector, u(k) ∈ R m is the input vector, y(k) ∈ R p is the output vector, and (A n ,B n ,C n ) represents the system in the nominal behavior and (A f ,B f ,C f ) represents the system in the faulty case.

First of all, it is supposed that the system is stable and controlled in three operating modes, and the structures (S nom and S f ) are known in nominal and faulty cases. It is also assumed that the state matrices (A n ,B n ,C n )a n d (A f ,B f ,C f ) are detectable in both cases. Thus, we can define the following sets: the set of the nominal controls, U nom , which correspond to the normal operating mode, and the set of the faulty controls, U f , which lead to abnor-mal functioning. The two sets can be written as

U nom u nom : u nom j,min ≤ u nom j ≤ u nom j,max , 1 ≤ j ≤ m, o ∈ O nom } , (13)
U f u f ∈ U nom : u f j,min ≤ u f j ≤ u f j,max , 1 ≤ j ≤ m, o ∈ O f , ( 14
)
where m is the number of actuators, (u nom j,min ,u nom j,max ) stand for the pair of upper and lower control bounds for each actuator j in the nominal mode, (u f j,min ,u f j,max ) stand for the pair of upper and lower control bounds for each actuator j in the faulty mode, o is the objective, O nom denotes the objective set in the nominal mode, O f means the objective set in the faulty mode.

For the sake of simplicity, only the static operating condition is considered. Without loss of generality, when this condition is established, it is assumed that y ⋍ y ref where y represents the output trajectory and y ref the reference trajectory that the system is supposed to be able to track. In the following, the output trajectory expresses explicitly a performance level and the reference trajectory reflects the objective that the system has to reach.

According to the previous set definitions, the following sets of nominal (respectively faulty) trajectories are defined:

Y nom y nom : u ∈ U nom ,y min j ≤ y ≤ y max j , 1 ≤ j ≤ m} , (15)
Y f y f : u ∈ U f ,y min j ≤ y ≤ y max j , 1 ≤ j ≤ m . (16)
In order to define the reference trajectories, the following sets are constructed:

Y nom ref {y nom ref : y ref ∈ D nom C } , (17)
Y f ref y f ref : y ref ∈ D f C , (18)
where y nom ref and y f ref denote the reference trajectory in nominal and faulty cases, respectively.

In the following study and for clearance, we suppose that the faults induce a slight shift of the performance coverage disc, which means that either D nom C or D deg C is nonempty. In addition, it is supposed that the degraded objectives are accepted and the corresponding controls are called admissible controls. The corresponding measured performances are also called degraded performances or degraded outputs. Hence, the set of these admissible controls, U adm , which satisfies the degraded objectives/outputs, Y deg , can be written as

U adm u adm ∈ U nom : u adm j,min ≤ u adm j ≤ u adm j,max , 1 ≤ j ≤ m, o ∈ O deg ,} (19) Y deg y deg : u ∈ U adm ,y min j ≤ y ≤ y max j , 1 ≤ j ≤ p} , (20)
where (u adm j,min ,u adm j,max ) is the pair of upper and lower admissible control bounds for each actuator j,a n dO deg is the degraded objective set.

Hence, the main idea is to determine a set of reference trajectories that allow the system under faulty conditions to operate with limited/degraded performances. Basically, the reference trajectories should be included in the degraded performance coverage domain in order to ensure the attainability of the objectives/performances. Thus, the degraded reference trajectory set, Y deg ref can be defined as

Y deg ref y deg ref : y ref ∈ D deg C . (21)

Fault tolerant controller design

4.1. FTC mechanism. The reconfiguration strategy of the proposed fault tolerant controller is described in Fig. 3, which includes modules of reference management, the MPC controller, FDD and the reconfiguration mechanism.

Plant Reference Management

Fault model.

In this study only actuator faults with the reduction of effectiveness are considered. In fact, such a fault can be considered a parametric fault modulated by a static or variable coefficient γ comprised between 0 and 1,

u f j (k)=(1-γ k j )u j (k), (22)
where γ k j ∈ [0, 1] and j is the actuator, j ∈{ 1,...,m}. In ( 22), the additive fault component which denotes a constant offset, occurs after actuator jamming is ignored. We only consider the loss of actuator effectiveness after fault occurrence.

Let Γ k ∈ R m × R m be the distribution matrix of actuator faults:

Γ k = ⎡ ⎣ γ k 1 00 0 γ k j 0 00γ k m ⎤ ⎦ .
Equation ( 22) can be written as

u f (k)=( I m - Γ k )u(k).W es e tB f = B n (I m -Γ k )=B n -B n Γ k and get x(k +1)=A f x(k)+B f u(k) =A n x(k)+B n u(k) -B n Γ k u(k) =A n x(k)+B n u(k)+f γ k ,
where f γ k ∈ R m is the vector of faults. The FDD module should estimate the fault magnitude f γ k , and then we deduce γk j for each actuator j.N otice that the FDD module is beyond the scope of our study.

Reconfiguration mechanism.

The reconfiguration module should find the pair (u,o) that makes the system operate in the nominal performance region or in the degraded performance region if the nominal ones are unachievable. Let us consider the worst case, where only the degraded performances could be met, and there exists ap a i r( u adm ,o deg ) such that Π deg is nonempty. In fact, u adm ∈ U adm and o deg ∈ O deg . Thus, the reconfiguration problem is reduced to determine the two sets U adm and O deg .

Admissible control set.

Let us start from the nominal control for an actuator j in the fault-free case:

u nom j,min ≤ u nom j ≤ u nom j,max , (23)
where the upper and lower nominal control limits are completely known after the design of the nominal controller.

Assuming that the fault occurs at the time k f , and is detected at the time k d with the estimated fault magnitude γk d ,fork>k d ,wehave

(1 -γk d j )u nom j,min ≤ (1 -γk d j )u nom j ≤ (1 -γk d j )u nom j,max , (24
) where γk d j is a positive real number and γk 23) and ( 24), we get

d j ∈ [0, 1] so that (1 -γk d j ) ∈ [0, 1]. Due to (
(1 -γk d j )u nom j,min ≤ u j ≤ (1 -γk d j )u nom j,max . (25)
We use the following notation

⎧ ⎨ ⎩ u adm j,min =(1-γk d j )u nom j,min , u adm j,max =(1-γk d j )u nom j,max .

B. Boussaid et al.

Thus ( 25) can be re-written as

u adm j,min ≤ u j ≤ u adm j,max . (26)
Besides, the admissible control should avoid any actuator saturation. Hence

u adm j = σ j (u j )
such as

σ j (u j )= ⎧ ⎪ ⎨ ⎪ ⎩ u adm j,max if u adm j,max ≤ u j , u j if u adm j,min ≤ u j ≤ u adm j,max , u adm j,min if u j ≤ u adm j,min . (27)
The admissible control set is determined by the definition of the upper and lower control limits. These constraints should be included in the MPC problem formulation in order to avoid any actuator saturation and to let the MPC optimizer find the appropriate control signal u adm .

Degraded objective set.

On the analogy of ( 7), the following index can be used, where o represents the objective:

ε o = o -o 2 o nom 2 . ( 28
(1 -ε deg π ) o nom 2 ≤ o 2 ≤ (1 + ε deg π ) o nom 2 (29)
Finally, the degraded objective set, O deg ,isgivenby

O deg o deg :(1 -ε deg π ) o nom 2 ≤ o deg 2 ≤ (1 + ε deg π ) o nom 2 . (30)
The optimum degraded objective should be as closed as possible to the nominal objective or the required performance. Then, it can be obtained by minimizing the following quadratic function under constraints on the control and objective:

o * deg =arg min u adm ∈U adm o deg ∈O deg o deg -o nom 2 . (31)
4.4. Reference management. As mentioned previously, the following function is defined to generate the reference trajectory input:

y ref (k)=o(k) -τ (o(k) -y ref (k -1)) . (32)
where y ref (0) means the initial condition and τ represent the system dynamics, with τ ∈ [0, 1] , a positive real number. After fault detection and diagnosis, the reconfigured reference trajectory, y deg ref ,is

y deg ref (k)=o * deg (k) -τ o * deg (k) -y deg ref (k -1) . (33)
4.5. Fault tolerant MPC controller. The control signal delivered by the MPC controller is calculated by minimizing the following cost function J with the prediction horizon h p :

J =m i n u(k),...,u(k+hp-1) hp-1 i=0 (y ref (k + i) -y(k + i)) T Q (y ref (k + i) -y(k + i)) + u(k + i) T Ru(k + i) , (34)
subject to the constraints

y min (k + i) ≤ y(k + i) ≤ y max (k + i), u min (k + i) ≤ u(k + i) ≤ u max (k + i), |∆u(k + i)|≤∆u max (k + i),
where y(k + i) is the output, u(k + i) the control and ∆u(k + i) the change in the control action from one instant to another, for each discrete time (k + i).T h et e r m s max and min refer to upper and lower limits, respectively. The weighting matrices Q (Q = Q T )andR (R = R T )of the cost function are positive definite and they are used to tune the MPC control law.

After fault occurrence, the cost function becomes

J deg =m i n u adm (k),...,u adm (k+hp-1) hp-1 i=0 y deg ref (k+i)-y(k+i) T Q y deg ref (k + i) -y(k + i) +u adm (k + i) T Ru adm (k + i) , (35)
subject to the constraints

y min (k + i) ≤ y(k + i) ≤ y max (k + i), u adm min (k + i) ≤ u adm (k + i) ≤ u adm max (k + i), |∆u(k + i)|≤∆u max (k + i).

Numerical example

The system to be studied is a hydrothermal process which consists of a tank with natural racking and electric heating [START_REF] Boussaid | A fault tolerant predictive control applied to hydro-thermal system[END_REF]. In this plant, the fluid is introduced into the tank with a controlled flow q e ,u s i n gt h e control input u q . The fluid is evacuated with the flow q s through a manual valve V m . The electric heating system allows the temperature of the fluid to increase by modulation of the current in resistance R using the order u e .A n agitator is used to homogenize the temperature of the fluid in the tank. This system is shown in Fig. 4. The characteristics of the installation are as follows:

HIC

• the tank has a section S of 0.5 m 2 ,

• the resistance R is 2.42 Ω,

• the level sensor used delivers 1 V/m,

• the temperature gauge used delivers 10 mV/ • C,

• the gain of the electro-valve K e is 2 l/mn/V,

• the manual valve V m is partially open, q s =0.2L,

• the sampling period T ech is 5 s,

• L denotes the level and T denotes the temperature.

The nominal system S nom :( A n ,B n ,C n ) is fully described by the following matrices:

A n = 0.9967 0 00 .1353 , B n = 0.0511 -0.0256 00 .0169 , C n = 0.0078 0 00 .0085 .

The control and output vectors are u = u e u q T and y = T s L T .

The calculation of the two controls u e and u q is performed by the MPC controller. The parameters of the function cost include

• the prediction horizon: h p =2,

• the output constraints:

{0 ≤ T s (k) ≤ 100} [ • C] and {0 ≤ L(k) ≤ 10}[m],
• the control constraints:

{0 ≤ u e (k) ≤ 10}[V], {0 ≤ u q (k) ≤ 10} [V] and |∆u(k)|≤0.5 [V],
• the relaxation matrices:

Q =10 -3 10 01 , R =10 -1 10 01 ,
• the nominal objectives (required performances):

o nom T =20[ • C] and o nom L =8[m].
In order to define the different performance regions, the following performance index thresholds are considered: The following figures show the results of simulations performed in MATLAB. Figure 5 presents the response of the nominal plant (fault-free system). It shows the two outputs, the tank temperature [ • C] and level [m], and their corresponding control signals, i.e., the electric valve and resistance control signals [V]. Notice that the nominal reference trajectories and the control limits are represented here with dashed lines.

ε nom π =5%,ε deg π = 20%, O nom o nom = o nom T o nom L T : 19 ≤ o nom T 2 ≤ 21, 7.6 ≤ o nom L 2 ≤ 8.4 , O deg o deg = o deg T o deg L T : 16 ≤ o deg
In the experiment, a fault on the controlled valve actuator occurs at time 50. It is assumed that the reconfiguring action starts three seconds after the fault occurrence. To study the impact of the fault on the system performance, three cases are considered as follows. The degraded objective is The effects of the fault on system performance are shown in Fig. 6. After fault occurrence, the valve control remains within the admissible limits and the performance index slightly deviates and returns to its nominal value. The large magnitude of the performance index between time 0 and time 20 corresponds to a transient response of the system.

o deg * T =arg min u adm ∈U adm o deg T ∈O deg o deg T -20 2 =20, o deg * L =arg min u adm ∈U adm o deg L ∈O deg o deg L -8 2 =8.
o deg * L =arg min u adm ∈U adm o deg L ∈O deg o deg L -8 2 =7.3.
In this case, the fault is more severe than in Case 1. It can be noticed that the performance index (shown in Fig. 7) does not return to the nominal value. Then, the nominal objective cannot be achieved. Nevertheless, the new degraded reference applied to the system allows the valve control to remain within the admissible limits and the level value follows the new set-point. The deterioration of the performance is critical since the performance index, shown in Fig. 8, is outside the degraded performance region. It is clear that the controls still operate on the admissible interval values. The evolution of the performance index with respect to the actuator fault magnitude, shown in Fig. 9, is valuable information which helps the operator in decision making. For this specific case, the degraded performance region corresponds to the fault magnitude comprised between 0.24 and 0.36.

Conclusion

This paper presented an active fault tolerant control system design strategy. The strategy is based on the computation of a degraded reference trajectory with respect to an admissible control set. The new reference trajectory is determined on the basis of an optimization algorithm which ensures the lower possible performance degradation. The impact of the faults on the system performances is evaluated by a specific index. This index provides useful information to the operator for the supervision of the process.

Fig. 1 .

 1 Fig. 1. Fault tolerant system closed loop.

Fig. 2 .

 2 Fig. 2. Domain of performance coverage in nominal and faulty mode functioning.

  Case 2: D nom C = ∅ and D deg C = ∅. This means that the nominal objectives are not achievable. Only the degraded objectives are met if they are accepted by the supervision instructions.

Case

Fig. 3 .

 3 Fig. 3. Proposed FTC reconfiguration mechanism.

  ) with o ∈ O deg , which means that ε o ≤ ε deg o where ε deg o is a threshold fixed by the operator. In practice, the value of the threshold is defined as ε deg π . ε o ≤ ε deg π leads to the following bounds:

  Fig. 4. System description.

  ,k ≥ 50. The admissible control set is U adm = u adm = u adm e u adm q T : 0 ≤ u adm e ≤ 10, 0 ≤ u adm q ≤ 8, 4 .

Fig. 5 .

 5 Fig. 5. Nominal plant outputs and controls.

Fig. 6 .

 6 Fig. 6. Low severity fault effect.

  Case 2: γk 1 =0for all k andγk 2 = 0,k < 50, 0.27,k ≥ 50.The admissible control set isU adm = u adm = u adm e

  Fig. 7. Medium severity fault effect. Case 3: γk 1 =0for all k and

Fig. 8 .

 8 Fig. 8. High severity fault effect.

Fig. 9 .

 9 Fig. 9. Evolution of the performance index.

 Noura, H., Theilliol, D., Ponsart, J. and Chamssedine, A. (2009). Mohamed Koni Ben Gayed received a Ph.D. in automatic and industrial data processing from the University of Lille, France, in 1987. He is currently a member of the research team on Modeling, Analysis and Control of Systems (MACS). His field of research is automation, including control, supervision, safety as well as man/machine interaction. He also teaches with the ISET of Gabès, Tunisia.