Nicolas Obin

Xavier Rodet

Anne Lacheret-Dujour

A Syllable-Based Prominence Detection Model Based on Discriminant Analysis and Context-Dependency

On the basis of our previous work, we propose a syllablebased prominence detection model within the framework of exploratory data analysis and discriminant learning in the acoustic domain. This paper investigates two hypothesis on the acoustic data processing: a linear discriminant analysis in which the relative discriminant ability of single prosodic cues are combined into prosodic patterns and a context-dependant model that accounts for phonological dependencies (phonetic intrinsic properties and coarticulation effect). The proposed approach significantly outperforms a baseline method on a corpus of French read speech with a performance of 87.5% in f-measure for the prominent syllables (respectively 90.4% in global accuracy).

Introduction

To enable the development of linguistic study and prosodic models over large corpora, a number of tasks should be automatically preprocessed such as phonetic decoding and segmentation, detection of para-verbal phenomena (fillers), prosodic saliencies (prominences) and disfluencies. Research on prominence (syllable or word based) has received a significant amount of attention over past years both in the study of its perception, its acoustic correlates as its automatic detection. We will focuse here mostly on this acoustic dimension. The perception of prominence is complex: it results from the interaction between acoustic perceptual salience (a bottom-up process involving significant acoustic deviation from more or less abstract references) and cognitive representations and associated expectations (a top-down process which involves lexical feature, syntactical structure, and informational structure) [START_REF] Erikson | Perception of syllable prominence by listeners with and without competence in the tested language[END_REF][START_REF] Wagner | Great expectations -introspective vs. perceptual prominence ratings and their acoustic correlates[END_REF]. Most of the studies on prominence, if they do not deny its topdown dimension, have nevertheless focused on its acoustic correlates [START_REF] Terken | Fundamental frequency and perceived prominence[END_REF][START_REF] Mertens | Local prominence of acoustic and psychoacoustic functions and perceived stress in french[END_REF][START_REF] Campbell | Loudness, spectral tilt, and perceived prominence in dialogues[END_REF][START_REF] Sluijter | Spectral balance as an acoustic correlate of linguistic stress[END_REF][START_REF] Gussenhoven | The perceptual prominence of fundamental frequency peaks[END_REF][START_REF] Fant | Prominence correlates in swedish prosody[END_REF][START_REF] Erikson | Syllable prominence: A matter of vocal effort, phonetic distinctness and topdown processing[END_REF][START_REF] Heldner | On the reliability of overall intensity and spectral emphasis as acoustic correlates of focal accents in swedish[END_REF]. These studies have pointed out that the perception of prominence is based on a set of acoustic correlates that seems to be independent of language (variations in duration, pitch [START_REF] Terken | Fundamental frequency and perceived prominence[END_REF][START_REF] Gussenhoven | The perceptual prominence of fundamental frequency peaks[END_REF], intensity, and vocal quality (spectral emphasis, vocal effort) [START_REF] Sluijter | Spectral balance as an acoustic correlate of linguistic stress[END_REF][START_REF] Campbell | Loudness, spectral tilt, and perceived prominence in dialogues[END_REF][START_REF] Erikson | Syllable prominence: A matter of vocal effort, phonetic distinctness and topdown processing[END_REF][START_REF] Heldner | On the reliability of overall intensity and spectral emphasis as acoustic correlates of focal accents in swedish[END_REF]), but whose relative contributions is language-specific [START_REF] Mertens | Local prominence of acoustic and psychoacoustic functions and perceived stress in french[END_REF]. Another property is the acoustical context-dependencies which account for the local acoustic reference baseline used by listeners to judge perceptual saliency [START_REF] Terken | Fundamental frequency and perceived prominence[END_REF][START_REF] Mertens | Local prominence of acoustic and psychoacoustic functions and perceived stress in french[END_REF][START_REF] Gussenhoven | The perceptual prominence of fundamental frequency peaks[END_REF]. Such findings, while they are fundamental to understand the mechanisms that underlay the perception of prominence in speech, could be of little interest in an automatic detection task. There are two reasons for this: 1) a significant difference observed between acoustic means does not ensure a good separability for each particular observation. 2) the use of controlled utterances neutralizes the phonetic contextual factor which is a source of variability in the acoustic observations (acoustic properties intrinsic units observed and coarticulation) in speech. These two reasons combined together explain why a relevant acoustic cue in the case of laboratory speech may have a low discriminating ability or difficult to integrate in a prominence automatic detection task [START_REF] Rosenberg | Detecting pitch accent using pitch-corrected energy-based predictors[END_REF]. In parallel, methods have been proposed for the automatic detection of prominence [START_REF] Jenkin | Development and comparison of three syllable stress classifiers[END_REF][START_REF] Tamburini | Automatic detection of prosodic prominence in continuous speech[END_REF][START_REF] Sun | Pitch accent prediction using ensemble machine learning[END_REF][START_REF] Kochanski | Loudness predicts prominence: fundamental frequency lends little[END_REF][START_REF] Ananthakrishnan | An automatic prosody recognizer using a coupled multi-stream acoustic model and a syntactic-prosodic language model[END_REF][START_REF] Avanzi | Méthodologie et algorithmes pour la détection automatique des syllabes proéminentes dans les corpus de franc ¸ais parlé[END_REF][START_REF] Rosenberg | Detecting pitch accent using pitch-corrected energy-based predictors[END_REF][START_REF] Obin | Prominence model: a probabilistic framework[END_REF]. Research in the field of automatic detection of prominence has mainly concentrated on the search and integration of discriminant acoustic features [START_REF] Tamburini | Automatic detection of prosodic prominence in continuous speech[END_REF][START_REF] Kochanski | Loudness predicts prominence: fundamental frequency lends little[END_REF][START_REF] Rosenberg | Detecting pitch accent using pitch-corrected energy-based predictors[END_REF][START_REF] Obin | Prominence model: a probabilistic framework[END_REF], relatively few on the choice of classifiers. On the acoustic level: pitch (f0) and duration feature are known and commonly used. If these correlates show robust performance in a prominence detection task, it seems that they have limited discriminant ability and are not sufficient in a significant number of cases. Therefore, more recent studies have tried to investigate some new acoustic correlates such as energy (acoustic and perceived) [START_REF] Kochanski | Loudness predicts prominence: fundamental frequency lends little[END_REF], energy bands [START_REF] Rosenberg | Detecting pitch accent using pitch-corrected energy-based predictors[END_REF], and spectral emphasis. These methods are based on the discriminant ability of a set of single acoustic cues; [START_REF] Tamburini | Automatic detection of prosodic prominence in continuous speech[END_REF] has proposed to combine a set of acoustic cues in order to cumulate their relative contributions in the achievement of an acoustic salience. However, the proposed method is based on some heuristics postulated solely on the basis of expert knowledge. On the classifier level: if the decision tree is referred as a baseline classifier, some other classifiers have been proposed through literature (neural networks, HMM, rule-based [START_REF] Jenkin | Development and comparison of three syllable stress classifiers[END_REF], bagging and boosting approach [START_REF] Sun | Pitch accent prediction using ensemble machine learning[END_REF], coupled-HMM [START_REF] Ananthakrishnan | An automatic prosody recognizer using a coupled multi-stream acoustic model and a syntactic-prosodic language model[END_REF], voting classifiers [START_REF] Rosenberg | Detecting pitch accent using pitch-corrected energy-based predictors[END_REF]). Following our previous work [START_REF] Obin | Prominence model: a probabilistic framework[END_REF], our approach is based on exploratory data analysis and discriminant learning in the acoustic domain. Our method aims to automatically determine from a large dimensional acoustic space a set of discriminant prosodic cues and/or patterns from a large dimensional acoustic space in a prominence detection task. In the present study, the authors propose to investigate two new assumptions:

