A Leygue

C Bailly
email: bailly@poly.ucl.ac.be

R Keunings
email: roland.keunings@inma.ucl.ac.be

A differential tube-based model for predicting the linear viscoelastic moduli of polydisperse entangled linear polymers

Keywords: Entangled linear polymers, Linear visco-elasticity, Polymer mixtures

may come

A differential tube-based model for predicting the linear viscoelastic moduli of polydisperse entangled linear polymers

Adrien Leygue, Christian Bailly, Roland

Introduction

Tube theories addressing the quantitative prediction of the linear viscoelastic properties of entangled systems based on their microstructure have now reached a high level of maturity. Following the pioneering work of de Gennes [START_REF] De | Reptation of a polymer chain in the presence of fixed obstacles[END_REF] on reptation, numerous tube models are now capable of a quantitative description of the linear viscoelastic properties of linear entangled polymers (see e.g. [START_REF] Pattamaprom | Quantitative predictions of linear viscoelastic rheological properties of entangled polymers[END_REF][START_REF] Léonardi | Rheological models based on the double reptation mixing rule: the effects of a polydisperse environment[END_REF][START_REF] Van Ruymbeke | Evaluation of reptation models for predicting the linear viscoelastic properties of entangled linear polymers[END_REF][START_REF] Likhtman | Quantative theory for linear dynamics of linear entangled polymers[END_REF][START_REF] Milner | Reptation and contour-length fluctuations in melts of linear polymers[END_REF]). The success of these theories is such that the focus is now on the prediction of the linear viscoelastic response of systems with a much more complex architecture such as mixtures of star and linear polymers [START_REF] Milner | Dynamic dilution, constraint-release, and star-linear blends[END_REF][START_REF] Lee | Stress relaxation of star/linear polymer blends[END_REF][START_REF] Van Ruymbeke | Prediction of linear viscoelastic properties for polydisperse mixtures of entangled star and linear polymers: modified tube-based model and comparison with experimental results[END_REF]. For systems involving more complex architectures such as long chain branching, a hierarchical algorithm has been proposed by Larson [START_REF] Larson | Combinatorial rheology of branched polymer melts[END_REF]. Park et al. [START_REF] Park | A hierarchical algorithm for predicting the linear viscoelastic properties of polymer melts with long-chain branching[END_REF] have modified it to incorporate recent theoretical developments. Despite their focus on the linear rheology of entangled systems, a striking fact about all those theories is that they are mathematically non-linear.

In order to predict the relaxation modulus of an entangled polymer, one must describe how the polymeric chains escape from their constraining tubes and how this relaxation process affects the tubes themselves. For linear polymers, the escape of a chain can occur through reptation or contour-length fluctuations [START_REF] Doi | The Theory of Polymer Dynamics[END_REF]. Although they are often addressed separately, these two phenomena can be viewed as different modes of a Rouse chain trapped in a tube potential. For example, the stochastic description of linear chains [START_REF] Likhtman | Quantative theory for linear dynamics of linear entangled polymers[END_REF][START_REF] Hua | Segment connectivity, chain-length breathing, segmental stretch, and constraint release in reptation models. i. theory and single-step strain predictions[END_REF][START_REF] Masubuchi | Brownian simulations of a network of reptating primitive chains[END_REF] does not consider reptation and fluctuations separately as fluctuations correspond to the breathing modes of a one-dimensional Rouse chain escaping from a tube. Different authors have proposed coupled descriptions of reptation and contour-length fluctuations by either considering a position-dependent [START_REF] Cloizeaux | Relaxation of entangled polymers in melts[END_REF][START_REF] Graham | Microscopic theory of linear entangled polymer chains under rapid deformation including chain stretch and convective constraint release[END_REF] or a time-dependent [START_REF] Cloizeaux | Relaxation and viscosity anomaly of melts made of long entangled polymers[END_REF] diffusion process along the primitive path.

Thermal Constraint Release [START_REF] Marrucci | Relaxation by reptation and tube enlargement: A model for polydisperse polymers[END_REF] theories address the effects, on the tube, of the relaxation of the surrounding chains. In a sense, constraint release is a closure to the tube potential. Viovy et al. [START_REF] Viovy | Constraint release in polymer melts: tube reorganization versus tube dilation[END_REF] have proposed a rigorous description of constraint release and have shown when Rouse-like tube motion induced by thermal constraint release can be approximated by a dilation of the tube. The double reptation picture [START_REF] Tsenoglou | Viscoelasticity of binary polymer blends[END_REF][START_REF] Cloizeaux | Double reptation versus simple reptation in polymer melts[END_REF] provides a simple and successful alternative to complex constraint re-lease theories. Double reptation actually is a fair approximation of Rouse tube motion. Likthman and McLeish have recently shown that generalized double reptation is essentially correct down to the Rouse time of the chain [START_REF] Likhtman | Quantative theory for linear dynamics of linear entangled polymers[END_REF]. For polydisperse systems, double reptation induces a non-linear mixing rule which has actually been recovered in the stochastic reptation model proposed by Hua et al. [START_REF] Hua | Segment connectivity, chain-length breathing, segmental stretch, and constraint release in reptation models. i. theory and single-step strain predictions[END_REF]. A limitation of double reptation appears for polydisperse systems where different chains relax on well separated timescales. In that case, the random pairwise interactions picture of double reptation breaks down, as the fast dynamics of the short chains would induce a very fast relaxation of the slower chains, which is not observed experimentally.

In this paper, we present a simple tube model for mixtures of linear polymers. This theory, valid in the rheologically linear regime only, proposes a coupled description of reptation, contour-length fluctuations and constraint release. Reptation and contour-length fluctuations are accounted for through a position dependent diffusion mechanism proposed by Graham et al. [START_REF] Graham | Microscopic theory of linear entangled polymer chains under rapid deformation including chain stretch and convective constraint release[END_REF], while thermal constraint release is introduced using an extension of a new differential formulation of double reptation [START_REF] Leygue | A differential formulation of thermal constraint release for entangled polymers[END_REF]. Our theory allows one to compute directly the tube survival probability along the primitive path. Furthermore, the extension to the polydisperse case yields a perfectly linear mixing rule. Then we demonstrate the predictive capabilities of our approach on a large set of polystyrene samples. Finally, we show how it is possible to build a mathematically linear approximation of our theory that retains both the structure of the model and the physical meaning of the variables. The resulting model is therefore mathematically linear (linear equations and mixing rule) and is able to predict quantitatively the linear viscoelastic properties of entangled linear polymers. The next step is of course to build a constitutive equation for the non-linear rheological regime whose linearization would be the linear model presented in this work [START_REF] Leygue | A tube-based constitutive equation for polydisperse entangled linear polymers[END_REF].

Tube theory for a monodisperse system

In this section we present a tube theory for entangled linear systems. We first consider the monodisperse case and then extend the theory to account for polydispersity.

Let s (-1 ≤ s ≤ 1) be a parametric coordinate along the primitive path of a polymer chain in a monodisperse environment. Let us define P f γ (t, s) the probability for a tube segment with position s to survive between the arbitrary initial time 0 and time t. The superscript f denotes that we do account for contour-length fluctuations, while the parameter γ is representative of the effectiveness of constraint release. Two limiting cases are P f 0 (t, s) and P γ (t, s). The quantity P f 0 (t, s) is representative of a chain relaxing in a fixed network of entanglements where no constraint release occurs, while P γ (t, s) properly describes a very long chain where the length-fluctuations become negligible. Reptation theory [START_REF] De | Reptation of a polymer chain in the presence of fixed obstacles[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF] suggests that the time evolution of P f 0 (t, s) is governed by a diffusion operator along the coordinate s. This operator represents the diffusive motion of the chain inside its tube that yields the destruction/renewal of a tube segment that reaches the chain ends. Following Des Cloizeaux [START_REF] Cloizeaux | Relaxation of entangled polymers in melts[END_REF] and Graham [START_REF] Graham | Microscopic theory of linear entangled polymer chains under rapid deformation including chain stretch and convective constraint release[END_REF], we consider that the effect of contour-length fluctuations can be accounted for through a position-dependent diffusion constant. The simplest form for the evolution of P f 0 (t, s) is therefore [START_REF] Doi | The Theory of Polymer Dynamics[END_REF]:

∂ ∂t P f 0 = ∂ ∂s α f d (s) ∂ ∂s P f 0 , ( 1
)
P f 0 (t, -1) = 0, (2)
P f 0 (t, 1) = 0 for t > 0, (3)
P f 0 (0, s) = 1 for -1 < s < 1,
where Eqs. ( 2) and ( 3) come from the assumption that chain ends are fully relaxed.

As suggested in [START_REF] Graham | Microscopic theory of linear entangled polymer chains under rapid deformation including chain stretch and convective constraint release[END_REF] we consider α f d (s) to be constant along the primitive path, except at the chain ends in two regions of size proportional to √ M/M e , where M is the molecular weight of the polymeric chain and M e is the molecular weight between entanglements. In those regions, where contour-length fluctuations allow for a faster renewal of the tube segments, the diffusion coefficient is proportional to the inverse of the quadratic distance to the chain end. As the molecular weight of the polymer increases, the influence of contour-length fluctuations becomes more and more negligible. In that limiting case, the reptation picture [START_REF] De | Reptation of a polymer chain in the presence of fixed obstacles[END_REF] tells us that the diffusion constant should scale like the inverse of the cube of the molecular weight of the polymer. The full expression for α f d (s) becomes:

α f d (s) = 4 K d π 2 M 3 K f 2 M e M(1 -s) 2 if s > 1 -K f M e M , ( 4
) 4 K d π 2 M 3 K f 2 M e M(1 + s) 2 if s < K f M e M -1 , (5) 4 K d π 2 M 3 otherwise, (6)
where K d is a material parameter and the M -3 dependence of α f d (s) is due to the non-dimensional nature of the spatial parametric coordinate s. The adjustable parameter K f is close to unity and controls the depth of the contour-length fluctuations within the model. We do not claim any universal value for K f and consider it as a parameter that has to be identified for each material. This is similar to the parameter M of Des Cloizeaux [START_REF] Cloizeaux | Relaxation and viscosity anomaly of melts made of long entangled polymers[END_REF] that also controls fluctuations. Taking K f = 0 actually suppresses the fluctuations as the diffusion constant becomes constant along the chain. In that special case, one recovers the simple reptation picture of Doi and Edwards [START_REF] Doi | The Theory of Polymer Dynamics[END_REF]. The reptation time τ d (M) can then be defined as:

τ d (M) = K d M 3 . ( 7
)
At this point, it is worth noticing that the simple scaling we adopt for the width of the fluctuation zone cannot be valid for poorly entangled systems. When M is only a few times M e , the reptation picture tends to break down as fluctuations become the dominant relaxation mechanism.

Finally we account for constraint release using the approach we developed in [START_REF] Leygue | A differential formulation of thermal constraint release for entangled polymers[END_REF]. This is done by adding a local relaxation term to Eq. ( 1) representative of relaxation through constraint release. Following the double reptation picture, we postulate that the rate of constraint release is proportional to the rate of relaxation through reptation and contour-length fluctuations, and that constraint release events occur with the same probability on all unrelaxed segments. The first assumption translates the fact that constraint release is actually driven by the other relaxation mechanisms. The second one states that the probability for a segment to relax through constraint release does not depend on its position along the chain. At a given position along the primitive path, the rate of relaxation through constraint release must therefore be proportional to the local tube survival probability. The resulting differential problem for P f γ (t, s) reads:

∂P f γ ∂t = ∂ ∂s α f d (s) ∂ ∂s P f γ + βP f γ , ( 8
)
P f γ (t, -1) = 0, P f γ (t, 1) = 0 for t > 0, P f γ (0, s) = 1 for -1 < s < 1,
where β(P f γ ) is the instantaneous rate of thermal constraint release. We now extend the expression for β(P f γ ) proposed in [START_REF] Leygue | A differential formulation of thermal constraint release for entangled polymers[END_REF] to account for both reptation and contour-length fluctuations:

β(P f γ ) = γ 1 -1 ∂/∂s(α f d (s)∂/∂s(P f γ )) ds 1 -1 P f γ ds , ( 9
)
where γ is an adjustable parameter of order unity that controls the amount of constraint release in the system. It is important to notice that β is not constant in time and therefore the local relaxation does not obey a simple first-order decay.

The relaxation modulus G(t) is simply defined as being proportional to the average tube survival probability along the chain:

G(t) = G 0 N 1 2 1 -1 P f γ (t, s) ds, (10)
G 0 N being the plateau modulus.

In order to predict the high frequency response, the expression (10) for the relaxation modulus has to be modified to account for the Rouse modes of the chain. In practice, this is simply done by superposing a Rouse relaxation modulus G r (t) on the relaxation modulus coming from reptation. Following [START_REF] Van Ruymbeke | Evaluation of reptation models for predicting the linear viscoelastic properties of entangled linear polymers[END_REF], the expression for the Rouse relaxation modulus for a chain of mass M writes:

G r (t) = G 0 N   ∞ p=Z+1 1 Z exp - p 2 τ r + 1 3 Z p=1 1 Z exp - p 2 τ r   , ( 11
)
where τ r = K r M 2 and K r is an additional linear material parameter. The integer Z is defined as the closest integer to the ratio M/M e . The material parameters K d , K r and M e are not independent as τ d (M) is linked to τ r (M) through the underlying segmental dynamics, characterized by the equilibration time τ e . Following Larson et al. [START_REF] Leal | Definitions of entanglement spacing and time constants in the tube model[END_REF] we have:

K d = 3 τ e M 3 e , ( 12
) K r = τ e M 2 e . ( 13
)
The independent parameters of the model are therefore M e , τ e , G 0 N , K f and γ.

Extension to the polydisperse case

We now extend the theory presented in the previous section to the polydisperse case, without adding any adjustable parameter. The extension procedure is similar to the one we proposed in [START_REF] Leygue | A differential formulation of thermal constraint release for entangled polymers[END_REF], where contour-length fluctuations were neglected. Let us consider a polydisperse system of polymers with N different masses M (i) . Let P f (i) γ (t, s) be the tube survival probability at time t and parametric position s for the chains of mass M (i) . We assume that all chains reptate and fluctuate independently of the polydispersity of the system. The rate of constraint release is however the same for all chains and is computed from the relaxation of all masses. The differential problem governing the evolution of P f (i) γ reads:

∂P f (i) γ ∂t = ∂ ∂s α f (i) d (s) ∂ ∂s P f (i) γ + βP f (i) γ , (14)
P f (i) γ (t, -1) = 0, P f (i) γ (t, 1) = 0 for t > 0,

P f (i) γ (0, s) = 1 for -1 < s < 1.
In view of Eq. ( 9), the rate of constraint release β is consistently defined as:

β = γ i φ (i) 1 -1 ∂/∂s(α f (i) d (s)∂/∂s(P f (i) γ )) ds i φ (i) 1 -1 P f (i) γ ds , ( 15
)
where φ (i) is the volume fraction of species (i) and γ keeps the same meaning as in the monodisperse case. The diffusion coefficient α f (i) d (s) is computed from Eq. ( 4), where M (i) is substituted for M. Similarly to Eq. ( 10), we define the relaxation modulus as the volume average over the different species of the tube survival probability:

G(t) = G 0 N N i=1 φ (i) 2 1 -1 P f (i) γ (t, s) ds. ( 16
)
By directly computing the relevant quantity i.e. the tube survival probability under reptation, contour-length fluctuations and thermal constraint release, we alleviate the need for a non-linear mixing rule.

In order to correctly model the short-time behaviour, the Rouse relaxation modulus (11) can be added for each species using a linear mixing rule.

Numerical solution of the model

Although the above theory is quite simple from a mathematical point of view, a closer look at Eq. ( 14) shows that the differential problems for all masses are coupled through the constraint release term. This is numerically expensive for a large number of masses. The numerical cost can however be reduced dramatically through the use of a semi-analytical solution that decouples the different masses.

In the absence of constraint release, i.e. when γ = 0, the partial differential Eq. ( 14) become decoupled. The functions P f (i) 0 (t, s) can therefore be computed numerically at a reasonable cost using a finite difference scheme. Through direct substitution, one can verify that the following expression for P f (i) γ (t, s) satisfies Eq. ( 14):

P f (i) γ (t, s) = P f (i) 0 (t, s) N i=1 φ (i) 2 1 -1 P f (i) 0 (t, s) ds γ . ( 17
)
The computation of P f (i) γ

with this method has a truly negligible numerical cost as, for a large number of masses, the algorithmic complexity is reduced from quadratic to linear in N.

Model predictions

We now present predictions of the model in the monodisperse case. Next, we apply the model for the prediction of linear viscoelastic moduli of different samples of polydisperse polystyrene melts. For all predictions, a Rouse spectrum described by Eq. ( 11) has been superposed using a linear mixing rule, as done by van Ruymbeke et al. [START_REF] Van Ruymbeke | Evaluation of reptation models for predicting the linear viscoelastic properties of entangled linear polymers[END_REF].

In Fig. 1, we show the molecular weight dependence of the predicted zero-shear viscosity of a monodisperse system. In the limit of highly entangled chains, the viscosity scales like the cube of the molar mass. As the number of entanglements gets lower, this scaling gradually changes to an exponent 3.4. In this example we chose to set γ = 1 and K f = 0.6, as these are the values we identified for the polystyrene melts discussed later in this section. The value of the other parameters does not change the shape of the viscosity curve but merely translates it vertically and horizontally. For the sake of completeness, we compare our predictions with those of Likthman and McLeish (C ν = 1) [START_REF] Likhtman | Quantative theory for linear dynamics of linear entangled polymers[END_REF]. The models predicts similar trends and apparently the same asymptotic value. The differences showing up for less entangled systems are mostly due to differences in the relative importance of contour-length fluctuations between the two models. With γ = 1.2 and K f = 2 our model predicts a viscosity curve very close to that of Likhtman and McLeish.

To evaluate the model further, we test its ability to correctly predict the linear viscoelastic properties of six polydisperse polystyrene samples. The parameters of the model are estimated using the molecular weight distribution and the viscoelastic moduli of a single sample. The value of the parameters is then frozen for the subsequent predictions on the remaining five samples. The data of the six samples named PS1, PS2, PS3, PS60, PS275 and PS330 were kindly provided to us by BASF. We also received the molecular weight distribution and the viscoelastic moduli at 170 • C from BASF. The characteristics of the various samples are provided in Table 1. Samples PS1, PS60 and PS275 are quite monodisperse, while PS2 and PS330 have a much broader molecular weight distribution. PS3 is a trimodal sample with a fraction of very high masses. The parameter estimation is actually quite a delicate task for two reasons. In the first place, the sensitivity of the model to the different parameters can be very uneven. Second, not all samples are suited to perform this estimation, if they only provide information in a narrow molecular weight range which is not representative of the physics built in the theory. Therefore, we chose to use the broad sample PS3 to fit the parameters of the model. In order to account for the whole molecular weight distribution, we discretize it with a small fixed step in logarithmic space. The efficient semi-analytical method described in Section 4 still allows us to compute quickly the model predictions. The parameters resulting from the fitting procedure are reported in Table 2. The molecular weight distribution as well as the measured and predicted viscoelastic moduli for PS3 are reported in Fig. 2(a andb). We see there a good quantitative agreement between the data and the predictions, over a wide range of frequencies.

In Figs. 234we show, for a wide range of molecular weight distributions, the quantitative agreement between experimentally measured viscoelastic moduli and the model's predictions. All these predictions were obtained using the single set of parameters from Table 2. The main discrepancies between predictions and measurements are found at low and high frequencies. At low frequency, the model sometimes fails to correctly predict the storage modulus. In this frequency range the dominating mechanisms are reptation and constraint release. Furthermore, the parameter γ, which controls constraint release, has the main influence on the slope of the storage modulus in this region. We would therefore suggest that a more accurate treatment of constraint release would help resolve these discrepancies. At intermediate and high frequencies the loss modulus tends to be underestimated. We connect this feature to the crude description of contour-length fluctuations we implemented.

Linear approximation of the theory

Despite the relative simplicity of our theory, it remains mathematically non-linear. The non-linear operator on P f γ (t, s) found in the definition makes the rate of constraint release β time dependent. This non-linearity is not a problem per se but having to rely on a mathematically non-linear theory to describe linear rheological properties of a system prevents us from building a theory for the rheologicaly non-linear regime on top of the previous one. In the present section we propose a simple ap- proximation of our theory that preserves its general mathematical structure, but where both the mixing rule and the operators acting on the variables of the model are mathematically linear.

The non-linearity of our theory comes from the timedependent relaxation induced by constraint release. Therefore, we propose to approximate constraint release through a superposition of few linear modes. Let us call Q CR (t) or constraint release kernel the second factor of the right-hand-side of Eq. ( 17). This factor represents the effects of constraint release on P f (i) 0 (t, s):

Q CR (t) = N i=1 φ (i) 2 1 -1 P f (i) 0 (t, s) ds γ . ( 18
)
Let us assume that the constraint release kernel Q CR (t) can be approximated by a finite sum of decreasing time-exponentials:

Q CR (t) ≈ N CR j=0 w j exp(-t/τ j ), (19)
where w j and τ j are adjustable parameters. For a given of N, we compute w j and τ j using a minimisation procedure for the the quadratic error between the exact constraint release kernel [START_REF] Cloizeaux | Relaxation and viscosity anomaly of melts made of long entangled polymers[END_REF] and its approximation [START_REF] Marrucci | Relaxation by reptation and tube enlargement: A model for polydisperse polymers[END_REF]. We can now define P f (i,j) γ (t, s) as the solution of the following linear differential problem:

∂P f (i,j) γ t = ∂ ∂s α f (i) d (s) ∂ ∂s P f (i,j) γ - 1 τ j P f (i,j) γ , ( 20
) P f (i,j) γ (t, -1) = 0, P f (i,j) γ (t, 1) = 0 for t > 0, P f (i,j) γ (0, s) = 1 for -1 < s < 1.
As for Eq. ( 14), the analytical solution of Eq. ( 20) can be expressed as a function of P f (i) 0 :

