Ehsan Enferad

Christophe Giraud-Audine
email: christophe.giraud-audine@ensam.eu

Frédéric Giraud

Michel Amberg

Betty Lemaire Semail

Generating controlled localized stimulations on haptic displays by modal superimposition

Keywords: Modal Analysis, Modal synthesis, Human Machine Interaction, Haptic, Focusing

In the field of tactile feedback, researchers try to generate localized stimulations on screens. Some solutions such as time reversal or phased array use vibration induced in the screen equipped with piezoelectric ceramics. We propose to use the modal basis to reproduce a specified velocity field on such devices. We explain the theory and propose a methodology to practically synthesize the voltages to achieved a controlled focusing in a given time. Experiments on a simple demonstrator are in good agreement with the theory for various velocity fields and a reduced number of modes.

Introduction

Tactile displays are now ubiquitous, the most common example being smartphones or tablets. They also have a great potential to realize virtual control panel for e.g. machine tools or automotive applications. However, one important limitation is that the main feedback channel toward the user is visual. Tactile or auditive feedback are at the moment crude, although it is commonly admitted and demonstrated that multiple sensory channels can enhance task performances and usability [START_REF] Kim | Interaction with a display panel -an evaluation of surface-transmitted haptic feedback[END_REF][START_REF] Prewett | A Meta-Analysis of Vibrotactile and Visual Information Displays for Improving Task Performance[END_REF]. Hence, there are active researches to reduce the cognitive overload induced by visual confirmation during human machine interaction by stimulating the fingertip [START_REF] Pasquero | Tactile feedback can assist vision during mobile interactions[END_REF]. In this paper, a solution that can be simply implemented on a screen is examined. Unlike pin array technologies which require specific designs [START_REF] Yang | Development of quantitative tactile display device to provide both pin-array-type tactile feedback and thermal feedback[END_REF], friction modulation is an approach where, by some mechanical action, the friction coefficient perceived by the user is modified. Two main techniques fall into this category: electro-vibration (EV) and vibro-tactile (VT) devices. Electro-vibration uses electrostatic force to attract the skin [START_REF] Agarwal | A hybrid natural/artificial electrostatic actuator for tactile stimulation[END_REF][START_REF] Bau | TeslaTouch: electrovibration for touch surfaces[END_REF]or a slider [START_REF] Yamamoto | Electrostatic tactile display for presenting surface roughness sensation[END_REF] on a surface submitted to a high voltage. Since the normal force at the contact area is modified [START_REF] Vezzoli | Electrovibration modeling analysis[END_REF][START_REF] Shultz | Surface haptics via electroadhesion: Expanding electrovibration with johnsen and rahbek[END_REF], the friction force is modulated simultaneously. This technology is versatile, however the effect can be mitigated by environmental conditions such as moisture [START_REF] Tang | A microfabricated electrostatic haptic display for persons with visual impairments[END_REF]. Vibro-tactile is already widely used in current devices where the vibrations produced by an eccentric-mass motor convey basic pieces of information. Elastic waves are well suited to realize tactile feedback, they require minimal modification thanks to the use of piezoelectric patches, and can provide interaction at distance. Hence, more sophisticated approaches have been proposed to render high fidelity vibrations. Out of plane vibrations can produce friction modulation due to intermittent contact [START_REF] Vezzoli | Role of fingerprint mechanics and non-coulombic friction in ultrasonic devices[END_REF] and squeeze film effect [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF][START_REF] Biet | Squeeze film effect for the design of an ultrasonic tactile plate[END_REF][START_REF] Wiertlewski | Partial squeeze film levitation modulates fingertip friction[END_REF]. Combination of in plane and out of plane vibration can also be used to increase the velocity at the screen/finger interface [START_REF] Dai | LateralPaD: A surfacehaptic device that produces lateral forces on a bare finger[END_REF][START_REF] Chubb | ShiverPaD: A glass haptic surface that produces shear force on a bare finger[END_REF]. Usually, resonance is taken advantage of to reduce the power requirement and achieve usable displacement. Therefore, the whole surface is active, and to render textures, the vibration must be modulated with respect to the finger position [START_REF] Amberg | STIMTAC: A Tactile Input Device with Programmable Friction[END_REF][START_REF] Meyer | high-fidelity texture display via programmable friction[END_REF][START_REF] Ben Messaoud | Textile fabrics texture: From multi-level feature extraction to tactile simulation[END_REF]. By principle, these techniques are not suited for local stimulation or to stimulate several fingers independently. Besides, as it relies on the friction coefficient, the finger must move on the surface. To address this shortcoming, phased array and time reversal of elastic waves have been proposed to achieve localized tactile stimulation [START_REF] Kim | Traveling vibrotactile wave -a new vibrotactile rendering method for mobile devices[END_REF][START_REF] Woo | Vibration rendering on a thin plate with actuator array at the periphery[END_REF]. Unlike the previous techniques, elastic waves are propagated in the screen before interfering constructively at given time and position. To do so, the transducers composing phased arrays produce the same vibration delayed by a time interval proportional to the distance separating them from the focus target. Time reversal allows to improve the focusing accuracy and resolution [START_REF] Bai | Impact localization combined with haptic feedback for touch panel applications based on the time-reversal approach[END_REF][START_REF] Hudin | Localized tactile feedback on a transparent surface through time-reversal wave focusing[END_REF]. This implies an additional knowledge resulting from an identification procedure. One limitation is that the focus pattern and its transient are not controlled, and external parameters can mitigate the effects [START_REF] Hudin | Spatial, temporal, and thermal contributions to focusing contrast by time reversal in a cavity[END_REF]. As a matter of fact, this brief overview shows that the two last techniques actually rely on similar devices, and so they could be combined. To address the robustness issue, we propose an alternative approach to localized stimulation that can be combined with friction control. It is based on a truncated modal decomposition to design a focusing shape and control its position as suggested in [START_REF] Woo | Vibration rendering on a thin plate with actuator array at the periphery[END_REF]. This paper is in the continuity of our previous work that demonstrated that the resulting focus was viable to generate stimuli [START_REF] Enferad | Differentiated haptic stimulation by modal synthesis of vibration field[END_REF]. First we recall the basic theory that links the piezoelectric patches placement to the vibration problem then the methodology to specify the time evolution. Optimal voltage distribution to avoid saturation is also discussed. In the next part, we demonstrate the validity through various experiments and address some specific issues. Finally the last sections are devoted to discussion and conclusion.

Operation principle

Assumptions

In order to illustrate the idea presented, a special case is considered, but the method actually uses the orthogonality of the elastic vibration mode shapes [START_REF] Meirovitch | Fundamentals of vibrations[END_REF] and therefore it can be generalized to other geometries. Since the point of this section is to explain how several modes are excited simultaneously in a controlled manner, a simple model involving a rectangular plate (the substrate) with piezoelectric patches glued symmetrically (bimorph regions) is considered. The volume occupied by the substrate is Ω and its surface in the mid-plane is Γ. The N b bimorph domains will be denoted Ω i and Γ i (i ∈ {1 . . . N b }) respectively for their volume and their surface in the mid-plane. Finally, the remaining substrate without patches is in

