Daniele Astolfi

Laurent Praly

Lorenzo Marconi

Francis-Wonham nonlinear viewpoint in output regulation of minimum phase systems

 showing that asymptotic regulation can be achieved with an infinite dimensional regulator embedding a linear internal model copying all the harmonics that are multiple of the one associated to the exosystem. The regulator follows a post-processing structure in which the (infinite dimensional) internal model is driven by the error and a static stabiliser is considered for the cascade extended system. The post-processing structure traces the one proposed by Francis and Wonham into a linear framework. The presented analysis is limited to nonlinear systems that have unitary relative degree and that are minimum-phase.

INTRODUCTION

The problem of asymptotically offsetting the effect of an exosystem-generated signal on some regulated outputs of a system, typically known as output regulation problem, is fundamental in control theory. A landmark in the field was undoubetly given by the set of papers by [START_REF] Francis | The linear multivariable regulator problem[END_REF]Wonham (1975, 1976), [START_REF] Francis | The linear multivariable regulator problem[END_REF], [START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF], in which the problem was fully addressed for linear systems and exosystems. In those papers the necessity of internal model-based controllers to achieve asymptotic regulation was for the first time proved and fully characterized also in terms of robustness to possible parametric uncertainties affecting the controlled system, leading to the celebrated internal model principle. The results obtained for linear systems then pushed research activities in the nonlinear field with the first pioneering results that came out in [START_REF] Huang | On a nonlinear multivariable servomechanism problem[END_REF], [START_REF] Isidori | Output regulation of nonlinear systems[END_REF]. In those papers the center manifold theory was the key tool to prove the necessity of the internal model whenever the problem of output regulation must be solved in spite of (small) parameter uncertainties. Those papers opened a very rich research period in nonlinear output regulation in which many attempts have been done to make the nonlinear framework even more general and constructive. Milestone contributions within the last 20 years of research in nonlinear output regulation were certainly the paper by [START_REF] Serrani | Semiglobal nonlinear output regulation with adaptive internal model[END_REF], in which the case of uncertain exosystems was for the first time addressed opening the field of adaptive output regulation, the paper by [START_REF] Byrnes | Limit sets, zero dynamics, and internal models in the problem of nonlinear output regulation[END_REF], in which a "non-equilibrium theory" for output regulation was laid, the paper by [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF] in which a clear link between the design of nonlinear internal model-based regulator and nonlinear (high-gain) observers was established, and the work of [START_REF] Marconi | Output stabilization via nonlinear luenberger observers[END_REF] in which a complete theory of nonlinear output regulation without the so-called immersion assumptions was proposed. All the efforts were essentially restricted to consider single input-single error regulated systems with a well defined normal form between the control input and the regulation error and with minimum-phase features The noteworthy research attempts in the last 5 years then shifted the attention to deal with more general classes of multivariable systems without necessarily well defined relative degree and minimum-phase behaviours, with the final goal of extending in its full generality the linear theory. The search of design methods able to cope with more general nonlinear systems immediately revealed a substantial difference between the regulator structure originally proposed in the linear framework by Francis, Wonham and Davison and the ones that implicitly came out in the nonlinear framework. In the former, in fact, the regulator is given by a copy of the exogenous dynamics (acting as internal model) processing the regulation error and by a stabiliser designed to stabilise the extended system given by the regulated system driving the internal model. The structure that was used in all (to the best of authors' knowledge) the nonlinear literature, on the other hand, has a "swapped" structure with an internal model acting on the input of the regulated plant and a stabilisation unit stabilising the cascade of the internal model driving the regulated plant. This justified the terminology of post-processing and pre-processing internal model-based regulators, used to classify respectively the two previous structures, proposed in [START_REF] Isidori | Shifting the internal model from control input to controlled output in nonlinear output regulation[END_REF]. In that work, moreover, it was shown that a post-processing structure is a more natural and general choice to deal with multivariable (not necessarily square) systems. This fact, in turn, justified a number of attempts that have been recently done in nonlinear output regulation to adapt existing (pre-processing) internal model structures by "shifting" from the input to the error the internal model. This is, for instance, the case of [START_REF] Isidori | Shifting the internal model from control input to controlled output in nonlinear output regulation[END_REF] and [START_REF] Astolfi | Nonlinear output regulation by post-processing internal model for multi-input multi-output systems[END_REF] in which the pre-processing structures respectively presented in [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF] and [START_REF] Marconi | Output stabilization via nonlinear luenberger observers[END_REF] were shifted on the output by thus obtaining post-processing internal model-based regulators. Those papers, though, are still limited to the class of nonlinear systems having normal forms and minimum-phase. A more general class of regulated systems was addressed in [START_REF] Poulain | Robust asymptotic stabilization of nonlinear systems by state feedback[END_REF] and [START_REF] Astolfi | Integral action in output feedback for multi-input multi-output nonlinear systems[END_REF], under the assumption that the exosystem is a simple integrator (namely that the desired steady state is an equilibrium point). This framework was then extended in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF] by considering a linear periodic exosystem and showing how practical regulation can be achieved (under the assumption of "small" exogenous signals) by means of a post-processing regulator structure having an internal model copying the exosystem frequencies and a stabiliser stabilising the nonlinear system and the internal model via forwarding. The remarkable (as far as the present work is concerned) result proved in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF] is that the steady state error is a periodic signal (with the same period of the exosystem) and that its harmonic components at the frequencies that are copied in the internal model are in fact zero. Practical and not asymptotic regulation is, however, the only result that can be claimed for general nonlinear systems due to the fact that the steady state signals to cope with have in general an infinite number of harmonics. To the best knowledge of the authors there are not in literature systematic results for designing postprocessing asymptotic regulators for fairly general classes of regulated systems and exosystems. In this respect, in fact, the difficulties showing up are summarised in Bin et al. (2018) and pertain an intertwining between the design of the internal model and the design of the stabiliser that makes very difficult the systematic design and lead to a "chicken-egg" dilemma, by using the terminology of Bin and Marconi (2018), in the regulator construction. This paper is strongly motivated by the results obtained in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF] and, specifically, by the fact that, in that "periodic" framework, a post-processing regulator embedding a finite-dimensional internal model copying the harmonics of the exosystem generates a steady state error whose Fourier coefficients associated at that harmonics are precisely zero. We are thus naturally interested to answer the question: can we succeed in achieving asymptotic regulation in the nonlinear case if the internal model is infinite dimensional and containing all the harmonics multiple of the basic one associated to the exosystem period? The answer is affirmative and the goal of the paper is to present the main technicalities to prove it still in the case the system has a well defined relative degree and it is minimum-phase.

Notation R is the set of real numbers, R >0 := (0, ∞), R ≥0 := [0, ∞), N the set of non-negative integer numbers and N >0 the set of strictly positive integer numbers. We denote with | • | the standard Euclidean norm. For x ∈ R n , we denote with x its transpose. For x ∈ R n , y ∈ R m , we compactly denote (x, y) := (x , y ) . We define with P 1 T (X) the set of C 1 continuous functions R → X which are T -periodic. Recall, that by standard results on Fourier theory, if ε ∈ P 1 T (R), then ε is equal to its Fourier series, namely ε

(t) = 0 + ∞ k=1 c k cos(k 2π T t) + s k sin(k 2π T t) where 0 = 1 T T 0 ε(t)dt, c k = 2 T T 0 e(t) cos k 2π T t dt, s k = 2 T T 0 e(t) sin k 2π T t dt.
Given an infinite dimensional vector of the form z = (z 0 , z 1 , . . . , z k , . . .) with z 0 ∈ R and z k ∈ R 2 , we define the associated norm • 2 as z 2 2 := ∞ k=0 z k z k , and we define with 2 the set of vectors z satisfying z 2 < ∞. The set 2 so defined is a Banach space. Finally, throughout the paper we will use the following definitions

S = 0 ω -ω 0 , G = 0 1 , ω = 2π T . (1)
2. REVIEW OF LINEAR AND NONLINEAR OUTPUT REGULATION

Nonlinear Byrnes-Isidori Solution

Consider single-input single output strongly minimumphase systems with unitary relative degree of the form ż = f (w, z, e) ė = q(w, z, e) + u

where

(z, e) ∈ R n × R is the state, u ∈ R is the control input. The initial conditions of (2) range in a compact set Z × E ⊂ R n × R.
By following [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF], w ∈ R nw is an exogenous signal which is supposed to be generated by an exosystem of the form ẇ = s(w).

(3)

The exosystem (3) is assumed to be Poisson stable (see [START_REF] Byrnes | Limit sets, zero dynamics, and internal models in the problem of nonlinear output regulation[END_REF]), and in particular, its solutions live in some compact set W ⊂ R nw for all time. Next, we suppose the existence of a differentiable function π : R → R n solution of

L s π(w) = f (w, π(w), 0) ∀w ∈ W .
With the map so defined, it can be verified that the set A ⊂ W ×R n defined as A := {(w, z) ∈ W ×R n : z = π(w)} is invariant for the zero dynamics of system (2). System (11) is also supposed to minimum phase as stated by the next assumption. Assumption 1. The set A is asymptotically and locally exponentially stable for system (2) with a domain of attraction of the form

W × D where D is an open set of R n such that Z ⊂ D.
Previous assumption allows to properly define the friend ψ, namely the steady-state input which is needed to steer constantly to zero the output, as

ψ(t) := -q(w(t), π(t), 0), ∀ t ≥ 0. ( 4
)
As a matter of fact, when selecting u = ψ, it can be verified that z = π(w) and e = 0 is an equilibrium of (2) for any w ∈ W . As a consequence, in order to build a dynamical system which is able to generate the function ψ in (4) (the internal model unit), we review the solution proposed in [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF]. The starting point is to suppose the existence of a regression law characterizing the friend ψ, as stated by the next assumption.

Assumption 2. There exist a number d ∈ N >0 and a locally Lipschitz function ϕ : R d → R, so that

ψ (d) (t) = ϕ(ψ(t), ψ (1) (t), . . . , ψ (d-1) (t)), ∀ t ≥ 0. (5)
Then, the proposed dynamical regulator solving the problem of output regulation is given by η

= Φ(η) -σD κ Γe u = -σe + Cη (6) where η ∈ R d is the state with initial conditions ranging in a compact set Ξ ⊂ R d , Φ(η) =     η 2 . . . η d ϕ s (η)     , Γ =    γ 1 . . . γ d    , C = (1 0 • • • 0), (γ 1 , . . . , γ d ) are selected so that the polynomial λ d + γ 1 λ d-1 +. . .+γ d-1 λ+γ d is Hurwitz, D κ = diag(κ, . . . , κ d ),
ϕ s is a globally Lipschitz function that agrees with ϕ on the final attractor, and σ, κ ∈ R >0 are high-gain parameters that, if chosen large enough, ensure asymptotic output regulation as stated in the following theorem. Theorem 1. [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF]). There exists κ > 0, and for any κ ≥ κ there exists σ > 0, such that for any σ ≥ σ the closed-loop trajectories (2), ( 3), ( 6) starting in Z × E × W × Ξ are bounded for all t ≥ 0 and moreover lim t→∞ e(t) = 0.