• H1: we support that the phenomenon of prominence results from a complex mechanism of interaction of different prosodic cues. This interaction takes place both in the combination of heterogeneous acoustic features (pitch, duration, vocal quality, etc.) and in the acoustic context-dependencies (combination of multiple reference floors according to multiple context-window sizes). • H2: we suggest that taking into account the phonetic-context, by neutralizing the variations related to the contextual pho-netic properties (acoustic intrinsic properties and coarticulation effect), will improve the acoustics separation of prominent classes (P) and non-prominent (NP) syllables;

The article is organized as follows: first, we briefly describe the corpus used for the experiment; then, we describe the scheme of our proposed method; finally, we introduce our evaluation protocol and present the obtained results.

Speech Material

Unlike for English studies (Boston University Radio Speech Corpus [START_REF] Ostendorf | The boston university radio news corpus[END_REF], and Boston Directions Corpus [START_REF] Nakatani | Discourse structure in spoken language: studies on speech corpora[END_REF]) there is currently no reference corpus available for prominence detection in French. Therefore a home-made corpus of a single speaker French read speech has been used for this preliminary study.

The corpus has been manually annotated for prominence by 2 non-expert annotators with an agreement score of 78% of fmeasure (respectively 80% of accuracy). This score is in concordance with the inter-annotator agreement usually observed in the literature (80-90% accuracy). For the detection task, only syllables for which agreement occurs were considered as being prominent; others have been declared non-prominent. This results in 1670 prominent syllables and 4635 non-prominent syllables.

Exploratory and Discriminant Learning Framework

Feature extraction

In [START_REF] Obin | Prominence model: a probabilistic framework[END_REF], we have proposed a feature extraction framework for the prominence detection task. We outline here the main principles of the feature extraction step: we represent the acoustic space into 5 of his prosodic dimensions: pitch (f0), durational (syllable duration and local speech rate), intensity (absolute and perceived), spectral (mel-frequency energy and loudness), and vocal quality (spectral slope, spectral emphasis, open quotient of the glottal source) features. On the basis of these prosodic dimensions, low-level statistical syllable-based features are extracted: minimum, maximum, mean, slope, glide, excursion, range, and standard deviation. These features are then used to give information on 2 different levels:

• intra-syllabic features: describe prosodic cues that occurs inside a given syllable;

• extra-syllabic features: describe the acoustic-context dependencies such as context-window and reference floor that account for acoustic-context.

In this experiment, two syllable-based units (syllable and nucleus of the syllable), several contextual-windows (none, {previous; next; surrounding } syllable, current breath group, current utterance) and two reference floors (minimum and mean value on the considered context-window) have been used.

Discriminant Analysis

The aim of our proposed method is to automatically estimate and accumulate the discriminant ability of the single prosodic cues of the considered acoustic space in order to improve performance in detecting prominence (H1). To formalize this hypothesis, we set our approach within the framework of discriminant analysis. However, as the feature extraction step described above leads to a high dimensional feature space, this could causes several problems in a classification task. In particular: poor classification performance due to the fact that some features are irrelevant for the task and therefore introduce noise in the learning procedure; over-fitting of the model to the learning set, especially in the case of dimension reduction algorithm such as LDA.

To prevent these bad properties, we propose the following learning scheme: 1) the initial feature space is reduced by feature selection into a subset of uncorrelated and discriminant features;

2) the resulting feature sub-space is transformed by means of linear discriminant analysis (LDA) in order to accumulate the discriminant ability of each single prosodic cue.

Let us introduce some notations and definitions. Let K be the total number of classes, N k the number of total feature vectors accounting for the training data from class k and N the total number of feature vectors, µ k and Σ k respectively the mean vector and the covariance matrix of the class k, and µ the overall mean vector. The within and between class scatter is then defined as follows:

B = P K k=1 N k (µ k -µ)(µ k -µ) T W = P K k=1 P N k n k =1 (xn k -µ k )(xn k -µ k ) T
And the inertia ratio is defined as:

r = B W

Feature Selection

The method used for the feature selection step is based on Inertia Ratio Maximization using Feature Space Projection (IRMFSP ) [START_REF] Peeters | Automatic classification of large musical instrument databases using hierachical classifiers with inertia ratio maximization[END_REF][START_REF] Obin | Prominence model: a probabilistic framework[END_REF]. Let F be the acoustic feature space with respective feature dimensions (fi) 1≤i≤D where D is the rank of the features space. At each iteration of the process, the best feature fi is selected as the one maximizing ri (inertia ratio along dimension i) and then the feature space is projected along the dimension fi thus defining new dimensions fi that verify:

f j = f j -(f j . f i f i ) f i f i , ∀f j ∈ F
This projection step ensures that selected features are little corellated thus avoiding redundancy in the resulting features set. Iteration is processed until the ratio of the inertia ratio at step j to the inertia ratio of the first selected feature drop below a given threshold or when a desired number of features is reached. This leads to select a set of little-correlated features that individually maximizes the classes separability conditionally to each previous subset of selected features.

Feature Transform

The selected features estimated in the previous step are then used as input features of a linear discriminant analysis. The aim of this step is to combine the discriminant ability of the feature, by estimating the linear combinations of single features of the original space that maximize the class separation (following the same criterion mentioned above) This problems turns into determining the K-1 linear combinations (α) 1≤k≤K-1 that maximize the objective function:

J(α) = α T Bα α T W α
Once the parameters α have been estimated, the original features are projected along these directions into a transformed space of rank K-1.

Context-Dependant Analysis

Even if context-dependent analysis has already been used by [START_REF] Jenkin | Development and comparison of three syllable stress classifiers[END_REF] in prominence detection, we would like to assume it as a full hypothesis in the acoustic modeling of prominence. In order to account for the effect of phonetic-context (intrinsic properties and coarticulation) on the observed prosodic parameters (H2), we introduced syllable-based left-to-right contextual factors. The phonologic contextual factors used in this experiment are:

• type of {previous, current, next} syllable;

• number of phonemes in {previous, current, next} syllable;

• phonetic label of the current syllable's nucleus;

• phonetic class of the current syllable's nucleus;

• number of phonemes of {previous; current; next} syllable onset;

• number of phonemes of {previous; current; next} syllable coda;

• phonetic class of the current syllable onset & coda (i.e glide, occlusive, fricative, liquid, nasal);

Models evaluation 4.1. Compared Models

In this study, we compared 5 types of models of prominence detection. In order to remove any influence of the classifier type on the performance measurement, each experiment was performed using weka's J48 decision tree [START_REF] Witten | Weka: Practical machine learning tools and techniques with java implementation[END_REF]. The compared models are:

• a baseline model (BL): decision tree learned on the whole acoustic space.

• a Context-Dependent model (H1): decision tree learned on the whole acoustic space + phonetic-context features.

• an IRMFSP model: decision tree learned with the N first discriminant features selected with IRFMSP.

