Benjamin Trévisan
email: benjamin.trevisan@insa-lyon.fr

Kerem Ege

Bernard Laulagnet

Vibroacoustics of orthotropic plates ribbed in both directions: application to stiffened rectangular wood panels

published or not. The documents may come

Vibroacoustics of orthotropic plates ribbed in both directions: Application to stiffened rectangular wood panels

I. Introduction

With their reduced weight and their high rigidity, stiffened panels are increasingly used in many industrial domains as aircraft manufacturing, naval (submarine shells) and automotive industries (rear and side walls of trucks cabs) but also in musical instrument making (piano or guitar soundboards for example). Contrary to more simple structures like beams or plates for which the vibro-acoustic behavior is now well described 1;2;3;4;5;6;7;8;9;10 , the design and modeling of ribbed structure with optimized vibroacoustic performances remains a major industrial challenge and a hot topic of research.

Since last decades several studies in various domains have concerned the vibroacoustic behavior of ribbed structures; most of them have been interested in periodically isotropic ribbed flat plates 11;12;13;14;15;16;17;18;19;21 , curved panels [START_REF] Chang | Noise transmission into semi cylindrical enclosures through discretely stiffened curved panels[END_REF] , stiffened cylindrical shells 23;24;25;26;27;28;29;30 . Nowadays, the development of composite materials offers new perspectives of improvements and a resurgence of interest for these structures. Moreover, the actual numerical resources make possible to describe many kinds of geometry and to include new material aspects (orthotropy by example) in order to have the more accurate results (see [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] for example). However, the effects of the broadband vibroacoustic behavior of these structures is not well known yet. In addition of the geometrical aspects, the current studies include the special characteristics of materials like interactions between the B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 4 different layers in the case of laminated composite panels 32;33 or the influence of air cavities on wave propagation in the case of sandwich structures 34;35;36;37 . On the particular case of ribbed wood panels, the pseudo-periodic parallel stiffening creates many localized phenomena in some specific frequency range (well described for the piano soundboard case 38;39;40 for example).

In order to understand the vibro-acoustic phenomena of this kind of ribbed structure, we develop a quasi-analytical model based on a variational approach and the judicious use of the associated orthotropic plate modes. We consider in this paper a simple geometry with: a thin flat plate made of wood (spruce) with a special orthotropy, pseudo-periodically ribbed in the direction of weak Young's modulus by eleven ribs and by one straight rib in the perpendicular direction (strong Young's modulus direction). The geometry is described in Fig. 1. This approach is an alternative to a complete and more accurate numerical modeling method, i.e the Finite Element Method for the structure and the Boundary Element Method for the acoustical radiation. In order to understand the vibroacoustic phenomena of stiffened structures, we assume that a method that avoids any discretization of the object would contribute to knowledge on the subject and would allow for easy parametric studies (influence of the ribs, their number, their width/length).

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 5

The first part of the article presents the quasi-analytical model and discusses modal aspects due to stiffening on the whole frequency range [0,5000] Hz. Then in a second part, the modal radiation coefficients of a ribbed structure that have never been shown before are given. Finally, in the last part of the article the numerical results show that the apparent critical frequency of a heterogeneous structure is very sensitive to the excitation point of the effort. In particular, when applied on a hard point of the structure, i.e. on the straight rib, it decreases the apparent critical frequency, selecting the more radiative modes.

II. Vibratory modeling A. Context of the study

In this study, we focus on a rectangular ribbed wood panel (see Fig. 1). We consider a special orthotropic plate (axis of orthoropy parallel to the edges, see Fig. 1). The boundary conditions are simply supported and the plate has one straight rib in the direction of the strong Young's modulus x and eleven pseudo-periodic smaller ribs in the perpendicular direction y. These superstructures are also made of spruce.

The eccentricity of the different superstructures is taken into account. Moreover, the thickness of the plate, width and height of the superstructures are constant along their length.

The analytical model developed is based on a variational approach that draws its B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 6 inspiration from Laulagnet & Guyader's works 23;41 . Moreover, the traversal displacement is developed on the basis of simply supported plate modes.

It is necessary to calculate the kinetic and strain energies of the whole structure to calculate the Hamilton function of the problem, and develop it on the basis mentioned above in order to write a matrix problem and, finally, determine a new specific basis for the ribbed structure.

B. Hamilton function for eleven stiffeners and one straight rib

The modeling is based on a variational approach. So, we want to express the Hamilton functional of the whole problem, which is defined as:

H = t 1 t 0 (E c -E d + E w ) dt (1)
with [t 0 ; t 1 ] an arbitrary time interval, E c and E d the kinetic and strain energies and E w , the work done by the external force or sound pressure. In order to solve the eigenvalue problem, E w is equal to zero for the moment.

Kinetic hypothesis and boundary conditions at the interface plate / superstructures

Because the thickness of the plate is very small compared to the other dimensions, the Love-Kirchhoff hypothesis referring to thin plates was adopted. Thus pumping effects and shearing in the two bending planes are ignored Taking into account these hypotheses and writing the motion as linear into the thickness (limited development at order 1), leads to the following motion field:

               u(x, y, z, t) = -z ∂w ∂x v(x, y, z, t) = -z ∂w ∂y w(x, y, z, t) = w(x, y, t) (2)
with z ∈ -h 2 ; h 2 (h being the plate thickness) and where u(x, y, z, t) and v(x, y, z, t) are the motions in directions x and y and w(x, y, t) is the plate transversal displacement.

The different superstructures are considered as beams driven by bending and torsion coupled to the plate. Nevertheless, superstructures cannot be approached as simple beams due to their coupling with the plate. Indeed, boundary conditions are applied to the plate, thus its motion governs all the elements attached to it, as shown in Fig. 2 which illustrates this phenomenon in bending and torsion motions.

Therefore displacements and rotations are considered as continuous at each plate / stiffener interface. Moreover, the negligible shearing effects of the Love-Kirchhof plate theory are propagated to the ribs through the continuities mentioned above. Likewise for the other kinetic hypothesis, like the pumping effect and shell assumptions applied to the B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 8 displacement field of the superstructure. This hypothesis is in agreement with the idea of considering the ribs and straight rib as beams.

After performing calculus, the following displacement field is obtained for ribs oriented in direction y at position x = x r :

               u r (y, z, t) = -z w ,x (x r , y, t) v r (y, z, t) = -z w ,y (x r , y, t) w r (x, y, t) = w(x r , y, t) + (x -x r ) w ,x (x r , y, t) (3)
where x ∈ x r -b 2 ; x r + b 2 and z ∈ h 2 ; h 2 + H and where h and H are respectively the plate thickness and rib height, and b its width.

Expression of the Hamilton function as a function of the kinetic field of the plate

We express the Hamilton functions of each sub-part on the basis of the different fields of motion. These functions are expressed by the integral of the difference between the kinetic and strain energies on an arbitrary time interval (see Appendix I for more information on this calculus).

Thus we write successively the actions of the plate, of the bending and torsion of a rib in direction y at position x = x r and of the bending and torsion of the straight rib in direction B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 9

x at position y = y c by:

H p = t 1 t 0 1 2 l x=0 L y=0 ρh ẇ2 -D 1 w 2 ,xx + D 3 w 2 ,yy + D 2 w ,xx w ,yy + D 4 w 2 ,xy dxdydt (4)
H rb = t 1 t 0 1 2 l x=0 L y=0 ρ r I f ẇ2 ,y + bH ẇ2 -E r I f w 2 ,yy δ(x -x r ) dxdydt (5)
H rt = t 1 t 0 1 2 l x=0 L y=0 ρ r I g w 2 ,x -G r I g w 2 ,xy δ(x -x r ) dxdydt (6)
H bb = t 1 t 0 1 2 l x=0 L y=0 ρ c I f c ẇ2 ,x + b c H c ẇ2 . -E c I f c w 2 ,xx δ(y -y c ) dxdydt (7)
H bt = t 1 t 0 1 2 l x=0 L y=0 ρ c I gc w 2 ,y -G c I gc w 2 ,xy δ(y -y c ) dxdydt (8)
where ρ and ρ r are respectively the mass density of the plate and the rib ; h, l and L the thickness, the length and the width of the plate ; The Table 1 recaps the different constants used in the previous equations.

D 1 = Exh 3 12(1-vxyvyx) , D 3 = Eyh 3 12(1-vxyvyx) , D 2 = vyxExh
Note that all these functions are defined on the intervals x ∈ [0; l] and y ∈ [0; L]. This formulation does not mean that the superstructures occupy the whole area of the plate.

Indeed, they are considered as punctual in their widths, justifying the Dirac distributions

Because the rib height is low compared to its width, the torsional inertia does not take warping into account. Moreover, because the plate controls the rib motion, the action of the whole system is expressed as a function of the plates transversal displacement and its spatial and temporal derivatives.

Finally, for a structure with a determined number N r of ribs in direction y and a straight rib in direction x, we express the entire function by:

H System = H p + Nr i=1 H (i) rb + H (i) rt + (H bb + H bt ) (9)

C. Hamilton function as a function of modal amplitudes of unribbed plate modes

We develop these energies on the basis of simply supported unribbed plate modes.

This basis, currently used in the area of vibrations 23;41 is particularly adapted for an analytical approach. Moreover, the radiation of the plates was studied specifically for these boundary conditions.

Modal decomposition on the basis of simply supported plate modes

Through modal decomposition, the plate transversal displacement is written as a B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 11 linear combination of unribbed plate modes weighted by modal amplitudes a mn (t):

w(x, y, t) = M m=1 N n=1 a mn (t)φ mn (x, y) ∀x ∈ [0; l] & y ∈ [0; L] (10)
where φ mn (x, y) = sin mπ l x sin nπ L y .

It should be recalled that these modal shapes have orthogonal properties:

l 0 L 0 φ mn (x, y)φ pq (x, y) dx dy =        l * L 4 if m=p & n=q 0 if m =p || n =q (11)
The modal decomposition (10) is finite; there are truncations on the orders (M,N).

The further the problem is from the initial basis, meaning without superstructures, the higher the order of the truncations must be to guarantee the convergence of the numerical results. In particular, we know that the strong inter-modal coupling introduced by the superstructures will lead us to high orders (M,N).

New Hamilton function as a function of the modal amplitudes of simply supported plate modes

By injecting the modal decomposition Eq. ( 10) into the Hamilton function expressions Eq. ( 4) to (8) and using the orthogonal properties Eq. ( 11), it is possible to analytically calculate the surface integral in the Hamilton function.

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 12

In this case, the functional depends on the couples of variables (a mn (t), ȧmn (t),a pq (t), ȧpq (t)) and no longer on the transversal displacement w(x, y, t) and its space and temporal derivatives w ,x , w ,xx , w ,yy , w ,xy and ẇ. Thus, we have:

H System (a mn (t), ȧmn (t), a pq (t), ȧpq (t)) = t 1 t 0 L(a mn (t), ȧmn (t), a pq (t), ȧpq (t)) dt (12)
where L(a mn (t), ȧmn (t), a pq (t), ȧpq (t)) is termed the Lagrangian of the system.