The main drawbacks of such approach can be summarized as follows.

• In practical applications, the function ϕ in ( 5) can be very difficult to obtain in closed form. • The regression law (5) implies smoothness of the function q, which has to be at least C d with respect to his arguments. • As shown in Bin et al. (2018), this solution is not robust with respect to C1 perturbations of the nominal function q since the function ϕ may vary according to such perturbations. This changes are not easy to map in ϕ s and, in practice, only practical regulation is obtained, see Astolfi et al. (2017). • Since the lower bounds g and σ depend on the Lipschitz constants of q, ϕ, g and σ can be very large. However, high-gain controller may be undesired in practical applications.

Note that also the alternative solution proposed in Marconi et al. ( 2007) presents similar drawbacks. The main limitation in particular is the design of the regulator in practical applications as it is based on the exact solution of a PDE. An alternative version, in which such PDE does not need to be solved, has been proposed in [START_REF] Marconi | Uniform practical nonlinear output regulation[END_REF], although only practical regulation is obtained and similar drawbacks hold.

The Linear Lesson

As explained in the introduction, the solution proposed in [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF] for the class of nonlinear system (2)-( 3) can be considered as a pre-processing design which differs from classical post-processing solution proposed in linear output regulation. In order to highlight this fact, we recall here the main results of [START_REF] Francis | The linear multivariable regulator problem[END_REF]Wonham (1975, 1976); [START_REF] Francis | The linear multivariable regulator problem[END_REF]; [START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF] by specializing their results to the class of minimum-phase systems with unitary relative degree. Although such restriction is not necessary, it is instrumental in this context to highlight the main differences between linear and nonlinear output regulation frameworks. In particular, consider a linear system of the form ż = F z + N e + P w ė = Lz + He + Qw + u (7)

with z ∈ R n , e ∈ R, u ∈ R and for some matrices F, N, P, L, H, Q of suitable dimension. In such context, also the signal w is supposed to be linear, which is therefore generated by a neutrally stable system. In order to simplify the rest of this section, we restrict here to the case in which w is a sinusoidal signal determined by only one frequency, namely it can be though as generated by an autonomous systems of the form ẇ = Sw (8) where w ∈ R 2 , S is defined in (1) for some ω ∈ R >0 . The more general case in which various frequencies and constant signals contribute to w can be easily dealt by superimposition of the effects. By following [START_REF] Francis | The linear multivariable regulator problem[END_REF]Wonham (1975, 1976); [START_REF] Francis | The linear multivariable regulator problem[END_REF]; [START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF], when the matrix F is Hurwitz, the output regulation problem can be addressed for system ( 7)-( 8) with the following two step procedure: P1) Extend the system (7) with the following internal model in which η ∈ R 2 and S, G selected as in (1), η = Sη + Ge .

(9) P2) Select σ, K so that the matrix

A := F N 0 L H -σ K 0 G S is
In order to show that output regulation can be achieved, the following two arguments can be used.

A1) Since the matrix A, which characterizes the unforced (namely w = 0) closed-loop system (7), ( 9), (10), is Hurwitz, the forced closed-loop system (7)-( 10) has an asymptotically stable steady-state solution (z, ē, η) which is fully characterized the Sylvester equation ΠS = AΠ + B where B = (P, Q, GQ). In particular, (z, ē, η) = (Π z w, Π e w, Π η w), where Π = (Π z , Π e , Π η ). A2) Geometric arguments, briefly reviewed in [START_REF] Byrnes | Output regulation of uncertain nonlinear systems[END_REF], can be used to exploit the properties of the internal model ( 9) to show that since S, G is a controllable pair, then Π η S = SΠ η + GΠ e implies necessarily Π e = 0, which implies ē = 0, namely asymptotic output regulation is achieved.

Opposed to the nonlinear case reviewed in Section 2.1, such approach is post-processing 1 In particular, it is readily seen that the design of the regulator ( 9) is a copy of the exosystem (8) driven by the regulated output e. This design is selected independently of the friend (4), namely, no matter what the right steady-state input q is needed, the internal model is not changed. This is a substantial difference between linear and nonlinear output frameworks which allows to conclude also robustness of the regulator ( 9), ( 10), with respect to parameter uncertainties of the matrices (F, N, P, L, H, Q), see Chapter 1.3 in [START_REF] Byrnes | Output regulation of uncertain nonlinear systems[END_REF] and Bin et al. (2018). The main obstructions that prevent the extension of this linear paradigm to the nonlinear framework (2) can be therefore summarized by the following two questions:

Q1) Can we mimic the two-step procedure P1) -P2) in the nonlinear framework (2)-( 3), in which the design of the internal model is not affected by the selection of the stabilizer and does not depend on the friend ψ? Q2) Can we find alternative arguments to the geometric analysis in A2) that can be extended to the nonlinear framework ( 2)-( 3) to conclude that output regulation is achieved, once boundedness of trajectories is established?

3. MAIN RESULT

Main Idea

In order to try to reply to the second question Q2) of previous section, we recall the main result established in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF], which motivates this work. Lemma 1. [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF]). Consider system η = Sη + Ge where S, G are defined in (1), η ∈ R 2 , e ∈ R, and suppose it admits a solution η ∈ P 1 T (R 2 ), ē ∈ P 0 T (R), with T = 2π/ω. Then, the Fourier coefficients of ē associated to ω are zero, namely 2

T T 0 ē(t) cos (ωt) dt = 2 T T 0 ē(t) sin (ωt) dt = 0.
Lemma 1 establishes that if (η, e) are T -periodic bounded functions, then the Fourier coefficients of e associated to the frequency ω = 2π/T must necessarily be zero. This fact has a certain number of interesting consequences/properties that can be highlighted.

The first consequence is that if the spectral content of e coincides with one single frequency component, then e must be necessarily zero. In particular, consider again the linear framework of Section 2.2. In light of ( 8), and by definition of S in (1), the signal w is a T -periodic trajectory containing only frequency contents at ω = 2π/T . As a consequence, in light of stability of the matrix A and by using linearity, it can be proved that the forced closed-loop system ( 7)-( 10), admits a unique Tperiodic solution (z, ē, η) ∈ P 1 T (R n × R × R 2 ) containing frequency components only at ω. As a consequence, direct application of Lemma 1 allows to conclude that ē = 0.

An second property of Lemma 1 to be highlighted is that such results hold as long as ē, η are bounded T -periodic functions, no matter how they are generated. This means that such result can be used also in the framework of nonlinear systems. In other words, Lemma 1, provide an answer to the question Q2) of Section 2.2 for the specific case of periodic exogenous systems. In particular, consider the nonlinear system (2), extended with the internal model ( 9), and suppose that we are able to design a controller of the form u = φ(e, η) able to guarantee the existence of an attractive and stable T -periodic solution (z, ē, η) for the closed-loop system when w is present. As a consequence, for such periodic solution, we can apply Lemma 1, to establish that the Fourier coefficient of ē corresponding to ω is zero. As shown in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF], however, such result is not enough to conclude that ē = 0, since ē could contain higher order Fourier coefficients which are caused by higher order deformations induced by the nonlinear functions f, q.

Finally, a third property of Lemma 1 is that superposition of effects can be used in the in which many systems of the form ηk = kSη k + Ge with k ranging in {1, . . . , N }, are considered, for which a unique solution (η 1 , . . . , ηk , ē)

∈ P 1 T (R 2 × • • • × R 2 × R) exists.
In this case, it can be shown, see [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF], that all the Fourier coefficients of ε corresponding to ω, . . . , N ω, are zero.

Francis-Wonham Nonlinear Viewpoint

The main idea of this work, is therefore that of building an internal model unit composed by an infinite number of oscillator of the form ( 9), so that we can apply Lemma 1 to any of the Fourier coefficients characterizing the Fourier developments of a periodic steady-state solution of (2). A first observation however is that in doing so, the model of the exosystem (3) generating w is no more needed. As a matter of fact, in view of previous consideration, we only need that the steady-state solution of the zero-dynamics of (2) are periodic, and so is the friend (4). For this, in the following, all the instances of w are replaced by t to point out that we only need functions which are periodic. In particular, the following system is considered ż = f (t, z, e) ė = q(t, z, e) + u (11)

where z ∈ R n , e ∈ R, u ∈ R. As in Section 2.1, we suppose that the initial conditions of (11) range in some given compact sets Z × E ⊂ R n × R, and the following assumptions, which replace Assumptions 1 and 2, are stated. Assumption 3. The functions f is continuous with respect to time and locally Lipschitz with respect to (z, e). The function q is C 1 with respect to (t, z, e). The functions f, q are T -periodic in the argument t, namely f (t + T, z, e) = f (t, z, e) and q(t+T, z, e) = q(t, z, e) for any t ≥ 0, z ∈ R n and e ∈ R. Assumption 4. The zero-dynamics ż = f (t, z, 0) admits a unique solution z ∈ P 1 T (R n ) which is asymptotically (and locally exponentially) stable with domain of attraction an open set A ⊆ R n containing Z.

Under previous assumptions, by following the two step procedure for linear systems in Section 2.2, we propose here a two-step procedure for the nonlinear system (11): P1) Design the internal model as a bunch of an infinite number of oscillators, namely η0

= e η 0 ∈ R ηk = kSη k + Ge η k ∈ R 2 (12)
where k ∈ N >0 , η = (η 0 , η 1 , . . .) = (η k ) k∈N is the state with initial condition η(0) ranging in some given closed set Ξ ⊂ 2 , and S, G are defined in (1).

Note that η 0 is an integral action needed to remove constant bias, see [START_REF] Astolfi | Integral action in output feedback for multi-input multi-output nonlinear systems[END_REF]. P2) Select the following regulator

u = -σe + σ ∞ k=0 M k (η k -M k e) (13)
where σ > 0 is a parameter to be chosen, and

M 0 ∈ R, M k ∈ R 2 , with k ∈ N >0 , are defined so that -σM 0 = 1, -σM k = k S + G, k ∈ N >0 . ( 14
)
Similarly to the nonlinear framework in Section 2.1, by using the compact notation Φ = blckdiag(0, ωS, 2ωS, . . .),

Γ = (1, G, G, . . .), M = (M 0 , M 1 , M 2 , . . .), (15)
the internal model unit ( 12) and the control law ( 13) can be compactly rewritten as η = Φη + Γe, (16a)