Ω 0 = Ω \ ∪ N b i=1 Ω i , defining the set Γ 0 = Γ \ ∪ N b
i=1 Γ i as it intercept the mid-plane. For each sub-domain, we consider a Kirchhoff-Love displacement field [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. That is, the components of the displacement field for sub-domain i ∈ 0 . . . N p are:

u iα (x 1 , x 2 , x 3 , t) = -x 3 w i,α (x 1 , x 2 , t) for α = {1, 2} (1)
u i3 (x 1 , x 2 , t) = w i (x 1 , x 2 , t) (2)
In this paper, partial spatial derivatives are denoted ∂f (x 1 ,x 2 ,t)

∂x 1 = f ,1 (x 1 , x 2 , t),
time derivative by dots e.g ∂f (x 1 ,x 2 ,t) ∂t = ḟ (x 1 , x 2 , t). Vectors are denoted in bold and v is the transpose of a vector v. Matrices are in bold between brackets. In the bimorph regions, the bonding is supposed to be perfect, hence the piezoelectric ceramics and the substrate share the same displacement field. Using engineering notation, the non-nil components of the strain tensor for the sub-domain i are:

@Γ Γ 0 @Γ i Ceramic Mid-plane U i h c h s x 1 x 2 x 3
S i1 = -x 3 w i,11 S i2 = -x 3 w i,22 S i6 = -2x 3 w i,12 i ∈ {0 . . . N b } (3)
The substrate is supposed to be linear hence the stresses are given by T si = [c]S i , where [c] is the stiffness tensor.

Considering now the piezoelectic ceramic, linear piezoelectric constitutive laws relating strain and stresses to the electric field E and the electric displacement vector D are used. Considering the electric enthalpy, they are given, for the considered subdomain i ∈ {1, . . . , N b }by [START_REF] Tiersten | Linear Piezoelectric Plate Vibrations Elements of the Linear Theory of Piezoelectricity and the Vibrations Piezoelectric Plates[END_REF][START_REF] Ikeda | Fundamentals of piezoelectricity[END_REF][START_REF] Mitchell | A refined hybrid plate theory for composite laminates with piezoelectric laminae[END_REF]:

D pi = [e]S i + [β S ]E i (4)
T pi = [c E ]S i -[e] E i (5)
for the direct and converse effects respectively. In practice, the piezoelectric patches used to exploit the transverse effect are thin and electroded on the upper and lower surfaces. Since the electric potential is usually prescribed on the electrodes, and assuming that thickness ratio between ceramics and substrate is small, it is commonly admitted that the electric field can be considered to be imposed by the voltage U i applied across the electrodes with a good approximation [START_REF] Heyliger | Exact freevibration analysis of laminated plates with embedded piezoelectric layers[END_REF][START_REF] Mitchell | A refined hybrid plate theory for composite laminates with piezoelectric laminae[END_REF]. Hence, it will be considered that

|E i3 | = |U i | hp .
As depicted in fig. 1, the polarizations of the patches and the voltage must be combined so that the piezoelectrically induced stress in each patch results in a moment. Finally, the displacement field is continuous, that is the out of plane displacement and the rotation at a boundary between two sub-domains i and j should be equal that is:

w i = w j ( 6
)
w i,α = w j,α α ∈ {1, 2} (7)
As for the excitation force, only the effects of the piezoelectric patches are included in the following.The excitation forces are supposed to be produced only by the piezoelectric patches.

Equations of vibrations

In the frame of the assumptions above, the equations of motion can be derived by applying Hamilton principle and are similar to the classical plate theory with segmented piezoelectric patches [START_REF] Vel | Analysis of piezoelectric bimorphs and plates with segmented actuators[END_REF]. The dynamic equilibrium of a subdomain is given by :

M i ẅi + M i αβ,αβ = 0 in Ω i (8)
where M i the masses per surface unit and M i αβ,αβ are the moment in the plate, which expressions are given in the appendices. The overall displacement field is continuous, hence on a boundary ∂Γ i , i ≥ 0 that defines a bimorph domain, the following kinematic constraint is verified:

w i -w 0 = 0 (9)
w i,n -w 0,n = 0 [START_REF] Tang | A microfabricated electrostatic haptic display for persons with visual impairments[END_REF] n indicating the normal direction of the local frame defined on the path Γ i .

To fulfil the equilibrium on the boundary the additional conditions must hold

M 0 nn,n + 2M 0 nt,t = M i nn,n + 2M i nt,t (11)
M 0 nn = M i nn -M U i (12)
t denoting the tangential direction of the local frame defined along the boundary. In Eq. 12, the moment M U i is the moment produced by the ceramics under the voltage U i , which results from the variational procedure [START_REF] Benjeddou | A Unified Beam Finite Element Model for Extension and Shear Piezoelectric Actuation Mechanisms[END_REF]:

M U i = κ b U i with κ b = e 31 d (13)
At the discontinuity points of Γ i , the conditions (11) degenerate, yielding the Kirchhoff's corner load. Assuming that the material axes are parallel to the structural axes, this equivalent corner load is:

M 0 12 (x c ) = M i 12 (x c ) (14)
at the corners at coordinates x c . It can be noted that, as a result of the previous assumption, there is no contribution of the piezoelectric effect [START_REF] Deraemaeker | Equivalent loads for twodimensional distributed anisotropic piezoelectric transducers with arbitrary shapes attached to thin plate structures[END_REF].

On the external boundary Γ, one has:

w 0 = 0 or M 0 nn,n + 2M 0 nt,t = 0 ( 15
)
w 0,n = 0 or M 0 nn = 0 [START_REF] Chubb | ShiverPaD: A glass haptic surface that produces shear force on a bare finger[END_REF] and at the corners of Γ, the condition (14) becomes:

M 0 12 (x c ) = 0 (17)

Projection

The modal analysis suppose to cancel the exogenous sources in the previous equation, hence U i = 0. This defines the Sturm-Liouville problem of the device. It turns out that the problem is merely a laminated and staged elastic plate for those conditions. It is well known that such problems have eigenvalues, the resonant frequencies ω k , and eigenvectors, the mode shapes that will be denoted ϕ k .for the modal analysis the exogenous sources are cancelled, that is U i = 0. So the problem is merely a laminated and staged elastic plate, for which it is possible to find resonant frequencies ω k and mode shapes denoted ϕ k . If the latest are normalized according to the vibrating mass, they verify the orthogonality condition:

N i=0 Ω i M i ϕ k ϕ l dΩ = δ kl ( 18
)
where δ kl is the Kronecker symbol, and ϕ k denotes the normalized mode k.

Let M ik αβ be the moments for the modecalculated with the mode shape ϕ k on the domain Γ i , it follows from the previous property that

N i=0 Ω i M ik αβ ϕ l,αβ dΩ = δ kl ω 2 k ( 19
)
The mode shapes ϕ k constitute a basis, and therefore any shape w(x 1 , x 2 , t) respecting the kinematic conditions on the boundaries can be decomposed as a weighted superimposition on the mode shapes i.e.

w(x 1 , x 2 , t) = ∞ k=1 ϕ k (x 1 , x 2 )η k (t) (20)
where the modal coordinates η k depend solely on time. Projecting the equations (8) verified by w as follows

N b i=0 Ω i M i ẅ + M i αβ,αβ ϕ l dΩ = 0 (21)
Performing the integration by parts and using properties ( 18), [START_REF] Ben Messaoud | Textile fabrics texture: From multi-level feature extraction to tactile simulation[END_REF], results in a set of decoupled dynamic equations that govern the modal coordinates where the contribution of the piezoelectric patches to the mode appears as the work produced by the piezoelectric moment for the rotation around the tangential axis along the periphery of the ceramic. It follows from this equation that several voltages applied to several ceramics can participate to the excitation of a given mode, and reciprocally, that a ceramic can contribute to excite several modes.Indeed, since U i is constant in a ceramicthe contribution to the modal force exciting a mode k is

ηk + ω 2 k η k = N b i=1 ∂Γ i M U i ϕ k,n dΓ (22)
∂Γ i M U i ϕ k,n dΓ = κ b ∂Γ i ϕ k,n dΓ U i = α ik U i (23)
where α ik , the participation coefficient of a ceramic i to a mode k per unit voltage, is fixed by the design, in other words the location, the dimensions -and ideally the shape -of the patches.

Strategy

The two results given by ( 20) and ( 23) are now exploited to specify a velocity field and a transient. The problem can be stated as sketched in fig. 2. The first step is to define the modal coordinates. Therefore, given a reference velocity field ẇ that should be realized at a given focusing time t f , we first compute its velocities spatial spectrum by:

η k = N i=0 Ω i M i ϕ k ẇ dΩ (24)
This specifies only the spectrum to be realized at the specified focusing time t f , that is ηk (t f ) = η k , but the time evolution is still to be defined. Arbitrarily setting the transients for each modal coordinates is however difficult as it would probably result in unrealistic dynamic involving very high voltages. In this paper, we adopt a different strategy based on the following observations:

1. Exciting an oscillator at its resonance frequency requires a minimal excitation force; 2. Rewriting the eq. ( 22) as a function of the modal velocities coordinates and taking the Laplace transform gives

L [ ηk ] (s) = s s 2 + ω 2 k L [ψ k ] (s)
where L is the Laplace transform with respect to time, and s the Laplace variable and the modal force ψ k is defined according to eq. 22

ψ k = N b i=1 α ik U i
This last equation shows that the modal velocity of a mode k results from applying a band-pass filter to the modal force that is produced from a linear combination of voltages. Therefore, it is to possible to selectively excite a mode by using a time-varying modal force such that its spectrum includes the relevent frequency ω k . For this approach to be valid however, the modes should not be too close in the frequency domain or coupled because of the damping. Hence, the following discussion holds for weakly damped structures.

To exploit these consideration, a modal force is now researched as a function modulated by an harmonic function at ω k as sketched in fig. 2. The purpose is to impose the envelop of the excitation of each mode independently to gain some control on their transient. Supposing that the frequency decoupling assumption holds, consider the modal dynamic equation where the damping ξ k of the mode is now introduced:

ηk + 2ξ k ω k ηk + ω 2 k η k = ψ k (25)
and suppose, as an illustration, that the desired time evolution of the modal velocity is ηk

(t) = A k e γ k a k t sin (b k t)
where

a k = ξ k ω k , b k = 1 -ξ 2
k ω k are the damping factor and the natural frequency of the mode, A k and γ k ≥ 0 are some factor that scales the transient. Using the Laplace transform for instance, one calculates the following expression for the modal force:

ψ k (t) = Ψ k e γ k ξ k ω k t [B k sin (b k t) + C k cos (b k t)] + D k ( 26
)
The details of the various constants B k , C k and D k can be found in the appendices.

To define the transient, γ k and A k still have to be fixed. This is done in this work by setting the time ∆t k = t f -t k , where t k is the starting time of the transient, allocated to each mode to reach its final value. It is also imposed that

ω k ∆t k = π 2 + 2πN k (N k ∈ N).
This ensures that all the modal coordinates reach their maximum value at the end of the focalisation (this constraint is actually a simple way to ensure that if only one ceramic is used, the maximum voltage is easily handled). The parameter γ k is used to ensure that the transient is below some level during a given time in view of the haptic interaction. Finally A k is calculated by:

A k = η k exp(γ k a k ∆t k ) (27)
In fig. 3, the case N 1 = 6 for γ k ∈ {1, 10, 50} is presented in a simulation (unit are normalized with regard to the mass i.e speed are in kg -0.5 mm s -1 and force in kg -1/2 N. The dynamic model is constituted by two modes having the same modal damping ξ 1 = ξ 2 = 0.01 and such that ω 2 = 2ω 1 .

The simulation is ran applying ψ 1 to the first mode. The velocities for each modes and the input are depicted. It can be noted that the first mode is indeed controlled (fig. 3-a) in every case. In this example, for γ 1 = 1, the offset D 1 is large (fig. 3-c), and produces an undesirable transient effect on the second mode (fig. 3-b), which would be a problem if both modes where to be controlled simultaneously. For γ 1 = 10, the transient is modified in a such a way that the amplitude remains below half the amplitude at focusing time for nearly 95 %, with limited effect on ϕ 2 . Finally, if the first mode is too strongly boosted, the input is very rapidly increasing at the end of the transient, meaning that its spectrum is spread. This induces a "spill-over" of the excitation that now also energize the second mode. Therefore, caution is required when accelerating a mode for this reason, but also because the modal force -that is the voltage required-increases. Indeed, practically, the voltage applied should not exceed some limits. One possibility is to take advantage of the multiplicity of the actuators to distribute the voltages as evenly as possible among them. To do so, the voltages U i are rewritten as a linear combination of the modal forces calculated above:

U i (t) = nm k=1 m ik ψ k (t)
where the weights m nk are used to allocate the desired modal force among the ceramics. Thus, the modal dynamic equation is now:

ηk + ω 2 k η k = N b i=1 α ki nm k=1 m ik ψ k (t) (28)
Assuming that the modes can be excited separately owing to the frequency selectivity, we require the conditions:

N b i=1 α ki m ik = 1 for k = {1, . . . , n m } (29)
to ensure that the necessary modal force is are indeed applied. Then, the weight can be tuned to reduce the voltages U i actually applied to the ce-ramics by stating the following optimization problem: minimize

m ik M Ψ subject to |U i | ≤ U max , i = 1, . . . , N b .
and to

N b i=1 α ki m ik = 1 for k = {1, . . . , N m } ( 30
)
where M is the matrix of the weights (which coefficients are the m ik , dim(M ) = N b × N m ). Ψ is a matrix which k-th row vector is constituted by ordered samples of the modal forces ψ k (nT s + t 0 ). The number of samples N s and the sampling should be chosen such that N s T s + t 0 = t f , and to cover a time interval where the modal forces ψ k become large with respect to the voltages rating. In this work, the Frobenius norm was used.

The results in fig. 4 illustrate the procedure. In this simulation, the voltages required with and without optimization are compared in fig. 4-a, together with the related transient are presented fig. 4-b. The same system as the previous simulation is considered, but this time, to model two ceramics, different modal piezoelectric coefficients are considered. In this scenario, the voltage limitation is set to 20 V, the two modes dynamics are accelerated by choosing γ 1 = 4 and γ 2 = 2. It is also required that N 1 = 10 and N 2 = 15, which results in a large offset for the first mode especially (dashed-black circles curve), the second mode (dashed-white diamonds) is less impacted because of the larger number of cycles allowed. The optimisation is performed considering the two last oscillations of the first mode (which exceed 40 V peak to peak). The simulation shows that the procedure indeed imposes new voltages that respect the given limitation almost everywhere during the transient. The transient (bottom) are quasi identical as expected.

Experimental results

Set-up

The experimental set-up is depicted in fig. 5. The device is fixed on a x-y table in order to position it under the vibrometer laser head (OFV 505 Polytec) to measure the vibration velocity at various locations. All measurements are monitored and recorded using an oscilloscope (picoscope 3403D pico technolgies LTD). The voltages are pre-calculated using the procedure described earlier, and sent to a two channel arbitrary waveform generator (Tektronix AFG3022B). These voltage references are then amplified thanks to two analog amplifiers supplying bipolar voltages up to 150 V (NF HSA

4051

). Since the phenomenon under investigation is essentially a transient, measurements and waveform generation are synchronized using an external synchronization clock. It generates a square waveform at a low frequency that is chosen to let the whole process of focalisation then decaying happen, in order to repeat the measurements starting from the same initial condition. The waveforms are acquired using the same signal and to ensure that the lag between the trigger signal and the actual generation is known, the clock signal was also measured. In practice, this lag is negligible. The voltages featuring multiple frequencies up to 100 kHz, it is required to respect the Shannon frequency. In the case of the generator used, for a focalisation time of 50 ms, the equivalent sampling frequency with the 2 16 points available is 1.3 MHz which verifyies the criterion. The waveforms specification and the measurements processing are performed by homemade software using matlab and the libraries provided by pico technology and tektronix. The demonstrator is a clamped staple shaped aluminium beam equipped with piezoelectric patches Noliac NC51 with dimensions (9 × 5 × 0.5) mm 3 as visible on the picture in fig. 6. This peculiar shape was chosen because i) it reduces the surface to be scanned hence the velocities cartography are considerably faster to realize; ii) it demonstrates that the method is usable even if actual solutions to the eigen problem are not known. As discussed above, the ability of a voltage U i to excite the k-th mode depends on the α ik coefficient eq [START_REF] Bai | Impact localization combined with haptic feedback for touch panel applications based on the time-reversal approach[END_REF] which depends on the geometry and location of the patchess. To avoid to use different ceramics, the electrodes were actually splitted (e.g second ceramic from the left). By properly chosing these cuts, it is possible to evenly excite the modes. The identification was realized using a white noise identification procedure. The velocity vs input voltage frequency response fonction was recorded at points defined on a regular grid on the upper part of the beam. The mapped region was defined between the limits of the curved parts, with a step small enough with regard to the smallest spatial wavelength. During the process the end ceramics were used in conjunction or in opposition, to excite the even and odd modes respectively, because they are able to excite all the bending modes in the 100 Hz-80 kHz band. During these measurements, the other ceramics are in short-circuit conditions. Modes were selected using the resonance peaks. Knowing the phase and amplitude at resonant frequencies at each measured point, the mode shapes were reconstructed. When the signal to noise ratio was too low, the mode shapes were more precisely mapped at resonance choosing ceramics with more effective participation coefficients.