P f (i,j) γ (t, s) = P f (i) 0 (t, s) exp(-t/τ j ). (21)
Making use of Eq. ( 19), one can verify that the following expression for the relaxation modulus correctly approximates Eq. ( 16) of the original theory:

G(t) = G 0 N N i=0   φ (i) 2 1 -1   N CR j=1 P f (i,j) γ (t, s)   ds   . ( 22
)
In the previous expression, the sub-expression contained within the innermost parenthesis actually is an approximation of P f (i) γ (t, s). The quality of this linear approximation depends only on the quality of the approximation of the constraint release kernel (18) by a sum of decreasing exponentials. In Fig. 5, we show for sample PS2 that, with seven modes, the linear approximation is almost indistinguishable from the original model.

Conclusions

We have presented a simple tube theory for linear entangled systems. The specificity of this theory is that, for a monodisperse system, a single partial differential equation is used to compute the tube survival probability under reptation, contourlength fluctuations and thermal constraint release. Reptation and contour-length fluctuations are modelled simultaneously by position-dependent diffusion along the primitive path. Thermal constraint release is introduced through a time dependent relaxation term. The intensity of the thermal constraint release is controlled by reptation and contour-length fluctuations. In a polydisperse environment, we directly compute the tube survival probability under the influence of all masses with a set of coupled partial differential equations. Additionally, we proposed a semi-analytical method for the numerical solution of our model. This method reduces the algorithmic complexity from quadratic to linear in terms of the number of masses, thus making the computational cost negligible even for a broad molecular weight distribution.

Although some discrepancies remain, we have shown the good predictive properties of our theory on a wide range of molten polystyrene samples.

Finally, we have proposed an approximation of our theory, that enables us to address the linear viscoelastic properties of entangled linear polymers with a mathematically linear theory. This approximation is suited as a starting point for building constitutive equations for the non-linear regime. Work on this topic is in progress.

Fig. 1 .

 1 Fig. 1. Dependence of the zero-shear viscosity η 0 with the molar mass for a strictly monodisperse sample. The plain line is the prediction of our model (γ = 1 and K f = 0.6), while the dashed line is the prediction of the Likhtman and McLeish model [6]. The dotted line is the prediction of our model where we increased the depth of contour length fluctuations and the effectiveness of constraint release (γ = 1.2 and K f = 2).

Fig. 2 .

 2 Fig. 2. Molecular weight distribution and linear viscoelastic moduli of PS3 (a and b) and PS1 (c and d). Plain lines are the model's output, while the dots are the experimental data.

Fig. 3 .

 3 Fig. 3. Molecular weight distribution and linear viscoelastic moduli of PS2 (a and b) and PS60 (c and d). Plain lines are the model's output, while the dots are the experimental data.

Fig. 4 .

 4 Fig. 4. Molecular weight distribution and linear viscoelastic moduli of PS275 (a and b) and PS330 (c and d). Plain lines are the model's output, while the dots are the experimental data.

Fig. 5 .

 5 Fig. 5. Comparison between the viscoelastic moduli of the full model (-) and its linear approximation(-•-•-). Seven modes were used to approximate the constraint release kernel of the full model.

Table 1

 1 Characteristics of the six PS samples provided by BASF

	Name	M w (kg/mol)	M n (kg/mol)	M w /M n
	PS1	320	270	1.18
	PS2	274	101	2.72
	PS3	407	143	2.83
	PS60	69.6	6 2 .7	1.11
	PS275	290	253	1.07
	PS330	324	112	2.89

Table 2

 2 Parameters resulting from the fitting of the model on the PS3 data

	G 0 N τ e	2 × 10 5 Pa 5.0 × 10 -4 s
	M e	15700 g/mol
	K f	0.6
	γ	1.0

Acknowledgements

This work is supported by the "Action de recherche concertée" programme of the Communauté Franc ¸aise de Belgique and by "BASF Aktiengesellschaft". The authors wish to thank Dr. Martin Laun for making the PS data available.


Le Ha

Vy Nguyen
email: nguyenlehavy@gmail.com

Catherine Bonnet
email: catherine.bonnet@inria.fr

L H V Nguyen

Stabilization of MISO fractional systems with delays

Keywords: Coprime factorization, H-infinity, stabilizing controllers, time delay, MIMO

published or not. The documents may come

Introduction

Fractional systems are systems described by differential equations involving non-integer order derivatives and/or integrals. Consequently, in the frequency domain, their transfer functions contain non-integer powers of the Laplace variable s. This kind of models has become more popular in many fields in the past two decades since it provides a better fit to data being then more succinct than a standard model. Refer, for example, to [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] for basic backgrounds on fractional calculus and to [START_REF] Sabatier | Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering[END_REF][START_REF] Freeborn | A survey of fractional-order circuit models for biology and biomedicine[END_REF] for its recent applications on modeling.

Delays are encountered almost everywhere due, for example, to distance of transmission and it is well-known that they have important effects on the stability of the systems [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF].

While integer-order systems with delays have been intensively studied [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF], the literature on fractional systems with delays is still quite small. Particularly, the stabilization problem of fractional systems with delays has rarely been addressed. Some available studies are classical [START_REF] Özbay | PID controller design for fractional-order systems with time delays[END_REF] and fractional PID controller design [START_REF] Hamamci | An algorithm for stabilization of fractional-order time delay systems using fractional-order PID controllers[END_REF], fractional sliding mode control [START_REF] Si-Ammour | A sliding mode control for linear fractional systems with input and state delays[END_REF], factorization approach to control synthesis [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF][START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF].

In the framework of fractional representation approach to synthesis problems [START_REF] Vidyasagar | Control System Synthesis: A Factorization Approach[END_REF], SISO fractional delay systems was considered in [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF][START_REF] Bonnet | Stabilization of some fractional delay systems of neutral type[END_REF] and coprime factorizations together with the corresponding Bézout factors of the transfer function of these systems have been derived. In [START_REF] Curtain | Coprime factorization for regular linear systems[END_REF], coprime factors were presented for a large class of MIMO infinite-dimensional systems which include delay systems. The factors were determined from a state-space realization of the (regular) system which was given in terms of the semigroup of the system. Such realizations are not much considered for fractional systems.