u = -σ(1 + M M )e + σM η, (16b
) where in this case the regulator is linear although of infinite dimension. It is readily seen that in the control law (16b) we have two terms: the first, -σe, is the stabilizing law, while the second, M (η -M e) is a term obtained by using the forwarding approach 2 and applied to the specific cascade system (11), ( 12). This term is therefore needed to stabilize the internal model unit ( 12). Note that with respect to the standard forwarding procedure, reviewed for instance in [START_REF] Poulain | Robust asymptotic stabilization of nonlinear systems by state feedback[END_REF], the control law ( 13) is here obtained by ignoring the term q(t, z, e) in the equation (11). By solving ( 14) we can also explicitly compute the matrices M k as follows

M 0 = - 1 σ , M k = 1 (kω) 2 + σ 2 kω -σ . ( 17
)
The following lemma state that the term M ∈ 2 , namely the term -σM M e in the control law (16b) is not diverging. Lemma 2. Let M be given by (17). Then, for any σ ≥ 1,

M 2 ≤ 1 σ π 2 . Proof. The result follows by computing M 2 2 = ∞ k=0 M k M k = 1 σ 2 + ∞ k=1 1 (kω) 2 + σ 2 = 1 σ 2 ∞ k=0 1 ( kω σ ) 2 + 1 ≤ 1 σ 2 ∞ 0 1 s 2 + 1 ds = 1 σ 2 π 2 .
Finally, we can state the following theorem, which is the main result of this work, claiming that the regulator (16) solves the robust asymptotic output regulation problem for system (11). Theorem 2. For any pair of functions f, q satisfying Assumptions 3-4 hold. there exist a σ ≥ 1 such that, for any σ > σ , any trajectory of the the closed-loop system (11), ( 16), starting inside Z × E × Ξ, is such that (z(t), e(t)) is bounded for all t ≥ 0, η(t) ∈ 2 for all t ≥ 0, and moreover lim t→∞ e(t) = 0.

2 See [START_REF] Mazenc | Adding integrations, saturated controls, and stabilization for feedforward systems[END_REF] and [START_REF] Poulain | Robust asymptotic stabilization of nonlinear systems by state feedback[END_REF], [START_REF] Astolfi | Integral action in output feedback for multi-input multi-output nonlinear systems[END_REF], for a review of forwarding approach.

Proof of Theorem 2

With the function z defined by Assumption 4 , let consider the following error coordinate system ζ := z -z whose dynamics are given by ζ

= F (t, ζ), F (t, ζ) := f (t, z + z, 0) -f (t, z, 0) . (18)
We have the following result, which can be inherited by standard converse Lyapunov theorems, see, for instance, in [START_REF] Marconi | Output stabilization via nonlinear luenberger observers[END_REF]. Lemma 3. There exist an open set A ⊂ R n , a C 1 function V : R × R n and class K ∞ functions α, α, ᾱ such that the following holds.

(1

) {ζ ∈ R n : ζ = z -z, z ∈ Z} ⊂ A. (2) The function V satisfies α(|ζ|) ≤ V (t, ζ) ≤ ᾱ(|ζ|) (19a) lim ζ→∂A α(|ζ|) = +∞ (19b) ∂V (t, ζ) ∂ζ F (t, ζ) + ∂V (t, ζ) ∂t ≤ -α(|ζ|) (19c)
for all t ∈ R ≥0 and z ∈ A.

(3) The functions α, α, ᾱ are quadratic around the origin and there exists a class K ∞ function β, linear around the origin, such that

∂V (t, ζ) ∂ζ ≤ β(|ζ|)
for all t ∈ R ≥0 and z ∈ A.

Now, let us introduce the following lemma, instrumental to the proof of Theorem 2. Lemma 4. For any ψ ∈ P 1 T (R) there exists b > 0, and, for any σ ≥ 1, there exist η ∈ 2 , satisfying η 2 ≤ b, such that the solution to η = Φη, η(0) = η (20a) satisfies ψ(t) = σM η (20b) for any t ≥ 0, with Φ, M defined by (15).

Proof. Given any function ψ : [0, T ] → R, let us introduce the following compact notation

F 0 (ψ) := 1 T T 0 ψ(t)dt, F c k (ψ) := 2 T T 0 ψ(t) cos 2kπ T t dt, ∀ k ∈ N >0 , F s k (ψ) := 2 T T 0 ψ(t) sin 2kπ T t dt, ∀ k ∈ N >0 .
Since ψ ∈ P 1 T (R), we can express it by a Fourier series. In particular we have that

ψ(t) = ψ 0 + ∞ k=1 ψ c k c k + ψ s k s k , (21)
where we compactly denoted ψ 0 = F 0 (ψ), and, for any

k ∈ N > 0, ψ c k = F c k (ψ), ψ c k = F s k (ψ)
, c k := cos(kω t), s k := sin(kω t), ω = 2π/T . Equation (20b) can be written by explicating each component as

ψ(t) = σM 0 η 0 (t) + σ ∞ k=1 M k η k (t) .
Note that for k = 0, equation (20a) reads η0 = 0, η 0 (0) = η0 , which implies η 0 (t) = η0 for all t ≥ 0. By identifying the latter with expression (21) and the definition of M 0 given in (17), we select η0 = -ψ 0 . For k > 0 we denote ψ k = (ψ c k , ψ s k ). The k-th equation of (20a) can be written as ηk = kωSη k , η k (0) = ηk . Its solution is given by

η k (t) = exp(kωS t)η k = c k s k -s k c k ηk,c ηk,s ,
We have also

M k η k (t) = kω ηk,c -σ ηk,s (kω) 2 + σ 2 c k + σ ηk,c + kω ηk,s (kω) 2 + σ 2 c s .
By identifying all the terms of the series ( 20) and ( 21) we obtain

σ kω ηk,c -σ ηk,s (kω) 2 + σ 2 = ψ c k , σ σ ηk,c + kω ηk,s (kω) 2 + σ 2 = ψ s k , which can be written in the compact notation σ (kω) 2 + σ 2 kω -σ σ kω ηk = ψ k . As a consequence we can select ηk = 1 σ kω σ -σ kω ψ k ,
for any k ∈ N >0 . In order to prove that η ∈ 2 , we compute

η 2 2 ≤ ψ 2 0 + ∞ k=1 1 σ 2 kω σ -σ kω 2 |ψ k | 2 ≤ ψ 2 0 + 2 ∞ k=1 kω σ 2 + 1 |ψ k | 2 .
Since ψ ∈ P 1 T (R), dψ dt is square integrable and there exists real numbers d 1 , d 2 > 0, independent of σ, such that

∞ k=1 |ψ k | 2 ≤ d 1 , ∞ k=1 k 2 |ψ k | 2 ≤ d 2 ,
and therefore

η 2 2 = (η k ) 2 + ∞ k=1 (η k ) 2 ≤ 2 1 2 ψ 2 0 + d 1 + ω σ 2 d 2 .
The proof concludes by selecting

b = √ 2 1 2 ψ 2 0 + d 1 + ω 2 d 2 1 2
and by noting that this number does not depend on σ.

With the previous lemma in mind, let ψ be defined as ψ(t) := -q(t, z(t), 0) . ( 22) with z given in Assumption 4. Since z ∈ P 1 T (R n ), the function q is C 1 and T -periodic in t, we conclude that ψ ∈ P 1 T (R). As a consequence, let ζ ∈ P 1 T ( 2 ) be solution defined in (20) with initial conditions η, and consider the following changes of coordinates where we used the following notation

z → ζ := z -z, η → ξ := η -M e -η, (
F (t, ζ) := f (t, ζ + z, 0) -f (t, z, 0) , N (t, ζ, e) := f (t, ζ + z, e) -f (t, z + z, 0) , p(t, ζ, e) := q(t, ζ + z, e) + σM ζ = q(t, ζ + z, e) -q(t, z, 0) , y := M ξ (25) Note that N (t, ζ, 0) = 0, p(t, ζ, 0) = 0, ∀ t ∈ [0, T ], ζ ∈ A.
Note that in view of the change of coordinates ( 23), the initial conditions of ξ range in some closed set X ⊂ 2 which is independent of σ. This can be proved by noting recalling that z(0) ∈ Z ⊂ 2 , e(0) ∈ E with E ⊂ R some given compact set, M ∈ 2 , M 2 does not depend on σ, and η 2 ≤ b with b independent of σ in view of Lemma 4. Now consider the Lyapunov function U (t, ζ, e, ξ) := V (t, ζ) + e 2 + ξ 2 2 , with V (t, ζ) given by Lemma 3. In view of (19b), there exists ĉ > 0, independent of σ, such that -2σe 2 -2σy 2 + 2σey where we used the definition of y in (25) and the fact that ξ Φξ = 0. Note that the latter can be verified component wise by recalling that ξ k Sξ k = 0. Therefore, we obtain

Z × E × X ⊂ Ω ĉ, with Ω ĉ = {(t, ζ, e, ξ) ∈ R ≥0 × R n × R × 2 : U (t, ζ, e, ξ) ≤ ĉ}. (26) Moreover, since the projection of Ω ĉ on (ζ, e) is compact, in light of Assumption 3, there exist L q , L N > 0 such that |p(t, ζ, 0)| = |q(t, ζ + z, e) -q(t, ζ + z, 0) +q(t, ζ + z, 0) -q(t, z, 0)| ≤ L q |e| + L q |ζ| |N (t, ζ, e)| ≤ L N |e| for all (t, ζ, e) ∈ Ω ĉ.
U ≤ -α(|ζ|) + 1 ν β 2 (|ζ|) + 2 ν |ζ| 2 -(σ -3L q -L 2 q ν -L 2 N ν)e 2 -(σ -L q -L 2 q
ν)y 2 where ν > 0 is a parameter to be chosen large enough. Since α is quadratic around the origin and β is linear around the origin, there exists ν > 0 large enough such that the following inequality holds να

(|ζ|) -2β 2 (|ζ|) -4|ζ| 2 ≥ 0 for all ζ ∈ Ω ĉ. Now let σ * = 3L q + L 2 q ν + L 2 N ν. As a consequence, there exists a ρ > 0, independent of σ, such that U ≤ - 1 2 α(|ζ|) -ρe 2 -ρy 2 (27)
for any σ > σ * . Previous inequality implies

U (t) ≤ U (0) ∀ t ∈ [0, τ ] . ( 28
) This says that the solution (ζ, e, ξ) remains in Ω ĉ for all time. So, from its definition τ is infinite. As the projection of Ω ĉ on (ζ, e) is compact, solutions t → (ζ(t), e(t)) are bounded for all t ≥ 0. Moreover, in light of ( 27), (28), we have that also that t → y(t) is bounded. Boundedness of y implies, by using its definition in (25) and in light of the definitions of ξ in (23) and Ω ĉ in (26), that η(t) ∈ 2 for all t ≥ 0. Next, by using (28) We discuss now the result of Theorem 2. First of all, we remark that asymptotic regulation is obtained robustly to any (small) C 1 model perturbation3 of the functions f, q as long as the resulting closed-loop trajectories of ( 11), ( 16), are bounded and T -periodic. As a matter of fact, the design of the regulator ( 16) does not make any explicit use of the functions f, q, since we only need to select the parameter σ large enough to achieve asymptotic regulation. His magnitude, in particular, has to be selected larger than the Lipschitz constants of q on the compact set of interest and so that stabilization of the zero-dynamics is achieved4 .

A second remark is that obviously the controller ( 16) is not implementable. Nevertheless, from a practical point of view, a finite version of the internal model ( 12), namely with k ∈ {1, . . . , N } for some N < ∞, is still very interesting if the closed-loop trajectories are bounded and T -periodic, since the result of Lemma 1 holds for these solutions, as showed in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF]. Moreover, recall that higher frequencies are naturally attenuated by low-pass behavior of system (11).

