text:
I Pontes Duff

C Poussot-Vassal

C Seren

Optimal H 2 model approximation based on multiple input/output delays systems

Keywords: Model reduction, time-delay systems, large-scale systems, linear systems. 1

 for the delay-free case, which is the main contribution of this paper. Secondly, a two stage algorithm is proposed in order to practically obtain such an approximation.

Introduction

Model approximation plays a pivotal role in many simulation based optimization, control, analysis procedures. Indeed, due to memory and computational burden limitations working with a reduced order model in place of the original one, potentially large-scale, might be a real advantage. To this aim, most of the results presented in the literature address the linear dynamical systems approximation problem in the delay-free case 1 . More specifically, this problem has been widely studied using either Lyapunov-based methods [START_REF] Spanos | A new algorithm for l 2 optimal model reduction[END_REF][START_REF] Hyland | The optimal projection equations for model reduction and the relationships among the methods of Wilson, Skelton, and Moore[END_REF][START_REF] Wilson | Optimum solution of model-reduction problem[END_REF], interpolation-based algorithm [START_REF] Meier | Approximation of linear constant systems[END_REF][START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF][START_REF] Van Dooren | H 2 -optimal model reduction of MIMO systems[END_REF][START_REF] Beattie | A trust region method for optimal H 2 model reduction[END_REF], or matching moments approaches [START_REF] Grimme | Krylov projection methods for model reduction[END_REF][START_REF] Astolfi | Model reduction by moment matching for linear and nonlinear systems[END_REF], leading to a variety of solutions and applications. Recent surveys are available in [START_REF] Antoulas | A survey of model reduction methods for large-scale systems[END_REF][START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF][START_REF] Opmeer | Optimal model reduction for non-rational functions[END_REF]. The presence of input/output delays in the approximation model was tackled in [START_REF] Halevi | Reduced-order models with delay[END_REF] (exploiting both Lyapunov equations and grammians properties derived in [START_REF] Hyland | The optimal projection equations for model reduction and the relationships among the methods of Wilson, Skelton, and Moore[END_REF] for the free-delay case). The bottleneck of this approach is that it requires to solve Lyapunov equations which might be costly in the large-scale context. From the moment matching side, [START_REF] Scarciotti | Model reduction by moment matching for linear time-delay systems[END_REF] proposed a problem formulation that enables the construction of an approximation which contains very rich delay structure (including state delay), but where the delays and the interpolation points are supposed to be a priori known. From the Loewner framework side, [START_REF] Duff | Realization independent single time-delay dynamical model interpolation and H 2 -optimal approximation[END_REF] and after [START_REF] Schulze | Data-driven interpolation of dynamical systems with delay[END_REF] generalizes the Loewner framework from [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF] to the state delay case enabling data-driven interpolation. However, as for the moment matching case, the delays and the interpolation points are supposed to be a priori known.

In this paper, the problem of approximating a given large-scale model by a low order one including (a priori unknown) I/O delays using the interpolatory framework, is addressed. An alternative "poles/residues"-based approach is developed, which enables to reach the H 2 optimality conditions, treated as interpolation ones. Then, the main contribution of this paper consists in extending the interpolation results of [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF] to the case of approximate models with an extended structure, namely, including non-zero input(s)/output(s) delays. Last but not least, H 2 optimality conditions for such cases are also elegantly derived together with a single numerical procedure.

The paper is organized as follows: after introducing the notations and the mathematical problem statement in Section 2, Section 3 recalls some necessary preliminary results related to the computational aspects of the H 2 inner product and H 2 norm when the calculations are based on the poles/residues decomposition of a transfer function. Section 4 establishes the H 2 optimality conditions solving the input/output delay dynamical model approximation problem. It also proposes an algorithm which permits to practically compute such an approximation. Section 5 details the results obtained after treating an academic example. Conclusions and prospects end this article in Section 6.

Notations and problem statement

Notations. Let us consider a stable Multiple-Input/Multiple-Output (MIMO) linear dynamical system, denoted by G in the sequel, with n u (resp. n y ) ∈ N * input(s) (resp. output(s)), represented by its transfer function G(s) ∈ C n y ×n u . Let H n y ×n u 2 be the Hilbert space of holomorphic functions F : C → C n y ×n u which are analytic in the open right-half plane and for which

+∞ -∞ trace F(iω)F T (iω) dω<+∞. For given G, H ∈ H n y ×n u 2
, the associated inner-product reads:

G, H H 2 = 1 2π +∞ -∞ trace G(iω)H T (iω) dω, (1)
and the H n y ×n u 2 induced norm can be explained:

G H 2 = 1 2π +∞ -∞ G(iω) 2 F dω 1/2 = G, G H 2 , (2)
where G 2 F = G, G F and G, H F = trace(GH T ) are the Frobenius norm and inner-product, respectively. Dynamical system H will be said real iff. ∀s ∈ C,

H(s) = H(s). It is noteworthy that if G(s), H(s) ∈ H n y ×n u 2 are real, then G, H H 2 = H, G H 2 ∈ R + .
Besides, any dynamical matrix ∆(s) will belong to H n y ×n u ∞ iff. sup{σ max (∆(iω))/ω ∈ R} < +∞. σ max (∆(iω)) refers to the largest singular value of matrix ∆(iω).

Followingly, let Ĥd be a multiple-input/output delays MIMO system s.t. Ĥd (s) ∈ H n y ×n u 2 and represented by: Ĥd :

Êẋ (t) = Âx(t) + B∆ i (u(t)) ŷ(t) = ∆ o ( Ĉx(t)) , (3)
where Ê, Â ∈ R n×n (with state dimension n ∈ N * ), B ∈ R n×n u , Ĉ ∈ R n y ×n and ∆ i and ∆ o are delay operators. The matrix transfer functions ∆i (s) and ∆o (s) defined in [START_REF] Wilson | Optimum solution of model-reduction problem[END_REF] represent the frequency behavior of the delays operators ∆ i and ∆ o , receptively. The transfer function of the underlying system (3) from input û(t) to output ŷ(t) vectors is given by:

Ĥd (s) = ∆o (s) Ĥ(s) ∆i (s) ∈ H n y ×n u 2 , (4)
where:

           Ĥ(s) = Ĉ( Ês -Â) -1 B ∈ H n y ×n u 2 ∆i (s) = diag(e -sτ 1 . . . e -sτ nu ) ∈ H n u ×n u ∞ ∆o (s) = diag(e -sγ 1 . . . e -sγ ny ) ∈ H n y ×n y ∞ . (5)
From this point, we will denote by Ĥd = ( Ê, Â, B, Ĉ, ∆i , ∆o ) a MIMO input/output delayed system of the form (4). Ĥd will also be said to have order n N (where N is the original model order).

Problem statement.

dim( Ĥd ) ≤ n G -Ĥd H 2 ,
where Ĥd = ∆o Ĥ ∆i as in, (

This search for an optimal solution will be carried out assuming that both G and Ĥ from Eq. ( 5) have semi-simple poles i.e., s.t. their respective transfer function matrix can be decomposed as follows:

G(s) = N j=1 l j r T j s -µ j and Ĥ(s) = n k=1 ĉk bT k s -λk , (6)
where ∀ j = 1 . . . N, ∀k = 1 . . . n, r j , bk ∈ C n u ×1 and l j , ĉk ∈ C n y ×1 . The poles µ j , λk are elements of C -so that G and Ĥ belong to H n y ×n u 2 .

Preliminary results

In this section, some elementary but important, results, which will be useful along this paper, are recalled and generalized.

First of all, a fundamental result dealing with the H 2 norm invariance in case of input/output delayed systems is presented.

Proposition 3.1. (H 2 norm invariance) Let Ĥ ∈ H n y ×n u 2
be a stable dynamical system and

M ∈ H n u ×n u ∞ , N ∈ H n y ×n y ∞ s.t.: ∀ω ∈ R, M(iω)M(iω) T = I n u , N(iω) T N(iω) = I n y . ( 7
)
If Ĥd = N ĤM then Ĥd H 2 = Ĥ H 2 .
Proof. If Ĥd = N ĤM, the scaled term 2π Ĥd 2 H 2 will then read by definition:

+∞ -∞ trace N(iω) Ĥ(iω)M(iω)M T (iω) ĤT (iω)N T (iω) dω = +∞ -∞ trace N(iω) Ĥ(iω) ĤT (iω)N T (iω) dω = +∞ -∞ trace Ĥ(iω) ĤT (iω)N T (iω)N(iω) dω = +∞ -∞ trace Ĥ(iω) Ĥ(iω) T dω = 2π Ĥ 2 H 2 .
One can easily check that condition [START_REF] Beattie | A trust region method for optimal H 2 model reduction[END_REF] appearing in Proposition 3.1 is satisfied by the delays matrices of the two last lines of ( 5) when M = ∆i and N = ∆o . In other words, the H 2 norm does not depend on the input, nor output delays. The following proposition makes now explicit the calculation of the H 2 norm associated with the dynamical mismatch gap G -Ĥd , which conditions Problem 2.1 criterion.

Proposition 3.2. Let G, Ĥd ∈ H n y ×n u 2 s.t.
Ĥd is given by Eq. ( 4). The H 2 norm of the approximation gap (or mismatch error), denoted by J, can be expressed as:

J = G -∆o Ĥ ∆i 2 H 2 = G 2 H 2 -2 G, ∆o Ĥ ∆i H 2 + Ĥ 2 H 2 . (8)
Proof. Simply develop the H 2 norm using the inner product definition and exploit the previous result ∆o Ĥ ∆i

H 2 = Ĥ H 2 .
Obviously, regarding Eq. ( 8), minimizing J is equivalent to minimize -2 G, ∆o Ĥ ∆i

H 2 + Ĥ 2
H 2 and thus to look for the optimal values of the decision variables contained in both the realization Ĥ ∈ H n y ×n u 2 and the delay blocks ∆i , ∆o ∈ H n y ×n u ∞ . At this point, it could be profitable to derive suitable analytical expressions for the inner-product and the H 2 norm of Ĥ in order to define more precisely the aforementioned H 2 gap between the two transfer functions. To this aim, the previous assumption made for both G and Ĥ systems (see Eq. ( 6)) will be essential to obtain the following results. Proposition 3.3. (H 2 inner product computation with input/output delays) Let G, Ĥ be two systems ∈ H n y ×n u 2 whose respective transfer functions G(s) and Ĥ(s) can be expressed as in [START_REF] Meier | Approximation of linear constant systems[END_REF]. Let ∆i , ∆o be real, H n u ×n u ∞ and H n y ×n y ∞ respectively, models satisfying sup{ ∆o (s), ∆i (s) /s ∈ C -} = M < +∞. By denoting Ĥd = ∆o Ĥ ∆i , the inner product Ĥd , G H 2 reads:

Ĥd , G H 2 = N j=1 trace Res Ĥd (-s)G T (s), µ j = N j=1 l T j ∆o (-µ j ) Ĥ(-µ j ) ∆i (-µ j )r j . (9)
with:

Γ I = {s ∈ C/s = iω and ω ∈ [-R; R], R ∈ R + } Γ R = {s ∈ C/s = Re iθ where θ ∈ [π/2; 3π/2]} .
Thus, for a sufficient large radius value R, the Γ C contour will contain all the poles of the transfer function G(s) i.e., µ 1 , µ 2 , . . . , µ N . Thus, by applying the residues theorem, it follows that:

Ĥd , G H 2 = 1 2π +∞ -∞ trace Ĥd (iω)G T (iω) dω = lim R→+∞ 1 2iπ Γ C Ĥd (-s)G(s)ds = N j=1 trace Res Ĥd (-s)G T (s), µ j .
where Res(.) denotes the residue operator. The second equality line holds true since:

Γ R Ĥd (-s)G(s)ds ≤ M 2 Γ R Ĥ(-s)G(s)ds → 0 + , when R → +∞.
One may note that Proposition 3.3 is a generalization of Lemma 3.5 appearing in [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF] in the case of MIMO systems with multiple-input/output delays. It is noteworthy that the ∆i , ∆o matrices defined by [START_REF] Wilson | Optimum solution of model-reduction problem[END_REF] clearly verifies the hypothesis Proposition 3.3. Remark 3.1. (Delay-free case "symmetry") An equivalent proposition was derived in the delay-free case [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF]. It can be recovered from Proposition 3.3 by taking ∆i = I n u and ∆o = I n y . The result corresponds to the symmetric expression of the inner product i.e., the evaluation of G in the poles of Ĥ and its associated residues ĉk and bk s.t.:

G, Ĥ H 2 = n k=1 ĉT k Ĝ(-λk ) bk = N j=1 l T j Ĥ(-µ j )r j = Ĥ, G H 2 .
In the presence of input/output delays, since the H 2 norm cannot be approximated using one contour containing the poles of Ĥd only, this result is no longer true. Indeed, it can be easily shown that in this case, the integral on Γ R will depend on a positive exponential argument which will not converge to 0 + when R → +∞. This justifies the assumption that sup{ ∆o (s), ∆i (s) /s ∈ C -} = M < +∞ and relevance of Proposition 3.3.

Finally, let us recall the pole(s)/residue(s) H 2 norm formula.

Corollary 3.1. (Poles/residues H 2 norm [1]) Assume that Ĥd (s), Ĥ(s) belong to H n y ×n u 2
and that Ĥd = ∆o Ĥ ∆i . Besides, suppose that Ĥ can be expressed such as in (6), then,

Ĥd 2 H 2 = n k=1 ĉT k Ĥ(-λk ) bk .
Proof. See [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF].

In the next section, the main result, namely H 2 optimality conditions related to Problem 2.1, are firstly established and an interpolation-based algorithm is proposed to numerically compute the approximation Ĥd .

Approximation by multiple I/O delays MIMO systems: H 2 optimality conditions

Considering the mathematical formulation of Problem 2.1 and the reduced order system structure Ĥd = ∆o Ĥ ∆i , where Ĥ(s) is given as in [START_REF] Meier | Approximation of linear constant systems[END_REF], the underlying optimization issue that must be solved is parameterized by (k = 1, . . . , n): (i) the n pole(s) λk ∈ C -; (ii) the n bi-tangential directions ( bk , ĉk ) ∈ C n u ×1 × C n y ×1 ; and (iii) the n u + n y delay values (τ l , γm ), l = 1 . . . n u , m = 1 . . . n y . Our primary objective consists in rewriting the expression of the H 2 gap J as a function of these latter parameters which will subsequently facilitate the derivation of the H 2 optimality conditions for Problem 2.1. This forms the topic of the three following propositions and of Theorem 4.1, which stands as the main result of the paper. Proposition 4.1. From the preliminary results, the mismatch H 2 gap defined previously in Proposition 3.2 can be equivalently rewritten as:

J = G 2 H 2 + n k=1 ĉT k Ĥ(-λk ) bk . . . -2 N j=1 l T j ∆o (-µ j ) Ĥ(-µ j ) ∆i (-µ j )r j . (10)
Proof. The result is immediate. To be established, it requires to develop the H 2 norm expression showing the inner product and then to use both Proposition 3.3 and Corollary 3.1 results.

From the previous equation ( 10), the first-order optimality conditions related to the minimization of J can be analytically computed. The gradient expressions of the H 2 gap w.r.t. each parameters (delays, tangential directions and poles) are detailed in the two following propositions. Starting with the simplest calculations, we first derive the gradient of J w.r.t. the delays since the second term of the right-hand side part of ( 10) is delay-dependent, only.

Proposition 4.2. The gradients of the H 2 gap J with respect to the delays read ∀l = 1 . . . n u , ∀m = 1 . . . n y :

                                           ∇ τl J = -2 ∂ Ĥd , G H 2 ∂τ l = -2 N j=1 µ j l T j ∆o (-µ j ) Ĥ(-µ j )D l ∆i (-µ j )r j , ∇ γm J = -2 ∂ Ĥd , G H 2 ∂γ m = -2 N j=1 µ j l T j D m ∆o (-µ j ) Ĥ(-µ j ) ∆i (-µ j )r j ,
where elements of D l ∈ R n u ×n u , D m ∈ R n y ×n y , are defined as:

[D k ] i j = δ i jk = 1 if i = j = k 0 otherwise .
Proof. The proof is straightforward to establish since both ∆i and ∆o terms are diagonal matrices and the exponential derivative function is obvious.

Proposition 4.3. The gradients of the H 2 gap J with respect to parameters ĉk , bk and λk , ∀k = 1 . . . n read:

                         ∇ ĉk J = -2 ∂ Ĥd , G H 2 ∂ĉ k + ∂ Ĥ 2 H 2 ∂ĉ k = -2 bT k G(-λk ) -Ĥ(-λk ) T , ∇ bk J = -2ĉ T k G(-λk ) -Ĥ(-λk ) , ∇ λk J = 2ĉ T k G (-λk ) -Ĥ (-λk ) bk ,
where:

G(s) = N j=1 ∆o (-µ j ) l T j r j s -µ j ∆i (-µ j ). ( 11
)
and where G and Ĥ are the Laplace derivative of G and Ĥ, respectively.

Proof. By defining r j = ∆i (-µ j )r j and lT j = l T j ∆o (-µ j ) with j = 1 . . . N, the H 2 gap can be written as:

J = G 2 H 2 -2 N j=1 lT j n m=1 ĉm bT m -µ j -λm r j + n k=1 ĉT k n m=1 ĉm bT m -λk -λm bk .
Then, calculating the gradients w.r.t. bl , ĉl and λl (l = 1 . . . n) gives:

∇ bl J = -2 ∂ Ĥd , G H 2 ∂ bl + ∂ Ĥ 2 H 2 ∂ bl
Thus, by computing both terms on this expression

∂ Ĥ 2 H 2 ∂ bl = n k=1 n m=1 (ĉ T k ĉm ) -λk -λm ∂ ∂ bl bT m bk = 2 n k=1 ĉT l ĉk bT k -λk -λl = 2ĉ T l Ĥ(-λl ) and ∂ Ĥd , G H 2 ∂ bl = N j=1 n m=1 ( ˜lT j ĉm )r T j -µ j -λm ∇ bl bm = ĉT l N j=1 ˜l j rT j -µ j -λl = ĉT l G(-λl ).
one obtains the gradient. It is noteworthy that ∇ ĉl J can be obtained in the same way as ∇ bl J. The calculation of ∇ λl J is straightforwardly derived as follows:

∇ λl J = -2 N j=1
lT j ĉl bT l r j (-λlµ j ) 2 -ĉT l Ĥ (-λl ) bl . . .

+ n k=1 ĉT k ĉl bT l bk (-λl -λk ) 2 = 2ĉ T l G (-λl ) -Ĥ (-λl ) bl .
           Ĥ(-λk ) bk = G(-λk ) bk , ĉT k Ĥ(-λk ) = ĉT k G(-λk ), ĉT k Ĥ (-λk ) bk = ĉT k G (-λk ) bk , (12)
                     N j=1 µ j l T j ∆o (-µ j ) Ĥ(-µ j )D l ∆i (-µ j )r j = 0, N j=1 µ j l T k D m ∆o (-µ j ) Ĥ(-µ j ) ∆i (-µ j )r j = 0, (13)
for all k = 1 . . . n, l = 1 . . . n u and m = 1 . . . n y where G(s) is given by [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF].

Proof. The interpolation conditions gathered in [START_REF] Opmeer | Optimal model reduction for non-rational functions[END_REF] are deduced by taking ∇ ĉl J = 0, ∇ bl J = 0 and ∇ λl J = 0. Conditions (13) are obtained similarly by taking ∇ τl J = 0 and ∇ γm J = 0. Theorem 4.1 asserts that any solution of the H 2 model approximation Problem 2.1, denoted by Ĥd = ∆o Ĥ ∆i is s.t. Ĥ satisfies, at the same time, a set of 3n bi-tangential interpolation conditions detailed in [START_REF] Opmeer | Optimal model reduction for non-rational functions[END_REF] and another set of n u + n y relations on the delays contained in the ∆i and ∆o diagonal matrices [START_REF] Halevi | Reduced-order models with delay[END_REF].

Remark 4.1. (H 2 optimality conditions in the SISO case) In the SISO case, all the conditions provided in Theorem 4.1 appear much simpler and can be stated as follows. Considering:

G(s) = N j=1 ψ j s -µ j , Ĥd (s) = n k=1 φ k e -τs s -λk , s.t.
Ĥd is a local optimum of Problem 2.1, then the following conditions hold:

Ĥ(-λk ) = G(-λk ), Ĥ (-λk ) = G (-λk ), ( 14
) N j=1 µ j ψ j        n k=1 φ k µ j + λk        e τµ j = 0. ( 15
)
for all k = 1 . . . n, and where G is as in (11):

G(s) = N j=1 ψ j s -µ j e τµ j .
Remark 4.2 (Impulse response of G(s) and advance effect). The H 2 -optimality conditions given in Theorem 4.1 involves a model G(s) which has a pole-residue decomposition defined by [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF].

For simplicity, let us consider the SISO case where G and G is given by

G(s) = N j=1 ψ j s -µ j G(s) = N j=1 ψ j s -µ j e µ j τ .
Thus, the the impulse response of G(s) is

g(t) = N j=1
ψ j e µ j t e µ j τ 1(t) = N j=1

ψ j e µ j (t+τ) 1(t)

= g(t + τ)1(t), t ∈ R
where 1(t) corresponds to the Heaviside step function and g(t) is the impulse response of model G(s). Therefore, G(s) behaves as a time advance of G(s) and correspond to the "causal part" of the model G(s)e sτ .

Practical considerations

In this subsection, three considerations about Problem 2.1 and Theorem 4.1 are discussed. These latter are relevant to sketch an algorithm which enables the computation of model ∆o Ĥ ∆i satisfying the optimality conditions of Theorem 4.1. Let us consider that Ĥd = ∆o Ĥ ∆i is a local minimum of the H 2 optimization Problem 2.1 where Ĥ is given by ( 6), then:

• Consideration . If the matrices ∆o , ∆i and the reduced order model poles λ1 , λ2 , . . . , λn are assumed to be known, Problem 2.1 is reduced to a much simpler problem that can be solved, for example, by using the well-known Loewner framework such as in [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF];

• Consideration . If the delay matrices ∆o , ∆i are known, then Problem 2.1 can be solved by finding a model realization Ĥ which satisfies the interpolation conditions (12) of Theorem 4.1, only. This can be done using, for instance, a very efficient iterative algorithm, e.g., IRKA (see [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF]);

• Consideration . Assume that the system realization Ĥ has already been determined. It follows that Problem 2.1 is equivalent to look for optimal delays matrices (

∆ o , ∆ i ) ∈ H n y ×n y ∞ × H n u ×n u ∞ s.t.: ( ∆ o , ∆ i ) = argmax ( ∆o , ∆i ) ∆o Ĥ ∆i , G H 2 . (16)
Interestingly, since ∆o Ĥ ∆i , G H 2 → 0 when the delays go to infinity, this problem can be restricted to a compact set and thus a global solution exists.

Computational considerations

An algorithm which allows to numerically compute a model Ĥd satisfying the previous H 2 optimality conditions is proposed in this subsection. It relies on the considerations above discussed (Section 4.1). Therefore, the proposed approach corresponds to an iterative algorithm in which each iteration can be decomposed in two steps. The first one aims at computing a realization Ĥ which satisfies the interpolation conditions [START_REF] Opmeer | Optimal model reduction for non-rational functions[END_REF] while fixing the matrices ∆o , ∆i at their values obtained from the previous iteration. This can be done using, for instance, the IRKA algorithm (Step 4). In the second step, the resulting Ĥ is then exploited to determine the n u + n y optimal values for the ∆o , ∆i matrices elements (Step 5). This step is achieved by solving the nonlinear optimization problem defined in ( 16) using an appropriate solver. Then, the whole process is repeated and these two steps performed again until the convergence 2 . At the end of the procedure, the model built will satisfy the H 2 optimality conditions on which Theorem 4.1 relies. This sequential procedure can be summarized such as in Algorithm 1, and referred to as MIMO IO-dIRKA.

Algorithm 1 MIMO IO-dIRKA (MIMO Input Output delay IRKA)

Require: A N th -order model G ∈ H n y ×n u 2 , dimension n ∈ N * (n N) and initial guesses for both ∆it=0 i , ∆it=0 o . 1: while not converged do 2: Set it ← it + 1 3:
Build Git as in [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF] 4:

Build Ĥit satisfying the bi-tangential interpolation conditions (12) using IRKA [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF] on Git

Structured input/output delays

All the previous results are left unchanged in the case of structured input/output delays i.e., if, for example, delays does not apply on given input(s) and/or output(s) of Ĥd . The results can be derived in a straightforward way, without any loss of generality, just by considering the following ordered delays matrices (where delays are present on the first n d1 < n u inputs and n d2 < n y outputs): ∆i (s) = diag(e -sτ 1 , e -sτ 2 , . . . , e -sτ n d1 , 1, . . . , 1) ∆o (s) = diag(e -sγ 1 , e -sγ 2 , . . . , e -sγ n d2 , 1, . . . , 1).

One can easily note that the preliminary results from Sections 3 and 4 still remain true when introducing these matrices. The main result stated in Theorem 4.1 thus remains unchanged.

Numerical application

This section is dedicated to the application of the results obtained in Sections 4, namely, the input/output-delay optimal H 2 model approximation and its first -order optimality conditions. We will emphasize the potential benefit and effectiveness of the proposed approach.

Let us consider a model G of order N = 20, given by the following transfer function

G(s) = N j=1 µ j s -µ j , (17)
where µ j ∈ R -( j = 1, . . . , N) are linearly spaced between [-2 -1]. The impulse response of G is given by the solid dotted blue line in Figure 1. Interestingly, it behaves like a system with an input delay. In order to fit the framework proposed in this paper, input-delay H 2 optimal model Ĥd = ∆o Ĥ ∆i of order n = 2 (solid red) was obtained by applying Theorem 4.1 and IO-dIRKA, as described in Section 4. The obtained delay model is compared with delay-free approximations of order n = {2, 3, 4}, obtained with IRKA3 . All the results are reported on Figure 1. As clearly shown on Figure 1, the proposed methodology allows to obtain an input-delay H 2 approximation of model G that clearly provides a better matching than the delay-free cases, even for higher orders (here, IRKA with n = 4 still have a bad matching and exhibits difficulties in accurately catching the delay and main dynamics). Indeed, the delay-free cases exhibits an oscillatory behaviour during the first seconds while the input-delay model Ĥd takes benefit of the delay structure to focus on the main dynamical effect. Moreover, the approximation model of Ĥd satisfies the conditions given in Theorem 4.1.

Remark 5.1 (Numerical results (SISO case, n = 2)). For sake of completeness, the optimal numerical values obtained with MIMO IO-dIRKA are: λ1,2 = -2.0320 × 10 -1 ± i 2.0700 × 10 -1 , φ1,2 = 1.5713 × 10 -3 ± i 1.8691 × 10 -1 and the optimal delay τ = 8.7179. The interpolation conditions can then easily be checked:

• Condition (14) leads to Ĥ(-λ1,2 ) = G(-λ1,2 ) = 2.3567 × 10 -1i 2.3614 × 10 -1 and Ĥ (-λ1,2 ) = G (-λ1,2 ) = 5.6466 × 10 -1 ± i 1.1465.

• When evaluating

N j=1 µ j ψ j        n k=1 φ k µ j + λk        e τµ j ,
one obtains 9.7284 × 10 -5 , which is close to zero, as stated by condition [START_REF] Duff | Realization independent single time-delay dynamical model interpolation and H 2 -optimal approximation[END_REF]. With reference to Figure 2, similar results are obtained in the case of an input delay-dependent approximation of order n = 4 (using IO-dIRKA) and delay-free approximation of order n = {4, 5, 6} (using IRKA). Then, Figure 3 shows the impulse response mismatch error for these different configurations. For each reduced order models, the mean square absolute error ε of the impulse response are computed. The main observation that can be made is that the mismatch error obtained for Ĥd of order n = 4 is lower that the one obtained by a delay-free model Ĥ of order n = 6 (a better result is obtained for a delay-free model with an order n = 7). This motivates the use of the specific approximation model delay structure. Figure 3: Impulse response error between the original model H of order N = 20 and the input-delay H 2 -optimal model Ĥd of order n = 4 (solid red line) and the delay-free H 2 -optimal models Ĥ of order n = {4, 5, 6} (dashed dark green, light green and yellow lines).

Conclusion

The main contribution of this paper is the derivation of the first-order H 2 optimality conditions for Problem 2.1. It forms a direct extension of the bi-tangential interpolation conditions of the delay-free case derived in [START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF][START_REF] Van Dooren | H 2 -optimal model reduction of MIMO systems[END_REF]. Theorem 4.1 establishes that if Ĥd = ∆o Ĥ ∆i is a local optimum, then the parameters of this latter verify an extended set of matricial equalities. These ones are of two types: first, (i) a subset of interpolation conditions [START_REF] Opmeer | Optimal model reduction for non-rational functions[END_REF] satisfied by the rational part Ĥ of Ĥd , which generalizes the delay-free case; secondly, (ii) a subset of matricial relationships [START_REF] Halevi | Reduced-order models with delay[END_REF] focussing on the input/output delay blocks ∆o , ∆i . These conditions all are dependent on the reduced order model parametrization described by bk , ĉk , λk , τl and γm , and solving Problem 2.1 requires to tackle a non-convex optmization problem. Nevertheless, an algorithm referred to as IO-dIRKA, has been proposed to practically address this issue. This latter decorrelates the decision variables between them by solving, firstly for given ∆i , ∆o matrices, an optimal H 2 approximation problem, and then, in a second stage, a nonlinear maximization problem [START_REF] Schulze | Data-driven interpolation of dynamical systems with delay[END_REF] to determine the optimal values of the delays. Both optimizations rely on descent methods, taking benefits from the analytical expressions of the gradients of the H 2 mismatch gap ∇J. Numerical experiment have also been presented, illustrating the benefit of the proposed approximation delay structure with respect to standard delay-free approximation methods.

Theorem 4 .

 4 1 gathers all the first-order optimality conditions related to Problem 2.1 and stands as the main result of the paper.

Theorem 4 .

 4 1. (Delay model approximation first-order H 2 optimality conditions) Let us consider G ∈ H n y ×n u 2 whose transfer function is G(s) ∈ C n y ×n u . Let Ĥd = ∆o Ĥ ∆i be a local optimum of Problem 2.1. It is assumed that Ĥ ∈ H n y ×n u 2 corresponds to a model with semi-simple poles only and whose transfer function is denoted by Ĥ(s) = Ĉ(s Ê -Â) -1 B ∈ C n y ×n u . Let ∆i , ∆o be elements of H n u ×n u ∞ and H n y ×n y ∞ , respectively, s.t. Propositions 3.1 and 3.3 are verified. Then, the following equalities hold:

5 :

 5 Determine ( ∆ i , ∆ o ) which solve (16Ĥd satisfies the interpolation conditions of Theorem 4.1.

Figure 1 :

 1 Figure1: Impulse response of the original model H of order N = 20 (solid dotted blue line), the input-delay H 2 -optimal model Ĥd of order n = 2 (solid red line) and the delay-free H 2 -optimal models Ĥ of order n = {2, 3, 4} (dashed dark green, light green and yellow lines).

Figure 2 :

 2 Figure2: Impulse response of the original model H of order N = 20 (solid dotted blue line), the input-delay H 2 -optimal model Ĥd of order n = 4 (solid red line) and the delay-free H 2 -optimal models Ĥ of order n = {4, 5, 6} (dashed dark green, light green and yellow lines).

  Ĥd, n = 4, τ = 6.4103 with mean error ε =2.55906e-07 (IO-dIRKA) Ĥ, n = 4 with mean error ε =3.75348e-05 (IRKA) Ĥ, n = 5 with mean error ε =4.94274e-06 (IRKA) Ĥ, n = 6 with mean error ε =4.56304e-07 (IRKA)

  The main objective addressed in this paper is to solve the following H 2 approximation problem:

	Problem 2.1. (Delay model H 2 -optimal approximation) Given a stable N th order system G ∈ H n y ×n u 2 , find a reduced n th order (s.t. n N) multiple-input/output delays model Ĥ d = ( Ê, Â, B, Ĉ, ∆i , ∆o ) s.t.:
	Ĥ	d = argmin ny ×nu Ĥd ∈ H 2

Proof. Observing that the poles of the complex function Ĥd (-s)G(s) are µ 1 , µ 2 , . . . , µ N ∈ C - and -λ1 , -λ2 , . . . , -λn ∈ C + , let us consider the following semi-circular contour Γ C located in the left half plane s.t.:Γ C = Γ I ∪ Γ R

In practice, different stopping criteria might be considered, e.g. (i) the variation of the interpolation points materialized by λk (k = 1, . . . , n), as in[START_REF] Gugercin | H 2 model reduction for large-scale linear dynamical systems[END_REF], (ii) the interpolation conditions check (Theorem 4.1) or (iii) the mismatch H 2 error check (if the order N of the original system is reasonably low).

Using the implementation available in the MORE toolbox[START_REF] Poussot-Vassal | Introduction to MORE: a MOdel REduction Toolbox[END_REF], http://w3.onera.fr/more/.


T Virin

J Koko
email: koko@isima.fr

E Piron

P Martinet
email: martinet@lasmea.univ-bpclermont.fr

M Berducat
email: michel.berducat@cemagref.fr

Spatial Optimization of Fertilizer Application by Centrifugal Spreading

Keywords: Augmented Lagrangian, Centrifugal spreaders, Discretization, Optimization

We study an approach for the optimization of spreading performed by centrifugal spreaders in order to minimize environmental effects due to application errors. Faced with a large scale problem, we divide the domain of study into subdomains so that each tramline is individually dealt with. In order to take into account the mechanical limits of the device, some inequality constraints are introduced. After cost function discretization, we use an augmented lagrangian algorithm associated with a L-BFGS technique to solve the problem. Results are presented for parallel and non parallel paths.

I. INTRODUCTION

Fertilization is a regular agricultural practice which is carried on in order to apply plant nutrients to satisfy the demand of crops. In most cases, this operation is performed thanks to a centrifugal spreader because of its low cost and robustness. Unfortunately, by using unsuitable strategy, this practice can result in ill environmental effects. Indeed, an over application mostly involves an over nitrate enrichment of surface waters leading to an excessive multiplication of algae which can result in fishes dying. In the contrary, when under-dosages occur, yield losses can be very important. Therefore, mineral fertilization is today pointed out as one of the major causes of the pollution of groundwaters and watercourses. Thus, the new rules imposed by the different governments in North America and Northern Europe (Bruxelles, 2005), lead us to investigate a new strategy to improve the quality of the distribution achieved by centrifugal spreaders. This paper presents an approach for optimization of fertilizer application achieved during spreading. Instead of using the traditional method based on the best arrangement of the transverse distribution, we consider an optimization method relying on a spread pattern model developed by [START_REF]Etude du procédé d'épandage centrifuge d'engrais minéraux[END_REF], [START_REF] Olieslagers | Fertilizer distribution modelling for centrifugal spreader design[END_REF]. Faced with a large scale problem we decompose the problem by using "sliding windows" and taking advantage of the symmetry properties of centrifugal spreading. Simulation results are presented for parallel and non parallel tramlines.

II. PROBLEM STATEMENT

Spreading regularity is achieved when the actual distributed amount of fertilizers in the field is equal to the prescribed dose determined by agronomical and pedological reasoning with respect to the crops. Here, we consider centrifugal spreaders equipped with two spinning discs which eject nutrients along each tramline in the field. The tractor-spreader combination is commonly tooled up with a GPS antenna, a radar speed sensor, an embedded computer, and an actuator. So, according to the location and speed of the tractor measured by the two first tools previously mentioned, the computer determines the actuator adjustments to control the dosage rate. The amount of applied fertilizers currently called spread pattern, has an irregular distribution shape which is often highlighted by the transverse distribution curve calculated by summing the amounts along each travel direction illustrated by the red curve in Figure (1(a)). Because of the spatial distribution heterogeneousness, the tractor follows back and forth paths to obtain an uniform deposit from transverse distributions summation for each successive travel shown in Figure (1(b)). As we can notice, transverse distributions play an important role in the fertilization strategy. Indeed, the device settings relie on their best arrangements according to the different trajectories that is to say working widths. The working width is a concept widely used in the fertilization community today and corresponds to the distance between two consecutive overlapping lines. Furthermore, when overlapping is optimal, these lines stand for symmetry axes that make two successive transverse distribution coincide. This variable is also linked to the applied mass flow rate during spreading determined by the mathematical relation m = (Q * × W × S)/600 where m is the mass flow rate (kg/min), Q * the prescribed dose (kg/Ha), W the working width (m) and S the speed of tractor (km/h). With this relation, we reason as if nutrients were homogeneously distributed by the machine onto a rectangular area which length is equal to the wished working width. In this case, we don't take into account the actual occurring phenomenon. Indeed, the true fertilizers global deposit in the field is due to spread pattern overlappings, and then is the result of the heterogeneous spatial distribution accumulation at each position of the spreader. Therefore, this practice can give satisfying results with regularly spaced parallel travel direction but is quickly inefficient when geometrical singularities are met in the field (end of field, irregularly spaced parallel paths, start and end of spreading) and then can produce local application errors as illustrated in Figure ( 2). Then, we can have over-dosage that can reach +80% and sub-dosage equivalent to -85%. Since the spreading strategy is traditionally based on the best arrangement of transverse distributions, some works were led to find appropriate trajectories with respect to these ones (Dillon et coll., 2003). Unfortunately, this strategy is unsuitable when tramlines are already fixed by others agricultural practices like sowing. Therefore, it is important to know how best arrange the shape and the placement of the true distributions, that is to say spread patterns, during spreading with respect to the imposed geometrical constraints in the field. This adjustment should be continuously performed for each GPS position of the machine by changing its settings. Then, in this paper, we focus on a method which computes optimal parameters that permit to obtain best spread patterns arrangement in the presence of imposed paths. Moreover, this study deals with only optimization along the trajectories in the field and not on the boundaries of this one.

III. COST FUNCTION AND OPTIMIZATION METHOD

In order to formalize a suitable optimization criterion, it is necessary to define at first the parameters used for the spread pattern model. So, we consider the field as a polygonal domain

Ω in R 2 and s : (0, T ) -→ R 2 a path in Ω that is to say s(t) = (s 1 (t), s 2 (t)) ∈ R 2 . s(t)
is assumed to be rectilinear. The mass flow rates at time t for the spinning disc on the right and the left of the spreader are respectively the functions m : (0, T ) -→ R and d : (0, T ) -→ R. The spread pattern shown in Figure 1(a) is in most cases described by its medium radius r(x, s(t)) and medium angle θ(x, s(t)). r(x, s(t)), varying with the speed of disc, is assumed to be the distance between the point x and s(t), while θ(x, s(t)), modifiable the fertilizers dropping point on the disc, corresponds to the angle between ---→ s(t)x and s(t). To simplify notations, we will consider r(x, t) and θ(x, t) instead of r(x, s(t)) and θ(x, s(t)). Besides, we consider the following functions ρ and ξ, the medium radius concerning respectively the disc on the right and the left of the device, defined by ρ : (0, T ) -→ R + , and ξ : (0, T ) -→ R + . Moreover, let us assume the functions ϕ : (0, T ) -→ R, the medium angle concerning the right spinning disc, and ψ : (0, T ) -→ R, the medium angle concerning the left spinning disc. According to [START_REF]Etude du procédé d'épandage centrifuge d'engrais minéraux[END_REF], the spatial distribution pattern for the right and left spinning disc, at time t, are respectively given by the functions q r : Ω × (0, T ) -→ R and q l : Ω × (0, T ) -→ R defined by:

q r (x, m(t), ρ(t), ϕ(t)) = τ • exp -(r(x, t) -ρ(t)) 2 2σ 2 r • exp -(θ(x, t) -ϕ(t)) 2 2σ 2 θ (1) q l (x, d(t), ξ(t), ψ(t)) = κ • exp -(r(x, t) -ξ(t)) 2 2σ 2 r • exp -(θ(x, t) -ψ(t)) 2 2σ 2 θ (2) with τ = m(t) 2πσ r σ θ and κ = d(t) 2πσ r σ θ
In (1) and (2), σ r and σ θ are known parameters and stand for the standard deviations concerning respectively the medium radius and the medium angle. If we define

M (t) = (m(t), d(t)) ∈ R 2 , R(t) = (ρ(t), ξ(t)) ∈ R 2 and Φ(t) = (ϕ(t), ψ(t)) ∈ R 2 ,
the distribution pattern is then obtained by the function q tot : Ω × (0, T ) -→ R defined by:

q tot (x, M (t), R(t), Φ(t)) = q r (x, m(t), ρ(t), ϕ(t)) + q l (x, d(t), ξ(t), ψ(t)) (3)
From this model, the actual distributed dose Q during the interval of time (0, T ) for a single tramline is given by the function Q : Ω × (0, T ) -→ R defined by:

Q(x, M, R, Φ) = T 0 q tot (x, M (t), R(t), Φ(t))dt (4)
For our problem, the objective is to reduce harmful fertilization effects by minimizing application error that is to say the difference between actual and desired dose for all points in the field. Then, we have to compute the optimal functions M , R and Φ, which minimize the functional:

F (M, R, Φ) = Ω [Q(x, M, R, Φ) -Q * ] 2 dx (5)
where Q * is the prescribed dose. It is well known that (5) cannot be analytically calculated. Thus, we use an approximative integration method. So, let us divide interval (0, T ) into n elements with equal length δ = T n so that t j = jδ, j = 0, 1, ..., n, with t 0 = 0 and t n = T . Then, we define

M j = M (t j ), R j = R(t j ), Φ j = Φ(t j ) and M =    M 0 . . . M n    , R =    R 0 . . . R n    , Φ =    Φ 0 . . . Φ n    .
Thus, by using the trapezoidal rule, (4) can be approximated by:

Q(x, M , R, Φ) = δ χ + n-1 i=1 q tot (x, M i , R i , Φ i (6)
with χ = q tot (x, M 0 , R 0 , Φ 0 ) + q tot (x, M n , R n , Φ n ) 2 (7)
The functional to be minimized is then given as:

F (M , R, Φ) = Ω [Q(x, M , R, Φ) -Q * ] 2 dx (8)
In order to take into account the mechanical limits of the device and not to untimely solicit actuators, the functions M , R, and Φ and their time derivative are subject to bound constraints. Then, the set of constraints is defined as:

S = (M , R, Φ) ∈ R 6(n+1) | M min ≤ M ≤ M max , R min ≤ R ≤ R max , Φ min ≤ Φ ≤ Φ max , |M i+1 -M i | ≤ αδ, |R i+1 -R i | ≤ αδ, |Φ i+1 -Φ i | ≤ βδ}
where α, β and γ are known parameters. Therefore, we obtain the nonlinear programming problem (P) given by:

(P) min (M ,R,Φ)∈S F (M , R, Φ) (9)
S is a bounded closed set, then it is a compact set in R 6(n+1) . Thus, according to the Weierstrass theorem, the problem (P) has at least one local minimum. In most cases, there exist several tramlines and the actual distributed dose for all trajectories is then obtained by the summation of the applied dose for each k indexed path:

Q(x, M, R, Φ) = w k=1 t k f t k i q(x, M (t k ), R(t k ), Φ(t k ))dt ( 10
)
where w is the number of paths and the trajectories s k (t) are defined in the interval (t k i , t k f ). If we assume that M k j = M (t k j ), R k j = R(t k j ) and Φ k j = Φ(t k j ), we can use the same discretization scheme as before. For practical reasons, like the prescribed dose map in Figure (2), the domain Ω is 1 mgridded. In order to lose informations as little as possible, at least 2 samples per mesh must be computed. Thus, if we take the example of a farmland with only 3 tramlines 100 m long, 3600 variables should be considered. It is then clear that we cannot directly solve the problem for the entire field. Faced with this difficulty, we decompose the problem (P) into sub-problems in order to deal with each trajectory individually. Given that the arable land contains w tramlines, Ω is at first divided into w subdomains. We use the following notations:

K 1 = {k ∈ N| 1 ≤ k ≤ w} , L 2 = {l ∈ N| ∀z ≥ 2 ∈ N, 2 ≤ l ≤ z} , K 2 = {k ∈ N| 1 ≤ k ≤ w -1} , Ω = k∈K 1 Ω k , K 3 = {k ∈ N| 2 ≤ k ≤ w} , Ω k = l∈L 1 Ω k l , L 1 = {l ∈ N| ∀z ≥ 2 ∈ N, 1 ≤ l ≤ z} ,
and the following definitions:

-

Ω k ∈ R 2 is the k th subdomain of Ω, -Ω k l ∈ R 2
is the l th subdomain of Ω k . The subdomains Ω k are defined so that:

-

∂Ω k ∩ Ω k+1 = s k+1 (t), ∀ k ∈ K 2 and ∀t ∈ (t k+1 i , t k+1 f ), -∂Ω k ∩ Ω k-1 = s k-1 (t), ∀ k ∈ K 3 and ∀t ∈ (t k-1 i , t k-1 f ).
To make easier to understand the spatial decomposition, the figure (3) illustrates the example of three parallel tramlines in a domain Ω with a rectangular geometry.

k l , 1 ≤ l ≤ 3, 1 ≤ k ≤ 3.
From this decomposition, if the vectors M k l , R k l , Φ k l are assumed to be respectively the restriction of M , R and Φ in Ω k l , we can also define S k l as the restriction of S in the same subdomain. Thus, by considering the symmetry properties of spreading exposed in the section II, we can then solve the problem (P) by defining a new formulation (P ) of this one given by:

(P ) min z l=1 w k=1 J k l (x, M k l , R k l , Φ k l ) (M k l , R k l , Φ k l ) ∈ S k l , (l, k) ∈ L 1 × K 1 (11) with J k l (x, M k l , R k l , Φ k l ) = Ω k l [Q k l -Q * ] 2 dx ( 12
)
where Q k l stands for the applied dose in the subdomain Ω k l . Then, solving the problem (P ) is equivalent to optimize the distributed dose in each subdomain Ω k l , called "window", by taking into account not only the dose already applied in the domains Ω k l-1 and Ω k-1 l , but also the future distributed amounts in Ω k+1 l predicted so that they respect the previously exposed symmetry relations. Therefore, the notion of "sliding windows" can be explained by the fact of optimizing fertilization accuracy by solving in a sequential way the problem (P). Given that the problem (P ) is an optimization problem subject to inequality constraints, minimizing the functional J k l for (l, k) ∈ L 1 × K 1 is equivalent to consider the following problem:

P ineq      min J k l (M k l , R k l , Φ k l ) min j ≤ h j (M k l , R k l , Φ k l ) ≤ max j , j = 1, 2, ...6 × dim(M k l ) -1 (13)
where h j denotes the j th double inequality, min j and max j its lower and upper bound. To avoid to deal with the problem which consists in looking for saturated constraints, we use an augmented lagrangian algorithm which permits to penalize severely unacceptable solutions [START_REF] Bertsekas | Constrained Optimization and Lagrange Multipliers Methods[END_REF]. Thus, let us assume:

P r                    min X∈S k l J k l (X) + j c j (h j (X), λ j , r j ) X = M k l , R k l , Φ k l ¡ ∈ S k l λ j ∈ R r j ∈ R j = 1, 2, ...6 × dim(M k l ) -1 (14)
where λ j is the j th lagrange multiplier, r j the j th penalty coefficient and c j defined by:

c j (h j (X), λ j , r j ) =      λ j (h j (X) -max j ) + r j 2 (h j (X) -max j ) 2 if Z(h j (X), λ j , max j ) > 0 λ j (h j (X) -min j ) + r j 2 (h j (X) -min j ) 2 if Z(h j (X), λ j , min j ) < 0 - λ 2 j 2r j otherwise with Z(h j (X), λ j , b(h j )) = λ j + r j (h j (X) -b(h j
) and b(h j ) the bounds on h j . Faced with a cost function and gradient which evaluations present a high computational time, we choose to apply a L-BFGS algorithm shown to be efficient in this case (Byrd et coll., 1994). Then, we obtain the following algorithm:

Step 0: k = 1, l = 1, X in S 1 1 , Step 1: if l ≤ z,

IV. RESULTS

In this study, we only interest in rectangular domain. Nowadays, even in the cases where uniform rate application are considered, over and under-application are often observed within farmlands. So, here we define a constant prescribed dose fixed at 100 Kg/Ha. The speed of tractor commonly evolving around 10 Km/h, we use this average value for the machine displacement velocity. We consider two scenarios: the case of three parallel tramlines with a working width equal to 24 m and the case of two consecutive 24 m spaced parallel paths preceded by an other presenting a travel direction shift.

In the first case, the domain Ω is defined by Ω = (-24, 24) × (0, 44). We choose z = 4, w = 3 and thus divide the domain into 12 windows. According to the discretization technique, we obtain dim(M k l ) = 81 for l = 1, 2, 3, 4 and k = 1, 2, 3. The number of variables is then 486 per window. After optimization with Matlab environment, we obtain in figure (4(a)) an absolute error inferior to 0.1%, wich is very satisfying. Since we study the case of 3 tramlines and because of the symmetry properties previously explained, we choose to present only the parameters for the outward path which is between the two other trajectories. By observing the optimal parameters in figure (4(b)), we can notice some oscillations for the medium radius and the medium angle. This phenomenon can be explained if we make reference to the natural circular movement which is achieved when trying to uniformly paint a surface with a circular brush. Furthermore, the mass flow rate has a mean closed to 20 Kg/min wich is the value found by applying the relation m = (Q * ×W ×S)/600. So, the computed values for the other parameters are coherent with the actual ones determined by the manufacturers settings. For the second case, Ω = (-22, 24) × (0, 15) and 3 windows are considered (w = 3, z = 1). By using the same discretization techniques, v k l = 79 for l = 1 and k = 1, 2, 3. The total number of variables per window raises then 474. The non parallel tramline has an angle equivalent to 81 degrees which is current in fields which boundaries are not orthogonal between them. The error after optimization and the optimal parameters for the outward path are represented by Figure ( 5). In this case, the absolute value of the error is slightly superior to the one mentioned above but is all the same inferior to 0.5%. Because of the narrowing occuring on the left of the outward path, we can notice a little increase in the mass flow rate while the other variables seem not to be affected by this geometrical singularity. This result suggests then that the direction shift is not large enough to cause a medium radius and angle adjustment. The results shown in Figure 4(b) and 5(b) respecting the mechanical constraints of the machine, the set of computed parameters should be stored in the embedded computer and used as reference variables for a future control of the spreader. Therefore, the speeds of discs and granulars dropping points values could be continuously adjusted during spreading process as recommended in (Olieslagers et coll., 1997) and maintain a correct distribution. V. CONCLUSION Because of environmental and economic pressures, an optimization method has been developed for minimizing application errors due to centrifugal spreading. Considering the great size of problem, we have divided the global spatial domain into subdomains with rectangular geometry to deal with each trajectory individually and to reduce the number of decision variables. Faced with a problem subject to inequality constraints, an augmented lagrangian algorithm combined with a l-bfgs technique has been implemented. The solutions have been computed for parallel and non parallel tramlines and permit to obtain an error inferior to 1% in either case. From these results, we could deduce the values of the disc speeds and fertilizers dropping points during time and consequently carry out the control of the applicator. This work suggests a study is needed about spreading on the boundaries in order to improve application accuracy in the whole arable land. This future study should enable to limit distributed amounts outside the field and reduce waste of fertilizers.

Figure 1 :

 1 Figure 1: Centrifugal spreading characteristics: (a) Spread pattern and its tranverse distribution -(b) Fertilization strategy based on transverse distribution summation.

Figure 2 :

 2 Figure 2: Application errors resulting from the reasoning based on the best transverse distribution investigation and not on spread patterns overlappings.

Figure 3 :

 3 Figure 3: Rectangular domain Ω divided into 9 subdomains Ω k l , 1 ≤ l ≤ 3, 1 ≤ k ≤ 3.

  minimize J k l by using the augmented lagrangian algorithm associated with the L-BFGS minimization technique, otherwise goto Step 3, Step 2: l ← l + 1, goto Step 1 Step 3: k ← k + 1, X in S k 1 , goto Step 1.

  (a) Error in % (b) Optimal parameters for the outward path

Figure 4 :

 4 Figure 4: Error and optimal parameters for parallel tramlines, CPU=100h.

  (a) Error in % (b) Optimal parameters for the outward path

Figure 5 :

 5 Figure 5: Error and optimal parameters for non parallel tramlines, CPU=26h.