For the particular class of MIMO (integer-order) systems with I/O delays, the problem of parametrization of stabilizing controllers was solved in [START_REF] Mirkin | State-space parametrization of all stabilizing dead-time controllers[END_REF][START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple I/O delays[END_REF]. The idea was to reduce the problem to an equivalent finite-dimensional stabilization problem by involving an unstable finitedimensional system and a stable infinite-dimensional system (FIR filter). In [START_REF] Mondié | Finite eigenstructure assignment for input delay systems[END_REF], a procedure to compute right coprime factorizations over a Bézout domain was proposed for spectrally controllable MIMO (integer-order) systems with input delays. For MISO structure, a class of (integer-order) systems with multiple transmission delays was studied in [START_REF] Bonnet | PID stabilization of SISO delay systems and robust stabilization of systems with multiple transmission delays[END_REF] and coprime factorizations and associated Bézout factors over H ∞ were derived.

In this paper, we are interested in the stabilization problem of MISO fractional systems with different I/O delays which are not necessarily commensurate. This MISO structure appeared, for example, in communication systems [START_REF] Quet | Rate-based flow controllers for communication networks in the presence of uncertain time-varying multiple time-delays[END_REF]. We would like to obtain the set of all stabilizing controllers by determining a doubly coprime factorization over H ∞ of the transfer matrix and the associated Bézout factors, which allow the construction of the Youla-Kučera parametrization [START_REF] Vidyasagar | Control System Synthesis: A Factorization Approach[END_REF]. As in the finitedimensional case, the Youla-Kučera parametrization gives the set of all H ∞ -stabilizing controllers in terms of one free parameter. Note that in [START_REF] Quadrat | On a generalization of the Youla-Kučera parametrization, Part II: The lattice approach to MIMO systems[END_REF], a parametriza-tion of the set of all stabilizing controllers is given in terms of two free parameters for MIMO systems once we already know a particular stabilizing controller. Our strategy here is to work directly on the Bézout identity in order to get explicit expressions of Bézout factors in terms of the matrix transfer function. Such explicit expressions could not be easily derived in [START_REF] Mirkin | State-space parametrization of all stabilizing dead-time controllers[END_REF][START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple I/O delays[END_REF][START_REF] Mondié | Finite eigenstructure assignment for input delay systems[END_REF] even in the case of standard delay systems. We hope that the explicit form will facilitate the use of these factors in controllers design while the use of the frequency domain representation of the systems agrees well with the modeling practice of fractional systems [START_REF] Sabatier | Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering[END_REF].

The paper is organized as follows. In Section 2, the class of systems of interest and some background are presented. The results are stated in Sections 3 and 4. We give in Section 3 explicit expressions of left coprime factorizations and associated Bézout factors over H ∞ of the transfer function of the systems under study. Right coprime factorizations and right Bézout factors are given in Section 4 for a large subclass of the class of systems considered. Examples are provided to illustrate the results. Finally, Section 5 gives conclusions and perspectives.

A class of MISO fractional time-delay systems

We consider systems described by transfer matrices of the form

G(s) = e -sh1 R 1 (s α ), . . . , e -shn R n (s α ) , (1)
where 0 ≤ h k ∈ R for k = 1, . . . , n are the delays; α ∈ R, 0 < α < 1; R k (s α ) = q k (s α )/ p k (s α ), where p k and q k are polynomials of integer degree in s α , p k (s α ) and q k (s α ) have no common roots, and deg p k ≥ deg q k for k = 1, . . . , n; d k is the degree in s α of p k ; s is in the principle branch C\R -, that is arg(s) ∈ (-π, π), in order to guarantee a unique value of the transfer function involving s α with α ∈ (0, 1).

Some notations used are

C + = {s ∈ C | Re(s) > 0}, Z + = {p ∈ Z | p > 0}, Z + = {p ∈ Z | p ≥ 0}.
We are interested in H ∞ -stability, i.e. a SISO system is stable if its transfer function K(s) belongs to the H ∞ space of analytic and bounded functions in

C + with K H∞ = sup s∈C+ |K(s)|. Let us denote M(H ∞ ) the set of matrices whose components are in H ∞ .
The following notion of coprimeness is considered.

A system G is said to have a right coprime factorization (r.c.f.) (N, M ) over

H ∞ if G = N M -1 , det M = 0, N , M ∈ M(H ∞ ) and there exist X, Y ∈ M(H ∞ ) such that XM + Y N = I. Then X, Y are called right Bézout factors.
A system G is said to have a left coprime factorization (l.c.f.) ( M , N ) over If G has an r.c.f. (N, M ) and an l.c.f. ( M , N ), then the set of all controllers guaranteeing the internal stability of the closed-loop system is given by the Youla-Kučera parametrization

H ∞ if G = M -1 N , det M = 0, M , N ∈ M(H ∞ ) and there exist X, Y ∈ M(H ∞ ) such that M X + N Y = I. Then X, Y are called left Bézout factors. For M , N ∈ M(H ∞ ), there exist X, Y ∈ M(H ∞ ) such that M X+ N Y = I if and only if inf s∈C+ σ m ([ M , N ] T ) > 0,
C(G) = {(X -R N ) -1 (Y + R M ) | R ∈ M(H ∞ ), det(X -R N ) = 0} = {( Y + M R)( X -N R) -1 | R ∈ M(H ∞ ), det( X -N R) = 0},
where X, Y and X, Y are respectively the corresponding right and left Bézout factors [START_REF] Vidyasagar | Control System Synthesis: A Factorization Approach[END_REF]. For R = 0, we obtain two particular stabilizing controllers

C = X -1 Y and C = Y X -1 .
Poles (resp. roots) in the closed right half-plane C + are referred to as unstable poles (resp. roots).

The following notations will be of intense use later. Denote p(s α ) the lowest common denominator of R k (s α ) for k = 1, . . . , n; d the degree in s α of p(s α ). Then rational transfer functions R k (s α ) can be rewritten as

R k (s α ) = q k (s α ) p(s α ) ,
where q k are polynomials in s α .

We can decompose

p(s α ) = (s α ) m0 N i=1 (s α -b i ) mi   N j=1 (s α -c j ) m j   , where b i ∈ D := {σ ∈ C\{0} | -πα/2 ≤ Arg(σ) ≤ πα/2}; c j ∈ C\{D ∪ {0}}; m 0 ∈ Z + , m i , m j ∈ Z + for i = 1, .
. . , N and j = 1, . . . , N . Hence

s i = b 1/α i
are the non-zero unstable roots in s of p(s α ).

Similarly, we write

p k (s α ) = (s α ) m 0k N i=1 (s α -b i ) m ik   N j=1 (s α -c j ) m jk   ,
where m 0k , m ik , m jk ∈ Z + for i = 1, . . . , N , j = 1, . . . , N and k = 1, . . . , n. It is obvious that m 0k ≤ m 0 , m ik ≤ m i , and m jk ≤ m j .

Left coprime factorizations and Bézout factors

In this section, we present l.c.f.'s and Bézout factors for the transfer matrix (1).

Left coprime factorizations

Due to the dimension of the transfer matrix, finding an l.c.f. is straightforward. The next proposition was presented in [START_REF] Nguyen | Coprime factorizations of MISO fractional time-delay systems[END_REF] and is recalled here together with its proof for the paper to be self-contained.

Proposition 1 ( [START_REF] Nguyen | Coprime factorizations of MISO fractional time-delay systems[END_REF]) Let G(s) be given by [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF]. Then ( M (s), N (s)), where

M (s) = p(s α ) (s α + 1) d , N (s) = 1 (s α + 1) d e -sh1 q 1 (s α ), . . . , e -shn q n (s α ) , (2)
is an l.c.f. over H ∞ of G(s). PROOF. It is obvious that M (s) -1 N (s) = G(s).
We see that M (s) ∈ H ∞ . Also, each component of N (s) is in H ∞ , and then N (s) ∈ M(H ∞ ).

For all roots σ of p(s α ), there exists at least one 1 ≤ k ≤ n such that q k (σ) = 0. Thus inf s∈C+ (

n k=1 | N k |+| M |) > 0 which ensures that ( M , N ) is an l.c.f. over H ∞ of G.

Bézout factors

Now we propose left Bézout factors corresponding to the l.c.f. obtained above. From the left Bézout identity, we derive that

X = M -1 (1 -N Y ). The idea is to choose Y ∈ M(H ∞ ) such that X ∈ H ∞ . To achieve that, we interpolate (1 -N Y ) at unstable zeros of M .
For some simple classes of systems (1), Bézout factors were derived in [START_REF] Nguyen | Coprime factorizations of MISO fractional time-delay systems[END_REF]. In this paper, we consider the most general case of systems [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF] where techniques developed in [START_REF] Bonnet | PID stabilization of SISO delay systems and robust stabilization of systems with multiple transmission delays[END_REF][START_REF] Nguyen | Coprime factorizations of MISO fractional time-delay systems[END_REF] are no longer applicable and we need here to include non-commensurate powers of s in X and Y .

First, let us denote

k 0 := min{k | k ∈ {1, . . . , n}, m 0k = m 0 }, (3)
k i := min{k | k ∈ {1, . . . , n}, m ik = m i }, i = 1, . . . , N, (4)
f k := i∈{1,...,N }:ki=k m i , k = 1, . . . , n, L(m 0 , α) := {x ∈ R | x = a + bα < m 0 α, a, b ∈ Z + }.
The latter set will be of use to describe the power in s of the product s a (s α ) b with a, b ∈ Z + .

Proposition 2 Let G(s) be given by [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF]. Then Bézout factors corresponding to the l.c.f. ( 2) are given by

X(s) = (s α + 1) d u(s α ) - n k=1 e -sh k q k (s α )µ k (s) p(s α )u(s α ) , Y (s) = µ 1 (s) u(s α ) , . . . , µ n (s) u(s α ) T ,
where u(s α ) is a polynomial in s α of degree greater or equal to d whose zeros are stable, and for k = 1, . . . , n the polynomials µ k (s) have the following form

µ k (s) =                λ∈L(m0,α) β λk s λ + m0+f k -1 j=m0 β (jα)k (s α ) j if k = k 0 , f k -1 j=0 β (jα)k (s α ) j if k = k 0 ,
and satisfy

(s α + 1) d u(s α ) - n k=1 e -sh k q k (s α )µ k (s) = O(s m0α )
(5) as s → 0 and

(s α + 1) d u(s α ) - n k=1 e -sh k q k (s α )µ k (s) (l) = 0 (6) for each non-zero unstable root s = b 1/α i , i = 1, . . . , N , of p(s α ) and for 0 ≤ l ≤ m i -1. Remark 1 If f k = 0, then µ k (s) = λ∈L(m0,α) β λk s λ when k = k 0 and µ k (s) = 0 otherwise.
PROOF. It is easy to verify that X(s) and Y (s) satisfy the left Bézout identity.

The degree of µ k (s) is smaller than or equal to the degree of u(s α ), and so Y (s) ∈ M(H ∞ ).

We see that X(s) is bounded at ∞ in C + . Moreover, due to (6), the numerator of X(s) has the same non-zero unstable roots as the denominator, which assures that X(s) is analytic in C + . The boundedness of the function at s = 0 is satisfied by [START_REF] Freeborn | A survey of fractional-order circuit models for biology and biomedicine[END_REF]. Then X(s) ∈ H ∞ . Now it remains to prove the existence of µ k (s) satisfying the two conditions ( 5) and ( 6). The details are reported in Appendix A.

Remark 2 If m 0 α ≤ 1 (this case has been considered in [START_REF] Nguyen | Coprime factorizations of MISO fractional time-delay systems[END_REF]) or α = 1/m with m ∈ Z + \{0, 1}, then λ are multiples of α and we obtain an elegant formula of µ k0 which only contains the terms in s α . More generally, if α is rational, then µ k0 contains powers of s of commensurate exponents.

The following example illustrates the case where α is irrational.

Example 1

G(s) = e -s s π/2 (s π/4 -1) 2 , e -3s s π/4 -1
We have p(s π/4 ) = s π/2 (s π/4 -1) 2 with degree d = 4 in s π/4 . Its unstable roots are b 0 = 0 and b 1 = 1 with multiplicity m 0 = 2 and m 1 = 2 respectively. Obviously, q 1 (s π/4 ) = 1 and q 2 (s π/4 ) = s π/2 (s π/4 -1). Then from Proposition 1, we obtain an l.c.f. as follows

M (s) = s π/2 (s π/4 -1) 2 (s π/4 + 1) 4 ,
N (s) = 1 (s π/4 + 1) 4 e -s , e -3s s π/2 (s π/4 -1) .

To complete the expressions of the Bézout factors given in Proposition 2, we now choose u(s π/4 ) and search for µ 1 (s) and µ 2 (s) by solving the equations imposed by the two conditions ( 5), [START_REF] Hamamci | An algorithm for stabilization of fractional-order time delay systems using fractional-order PID controllers[END_REF].

Here, we choose u(s π/4 ) = (s π/4 + 1) 4 . It is easy to see that

k 0 = 1, k 1 = 1, f 1 = 2, f 2 = 0 and L(m 0 , α) = L(2, π/4) = {0, π/4, 1}. Therefore, µ 1 (s) and µ 2 (s) have the forms µ 1 (s) = β 01 + β (π/4)1 s π/4 + β 11 s + β (π/2)1 s π/2 + β (3π/4)1 s 3π/4 , µ 2 (s) = 0.
The numerator of X(s) is then (s π/4 + 1) 8 -e -s µ 1 (s).

Its development around zero is

(1 -β 01 ) + (8 -β (π/4)1 )s π/4 + (β 01 -β 11 )s + O(s π/2 ).
The condition [START_REF] Freeborn | A survey of fractional-order circuit models for biology and biomedicine[END_REF] implies that all powers of s with degree smaller than π/2 vanish, thus leads to β 01 = 1, β (π/4)1 = 8, and β 11 = 1.

The other coefficients are derived from the condition [START_REF] Hamamci | An algorithm for stabilization of fractional-order time delay systems using fractional-order PID controllers[END_REF], which is represented by

(s π/4 + 1) 8 -e -s µ 1 (s) = 0, [(s π/4 + 1) 8 -e -s µ 1 (s)] = 0 at s = 1.
The unique solution of these two equations is

β (π/2)1 = -2(11π + 128eπ -2 + 512e)/π and β (3π/4)1 = 4(3π + 128eπ -1 + 256e)/π.
Hence, the Bézout factors are

X(s) = (s π/4 + 1) 8 -e -s µ 1 (s) s π/2 (s π/4 -1) 2 (s π/4 + 1) 4 , Y (s) = µ 1 (s) (s π/4 + 1) 4 , 0 T , where µ 1 (s) = 1 + 8s π/4 + s - 2(11π + 128eπ -2 + 512e) π s π/2 + 4(3π + 128eπ -1 + 256e) π s 3π/4 .

Right coprime factorizations and Bézout factors

The previous section showed that the systems G(s) under study admit l.c.f.'s over H ∞ . Since H ∞ is a Hermite ring, then there exist r.c.f.'s for G(s) [START_REF] Smith | On stabilization and the existence of coprime factorizations[END_REF].

While l.c.f.'s and left Bézout factors of our transfer matrices are vectors and scalars, the right counterparts are vectors and square matrices with a lot more entries to be determined. Regarding calculation complexity, we separate the study into two classes of systems. The first class consists of systems with distinct unstable poles, i.e. p k (s α ) and p k (s α ) have no common unstable roots if k = k , for which M (s) can be chosen to be of diagonal form. This simple form cannot be applied to the second class of systems with identical unstable poles, for which our first attempt is to consider a particular case.

Some particular cases of the systems considered here have been studied in [START_REF] Nguyen | Right coprime factorizations of MISO fractional time-delay systems[END_REF].

Distinct poles

In this context of determining right Bézout factors, the choice of µ k for k = 1, . . . , n is not very different from that presented in Proposition 2.

Proposition 3 Let G(s) be given by [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF]. Suppose that all (zero and non-zero) unstable roots of p k (s α ) for k = 1, . . . , n are distinct. Then one r.c.f. and associated Bézout factors are given by

N (s) = [N 1 (s), . . . , N n (s)], M (s) =      M 11 (s) • • • 0 . . . . . . . . . 0 • • • M nn (s)      , X(s) =      X 11 (s) • • • X 1n (s) . . . . . . . . . X n1 (s) • • • X nn (s)      , Y (s) = [Y 1 (s), . . . , Y n (s)] T ,
where for k, k ∈ {1, . . . , n} and k = k

N k (s) = e -sh k q k (s α ) (s α + 1) d k , (7)
M kk (s) = p k (s α ) (s α + 1) d k , (8)
Y k (s) = µ k (s) u(s α ) 1≤j≤n,j =k (s α ) m0j N i=1 (s α -b i ) mij , X kk (s) = 1 -Y k (s)N k (s) M kk (s) , X kk (s) = -Y k (s)e -sh k q k (s α ) p k (s α ) ,
where d k is the degree of p k in s α ; u is a polynomial of degree d in s α that has no unstable zeros; and µ k (s) have the following form

µ k (s) = λ∈L(m 0k ,α) β λk s λ + m 0k + N i=1 m ik -1 j=m 0k β (jα)k (s α ) j , satisfying u(s α )(s α + 1) d k -e -sh k µ k (s) q k (s α ) × 1≤j≤n,j =k N i=1 (s α -b i ) mij = O(s m 0k α ) ( 9
)
as s → 0 if p k (s α ) has a root at zero, and

u(s α )(s α + 1) d k -e -sh k µ k (s) q k (s α ) × 1≤j≤n,j =k (s α ) m0j N i=1 (s α -b i ) mij   (l)
= 0 [START_REF] Mondié | Finite eigenstructure assignment for input delay systems[END_REF] for each non-zero unstable root of p k (s α ), i.e. s = b

1/α i with m ik = 0 for i = 1, . . . , N , where l = 0, . . . , m ik -1.

PROOF. It is obvious that N k (s), M kk (s), Y k (s), X kk (s) ∈ H ∞ .
The two conditions ( 9) and ( 10) guarantee X kk (s) ∈ H ∞ for k = 1, . . . , n. We also see that G(s) = N (s)M (s) -1 and the right Bézout identity X(s)M (s) + Y (s)N (s) = I is verified.

To complete the proof, we prove the existence of µ k (s) satisfying the two conditions in Appendix B.

Here is a numerical example of the use of the proposition.

Example 2

G(s) = e -s s π/2 (s π/4 -1) 2 , e -3s s π/4 -2
By applying Proposition 3 we obtain the r.c.f. as follows

N (s) = e -s (s π/4 + 1) 4 , e -3s s π/4 + 1 , M (s) =   s π/2 (s π/4 -1) 2 (s π/4 +1) 4 0 0 s π/4 -2 s π/4 +1   .
We choose u(s π/4 ) = (s π/4 + 1) 5 , then Y (s) has the form

Y (s) = µ 1 (s)(s π/4 -2) (s π/4 + 1) 5 , µ 2 (s)s π/2 (s π/4 -1) 2 (s π/4 + 1) 5 T
where µ 1 (s) = β 01 + β (π/4)1 s π/4 + β 11 s + β (π/2)1 s π/2 + β (3π/4)1 s 3π/4 and µ 2 (s) = β 02 .

The condition ( 9) is only applied for k = 1. We develop the numerator of X 11 (s) around zero as follows

(s π/4 + 1) 9 -e -s µ 1 (s)(s π/4 -2) = (1 + 2β 01 ) + (9 -β 01 + 2β (π/4)1 )s π/4 + 2(β 11 -β 01 )s + O(s π/2 ).
It turns out that all the terms with order smaller than π/2 in the development have to be zero, thus giving β 01 = -1/2, β (π/4)1 = -19/4, and

β 11 = -1/2.
Other unknown coefficients are deduced from applying the condition [START_REF] Mondié | Finite eigenstructure assignment for input delay systems[END_REF]. We have (s π/4 + 1) 9 -e -s µ 1 (s)(s π/4 -2) = 0, [(s π/4 + 1) 9 -e -s µ 1 (s)(s π/4 -2)] = 0 at s = 1 and (s π/4 + 1) 6 -e -3s µ 2 (s)s π/2 (s π/4 -1) 2 = 0 at s = 2 4/π , and the unique solution of the above equations is β (π/2)1 = (25π + 2560eπ -4 + 4096e)/(2π), β (3π/4)1 = -(27π +7168eπ -8+8192e)/(4π), and β 02 = (729(e 2 4/π ) 3 )/4.

Finally, X(s) = X 11 (s) X 12 (s)

X 21 (s) X 22 (s)
where

X 11 (s) = (s π/4 + 1) 9 -e -s µ 1 (s)(s π/4 -2) s π/2 (s π/4 -1) 2 (s π/4 + 1) 5 , X 12 (s) = -e -3s µ 1 (s) (s π/4 + 1) 5 , X 21 (s) = -e -s µ 2 (s) (s π/4 + 1) 5 , X 22 (s) = (s π/4 + 1) 6 -e -3s µ 2 (s)s π/2 (s π/4 -1) 2 (s π/4 -2)(s π/4 + 1) 5 , µ 1 (s) = - 1 2 - 19 4 s π/4 - 1 2 s + 25π + 2560eπ -4 + 4096e 2π s π/2 - 27π + 7168eπ -8 + 8192e 4π s 3π/4 , µ 2 (s) = 729(e 2 4/π ) 3 4 .

Identical unstable poles

In the case of identical unstable poles, the following lemma characterizes the form of the matrix M (s).

Lemma 4 Let G(s) be given by [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF]. Suppose that there exist k, k ∈ {1, . . . , n}, k = k such that p k (s α ) and p k (s α ) have common unstable root(s). Then the factor M (s) of an r.c.f. of G(s) cannot be a diagonal matrix.

PROOF. Let us prove the fact by contradiction. As

- sume that M (s) is diagonal, then G k (s) = N k (s)/M kk (s), k = 1, . . . , n. Since N k (s), M kk (s) ∈ H ∞ ,
X kk (s) = 1 -Y k (s)N k (s) M kk (s) , X kk (s) = -Y k (s) N k (s) M k k (s) .
Since X kk (s) ∈ H ∞ , then all unstable zeros of M k k (s) are zeros of Y k (s). Now, if G k (s) and G k (s) have a common unstable pole s 0 , then s 0 is a common zero of M kk (s) and M k k (s). Also Y k (s 0 ) = 0, and thus X kk (s 0 ) = ∞, which contradicts the fact that X kk (s) ∈ H ∞ .

In the following part, we derive r.c.f.'s and Bézout factors for a particular system whose entries (which only have one simple pole) may involve identical poles. To help clarify the demonstration of those results, we will begin with a lemma who derives the inverse of a particular upper triangular matrix.

We consider sparse matrices with the following condition imposed on the entries above the main diagonal: if any entry on the k-th column is non-zero, then all other entries on the k-th column as well as those on the k-th row must be zeros.

Lemma 5 Let M ∈ R n×n be an upper triangular matrix whose entries on the main diagonal are not equal to zero and the entries above the main diagonal satisfy (A) for k = 1, . . . , n, if there exists l ∈ Z, l ∈ [1, k) such that M l k = 0 then M lk = 0 for l ∈ {1, . . . , k -1}\{l } and M kl = 0 for l ∈ {k + 1, . . . , n}.

Then its inverse is given by M -1 = M inv , where M inv is an upper triangular matrix whose entries on and above the main diagonal satisfy

M inv kk = 1 M kk , (11)
M inv kk = - M kk M kk M k k ( 12
)
for k, k ∈ {1, . . . , n} and k < k .

PROOF. It is obvious that the entries below the main diagonal of the product M M inv are all zero, and the entries on the main diagonal are all one. Now we consider the entries above the main diagonal of the product. For i < j,

(M M inv ) ij = n k=1 M ik M inv kj = j k=i M ik M inv kj .
Considering i < k < j, if M ik = 0, then M kj = 0 under the assumption (A), and thus M inv kj = 0 due to [START_REF] Nguyen | Right coprime factorizations of MISO fractional time-delay systems[END_REF]. Therefore,

(M M inv ) ij = M ii M inv ij + M ij M inv jj .
By replacing M inv ij and M inv jj with ( 12) and ( 11) respectively, we obtain (M M inv ) ij = 0.

Example 3 M =          M 11 0 M 13 0 0 0 M 22 0 M 24 M 25 0 0 M 33 0 0 0 0 0 M 44 0 0 0 0 0 M 55         
with M ii = 0, i = 1, . . . , 5, and M 13 , M 24 , M 25 = 0. This matrix satisfies all the conditions in Lemma 5 and thus its inverse is

M -1 =          1 M11 0 -M13 M11M33 0 0 0 1 M22 0 -M24 M22M44 -M25 M22M55 0 0 1 M33 0 0 0 0 0 1 M44 0 0 0 0 0 1 M55          .
In the following proposition, we consider G(s) with one pole for each of its entries and some entries may have the same pole. To simplify the presentation, we assume that the delays are ordered. A discussion on how to apply the next result to the case of unordered delays will follow the proposition.

Proposition 6 Let G(s) be given by ( 1) with h 1 ≤ . . . ≤ h n and

R k (s α ) = a k s α -σ k ,
where a k , σ k ∈ R for k = 1, . . . , n. We denote I 1 := ∅ and

I k := {j | j ∈ {1, . . . , k -1}, σ j = σ k } for k = 2, . . . , n.
One r.c.f. and associated Bézout factors are given by

N (s) = [N 1 (s), . . . , N n (s)], (13)
M (s) =      M 11 (s) • • • M 1n (s) . . . . . . . . . 0 • • • M nn (s)      , (14)
Y (s) = [Y 1 (s), . . . , Y n (s)] T , X(s) = M -1 (s) -Y (s)G(s),
where for k, k ∈ {1, . . . , n} and k = k

N k (s) = 0 if I k = ∅, a k e -sh k s α +1 otherwise, (15)
M kk (s) = 1 if I k = ∅ s α -σ k s α +1 otherwise, (16)
M k k (s) =    -a k e -s(h k -h k ) a k if k = min I k , 0 otherwise, (17)
Y k (s) = 0 if I k = ∅, β k pk (s α ) u(s α )
otherwise, where u(s α ) is a polynomial of degree d in s α that has no unstable zeros; pk (s α ) = p(s α ) s α -σ k ; the coefficients β k , for those k such that I k = ∅ and σ k ≥ 0, are given by

β k = u(σ k )(σ k + 1)e σ 1/α k h k a k pk (σ k ) , ( 18
)
while β k for other k can be chosen arbitrarily, and

M -1 (s) = M inv (s)
, where M inv (s) is an upper triangular matrix whose entries on and above the main diagonal are given by

M inv kk (s) = 1 M kk (s) , M inv kk (s) = - M kk (s) M kk (s)M k k (s) (19)
for k, k ∈ {1, . . . , n} and k < k .

PROOF. Let us prove that M (s) given by ( 14), [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF], and ( 17 Consequently, due to Lemma 5, the inverse of M (s) is M inv (s) whose entries are given by [START_REF] Smith | On stabilization and the existence of coprime factorizations[END_REF].

We now prove that N (s)M -1 (s) = G(s).

For k ∈ {1, . . . , n}, we have

(N (s)M (s) -1 ) k = n l=1 N l (s)M inv lk (s) = k l=1 N l (s)M inv lk (s) • If I k = ∅, then M inv lk (s) = 0 for l = 1, . . . , k -1, and (N (s)M (s) -1 ) k = N k (s)M inv kk (s) = e -sh k R k (s α ). • If I k = ∅, then M inv lk (s) = 0 for l ∈ {1, . . . , k-1}\{k } where k = min I k and N k (s) = 0. Therefore, (N (s)M (s) -1 ) k = N k (s)M inv k k (s) = -N k (s) M k k (s) M k k (s)M kk (s)
.

Note that I k = ∅ since k = min I k . By replacing the above terms with appropriate expressions in ( 15), ( 16) and ( 17) and by noting that

σ k = σ k , we get (N (s)M (s) -1 ) k = e -sh k R k (s α ). It is obvious that N k (s), M kk (s), M k k (s), Y k (s) ∈ H ∞ .
Let us now prove that X(s) ∈ M(H ∞ ).

For k, k ∈ {1, . . . , n}, we have

X k k (s) = M inv k k (s) - Y k (s)e -sh k R k (s α ). • If I k = ∅, then Y k (s) = 0, and thus X k k (s) = M inv k k (s). Now, for k > k, M inv k k (s) = 0. For k = k, M inv k k (s) = 1. For k < k, from the fact that I k = ∅, we deduce that k = min I k , thus M inv k k (s) = 0. • If I k = ∅, then Y k (s) involves pk (s α ).
For k > k, the fact that

I k = ∅ leads to σ k = σ k . Therefore, Y k (s)R k (s α ) = a k β k p(s α ) (s α -σ k )(s α -σ k )u(s α ) belongs to H ∞ since (s α -σ k )(s α -σ k ) is eliminated by the same term in p(s α ). It is also obvious that M inv k k (s) = 0. Therefore, X k k (s) ∈ H ∞ . For k = k, we have X k k (s) = u(s α )(s α + 1) -β k pk (s α )a k e -sh k u(s α )(s α -σ k ) . If σ k < 0, then X k k (s) ∈ H ∞ for all β k . If σ k ≥ 0,
since β k given by ( 18) makes the numerator of

X k k (s) vanish at s = σ 1/α k , then X k k (s) ∈ H ∞ . For k < k, suppose first that k = min I k . From this fact together with I k = ∅ we deduce that σ k = σ k , and thus Y k (s)R k (s α ) ∈ H ∞ . We also have M inv k k (s) = 0, leading to X k k (s) ∈ H ∞ . In the case where k = min I k , we have σ k = σ k and X k k (s) = a k e -sh k u(s α )(s α + 1)e sh k -β k pk (s α )a k u(s α )(s α -σ k )a k .
By the same argument as in the case where k = k, we conclude that

X k k (s) ∈ H ∞ .
The right Bézout identity X(s)M (s) + Y (s)N (s) = I is clearly satisfied.

Remark 3 A transfer matrix G given by (1) with the delays of its elements not in order can be transformed to a transfer matrix G 0 with ordered delays by multiplying G by an appropriate permutation matrix P . It is well known that this matrix P is orthogonal and its inverse is P T . Assume that (M 0 , N 0 ) is an r.c.f. over H ∞ of G 0 and X 0 , Y 0 are the corresponding right Bézout factors. We have then G = G 0 P -1 = N 0 M -1 0 P -1 = N 0 (P M 0 ) -1 . It is obvious that P M 0 ∈ M(H ∞ ). Furthermore, X 0 P -1 P M 0 + Y 0 N 0 = I and X 0 P -1 ∈ M(H ∞ ). Hence, (P M 0 , N 0 ) is an r.c.f. of G and X 0 P -1 , Y 0 are the corresponding Bézout factors.

The next example illustrates the proposition.

Example 4

G(s) = e -s √ s , e -s √ s -1 , e -3s √ s -1
From ( 13), ( 14), ( 15), [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF], and ( 17), we obtain

N (s) = e -s √ s + 1 , e -s √ s + 1 , 0 , M (s) =     √ s √ s+1 0 0 0 √ s-1 √ s+1 -e -2s 0 0 1     .
The inverse of M (s) is

M -1 (s) =     √ s+1 √ s 0 0 0 √ s+1 √ s-1 e -2s ( √ s+1) √ s-1 0 0 1     .
The least common denominator of the entries of G(s)

is p( √ s) = √ s( √ s -1). Then p1 ( √ s) = √ s -1, and p2 ( √ s) = √ s. We choose u( √ s) = ( √ s + 1)
2 , which has no unstable roots. We have then

Y (s) = β 1 ( √ s -1) ( √ s + 1) 2 , β 2 √ s ( √ s + 1) 2 , 0 T , X(s) = M -1 (s) -Y (s)G(s) =     X 11 (s) X 12 (s) X 13 (s) X 21 (s) X 22 (s) X 23 (s) 0 0 1     ,
where

X 11 (s) = ( √ s + 1) 3 -β 1 ( √ s -1)e -s √ s( √ s + 1) 2 , X 12 (s) = - β 1 e -s ( √ s + 1) 2 , X 13 (s) = - β 1 e -3s ( √ s + 1) 2 , X 21 (s) = - β 2 e -s ( √ s + 1) 2 , X 22 (s) = ( √ s + 1) 3 -β 2 e -s √ s ( √ s -1)( √ s + 1) 2 , X 23 (s) = e -2s ( √ s + 1) 3 -β 2 e -3s √ s ( √ s -1)( √ s + 1) 2 .
We see that X 12 (s), X 13 (s), X 21 (s) ∈ H ∞ . From [START_REF] Si-Ammour | A sliding mode control for linear fractional systems with input and state delays[END_REF], we obtain β 1 = -1 and β 2 = 8e, which make X 11 (s), X 22 (s), and X 23 (s) be in H ∞ respectively.

Conclusion

In this paper, we have considered MISO fractional systems with input or output delays. Explicit expressions of an l.c.f. over H ∞ of the transfer matrices as well as the corresponding Bézout factors are given. Right coprime factorizations and right Bézout factors are also found for systems with entries of the transfer matrix containing different poles. In the case of identical poles, the right factors are primarily found for some simple classes of systems. Hence, in conclusion, we can have Youla-Kučera parametrization of stabilizing controllers for all systems with distinct unstable poles and a class of systems with identical unstable poles where each element of the transfer matrix involves one pole since for these systems both l.c.f.'s and r.c.f.'s and Bézout factors are available.

Even if similar techniques may in principle be used for MIMO systems, the analysis appears to be much more involved and requires a full separate study.

The proposed controllers in this paper given in terms of proper transfer functions containing only standard polynomials, fractional polynomials and delays are realizable. Their practical implementation is the objective of a forthcoming work.

A Proof of Proposition 2

First, we consider the condition (5) on the poles at zero.

If the system has no zero pole, then the condition is satisfied. Otherwise, the condition imposes that all powers of s whose order is smaller than m 0 α are eliminated in the development of the numerator of X(s) around zero. When s → 0, we have

(s α + 1) d u(s α ) = λ∈L(m0,α) a λ s λ + O(s m0α ), e -sh k q k (s α ) = λ∈L(m0,α) γ λk s λ + O(s m0α ), (A.1) µ k (s) = λ∈L(m0,α) β λk s λ + O(s m0α ).
Then the development of the numerator of X(s) around zero can be represented in matrix form as follows

(s α + 1) d u(s α ) - n k=1 e -sh k q k (s α )µ k (s) = [s λ N -1 , s λ N -2 , . . . , s λ0 ] A - n k=1 Γ k B k + O(s m0α ),
where N = card(L(m 0 , α)) and 0 = λ 0 < . . .

< λ N -1 ∈ L(m 0 , α), A = [a λ N -1 , . . . , a λ0 ] T , B k = [β λ N -1 k , .
. . , β λ0k ] T , and Γ k ∈ R N ×N are upper triangular matrices which contain γ λk and whose entries on the main diagonal are all γ λ0k . The condition ( 5) is then equivalent to

n k=1 Γ k B k = A.
For k = k 0 with k 0 defined by (3), since q k0 (s α ) does not have roots at zero, then (A.1) can be used and implies that γ λ0k0 = 0. Hence det Γ k0 = 0 and thus B k0 admits a unique solution for any values of β λk with λ < m 0 α and k ∈ {1, . . . , n}, k = k 0 .

Next, we analyze the condition (6) on non-zero poles.

We first examine the system of equations obtained by replacing s by a non-zero unstable pole b 1/α i and study the existence of µ

(l) ki (b 1/α i ),
where k i is defined by ( 4) and l = 0, . . . , m i -1, satisfying the equations. The first equation which corresponds to l = 0 contains q ki (b i )µ ki (b

1/α i ).
The second equation, i.e. l = 1, contains a linear sum of q ki (b i )µ ki (b

1/α i ) and q ki (b i )µ ki (b 1/α i
). Generally, the equation corresponding to the l-th derivative contains a linear sum of q ki (b i )µ ) with l 2 = 0, . . . , l -1. Since q ki (b i ) = 0, the system of m i equations can be recursively solved for m i unknowns µ ) for k = 1, . . . , n, k = k i and l = 0, . . . , m i -1 and admits a unique solution.

Hence, for each k ∈ {1, . . . , n} and k = k 0 , the coefficients β (jα)k , j = 0, . . . , f k -1 of µ k (s) satisfy µ (l) k (b 1/α i ) = a k,i,l for i = 1, . . . , N such that k i = k and l = 0, . . . , m i -1. This is the problem of Hermite interpolation and there exists a unique solution.

For k = k 0 , the coefficients of µ k0 (s), i.e. β λk0 with λ ∈ L(m 0 , α) and β (jα)k0 with j = m 0 , . . . , m 0 + f k0 -1, satisfy the equations µ (l) k0 (b 1/α i ) = a k0,i,l for i = 1, . . . , N such that k i = k 0 and l = 0, . . . , m i -1. We can write µ k0 (s) = ν k0 (s) + s m0α η k0 (s α ), where ν k0 (s) = λ∈L(m0,α) β λk0 s λ and η k0 (s α ) = m0+f k 0 -1 j=m0 β (jα)k0 s (j-m0)α . For arbitrary values of the coefficients β λk0 with λ ∈ L(m 0 , α), we can derive the values of η (l) k0 (b i ) for i = 1, . . . , N such that k i = k 0 and l = 0, . . . , m i -1. Note that the numbers of unknowns and of equations are the same and are equal to f k0 . This returns to the problem of Hermite interpolation and there exists a unique solution.

B Proof of Proposition 3

First, we consider the condition (9) on zero pole. For k ∈ {1, . . . , n} such that m 0k > 0, we develop the terms of the numerator of X kk (s) around zero as follows u(s α )(s α + 1)

  where σ m (•) is the smallest singular value of a matrix [20, Lemma 8.1.13 and Example 8.1.15].

  then the unstable poles of G k (s) are the unstable zeros of M kk (s). Now, from the right Bézout identity X(s)M (s) + Y (s)N (s) = I we obtain for k, k ∈ {1, . . . , n} and k = k

  ) satisfies the assumptions in Lemma 5. Let k ∈ {1, . . . , n}, if M k k (s) = 0 then due to (17) k = min I k , and M k k (s) = 0 for k = k . Also, since k = min I k , then k = min I k for k > k, and thus M kk (s) = 0. Hence the assumption (A) in Lemma 5 is satisfied.

  . . , m i -1 for arbitrary values of µ

  d k = λ∈L(m 0k ,α) a λ s λ + O(s m 0k α ), µ k (s) = λ∈L(m 0k ,α) β λk s λ + O(s m 0k α ), e -sh k q k (s α ) 1≤j≤n,j =k N i=1 (s α -b i ) mij = λ∈L(m 0k ,α) γ λk s λ + O(s m 0k α ). (B.1)Then the development of the numerator of X kk (s) when s → 0 can be written asu(s α )(s α + 1) d k -e -sh k µ k (s) q k (s α ) × 1≤j≤n,j =k N i=1 (s α -b i ) mij = [s λ N k -1 , s λ N k -2 , . . . , s λ0 ](A k -Γ k B k ) + O(s m 0k α ),whereN k = card(L(m 0k , α)), λ j ∈ L(m 0k , α), 0 = λ 0 < . . . < λ N k -1 ; A k = [a N k -1 , . . . , a 0 ] T ; B k = [β λ N k -1 k , .. . , β λ0k ] T ; and Γ k is an upper triangular matrix whose entries on the main diagonal are all γ λ0k . The condition (9) implies the elimination of the powers of s with order smaller than m 0k α in the numerator, which leads to Γ k B k = A k . From (B.1), we see that γ λ0k = 0, which implies Γ k is invertible and B k admits a unique solution. Now, we analyze the second condition (10) related to non-zero unstable poles. We can writeµ k (s) = ν k (s) + s m 0k α η k (s α ), where ν k (s) = λ∈L(m 0k ,α) β λk s λ and η k (s α ) = m 0k + N i=1m ik -1 j=m 0k β (jα)k s (j-m 0k )α . Then for each k ∈ {1, . . . , n}, we can derive the values of η