• an IRMFSP + LDA model (H2): decision tree learned on the optimal linear combination of the N first discriminant features selected with IRMFSP.

• an IRMFSP + LDA + context-dependent (H1 + H2): decision tree learned on the optimal linear combination of the N first discriminant features selected with IRFMSP + phoneticcontext feature

Evaluation scheme

The evaluation was conducted on the read speech corpus presented in section 2 within a 10-fold cross validation. Unlike other studies that use global accuracy as a measure of prominence detection performance, in this study we used the f-measure of the prominence class only. This choice is motivated by the following reasons: 1) prominence classification task is a detection task, i.e. a classification task in which only one of the class is of interest. This choice is all the more justified in the prominence detection task since the class of interest has much less observations than the other. 2) the f-measure provides a good compromise between recall and accuracy measures (insertion and deletion rates in the case of a two-class classification task).

Results and Discussion

In table 1, we summarize the results obtained with the 5 proposed models. Our proposed hypothesis (H1: phonetic contextdependency and H2: prosodic cues combination) leads to both significant improvement compared to the baseline method, with respective performance of 84.1 and 87.5% f-measure on the prominence class (this represents for comparison with standard performance measure respectively an accuracy of 86.5 and 90.4%). From the prosodic cue combination hypothesis (H2), we have investigated several properties compared to the baseline as well as the IRMFSP models. First we have investigated the discriminant ability of the most discriminant prosodic feature (syllable duration) compared to the combination of the N most discriminant features (H2). An anova analysis reveals that combination of prosodic features is significantly more discriminant for any number of selected features (Fstats = 385 for the most discriminant features vs. 511 ≤ F -stats ≤ 990 for the combination of the N most discriminant features N ≥ 2 on the test set) . A comparison of the IRMFSP and IRMFSP+LDA models by mean of prominence detection performance is presented in figure 1. Both models outperform the baseline model with a very small feature space dimension. The IRMFSP+LDA model outperforms the IRMFSP for almost all number of selected features (N > 2). This means that our proposed model succeed in finding discriminant prosodic patterns that are more discriminant than single cues estimated with the IRMFSP method. These results support the hypothesis that prosodic cues combine their contributions to achieve an acoustic saliency.

Models

From the context-dependant hypothesis, it can be seen that the increase in performance observed with our H1 model compared to the BL model and the combined H1 + H2 model compared to the H2 model is relatively small or even null. This means that the proposed model somehow failed in catching discriminant phonetic distinctions in prominence structure. This suggests that the phonetic distinction does not affect overall prosodic features homogeneously; phonetic-context rather affects each particular prosodic cue depending on his nature (f0, duration, intensity, spectral feature, voice quality). Thus, distinct context-dependent models should be learned for each prosodic cue or at least for each prosodic dimension. This remains to be investigated in a further study.

Conclusion and future works

In this paper we have proposed the investigation of two hypotheses for a syllable-based prominence detection model within a discriminant-based exploratory data analysis framework. These hypotheses were phonetic context-dependencies and prosodic cues combination. They both significantly outperform a baseline decision tree model in prominence detection. The improvement obtained by combining prosodic cues gives strong support to the hypothesis that prominence perception results from the combination of multiple prosodic cues contributions. However, our context-dependant model could be greatly improved by learning different context-dependant models according to the considered prosodic feature. Moreover, since the discriminant ability of a given prosodic cue is somehow related to his phonetic context, the context-dependant model should be estimated before estimating the discriminant ability of any prosodic cue; thus the context-dependency model and the feature selection should be merged into a single process. This preliminary study was achieved on read speech, further experiments will now be carried out on spontaneous speech corpora in order to propose a model robust to speaker, speaking style, expressivity and recording conditions. This will be carried out via the elaboration of a corpus of six hours of spoken French sampled into different speaking styles with manual prosodic annotations (Rhapsodie project). Furthermore, our proposed approach can be applied to any language since our model automatically estimates and combines relevant prosodic cues used for prominence achievement without any language-specific assumption. Thus it could be used in other languages to compare the language-specific prominence acoustic correlates and their relative weights in the achievement of an acoustic saliency.

Figure 1 :

 1 Figure 1: Performance of the IRMFSP (white bars) and the IRMFSP/LDA (gray bars) models as a function of the number of selected features in the feature selection step. Performances are presented as mean and standard deviation of the f-measure on the prominence class. The dotted line gives the mean performance of the baseline model. For convenience, performance is shown only for the first 15 selected features.

Table 1 :

 1 Performance of the compared models. In the case of the IRMFSP-based model, only the optimal model's performance is referred.

						Mean Performance ( % f-measure)
			BL			81.5
			H1			84.1
			IRMFSP		84.5
			H2			87.5
			H1 + H2	87.4
		88
		86
		84
	% performance	78 80 82
		76
		74
		72
		1	2	3	4	5 # of selected features 6 7 8 9 10 11	12	13	14	15

Acknowledgments

This study was supported by: • ANR Rhapsodie 07 Corp-030-01; reference prosody corpus of spoken French; French National Agency of research (ANR); 2008-2012., • Programmes Exploratoires Pluridisciplinaires (PEPS), CNRS/ST2I, 2008-2010.


Nicolas Obin
email: nobin@ircam.fr

Xavier Rodet
email: rodet@ircam.fr

Anne Lacheret-Dujour

A Multi-Level Context-Dependent Prosodic Model Applied to Durational Modeling

Keywords: speech synthesis, prosody, multi-level model, context-dependent model

We present in this article a multi-level prosodic model based on the estimation of prosodic parameters on a set of well defined linguistic units. Different linguistic units are used to represent different scales of prosodic variations (local and global forms) and thus to estimate the linguistic factors that can explain the variations of prosodic parameters independently on each level. This model is applied to the modeling of syllablebased durational parameters on two read speech corpora -laboratory and acted speech. Compared to a syllable-based baseline model, the proposed approach improves performance in terms of the temporal organization of the predicted durations (correlation score) and reduces model's complexity, when showing comparable performance in terms of relative prediction error.

Introduction

Research on speech synthesis has brought significant improvements over the past decade that makes possible to generate natural speech from text. However, if the synthesized speech sounds acoustically natural, it is often considered poor accordingly to the way of speaking (prosodic artefacts and monotony). Now modeling the variability in the way of speaking (variations of prosodic parameters) is required to provide natural expressive speech in many applications of high-quality speech synthesis such as multi-media (avatar, video game, story telling) and artistic (theater, music) applications. Despite growing attention to prosody modeling over the past few years, one of the major drawback of actual prosody models remains the monotony of the generated prosodic parameters. That problem (averaging problem) appears mainly due to two inter-related causes : a lack of linguistic knowledge extracted from the text that could be used to explain more accurately the observed prosodic variations and a representation of the prosodic parameters which is essentially a superposition of acoustic forms observed on a set of different + / -well defined linguistic units related to different communication functions.

A short overview of linguistic units and associated +/-linguistic parameters affecting these units leads to distinguish at least three main classes of prosodic scales and associated linguistic parameters : 1) global variations : overall properties of prosodic parameters (mean, variance). This level is associated with global characteristics of a speaker and to speaking styles associated with specific discourse genres. 2) local variations : variations on the smallest linguistic units (sub-state of a phoneme, syllable or phoneme). These variations are associated with phonological properties of these levels (coarticulation, syllabic structure, accentuation).