In fig. 7, the modes ϕ 1 to ϕ 15 (up to 30.9 kHz) are depicted in grey. The orthogonality of the measured modes, however, is only approximatively verified, as shown in fig. 8 on the left, where the result of the normalized projection of modes on each other reveals cross talk (e.g between 16 th and 17 th mode). This results from the mapping that does not cover the whole device, but only a rectangular area on the top. Indeed, finite element simulation showed that for low frequency modes (having large wavelength) the corners of the staple where acting as clamps, i.e. the modes were located on the top of the staple only. However, as the frequency increases, this is not true any more, and modes tend to occupy the whole staple. This is also verified experimentally. To address this problem, the modes are orthogonalized using a Gram-Schmidt procedure which gave significant improvement (fig. 8 right) although the shapes are only slightly modified (depicted in color in fig. 7). The time synthesis of the voltages discussed previously requires a fair knowledge of the damping of a mode, which was done by measuring the -3 dB bandwith. Finally, the participation factors were estimated from the measurement of the velocity at resonance obtained when supplying each ceramic individually (the other being short-circuited).

Protocol

The main experiments presented here aimed at testing the viability of the approach to reproduce a specified shape on the top of the "staple". A sinc function, where position, width and amplitude of a relatively homogeneous high velocity region can be easily specified, was adopted. From fig. 7, one can observe that the first modes are relatively similar to harmonic functions. The clamp-clamp modal shapes of a beam can be well approximated by cosine function from the third mode. To understand the trend of the relationship between mode amplitude and the reference shape, we assimilate all the modal shape to cosine function ϕ k (x) = cos (β k x) for β k = kπ L (L:length of the beam). Here, without loss of generality, the normalization factor is omitted as it is the unity for all modes of a uniform beam . This approximation leads to overestimate the contribution of the first two modes, but provides insight of the general trends found later. Under those assumptions, the projection of the reference function is