Furthermore, it is worth also noticing that asymptotic regulation could be obtained by means of a finite dimensional regulator if we are able to identify the spectral content of friend ψ(t) in (4). In particular, by developing the Fourier series of ψ as ψ(t) = ψ 0 + ∞ k=1 ψ c k cos(kωt) + ψ s k sin(kωt), suppose that we know the existence of a finite set of integers I ⊂ N so that ψ c k = ψ s k = 0 for all k ∈ N ⊂ I. In such case, we may design the regulator (16) with a finite number of states η k k ∈ I. Asymptotic regulation can be then proved by combining the technicalities of Lemma 1, Theorem 2, and the results presented in [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF].

Finally, we remark the following main differences that with respect to the methodology proposed by [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF] and reviewed in Section 2.1: we do not need the exact knowledge of the friend ψ in (4) and the regression law ϕ in Assumption 2; the function q is allowed to be only C 1 with respect to his arguments, while Assumption 2 implies the signals to be C d ; we restrict to the framework of T -periodic exosignals.

CONCLUSIONS

We showed that robust asymptotic output regulation can be obtained by means of a infinite dimensional internalmodel based regulator which follows the linear paradigm proposed by Francis, Wonham and Davison in the 70's, when periodic exosignals are considered. This work can be viewed as a first step to be developed in the following directions: a full characterization of its finite-dimensional version in the case of minimum phase systems by combining the results developed here with [START_REF] Astolfi | Approximate regulation for nonlinear systems in presence of periodic disturbances[END_REF]; the use of the proposed approach for more general classes of nonlinear systems, with a particular focus on the nonminimum phase case; the use of adaptive or identification techniques to address the case in which the period T is unknown.

  Hurwitz and select the regulator u = -σe + Kη .

  23) by which we obtain the following closed-loop system ζ = F (t, ζ) + N (t, ζ, e) ė = p(t, ζ, e) -σ e + σy ξ = Φξ -M p(t, ζ, e) -σM y (24)

  On the left hand side each term in the integrand is a function with non negative values, and, since, according to (24), the functions t → ė(t) and t → ζ(t) are bounded, the functions t → e(t) and t → ζ(t) are Lipschitz (and therefore uniformly continuous). Hence by applying Barbalat's lemma we get lim

	lim τ →∞ α(|ζ(t→∞ τ 0 1 2 |e(t)| = 0 , lim t→∞ |ζ(t)| = 0 , which concludes the proof. 2 3.4 Discussion of the Main Result	along solutions, we obtain

t)|) + ρe 2 (t) dt ≤ U (0) .

We refer also toBin and Marconi (2018) for a more detailed discussion about differences between pre and post-processing solutions.

In the sense defined in[START_REF] Astolfi | Integral action in output feedback for multi-input multi-output nonlinear systems[END_REF],Bin et al. (2018).

In the same spirit of semi-global stabilization problems of minimum-phase systems, see Chapter 9.3 in[START_REF] Isidori | Nonlinear control systems[END_REF].


Giacomo Casadei

Daniele Astolfi

Angelo Alessandri

Luca Zaccarian

Synchronization of interconnected linear systems via dynamic saturation redesign

come

Synchronization of interconnected linear systems via dynamic saturation redesign

Giacomo Casadei * Daniele Astolfi * * Angelo Alessandri * * * Luca Zaccarian * * * * , †

INTRODUCTION

Networks of dynamical systems and synchronization are omnipresent in our daily lives. Power networks [START_REF] Dörfler | Synchronization in complex oscillator networks and smart grids[END_REF]), flow networks [START_REF] Bürger | Dynamic coupling design for nonlinear output agreement and time-varying flow control[END_REF]), robot or vehicle fleets [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: Algorithms and theory[END_REF]), sensor networks [START_REF] Sivrikaya | Time synchronization in sensor networks: a survey[END_REF]), and social networks [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: Convergence conjectures and theorems[END_REF]) are some of the most popular examples that can be framed in this same context.

From the seminal works on consensus [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF]) and synchronization of identical linear systems [START_REF] Scardovi | Synchronization in networks of identical linear systems[END_REF]), in the last decade researchers have considered many aspects of networks control such as heterogeneous networks [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF]), switching networks [START_REF] Lu | Exponential synchronization of linearly coupled neural networks with impulsive disturbances[END_REF]), and high-order systems [START_REF] Seo | Consensus of high-order linear systems using dynamic output feedback compensator: Low gain approach[END_REF]). In more recent years the activity has focused on the more challenging problem of nonlinear networks. Different approaches such as passivity [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF]), dissipativity [START_REF] Stan | Analysis of interconnected oscillators by dissipativity theory[END_REF]) and high gain [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]) have been considered and the literature on the subject is now vast. Indeed, most of the results of the linear framework have been extended succesfully to the nonlinear context [START_REF] Abdessameud | Synchronization of nonlinear systems with communication delays and intermittent information exchange[END_REF], [START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF]).

These different approaches have in common the possibility to rewrite the synchronization has a stabilization problem of a certain set. For instance, the synchronization error with respect to the origin in [START_REF] Scardovi | Synchronization in networks of identical linear systems[END_REF], or passivity with respect to a desired target set in [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF]. While the stability of the aforementioned set and the impact of noise [START_REF] Wells | Control of stochastic and induced switching in biophysical networks[END_REF]), perturbations [START_REF] He | Accurate clock synchronization in wireless sensor net-works with bounded noise[END_REF]), and hybrid phenomena

The research leading to these results is partially funded by ANR via project HANDY, number ANR-18-CE40-0010.

(such as switching topology [START_REF] Casadei | About disconnected topologies and synchronization of homogeneous nonlinear agents over switching networks[END_REF] and open networks [START_REF] Hendrickx | Open multi-agent systems: Gossiping with random arrivals and departures[END_REF]) have been considered, novel approaches to the control design that allow mitigating the effect of these perturbations over the networks have still to be explored. Some authors have adapted classical H ∞ techniques to the case of networks of identical linear systems [START_REF] Huang | Robust h∞ synchronization of chaotic lur'e systems[END_REF]Feng (2008), Dal Col et al. (2018)). In the context of heterogeneous networks, [START_REF] Khong | A unifying framework for robust synchronization of heterogeneous networks via integral quadratic constraints[END_REF] studied the problem of synchronization by means of integral quadratic constraint (IQC). However, these techniques are often hard to extend to the nonlinear framework.

Recent works concerning the design of observers have suggested that the use of nonlinear functions such as saturations and deadzones might help in reducing the effect of perturbations affecting the measurements. Both in the context of linear [START_REF] Alessandri | Stubborn state observers for linear time-invariant systems[END_REF], [START_REF] Cocetti | On dead-zone observers for linear plants[END_REF]) and nonlinear systems [START_REF] Astolfi | Stubborn iss redesign for nonlinear high-gain observers[END_REF]; [START_REF] Cocetti | High-gain dead-zone observers for linear and nonlinear plants[END_REF]), the improvements obtained are significant. Inspired by these results and their applicability to the nonlinear framework, in this paper we investigate the use of saturations in networks with the aim of reducing the effect of impulsive perturbations acting over the systems communication lines. The perturbations may represent external/exogenous attacks to the network, commutations in the topology or the emergence of new nodes. First, we prove that the nominal behavior of the networks, namely the synchronization, is not compromised. Then, with the aid of simulation results, we show that the saturations helps in reducing the propagation of impulsive perturbations over the network. The paper is organized as follows.

In Section 2, we briefly review some basic results in the literature of synchronization. In Section 3, we formulate the problem at stake. In Section 4, the main result of the paper is stated. Finally, in Section 5, simulation results show the effectiveness of the proposed approach.

Notation R is the set of real numbers. On R n , we define

1 n = [1, • • • , 1] . On R n , we define the row-vector basis b 1 = [1, 0, • • • , 0], b i = [0, • • • , 0, 1, 0, • • • , 0] with all zeros except 1 in the i-th position, b n = [0, • • • , 0, 1], so that I n = col(b 1 , . . . , b n ). Let (x, y) := [x ,
y ] for any column vectors x and y. Given a symmetric matrix P ∈ R n×n , let us denote by λ m (P ) and λ M (P ) the minimum and maximum eigenvalues of P , respectively. Given y ∈ R and σ ∈ R ≥0 , we define sat σ (y) := max{-σ, min{σ, y}} and dz σ (y) := y -sat σ (y).

Graph theory In a general framework, a communication graph is described by a triplet G = {V, E, A} in which V is a set of n nodes V = {v 1 , v 2 , . . . , v n }, E ⊂ V × V is a set of edges e jk that models the interconnection between nodes with the flow of information from node j to node k weighted by the (k, j)-th entry a kj ≥ 0 of the adjacency matrix A ∈ R N ×N . We denote by L ∈ R n×n the Laplacian matrix of the graph, defined as

kj = -a kj for k = j, kj = n i=1 a ki for k = j .
For a time-invariant graph, the following result holds, see [START_REF] Godsil | Algebraic graph theory[END_REF]. Lemma 1. A time-invariant graph, is connected if and only if L has only one trivial eigenvalue λ 1 (L) = 0 and all other eigenvalues λ 2 (L), . . . , λ n (L) have positive real parts.

PRELIMINARIES

Consider a network of n identical agents described by ẋj = Ax j + u j y j = Cx j + Dw j (1) for j = 1, . . . , n, where x j ∈ R d , is the state, y j ∈ is the output, w j ∈ R is an impulsive disturbance acting over the output through the gain D ∈ R, (A, C) is a detectable pair matrices, and u j ∈ R d is the diffusive coupling control input to be defined.

Agents (1) are connected according to a directed graph G = {V, E, A}, fulfilling the following assumption. Assumption 1. The graph G contains at least a spanning tree and its Laplacian L is diagonalizable. As a consequence (see [START_REF] Godsil | Algebraic graph theory[END_REF]), there exists µ > 0 such that, for all i = 2, . . . , n the following holds:

Reλ i (L) ≥ µ .
In the unperturbed case, i.e. when w j = 0 for all j = 1, . . . , n, it is well known (see, e.g., [START_REF] Scardovi | Synchronization in networks of identical linear systems[END_REF] and references therein) that by assigning the input of all the agents j = 1, . . . , n as

u j = -K n i=1 ji y i , (2)
where ki denotes the (k, i) entry of the Laplacian matrix L, synchronization, namely asymptotic convergence to zero of x i -x k for all i, k ∈ {1, . . . , n}, is ensured for suitable selections of K (a possible one is described in Theorem 1 below).

To characterize the collective behavior, and as a guideline for selecting K, it is customary to write in compact form the network interconnection of (1) via ( 2)

as ẋ = [(I n ⊗ A) -(L ⊗ KC)] x
(3) where x := (x 1 , . . . , x n ). Following the same approach by [START_REF] Fax | Information flow and cooperative control of vehicle formations[END_REF], [START_REF] Seo | Output feedback consensus for high-order linear systems having uniform ranks under switching topology[END_REF], [START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF], the compact form (3) can be conveniently manipulated by introducing the transformation T ∈ R n×n defined as

T := 1 0 1×(n-1) 1 n-1 I n-1 , T -1 = 1 0 1×(n-1) -1 n-1 I n-1 , (4) which (from Assumption 1) satisfies L = T -1 LT = 0 L 12 0 (n-1)×1 L 22 , (5)
where eig{L 22 } = {λ 2 , . . . , λ n }.

Consider then the change of variables x = (T -1 ⊗ I d )x, with T as in (4). From the structure of T , one immediately realizes that x = (x 1 , x 2 -x 1 , . . . , x n -x 1 ) and the system in the new coordinates reads ẋ

= (T -1 ⊗ I d ) [(I n ⊗ A) -(L ⊗ KC)] (T ⊗ I d )x = (I n ⊗ A) -( L ⊗ KC) x .
Moreover, introducing e k = x k -x 1 , for k = 2, 3, . . . , n, and e := (e 2 , e 3 , . . . , e n ) , yielding x = col(x 1 , e), it is readily seen that the system above exhibits a triangular structure of the form ẋ1

= Ax 1 -(L 12 ⊗ KC)e (6a) ė = (I n-1 ⊗ A) -(L 22 ⊗ KC) e. (6b)
With this triangular structure, it is well known that assessing global exponential stability of e (namely exponential convergence to zero of the error coordinates) guarantees asymptotic synchronization of the agents to the perturbed trajectory of x 1 (which may be bounded or unbounded).

Theorem 1. Under Assumption 1, the following dual algebraic Riccati equation admits a solution P = P T > 0 P A + AP -2µP C CP + aI = 0 for any a > 0 and µ given in Assumption 1. Moreover, selecting K = P C in (2) ensures that the e-dynamics in (6b) are globally exponentially stable.

Proof. The proof follows from the construction in (Isidori et al., 2014, App. A). In particular, consider the invertible complex matrix M ∈ C n×n such that M L 22 M -1 = diag(λ 2 , . . . , λ n ). Then, define the quadratic function V (e) := e (M * M ⊗ P -1 )e = e (H R ⊗ P -1 )e, (7) where M * is the conjugate transpose of M and H R = Re(M * M ) is Hermitian and positive definite by construction. From the property of the Kronecker product, we get α 2 e 2 ≤ V (e) ≤ ᾱ2 e 2 , (8) with the positive scalars α 2 := λ m (H R )λ M (P ) -1 > 0 and ᾱ2 := λ M (H R )λ m (P ) -1 > 0. Moreover, from the calculations in [START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF](Isidori et al., , page 2689)), there exists α > 0, possibly depending on P , such that, for all i = 2, . . . , n, P -1 (A -λ i KC) + (A -λ i KC) * P -1 < -αI .

Then, according to (Isidori et al., 2014, eqn. (36)), we get

V = 2e (H R ⊗ P -1 ) ((I n-1 ⊗ A) -(L 22 ⊗ KC))e ≤ -αe (H R ⊗ I d )e , (9)
which concludes the proof.

PROBLEM FORMULATION