3) intermediate variations : variations on a set of units larger than the syllable and + / -linguistically well defined (accentual group, interpausal group, prosodic group, intonational phrase, period, verbal construction, discourse sequence, ...) and associated with + / -linguistic factors : physiological (f0 declination), modalities (questions, ...), syntactical (prosodic contrasts related to some specific syntactical sequence), semantic (informational structure) and discursive.

In order to model the variability of these parameters, it is necessary to determine an appropriate representation of prosodic parameters and then to extract and estimate the effects of high-level linguistic features (syntactic, semantic, discursive) on the observed prosodic parameters. At the signal level several approaches have been proposed to represent the variations of the fundamental frequency, HMM phone state-based [START_REF] Yoshimura | Simultaneous modeling of spectrum, pitch and duration in hmm-based speech synthesis[END_REF], syllablebased [START_REF] Boidin | Generating intonation from a mixed cart-hmm model for speech synthesis[END_REF], a set of well defined linguistic units [START_REF] Latorre | Multilevel parametric-base f0 model for speech synthesis[END_REF], or a set of units that are estimated with unsupervised statistical methods [START_REF] Morlec | Génération multiparamétrique de la prosodie du français par apprentissage automatique[END_REF][START_REF] Holm | Sfc : un modèle de superposition de contours multiparamétriques pour la génération automatique de la prosodie -apprentissage automatique et application à l'énonciation de formules mathématiques[END_REF]. For duration, a phone-based [START_REF] Zen | Hidden semi-markov model based speech synthesis[END_REF][START_REF] Gao | Duration refinement by jointly optimizing state and longer unit likelihood[END_REF][START_REF] Yamagishi | Phone duration modeling using gradient tree boosting[END_REF] or a syllable-based [START_REF] Chen | A new duration modeling approach for mandarin speech[END_REF][START_REF] Rao | Modeling durations of syllables using neural networks[END_REF] representation. If the linguistic studies define the syllable or syllable-like [START_REF] Barbosa | Generating segmental duration by p-centers[END_REF] unit as the minimal unit of prosodic rhythm, some studies have shown that rythm variations were also used as prosodic cues by speakers on larger units such as speech rate variations on some specific verbal constructions (oral parenthesis [START_REF] Gachet | Les parenthèses en français : Etude prosodique[END_REF]) as well as informational and discourse structure [START_REF] Koopmans-Van Beinum | Relationship between discourse structure and dynamic speech rate[END_REF]. Thus in the similar manner as for fundamental frequency variations, it could be usefull to represent durational variations on different linguistic units, thus to estimate the linguistic factors that affect these variations on each unit.

We propose in this article a syllable-based duration model based on multi-level context-dependent analysis. This paper is organized as follow : in section 2 our proposed approach is presented ; in section 3 speech material and evaluation scheme are presented ; results are discussed in section 4.

A Multi-Level Context-Dependent Model

Let L = {l1, l2, ..., lN } a set of continuous and nonrecursive linguistic units (i.e. phonem, syllable, prosodic group, period, ...), θ = {θ l 1 , θ l 2 , ..., θ l N } the durational features related to the linguistic levels li and F = {F l 1 , F l 2 , ..., F l N } a set of linguistic features that describe linguistic level li.

Accordingly to the General Superpositional Model [START_REF] Mishra | Decomposition of pitch curves in the general superpositional intonation model[END_REF], the duration θ(n) of the minimal linguistic unit l1 could be written as :

θ(n) = ⊕ l∈L ⊕ k∈l θ l,k (n)
For each linguistic level li, a model typology of θ l i conditionally to F l i is estimated using the Minimum Description Length (MDL) criterion with normal distributions [START_REF] Dell | The use of context in large vocabulary speech recognition[END_REF]. Then given a sequence of linguistic features F l i , the optimal sequence of observations Θl i is estimated as the one that maximizes the conditional probability of observations Θ l i conditionnaly to F l i :

Θl i = argmax Θ l i P (Θ l i |F l i )
which is according to the independance hypothesis the sequence of means :

Θl i == [µ l i ,1 ; ...; µ l i ,K ]
where :

µ l i ,k = argmax θ l i ,k P (θ l i ,k |F l i ,k )
In contrast to models based on modeling durational features on a single linguistic unit (phoneme, syllable), the proposed approach shows several advantages : 1) distinguishing several linguistic units in the representation of durational features variations enables to explicite the superposition of prosodic forms jointly observed on a given unit, 2) each prosodic level (speech rate, duration syllabic residual, ...) can be modeled and controlled independently from each other, 3) estimate the set of linguistic parameters affecting each linguistic unit independently.

Evaluation

Material

The proposed model has been evaluated on two read speech +/-neutral described in table 1. All analysis were conducted within the ircamCorpusTools [START_REF] Veaux | Ircamcorpustools : an extensible plateform for speech corpora exploitation[END_REF] framework. Corpora were segmented using ircamAlign [START_REF] Lanchantin | Automatic phoneme segmentation with relaxed textual constraints[END_REF], syllabified on each interpausal group with Liaphon [START_REF] Béchet | Liaphon : un système complet de phonétisation de textes[END_REF] and syllable prominences have been automatically estimated using ir-camProm [START_REF] Obin | A syllablebased prominence model based on discriminant analysis and context-dependency[END_REF].

Model's parameters

Prosodic features and Linguistic units

The proposed model in this study is a syllable-based duration model. In this experiment durational features were chosen as being speech rate on interpausal group and residual syllable duration.

Linguistic factors

In this experiment low-levels linguistic parameters were used, such as locational features (position of a given unit within higher level units), weight features (number of observations of a given linguistic unit within higher level units) and phonological features (syllabic structure and prominence). The used linguistic parameters on each linguistic unit are presented in table 2.

Syllable level Phonological features

• {preceding, current, succeeding} syllable type, • {preceding, current, succeeding} syllable structure sequence (onset, nucleus, coda),

• phone label of {preceding, current, succeeding} syllable nucleus,

• phone class in {preceding, current, succeeding} syllable nucleus,

• phone label sequence in {preceding, current, succeeding} syllable,

• phone class sequence in {preceding, current, succeeding} syllable,

• phonem sequence in {preceding, current, succeeding} syllable onset,

• phonem sequence in {preceding, current, succeeding} syllable coda,

• prominence state of {preceding, current, succeeding} syllable.

Locational features

• syllable ordinal position in current prosodic group,

• syllable ordinal position in sentence,

• current prosodic group ordinal position in sentence,

• syllable caterogical position (begin ; middle ; end) in current prosodic group,

• syllable categorical position (begin ; middle ; end) in sentence,

• current prosodic group categorical position (begin ; middle ; end) in sentence,

Weight features

• phone number in {preceding, current, succeeding} syllable,

• phone number in {preceding, current, succeeding} syllable onset,

• phone number in {preceding, current, succeeding} syllable coda,

• syllable number in {preceding, current, succeeding} prosodic group,

• syllable number in {preceding, current, succeeding} sentence,

Interpausal group level Locational features

• current prosodic group ordinal position in sentence,

• current prosodic group categorical position (begin ; middle ; end) in sentence.

Weigth features

• number of syllables in {previous ; current ; next} prosodic group,

• number of syllables in sentence,

• number of prosodic groups in sentence, TAB. 2 -List of the contextual features used for each linguistic level.

Evaluation scheme