k

  (b i ) for i ∈ {1, . . . , N } such that m ik > 0 and l = 0, . . . , m ik -1. If we consider the coefficients β (jα)k , j = m 0k , . . . , m 0k + N i=1 m ik -1 as unknowns, then the number of unknowns is the same as the number of equations and is equal to N i=1 m ik . This is the problem of Hermite interpolation and there exists a unique solution.


Mercedes García-Martínez

Loïc Barrault

Fethi Bougares

Neural Machine Translation By Generating Multiple Linguistic Factors

Keywords: Machine Translation, Neural Networks, Deep Learning, Factored representation

Factored neural machine translation (FNMT) is founded on the idea of using the morphological and grammatical decomposition of the words (factors) at the output side of the neural network. This architecture addresses two well-known problems occurring in MT, namely the size of target language vocabulary and the number of unknown tokens produced in the translation. FNMT system is designed to manage larger vocabulary and reduce the training time (for systems with equivalent target language vocabulary size). Moreover, we can produce grammatically correct words that are not part of the vocabulary. FNMT model is evaluated on IWSLT'15 English to French task and compared to the baseline word-based and BPE-based NMT systems. Promising qualitative and quantitative results (in terms of BLEU and METEOR) are reported.

Introduction and related works

In contrast to the traditional phrased-based statistical machine translation [START_REF] Koehn | Moses: Open source toolkit for statistical machine translation[END_REF] that automatically translates subparts of the sentences, standard Neural Machine Translation (NMT) systems use the sequence to sequence approach at word level and consider the entire input sentence as a unit for translation [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF][START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF][START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF].