In this paper we are interested in studying the case in which impulsive perturbations act on the communication among agents (1), namely w j = 0 for some j. As mentioned in the introduction, w j can arise from hybrid phenomena (such as a switching topology) or glitches in the communication among agents. By changing coordinates according to (4), the network of ( 1), ( 2

T -1 L = N 12 N 22 ,
and w := (w 1 , . . . , w n ). The forthcoming results follows from Theorem 1 for the perturbed system (10).

Corollary 1. Consider the perturbed e-dynamics in (10b).

If K is selected according to Theorem 1, there exist c, γ > 0 such that the following bound holds

∇V (e), ė ≤ -c e 2 + γ w 2 (11)
for all e ∈ R (n-1)×d ,

x 1 ∈ R d , w ∈ R n .
Proof. Following the steps of the proof of Theorem 1, we obtain, along (10b), V = 2e (H R ⊗ P -1 ) Then,(11) follows from the Young inequality, with simple computations.

((I n-1 ⊗ A) -(L 22 ⊗ KC))e + (N 22 ⊗ KD)w ≤ -αe (H R ⊗ I d )e + e H R P -1 N 22 K w.
Corollary 1 guarantees that a bounded exogenous perturbation w always produces a bounded synchronization error e, namely practical synchronization is guaranteed. If now we are interested in mitigating the effect of w over the network, one could think of designing K to both guarantee synchronization and attenutation of perturbations (for instance, exploiting H ∞ methods as in [START_REF] Khong | A unifying framework for robust synchronization of heterogeneous networks via integral quadratic constraints[END_REF] or Dal [START_REF] Dal Col | H ∞ control design for synchronization of identical linear multi-agent systems[END_REF]). However, it is worth pointing out that both α and γ in (11) are functions of the degree of freedom K and structural limitations arise when a control input is restricted within the class of linear regulator, see [START_REF] Seron | Fundamental limitations in filtering and control[END_REF].

In this paper we propose a redesign technique that maintains the design of K proposed in Theorem 1 and adds a dynamic saturation effect over the diffusive coupling (2). First we show that the proposed technique still guarantees synchronization of the unperturbed network. Then, simulation results show that the impact of the disturbance w over the network is significantly reduced.

MAIN RESULT

In dynamic saturation redesign we introduce a saturation in the communication among agents. In particular, we replace the linear selection in (2) by

u j = -K sat √ σj n i=1 ji y i (12)
where the saturation level σ j ∈ R obeys the following dynamic depending on design parameters θ, r ∈ R

σj = -θσ j + r n i=1 ji y i 2 . ( 13
)
Note that the σ j dynamics in (13) can further expressed as σj = -θσ j + r ((b j Ly)

2 = -θσ j + r ((b j L ⊗ C)x) 2 = -θσ j + (b j L ⊗ C)xr((b j L ⊗ C)x) = -θσ j + (b j L ⊗ C)xrx (L b j ⊗ C ))
. As a consequence, the network (1) interconnected via ( 12) and ( 13) can be written in the following compact form ẋ = (

I n ⊗ A)x -(I n ⊗ K) sat √ σ ((L ⊗ C)x) σ = Θσ + (L ⊗ C)xRx (L ⊗ C ) ( 14
)
with Θ := diag{-θ, . . . , -θ}, R := diag{r, . . . , r}, and sat √ σ (•) defined as

sat √ σ ((L ⊗ C)x) :=      sat √ σ1 ((b 1 L ⊗ C)x) sat √ σ2 ((b 2 L ⊗ C)x) . . . sat √ σn ((b n L ⊗ C)x)      .
Theorem 2. Consider the network of n agents ( 14), under Assumption 1 with K designed according to Theorem 1. Then, for any θ there exists r > 0 such that for any r > r the set X = {(x 1 , x 2 , . . . , x n ) :

x 1 = x 2 = • • • = x N } (15) is globally exponentially stable.
Proof. We start by suitably rewriting the x-dynamics in (14

) as ẋ = (I n ⊗ A)x -(I n ⊗ K)(L ⊗ C)x + (I n ⊗ K)(L ⊗ C)x -(I n ⊗ K)sat √ σ ((L ⊗ C)x) = [(I n ⊗ A) -(L ⊗ KC)] x + (I n ⊗ K)dz √ σ ((L ⊗ C)x) , where dz √ σ (•) is defined as dz √ σ ((L ⊗ C)x) =      dz √ σ1 ((b 1 L ⊗ C)x) dz √ σ2 ((b 2 L ⊗ C)x) . . . dz √ σn ((b n L ⊗ C)x)     
.

We now apply the change of coordinates

x = (T -1 ⊗ I d )x with T as in (4) to obtain ẋ = [(I n ⊗ A) -( L ⊗ KC)]x +(T -1 ⊗ K)dz √ σ ((LT ⊗ C)x) σ = Θσ + (LT ⊗ C)xRx (T L ⊗ C ) (16)
Again, just as in ( 6), we want to exploit the structure of x in order to express in a more suitable form the dynamics (16). To this end, recall that we have e = (e 2 , e 3 , . . . , e n ) ,

x = (x 1 , e). ( 17
)
By definition of L and the structure of T in (4), we have

LT = 0 | M =     0 M 1 0 M 2 . . . . . . 0 M n     , (18)
where M ∈ R n× (n-1) . By using the structure in ( 18), we also obtain (LT ⊗ C)x = (M ⊗ C)e. As a consequence, in order to simplify the forthcoming computations, let us introduce the following notation

ζ := (LT ⊗ C)x = (M ⊗ C)e, ζ i := b i (LT ⊗ C)x = (M i ⊗ C)e ∀ i = 1, . . . , n. (19)
First of all, note that we can take advantage of ( 18) and ( 19), to get

dz √ σ ((LT ⊗ C)x) = dz √ σ (ζ) =    dz √ σ1 (ζ 1 ) . . . dz √ σn (ζ n )    . ( 20
)
As a consequence, by using (6b), ( 17), ( 19) and ( 20), dynamics ( 16)

reads ẋ1 = Ax 1 -(L 12 ⊗ KC)e + K dz √ σ1 (ζ 1 ) (21a) ėi = (A -b i L 22 KC)e i + K[dz √ σi (ζ i ) -dz √ σ1 (ζ 1 )] = Ae i + K[sat √ σi (ζ i ) -sat √ σ1 (ζ 1 )], i = 2, . . . , n, (21b) σi
= -θσ i + rζ 2 i , i = 1, . . . , n, (21c) Moreover
, according to (4), (17) (19), we have the following identities b

1 (T -1 ⊗ K)dz √ σ (ζ) = K dz √ σ1 (ζ 1 ), b i (T -1 ⊗ K)dz √ σ (ζ) = K[dz √ σi (ζ i ) -dz √ σ1 (ζ 1 )
], for all i = 2, . . . , n. Then the e-dynamics (21b) can be written in compact form, by also using (20), as ė = ((

I n-1 ⊗ A) -(L 22 ⊗ KC))e + K e dz √ σ (ζ) = (I n-1 ⊗ A)e + K e sat √ σ (ζ) (22)
where

K e := col(b 2 , . . . , b n )(T -1 ⊗ K).
As a consequence, similarly to the proof of Theorem 1, we consider a Lyapunov function for the e, σ dynamics which is independent of x 1 . In particular, consider the Lyapunov function

W (e, σ) = V (e) + n i=1 r -2 σ i + max{ζ 2 i -σ i , 0} (23)
with V defined as in (7). From ( 7) and ( 8) we have that there exist w > w > 0 satisfying w( e 2 + σ ) ≤ W (e, σ) ≤ w( e 2 + σ ). Now let I ⊂ {1, . . . , n} be the subset of indexes i for which σ i ≥ ζ 2 i and J = {1, . . . , n} \ I be the set of indexes j for which σ j < ζ 2 j . Then we have

σ i ≥ ζ 2 i ⇒ dz √ σi (ζ i ) = 0, ∀ i ∈ I σ j < ζ 2 j ⇒ dz √ σj (ζ j ) = 0, ∀ j ∈ J . (24)
Now, we compute the time derivative of W defined in (23). We obtain, using the definitions in ( 19),

Ẇ = V + i∈I r -2 σi + j∈J (r -2 -1) σj + j∈J 2ζ j (M j ⊗ C) ė
and therefore, using (21c), ( 22) and also (9),

Ẇ ≤ -αe (H R ⊗ I d )e + 2e (H R ⊗ P -1 )K e dz √ σ (ζ) + r -2 i∈I -θσ i + rζ 2 i + (r -2 -1) j∈J -θσ j + rζ 2 j + j∈J 2ζ j (M j ⊗ CA)e (25)
+ j∈J 2ζ j (M j ⊗ C)K e sat √ σ (ζ).
Let us now introduce the following positive scalars α and κ, independent of r, to simplify notation:

α(H R ⊗ I d ) ≥ α I > 0, κ := max j∈{1,...,n} { (H R ⊗ P -1 )K e 2 , (M j ⊗ CA) 2 , (M j ⊗ C)K e 2 (M ⊗ C) 2 },
and exploit the following bounds, resulting from the properties in (24) and definition ( 19):

dz √ σ (ζ) 2 ≤ j∈J ζ 2 j , sat √ σ (ζ) ≤ ζ ≤ (M ⊗ C) e . ( 26
)
Then we may refine the upper bound (25) for Ẇ (for reading convenience we preserve the position of each term in ( 25)), where we impose r > 1 (so that 1 -r -2 > 0) and apply standard Young inequalities1 in lines 1, 3 and 4, and use the properties of sets I, J in (24

): Ẇ ≤ -α e 2 + ν e 2 + κ ν j∈J ζ 2 j -r -2 θ i∈I σ i + n r ζ 2 -1- 1 r 2 j∈J r 2 -θ σ j + r 2 ζ 2 j + j∈J ν e 2 + κ ν ζ 2 j (27) + j∈J ν e 2 + κ ν ζ 2 j .
Finally, combining the different terms in (27) and using again the second bound in (26), we obtain

Ẇ ≤ -e 2 α -ν - n r (M ⊗ C) 2 - j∈J 2ν - r -r -1 2 - κ ν -2 κ ν j∈J ζ 2 j -r -2 θ i∈I σ i -1 - 1 r 2 r 2 -θ j∈J σ j . ( 28
)
Since κ is indepedent of r and ν, for any θ > 0 we can first select r ≥ r * ≥ ν -1 and ν small enough so that

α -ν -nν (M ⊗ C) 2 - j∈J 2ν > 0,
and finally fix r ≥ r * > 0 large enough so that r

-r -1 2 - κ ν -2 κ ν > 0, 1 - 1 r 2 r 2 -θ > 0.
By recalling that I ∪J = {1, . . . , n}, from (28) there exists a small enough ε > 0 satisfying Ẇ ≤ -ε( e 2 + σ ).

The last inequality, combined with the definition of W and (e, σ), completes the proof. The effect of a perturbation over the output of an agent has a variable impact depending on what output is affected. With the dynamic saturation design, this impact is highly reduced independently of the affected node.

In the subfigure 1b, norm u j of of each input without redesign (upper plot) and with saturation redesign (lower plot). Every time an impulsive perturbation occurs, the classic design provides and input about twice as large than the saturation redesign.

SIMULATION RESULTS

We consider a network of n = 6 linear oscillators, whose dynamics are ẋj1 = x j2 + u j1 , ẋj2 = -x j1 + u j2 y j = x j1 + w j (29) for j = 1, . . . , 6, exchanging their output information y j with their neighbors with control inputs u j = (u j1 , u j2 ) defined according to (2) and K chosen according to Theorem 1. An impulsive perturbation δ j acting over the communication is considered. In particular, we consider the case where the perturbation δ j acts on only one agent at a time. The same network is then considered with input designed according to (12)-( 13), with θ = 1 and r = 10. The network considered in this example is described by the Laplacian matrix

L =        -2 1 1 0 0 0 1 -4 1 0 1 1 0 1 -2 0 1 0 0 0 1 -2 1 0 0 1 0 0 -1 0 1 1 0 0 0 -2        . ( 30
)
Simulation results are shown in Figure 1, both for the linear design (upper plot) and the redesign proposed in Section 4 (lower plot): after reaching synchronization, the output of one of the oscillators is perturbed at t = 15 s. Then, each oscillator's output is perturbed with an interval of 5 s. First, in Figure 1a we can observe that the redesign improves the performances of the network importantly. While the nominal linear design fails to return rapidly to synchronization, the redesign helps to make the network more robust. Second, in Figure 1b we can clearly observe that the norm of the inputs u j without redesign and with saturation redesign is significantly different: thanks to the dynamic saturations not only the impact of the impulsive perturbations is highly reduced, but also the control effort is reduced to less than one half. This fact is a direct consequence of the evolution of the saturation levels σ j in Figure 2.

CONCLUSIONS

In this paper, a preliminary analysis of the use of dynamic saturations to reduce the effect of impulsive perturbations in networks has been developed. In the context of linear systems, we have shown that it is possible to redesign the classical linear diffusive coupling by adding a dynamic saturation, without compromising the synchronization property. Then, with the aid of simulations, we have verified that the dynamic saturations reduce the impact of impulsive perturbations over the communication among agents. In the near future, we would like to consider the possibility of decentralizing design of the dynamic saturations, namely each agent may have its own design parameters. More importantly, we aim to quantifying analytically the gain in terms of performances that we obtain with respect to the standard linear design. Another remarkable aspect of this approach is its applicability to the nonlinear framework. With respect to this, the research activity has already started and the first promising results are under investigation.

  ), reads ẋ1 = Ax 1 -(L 12 ⊗ KC)e + (N 12 ⊗ KD)w (10a) ė = (I n-1 ⊗ A) -(L 22 ⊗ KC) e + (N 22 ⊗ KD)w(10b) where N 12 , N 22 are matrices of appropriate dimensions satisfying

Fig. 1 .

 1 Fig. 1. In subfigure 1a, synchronization without coupling redesign (upper plot) and with coupling redesign (lower plot).The effect of a perturbation over the output of an agent has a variable impact depending on what output is affected. With the dynamic saturation design, this impact is highly reduced independently of the affected node.In the subfigure 1b, norm u j of of each input without redesign (upper plot) and with saturation redesign (lower plot). Every time an impulsive perturbation occurs, the classic design provides and input about twice as large than the saturation redesign.

Fig. 2 .

 2 Fig. 2. The evolution of the local saturation values σ j .

The Young inequality is the well-knonw upper bound

2ab ≤ νa 2 + b 2 /ν, holding for any a, b ∈ R and any ν > 0.


Aurelien Grolet
email: aurelien.grolet@ec-lyon.fr

Fabrice Thouverez
email: fabrice.thouverez@ec-lyon.fr

ON THE USE OF THE PROPER GENERALISED DECOMPOSITION FOR SOLVING NONLINEAR VIBRATION PROBLEMS

This paper presents the use of the so called Proper Generalized Decomposition method (PGD) for solving nonlinear vibration problems. PGD is often presented as an a priori reduction technique meaning that the reduction basis for expressing the solution is computed during the computation of the solution itself. In this paper, the PGD is applied in addition with the Harmonic Balance Method (HBM) in order to find periodic solutions of nonlinear dynamic systems. Several algorithms are presented in order to compute nonlinear normal modes and forced solutions. Application is carried out on systems containing geometrical nonlinearity and/or friction damping. We show that the PGD is able to compute a good approximation of the solutions event with a projection basis of small size. Results are compared with a Proper Orthogonal Decomposition (POD) method showing that the PGD can sometimes provide an optimal reduction basis relative to the number of basis components.

INTRODUCTION

The objective of this paper is to evaluate the efficiency of the proper generalized decomposition (PGD) for solving nonlinear vibration problems in the context of reduced order modeling. Most of the reduction method (e.g. krylov subspace [START_REF] Bai | Soar: A second order arnoldi method for the solution of the quadratic eigenvalue problem[END_REF][START_REF] Bai | Strcuture preserving model re-duction using a krylov subspace projection formulation[END_REF], proper orthogonal decomposition [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterisation and order reduction of mechanical systems: an overview[END_REF], dual modes [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF]) aim at reducing the model by first computing a set of vector which will constitute a reduced basis for solving the problem. This basis is computed relatively to some information already known about the system behaviour. For instance, the Proper Orthogonal Decomposition (POD) uses preliminary simulations (or mesurements) to extract a basis by using data analysis tools such as singular value decomposition. Thus those methods are often called a posteriori reduction methods. Contrary to those methods, the PGD aims at finding a reduced basis during the computation of the solution itself without any previous information on the system, thus making it an a priori reduction method. PGD was introduced under a different name (radial time-space approximation) by P.Ladeveze in the LATIN method [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF]. In the last decade, A.Nouy wrote a series of paper on the use of PGD for solving linear stochastic problems [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition for stochastic nonlinear problems[END_REF] and first order (in time) linear partial differential equations [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. The core of the method resides in separation of variable (also called radial decomposition by P.Ladeveze [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF]), alternate Galerkin projections in time and space and fixed point computation. In this paper the PGD is used for finding steady state solutions of nonlinear second order differential equation arising for instance when modeling nonlinear structural dynamic problems with finite element methods. We will show how the harmonic balance method (HBM) can be used in conjunction with PGD for finding periodic solution of the problem. A comparison will be drawn between PGD and POD on a numerical example. The paper is organized as follows: the first section presents generically the kind of nonlinear problems studied in this paper along with the description of Galerkin projections and the particular case of the Harmonic Balance Method. In the second section the POD and the PGD are presented in detail along with the modifications used to find periodic solutions. The last section is dedicated to numerical examples and to the comparison between POD and PGD.

NONLINEAR VIBRATION PROBLEMS AND GALERKIN PROJECTIONS

In this section we presents the kind of nonlinear problems studied in the paper and the application of Galerkin projection on such problems.

Nonlinear vibrations

Problems considered in this paper are related to the study of nonlinear vibration of structural systems. After discretization of the problem (e.g. by finite elements methods) one have to solve a systems of n second order nonlinear differential equations, where n is the number of degree of freedom (dof) of the model. Generally speaking the system can be put in the following form:

M ü u u + C u u u + Ku u u + g g g(u u u, u u u) = f f f (t) (1)
where M, C, K are respectively the mass, damping and stiffness matrices of size n × n, u u u is the vector of dof, f f f (t) is the vector of excitation forces which will be assumed to be harmonic with period T , and g g g is the vector of nonlinear forces acting on the system. In this paper we will consider a polynomial forces of maximum degree 3 corresponding to geometric nonlinearity, this force can be expressed as follows [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF]:

g g g i (u u u) = ∑ j k (1) i j u j + ∑ j,k k (2) i jk u j u k + ∑ j,k,l k (3) i jkl u j u k u l , 1 ≤ i ≤ n (2)
where k (1) , k (2) , k (3) are tensors of coefficient computed by the finite element model.

Galerkin projections

Here we briefly recall the results given by Galerkin projection applied to the system in Eq.( 1) as they will be used through the paper. Galerkin projection assumes a separation of variables such that a solution u u u can be expressed under the form u u u = Pq q q(t) = ∑ r i=1 p p p i q i (t) with P a matrix of size n × r and q q q(t) a vector of size r. A time (resp. space) projection basis of size r is then set and projections are computed in order to derive a purely spatial (resp. temporal) system.

Time projections

In the case of a projection in time, the basis q q q is set and we search for P such that the following holds for 1 ≤ i ≤ r: r ∑ j=1 [ I q i q j M + I q i q j C + I q i q j K]p p p j + I q i g g g(Pq q q, P q q q) = I q i f f f

(3) which can be rewritten as:

(I 2 ⊗ M + I 1 ⊗ C + I 0 ⊗ K) p p p + g g g( p p p) = f f f (4)
where p p p = [p p p T 1 , . . . , p p p T r ] T , g g g( p p p) = [ I q i g g g] 1≤i≤r and f f f = [ I q i f f f ] 1≤i≤r are vector of size nr and I k = [ I q i d k q j dt k ] are matrices of size r × r. The approximated solution u u u a (t) = Pq q q(t) is then computed by solving the nr nonlinear algebraic equations Eq.( 4). A typical use of Galerkin projection in time is the harmonic balance method described later in the paper.

Space Projections

In the case of a projection in space, the basis P is set and we search for q q q(t) such that the following holds: P T (MP q q q + CP q q q + KPq q q + g g g(Pq q q, P q q q)) = P T f f f [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF] or M r q q q + C r q q q + K r q q q + g g g r (q q q, q q q) = f f f r [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF] where A r = P T AP, A = M, C, K are reduced matrices of sizes r × r and g g g r , f f f r reduced vector of size r. The approximated solution is then computed by solving the r non linear differential equations Eq.( 6).

Harmonic Balance Method

The harmonic balance method is a particular case of Galerkin projection in time, and a widespread way to find periodic solutions of nonlinear systems such as Eq.(1) [START_REF] Sanliturk | Modelling two-dimensional friction contact and its application using harmonic balance method[END_REF][START_REF] Groll | The harmonic balance method with arc-length continuation in rotor stator contact problems[END_REF][START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF].

The basis used for this projection is a truncated Fourier basis up to H harmonics T T T (t) = [1, cos(ωt), sin(ωt), . . . , cos(Hωt), sin(Hωt)] and the time interval for integration is reduced to I = [0, 2π/ω]. In this case the integrals I k Eq.( 4) are simplified into I 0 = I d (identity matrix),

I 1 = D = ω diag(0, 0 -1 1 0 , ..., H 0 -1 1 0 ), I 2 = D 2 leading
to the following algebraic system with n = 2H + 1 equations:

L(ω) p p p + g g g( p p p) = f f f (7) with L = D 2 ⊗ M + D ⊗ C + I d ⊗ K a matrix of size (2H + 1)n × (2H + 1
)n, g g g and f f f vectors of size (2H + 1)n corresponding respectively to the nonlinear forces and the excitation forces in the frequency domain. When no analytical expressions can be derived for the nonlinear forces in the frequency domain g g g, an alternate frequency time (AFT) procedure is used to approximate g g g [START_REF] Cameron | An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems[END_REF].

Let distinguish the two cases of forced and free solutions. In the case of forced solution, the frequency ω is taken equal to the frequency of the harmonic force f f f , that is ω = 2π/T . In the case of free solution, the excitation forces are null, and the frequency ω is considered as an unknown. In this case a phase equation is added (for instance null velocities at t = 0) along with an equation controlling the amplitude of the displacements in order to avoid obtaining the null solution.

MODEL REDUCTION METHODS: POD AND PGD Proper Orthogonal Decomposition

The proper orthogonal decomposition is a model reduction technique which has been employed in various domain such as fluids [START_REF] Sterian | Low pressure centrifugal compressor[END_REF] or structural [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterisation and order reduction of mechanical systems: an overview[END_REF] dynamics. It consist in extracting a reduction basis from data acquired during numerical simulations or measurements. From a solution u u u = u u u(x x x,t) the reduction basis V(x) is computed such that it maximize the average in time of the inner product u u u(x,t), V . The result is a set of proper orthogonal modes (POM) V V V i (x x x) associated to a set of proper orthogonal values (POV) λ i . The participation of POM number i in the solution is given by the ratio λ i / ∑ j λ j . A convenient way to compute the reduced basis in the discrete case is the so called snapshot methods [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterisation and order reduction of mechanical systems: an overview[END_REF]. From a set of data U = u u u(x x x i ,t j ), the POD is obtained by performing a singular value decomposition of the matrix U, i.e. U = VSW T where V is the matrix of proper orthogonal modes and S the matrix containing the singular values on its diagonal. When carrying the POD, one observe that the magnitude of the POV is rapidly decreasing by several order of magnitude. This fact allow for model reduction by building a reduction basis P composed of only the r most significant POMs. The reduction is then carried out by searching an approximated solution u u u a of Eq.( 1) under the form u u u a (t) = Pq q q(t), and performing a Galerkin projection in space as described in Eq.( 6). The resolution of the reduced r differential equations finally allow to find the approximated solution u u u a .

Although the POD can give good approximations, its majors drawback resides in the fact that data is needed to build the reduced model. If the data comes from numerical simulations, this mean that one has to solve the n dof system in Eq.( 1) (e.g. by numerical integration, or HBM) for various loading case f f f . Those simulations can time consuming and are only valid for an excitation force f f f lying in the neighborhoods of the load cases used to derived the reduced basis.

Proper generalized decomposition

The PGD is a numerical method allowing to reduce models during the computation of the solution itself. The method and the related algorithms have been described in numerous papers [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition for stochastic nonlinear problems[END_REF][START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF] and especially in [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. It is based on variable separation and on an iterative procedure including alternate Galerkin projection in space and in time. In this section we describe how we applied the PGD to find approximated solutions of system in Eq. [START_REF] Bai | Soar: A second order arnoldi method for the solution of the quadratic eigenvalue problem[END_REF]. Two variants are presented, namely optimal Galerkin PGD (oPGD) and progressive PGD (pPGD), along with their corresponding algorithms.

Optimal

Galerkin PGD During the optimal Galerkin PGD, solutions are sought under the form u u u m (t) = ∑ m k=1 p p p k q k (t) = Pq q q(t). The matrix P and the time evolution vector q q q(t) are computed such that the following criterion hold simultaneously: B(u u u m , P * q q q(t)) = L(P * q q q(t)), ∀P * ∈ E oPGD S (8a)

B(u u u m , Pq q q * (t)) = L(Pq q q * (t)), ∀q q q * (t) ∈ E oPGD T (8b)
where E oPGD S is the set of real-valued matrix of size n × m, E oPGD T is the set of periodic functions vector of size m with period T = 2π/ω continuously derivable and B, L are two applications defined by the following relations:

B(u u u, v v v) = I v v v T (M ü u u + C u u u + Ku u u) + I v v v T g g g(u u u, u u u) L(v v v) = I v v v T f f f (9)
We define the two applications S S S m : q q q → P such that Eq.(8a) holds (the space problem), and T T T m : P → q q q such that Eq.(8b) holds (the time problem). The couple (P m , q q q m (t)) is given when stationarity has been reached, i.e when P m = S S S m (q q q m (t)) and q q q m (t) = T T T m (P m ).

When applied to problem in Eq.( 1), this procedure reduces to perform Galerkin projections in time and in space alternatively. The space problem corresponds to solve equation Eq.( 4) and the time problem corresponds to solve Eq.( 6).

Algorithm 1 gives the algorithm corresponding to oPGD. At each step k, the computation of the couple (P m , q q q m (t)) requires solving nm nonlinear algebraic equation (space problem in Eq.( 4)) and m nonlinear differential equations (time problem in Eq.( 6)). The time problem is quite small (only m equations) and can solved by HBM efficiently. The resolution of the space problem can become very costly when dealing with large number of

Algorithm 1 ALGORITHM FOR OPTIMAL PGD

for m = 1 to m max do Initialize q q q m (t) for k = 1 to k max do Compute P m = S S S m (q q q m m m ) by solving Eq.( 4) Orthonormalize P m [optional] Compute q q q m (t) = T T T m (P m ) by solving Eq.( 6) end for Set u u u m (t) = P m q q q m (t) and check convergence end for dof n. To improve performances, solutions of the space problem can be searched for in a smaller space, for instance by searching solutions as a linear combination of linear modes shapes [START_REF] Chinesta | Non incermental strategies for simulating thermomechanical models with uncertainty[END_REF].

Checking convergence is done by monitoring the relative change ε between two solutions given by ε = u u u m+1u u u m / u u u m .

Progressive PGD During the progressive PGD, solutions are computed iteratively. At step m, the previously computed solution u u u m-1 is completed with a new term p p p m q m (t) such that u u u m (t) = u u u m-1 (t) + p p p m q m (t) with u u u m-1 = ∑ m-1 k=1 p p p k q k (t) = P m-1 q q q m-1 (t). The vector p p p and the time evolution q(t) are computed such that the following orthogonality criterion hold simultaneously:

B(u u u m , p p p * q(t)) = L(p p p * q(t)), ∀p p p * ∈ E pPGD S B(u u u m , p p pq * (t)) = L(p p pq * (t)), ∀q * (t) ∈ E pPGD T [START_REF] Sanliturk | Modelling two-dimensional friction contact and its application using harmonic balance method[END_REF] where E pPGD S is the set of real-valued vector of size n, E pPGD T is the set of periodic function of period T = 2π/ω, B(u u u, v v v) and L(v v v) are given in Eq.( 9). We define the two applications S m-1 : q → p p p such that Eq.(10.1) holds (the space problem), and T m-1 : p p p → q such that Eq.(10.2) holds (the time problem). The couple (p p p m , q m (t)) is given when stationarity has been reached, i.e when p p p m = S m (q m (t)) and q m (t) = T m (p p p m ).

This procedure again results in Galerkin projections, but the projection is carried out on only one vector of the basis (either p p p m or q m (t)), leading to the following equation for the time (Eq.( 11)) and the space problem (Eq.( 12)): p T m (Mp p p m qm + Cp p p m qm + Kp p p m q m ) + g g g(Pq q q + p p p m q m ) = p T m f f fp T m (MP q q q + CP q q q + KPq q q) (11) I q m qm M + I q m qm C + I q m q m K]p p p m + I q m g g g(Pq q q + p p p m q m ) =

I q m f f f - m-1 ∑ k=1 [ I q m qk M + I q m qk C + I q m q k K]p p p k ( 12
)
Algorithm 2 ALGORITHM FOR PROGRESSIVE PGD for m = 1 to m max do Initialize q m (t) for k = 1 to k max do Compute p p p m = S m (q m ) by solving Eq.( 11) Orthonormalize p p p m relative to p p p 1 , . .

. , p p p m-1 [optional]

Compute q m (t) = T m (p p p m ) by solving Eq.( 12) end for if update then Update time functions q q q m = T T T m (P m ) by solving Eq.( 6) Set u u u m (t) = P m q q q m (t) and check convergence else Set u u u m = u u u m-1 + p p p m q m (t) and check convergence end if end for

The algorithm used to implement pPGD is depicted in Algorithm 2. At each step m the computation of the new term to be added consist in solving k max times an algebraic problem with n unknowns and a differential equation with one variable, where k max is the number of iteration needed to achieve stationarity. Since we are searching for periodic solutions, HBM is used to solve time problems in Eq. [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF]. The number of pPGD step m max can be determined by monitoring the relative change in norm between two solutions (as in oPGD), or can be set a priori by the user. The results of the pPGD can be improved by updating the whole set of time functions q q q after a new couple (p p p m , q m ) has been computed. This update is done by solving the time problem associated with the mapping T T T : P → q q q(P), i.e. by solving Eq. [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF]. pPGD can also be improved by computing several term at the same times. This method is based on the progressive PGD, but instead of computing one term (p p p m , q m ) at each iteration we choose to compute r terms simultaneously (with r remaining relatively small) by solving spacial and time problems derived from oPGD.

Remarks about the time and space problems

In all PGD methods described above, time problems (Eqs.( 6) or ( 12)) are solved by HBM with H harmonics in order to find steady 6) or ( 12)) is transformed into a system of nonlinear algebraic equation having the form of Eq.( 7) Since time problems only consists in a small set of differential equations of size m (m ≤ n), a large number of harmonics H can be retained, leading to solve a system of algebraic equations with m(2H + 1) unknowns at each iteration of oPGD (or at each update of pPGD).

When building space problems (Eqs.( 4) or ( 11)) integrals involving q(t), q(t) and q(t) are easily computed in the frequency domain. Indeed, the solution of time problems by HBM gives q i (t) = T(t) q q q i where T(t) contains the 2H + 1 vectors of Fourier basis. Integrals can then be approximated by:

I k = I q i (t) dq k j dt k (t)dt ≈ q q q T i H(D k q q q j ) (13)
where

H = I T T T T (t)T T T (t)dt is a 2H + 1 × 2H + 1 matrix.

NUMERICAL APPLICATION Beam featuring geometric nonlinearities

The system considered here consists in an Euler-Bernoulli beam as depicted in Fig. 1. Since the beam is clamped/clamped, interactions can occur between axial and transverse displacements leading to geometric nonlinearities. The numerical values of physical and geometrical parameters are given in Table 1. The beam is discretized by mean of finite element method with three dof per node (axial displacement u, transverse displacement v, rotation θ ) and n e = 20 elements.

FIGURE 2. DISPLACEMENT OF THE CENTER OF THE BEAM AS A FUNCTION OF FREQUENCY FOR THE REFERENCE SOLU-TION

After assembling and adding a damping term C = 3M, the resulting vibration equation for the beam is in the form of Eq.( 1) with n = 60 dof and with the nonlinear term g g g arising from geometric nonlinearities given by Eq.( 2).

Forced Solution

In this section we consider that the beam is excited at its center by an purely transverse harmonic force f (t) = A cos(ωt) with A = 200N. The reference solution u u u re f is computed by HBM with H = 3 harmonics for ω 2π from 155 to 180 Hz (upward direction) and for ω 2π from 190 to 162 Hz(downward direction) see Fig. 2.

For POD, at each frequency step, a POD is computed from the reference solution, and an approximated solution u u u m POD is computed by projecting the system onto a basis with m elements (this ideal version of the POD have no practical interests but for comparison with other methods).

For optimal PGD the initialization of q q q ω 1 m (t) (Alg.1) is chosen randomly for the first frequency step, and for all subsequent steps it is initialized with the previously computed value q q q ω i-1 m . The same strategy is used for progressive PGD (Alg.2) along with the update step each time a new vector has been computed. Solutions computed by oPGD (resp.pPGD) with m terms are denoted u u u m oPGD (resp.u u u m pPGD ). Since no theoretical results are available for proving the convergence of our fixed point problem, we did some preliminary studies about the numerical convergence of the fixed point problem. It turns out that the convergence is actually very fast and can be achieve most of the time within 3 iterations. Consequently, we choose to set the number of maxi-mum iterations to k max = 3. Finally the threshold for the relative change is set to ε = 0.005.

Solution are compared relative to the reference solution by using the relative difference ε r = u u u re f (t)u u u(t) / u u u re f . Figure 3 depicts the error relative to the reference solution for various number of retained mode m ≤ 4 for the three methods POD, oPGD and pPGD. For POD the error is clearly decreasing each time a new mode is retained. For oPGD we can see that increasing the number of maximum mode does not always decrease the error since the number of mode is actually controlled by the relative change in solution ε. The same remark holds for progressive PGD. We can see that oPGD gives very similar result than POD in term of error, and can sometimes provides better result than POD for a fixed number of modes, for example oPGD with m = 2 produces a smaller error than POD with the same number of mode (this is also happening for m = 4 see Fig. 3) (a similar fact can be observed for pPGD with m = 2 which produces a smaller error than POD with m = 2). We believe this phenomenon is due to the fact that oPGD introduces axial displacements sooner that POD (recall that nonlinearity arise from axial/transverse coupling), axial displacements are present in the second mode of oPGD and in the third mode of POD (see Figs. [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF][START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF].

In term of performance, solving the spatial problem is very costly since one need to solve a nonlinear algebraic set of n × m equations with Newton Raphson method. One way to improve performance would be to consider an incremental or a Newton linearization of the space problem as proposed in [START_REF] Chinesta | Recent advance and new challenges in the use of the proper generalized decomposition for solving multi dimensional models[END_REF]. Another mean of improving performance is to seek the solution of the spatial problem in a smaller subspace as explained in [START_REF] Chinesta | Non incermental strategies for simulating thermomechanical models with uncertainty[END_REF].

Free Solutions

In this section we search for free solutions of Eq.( 1) i.e. nonlinear normal modes (NNM). We only concentrate on the first mode of vibration but this procedure can be applied for any mode. A reference solution is computed by HBM with H = 3 harmonics and with cosine term only, this last condition corresponds to a phase condition in which the velocities are set to zeros for t = 0. As in the forced case, a POD basis with m components is computed from the reference solution for each frequency step, and an approximated solution is obtained by projecting the system onto this basis. Finally, a PGD approximation of the NNM is computed by optimal PGD. For the first frequency step, the solution is initialized with the linear mode and then the amplitude of the time evolution of the first mode is controlled in order to avoid null solution. The parameter used for this computation with oPGD are k max = 3 and ε = 0.005.

Figure 7 shows the backbone curve of the first mode obtained for the various methods (HBM,POD,oPGD) and Fig. 8 show the relative error with respect to the reference solution. Once again we see that oPGD gives results very close to POD when using only one mode in the decomposition. For m = 2 we can see that oPGD produces a smaller error than POD, and that this solution is very close to the reference one.

CONCLUSION

In this paper proper generalized decomposition was applied to solve nonlinear vibration problem. During PGD the solution of the space problems was found by using the harmonic balance method thus allowing to search for periodic solutions. Several algorithm were applied through a numerical example consisting in beam featuring geometric nonlinearity, modeled by finite element method. Comparison were drawn between solution computed by HBM, POD and PGD. Results show that optimal Galerkin gives very similar result than POD and can sometimes produce a smaller error than POD for the same number of mode in the solution. Although PGD gives acceptable results, work has to be done on performance and particularly on solving the space problem which can be very costly for large dof systems.

FIGURE 1 .

 1 FIGURE 1. SCHEME OF THE CLAMPED/CLAMPED BEAM AND REFERENCE ELEMENT

FIGURE 3 .

 3 FIGURE 3. RELATIVE ERROR ε r FROM THE REFERENCE SO-LUTION FOR POD, oPGD AND pPGD FOR m = 1, 2, 3, 4

FIGURE 4 .

 4 FIGURE 4. DECOMPOSITION OF THE SOLUTION COMPUTED BY oPGD FOR m = 4 AT ω 2π = 180 (shape have been normed for the purpose of comparison)

FIGURE 5 .FIGURE 6 .

 56 FIGURE 5. DECOMPOSITION OF THE SOLUTION COMPUTED WITH pPGD FOR m = 4 AT ω 2π = 180

FIGURE 7 .FIGURE 8 .

 78 FIGURE 7. BACKBONE CURVE OF THE FIRST NONLINEAR MODE COMPUTED BY HBM, POD and oPGD

TABLE 1 .

 1 NUMERICAL PARAMETERS USED IN THE BEAM MODEL

	L (m)	h (m)	ES (Pa.m 2 ) EI (Pa.m 4 ) ρS (kg.m -1 )
	1	3 • 10 -2	1.89 • 10 8	1.41 • 10 4	7.02
	state solutions. This mean that the set of differential equations
	defining the time problem (Eqs.(

ACKNOWLEDGMENT

Thanks go to SNECMA for its technical and financial support.