We have compared our approach to a baseline syllablebased model in which all described linguistic features were projected on the syllable unit. Speaker-dependent models were estimated on a training corpus of variable size and performances were estimated on a test corpus of fixed size (30mn.) In order to evaluate the performance as a function of the amount of training data, variable size training sets were used according to the following sizes : {1mn ; 2mn ; 5mn ; 30mn ; 1hr ; 2hrs ; 5hrs ; [9hrs]}. Performance measures were chosen as being the relative prediction error and correlation score between observed and predicted syllable durations. Relative error measure removes the influence of observed durations on the error measure (observed syllable durations varie from 14 to 900 ms. for the laboratory speech corpus and from 17 to 910 ms. for the acted speech corpus -pauses included). Correlation score is a performance parameter used to estimate the goodness of the predicted temporal structure. This last parameters appears more adequate for measuring prosodic prediction performance since prosody is more related to temporal variations rather than absolute error measurement.

Results and Discussion

Results are presented in figure 1. An overall comparison of the obtained performances on both corpus shows better performances (relative errors and correlation scores) for the laboratory speech than for the acted speech in all cases. Even if these corpora shows comparable overall relative dispersion (table 1), such difference could be explained by the fact that laboratory speech is more prototypical than acted speech and then more predictable from linguistic features. Our proposed approach improves performance in terms of correlation score in almost all cases (with exception of the laboratory corpus with a 1mn. training size). Performance improvement is more significant on the acted corpus for which performance gain varies from 6 to 20% compared to the baseline model. However the proposed approach shows a little degradation in relative errors performance (worse for the laboratory speech and similar for the acted speech corpus)

Gain in correlation score could be explained by several properties of our proposed model : 1) modeling durational features on different linguistic units enables to explicite the observed acoustic forms on each level separately, thus removing from the syllabic level variations due to speech rate variations. This leads to a reduction of the dispersion of residual syllable durations.

2) the correlation score is little sensitive to offset errors caused by speech rate prediction errors. In contrast, the performance degradation in terms of relative error of prediction is due both to the fact that speech rate prediction errors cause a prediction bias on all syllables of the considered interpausal group, but also by the lack of linguistic features that can explain such speech rate variations.

If we look at the evolution of the performances regardless of the model used, we see a clear asymptotic behavior and the gain in performance is no more significant after 1or 2 hours training size. However, this asymptotic behavior in terms of performance is observed jointly to a significant increase in terms of complexity of the models (table 3). On the one hand this means that if low-level linguistic features could be used to estimate robust models with a relatively small amount of data, these remain insufficient to explain more accurately the observed variations. The increase in complexity jointly observed could be explained by the fact that increasement of the number of observations increases at the same time the number of competitive linguistic parameters at each step of the estimattion of the typology of the model. This tends to increase errors on the prioritization of these parameters. Then a larger number of linguistic parameters is needed to obtain similar results.

In terms of complexity (table 3), the proposed model shows a significant reduction in the complexity of the models obtained when a sufficient amount of data is available to learn robust models. (small reduction observed for the laboratory speech from 5hrs. and a significant reduction for the acted speech from 1hrs.). This is mainly due to the fact that the baseline model needs more linguistic parameters to explain variations in speech rate, especially when the model is mixed by trying to jointly explain syllable duration and speech rate variations.

Conclusion

A multi-level approach of prosodic parameters variations applied to durational features has been presented. This approach explicitly represents the different observation levels of prosodic forms and estimate the linguistic parameters that affects these forms separately on each linguistic unit. Compared to a baseline syllable duration model, the proposed approach improves the performance of the predicted temporal structure (correlation score) and reduces the complexity of the models while providing comparable performance in terms of relative prediction error .

In further works, we will focus on two aspects of prosody modeling. At the acoustic level, we will explicitely model the prosodic temporal structure by introducing dependencies between acoustic observations. At the symbolic level, we will extract high-level linguistic features (morpho-syntactic, syntactic, semantic, discursive) from text that will be used to improve the estimation of speech rate variations over intermediate linguistic units.

1 -

 1 Description of the used corpora.

  FIG. 1 -Performance as a function of the training set size. Top : laboratory speech. Bottom : acted speech. Baseline model is represented with dash-lines. Proposed model is represented with plain-lines. Relative prediction errors are plotted with triangles and correlation score with squares.

		laboratory	acted
		speech		speech
		syllable multi-level	syllable multi-level
	1 mn.	16	29	19	24
	2mn.	31	36	28	35
	10mn.	39	48	41	30
	30mn.	103	126	95	89
	1h.	151	167	143	116
	2h.	236	250	230	192
	5h.	488	463	425	331
	9h.	766	722	NA	NA
	TAB. 3 -Models complexity as a function of the training set
	size.

Acknowledgements

This study was supported by : • ANR Rhapsodie 07 Corp-030-01 ; reference prosody corpus of spoken French ; French National Agency of research (ANR) ; 2008-2012., • Programmes Exploratoires Pluridisciplinaires (PEPS), CNRS/ST2I, 2008-2010.


Christophe Louembet
email: christophe.louembet@laas.fr

Franck Cazaurang

Ali Zolghadri

Motion planning for flat systems using positive B-splines: An LMI approach

Keywords: Motion planning, Nonlinear differentially flat systems, B-splines positivity, LMI

In this paper, the motion planning problem is studied for nonlinear differentially flat systems using B-splines parametrization of the flat output history. In order to satisfy the constraints continuously in time, the motion planning problem is transformed into a B-splines positivity problem. The latter problem is formulated as a convex semidefinite programming problem by means of a non-negative piecewise polynomial functions description based on sum of squares decomposition. The contribution of the paper is thus a one-step design procedure for motion planning that satisfies constraints continuously in time where usual Bspline and collocation techniques need post-analysis. Finally, an example of flexible link manipulator motion is presented to illustrate the overall approach.

Introduction

In this paper, the motion planning problem is considered for a subset of nonlinear systems referred to as differentially flat systems. For such systems, introduced first in [START_REF] Fliess | On differentially flat nonliear systems[END_REF], there exists a minimal set of particular outputs (the so-called flat outputs) that characterize all the state space motions and the corresponding input history. Such systems arise in several relevant domains of engineering control (see [START_REF] Martin | Flat systems, equivalence and trajectography[END_REF] and the references therein). Contrary to classical methods for solving optimal control problem [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF][START_REF] Bryson | Dynamic optimization[END_REF], solving the motion planning problem through flatness avoids integrating the differential equations of the dynamics. In fact, since the flat output is the state of a trivial system [START_REF] Fliess | a Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF], the above optimal control problem boils down to find the best flat output motion that lies into the subspace defined by the path constraints expressed in terms of flat outputs and their derivatives and passes through a given set of points. In this context, a classical and tractable methodology relies on Bsplines based collocation [START_REF] Milam | A new computational approach to real-time trajectory generation for constrained mechanical systems[END_REF][START_REF] Faiz | Differentially flat systems with inequality constraints: an approach to real-time feasible trajectory generation[END_REF][START_REF] Louembet | Path planning for satellite slew maneuvers: a combined flatness and collocation based approach[END_REF]. However, this technique involves time sampling: there is no guarantee on constraints satisfaction between collocation points. This may lead to critical issues that need to be detected by an appropriate post-analysis. Thus, the B-splines col-location requires an interactive procedure between the trajectory synthesis and a specific post-analysis. In this paper, our goal is to design flat system trajectories that guarantee constraints continuously in time and without a post-analysis. Henrion and Lasserre tackled this problem in the case of linear systems in [START_REF] Henrion | LMIs for constrained polynomial interpolation with application in trajectory planning[END_REF]. They proved that motion planning under constraints can be recast as the inclusion of a univariate polynomial in a linear semi-algebraic subset. By using results on positive polynomials [START_REF] Nesterov | High Performance Optimization[END_REF], the late problem turns out to be a linear matrix inequalities (LMI) problem. The main result of this paper is obtained in two steps. First, We extend the results from [START_REF] Nesterov | High Performance Optimization[END_REF] to positive piecewise polynomials. This result is then used, in conjunction with Bsplines parametrization, to provide a new motion planning methodology that allows to design a trajectory fulfilling the constraints continuously in time. In this works, we choose to use the B-splines basis to represent the squared piecewise polynomials rather than the monomial representation and to develop a B-spline positivity concept instead of classical SOS methods used in [START_REF] Henrion | LMIs for constrained polynomial interpolation with application in trajectory planning[END_REF] for the following reasons: (i) B-splines provide a minimal basis for piecewise polynomials, (ii) B-splines basis contains explicitely information about continuity class at breakpoints contrary to others representations. In section 2, the optimal path planning problem for flat systems is described. Then, our contribution is detailed in two steps. First, in section 3.1 and 3.2, the results on positive polynomials [START_REF] Nesterov | High Performance Optimization[END_REF] are extended to positive piecewise polynomials by characterizing them over the cone of positive semidefinite matrices using the sum of squares formalism. Subsequently, in section 3.3, the constrained B-splines optimization is formulated as a convex optimization problem over linear matrix inequalities (LMI) for which efficient programming (SDP) solvers are available. In section 4, an example illustrates the methodology.

Problem Statement

Differential flatness

Before adressing the motion planning problem for flat systems, let us briefly recall the differential flatness concept. Differential flatness, or shortly flatness, has been introduced by Fliess et al. [START_REF] Fliess | On differentially flat nonliear systems[END_REF] in 1992. Consider a nonlinear system: ẋ = f (x, u), where x is the n-component state vector and u the m-component control with m ≤ n.

Definition 1 The nonlinear system is differentially flat if there exists an m-dimensional vector z, so-called flat outputs, such that:

z(t) = Φ x(t), u(t), u(t), . . . , u (η) (t) , and (1)
x(t) = Ψ x z(t), ż(t), . . . , z (β-1) (t) , u(t) = Ψ u z(t), ż(t), . . . , z (β) (t) , (2)
where Ψ x and Ψ u are smooth functions, z

i (t) denoting the k th order time derivative of the i th component of z(t) and η, β ∈ N.

Flatness is also defined as a Lie-Bäcklund equivalence between a nonlinear system and a trivial system in [START_REF] Fliess | a Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF]. As z represents the state of the equivalent trivial system, the m-components of z are differentially independent. Indeed, a z-space of dimension n z is considered with the coordinates z = {z, ż, z, . . . , z (p) } with p ∈ N where any curve is equivalent to a trajectory of the corresponding nonlinear system.

Optimal path planning for the flat system

Using the flatness properties in line with the works [START_REF] Milam | A new computational approach to real-time trajectory generation for constrained mechanical systems[END_REF][START_REF] Louembet | Path planning for satellite slew maneuvers: a combined flatness and collocation based approach[END_REF], the generation of an optimal and constrained trajectory, t → (x * (t), u * (t)), for nonlinear flat system is equivalent to the following integration-free optimization problem:

min z J(z) subject to: z(t i ) = z i , z(t f ) = z f , z(t) ∈ S, (3)
where z is flat space coordinate, J is an objective function, and the subset S, so-called the feasible region, is such that S = {z| γ(z(t)) ≤ 0}, where the inequality γ(z(t)) expresses the path and actuators constraints. We assume that S is a polytopic subset. If it is not the case, it will be replaced by a polytopic inner convex approximation (see [START_REF] Faiz | Differentially flat systems with inequality constraints: an approach to real-time feasible trajectory generation[END_REF] for details on the approximation technique and [START_REF] Louembet | Path planning for satellite slew maneuvers: a combined flatness and collocation based approach[END_REF] for the study of approximation conservatism). Let us describe the flat output t → z(t) components by univariate piecewise polynomials. Piecewise polynomials have been historically used in solving optimal control problem due to their ability to handle multiple convex or non-convex constraints from especially dynamics transcription [START_REF] Ascher | Numerical solution of boundary value problems in ordinary differential equations[END_REF][START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF][START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF]. Indeed, z(t) admits the k th order B-splines basis {B j,k } (see [5, Chap. IX]) as,

z i (t) = nB j=1 C i,j .B j,k (t), i = 1, . . . , m z (p) i (t) = nB j=1 C i,j .B (p) j,k (t), i = 1, . . . , m. (4)
with C = (C 1,1 , . . . , C 1,nB , . . . , C m,1 , . . . , C m,nB ) T as the coordinates of z(t) in the B-splines basis B(t). In the sequel, the control points C become the decision variables of the flat optimal control problem (3). It comes the parametrized optimal flat path planning:

min C J(z(C))
subject to:

       z(t i , C) = z i , z(t f , C) = z f , z(C) ∈ S.
.

The purpose of section 3 is to set constraint z(C) ∈ S as an inclusion problem of z(t) trajectory within the intersection of several half-spaces. In fact, it will be shown that positioning the trajectory z(t) in a half-space is equivalent to evaluate the sign of the piecewise polynomial gap function, κ(t), between z(t) and the hyperplane boundary. In subsections 3.1 and 3.2, the positivity of a piecewise polynomial function in terms of its coordinates in a B-splines basis v(t) is defined through the sums of squares formalism. Then, in subsection 3.3, we formulate the κ(t)-positivity in terms of the coefficients C using the LMI condition developped in section 3.2 in order to solve [START_REF] De Boor | A pratical guide to spline[END_REF].

3 Motion planning as B-splines positivity problem

B-splines representation of squared piecewise polynomials

Inspired by the seminal work of Nesterov [START_REF] Nesterov | High Performance Optimization[END_REF], let us describe here the decomposition in sums of squares of the non-negative piecewise polynomial functions (PP). First, the k v -order PP functions will be described as a weighted sum of squared k w -order PP functions. Let P kw ,ξ,νw be the linear subspace defined by the collection of the k w -order piecewise polynomial functions of P kw,ξ defined with breakpoints ξ = {ξ

F (w 2 ) = {P (t) = N i=1 d i τ 2 i (t)|τ i (t) ∈ P kw,ξ,νw , d i ∈ R}.
By construction, the description of F (w 2 ) depends only on the properties of the set of the squared basis functions:

w 2 = {w i (t)w j (t), i, j = 1, . . . , n w }. (6)
It can be proved that F (w 2 ) is included P kv ,ξ,νv with k v = 2k w -1, ν v = ν w and ξ as defined above. Using the Curry-Schoenberg theorem [START_REF] Curry | On pólya frequency functions: The fundamental spline functions and their limits[END_REF] it follows that F (w 2 ) admits a k v -order B-splines basis v(t) such that:

for all P (t) ∈ F(w 2 ), P (t) = nv i=1 µ i v i (t). ( 7
)
Let us now introduce the operators Λ and Λ * that determine the coordinates of a element of F (w 2 ) in the B-splines basis v(t). These definitions will be useful in section 3.2.

Definition 2 (Operator Λ) Let Λ a linear operator such that P nv kv,ξ,νv → P nw×nw kv ,ξ,νv : v(t) → Λ(v) ≡ w(t)w(t) T be defined by:

Λ(v) = [w i .w j ], i, j = 1, . . . , n w =      λ T 1,1 v(t) . . . λ T 1,nw v(t) . . . . . . . . . λ T nw ,1 v(t) . . . λ T nw ,nw v(t)      , λ i,j ∈ R n (8)
In the following, •, • denotes the bilinear form:

R α×β × R α×β → R, A, B → A, B = α i=1 β j=1 a i,j b i,j . (9)
Definition 3 (Operator Λ * ) The linear dual operator

Λ * : R nw×nw → R nv , Y → Λ * (Y ) is defined such that: Y, Λ(v(t)) ≡ Λ * (Y ), v(t) , ∀Y ∈ R nw×nw , v(t) ∈ R nv . ( 10
) Recalling that Y, Λ(v(t)) = nw i,j=1 Y i,j nv l=1 λ i,j,l v l (t),
the operator Λ * is defined as follow:

Λ * (Y ) = (Λ * 1 (Y i,j , λ i,j ), . . . , Λ * nv (Y i,j , λ i,j )) T , ( 11
)
where Λ * i are scalar functions linear in Y i,j .

The latter definition describes the coordinates in the Bsplines basis v(t) of a F (w 2 ) PP function. In fact, the matrix Y represents the weighting coefficients on w 2 (t) and Λ * (Y ) denotes the coordinates in the B-splines basis v(t).

Piecewise polynomials positivity

A piecewise polynomial has two representations: the Bsplines and the sums of squares. The sums of squared PP function representation is convenient since its positiveness only depends on the semi definite positiveness of the weighting matrix Y . Then, through the operators defined above, we are able to describe the set of the coefficients µ on basis v(t) that define positive PP function. In fact, this set is shown to be a linear image of the cone of the positive semidefinite matrices.

Theorem 4 Let µ be an element of the closed, pointed and convex cone K defined by:

K = {µ ∈ R nv : µ = Λ * (Y ), Y 0}. ( 12
)
where Y 0 define the semidefinite positiveness of Y .eigenvalues non negative. Each element µ of K describes a positive semidefinite polynomial on the basis v(t) so that

P (t) = nv i=1 µ i v i (t) ≥ 0. ( 13
)
Proof First, Let P = µ T v with µ = Λ * (Y ) and Y 0. Then, P can be written as:

P (t) = Λ * (Y ), v(t) = Y, Λ(v(t)) = Y, w(t)w(t) T . (14)
Considering Iw(t), w(t) = I, w(t)w(t) T with I the identity matrix, the previous relation becomes:

P (t) = Y w(t), w(t) = (Y w(t)) T w(t) = w(t) T Y T w(t) ≥ 0. (15)
Second, if P (t) ≥ 0, then there exits a set of vectors γ i ∈ R nw , i = 1, . . . N such that:

P (t) = N i=1 γ i , w(t) 2 = N i=1 γ i γ T i , w(t)w(t) T = N i=1 γ i γ T i , Λ(v(t)) = Λ * N i=1 γ i γ T i , v(t) . (16)
Thus, taking Y = N i=1 γ i γ T i 0 and µ = Λ * (Y ) allows us to establish [START_REF] Louembet | Path planning for satellite slew maneuvers: a combined flatness and collocation based approach[END_REF].

Remark For the sake of conciseness only the odd order case, k v = 2k w -1, has been exposed here. But the even order case k v = 2k w is very similar to the odd case: from the beginning, one need to replace the basis function w i by wi = τ (t)w i with τ (t) = tξ 1 + 1 strictly positive on [ξ 1 , ξ l+1 ]. In this way, operator Λ * can be redefined in order to describe the non negative even order PP function:

P (t) = nv i=1 µ i v i (t) = τ (t) N i=1 τ 2 i (t) ≥ 0

Path planning problem as a positivity problem

The next step of the contribution consists of setting problem [START_REF] De Boor | A pratical guide to spline[END_REF] as an SDP problem. This result, mainly based on theorem 4, is the description of the PP trajectory inclusion into a polytope as a B-spline positivity problem and consequently as an LMI problem. Let O z be the finite dimensional flat output space with the following coordinates:

z = (z 1 , . . . , z m , ż1 , . . . , żm , z (r) 1 , . . . , z (r) m ) T .
Recall that the flat trajectories [t 0 , t f ] → R n z , t → z(t) are parametrized on k-order B-splines basis {B k } (see equation ( 4)).

Let the feasible region S be an intersection of n c halfspaces of O z and H i be the i th half-space described by its Cartesian coordinates:

H i = {z ∈ R n z | a T i z ≤ b i }, (17)
where a i ∈ R n z and b i ∈ R with i = 1, . . . , n c . We note thatz(t) belongs to the half-space H i if and only if

a T i z(t) ≤ b i (18)
Theorem 5 Solving the path planning problem defined by [START_REF] De Boor | A pratical guide to spline[END_REF], is equivalent to solving the following SDP problem:

min C J(z(C))
subject to:

   α i C -b i = Λ * (Y i ) Y i 0 ΘC = θ , ∀i = 1, . . . , n c . (19
) with the objective function assumed to be linear in z and so in

C = (C 1,1 , . . . , C 1,n , C 2,1 , . . . , C 2,n , C 3,1 , . . . , C 3,n ), the control points.
Λ * is the dual operator defined by [START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF]. α i ∈ R nv×NC are linear matrix functions of a i , with a i and b i associated to the i th half-space H i (cf. equation [START_REF] Toh | Sdpt3 -a matlab software package for semidefinite programming[END_REF]). The equality constraint ΘC = θ represents the initial and final conditions. The proof is detailled in appendix A.

Example

To illustrate the developed methodology, the motion planning problem of a flexible link manipulator system [START_REF] Faiz | Real-time and optimal trajectory generation for nonlinear systems[END_REF] will be solved, subject to path constraints and a cost criterion. The dynamics is described by the following equations:

I 1 q1 + M gL sin q 1 + k(q 1 -q 2 ) = 0, (20) I 2 q2 -k(q 1 -q 2 ) = u, (21)
where x = [q 1 , q1 , q 2 , q2 ] T is the state vector and w the input. For this example, the following values (in S.I. units) are considered:

L = 0.1, M = 0.1, I 1 = 1, I 2 = 1, k = 1, g = 9.
8. The goal is to steer the system from initial state q(t 0 ) = [0.8, 0, 0.67, 0] to final state q(t f ) = [-0.8, 0, -0.67, 0] at final time t f = 5.35s. The operating constraints are

-π 3 ≤ q 1 ≤ π 3 , -π 4 ≤ q 2 ≤ π 4 , -π 12 ≤ q 2 -q 1 ≤ π 12 .
A flat output for this system has been described in [START_REF] Faiz | Real-time and optimal trajectory generation for nonlinear systems[END_REF]: z = q 1 . Thus, the flat parametrization of the state x is described by:

x = φ(z, ż, z, z (3) ) = z, ż, I1 k z + MgL k sin z + z, I1 k z (3) + MgL k ż cos z + ż T (22)
Using the flat parametrization ( 22), the path planning problem is translated in the following flat optimal control problem described by: min z J(z)

subject to:

       z(t i ) = z i , z(t f ) = z f -π 3 ≤ z ≤ π 3 , -π 4 ≤ I1 k z + MgL k sin z + z ≤ π 4 , -π 16 ≤ I1 k z + MgL k sin z ≤ π 16 . (23)
Note that the operating constraints depend only on z and z. Thus, the flat output space O z has the following coordinates z = {z, z}. We note that S is non convex. However as we said earlier, it can be replaced by a polytopic inner approximation S p calculated with Faiz' techniques [START_REF] Faiz | Differentially flat systems with inequality constraints: an approach to real-time feasible trajectory generation[END_REF]. The half spaces representation of S p = {z|Hz ≤ K} is given by the matrices

z(t i , C) = z i , z(t f , C) = z f , z(t) ∈ S p . ( 24
)
The trajectory t → z(t) is a 5 th order PP function defined on the sequence of equidistant knots ξ = {ξ 1 , . . . , ξ 10 }. Indeed z(t) as element of P k,ξ,ν admits B(t) as B-splines basis. The continuity parameter

         α 1 C -b 1 = Λ * (Y 1 ), Y 1 0, . . . α 4 C -b 4 = Λ * (Y 4 ), Y 4 0, ΘC = θ (25)
where C ∈ R 3n , α i ∈ R nv×3n . For the sake if simplicity, the cost is linear here, J(C) = z(t f , C). But one can choose any convex function J LMI (C) that is LMI representable. In this case, an additional decision variable τ is minimized instead of J LMI (C) such that τ ≥ J LMI (C). Problem (25) is solved using SDPT3 [START_REF] Toh | Sdpt3 -a matlab software package for semidefinite programming[END_REF] in 0.29 s with a 2.26GHz Intel Core 2 processor. The obtained trajectory is given in figure 1. Figure 1 shows that trajectory z is clearly included in S p as expected.

For sake of comparison, we solve problem (5) by means of flatness and collocation methods described in [START_REF] Louembet | Path planning for satellite slew maneuvers: a combined flatness and collocation based approach[END_REF][START_REF] Milam | A new computational approach to real-time trajectory generation for constrained mechanical systems[END_REF][START_REF] Faiz | Differentially flat systems with inequality constraints: an approach to real-time feasible trajectory generation[END_REF].

Recall that the constraints are checked in a finite number of time called collocation points. The method is programmed and solved with the linear solver Matlab linprog considering 10 equidistant collocation points.

Although an admissible solution for the collocation problem is obtained quickly (0.06 s), the trajectory does not respect the constraints between collocation points ♯2 and ♯3 then between ♯8 and ♯9 (see figure 1). Then, with the above collocation methods, an iterative process is needed to re-distribute the sequence or increase the number the collocation points, conversely to the one-run method exposed in this paper.

Concluding remarks

This paper considers motion planning for differentially flat systems. As opposed to most works on direct methods for optimal control problem reported in the literature, the developed methodology provides a new framework for satisfying constraints continuously in time. In this methodology, a crucial point is the approximation of non convex feasible regions since the accuracy of the results depends on its conservatism. Despite the method used above, the computation of convex (especially LMI representable) inner approximations remains an open research area.

A Proof of theorem 5

To apply the positivity theorem, the inequality a T i z(t) ≥ b i must be expressed in a B-splines basis (we flipped the inequality (18) without loss of generality). Thus, by using (4), equation ( 18) is equivalent to

  086 0.731 0.086 -0.731 -0.996 0.683 0.996 -0.683 T and K = [0.192 0.527 0.192 0.527]. After B-splines parametrization, the motion planning problem can be represented as min C

Fig. 1 .

 1 Fig. 1. Trajectories z(t) obtained by SDP (plain line) and by collocation (dotted line), the collocation points are the crosses, thin lines give actual constraints, and the polytopic approximation is the shaded zone.

  ,1 B j,k (t) + • • • + a i,(r-1)m+1 B (r) j,k (t))C 1,j . . . + (a i,2 B j,k (t) + • • • + a i,(r-1)m+2 B (r) j,k (t))C 2,j + . . . + (a i,m B j,k (t) + • • • + a i,rm B (r) j,k (t))C m,j ≥ b i . (A.1)In inequality (A.1), the PP function is composed of a P k,ξ,ν piecewise polynomial and its r first derivatives. Considering that for a B-spline B j,k ∈ P k,ξ,ν , one has Ḃj,k ∈ P k-1,ξ,ν⊖1 , . . . , B(r)j,k ∈ P k-r,ξ,ν⊖r where

  1 , • • • , ξ l+1 } whose first ν wi derivatives are continuous at ξ i (i.e. that are C νw i at ξ i with ν w = {ν w1 , • • • , ν w l+1 }); see [5, Chap. VIII] for more detailed definitions. Then, by virtue of the Curry-Schoenberg theorem [4], P kw,ξ,νw