0.5 -0.5 sinc (ax) cos (β k x) dx = Si 1 2 (kπ -a ] + Si 1 2 (kπ + a) a
where Si is the integral sine. Hence a narrow "high velocity" region will require more modes for a satisfactory approximation. On the other hand, their amplitude decreases so it should be possible to realize velocity fields with "sharp" variation despite the voltage limitation. Actually, the modal spectrum follows the same trend as the Fourier transform of sinc function. Five different focusing locations, equidistantly distributed on the beam, were tested. For the outermost points and large central lobe of the sinc function, the kinematic condition cannot be respected. So a Tukey window is used to respected the boundaries of the cartographies.

To assess the control of the mid height width of the central lobe, successive widths of 25, 18 and 10 mm were tested. Amplitude control was indirectly assessed by imposing 40 mm s -1 as the maximal velocity for every test (before applying the Kaiser window). For practical reason, the modal projection must be truncated. To select the modes, a criterion is defined as follow. Considering that twice the kinetic energy is given by:

K = Ω ρ ẇ ẇdΩ = Ω ρϕ ϕdΩ η2 k = ∞ k=1 η2 k (31)
one can define the τ -score criterion for a truncation order N m :

τ = 1 - Nm k=1 η2 k K 2 1 2 (32)
where K can be determined directly from the specified velocity field ẇ. The τ -score is thus a measure of the residue resulting from the truncation. For each test, the velocity field was measured, projected on the modal basis and compared to the specified modal coordinates resulting from the truncation. The criterionυ-score used to assess the "misalignement" of the obtained shape compared to the specified one is defined as:

υ = cov [ ẇ (x m ) , ẇ (x m , t f )] var [ ẇ (x m )] × var [ ẇ (x m , t f )] (33)
where var and cov are the variance and the covariance respectively, and x m are the location where the velocities are sampled. The criterion is chosen to evaluate the discrepencies between the truncated modal projection and the experimental modal projection, regardless of the error on the amplitude which can be due to a miscalibration of the amplification gain. In this way, one can appreciate whether the obtained shape is similar to the (truncated) reference up to a scaling factor.

Results

Voltages

The strategy proposed to comopute the voltages applied is examined in this section. As explained in section 2.4, the voltages deduced from the prescribed modal forces given by eq. ( 26) can be "mixed" between the patches using a linear combination. Here, the scenario examined is to control two modes simultaneously, and force them to achieve a given set of velocities Vel. Meas (V) Comparison between predicted and measured velocities at the center of the staple: measurements (grey), filtered measurements (plain line, black) and predicted (dashed line, black). The gain of the velocity measurement is 5 V s mm -1 . Zoom at focus time is depicted in the frame. of 20 mm s -1 for the modes ϕ 1 and ϕ 3 at the center of the experimental setup (for simplicity, the modes are renormalize in order to have 1 mm s -1 amplitude for the maximum velocity of a given mode shape). Therefore, the expected velocity at the center is 40 mm s -1 . Moreover, a speed up of γ 1 = 4 and γ 3 = 1 over the natural decaying speed of modes ϕ 1 and ϕ 3 is imposed. Finally, ϕ 1 and ϕ 3 should achieve the focusing within N 1 = 10 and N 3 = 25 oscillations respectively. For the purpose of the demonstration, two set of patches were selected. The first one (denoted set 1 in the sequel) has a strong effect on ϕ 1 , but a very weak effect on the second, while the second set is more balanced with a prominent effect on ϕ 3 . This results in the voltages depicted in fig. 10 (bottom). In fig. 9, a comparison of the simulated and measured velocity transients is presented. The agreement is good despite the uncertainties on the dynamic parameters. To assess the effect of "mixing" the voltages, a second run consisted in applying the optimisation procedure in order to conceal the voltages amplitude below 50 V as long as possible during the transient of the excitation signal is well controled since the second, 4th and 5th modes (at f =963, 3907 and 5719 Hz)are not excited. The transients are still very similar as expected, and the velocity at focusing time is still realized.

Location

The results from tests to control the location at x =-46, -24, 0, 24 and 46 mm of the high velocity zone (25 mm of width at mid-height of the main lobe) are presented in fig. 12. For each test, the specified, the approached and the measured velocities at focusing time along the center line in the y direction are presented (left column). The approached reference velocity field results from the truncation of the actual modal decomposition of the reference which was fixed here to 6 modes for all references in this case. As explained earlier, the reference field is windowed with a Tukey window Figure 11: Spectra of the voltages and velocity at the center for case 2: since set 1 of ceramics is poorly coupled to the 3rd mode, the optimization requires that the set 2 is used also to excite mode 1, hence reducing the voltages as shown in fig. 10 having a sharp transition to ensure a nil value at the boundary, without causing too much distortion of the initial reference. Similar test where ran for sharper pulse. The τ -score of the various truncation are summarized in Tab. C.1 in the appendix Appendix C. One can observe that the higher the truncation the better the approximation of the velocity field, the worst cases being on the ends. This is due to the need to "cut-off" the actual waveform, which, even with the use of the Tukey-window, requires a wider spatial spectrum. This is clear from the improvement of the score with the reduction of the width of the pulse. For each test, there is a good agreement of the approached and the measured velocity field, despite the imperfect knowledge of the damping and the participation factors. At the ends of the beam, it would require more modes to achieve the stiff edge of the reference. Hence, the assigned maximum velocity cannot be achieved, still the selected truncation order is high enough to achieve a correct reconstruction (see Tab. C.2).

Width

The results are summed up in fig. 13, with the same convention as before. As already explained, reducing the width of the main lobe of the sinc function, results in a smaller and wider spectrum in the modal basis and is therefore more challenging. One can observe on the spectrum that Figure 12: location test at -46, -24, 0, 24 and 46 mm, on the left velocity fields for the references, the approached ones resulting from the truncation and the measured velocities (resp. black, grey and white) and their respective modal projection on the right. Note that for the ends velocity fields a Tukey window is applied to respect the kinematic boundary conditions.

Figure 13: Width test results with the reference (black), the truncated reference (grey) and measured (white), the velocities and their modal projection are depicted in the left and right columns respectively. higher order modes are not correctly excited. This is partly explained by the high speed transient of these modes: the pre-calculated signals also include the lower frequency transient which are much slower, hence a long signal (compared to the actual transient time of the high frequency modes) is needed. This results in a coarse sampling given the limited number of points that can be programmed into the generators. However the effect is not too critical on the actual velocity field.

Transient

The transient can be independently controlled by imposing the transient of the voltage harmonics. In fig. 14, the transients of the five first modes are depicted in the case where ϕ 1 is strongly accelerated. These measurement are obtained by projecting the cartography of the transient on the modal basis. In the test presented, the references forced regime of the modal coordinate during the forced period are in gray, the theoretical decaying is in dashed gray, and the actual measurement are in black.

Because the first mode has a very high quality factor, its transient is much slower compared to the other modes. This interferes with the intended tactile effect, so it is desirable to apply the development previously discussed. In this case, the acceleration factor is 20 (see the spectrogram in fig. 15), which is close to the limit of the method (spill-over). Comparing the growth to the decay, there is a substantial effect for the first mode. Some spillover resulting from the fast dynamic excite the second and fourth mode. There is also some lag on the first mode that is the result of a mismatch of the identified and the actual resonant frequency. Since the quality factor is close to 190 (for comparison the 3rd and 5th modes quality factors are 47 and 63 repectively) a slight error on the resonant frequency (possibly due to variation of the environmental conditions) introduces a phase error. The benefit of the approach is more clear on the spectogram of the velocity at the point of focalisation. The modes are correctly excited; all the modes grow simultaneously until the focusing time at 25 ms. This also shows that for the application the decay -which were left here for demonstration-should also be controlled in practice.

Discussion

The price of the versatility of the approach is the identification step. However, the results above are encouraging as they demonstrate that, despite the uncertainties on the parameters of the modes and their modal shape, it is possible to realize a given velocity field with an acceptable precision. So the most salient difficulty is the number of modes required which increases the number of uncertainties as demonstrated by fig. 13 but also, makes it more challenging to remain within the voltage limits. This can be addressed by a proper choice of the actual velocity profile and the truncation. For instance, the graphs in fig. 16 depict another example where the pattern is defined as

ẇ(t f ) =    A exp - 1 1-( x-xc a ) 2 if x ∈ (-a + x c , a + x c ) 0 otherwise (34)
known as the bump function and where A is a scaling factor, x c the center coordinate, and a controls the width. This function has interesting continuity properties and its Fourier transform is rapidly decaying [START_REF] Johnson | Saddle-point integration of C ∞ "bump" functions[END_REF]. As mentioned earlier, there exists an analogy between the Fourier transform and the modal projection in this case, hence the result presented in fig. 16 was expected, despite the difference between the actual modes and the

Conclusion

This paper has discussed a method to synthesize and approximate velocity field in order to realize a specified velocity field at a given time. The scope of the paper is haptic application where differentiated stimuli are required. Given that such devices are generally realized using plates excited with piezoelectric patches, the equation have been recalled to demonstrate the capacity to excite specifically several modes by a applying well defined voltages to the different ceramics. We have proposed a method to calculate the voltages that must be applied and dispatch them among the piezoelec- So the modal force is:

ψ k (t) = Ψ k e a k ξ k ω k t [C k sin (b k t) + D k cos (b k t)] + E k with : Ψ k = A k a k (γ k + 1) a 2 k γ 2 k + b 2 k C k = a 2 k γ k (γ k + 1) + 2b 2 k D k = a k b k (γ k -1) E k = b k a 2 k + b 2 k a 2 k γ 2 k + b 2 k
Since in this work ω k ∆t k = π 2 + 2πN k is imposed, A k is such that:

η k = A k e γ k ξ k ω k ∆t k
Note that the differential equation could be directly used once the envelop is defined.

Appendix C. Experimental results

In this section, the τ and υ scores of the various shapes tested in this paper are summarized.

Figure 1 :

 1 Figure 1: Schematic of the modelled device with the electric connexions of a ceramic pair (above) and the mechanic layout (below). Only one ceramic pair is represented for clarity.

Figure 2 :

 2 Figure 2: Schematic of the method: the projections of the desired velocity field yield the final modal coordinates at focusing time in the modal space, their evolutions must then be specified by individual time functions between the rest and the focus configurations

Figure 3 :

 3 Figure 3: Example of transient in the case where mode ϕ1 (a) is excited by the modal force (c) to impose an transient which growth is controled by the factor γ. The neighboring mode ϕ2 (b) can be excited either by either by the offset (for γ = 1 or too fast a transient γ = 50.

Figure 4 :

 4 Figure 4: Effect of allocating the modal forces among two ceramics. Voltages U1 and U2 (top) needed in the case of assignement of one ceramic to excite one mode are in dashed, and in plain for the case of optimized allocation. The transient are reaching the desired modal speed in both case (bottom).

Figure 5 :

 5 Figure 5: Schematic of experimental setup.

Figure 6 :

 6 Figure 6: Picture of the experiment: A) Staple shaped aluminum beam, the velocity is controlled within a rectangular region on the upper part; B) Fixture of the beam to realize clamp-clamp boundary conditions; C) Mirror showing the ceramics; D) Reflection revealing some of the ceramics fixed beneath the beam. The second ceramic from the left is split as explained in the text; E) Electrical connections of transducers allowing the reconfiguration according to the set of ceramics used; F) Laser spot positioned at one of the corner of the studied region.

Figure 7 :

 7 Figure 7: Experimental mode shapes (grey) and reorthogonalized mode shapes (colored) in the frequency range 0 to 39 kHz

Figure 8 :

 8 Figure 8: Cross projection between experimental modes reveals their imperfect orthogonality (left) which can be rectified using an orthogonalization procedure (right)

Figure 9 :

 9 Figure9: Comparison between predicted and measured velocities at the center of the staple: measurements (grey), filtered measurements (plain line, black) and predicted (dashed line, black). The gain of the velocity measurement is 5 V s mm -1 . Zoom at focus time is depicted in the frame.

Figure 10 :

 10 Figure 10: Validation of the voltages allocation strategy: (a) velocity measurement with "one mode one patch" strategy (case 1 black) and mixed voltages between patches (case 2 dark grey). The corresponding voltages applied to the patches set are depicted in the case 1 (b) and case 2 (c).

Figure 14 :

 14 Figure 14: Transients of the five first excited modes used to synthetize a 25 mm wide sinclike velocity field. The theoretical forced regime (red) is followed by a free decay (green), compared to the actual measured transient obtained by (blue).

Figure 15 :

 15 Figure 15: Spectrogram of the velocity measured at the focusing point.

Figure 17 :

 17 Figure 17: Test results of approximated bump velocity profiles at location points -46, -24, 0, 24 and 46 mm, on the left velocity fields for the references, the approached ones and the measured velocities (resp. red, black and blue) and their respective modal projection (right). The truncation includes 7 modes, based on the distortion criteria.

  Table C.1: Effect of the truncation: τ -score for the various shapes studied. Table C.2: Quality of the focusing : υ-score for the various shapes studied.

	Pulse center	P 1	P 2	P 3	P 4	P 5
	Sinc 25 mm	0.6831 0.9418 0.9755 0.9359 0.7124
	Sinc 18 mm	0.8833 0.9868 0.9838 0.9887 0.9757
	Sinc 10 mm	0.9521 0.9948 0.9897 0.9928 0.9981
	Bump 25 mm 0.8285 0.9727 0.9544 0.9816 0.8717
	Pulse center	P 1	P 2	P 3	P 4	P 5
	Sinc 25 mm	0.9739 0.9739 0.9739 0.9739 0.9739
	Sinc 18 mm	0.9792 0.9885 0.9902 0.9611 0.9945
	Sinc 10 mm	0.8771 0.9768 0.9759 0.9386 0.9739
	Bump 25 mm 0.9877 0.9953 0.9911 0.9609 0.9973

Acknowledgements

This work has been carried out within the framework of the project StimTac of IRCICA (institut de recherche sur les composants logiciels et matériel pour la communication avancée).

tric patches based on the physical limitation of the ceramics. To do so, the problem is recasted in the modal basis, and the modal forces are specified based on the projection of the projection of the desired velocity field. The advantage of the method is the possibility to generate stimuli independently from the actual location of the actuators and to control the time needed for focusing by properly controlling the rise time for each mode. Experiments have confirmed that it is possible with a reduced number of modes. Several velocity fields have been tested for different velocity fields. Future work will address the problem of parametric variation by introducing a close loop control, and psycho-physical studies to figure out the relevant velocity distribution for haptic applications.

Appendix A. Reduction elements of the plate model

The model is based on the standard Euler-Bernouilli kinematic assumptions. Since the device is subdivided into domains including either the substrate (index or superscript 0), or the bimorph (index or superscript greater than 0, within {1 . . . N b } where N b is the number of the bimorph domain), the different expressions are recalled below. The mass per surface unit are:

the moments are given by:

2)

with the bending stiffnesses:

hp 2 is distance between the ceramics and the substrate midplane.

Appendix B. Computation of the modal forces

The procedure to calculate the modal force discussed supposes to specify the transient in the form:

f k (t) defines the envelop of ηk (t). Presently, a increasing exponential is chosen with the tuning parameters A k and γ k :

Using the Laplace transform, one can write:

yielding:


T Qoria

F Gruson

F Colas

X Guillaud

M-S Debry

T Prevost

Tuning of cascaded controllers for robust grid-forming Voltage Source Converter

Keywords: Transmission system, Grid-forming control, Parameters of Controllers, Tuning algorithm, Small-signals analysis

From the origin of the grid, energy has been delivered to electrical loads mainly by synchronous machines. All the main rules to manage the grid have been based on the electromechanical behavior of these machines which have been extensively studied for many years. Due to the increase of HVDC link and renewable energy sources as wind turbine and PV, power converters are massively introduced in the grid with a fundamentally different dynamic behavior. Some years ago, they were connected as simple power injector. Then, they were asked to provide some ancillary services to the grid, in the future, grid forming capability will be required. Even if gridforming converters had been extensively studied for microgrids and offshore grids, it has to be adapted to transmission grid where the topology may be largely modified. This paper presents an algorithm for calculating the controller parameters of a gridforming converter which guarantee a stable behavior for many different configurations of the grid.

I. INTRODUCTION

The rapid increase in the production of electricity from renewable energies (solar, wind) as well as the interconnection of synchronous or asynchronous zones by high voltage direct current (HVDC) links [START_REF] Rault | Introduction au contrôle commande des réseaux HVDC Multiterminaux[END_REF] leads the electric power transmission grid to undergo an unprecedented mutation. Indeed, photovoltaic or wind-based production systems are interfaced to the alternative transmission system through power electronics devices, and HVDC interconnections also use static converters. Then, the penetration rate of power electronics in transmission grid is steadily increasing, which leads to some fundamental modifications in the grid since their physical behavior is fundamentally different from the classical synchronous machine. When the renewable energy started to grow, these sources were considered as "negative load" meaning that the grid had to cope with new kind of loads which were producing more power than consuming. This is the base of grid following mode [START_REF] Ottersten | Vector current controlled voltage source converter-deadbeat control and saturation strategies[END_REF] [START_REF] Svensson | Vector current controlled grid connected voltage source converter-influence of nonlinearities on the performance[END_REF] where the converter injects power in the grid and needs a PLL to be synchronized to the grid voltage [START_REF] Zhang | Power-Synchronization Control of Grid-Connected Voltage-Source Converters[END_REF]. By increasing the power electronic converter penetration, a grid supporting capability has been asked for [START_REF] Allal | A survey on control of electric power distributed generation systems for microgrid applications[END_REF], mainly in term of voltage and sometime in term of frequency. The way to connect this kind of converter is still the same: an active and reactive power injector with an added function of voltage or frequency regulation as in classical electrical power sources. The major disadvantage of this mode of connection is the impossibility to operate in islanded mode and stability issues in case of a lower Short Circuit Ratio (SCR) mainly due to the PLL dynamics [START_REF] Zhou | Stability of DC-link voltage as affected by phase locked loop in VSC when attached to weak grid[END_REF] [START_REF] Zhou | Impact of Short-Circuit Ratio and Phase-Locked-Loop Parameters on the Small-Signal Behavior of a VSC-HVDC Converter[END_REF]. In the future, it is expected that large geographical areas will be fed at some time of the year only by power converters [START_REF] Debry | Maximizing the penetration of inverterbased generation on large transmission systems: the MIGRATE project[END_REF], so a power injector behavior for the converters, as it is now, will not be able to manage this kind of situation. In this case, the VSCs may need to operate autonomously, the power converter will have to generate the voltage waveform . This is the concept of grid forming converter [START_REF] Rocabert | Control of Power Converters in AC Microgrids[END_REF]. Due to the physical nature of the converter, it will not control only the voltage RMS value, but also the instantaneous waveform. This control technique has been used in UPS [START_REF] Guerrero | Control strategy for flexible microgrid based on parallel lineinteractive UPS systems[END_REF] and extensively studied in micro-grid [START_REF] Wu | A Virtual Inertia Control Strategy for DC Microgrids Analogized With Virtual Synchronous Machines[END_REF]- [START_REF] Katiraei | Microgrids management[END_REF]. But transmission grids have their own specificities compared to the distribution grid which induces specific studies to guarantee the same quality of service as today [START_REF] Denis | Review on potential strategies for transmission grid operations based on power electronics interfaced voltage sources[END_REF]. Among them, the robustness to various topologies is quite challenging since the transmission grids are subject to various events as loads connection, trip of lines, short circuits, etc. The high controllability of the power electronic converter may be an advantage but also a drawback since a high dynamic voltage control may deliver a very good voltage quality to the grid around an operating point but also may suffer of a poor robustness when subjected to different situations previously mentioned. The robustness improvement of the control of these converters is the main point addressed in this paper. This approach has to take into account the specificities of converters connected to transmission grid. Indeed, since the power rating of VSCs in transmission grid can be very high, the switching frequency is rather low to limit the losses, this implies a limitation of the control bandwidth. Furthermore, the converters are connected to the grid thanks to an LCL filter, the sizing of this filter has also to be taken into account. Various types of control algorithms could be considered but to ease the implementation in an industrial converter, a classical voltage and current loop has been kept. Two solutions had been proposed in the literature to deal with this issue, the first one is based on the choice of the LCL filter parameters, e.g., as increasing on the capacitor value [START_REF] D'arco | Automatic Tuning of Cascaded Controllers for Power Converters using Eigenvalue Parametric Sensitivities[END_REF]. The second one consists in acting on the feed-forward disturbances rejection in the inner loops at the control level [START_REF] Salha | Microréseaux îlotables: étude et coordination des protections des générateurs et du réseau[END_REF]- [START_REF] Coelho | Small-signal stability for parallel-connected inverters in stand-alone AC supply systems[END_REF]. In all these solutions, the system dynamics are not very well controlled. This paper exposes clearly the instability origin due to the interaction between system modes, therefore, time decoupling is not possible due to the bandwidth limitation. The objective of this paper is to take into account the overall dynamics of the system and to find a setting of controllers parameters satisfying a good stability, damping as well as a very high robustness with respect to the various grid topology. The determination of the parameters is done via an algorithm which responds to several imposed constraints. To formalize this problem and check the system stability, small-signal analysis and participation factors tool [START_REF] Danielsen | Use of participation factors and parameter sensitivities in the study and improvement of low frequency stability between electric rail vehicle and power supply[END_REF] have been used. The paper is organized as follows. In section II, modeling of a grid-forming control for transmission power systems and recall on classical controllers tuning are presented and tested in two different grid topologies in order to check the classical tuning method limitations. In section III, small-signal analysis for grid-forming converter is discussed in order to understand the instability origin. Finally, the proposed controllers tuning algorithm is presented and tested in section IV to highlight the very good dynamics and robustness obtained with the proposed method. Conclusions are drawn in Section V.

II. GRID-FORMING CONTROL FOR TRANSMISSION POWER

SYSTEM

Because of their capability to operate in different topologies, grid-forming voltage source converters will be used for future transmission grid (MIGRATE PROJECT [23]). In this article, the studied grid-forming VSC consists of a 2-Levels VSC supplied by a DC bus and connected to the AC side through an LCL filter. As explained in the introduction, this converter has to cope with various situations which are modeled here by two extreme cases : -Standalone mode where the converter and its filter is connected only to one load (K1 open in Figure 1). -Connection to an infinite grid modeled by a Thevenin voltage source and its equivalent impedance (K1 closed in Figure 1). An LCL filter is very often used in grid-feeding converter, so the main difference is not the topology of the converter itself but its control since an instantaneous voltage loop is implemented in grid-forming control whereas grid-feeding is mainly focused on the control of the current or the power. In case of an association of several grid forming converters, a droop control [START_REF] Mohamed | Adaptive Decentralized Droop Controller to Preserve Power Sharing Stability of Paralleled Inverters in Distributed Generation Microgrids[END_REF] is added to share the active and reactive power between them. The upper part of Figure 1 shows the general organization of the considered grid-forming control.

A. System Modeling

Following the notations in Figure 1, the accessible state variables are respectively the VSC current i s , the voltage across the filter capacitor e g and the current i g through the transformer inductance L c . To simplify the system study, the VSC is represented by its average model, and the state variables are represented in d -q frame in order to control DC variables as expressed by the dynamic equations below.

v md -e gd = L f w b . i sd + R f i sd -L f w 0 i sq (1)
v mq -e gq = L f w b . i sq + R f i sq + L f w 0 i sd (2)
i sd -i gd = C f w b . e gd -C f w 0 e gq (3)
i sq -i gq = C f w b . e gq + C f w 0 e gd (4)
e gd -v pcc d = L c w b . i gd + R c i gd -L c w 0 i gq (5) e gq -v pccq = L c w b . i gq + R c i gq + L c w 0 i gd (6)
The quantities are expressed in (P er -unit) except w b which is the base frequency in (rad/s).

B. Grid-forming Cascaded-Loops Control

The grid-forming control consists of three main stages as shown in Figure 1 : • The linearization stage that provides modulation signals.

• The inner control loops : current and voltage control.

• The primary control which is constituted of the active power droop that provides the control angle θ ref and the reactive power droop that provides the d-axis reference e g dref to the voltage loop. Using the two-axis theory, the injected instantaneous active and reactive power components, p mes and q mes , are given by the following equations : p mes = e gd i gd + e gq i gq (7) q mes = e gd i gq -e gq i gd [START_REF] Debry | Maximizing the penetration of inverterbased generation on large transmission systems: the MIGRATE project[END_REF] To ensure a good dynamic decoupling between the primary control and the inner control, a filter is added on the active and reactive power measurement. Where w c is the cutting frequency [START_REF] Mohamed | Adaptive Decentralized Droop Controller to Preserve Power Sharing Stability of Paralleled Inverters in Distributed Generation Microgrids[END_REF].

p f = w c w c + s p mes (9)
q f = w c w c + s q mes (10)

C. Recall on Classical Tunning of Controllers

As the states are transformed to the synchronous reference frame, they become DC signals under balanced sinusoidal grid and perfect synchronization. Under these conditions, PI controllers are sufficient to control the state variables without error in steady state. The control system requires a decoupled control of i sd and i sq . However, the two axes are coupled due to the cross terms L f .i sq and L f .i sd . These coupling terms are compensated by a feed-forward action, so that, the d-axis and q-axis components of currents can be controlled independently. The current loop is constituted of three mains parts (Fig. 2): -A PI controller.

v m e g i s v pcc L f , R f C f L c , R c VSC 2-Levels P * Q * kpc + k ic s - + -+ + + P -1 (θ) P (θ) P W M m abc i sd i sq v g d v g q L f w L f w w * U dc i sd i sq - + -+ + + e g d e g q i g d i g q C f w C f w mp s -+ + nq -+ + P mes Q mes P (θ) i g d i g q P (θ)
-An approximated transfer function of the PWM.

-A transfer function of the R f ,L f filter.

k pc + kic s 1 1 2fsw s+1 i sref i s i smes PI controller PWM VSC Current TF + - 1 Lf s+Rf Figure 2. Current control loop
Were k p , k i and f sw are respectively the PI controller gains,and the switching frequency. For the current controller design, the approximated PWM transfer function can be neglected if the current loop is at least 10 times slower than this latter, meaning that the response time of the current loop should be chosen higher than:

T rc ≥ 30 2f sw (11)
Where T rc is the current loop response time.

2) Voltage Controller Design:

Similar to what is done for the current loop, the voltage loop is constituted of three mains parts (Fig. 3) : -A PI controller.

-An approximated transfer function of the current loop.

-A transfer function of the capacitor. In this case, if the response time of the voltage loop is 10 times higher than the current loop response time, the effect of the current transfer function may be neglected in the calculation of the PI controller parameters. All the simulation parameters are given in Table I. The response time of the current and voltage loops are respectively set to 5ms and 50ms with damping ratio of ζ = 0.7.

D. Limitation of the classical controllers tuning

Given the controller parameters calculated previously, some time domain simulations are presented for two different cases: -K1 OFF: Standalone mode -K1 ON: Grid connected mode 1) Stand-alone mode: The load is assumed purely resistive. Two events are applied to the system in order to check its behavior and stability : -Load variation from P load = 0.1 pu to P load = 0.4 pu. -Variation of 0.1 on AC voltage step of E set = 0.1 pu. The dynamic response of the AC voltage and the power correspond to what was expected.

2) Grid-connected mode: The grid is assumed very stiff (SCR = 20). The control developed in stand-alone remains unchanged but the system is now largely modified. Following an operating point change (P 0 = 0pu to P 0 = 0.4pu), Figure 6 highlights that the power reference is not controlled any more since the system is unstable. Small-signal analysis and participation factors are used in the next section in order to analyze the dynamics of the system and find the origin of the instability.

III. SMALL-SIGNAL ANALYSIS FOR GRID-FORMING

CONVERTER

To analyze the problem previously described, a mathematical formulation is developed by combining the system dynamics and the control equations. These latter can be linearized around an operating point using small-signal analysis:

∆ . X = [A l ]∆X + [B l ]∆U (12)
∆X T = ∆θ ∆P ∆Q ∆C vdq ∆C cdq ∆i sdq ∆e gdq ∆i gdq ∆X represents all the state variables of the system when the converter is connected to the grid. To the 6 variables of the LCL filter, 7 others has been added :

• Four intermediate states of PI controllers of the current and voltage loop.

• Two outputs of the filter on P and Q. It can be noticed on table II that λ 5-6 have a positive real part which does not respect the stability condition as verified on the time domain simulation (Fig. 6). To define the link between the unstable modes and system states, the participation factors [START_REF] Danielsen | Use of participation factors and parameter sensitivities in the study and improvement of low frequency stability between electric rail vehicle and power supply[END_REF] P F actor are used to determine how much the i th state variable participates in the j th mode. The participation factor for λ 5-6 is presented on Fig. 7. It reveals an important link between the unstable modes and especially the currents i sd , i sq .

The most common solutions [START_REF] Salha | Microréseaux îlotables: étude et coordination des protections des générateurs et du réseau[END_REF]- [START_REF] Coelho | Small-signal stability for parallel-connected inverters in stand-alone AC supply systems[END_REF] consist to act on the feed-forward disturbances rejection in the inner loops, e.g., in [START_REF] Denis | Improving robustness against grid stiffness, with internal control of an AC voltage-controlled VSC[END_REF] H2 is set to zero, and in [START_REF] D'arco | Automatic Tuning of Cascaded Controllers for Power Converters using Eigenvalue Parametric Sensitivities[END_REF] H1 is set to zero. These solutions are considered as the references solutions which will be compared to the one developed in this paper.

IV. PROPOSED TUNING CONTROL ALGORITHM

The proposed method assesses the gains of the controller which gives the best placement of the poles of the system calculated on the 13 rd order system meaning that the converter is grid-connected mode. Thereafter, it is verified that it is also valid for the isolated mode.

A. Methodology

The determination of controllers parameters is done following three steps:

1-At first, the eigenvalues are calculated for a large range of the four controllers parameters. 2-Two constraints are imposed on the system poles:

λ i (A l ) < 0 (13) λ i (A l ) > w 0 (14)
Where w 0 = -800 rad/s. The choice of w 0 is linked to the limited bandwidth of converters which are connected to the transmission grids.

When the poles respect these two constraints, the set of the 4 parameters is stored.

3-The third step consists in calculating the damping of each pole for each solution and comparing them. The damping ratio is calculated for different poles for all the solutions which respect the constraints. For each solution, the smaller damping is isolated (ζ imin ). All the (ζ imin ) are compared. Finally, the solution which is chosen is the one which gives the maximum damping in all the (ζ imin ). This guarantee to have the optimal dynamic behavior among all the possible solutions with the chosen set of parameters for the 2 controllers. It can be noted from Table III that the obtained controllers parameters are very different from those calculated using the classical tuning method (Table I). The two parameters that had a lot of impact on pole placement are especially the current IV shows the resultant eigenvalues of the system. The minimum damping is associated to λ 1-2 , its value is 0.21. Since the system is non linear, the choice of the operating point is important. The gains have been calculated for the nominal operating point (P0 = 1 pu).

B. Comparison with other methods

A time domain simulation is presented in order to validate the parameters derived from the proposed algorithm and compare with the two other solutions proposed in the literature: The simulation is done on MATLAB/SimpowerSystem using the system parameters in Table I.

-In [START_REF] Denis | Improving robustness against grid stiffness, with internal control of an AC voltage-controlled VSC[END_REF], the proposed solution is to set H2 to zero keeping the classical tuning obtained in standalone mode. It is shown that it improves the stability. However the dynamics of the system is not controlled anymore.

-In [START_REF] D'arco | Automatic Tuning of Cascaded Controllers for Power Converters using Eigenvalue Parametric Sensitivities[END_REF], H1 is set to zero and the gain of the controller are calculated to get the best behavior as possible but the damping is not good enough. . In order to excite the modes linked to the inner loops, an event is applied on AC voltage E set = 1.1 pu at t = 2s. On Figure 8, the results show clearly the improvement brought by the proposed method. It illustrates that the proposed algorithm is effective in ensuring the stability and improving the dynamic performance of the system compared to what have been achieved using the two other approaches.

C. Robustness of the Grid-forming Voltage Source Converter

As mentioned in the introduction, the robustness of the proposed control is a major property to be tested since the converter is subjected to various situations due to the modification of the grid topology. Four studies are realized without changing any parameters of controllers presented in table III Figure 9 analyses the evolution of poles in respect with the active power operating point. It highlights that an operating point modification has a small impact on the system stability.

2) Variation of Short circuit Ratio:

The evolution of poles is now analyzed in respect with the Short Circuit Ratio from 20, very strong grid, to 1.2 which is a very weak grid. The poles have been calculated for the nominal operating point: P0 = 1 pu. Figure 10 highlights that the variation of the Short Circuit Ratio has a stronger impact particularly on the mode related to the angle. However the system remains stable for a very wide interval.

3) Standalone mode:

In figure 11, the same tests as in Section II-D-1 are achieved: a load variation from P load = 0.1 pu to P load = 1 pu and AC voltage variation E set = 1.1 pu at t = 0.2s. These results confirm that the control strategy, designed for a grid-connected mode, is also able to operate in stand-alone mode with a very good stability margin.

4) Multi-converter grid topology:

An extension of the stand alone configuration is the a connection of two grid forming converters through a transmission line. For this test, three events are realized: -At t = 1s a load variation is applied (from P= 0.2 pu to P = 0.6 pu).

-At t = 3s a power step is applied on VSC2 (P0 = 0.4 pu).

-At t = 5s, a line tripping is applied. Both converters are of identical rated power. The parameters of the controllers used are those in table II. From Fig. 13, it can be concluded that this control strategy is robust against events that may occur in a transmission grid (line opening, load variation, etc.). This paper presents an algorithm to determine the controllers parameters which consists in giving the pole of the system the best damping ratio using the linearized small signals model. The principle of the algorithm consists in varying the controllers parameters at each iteration by imposing a set of stability, dynamics and damping constraints on the system poles. This algorithm is applied to a gridforming voltage source converter with a cascaded multi-loop controller. The results obtained show that the proposed method provides a good dynamic behavior comparing to what has been achieved with the controllers tuning in the literature. More over, these set of parameters calculated for a given operating point and topology reveals a very good robustness in various topologies. The following works will be focused on the behavior in case of short circuit, which is another challenge to tackle for the grid forming converters. In addition, the interest in the Modular Multilevel Converters topology will increase, given its ability to generate an output AC voltage close to a sinusoid, therefore, no need to add an LCL filter which will greatly reduce dynamic interactions.

Figure 1 .

 1 Figure 1. Grid-forming Control

Figure 3 .

 3 Figure 3. Voltage control loop

Figure 4 .Figure 5 .

 45 Figure 4. Active power dynamics over load variation

Figure 6 .

 6 Figure 6. Active power dynamics over operating point change

State participation on λ 5- 6 λ 5 - 6 Figure 7 .

 6567 Figure 7. State participation in λ 5-6

Figure 8 .

 8 Figure 8. AC d-axis voltage and VSC output current

  : -Influence of the operating point. -Variation of Short Circuit Ratio. -Standalone mode. -Multi-converters grid topology. 1) Influence of operating point:

Figure 9 .

 9 Figure 9. Modes evolution in case of operating point change

-Figure 10 .

 10 Figure 10. Modes evolution in case of SCR variation.Left:All system poles. Right:Zoom on poles close to zero

Figure 11 .

 11 Figure 11. Left:Effect of the load variation on the active power. Right:Effect of step on d-axis AC Voltage

Figure 12 .

 12 Figure 12. Multi-converter topology

Figure 13 .

 13 Figure 13. Multi-converters topology

  ic and voltage controller proportional gain K pv . Using these new controllers parameters, The table

	TABLE IV
	SYSTEM EIGENVALUES
	λ 1-2 = -781.07± j3534	λ 3-4 = -754.64± j3011.8
	λ 5-6 = -36.65± j97.433	λ 7-8 = -12.98± j29.49
	λ 9 = -31.57	λ 10 = -2.3468
	λ 11 = -2.1055	λ 12 = -1.5605
	λ 13 = -1.6813
	controller integral gain K

ACKNOWLEDGMENT This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 691800. DISCLAIMER This paper reflects only the author's views and the European Commission is not responsible for any use that may be made of the information it contains.