Recently, NMT showed better accuracy than existing phrase-based systems for several language pairs. Despite these positive results, NMT systems still face several challenges. These challenges include the high computational complexity of the softmax function which is linear to the target language vocabulary size (Equation 1).

p i = e oi / N r=1
e or for i ∈ {1, . . . , N }

where o i are the outputs, p i their softmax normalization and N the total number of outputs.

In order to solve this issue, a standard technique is to define a short-list limited to the s most frequent words where s << N . The major drawback of this technique is the growing rate of unknown tokens generated at the output. Another work around has been proposed in [START_REF] Jean | On using very large target vocabulary for neural machine translation[END_REF] by carefully organising the batches so that only a subset K of the target vocabulary is possibly generated at training time. This allows the system to train a model with much larger target vocabulary without substantially increasing the computational complexity. Another possibility is to define a structured output layer (SOUL) to handle the words not appearing in the shortlist. This allows the system to always apply the softmax normalization on a layer with reduced size [START_REF] Le | Large vocabulary SOUL neural network language models[END_REF]. The problem of unknown words was addressed making use of the alignments produced by an unsupervised aligner [START_REF] Luong | Addressing the rare word problem in neural machine translation[END_REF]. The unknown generated words are substituted in a post-process step by the translation of their corresponding aligned source word or copying the source word if no translation is found. The translation of the source word is made by means of a dictionary.

Other recent work have used subword units instead of words. In [START_REF] Sennrich | Neural machine translation of rare words with subword units[END_REF], some unknown and rare words are encoded as subword units with the Byte Pair Encoding (BPE) method. Authors show that this can also generates words unseen at training time. As an extreme case, the character-level neural machine translation has been presented in several works [START_REF] Chung | A character-level decoder without explicit segmentation for neural machine translation[END_REF][START_REF] Ling | Character-based neural machine translation[END_REF][START_REF] Costa-Jussà | Character-based neural machine translation[END_REF] and showed very promising results. The character-level NMT architectures are composed of many layers, to deal with the long distance dependencies, increasing aggressively the computational complexity of the training process. In [START_REF] Sennrich | How grammatical is character-level neural machine translation? assessing mt quality with contrastive translation pairs[END_REF] has been shown that character-level decoders outperform subwords units using BPE method when processing unknown words, but they perform worse when extracting morphosyntactic information about the sentences, due to the long distances.

Among other previous works, our work can be seen as a continuation of [START_REF] García-Martínez | Factored neural machine translation architectures[END_REF]. Several works have used factors as additional information for the input words in neural language modelling with interesting results [START_REF] Niehues | Using factored word representation in neural network language models[END_REF][START_REF] Alexandrescu | Factored neural language models[END_REF][START_REF] Wu | Factored recurrent neural network language model in ted lecture transcription[END_REF]. More recently, factors have also been integrated into a word-level NMT system as additional linguistic input features [START_REF] Sennrich | Linguistic input features improve neural machine translation[END_REF]. Unlike these previous works, we are considering factors as translation unit. We refer to factors as some linguistic annotations at word level, e.g. the Part of Speech (POS) tag, number, gender, etc. The advantages of using factors as translation unit are two-fold: reducing the output vocabulary size and allowing to generate surface forms which are never seen in the training data.

Factors were first introduced for NMT at output side in [START_REF] García-Martínez | Factored neural machine translation architectures[END_REF] where two factored synchronous symbols are simultaneously generated. Authors presented an investigation of the architecture of their factored NMT system to show that better results are obtained using a feedback of the two generated outputs concatenation.

Our work is different from previous efforts in that we consider only the best type of feedback for the network. We also introduce an additional factor about the case information (lowercase, uppercase or in capitals) and evaluate using a different translation test. Moreover, we apply an unknown words (unk) replacement technique using the alignments of the attention mechanism to replace the generated unknown words in target side. For that, we make use of an unigram dictionary to find the translation of the source word corresponding to the generated unk.

We compare this architecture to the state of the art BPE approach and the classic word-level NMT approach on the English to French dataset from IWSLT'15 evaluation campaign. We provide, in addition a quantitative and qualitative study about the obtained results.

The remainder of this paper is organized as follows: Section 2 describes the attentionbased NMT system and Section 3 its extension using the factored approach. In Section 4, we describe the experiments and the obtained results. Finally, Section 5 concludes the paper and presents the future work.

Neural Machine Translation

The standard NMT model consists of a sequence to sequence encoder-decoder of two recurrent neural networks (RNN), one used by the encoder and the other by the decoder. The source language sequence is mapped into an embedded dimension in the encoder and the decoder maps the representation back to a target language sequence.

DECODER ENCODER Attention Mechanism x 1 x 2 x N • • • x 1 x 2 x N • • • • • • Annotations Bidirectional RNN GRU1 GRU2 X ↵ijai a 1 a N C j L C L R L O Softmax ⇣X⌘ Word Fig. 1. Attention-based NMT system
The architecture includes a bidirectional RNN encoder (see left part of Figure 1) equipped with an attention mechanism [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF]. Each input sentence word x i (i ∈ 1 . . . N with N the source sequence length) is encoded into an annotation a i by concatenating the hidden states of a forward and a backward RNN provided by a gated recurrent unit (GRU) [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF] to control the flow of information. These annotations a 1 . . . a N represents the whole sentence with a focus on the word being processed. One difference from the architecture of [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF] is that the decoder contains a conditional GRU [START_REF] Firat | Conditional gated recurrent unit with attention mechanism[END_REF] which consists of two GRUs interspersed with the attention mechanism (see right top part of the Figure 1). The first GRU combines the embedding of the previous decoded token and the previous hidden state in order to generate an intermediate representation which is an input of the attention mechanism and the second GRU. The attention mechanism (bottom yellow part of the Figure 1) computes a source context vector C j as a convex combination of annotation vectors, where the weights of each annotation are computed locally using a feed-forward network. These weights can be used to align the target words with the source positions. The second GRU generates the hidden state of the conditional GRU by looking at the output of the first GRU and the context vector C j . The decoder RNN takes as input the embedding of the previous output word (feedback of the network) in the first GRU, the context vector C j in the second GRU and its hidden state. The output layer L O is connected to the network through a hyperbolic tangent sum operation Φ( ) which takes as input the embedding of the previous output word as well as the context vector and the output of the decoder from the second GRU (both adapted with a linear transformation, respectively, L C and L R ). Finally, the output probabilities for each word in the target vocabulary are computed with a softmax function. The word with the highest probability is the translation output at each timestep. The encoder and the decoder are trained jointly to maximize the conditional probability of the reference translation.

Factored Neural Machine Translation

The Factored Neural Machine Translation (FNMT) [START_REF] García-Martínez | Factored neural machine translation architectures[END_REF] is an extension of the standard NMT architecture which allows the system to generate several output symbols at the same time. For the sake of simplicity, only two symbols are generated: the lemma and the concatenation of the different factors (verb, tense, person, gender, number and case information). The target words are then represented by a factored output: lemmas and factors. Factors may help the translation process providing grammatical information to enrich the output. The task of this work is English to French translation, English is a grammatically poor language and factors do not help for its translation, this has been tested in previous experiments. Therefore, we apply the factors only in the target side when translating to French which is a grammatically rich language. In the example shown in Figure 3, from the verbal form in French devient, we obtain the lemma devenir and its factors VP3#SL (Verb, in Present, 3rd person, no gender (#), Singular and Lowercased form). Moreover, we can see the word interéssant with the lemma interéssant and factors Adj##MSL (Adjective, no tense (#) and no person (#), Masculine gender, Singular number and Lowercased form). The morphological analyser MACAON toolkit [START_REF] Nasr | Macaon, an nlp tool suite for processing word lattices[END_REF] is used to obtain the lemma and factors for each word taking into account its context with nearly 100% accuracy. The first entry is used in the few cases that MACAON proposes multiple words (e.g. same word written in two forms).

The FNMT architecture is presented in Figure 2. The encoder and attention mechanism of Figure 1 remain unchanged. However, the decoder has been modified to get multiple outputs. The hidden state of the conditional GRU (cGRU) is shared to produce simultaneously several outputs. The output from the layer L O has been diversified to two softmax layers, one to generate the lemma and the other to generate the factors. An additional design decision is related to the decoder feedback. Contrary to the word based model, where the feedback is naturally the previous word (see Figure 1), we have multiple choices where multiple outputs are generated for each decoding time-step. We have decided to use the concatenation of the embeddings of both generated symbols based on the work [START_REF] García-Martínez | Factored neural machine translation architectures[END_REF].

The FNMT model may lead to sequences with a different length, since lemmas and factors are generated simultaneously but separately (each sequence ends after the generation of the end of sequence < eos > token). To avoid this, the sequences length is decided based on the lemmas stream length (i.e the length of the factors sequence is constrained to be equal to the length of the lemma sequence). This is motivated by the fact that the lemmas contain most of the information of the final surface form (word).

Once we obtain the factored outputs from the neural network, we need to combine them to obtain the surface form (word representation). This operation is also performed with the MACAON tool, which, given a lemma and some factors, provides the word. Word forms given by MACAON toolkit have a 99% success rate. In the cases (e.g. name entities) that the word corresponding to the lemma and factors is not found, the system outputs the lemma itself.

Experiments

We performed a set of experiments for Factored NMT (FNMT) and compared them with the word-based NMT and BPE-based NMT systems.

Data processing and selection

The systems are trained on the English to French (EN-FR) Spoken Language Translation task from IWSLT 2015 evaluation campaign 1 . We applied data selection using modified Moore-Lewis filtered by XenC [START_REF] Rousseau | XenC: An open-source tool for data selection in natural language processing[END_REF] to obtain a sub part of the available parallel corpora (news-commentary, united-nations, europarl, wikipedia, and two crawled corpora). The Technology Entertainment Design (TED) [START_REF] Cettolo | Wit 3 : Web inventory of transcribed and translated talks[END_REF] corpus has been used as in-domain corpus.

We preprocess the data to convert html entities and filter out the sentences with more than 50 words for both source and target languages. Finally, we obtain a corpus of 2M sentences with 147k unique words for the English side and 266k unique words for the French side. French vocabulary is bigger than English since French is more highly inflected language. Table 1 shows training, development and testing sets statistics.

Training

Models are trained using NMTPY [START_REF] Caglayan | Nmtpy: A flexible toolkit for advanced neural machine translation systems[END_REF], an NMT toolkit in Python based on Theano2 . The following hyperparameters have been chosen to train the systems. The embedding and recurrent layers have the dimensions 620 and 1000, respectively. The batch size is set to 80 sentences and the parameters are trained using the Adadelta [START_REF] Zeiler | Adadelta: an adaptive learning rate method[END_REF] optimizer. We clipped the norm of the gradient to be no more than 1 [START_REF] Pascanu | Understanding the exploding gradient problem[END_REF] and initialize the weights using Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. The systems are validated on dev15 dataset using early stopping based on BLEU [START_REF] Papineni | Bleu: a method for automatic evaluation of machine translation[END_REF]. The vocabulary size of the source language is set to 30K. The output layer size of the baseline NMT system is set to 30K. For the sake of comparability and consistency, the same value (30k) is used for the lemma output of the FNMT system. This 30K FNMT vocabulary includes 17k lemmas obtained from the original NMT vocabulary (30k word level gives 17k lemmas when all the derived forms of the verbs, nouns, adjectives, etc are discarded) increased with additional new lemmas to fit the 30K desired value. The factors have 142 different units in their vocabulary. When it comes to combining the lemmas and the factors vocabulary, the system is able to generate 172K different words, using the external linguistic resources, which is 6 times bigger than a standard word-based NMT vocabulary.

For BPE systems, bilingual vocabulary has been built using source and target language applying the joint vocabulary BPE approach. In order to create comparable BPE systems, we set the number of merge operations for the BPE algorithm (the only hyperparameter of the method) as 30K minus the number of character according to the paper [START_REF] Sennrich | Neural machine translation of rare words with subword units[END_REF]. Then, we apply a total of 29388 merge operations to learn the BPE models on the training and validation sets. During the decoding process, we use a beam size of 12 as used in [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF].

Quantitative results

The Factored NMT system aims at integrating linguistic knowledge into the decoder in order to overcome the restriction of having a large vocabulary at target side. We first compare our system with the standard word-level NMT system. For the sake of comparison with state of the art systems, we have built a subword system using the BPE method. Subwords were calculated at the input and the output side of the neural network as described in [START_REF] Sennrich | Neural machine translation of rare words with subword units[END_REF]. The results are measured with two automatic metrics, the most common metric for machine translation BLEU and METEOR [START_REF] Lavie | Meteor: an automatic metric for mt evaluation with high levels of correlation with human judgments[END_REF]. We evaluate on test15 dataset from the IWSLT 2015 campaign and results are presented in Table 2. As we can see from the Table 2 results, the FNMT system obtains better %BLEU and %METEOR scores compared to the state of the art NMT and BPE systems. An improvement of about 1 %BLEU point is achieved compared to the best baseline system (BPE). This improvement is even bigger (1.4 %BLEU point) when UNK replacement is applied to both systems. In a quest to better understand the reasons of this improvement, we also computed the %BLEU scores of each output level (lemmas and factors) for FNMT. Theses scores are presented in Table 2. The lemma and factors scores of NMT and BPE systems are obtained through a decomposing of their word level output into lemma and factors. We observe yet again that FNMT systems gives better score at both lemma and factors level. Replacement of unknown words has been performed using the alignments extracted from the attention mechanism. We have replaced the generated UNK tokens by translating its highest probability aligned source word. We see an improvement of around 1 point %BLEU score in both NMT and FNMT systems.

The last column of Table 2 shows, for each system, the the number of generated UNK tokens. As shown in the table our FNMT system produces half of the UNK tokens compared to the word-based NMT system. This tends to prove that the Factored NMT system effectively succeed in modelling more words compared to the word based NMT system augmenting the generalization power of our model and preserving manageable output layer sizes. Though we can see that BPE system does not produce UNK tokens, this is not reflected in the scores. Indeed, this can be due to the possibility of generation of incorrect words using BPE units in contrast to the FNMT system.

Qualitative analysis

The strengths of FNMT are considered under this qualitative analysis. We have studied and compared the translation outputs of NMT at word-level and BPE-level with the ones of FNMT systems. Two examples are presented in Figure 3 and Table 3. The reference translation of the source sentence presented in Figure 3 is "mais voilà où c ¸a devient intéressant". As we can see, contrary to the baseline NMT system, the FNMT system matches exactly the reference and thus produces the correct translation. An additional interesting observation is that the alignment provided by the attention mechanism seems to be better defined and more helpful when using factors. Also, one can notice the difference between the attention distributions made by the systems over the source sentence. The NMT system first translated "here" into "là", added a coma, and then was in trouble for translating the rest of the sentence, which is reflected by the rather fuzzy attention weights. The FNMT system had better attention distribution over of the source sentence in this case.

-p-l pct-l prep-l nc-f-s-l pct-l cln-1-s-l v-P-1-s-l pct-l v-P-3-p-l vppart-K-m-p-l
Table 3 shows another example comparing NMT, BPE and FNMT systems. The NMT system generated an unknown token (UNK) when translating the English word "baffled". We observe that BPE translates "baffled" to "bafés" which does not exist in French. This error probably comes from the shared vocabulary between the source and target languages creating an incorrect word very similar to its aligned source tokens. FNMT translates it to "déconcertés" which is a better translation than in the reference. One should note that it is not generated by the unknown word replacement method. However, for this particular example, an error on the factors leads to the word "sont" instead of "sommes", resulting in lower automatic scores for FNMT output.

Conclusion

In this paper, the Factored NMT approach has been further explored. Factors based on linguistic a priori knowledge have been used to decompose the target words. This approach outperforms a strong baseline system using subword units computed with byte pair encoding. Our FNMT system is able to model an almost 6 times bigger word vocabulary with only a slight increase of the computational cost. By these means, the FNMT system is able to halve the generation of unknown tokens compared to wordlevel NMT. Using a simple unknown word replacement procedure involving a bilingual dictionary, we are able to obtain even better results (+0.8 %BLEU compared to previous best system). Also, the use of external linguistic resources allows us to generate new word forms that would not be included in the standard NMT system shortlist. The advantage of this approach is that the new generated words are controlled by the linguistic knowledge, that avoid producing incorrect words, as opposed to actual systems using BPE. We demonstrated the performance of such a system on an inflected language (French). The results are very promising for use with highly inflected languages like Arabic or Czech.

Fig. 2 .

 2 Fig. 2. Detailed view of the decoder of the Factored NMT system

Fig. 3 .

 3 Fig. 3. Examples of NMT and FNMT outputs aligned against the source sentence

Table 1 .

 1 Datasets statistics

	Data	Corpus name	Datasets	# Sents # Words EN-FR
	Training	train15	data selection	2M	147k-266k
	Development	dev15	dev10+test10+test13 3.6k	7.3k-8.9k
	Testing	test15	test11+test12	1.9k	4.5k-5.4k

Table 2 .

 2 Results on IWSLT test15. %BLEU and %METEOR performance of NMT and FNMT systems with and without UNK replacement (UR) are presented. For each system we provide the number of generated UNK tokens in the last column

		%METEOR↑		%BLEU↑
	Model	word	word	lemma factors #UNK
	NMT / +UR 62.21 / 63.38 41.80 / 42.74 45.10 51.80 1111
	BPE	62.87	42.37	45.96 53.31	0
	FNMT / +UR 64.10 / 64.81 43.42 / 44.15 47.18 54.24 604

Table 3 .

 3 Examples of translations with NMT, BPE and FNMT systems (without unknown words replacement)

	Src			we in medicine , I think , are baffled
	Ref			Je pense que en médecine nous sommes dépassés
	NMT	Nous	,	en médecine ,	je	pense	,	sont	UNK
	BPE	nous	,	en médecine ,	je	pense	, sommes b@@ af@@ és
	FNMT	nous	,	en médecine ,	je	pense	,	sont	déconcertés
	Lemmas lui	,	en médecine ,	je	penser ,	être	déconcerter
	Factors pro-1

https://sites.google.com/site/iwsltevaluation2015

https://github.com/lium-lst/nmtpy

Acknowledgments

This work was partially funded by the French National Research Agency (ANR) through the CHIST-ERA M2CR project, under the contract number ANR-15-CHR2-0006-01.