admits a B-splines basis w(t) = (B 1,kw , . . . , B nw ,kw ) T = (w 1 (t), . . . , w nw (t)) T . The function subspace of sums of squared functions from P kw,ξ,νw is described with N ∈ N by:

Table 1 B

 1 -splines basis parameters vector ν are given in table1. The dimension of the basis B(t) is n = 14. In order to define Λ * , we need to characterize the gap function κ(t). Since the higher derivation order involved in O z is two, κ(t) belongs to P k,ξ,ν⊖2 (see appendix A) and thus it admits a B-splines basis v(t) of dimension n v = 32. Then, we calculate the corresponding basis w(t). Its dimension is n w = 12. The operator Λ

* is deduced from v(t), w(t) and definition 3. Coefficient matrices α i of problem (

19

) are then calculated for each half-space of polytope S p (see (A.3)). The problem equivalent to (24) is finally set using Yalmip

[START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in MATLAB[END_REF]

: min C J(C) subject to:

the operator ⊖ is defined by ν ⊖ r = (max{ν 1r, 0}, . . . , max{ν l+1r, 0}), with r ≤ ν i for i = 2, . . . , l. Then the sum (A.1), representing the gap, belongs to P k,ξ,ν⊖r and, by virtue of the Curry-Schoenberg theorem [START_REF] Curry | On pólya frequency functions: The fundamental spline functions and their limits[END_REF], admits the B-splines representation basis v(t) such that (A.1) becomes:

with α p,j ∈ R nv p = 1, . . . , m and j = 1, . . . , n. Vectors α l,i,j are identified using the following equalities system:

3) with l = 1, . . . , m and j = 1, . . . , n. E p (f ) = t f 0 x p f (x) dx denotes the p th order moment of the function f . The index p is chosen such that equation (A.3) leads to a square linear matrix equality to obtain α l,i,j . Thus, inequality (18) is equivalent to the following positivity problem:

where

Then, determining the operators Λ and Λ * is needed to recast the positivity problem (A.4) into an LMI problem by using theorem 4. These operators are built with a basis u(t) satisfying the following inequality: ν ⊖ r < k+1 2 if k is odd, or ν ⊖ r < k 2 if k is even. So, theorem 4 gives conditions on the κ coefficients so that inequality (A) holds:

(A.6) Finally, the inclusion of a trajectory t → z(t) to the intersection of n c half-spaces is written as the conjunction of the n c membership problem defined in the theorem 5 i.e.