Afterwards, the action of the system is minimized following the principle of less action. In practice, the differential form of Euler-Lagrange is used to determine the evolution of the system, leading to the following equation:

δH System = 0 ⇔ ∂L ∂a pq - d dt ∂L ∂ ȧpq = 0 ( 13
)
express for each mode of the unribbed plate and where p = 1 → M and q = 1 → N .

This minimization leads us to a single equation for a particular mode (p,q):

M m=1 N n=1 M(m, n, p, q).ä pq + M m=1 N n=1 K(m, n, p, q).a pq = 0 ( 14
)
where M(m, n, p, q) and K(m, n, p, q) are respectively mass and rigidity matrix terms depending on integers (m, n, p, q). The sums on the modal amplitude integers highlight the inter-modal coupling introduced by the superstructures. However, this expression does not take into account the external effort at this stage of the development.

Expressing the problem in this way makes it easier to show the coupling phenomena due to the presence of the superstructures. Finally, we obtain a homogenized problem for which the size is conditioned by an adapted truncation guaranteeing the convergence of the solution. Thus we obtain the following matrix formulation:

[M plate r ] + [M ribs rs + [M straight rib rs ] {ä p } + [K plate r ] + [K ribs rs + [K straight rib rs ] {a p } = 0 ( 15
)
where r = (m, n) and s = (p, q) and where M and K represent mass and rigidity matrices.

The system ( 15) is constituted by the sum of diagonal matrices for the plate and of full and symmetrical matrices for the superstructures. These introduce very strong coupling of the unribbed plate modes, as we will see in the following sections. Indeed, these matrices do not have a dominant diagonal.

Consequently, the size of the matrix is also important for the convergence and precision of the results. Indeed, since the wood panel studied is an heterogeneous structure, its modes are far from the initial basis (simply supported plate). Thus it is necessary to take into account modes with small wavelengths (high order modes) in order to describe its vibrational behavior correctly. In the present case, this leads us to a truncation on the B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 14 couple (p,q) around (40,30), i.e. 1200 simply supported plate modes. In the following, when the new basis of the ribbed structure is obtained, these 1200 modes will obtain around 350 panel modes in the frequency range [0;5000] Hz.

D. Calculus of ribbed panel modes

Ribbed modes are not as simple as simply supported plate modes, as we will see in the following. The first step to understand the behavior of such a structure is to determine and thus calculate the new basis specific to the panel.

Eigenvalue problem

In order to determine the new ribbed basis, we are interested in the canonical form of the eigenvalue problem:

N -λ Ī ā = 0 ( 16
)
where N is the dynamical rigidity matrix, also called core of the system 42 , with N = M -1 . K, where K and M are the mass and rigidity matrices of the system and with λ

the eigenvalues : λ = (2πf ) 2 .
The search for the eigenvalues and eigen vectors of the N matrix solving Eq. ( 16) leads to a diagonal matrix called D with square angular frequencies and to a matrix of eigen vectors called P whose terms are the modal amplitudes weighted by the eigen modes of the unribbed plate used to re-create ribbed plate modes.

D =             (2πf 1 ) 2 (0) (2πf 2 ) 2 . . . (0) (2πf i ) 2             ; P =             a (1
P Q a (2) P Q a (P * Q) P Q             (17)

Ribbed plate modes

Each ribbed plate mode is reconstituted by a linear combination of simply supported plate modes. Thus for each eigen frequency f i , we calculate the corresponding modal shape by:

φ (i) (x, y) = M m=1 N n=1 a (i) mn φ mn (x, y) ∀x ∈ [0; l] & y ∈ [0; L] (18)

Numerical results: panel modes

In the aim of applying in the future the model presented to musical instruments, and in particular to the piano soundboard, dimensions of superstructures are inspired from [START_REF] Ege | Vibroacoustics of the piano soundboard : (Non)linearity and modal properties in the low-and mid-frequency ranges[END_REF] (see the perspectives in the last section of the article).For this study, the frequency range of interest is [0;5000] Hz. The modal shapes are calculated for an order (p,q) equal to (41,27) to ensure the convergence of the solution.

The presence of only one straight rib in the direction x seems to separate the plate into two parts. A vibrating area and a non-vibrating area is often found on each side of this one (3 th , 4 th and 5 th mode in Fig. 3).

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 16

Under these conditions, only the two first modes seem to be similar to the homogeneous unribbed plate modes. The separating effect appears from the third mode (242 Hz), splitting the plate into 2 areas; for each side the vibration is approximately homogeneous.

With this consideration, the ribbed structure appears like an equivalent homogeneous unribbed plate but only over a reduced area delimited by the straight rib.

Similarly, the regular addition of stiffeners causes localized vibrations. However, due to these several additional ribs, the split areas are smaller than those with the straight rib along x and this phenomenon does not appear in the first modes. Thus Fig. 3 shows that when the frequency increases, modes can be seen for which the vibration is delimited over a small surface, for example the 17 th mode. At high frequency, the modes become very complex. A large number of localized vibratory phenomena can be seen such as oblique waves and wave reflections in the inter-rib spaces. See Fig. 3.

III. Acoustic radiation of ribbed panel

The behavior of such a structure is not limited to a purely vibratory aspect. Its radiation, which is complex due to its particular geometry, has rarely been studied with the Finite Element Method, due to high numerical and temporal costs (Chabassier & al. 31 in the case of a grand piano soundboard). This is why the analytical model presented here is a good alternative in order to perform parametric studies quickly and with minimum computer resources.

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 17

In the modeling, the baffled hypothesis is made. In real situation, a panel radiates in the two half spaces (up and down) which naturally communicate, tending at low frequency to considerable acoustical short cuts. Indeed, the radiation becomes dipolar as the wavelengths are large compared to the size of the structure. Thus it is better to assume no baffle in this frequency range. The influence of the baffle on the acoustic radiation at low frequency can be found in 3;9 .

Nevertheless, when focusing on a frequency higher than the first octave, the acoustic wavelengths fall and become smaller than the size of the panel and the baffled hypothesis becomes solid. Therefore only the radiation on one face is calculated.

We are interested in the acoustic radiation of the ribbed wood panel and more particularly in the general indicators such as the acoustical radiated power W (ω), the space average quadratic velocity < v 2 > (ω) and the radiation coefficient σ(ω). We invite the reader to refer to the Table 2 to read the rest of this paper.

A. Vibroacoustic problem: forced response

In the following, we considered a sinusoidal and harmonic excitation applied at (x e ,y e ) with an amplitude F of 1N. So, the Hamilton function now includes the work done by this external force minimized following the same approach that in the previous parts.

with a pq an arbitrary mode of the series. Thus, in the frequency domain, the problem with a second member and the vector of modal amplitudes can be written with the following expressions:

K * -ω 2 M ā = Fgen ( 20
)
⇔ ā = K * -ω 2 M -1 Fgen (21)
Note that the light fluid assumption is made. Therefore the Eq. (20 & 21) does not take into account the acoustical boundary pressure that we assume to be negligible at the first order. For a full calculus, see [START_REF] Laulagnet | Sound radiation by a simply supported unbaffled plate[END_REF] .

Moreover the structural damping is taken into account by making the stiffness matrix K * = K * (1 + jη) as a complex number with η = 2% for all the frequencies.

B. Expression of the ribbed panel radiated power

In the case of a harmonic problem, after application of the modal decomposition (10), we give the following expression of the radiated power [START_REF] Lesueur | Acoustic Radiation of Structures, vibroacoustic, fluid/structure interactions[END_REF] :

W (ω) = 1 2 -jω p q a * pq (ω)P a pq (ω) ( 22
)
where a * pq (ω) is the conjugated modal amplitude of a mode (p,q) and P a pq (ω) the B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 19 generalized boundary pressure defined by:

P a pq (ω) = Sp φ pq (x, y)p(x, y, z = 0) dS (23)
with p(x, y, z = 0), the boundary pressure.

Solving the Helmoltz equation in the wave number domain 3;43 and applying the boundary conditions (baffled plate), we obtain an expression of the pressure in the wave number domain:

p(k x , k y , z) = jρ 0 ω 2 ũz (k x , k y , z = 0) k z e -j.kz.z avec k x , k y ∈ 2 (24)
where

k z = k 2 -k 2 x -k 2
y is the reduced wave number in the direction z and ũz (k x , k y , z = 0) the particular acoustical displacement Fourier transform of the plate.

With the series development of the plate Eq. ( 10), we obtain:

ũz (k x , k y , z = 0) = m n a mn φmn (k x , k y ) (25)
where φmn (k x , k y ) refers to the Fourier transform of a simply supported mode. Returning to the spatial domain, using an inverse Fourier transform leads to the following expression of the pressure in the fluid:

p(x, y, z) = 1 4π 2 +∞ -∞ jρ 0 ω 2 k z m n
a mn φmn (k x , k y )e +jkxx e +jkyy e -jkzz dk x dk y

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 20

In this expression, note the undefined double integral in the wave number domain.

Therefore it is estimated numerically [START_REF] Laulagnet | Sound radiation by a simply supported unbaffled plate[END_REF] . By integrating this expression of the pressure in Eq. ( 23), we define the generalized acoustic pressure:

P a pq = jω m n a mn Z mnpq (ω) (27)
where Z mnpq (ω) is the inter-modal radiation impedance of modes (m,n) and (p,q) with:

Z mnpq (ω) = ρ 0 ω 4π 2 +∞ -∞ φmn φ * pq k 2 -k 2 x -k 2 y dk x dk y ( 28
)
where * is the conjugate complex number.

Neglecting inter-modal coupling, i.e. Z mnpq = 0 if nm = pq, thus the classical light fluid assumption (effects of inter-modal coupling can be consulted in 8;18;32;44 ), we can define a new expression of the acoustic radiated power of the panel in the frequency domain as [START_REF] Lesueur | Acoustic Radiation of Structures, vibroacoustic, fluid/structure interactions[END_REF] :

W (ω) = ω 2 2 m n |a mn | 2 R mnmn (ω) (29)
where R mnmn (ω) = (Z mnmn (ω)). The R mnmn (ω) calculus is performed numerically. We show that 3;6 :

R mnmn (ω) = ρ 0 ω 4π 2 k -k | φmn | 2 k 2 -k 2 x -k 2 y dk x dk y (30)
C. Expressions of the space average quadratic velocity < v 2 > (ω) and of the radiation coefficient σ(ω)

B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 21

The space average quadratic velocity expression is given by 6 :

< v 2 (ω) >= 1 2S l x=0 L y=0 |v(x, y, ω) 2 | dxdy (31)
where S = l * L is the plate surface. Applying the modal decomposition Eq. ( 10) and using the orthogonal properties of the modal shapes in the simply supported basis, leads to the following expression of the space average quadratic velocity in the frequency domain:

< v 2 (ω) >= ω 2 8 m n |a mn (ω)| 2 (32)
Finally, we define the radiation coefficient by [START_REF] Lesueur | Acoustic Radiation of Structures, vibroacoustic, fluid/structure interactions[END_REF] :

σ(ω) = 1 ρ 0 cS W (ω) < v 2 > (ω) (33)
where ρ 0 and c are respectively the mass density and the celerity of sound.

D. Expression of the radiated power in the panel basis

Calculating the radiated power in the ribbed modes basis is not necessary but allows defining the modal radiation coefficient of a ribbed structure, which has never been presented in the literature before. We give the expression of the power Eq. ( 22) in scalar form:

W (ω) = 1 2 -jω p q a * pq (ω)P a pq (ω) ⇔ W (ω) = 1 2 -jω ā * T (ω). Pa (34)
B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 22

where ā * T is the transposed vector of modal amplitudes taken from Eq. ( 21) and Pa is the vector of generalized boundary pressure that we defined in the same manner as:

P a pq = jω m n a mn (ω)Z mnpq (ω) ⇒ Pa = jω Za (ω) a(ω) ( 35
)
where Za is the inter-modal acoustic impedance matrix whose terms are Z mnpq (ω). This leads us to a compact matrix expression of the radiated power. Thus we obtain:

W (ω) = ω 2 2 ā * T (ω). Ra (ω).ā(ω) (36)
where Ra (ω) is the real part of the simply supported plate inter-modal impedance matrix and ā(ω) the vector of modal amplitudes for an excitation at a specific frequency obtained from solving Eq. ( 21).

As presented in the first part, solving the homogeneous problem Eq. ( 15) conferred the matrix of eigenvectors P , and consequently, the transfer matrix T to the new ribbed basis which is the orthonormed basis of eigenvectors. By definition, we have:

ā = T. b ( 37
)
where ā is the vector of modal amplitudes in the initial plate basis and b the same vector now in the ribbed basis. Therefore we give the acoustic radiated power expression in the ribbed basis:

W (ω) = ω 2 2 b * T (ω). T T . Ra (ω). T. b(ω) (38)
Rs (ω) = T T . Ra (ω). T

E. Radiation coefficient definition of a ribbed mode

We explained in the introduction that we neglected inter-modal coupling between simply supported plate modes, which means considering zero energy exchanges between two modes of the simply supported plate. From the matrix point of view, the inter-modal impedance matrix Ra will be assumed to be diagonal.

This simplification allows easily calculating the radiation of the structure in a light fluid and introducing the radiation coefficient of a particular mode (p,q). We now want to determine if this hypothesis can be generalized to a ribbed structure, at least in a frequency range.

This raises the question of whether extra-diagonal terms (inter-modal coupling) are negligible, which implies a dominant diagonal matrix. Indeed, diagonal terms represent the interaction of a ribbed mode with itself whereas extra-diagonal terms represent energy exchanges between modes.

F. Inter-modal coupling in the acoustic radiation of a ribbed structure

In this part, we focus on the normalized radiation impedance 8;11 r(ω) = R a||s (ω) interest in unstiffened coefficients or stiffened coefficients. The results will be presented as 10.log(|r ij (ω)|). When i = j, we obtain the real radiation impedance part of a ribbed panel mode which is what we will study in the following. When i = j, we talk about inter-modal radiation impedance between modes, an impedance which can be negative, contrary to the case i = j.

ρ 0 c S
In the same way as for the radiation impedance of simply supported plate modes, it is common to verify if interaction terms are negligible compared to diagonal terms:

r mnmn (ω) >> |r mnpq (ω)| (40)
Indeed, the radiation impedance notion of a mode (m,n) is meaningful only if these interaction terms (m,n),(p,q) between two modes of simply supported plate are negligible 3;6;8 .

Fig. 4-a and 4-c present the modal radiation coefficients of 1 st and 5 th ribbed panel modes (see Fig. 1) as well as examples of inter-modal coupling. Fig. 4-a shows that for the first mode, the extra-diagonal terms (dashed plots) are far below the diagonal term (full plot) with differences exceeding 10 dB.

On the other hand, for high order modes, we quickly see that energy exchanges are strong.

Even if it is negligibly higher than the coincidence frequency (here around 700/800 Hz), Fig. 4-c shows that these energy exchanges (dashed plots) have the same order of presented in Fig. 1 for which we remove the straight rib in the direction x. Thus the plate is pseudo-periodically ribbed in the direction x and each rib has a different height and width with mean values around 19 and 27mm [START_REF] Ege | Vibroacoustics of the piano soundboard : (Non)linearity and modal properties in the low-and mid-frequency ranges[END_REF] . Therefore the structure is quasi-homogeneous but where small variations of dimensions can lead to considerable asymmetry on the vibration 45;46 . This is the case of the two first modes. For a simply supported plate with the same dimensions, the first mode (1,1) presents a symmetry, contrary to our case study where the vibration is localized on the right half part (Fig. 5). This loss of symmetry implies a decrease of the modal radiation coefficient of the first ribbed mode compared to that of a simply supported plate mode (1,1) under the coincidence frequency. But both are close as shown by the small differences, with a maximum of 1.5 dB at 25 Hz.

The second ribbed plate mode is also strongly affected by small variations of rib dimensions. Although it resembles, at the first order, the mode (2,1) of a simply supported plate mode, we note a considerable loss of the anti-symmetry of the latter with different magnitudes at either side. Contrary to the first radiation coefficient, which is only weakly affected, this coefficient is very sensitive at low frequency, as shown in Fig. 6. Indeed, we see differences higher than 20 dB between the two coefficients at 25 Hz and, contrary to the first mode, this time the stiffeners increase the radiation coefficient. We assume that the B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 27 reduction of the vibrating anti-symmetry tends to make it more efficient at low frequency.

These first results show the influence of small variations of rib dimensions but they also show that it is impossible to predict if the radiation increases or decreases.

At a higher frequency, we see an alternation in most cases of quasi identical non ribbed plate modes and of modes with localized vibrations. Therefore most of the radiation coefficients are similar to those of a non-ribbed plate like the 5 th mode compared in Fig. 7 to the mode (5,1) of a rectangular plate with the same dimensions.

On the contrary, some modes present localized vibrations in inter-ribbed spaces. It would be pointless to compare them to the modes of an unribbed plate with the same dimensions.

For the 6th mode (Fig. 8), a comparison with the mode (1,2) of a rectangular plate reduced by 4 in the direction x slightly overestimates the radiation coefficient around 4 dB at 25 Hz but shows that for some modes, only a small area of the structure radiates.

Thus the radiation coefficient of a quasi homogeneous structure, pseudo-periodically ribbed in only one direction, can be close to those of an unribbed rectangular plate. In most cases, the modal shapes are similar to those of an unribbed plate and so the corresponding radiation coefficients are similar (Fig. 7). However, there are also localized modes for which the radiation coefficients are close to those of a reduced plate (Fig. 8). Finally, it is impossible to predict how the superstructures will modify the radiation because of the smallest modification of the geometry can have an considerable influence.

Increasing the frequency considerably changes the behavior of the ribbed plate in two directions in comparison to the previous case. Indeed, reduced homogeneous modal shapes are rare from the 7 th mode which is the last distinctly identifiable mode. Fig. 12 shows a comparison between its radiation coefficient and that of mode (2,2) of an unribbed plate with the same dimensions. Despite this relative agreement, this comparison shows that both coefficients are only comparable at a frequency higher than 100 Hz. Indeed, at low frequency, the differences become very significant, reaching 20 dB at 25 Hz.

For high order modes, the straight rib introduces geometric complexity, having considerable impact on the modal aspect of the system and making mode identification almost impossible. This is alternated with localized modes. The 10 th panel mode (Fig. 13) is a perfect example of this aspect. The vibration is split into two areas where we can see 4 maximum amplitudes of vibration below the straight rib and 3 above it. Consequently, its radiation coefficient is also unique and not comparable to the radiation coefficient of unribbed plate modes or that of the 271 th mode (Fig. 14) for which the vibration is extremely complex.

Finally, the high complexity of the vibratory aspect caused by structural heterogeneity is similar to that caused to the radiation. Thus the modes quickly become dependent on the geometry of the structure, removing the pertinence of making a comparison with a simple structure like the rectangular simply supported plate.

H. Influence of inter-modal coupling on the radiation coefficient σ(ω)

In view to providing an interpretation of ribbed mode radiation coefficients, we have assumed up to now that inter-modal coupling of the radiation impedance matrix was negligible. We now focus on the calculus of the radiation coefficient σ(ω) Eq. ( 33), plotted in 10.log(σ(ω)) for a sinusoidal excitation placed on the straight rib at the coordinates (x = 0.73 m ; y = 0.61 m). The structure modeled refers to Fig. 1, i.e. the ribbed wood panel.

Fig. 15 shows the influence of taking into account coupling terms. Globally, we observe very good agreement between the full and dashed plots, except for the frequency range [400;800] Hz. Indeed, in this range we tend to underestimate the radiation coefficient σ(ω) at a maximum of 2 dB.

Globally, this result shows the bias involved when couplings are assumed to be negligible. It is localized within a short frequency range and is rather weak.

IV. Numerical comparison between several excitation points

In this last part, we deal with the radiation of the same ribbed panel with different excitation point. Moreover, we now take into account the inter-modal coupling of the radiation impedance matrix in the ribbed basis. Since the radiation coefficient alone is B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 31 insufficient for understanding the influence of superstructures on radiation, we consider the acoustical radiated power, the space average quadratic velocity and the radiation coefficient Eq. ( 29), ( 32) and ( 33) in the following.

A. Influence of the excitation point: hard and weak point

For a ribbed structure, the excitation point has a considerable influence on acoustical radiation. To show this, we chose two characteristic points on the panel. These two points (see Fig. 16) are at the abscissa x e = 1m and at the ordinate:

• y e = 0.61m, on the largest superstructure (straight rib), so the hard point,

• y e = 0.47m, on the plate far from the first point, near the middle between the edge of the plate and the straight rib.

However, first let us define the critical frequency f c of a plate. At this frequency, the propagation velocity of flexural waves in the plate and in the fluid (air) is the same. Below this frequency, the bending wave is subsonic and poorly radiative; higher, it is supersonic and radiative.

For a simple case like an unribbed isotropic plate, we can easily determine the propagation velocity of the flexural wave in the structure. Thus we have the following expression of the critical frequency for an isotropic plate 2 :

f c = c 2 2πh 12ρ(1 -ν 2 ) E (41)
the plate mass density ρ, the Young's modulus E and the Poisson's coefficient ν of the plate.

For an orthotropic plate, there are two critical frequencies because of different propagation velocities due to the two Young's moduli. In our case, these frequencies are around 1450 Hz and 7250 Hz. For a structure with a complex geometry and superstructures, it is impossible to calculate f c . Nevertheless, we define the apparent critical frequency when the radiation coefficient is close to unity (0 dB in log scale).

Therefore there are two important results concerning the radiation coefficient (Fig. 17-c). Firstly, without mentioning the excitation point, we can see that the simple presence of the superstructures decreases the apparent critical frequency from 1450 Hz to less than 800 Hz. Thus the superstructures increase the radiation efficiency in the frequency range [800;1450] Hz. Obviously, above 800 Hz all the modes excited are radiative; on the contrary, at low frequencies below 400 Hz, the modes are weakly radiative, independently of the excitation point.

The second important result concerns the transitional frequency range [400;800] Hz which can be defined when the excitation position has an influence on the radiation efficiency.

Fig. 19-a shows that the radiated power is globally homogenous for an excitation on the straight rib. On the contrary, it is strongly affected for an excitation far from this superstructure (Fig. 19-b). As mentioned before, for a point placed directly on the plate, the average value of the acoustical radiated power is around 15 dB higher than for an excitation on the straight rib / bridge. However, we can also see that the value falls when the excitation moves to the ribs. This means that the response has a large number of discontinuities. Indeed, the ribs can be seen clearly in Fig. 19-b. These results show the importance of the straight rib / bridge in order to obtain the most homogenous acoustical response possible.

V. Conclusion

This article has presented the basis of a new analytical model for predicting the vibroacoustic behavior of orthotropic ribbed plates. Here, we have limited our study to the case of a wood panel: a special rectangular orthotropic plate with several ribs in one direction and one perpendicular straight rib.

The model developed allowed us to present the modes of a complex structure and focuses on the acoustic radiated power, the space average quadratic velocity and the radiation coefficients for different excitation points. It also allowed determining the influence of stiffeners on acoustic radiation.

Moreover we have introduced the notion of modal radiation coefficients for a ribbed B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 36 structure and have showed that it is meaningful. Namely the first modes of the stiffened panel exhibit radiation coefficients whose behaviour compares to the one of the unstiffened corresponding mode. On contrary, when the modal shape becomes complex, no similarity can be done with the radiation factor of existing modes and the tendency needs to be evaluated numerically, because of the specificity of the complex stiffened mode under consideration. One other conclusion is that the radiation coefficient of a stiffened panel increases sensibly, when the point excitation is located on a rib. Consequently, the apparent critical frequency decreases too in a sensible manner, rendering this critical frequency depending of the localization of the effort.

Acknowledgment

This work was performed within the framework of the Labex CeLyA of Université de Lyon, supervised by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).

The result is not immediate for the superstructures. We give the displacement field of a rib in direction y at position x = x r :

               u r (y, z, t) = -z w ,x (x r , y, t) v r (y, z, t) = -z w ,y (x r , y, t) w r (x, y, t) = w(x r , y, t) + (x -x r ) w ,x (x r , y, t) (45)
where x ∈ x r -b 2 ; x r + b 2 and z ∈ h 2 ; h 2 + H where h and H are respectively the plate thickness and the rib height.

Kinetic and strain energies are calculated on the entire rib volume placed at x = x r .

Let the kinetic energy be:

E C = 1 2 b 2 x=-b 2 L y=0 b 2 +H h= h 2 ρ r u2 rf t + v2 rf t + ẇ2 rf t dxdydz (46)
Replacing u2 rf t , v2 rf t and ẇ2 rf t by their expression from (45) and after making the integral on variables x and z, we have:

E C = 1 2 ρ r L y=0 I f ẇ2 ,y + bH ẇ2 + I g ẇ2 ,x dy (47)
where bending inertia is equal to

I f = b 3 h 2 + H 3 -h 2
3 and torsion inertia is

I g = Hb 3 12 + I f .
ḡrad( Ū ) + ḡrad( Ū ) T and the strain tensor ¯ = MS σ (Inverse Hook law). We find the following small deformation tensor and its vector form:

¯ = 1 2 0 -zw ,
The strains are linked to deformations through a flexibility matrix. The matrix MS for an orthotropic material is:

MS = 1 Ex -νxy Ex -νxz Ex -νyx Ey 1 Ey -νyz Ey (0) -νzx E Z -νzy E Z 1 Ez 1 2Gxy (0) 1 2Gxz 1 2Gyz (49)
B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 40 Thus we determine the three strains in the rib:

               σ yy = -zE y w ,yy (x r , y, t) σ xy = -zG xy w ,xy (x r , y, t) σ yz = xG yz w ,xy (x r , y, t) (50)
After which we calculate the strain energy:

E D = 1 2 b 2 x=-b 2 L y=0 b 2 +H h= h 2 σ : dxdydz (51)
with σ : = E y z 2 w 2 ,yy (x r , y, t) + (G xy z 2 + G yz x 2 ) w 2 ,xy (x r , y, t), let :

E D = 1 2 L y=0 E y I f w 2 ,yy (x r , y, t) + D g w 2 ,xy (x r , y, t) dy (52)
with 12 + G xy I f that we define as the torsion dynamic rigidity of an orthotropic beam in the plane ( x, z).

I f = b 3 h 2 + H 3 -h 2
For the sake of simplicity, we consider the beam as isotropic in the following. D g becomes D g = G Hb 3 12 + I f = G r I g . So, we express the Hamilton function of the rib H r :

H r = t 1 t 0 1 2 L y=0 ρ r I f ẇ2 ,y (x r , y, t) + bH ẇ2 (x r , y, t) + I g ẇ2 ,x (x r , y, t) -E r I f w 2 ,yy (x r , y, t) + G r I g w 2 ,xy (x r , y, t) dydt (53)
B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 41 that we gather on a line at x = x r :

H r = t 1 t 0 1 2 l x=0 L y=0 ρ r I f ẇ2 ,y (x, y, t) + bH ẇ2 (x, y, t) -E r I f w 2 ,yy (x, y, t) δ(x -x r ) dxdydt + t 1 t 0 1 2 l x=0 L y=0 ρ r I g ẇ2 ,x (x, y, t) + G r I g w 2 ,xy (x, y, t) δ(x -x r ) dxdydt = H rb + H rt (54)
For the straight rib, the functional H r is determined following the same approach so we obtain:

H b = t 1 t 0 1 2 l x=0 L y=0 ρ c I c ẇ2 ,x (x, y, t) + b c H c ẇ2 (x, y, t) -E c I f c w 2 ,xx (x, y, t) δ(y -y c ) dxdydt + t 1 t 0 1 2 l x=0 L y=0 ρ c I gc ẇ2 ,y (x, y, t) + G c I gc w 2 ,xy (x, y, t) δ(y -y c ) dxdydt = H bb + H bt (

3 6 ( 1 - 3 3

 613 vxyvyx) and D 4 = Gxyh the plate dynamic stiffness ; v xy and v yx = v xy Ey Ex Poisson's coefficients of the plate; E x and E y the two Young's modulus of the orthotropic plate ; E r the Young's modulus of a rib ; b and H respectively the width and the height of a rib ; I f and I g the rib momentum of bending and torsion inertia. Concerning the variables (.) c , they are referring to the straight rib oriented in the direction x.

  range[25;5000] Hz. Regarding R a||s , it could be R a () or R s () according to have B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 24

h= h 2 G xy z 2 +

 22 G yz x 2 dxdz = G yz Hb 3

42 Figure 1 :

 421 Figure 1: (a) top view of the stiffened rectangular wood panel ; (b) section view of the panel.

Figure 2 :

 2 Figure 2: Diagrams of displacement continuity at the interface plate / superstructure along

Figure 3 :

 3 Figure 3: Modes of the stiffened rectangular wood panel. (a) 1 st to 6 th modes (low frequen-

Figure 4 :

 4 Figure 4: Real parts of normalized radiation impedances of the 1 st (top: fig a) and 5 th

Figure 14 :

 14 Figure 14: Normalized radiation coefficient of the 271 th mode of a rectangular wood panel

6

  σ ( ω) as function of the frequency for an excitation on the strong beam (x=0.73 , y=0.Difference of radiation coefficient σ ( ω) as function of the frequency for an excitation on the strong beam (x=0.73 , y=0.61)

Figure 15 :

 15 Figure 15: Normalized radiation coefficient with and without taking into account crossmodal

Figure 16 :

 16 Figure 16: Drawing of the positions of the two excitation points.

Figure 19 :W

 19 Figure 19: Evolution of the acoustical radiated power L W = 10log(W/W ref ) [dB] as a

  It comes the vector of generalized effort Fgen and these components are defined by the following B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 18 ∂F m n a mn φ mn (x e , y e ) ∂a pq = F φ pq (x e , y e )

	relation:
	F pq =	∂E w ∂a pq	=

  xy (x r , y, t) 0 -2zw ,yy (x r , y, t) xw ,xy (x r , y, t)

				0
				-2zw ,yy (x r , y, t)
		⇐⇒ ¯ =	1 2	0 -zw ,xy (x
	(sym)	0

r , y, t) 0 xw ,xy (x r , y, t)

Table 2 :

 2 Acoustic radiation glossary

magnitude as the modal radiation coefficient of the mode coupled with itself (full plot) under 150 Hz and in the frequency range [300;500] Hz.

Nevertheless, these tendencies also exist for simply supported plate modes. Indeed, locally, the coupling terms of simply supported plate modes can be at the same order of magnitudes as coefficients with themselves, as show Fig. 4-b and4-d.

Indeed, for the mode (1,1) of a simply supported rectangular plate with the same dimensions, we also notice that differences between some couplings (dashed plots) and the radiation coefficient r 1,1;1,1 (full plot) are only around 5 dB under 100 Hz. See Fig. 4-b. For high order modes such as mode (5,5), presented in Fig. 4-d, we note that energy exchanges are sometimes higher than the exchange of the mode (5,5) with itself, from low to high frequencies. These less common results can be consulted in CVALOR [START_REF] Laulagnet | Validation commission of vibroacoustic calculation software packages: results of reference for the acoustic radiation of the B. Trévisan[END_REF] .

G. Radiation coefficient of ribbed modes

As recalled above, speaking about modal radiation coefficients supposes that modal interaction terms could be neglected initially. We make this assumption in the following, meaning we neglect the extra-diagonal terms of the matrix Rs (ω). This allows focusing on the notion of the radiation coefficient of a ribbed mode.

Plate ribbed by eleven stiffeners in one direction x

Although we can assume that a structure ribbed in both directions quasi homogeneously will have a similar behavior to the previous case, this is completely different for the ribbed panel. Indeed, the straight rib introduces a strong heterogeneity which splits the panel into two vibrating areas, as seen in section II.D.3. Thus we are now interested in the radiation coefficients of the ribbed wood panel (Fig. 1).

The presence of the straight rib along x tends to regularize the small variations of different heights / widths of each rib along y and implies a symmetrization effect on the two first modes. Therefore they are similar to modes (1,1) and (2,1) of a rectangular unribbed plate with the same dimensions 45;46 . Thus the conclusions are the same concerning the radiation coefficients, as shown in Fig. 9 and 10.

On the other hand, the heterogeneous aspect due to the straight rib appears from the 3th mode for which the panel is split into two. Thus the radiation coefficient is similar to that of a rectangular unribbed plate with reduced dimensions along y. Indeed, the radiation coefficient of the mode (3,1) of the reduced plate (dashed dark plot in Fig. 11) is a very good approximation of that of the 3 th mode (full plot). However, we note that at the first order the mode (3,1) of a plate with the same dimensions (light dash-dotted plot) is also a correct approximation, as shown by the small differences, with a maximum of 2 dB at 25 Hz.

Thus, excitation at a hard point of the structure (straight rib) tends to increase the radiation factor in this transitional frequency range: the apparent critical frequency decreases from 800 Hz (dashed plot) to 400 Hz (full plot). Indeed, far from the straight rib, the radiation factor is low under 800 Hz.

However, the conclusions are different for the two other indicators (Fig. 17-a and 17-b). Under 800 Hz the excitation point has a small influence. On the contrary, at a high frequency, we can observe that for both the acoustical radiated power and the space average quadratic velocity that the average values are around 15 dB lower for an excitation placed on a hard point. Obviously, the straight rib has a lower mobility which is why an excitation on the hardest point of this heterogeneous structure leads to decreasing the radiated power and the space average plate velocity.

Finally, these results show that an excitation on a hard point implies a decrease of the average velocity and of the acoustical radiated power in favor of better radiation efficiency at mid-frequency.

B. Perspectives: toward the piano soundboard

In this paper, we have mentioned the possibility to apply our model to the piano B. Trévisan, Kerem Ege and Bernard Laulagnet, JASA, p. 34 soundboard. The soundboard of a piano plays an essential role in the functioning of the instrument. Indeed, string sections are too thin to radiate on their own. Thus their vibrations are transmitted to the soundboard through the bridge which serves as an effective acoustic radiator. The vibroacoustic behavior of this ribbed structure has been studied by Suzuki, Conklin, Giordano, Boutillon With a work of mind, the bridge can be assimilated to the straight stiff rib of our semi-analytical model whereas the eleven smaller ribs can be assimilated to the stiffeners (See Fig. 18). Our results show that a smaller bridge could increase the acoustical power of high pitched notes but this also implies increased mobility and thus lower sustain of sound. This raises the question of the optimal dimensions of the bridge in order to find the best acoustical power/sustain trade-off.

An excitation placed on the straight rib / bridge provides another non negligible advantage for acoustical radiation of a stiffened structure: in other words the presence of this straight rib homogenizes the response of the ribbed structure. To demonstrate this, we consider an excitation moving along the straight rib and a second moving in the same direction at y e = 0.47m (see dashed lines on Fig. 16). Obviously, for the latter the excitation will sometimes occur on the plate, sometimes on the ribs. Fig. 19 presents the evolution of the acoustical radiated power as a function of the frequency along the x axis for these two cases.

Appendix I : Hamilton function calculus

In this Appendix, we give the main calculus steps for expressing the Hamilton function of the system without external force. The original expression is:

where E c refers to kinetic energy and E D the strain energy of the system. [t 0 , t 1 ] is an arbitrary time interval. We initially focus on the contribution of the plate and then the contributions of the superstructures.

Thus, we split the calculus into three sub-functions:

• H p for the plate

• H r for a rib in the direction y

• H b for the straight rib in the direction x with :

The contribution of the plate can be found in the literature 2 :


David Fernando

Novella Rodríguez
email: dnovellar@gmail.com

Del Muro Cuéllar

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Basilio Del Muro Cuéllar

Observer Based Scheme for the Control of High Order Systems with Two Unstable Poles Plus Time Delay

This paper deals with the problem of the stabilization and control of linear time invariant high order systems with two unstable real poles plus time delay. A simple observer based controller is designed in order to achieve a stable behavior of the closed loop system. Necessary and sufficient conditions for the existence of the proposed control structure are stated. Hence, only four proportional gains and the model of the plant are enough to obtain a stable response of the delayed system. Moreover, a robustness analysis is presented in order to compute the maximal uncertainty bound accepted for the delay term. In addition, a two degrees of freedom PI control action is implemented in order to track step references and to reject step disturbances. The achieved performance of the proposed control strategy is illustrated by mean of numerical simulations.

Introduction

Systems with significant delay term are very common, they appear in various systems as biological, ecological, social, industrial, etc., and are due to several mechanisms like material or energy transport, recycling loops, etc,. In addition, actuators, sensors and field networks that are involved in feedback loops usually introduce such delays, [START_REF] Richard | Time-Delay Systems:An Overview of Some Recent Advances and Open Problems[END_REF][START_REF] Zhong | Robust Control of Time-Delay Systems[END_REF]. Also delays can be used in model reduction where high-order (finite-dimensional) systems are approximated (in some norm sense) by low order systems with delays, [START_REF] Skogestad | Simple Analytic Rules for Model Reduction and PID Controller Tunning[END_REF]. It is known that time-delay is often a source of complex behaviors (oscillations, instability, bad performance), in many dynamic systems, and thus considerable attention has been paid on the stability analysis and controller design of time delay systems. Hence, there exists a great motivation to study delay effects on dynamical systems properties for two main reasons: first to understand how the delay presence may deteriorate the behavior of the system, and second to control their effects for better performance achievement on closed-loop systems, [START_REF] Gu | Stability of Time-Delay Systems[END_REF][START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF].

Several control strategies have been developed to deal with delayed systems. A common approach is to approximate the time-delay operator by means of Taylor or Padé series which could lead to a non minimum-phase system with rational transfer function representation. This method can be successfully used only when the delay term is small compared with the dominant time constant of the plant.

A different approach to deal with dead time systems is the classical Smith Predictor (SP), which consists in counteracting the time delay effects by mean of strategies intended to estimate the effects of current inputs over future outputs, [START_REF] Palmor | Time delay compensation smith predictor and its modifications[END_REF][START_REF] Smith | Close Control of Loops with Dead Time[END_REF]. The main limitation of the original SP is the fact that the prediction scheme does not have a stabilization step, which restricts its application to open-loop stable plants. However, open-loop unstable processes arise frequently in different dynamical systems and are fundamentally difficult to control. To overcome this problem, some modifications of the SP original structure have been proposed to deal with non-stable delayed process. For instance, [START_REF] Seshagiri | Simple analytical design of modified smith predictor with improved performance for unstable firstorder plus time delay (fodtp) processes[END_REF] has presented an efficient modification to the Smith predictor in order to control unstable first order system plus time delay, in this case the stability bound is determined by the relationship τ /τ un < 1.5, where τ is the process time delay and τ un the dominant unstable time-constant. With a different perspective, [START_REF] Normey-Rico | Dead-time compensators: A survey[END_REF] proposes a modification to the original Smith structure in order to deal with unstable first order delayed systems. Using a similar structure, the result is extended to delayed high order systems in [START_REF] Normey-Rico | Unified approach for robust dead-time compesator design[END_REF]. In both works, a robustness analysis is done concluding that for unstable dead time dominant systems, the closed-loop system can be unstabilized with an infinitesimal value of the modeling error, i.e., that robustness is strongly dependent on the relationship τ /τ un . For the control scheme proposed in this later work, it can be easily proven that in the case of unstable plants, the internal stability is not guaranteed. In fact it is obtained an unstable estimation error and, as a result, a minimal initial condition difference between the original plant and the model produces an internal unbounded signal.

The classical controllers Proportional-Integral-Derivative (P, PI, PID or PD) are also included between the tools used to deal with time delay systems. For instance, [START_REF] Visioli | Control of Integral Processes with Dead Time[END_REF] introduce a study about the control of integral systems with dead time which includes the design and tune of PI and PID controllers by mean of different methods. In [START_REF] Silva | PID Controllers for Time-Delay Systems[END_REF], the authors proposes a methodology to obtain a complete PI and PID parametrization for first order unstable delayed systems when τ < 2τ un . With a different perspective, in [START_REF] Middleton | On the Achievable Delay Margin Using LTI Control for Unstable Plants[END_REF] upper bounds on the delay size (τ < 2τ un ) are provided when using linear time invariant controllers on the stabilization of strictly proper delayed real rational plants. It is important to note that in general, the provided bounds are not tight and the authors prevent that the developed controllers are not intended as practical solutions and are only used as a tool to compute the achievable delay margin for some particular cases. In [START_REF] Michiels | Continuous pole placement for delay equations[END_REF], based in a numerical method, it is considered the stabilization of linear time-delay systems of order n. However, stability conditions with respect to time-delay and time constant of the process are not provided. The proposed method consists in shifting the unstable eigenvalues to the left half plane by static state feedback by applying small changes to the feedback gain. Using an observer based approach, in [START_REF] Del Muro Cuéllar | On the control of unstable first order linear systems with large time lag: Observer based approach[END_REF] a methodology is proposed in order to control unstable first order delayed systems allowing to deal with τ /τ un < 3 and with a similar approach but including PID controllers in the loop, in [START_REF] Márquez-Rubio | Observer PID stabilization strategy for unstable first-order linear systems with large time delay[END_REF] is stated a methodology to stabilize FOTD unstable Systems with τ /τ un < 4.

The problem of proposing a controller scheme to deal with high order unstable systems, providing specific conditions on the delay size is more complex and they are not many works on this problem. For instance, in the work presented by Madhuranthakam et al., [START_REF] Madhuranthakam | Optimal tuning of PID controllers for foptd, soptd and soptd with lead processes[END_REF] it is proposed a method to the optimal tuning of PID controllers for first and second order systems plus time delay, nevertheless the method is restricted to stable processes. With a different perspective, [18] presents a strategy based on a modified smith predictor in order to control second order system but only the case of one unstable pole is dealt. On the other hand, [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF] introduces a PID control strategy for second order unstable processes designed using the direct synthesis method. With a similar perspective, [START_REF] Shamsuzzoha | Design of IMC filter for PID control strategy of open-loop unstable processes with time delay[END_REF] presents a methodology in order to control second order delayed systems using the Internal Model Control guides for PID controllers. However, this works do not provide the stability bound in terms of the time delay size.

This paper is concerned with the stabilization problem of systems with two unstable poles and m stable poles plus time delay. The control scheme relies on an observer-based structure with a memory observer and a memoryless state feedback, so only two gains are enough to stabilize the observer scheme and two other to stabilize the open loop unstable plant. Necessary and sufficient conditions are stated to guarantee the existence of the proposed scheme in terms of the relation between the size of the delayed term and the system time constants. On the contrary of modified Smith predictors, the scheme only contains discrete time delay (and not distributed ones) which makes easy its practical implementation (see [START_REF] Zhong | Robust Control of Time-Delay Systems[END_REF] for details on numerical implementation of modified Smith predictor scheme). Moreover, it is implemented a Proportional-Integral controller which allows to track step references. It is worth stressing that, to the best of our knowledge, in the literature there is not reported similar results for this class of systems.

This paper is organized as follows; the Section 2 is dedicated to the problem formulation. The Section 3 yields the preliminaries results used to obtain the main result of this work. An observer based controller is proposed in the Section 4 in order to stabilize the unstable delayed system previously described, also the stability conditions of the proposed control structure are stated. In addition, it is illustrated how to deal with a zero in the left half plane in the delayed transfer function. Then, in Section 5 a robustness analysis for the observer based scheme with respect to the uncertainty in the delay operator is presented. Numerical simulations are presented in order to show the controller performance under different conditions in Section 6. Moreover, an application example related with a continuous stirred tank reactor (CSTR) is presented in order to refer the control strategy proposed with a practical framework. Finally some conclusions are drawn in Section 7.

Problem Formulation

Consider the following class of single-input single-output (SISO) linear systems with input delay:

Y (s) U (s) = N (s) D(s) e -τ s = G(s)e -τ s , (1)
where U (s) and Y (s) are the input and output signals respectively, τ ≥ 0 is the constant time delay, N (s) and D(s) are polynomials in the complex variable s and G(s) is the delay-free transfer function. Notice that with respect to the class of systems (1) a traditional control strategy based on an output feedback of the form:

U (s) = C(s)[R(s) -Y (s)], (2)
yields a closed-loop system given by:

Y (s) R(s) = C(s)G(s)e -τ s 1 + C(s)G(s)e -τ s , (3)
where the exponential term e -τ s located at the denominator of the transfer function (3) leads to a system with an infinite number of poles and where the closed-loop stability properties must be carefully stated. This work proposes an observer based control scheme in order to stabilize a system characterized by the following transfer function:

Y (s) U (s) = α (s -a)(s -b)(s + c 1 )(s + c 2 )...(s + c m ) e -τ s . (4)
Where τ ≥ 0, a, , b , c 1 , c 2 ...c m > 0 are positive constants, and without loss of generality, a ≥ b > 0. The proposed control scheme has been designed taking into account the traditional observer theory hence only the plant model and two static gains are enough to get an adequate estimation of an internal delay free variable which will be used in the final stabilizing control scheme.

For sake of simplicity, we denote the stable subsystem G stb (s) as follows

G stb (s) = α (s + c 1 )(s + c 2 )...(s + c m ) . ( 5
)

Preliminary Results

Preliminary results are presented, which will be used later in order to state the stability conditions of the proposed strategy in this work.

Considering the following unstable first order system plus time delay:

Y (s) U (s) = G(s)e -τ s = α s -σ e -τ s , (6)
with σ > 0, and a proportional output feedback control as follows:

U (s) = R(s) -kY (s), (7)
which produces a closed-loop system:

Y (s) R(s) = αe -τ s s -σ + kαe -τ s . ( 8
)
The following result has been widely studied in the literature and the proof can be easily obtained by considering different approaches as a classical frequency domain. An alternative simple proof based on a discrete time approach is shown in [START_REF] Del-Muro-Cuéllar | Stabilisation and control of unstable first-order linear delay systems[END_REF].

Lemma 1. Taking into account the delayed system (6) and the proportional output feedback [START_REF] Smith | Close Control of Loops with Dead Time[END_REF]. Then, there exists a proportional gain k such that the closed loop system [START_REF] Seshagiri | Simple analytical design of modified smith predictor with improved performance for unstable firstorder plus time delay (fodtp) processes[END_REF], is stable if and only if τ < 1 σ . Now, take into consideration the high-order unstable system characterized by:

G(s) = α (s -σ)(s + φ 1 )(s + φ 2 )...(s + φ n ) e -τ s . (9)
with σ, φ 1 , φ 2 ...φ n > 0. With the proportional output feedback [START_REF] Smith | Close Control of Loops with Dead Time[END_REF], the closed-loop system obtained is:

Y (s) R(s) = αe -τ s [(s -σ)(s + φ 1 )(s + φ 2 )...(s + φ n )] + kαe -τ s . ( 10
)
Lemma 2. Considering the delayed system (9) and the proportional output feedback [START_REF] Smith | Close Control of Loops with Dead Time[END_REF]. Then, there exists a proportional gain k such that the closed loop system [START_REF] Normey-Rico | Unified approach for robust dead-time compesator design[END_REF], is stable if and only if τ < 1 σ -

n i=1 1 φ i .
The proof of the previous lemma is presented in the Appendix A.

Control Strategy Proposed

Taking into account the class of systems studied in this work and characterized by the transfer function ( 4) with a, b, c 1 , c 2 , • • • , c m > 0 and the time delay τ ≥ 0, and assuming without loss of generality a ≥ b. An observer based control is designed in order to obtain an estimation of the internal states of the system to be used as control signals for the original process.

As a first step, the stability conditions for the controller and the observer systems are stated separately. This conditions will be used later in order to state the closed loop stability conditions for the proposed observer based controller.

Controller Scheme

First, taking into consideration the controller structure shown in Figure 1, with the control law u(t) = r(t) -k 1 w(t) -k 2 y(t), let us introduce the following result. 4), and the control law mentioned above. There exist gains k 1 and k 2 such that the closed-loop system is stable if and only if

τ < 1 b - m i=1 1 c i .
The aim of the foregoing proof is to apply the stability conditions given in Lemma 2 to the state feedback strategy shown in the Figure 1.

Proof. Sufficiency. Let us consider τ < 1 b - m i=1 1 c i . Then, τ = 1 b - m i=1 1 c i -β,
for some β > 0. Therefore, there exists

k 1 such that β > 1 k1-a > 0. Then, it is easy to determine τ < 1 b - m i=1 1 c i - 1 k 1 -a . (11)
Finally, we can conclude from Lemma 2, where σ = b, φ i = c i and φ m+1 = k 1 -a there exists k 2 such that the plant behavior is stable.

Necessity. Considering the delayed system (4), and the state feedback controller shown in Figure 3, with constant gains k 1 and k 2 such that the process is stable. The closed loop transfer function of the plant can be written as follows:

Y (s) R(s) = αe -τ s (s -b)(s + c 1 ) • • • (s + c m )(s + φ) + αk 2 e -τ s , (12)
with φ = k 1 -a. It is well known that a k 2 that stabilizes the delayed system (12) must also stabilize the delay free system (see for instance [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF] or [START_REF] Malakhovski | On Stability of Second-order Quasi-polynomials with a Single Delay[END_REF]), which

implies that φ > 0. Indeed, from Lemma 2, τ < 1 σ - n i=1 1 φ i with σ = b and φ i = c i , where n = m + 1 (note that φ > 0 is a free parameter function of k 1 ). Let us consider β > 1 φ > 0, denoting β = 1 σ - m i=1 1 φ i -τ , therefore: τ = 1 σ - m i=1 1 φ i -β < 1 σ - m i=1 1 φ i = 1 b - m i=1 1 c i Remark 1.
A simple methodology for choosing the values for the gains of the controller will be given below.

Step 1. In order to ensure the existence of a proportional gain k 2 such that the closed loop system is stable, from the proof of Lemma 3, eq.( 11) we obtain

k 1 > 1 1 b - m i=1 1 c i -τ + a. (13)
It is worth stress that if k 1 is selected near to the bound stated in [START_REF] Middleton | On the Achievable Delay Margin Using LTI Control for Unstable Plants[END_REF], the gain margin for the gain k 2 will be reduced.

Step 2. Once a gain k 1 is selected, we can compute the gain k 2 by means of a frequency domain analysis, Nyquist stability criterion for instance, for the auxiliary system [START_REF] Silva | PID Controllers for Time-Delay Systems[END_REF] such that the controller scheme shown in the Figure 1 is stable.

Observer Scheme

In most of the practical applications, the internal variables are not measured. Thus, an observer based on an output injection strategy is proposed, let us take into consideration the static output injection scheme shown in Figure 2, the stability of the observer can be tackled as follows.

Lemma 4. Considering the delayed system (4), and the static output injection scheme shown in Figure 2. There exist constants g 1 and g 2 such that the closedloop system is stable if and only if

τ < 1 a - m i=1 1 c i .
Proof. The proof can be easily derived from a dual procedure of the previous result. Remark 2. As in the controller design, the computation of the gains g 1 and g 2 can be obtained in a dual way to Remark 1, i.e.,

Step 1. In order to ensure the existence of a proportional gain g 1 such that the closed loop system is stable, from Lemma 4 we obtain

g 2 > 1 ( 1 a - m i=1 1 c i -τ ) + b.
Step 2. Once the gain g 2 is selected, we can compute the gain g 1 by means of a frequency domain analysis, Nyquist stability criterion for instance, such that the controller scheme shown in Figure 5 is stable.

Observer-Based Controller: the case without zeros

Finally, the main result of this work is presented, we propose an observed based controller as in the Figure 3, where the observer allows to estimate the state variables, to be used in state feedback controller. The authors would like to stress that, in the proposed scheme, only four proportional gains are enough to get a stable closed loop behavior. As a consequence of the previous results, the following lemma can be stated. Lemma 5. Consider the observer based controller scheme shown in Figure 3. There exist proportional gains k 1 , k 2 , g 1 and g 2 such that the closed-loop system is stable if and only if

τ < 1 a - m i=1 1 c i .
Proof. Consider a state space representation of the system (4) characterized by the following equation:

ẋ(t) = A 0 x(t) + A 1 x(t -τ ) + Bu(t) (14)
y(t) = Cx(t), with x(t) = w(t) x 1 (t) x 2 (t) • • • x m (t) z(t) T ,
where the states x i (t) represent the stable part of the system, w(t) and z(t) are linked with the unstable poles of the open loop system. Then we have,

A 0 =          a 0 0 • • • 0 0 0 1 -c 1 0 • • • 0 0 0 0 1 -c 2 • • • 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 • • • 1 -c m 0 0 0 0 • • • 0 0 b          , A 1 =          0 0 0 • • • 0 0 0 0 0 0 • • • 0 0 0 0 0 0 • • • 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 • • • 0 0 0 0 0 0 • • • 0 1 0          , B =          α 0 0 . . . 0 0          , C = 0 0 0 • • • 0 1 .
Note that the state state representation characterized by ( 14) can be returned to its transfer function representation by mean of:

Y (s) U (s) = C(sI -(A 0 + A 1 e -τ s )) -1 B, (15)
which brings us back to the delayed transfer function [START_REF] Gu | Stability of Time-Delay Systems[END_REF]. The dynamics of the estimated states and the control law can be described as follows.

ẋ(t) = A 0 x(t) + A 1 x(t -τ ) + Bu(t) -G(C x(t) -y(t)) (16)
Where x(t) is the estimated state of x(t), and the gain vectors K and G are defined by

K = k 1 0 0 • • • 0 k 2 , G = g 1 0 0 • • • 0 g 2 T .
Let e(t) := x(t) -x(t), then we have:

ė(t) = ẋ(t) -ẋ(t) = (A 0 -GC)e(t) + A 1 e(t -τ ), (17)
and the controlled system:

ẋ(t) = A 0 x(t) + A 1 x(t -τ ) -BK x(t). (18)
Noting x e = [x(t) e(t)] T and after a simple manipulation of variables we have the following closed loop system with the observer and the controller proposed in the Figure 3 ẋe

(t) = A 0 -BK BK 0 A 0 -GC x e (t) + A 1 0 0 A 1 x e (t -τ ) (19)
y(t) = C 0 x e (t).
It is easy to see that the observer based controller proposed satisfies the separation principle. Hence, the stability of the observer scheme is enough to assure an adequate error convergence, i.e. there exist proportional gains g 1 and g 2 such that lim

t→∞ [ ŵ(t) -w(t)] = 0 if and only if τ < 1 a - m i=1 1 c i ,
then, considering the fact of the observer and controller can be designed separately and reminding the stability conditions stated previously in Lemmas 3 and 4, is clear that the observer stability condition is more restrictive than the controller one, i.e., 1 a -

m i=1 c i < 1 b - m i=1 1 c i ,
therefore, there exist k 1 , k 2 , g 1 and g 2 such that the closed-loop system is stable if and only if

τ < 1 a - m i=1 1 c i .
Remark 3. Different works have been presented recently in order to obtain a maximum upper bound for the time delay. For instance, in [START_REF] Middleton | On the Achievable Delay Margin Using LTI Control for Unstable Plants[END_REF] the authors adopt a frequency domain approach to demonstrate that, for a strictly proper rational plant, there exists an uniform upper bound on the delay that can be tolerated when using a finite dimensional linear time invariant (FDLTI) controller. Thus, for plants with an unstable pole p, the maximum margin delay obtained is DM ≤ 2/p. Furthermore, the authors consider a controller (see [START_REF] Middleton | On the Achievable Delay Margin Using LTI Control for Unstable Plants[END_REF], Remark 8) to achieve this delay margin, although this class of controller may suffer difficulties such poor sensitivity, phase margin, etc., therefore, it is not intended as a practical controller. Also, the authors have the assumption that, for a plant with two unstable poles (the class of systems dealt in this work) the delay margin of the system is DM ≤ min{2/p 1 , 2/p 2 }, where p 1 and p 2 are the unstable poles of the system. However the authors do not propose a controller in order to obtain this delay margin.

Observer-Based Controller: the case with a left half plane zero

The proposed control strategy can be easily implemented when the system contains a left half plane zero. For this, the main idea is to cancel the zero dynamics for each case (controller and observer), placing an unstable pole in the position of the zero. For this consider the stable subsystem G stb (s) given by the following transfer function

G stb (s) = α(s + β) (s + c 1 )(s + c 2 )...(s + c m ) . ( 20
)
Remark 4. As was mentioned before, the main idea to deal with the case of systems with a left half plane zero is to cancel the zero dynamics by means the relocation of an unstable pole, from this we can state the following procedure to obtain the gains for the observer based structure in this case.

Step 1. In order to cancel the zero dynamics by means of the selection of k 1 in the controller scheme, we have the following equality

k 1 = β + a (21)
Step 2. From this, by means of a frequency domain analysis, Nyquist stability criterion for instance, we can compute the gain k 2 for the obtained system

Y (s) R(s) = αe -τ s (s -b)(s + c 1 ) • • • (s + c m ) , (22)
Step 3. In order to cancel the zero dynamics by means of the selection of g 2 in the observer scheme, we have the following equality

g 2 = β + b (23)
Step 4. From this, by means of a frequency domain analysis, Nyquist stability criterion for instance, we can compute the gain g 1 for the obtained system

Y (s) R(s) = αe -τ s (s -a)(s + c 1 ) • • • (s + c m ) , (24)

Step tracking reference

In order to solve the step tracking reference problem, the observed based structure can be complemented with a proportional-integral (PI) control action, shown in Figure 4. For high-order time delay systems, a PI controller design is proposed in [START_REF] Lee | Stabilization of all-pole unstable delay processes by simple controllers[END_REF]. As it is well known, conventional methods of PI controllers design involve a zero on the closed loop system which yields an undesirable overshoot on the response. Therefore, a two degrees of freedom PI controller will be considered to overcome this undesirable overshot response, [START_REF] Visioli | Control of Integral Processes with Dead Time[END_REF], [START_REF] Astrom | PID Controllers: Theory, Design, and Tuning[END_REF]. The PI controller is characterized as follows

U (s) = G f f (s)R(s) -G c (s) Ŷ (s), (25)
where

G c (s) = K p 1 + K I s (26)
and

G f f (s) = K p ǫ + K I s . (27)
Here, the main idea is to relocate the unstable pole a by means of the gain k 1 , then we have a open loop process with an unstable pole located at s = b and m + 1 stable poles, characterized by the equation [START_REF] Silva | PID Controllers for Time-Delay Systems[END_REF]. Then it is possible to use the methodology given in [START_REF] Lee | Stabilization of all-pole unstable delay processes by simple controllers[END_REF] in order to tune a PI controller. We can resume the procedure as follows (for more details see [START_REF] Lee | Stabilization of all-pole unstable delay processes by simple controllers[END_REF]); first the gain K I must be selected such that Φ(ω) > -π, for some w ⋆ > 0, with the function

Φ(ω) = -τ ω -π + arctan ω b -arctan K I ω - m+1 i=1 arctan ω c i . (28)
Then, the range of stabilizing gains K p is given by

(ω 2 c1 + b 2 ) m i=1 +1(ω 2 c1 + c 2 i ) 1 + KI ωc1 < K p < (ω 2 c2 + b 2 ) m i=1 +1(ω 2 c2 + c 2 i ) 1 + KI ωc2 ,
with ω c1 < ω c2 being the first two phase frequencies solved from

-τ ω + arctan ω b -arctan K I ω - m+1 i=1 arctan ω c i = 0.
Then, in order to reduce the undesirable overshoot in the response, the parameter ǫ can be chosen from 0 < ǫ < 1.

Robustness Analysis With Respect To Delay Uncertainty.

Here, we refer to robustness with respect to delay to the fact that the exact value of the delay is not known a priori. Therefore, we assume now that the delay is uncertain, i. e. the delay of the real system is τ = d(1 + δθ) and may be different from the one used in the observer (the nominal delay d). Consider the delay uncertainty in the original process, and the design of the observer taking into account the nominal time delay, i.e.

ẋ(t) = A 0 x(t) + A 1 x(t -τ ) + Bu(t) (29) ẋ(t) = Ax(t) + A 1 x(t -d) + Bu(t) -G(C x(t) -y(t)) u(t) = -K x(t)
Let us consider the observer-controller system with the form (29) and its extended closed loop system considering the delay uncertainty characterized as follows.

ẋe (t) =

A 0 -BK BK 0 A 0 -GC x e (t) + A 1 0 A 1 0 (x e ) τ (t -τ ) + • • • • • • + 0 0 -A 1 A 1 (x e ) d (t -d). ( 30
) = A 0 x e (t) + A τ (x e ) τ (t -τ ) + A d (x e ) d (t -d) (31)
Note that when the delay value of the actual process matches to the nominal one, the extended closed loop system (31) is equivalent to the nominal closed loop system [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF], [START_REF] Verriest | Robust Observer-Controller for Delay Differential Systems[END_REF]. Now modeling the time delay with a multiplicative uncertainty as τ = d(1 + δθ), with θ normalized, i.e. |θ| ≤ 1 we obtain.

e -sτ = e -sd(1+δθ) = e -sd e -sdδθ = e -sd (1 -∆) (32) With ∆(s) = 1 -e -sdδθ . Therefore the characteristic polynomial of the above system is given by:

Ψ(s) = det[Ψ 0 (s)]det[I n + Ψ -1 0 (s)A τ e -sd ∆(s)] ( 33
)
Where Ψ 0 (s) = I 2n -A 0 -(A τ + A d )e -sd . Now, when the design of the controller and the observer ensures that the nominal extended system is stable, i. e. det[Ψ 0 (s)] is stable, the perturbed closed loop system remains stable if det[I n + Ψ -1 0 (s)A τ e -sd ∆(s)] does not change sign when s sweeps the imaginary axis. Invoking Rouche's Theorem, it follows that the condition for stability is:

Q d (s)∆(s) ∞ < 1 ( 34
)
Where Q d (s) = Ψ -1 0 (s)A τ e -sd . Note that the term 1-e -sdδθ = sdδθe -sdδθ/2 | s=jω sinc ωdδθ 2π . Hence, on the imaginary axis, |1 -e -sdδθ | ≤ |sdδθ|, then if one can show stability for 1 -e sdδθ replacing it by sdδθ in the analysis, then the stability for the uncertain system will follow. Note that this amounts to replacing the delayed term by its first order Taylor expansion. This means that the maximal uncertainty bound δ that preserves stability for the worst case, i.e. θ = 1, is determined by:

δ max = 1 sde -sd Ψ -1 0 (s)A τ ∞ (35)
Then for all |θ| ≤ 1 and the uncertainty bound δ max the determinant has a fixed sign, implying the absence of zero crossings, and henceforth the stability of the perturbed system (provided the nominal one is stable) [START_REF] Verriest | Robust Observer-Controller for Delay Differential Systems[END_REF][START_REF] Sename | Is a Mixed Design of Observer Controllers for Time Delay Systems Interesting?[END_REF].

Remark 5. Note that the stabilizing gains computed for the observer based controller for the nominal time delay d are able to preserve stability of the closed loop system with a nominal delay belongs to the interval [0 d], the analysis of the stability sets and regions are dealt in [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF] and [START_REF] Malakhovski | On Stability of Second-order Quasi-polynomials with a Single Delay[END_REF]

Examples

Below, numeric simulations are performed in order to illustrate the observer based scheme performance. Moreover, the PI control action is implemented to achieve step reference tracking and step disturbance rejection.

Example 1. Consider the fourth order delayed system with two unstable poles characterized with the following transfer function:

Y (s) U (s) = 5 (s -0.7)(s -0.4)(s + 5)(s + 10) e -s (36)
Let be a = 0.7, b = 0.4, c 1 = 5, c 2 = 10, α = 1 and τ = 1 the parameters of the system, it is clear that the stability condition given in Lemma 5 is satisfied, therefore there exists an observer based structure with proportional gains k 1 , k 2 , g 1 and g 2 such that the resulting closed-loop system is stable due to:

τ = 1 < 1 a - m i=1 1 c i = 1.128.
Following the methodology given in remarks 1 and 2, the stability regions for the stabilizing gains are presented in Figure 5, from this the selected values for the control structure are shown in Table 1. Figure 6 illustrate the output performance of the observer based controller in numerical simulations. The continuous line indicates the output of the closed loop system with identical initial conditions between y(t) and ŷ(t). The dashed line point to the system performance whit different initial conditions, (ŷ(0) -y(0) = 0.01). As can be shown in Figure 5, the selection of a larger value of k 1 , (resp. g 2 ), gives a larger gain margin for k 2 , (resp. g 1 ). However, a selection of high gains for the observer based controller could give the problem of high control signals.

k 1 10.7 k 2 55 g 1 145 g 2 20.4
Remark 7. It is possible to implement a similar controller as the presented in [START_REF] Middleton | On the Achievable Delay Margin Using LTI Control for Unstable Plants[END_REF] (Remark 8) in order to obtain better stability conditions for the control strategy proposed in this work, however some properties of the system (phase margin, gain margin, etc.) can be affected.

Robustness Analysis

The Figure 7 shown the output signal of the closed loop system with different initial conditions between the plant and the observer system,(ŷ(t)-y(t) = 0.01), when there is not uncertainty in the delay value (continuous line), and the closed loop system behavior when the maximum uncertainty bound is presented in the system (dashed line). The following table shows the relation between the size of the time delay d of the nominal system and the maximal uncertainty bound δ max in the above example. The terms τ min and τ max indicate the minimal and maximal value respectively of the delay term τ , such that the closed loop system with nominal delay d remains stable. Example 2. Chemical reactors often have significant heat effects, so it is important to be able to add or remove heat from them. In a jacketed CSTR (continuously stirred tank reactor) the heat is added or removed by virtue of the temperature difference between a jacket fluid and the reactor fluid. Often, the heat transfer fluid is pumped through agitation nozzles that circulate the fluid through the jacket at high velocity, as is shown in Figure 8. Here we consider a CSTR carrying out the simple reaction A → B. The balance on component A is

V dC A dt = F C Af -F C A -V r A , (37)
where C A is the concentration of component A in the reactor and r A is the rate of reaction per unit volume. The Arrhenius expression is normally used for the rate of reaction. A first order reaction results in the following

r A = k 0 e ( -∆E RT CA) , (38)
where k 0 is the frequency factor, E a is the activation energy, R is the ideal gas constant, and T is the reactor temperature on an absolute scale. A detailed model of a CSTR includes the effect of cooling jacket dynamics. Then, the three modelling equations are

dC A dt = F V (C Af -C A ) -k 0 e ( -∆E RT ) C A , (39)
dT dt = F V (T f -T ) + -∆H RT k 0 e ( -∆E RT ) C A - U A V ρc p (T -T j ), ( 40
)
dT j dt = F jf V j (T jf -T j ) + U A V j ρ j c pj (T -T j ). ( 41
)
Where where -∆H is the heat of the reaction, U is the heat transfer coefficient, A is the heat transfer area, T f is the feed temperature, T j is the jacket temperature, T jf is the jacket feed temperature and F jf the jacket make-up flow rate, [START_REF] Bequette | Process Control. Modeling, Design and Simulation[END_REF].

The parameter values of the system are given as ∆E = 20000Btu/lbmol, k 0 = 16.96 × 10 7 hr -1 , -∆H = 19000Btu/lbmol, U = 75Btu/hr • F f t 2 , ρc p = 53.25Btu/ • F f t 3 , ρ j c pj = 55.6Btu/ • F f t 3 . The operating values of the reactor are given by A = 88f t 2 , C Af = 0.132lbmol/f t 3 , T f = 60 • F , the operating volume V = 85f t 3 , the flow rate F = 35f t 3 /hr, the jacket volume V j = 25f t 3 , the jacket temperature T jf = 0 • F and the jacket make-up flow rate is F jf = 28.75f t 3 /hr. A steady state operating point is C A = 0.166lbmol/f t 3 and T = 101.1 • F . Let consider the flow rate as the manipulated variable and the temperature of the CSTR as the controlled variable. Linearisation around this steady state operating point yields the following transfer function model (by assuming a measurement time delay of 54 minutes) T (s) F j (s) = -4.7475(s + 3.129) (s -0.5545)(s -0.09395)(s + 6.394) e -0.9s (42)

Figure 9 illustrates the output response of the closed loop system stabilized by means of the observer-based controller for a unit step reference in numerical simulations, the computed gains are shown in Table 3. Continuous line denotes the nominal system, dashed line denotes the closed loop system with an uncertainty of 10% in the real process.

For the present example, following the methodology proposed in [START_REF] Lee | Stabilization of all-pole unstable delay processes by simple controllers[END_REF], the parameters computed for the PI controller are shown in Table 3 and for the feed forward compensator G f f , the set-point weighting parameter ǫ = 0.01. The proportional gains k 1 , g 1 and g 2 are the same that the computed before. The closed loop response is shown in Figure 10 applied at t = 50hr, as it can be noticed, also the PI controller is able to reject step disturbances. The dashed line point to the behavior for the closed loop system with different initial conditions in the observer scheme (y(t) -ŷ(t) = 0.2). Remark 8. Numerical simulations above show an efficient performance of the control strategy proposed even in face of different initial conditions as well as parametric uncertainties in the system. Example 3. The following example is presented in [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF], it is a second order system plus time delay (a particular case of the class of systems dealt in this work) with the transfer function given by

k 1 3.6835 k 2 -1 g 1 -1.2 g 2 3.2230 K p -0.6 K I 0.2774
Y (s) U (s) = 2 (3s -1)(s -1) e -τ s (43)
For this example, the parameters of the system are a = 1, b = 0.333, α = 0.6667 and τ = 0.3.The gains computed for the observer and the PI controller 4. Figure 11 shows the response of the system of example 3. A step negative disturbance is presented in t = 30. Solid line point the performance for the control strategy proposed in this work and dashed line shows the performance of the methodology proposed in [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF]. It can be noticed that the observer based structure proposed in this work gives a better performance that the PID controller proposed in [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF]. In addition, the methodology proposed in [START_REF] Rao | Control of unstable processes with two rhp poles, a zero and time delay[END_REF] does not state bounds in terms of the time delay size, while our control structure can be implemented in the system with a larger delay as it is shown in Figure 12. The system of Example 3 is controlled with the observer based controller proposed, but with a time delay τ = 0.8. The computed parameters of the observer structure and the PI controller are shown in Table 4.

τ = 0.3 τ = 0.

Conclusions

An observer based controller is proposed in order to stabilize high order system with two unstable poles plus time delay. The necessary and sufficient conditions that ensure the existence of the stabilizing controller scheme are stated. The scheme is simple and may be easily implemented: only four proportional gains and the model of the process are enough to stabilize the system. The procedure to design the controller can be performed easily using well know analysis of classical control theory. Numerical examples illustrate the performance of the controller under ideal conditions. The examples show the closed loop system behavior working under different initial conditions between the plant and the observer, as well as parametric variations. Moreover, the introduction of a PI controller allows to one to solve the problem of regulation and step disturbance rejection. Also, a robustness analysis is presented to show the effects of the modeling errors of the delay operator.

First, analyzing the First Order Unstable Delayed System given by: G(s) = α s -σ e -τ s . (A.2)

Considering the Lemma 1, there exists a proportional gain k such that a closed loop system with a simple proportional output feedback is stable if and only if τ < 1

σ . An analysis in the frequency domain shall confirm this result. The Nyquist stability criteria establishes that, when closing the loop whit a proportional gain k, the system will be stable if 0 = N + P , with P being the number of poles in the right half plane "s" and N the numbers of clockwise round trips to the point -1 (N negative suggests round trips in the opposite direction) in the Nyquist diagram. The angle as a function of the frequency ω is given by:

∠G(jω) = -(180 • -arctan ω σ ) -(ωτ ) (A.3)
Taking into consideration that for small frequencies arctan ωϕ ≈ ωϕ, it can be shown that the condition τ < 1 σ is equivalent to ask that the angle path taps at least one point (for some frequency) with a value exceeding -180 • , that is ∠G(jω) > -180 • , i.e., one counter-clockwise round trip to the point -1 in the Nyquist diagram. Now, analyzing the system with n = 1 characterized by: G(s) = α (s -σ)(s + φ) e -τ s (A.4)

Then, there exists a proportional gain k such that a closed loop system whit a simple proportional output feedback is stable if and only if τ < 1 σ -1 φ . It can be easy to see that the Nyquist condition remains the same (one counterclockwise round trip to the point -1 in the Nyquist diagram), now the angle condition is:

∠G(jω) = -(180 • -arctan ω σ ) -(arctan ω φ ) -(ωτ ) (A.5)
For small frequencies arctan ωϕ ≈ ωϕ, and starting from ∠G(jω) > -180 • it is no difficult to conclude the relation τ < 1 σ -1 φ . Below considering the system whit n = 2 given by:

G(s) = α (s -σ)(s + φ 1 )(s + φ 2 ) e -τ s . (A.6)
It can be easy to see that the Nyquist condition remains the same (one counter-clockwise round trip to the point -1 in the Nyquist diagram), now the angle condition is:

∠G(jω) = -(180 • -arctan ω σ ) -(arctan ω φ 1 ) -• • • • • • -(arctan ω φ 2 ) -(ωτ ) (A.7)
Again considering that for small frequencies arctan ωϕ ≈ ωϕ, and starting from ∠G(jω) > -180 • it is no difficult to conclude the relation τ < 1 σ -1 φ1 -1 φ2 . This reflection can be generalized to any n ∈ R concluding that for the systems characterized by (A.1) with a proportional output feedback, there exist a constant gain k such that the closed loop system is stable if and only if:

τ < 1 σ - n i=1 1 φ i . (A.8)

Figure 1 :

 1 Figure 1: Control Scheme

Figure 2 :

 2 Figure 2: Observer Scheme.

Figure 3 :

 3 Figure 3: Control Strategy Proposed.

Figure 4 :

 4 Figure 4: Control Structure

1 Figure 5 :

 15 Figure 5: Stability regions for proportional gains in the observer based controller

Figure 6 :

 6 Figure 6: Numerical Simulation Results.

Figure 7 :

 7 Figure 7: Robustness Analysis with respect to Time Delay

Figure 8 :

 8 Figure 8: Continuos Stirred Tank Reactor

  Figure9illustrates the output response of the closed loop system stabilized by means of the observer-based controller for a unit step reference in numerical simulations, the computed gains are shown in Table3. Continuous line denotes the nominal system, dashed line denotes the closed loop system with an uncertainty of 10% in the real process.For the present example, following the methodology proposed in[START_REF] Lee | Stabilization of all-pole unstable delay processes by simple controllers[END_REF], the parameters computed for the PI controller are shown in Table3and for the feed forward compensator G f f , the set-point weighting parameter ǫ = 0.01. The proportional gains k 1 , g 1 and g 2 are the same that the computed before. The closed loop response is shown in Figure10, continuous line shows the performance of the observer based controller with the complement of the PI control action considering the nominal plant. A step disturbance with magnitude d = 0.01 is

Figure 9 :

 9 Figure 9: Control Strategy Performance

20

 20

Figure 10 :

 10 Figure 10: Step Tracking Performance

Figure 11

 11 Figure 11: Performance Comparison

Figure 12 :

 12 Figure 12: Closed Loop Performance

Table 1 :

 1 Observer Based Controller Gains: Example 1.

Table 2 :

 2 Robustness Analysis with different time delay size.

Table 3 :

 3 Observer Based Controller Gains: Example 2.

Table 4 :

 4 Observer Based Controller Gains: Example 3.

Appendix A. Proof of The Lemma 2.

Let us consider the High-Order Unstable System with time delay given by: