Emmanuel Schermann

Hassan Omran

Sylvain Durand

Renaud Kiefer

Stochastic Trajectory Optimization for Autonomous Soaring of UAV

Keywords: Thermal soaring, speed to fly, stochastic optimization

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Stochastic Trajectory Optimization for Autonomous Soaring of UAV

Emmanuel Schermann * , * * Hassan Omran * , * * * Sylvain Durand * , * * Renaud Kiefer * , * *

INTRODUCTION

Large birds and glider pilots commonly use motions of air to extract energy from the atmosphere, enabling them to fly with minimum energy consumption. This process of using motions of the air mass to gain altitude is usually called "soaring". Soaring can be classified into i) dynamic soaring, where the source of the energy is a wind gradient Lawrance andSukkarieh (2009, 2011); [START_REF] Bird | Closing the Loop in Dynamic Soaring[END_REF], and ii) static soaring -the focus of the present studywhere rising air is the source of the lift. Buoyant plumes of warm or moist air called updrafts or thermals are generally used in static soaring. This rising air usually comes from near the ground where the thermal energy of the sun transfers to the air by conduction. The heated and/or moisten air then becomes lighter and rises. The greater the temperature difference between the warm air and the surrounding air, the higher the rate of ascent.

Rising air is a very complex and turbulent phenomenon FAA (2013). If the air close to the ground rises, there will be no gas to take its place. Thus when the air is homogeneously heated, the air does not rise. Nothing can happen until an instability ultimately causes the release of warm air (like non-homogeneous ground surface, nonhomogeneous heating, air turbulence, etc). Realistic dynamic convective environment simulation has been developed based on Rayleigh-Bénard convection model (1916).

The speed distribution of a thermal tends to be stronger in the core and reduces progressively outwards. Through a deeper analysis of the behavior of the air mass, several factors can contribute to increasing the difficulty of modeling such a phenomenon. Ascending currents are very hard to predict since they are the result of instability. They have a limited life span and drift with the wind. Moreover, they are transparent and there is currently no reliable way to acquire information on air mass movements, except by flying through them. Therefore, it seems extremely difficult to give a reliable model of this uncertain and nondeterministic parameter.

Many soaring strategies have been developed over the years by gliding pilots to improve the energy harvesting, allowing their plane to stay aloft for several hours without landing. On the other hand, autonomous UAVs (Unmanned Autonomous Vehicles) have become a subject of keen interest in different fields (such as surveillance, atmospheric measurements, fire monitoring, etc). Their main drawback is endurance because of the limited embedded energy resources. Thanks to the increasing computational power of on-board computers, autonomous soaring algorithms have been developed in recent years, allowing considerable improvements in endurance, by using natural energy resources of the atmosphere [START_REF] Lee | Predictive control for soaring of unpowered autonomous UAVs[END_REF]; [START_REF] Kagabo | Trajectory determination for energy efficient autonomous soaring[END_REF]; [START_REF] Andersson | Thermal Centering Control for Autonomous Soaring; Stability Analysis and Flight Test Results[END_REF].

The trade-off with thermal soaring is that UAV has to stop its cruising progression while gaining altitude (the time increases without changing the covered distance).

The present study is devoted to an optimization approach for thermal soaring, which is stochastic by nature because sources of ascending energy appear spontaneously and randomly. In a seminal work, [START_REF] Cochrane | MacCready theory with uncertain lift and limited altitude[END_REF] proposed a mathematical solution to the problem of static soaring, minimizing the time of flight for a given distance to go given a statistical model of the thermals. Here, the dual problem of maximizing the distance covered in a limited time is addressed. In order to solve the stochastic optimization problem despite the thermal modeling difficulties, the problem is separated into two sub-problems: i) decide whether to exploit a thermal when the glider encounters one (depending on its strength) or to move forward, and ii) find the optimal speed of flight between the different thermals (note that their position and strength are not a priori known). The uncertain parameter is then how probable it is to find a thermal with a given strength in the future.

The proposal could be used to provide a stochastic optimal trajectory, for instance in soaring competitions. Typically, low-scaled radio-controlled gliders operated by ground pilots contest in the so-called GPS-triangle competition.

The aim is to fly in a circuit of 2 kilometers, defined by 3 checkpoints, and do as many rounds as possible in 30 minutes without using any engine (see http:// gps-triangle.net/ for the rules and further details). The gliders have to efficiently acquire energy using static soaring to win the competition. Pilots could be interested in confronting their performance with an optimal procedure. Note that the proposed algorithm will not have any a priori information on the field (except a probabilistic model) compared to other pilots who can gather information by looking at various indicators (like other glider behaviors, birds, ground texture, ground exposure to sunlight, etc).

The rest of the document is organized as follows. The fly mechanics context and original Cochrane's theory are introduced in section 2. The proposed method is developed in section 3, with simulation results in section 4. Section 5 finally concludes the paper and gives some perspectives.

FLY MECHANICS

Sink rate and glide ratio

An aircraft (heavier than air) requires a series of aerodynamic forces to be able to fly without falling. During the flight, a glider can be considered as being acted by three forces. The force generated by the air is decomposed as the drag force, which is caused by the dissipating of the aircraft energy through air friction, and the lift which is by convention the rest of the aerodynamic forces. Thus, drag acts in the same direction as the air mass motion in the craft reference, while lift acts in the plane orthogonal to drag. Weight acts vertically at the center of mass.

As a result to the aerodynamic forces, and more particularly to the lift-to-drag ratio (one of the major goals in aircraft design is to maximize this ratio), when a glider moves forward in a steady flight (constant speed in still air), it will also sink. The relationship between the forward and the downward traveled distances, or equivalently speeds, is represented using the so-called polar curve (see Fig. 1), which can be found experimentally for a specific glider by measuring the sink at different speeds in calm conditions (FAA ( 2013)). The sink rate is how much altitude the glider loses over time and the glide ratio is how much distance it can travel over time. Two points are particularly interesting for soaring: the minimum sink speed where the plane maximizes the endurance and the best glide speed where the plane maximizes the forward-to-downward traveled distance, as depicted in Fig. 1. Actually, a model of the glider is not needed since the trajectory optimization will be based on the polar curve.

Fig. 1. Typical polar curve for a glider.

Speed to fly theory and MacCready setting

Speed to fly theory FAA (2013) gives the plane the ability to adjust its velocity in an optimal manner, taking into account the real-time measurement of the air motion and the expected updraft velocity of the next thermal that the plane will encounter. The aim is to maximize the average cross-country speed by optimizing the airspeed in both rising and sinking air. This theory is a well-known tool used by sailplane pilots. Speed to fly depends on

(1) The rate of climb expected to achieve the next thermal or updraft (2) The rate of ascent or descent of the air mass through which the glider is flying (3) The glider's inherent sink ratio as a function of glide ratio

The calculation of the optimal speed to fly to obtain the best glide ratio is the horizontal speed of the glider V that solves

V dS(V ) dV -S(V ) = V c = W + S min (1)
where V c is the average climb rate in the thermal (in m.s -1 ), W the thermal strength, S min the minimum sink the plane is capable of and S(V ) the sink rate of the glider. In other words, the optimal speed V depends on the average climb rate and should be continuously adjusted to obtain the best glide ratio. The mathematical development to obtain (1) is not detailed here for the sake of space constraints. This theory is very useful for the present study, because it gives the possibility to fly at the optimal speed given a sequence of checkpoints and thermals to follow, under the assumptions of knowing strength and position of the next thermals and uniform updraft characteristics over altitude. However, these assumptions are far from being realistic since i) neither the strength nor the shape are uniform with respect to altitude and ii) the position of thermals can only be approximately known (e.g. drift and lifespan) if a mapping is performed before the experiment.

The volatile and elusive character of thermals make any deterministic analysis of little use and, therefore, probabilistic studies should be more promising.

With experience, sailplane pilots can have an idea of the typical strength of a thermal and its particular at-mospheric conditions, but it is not guaranteed that the best thermals will be encountered while gliding. Therefore, they make use the speed to fly theory by using an expected thermal strength, i.e. V c in (1), which is called the MacCready setting. If the strength is higher than this setting, the plane should climb the thermal, otherwise the plane should continue cruising. When the plane is at high altitude, it is more likely that it will meet a strong thermal before reaching the ground, and the probability decreases as altitude decreases. The MacCready setting is eventually reduced when losing altitude, that might lead to a conservative flight, but it lowers the threshold for acceptable thermals in order to be able to gain altitude and avoid premature landing.

Cochrane statistical approach and extension

The MacCready setting ( 1) is an intuitive and qualitative approach, which has been confirmed by a quantitative statistical analysis in Cochrane (1999). Cochrane's strategy depends on a dynamic programming algorithm, in which the distance to fly (in the context of a standard US crosscountry contest) is divided into one-mile steps. The lifts the plane could encounter are based on a probabilistic characterization of thermals. For each step and any altitude, the optimal speed to fly is computed to minimize the expected time to complete the race, taking into account all possible thermals. The model is solved numerically from the last mile where the state of the plane is known (boundary condition is zero altitude) backward to the current position of the plane. At each step, the optimal setting is determined with respect to altitude for the given thermal condition.

US soaring contest (which Cochrane's approach was designed for) and GPS-triangle contest (application for the present proposal) are different in their objectives: 1) the aim in US soaring contest is to minimize the time required to cover a given distance, whereas 2) in GPS-triangle contest it is to maximize the covered distance in a given time.

While it could seem interesting to use the information harvested in the previous rounds, in reality thermals move and can dissipate within minutes (the glider takes approximately 2 minutes to make a turn in the GPS-triangle competition, without stopping at any thermal). It is then assumed here that no information is known concerning the location of thermals on the triangle trajectory. Therefore, to simplify the problem, the triangle can be unwrapped and considered as a straight line. This does not affect the result, considering that the turning time at a checkpoint is negligible compared to the whole flight.

TRAJECTORY STOCHASTIC OPTIMIZATION FOR AUTONOMOUS SOARING OF UAV

In the present study, the speed to fly as well as the choice to climb a thermal or cruise is decided by solving a stochastic optimization problem, using a dynamic programming resolution similar to Cochrane's work. This requires to find a set of state variables that contain all the information the autopilot of the UAV needs to compute the optimal MacCready setting (1). Cochrane used altitude and distance left as the state, but in the present case, the altitude h and the remaining time t are enough since the thermal probability model and all other parameters are not changing during the flight (the current traveled distance x has no impact on the optimal strategy). Therefore, the objective is to maximize the expectation of X, which is the total traveled distance in a given time interval max (E(X)) where E(•) denotes the expectation operator that applies over the random variation of future thermals.

Probability function and airlift/sink

It is supposed that for all possible thermal strength W the probability of finding at least one thermal of this strength over a given distance d is known. These probabilities can be found based on pilot experiences, or by making some statistics flying over the field before the experiment. Based on this, it is then possible to find a function f ∆x (W ) which corresponds to the probability of finding at least one thermal of strength W and no thermal of higher strength over a distance ∆x. Because the state is discretized over time, the traveled distance during one time step while flying at speed

v is ∆x = v∆t. So f ∆x (W ) is a density function with W f ∆x (W ) = 1.
It is assumed not to have horizontal wind, but it can be interesting to take into account the vertical air motion during cruising. Between the time instants t and t + ∆t, if the plane meets and crosses a thermal of core strength W , it will take some benefit from it. Otherwise, the air mass will only sink to compensate from the other thermals (mass conservation). The model of the airlift/sink during cruising is given by l(W ) = a + W b where a is the constant air sink and b is a coefficient with value between 0 and 1 (0 means no benefit from the thermal). The lift l(W ) is abbreviated as l in the sequel.

Value function

The value function Φ(•) denotes the total distance X the plane will achieve after a certain interval time, given the current time t and altitude h, but also the lift l that the plane will meet during the next time step. At a given state (t, h), if the already traveled distance x increases the value function will also increase and, therefore, it also varies with respect to x. The value function is finally defined as Φ(h, x, t, l) = max E(X|h, x, t, l) where the vertical bar | denotes a conditional expectation, i.e. this is the best E(X) the plane can do knowing h, x, t and l. It is the best guess over the performance the plane is able to achieve, knowing information on the plane position and the air it is flying into. In order to find the value function whatever the next lift l is, it is possible to find Φ(h, x, t) by writing

Φ(h, x, t) = max E(X|h, x, t) = E(Φ(h, x, t, l)|h, x, t) = W Φ(h, x, t, l(W ))f ∆x (W ) (2)
Note that the function value is never known in the following, even at the boundary conditions, but in fact its knowledge is not necessary to solve the problem.

Optimal speed to fly

Suppose that the choice of cruising has been made. Because Φ(h, x, t) needs to always be maximized from one step to another, the optimal speed to fly between the time t and t + ∆t can be found by picking the speed that maximizes the value function at the end of the next time step. The optimal speed v * is thus the one that gives

max v Φ( h, x, t)
where h, x, t are the values of the variables h, x, t at the next time step. For a plane flying at speed v, they are defined by

h = h + (l -S(v))∆t x = x + v∆t t = t + ∆t (3)
where S(v) is the sink rate at speed v, i.e. the polar. By setting the derivative of the value function Φ( h, x, t) with respect to v to zero (to be optimal) gives

Φh d h dv + Φ x dx dv + Φt d t dv = 0
where

Φ x := ∂Φ( h,x, t) ∂ x , Φh := ∂Φ( h,x, t) ∂ h
and Φt := ∂Φ( h,x, t) ∂ t

for ease of reading. Then, substituting (3) yields

Φ x Φh = dS(v) dv := S (v) (4)
Note that S (v) is the derivative of S(v) with respect to v.

If Φx

Φh is known, ( 4) is easy to solve because dS (v) dv is an affine function of v. Its solution is the optimal speed to fly v * that maximizes the value function for this time step. In order to find v * at all other steps, it is necessary to find a relation between Φ x/Φh := ∂Φ( h,x, t)

∂ x / ∂Φ(
where " " can be replaced by h, x or t. Substituting (2), then (5) turns out

∂Φ(h, x, t) ∂ = E ∂Φ(h, x, t, l) ∂ |h, x, t = W ∂Φ(h, x, t, l) ∂ f ∆x (W ) (6)
By combining (4), ( 5) and ( 6), then any optimal speed can be recursively computed in all possible states (h, t) if the boundary conditions are known.

Thermals and climb-or-cruise decision

So far, only the cruising case has been explored. The possibility of using thermals has then to be also included in the value function. When the plane stays in a thermal, its altitude increases, but the time spent also increases because the plane does not improve the traveled distance.

For a given plane position (h, x), a loss of time ∆t > 0 costs ∂Φ(h,x,t,l) ∂t ∆t < 0 in expected total distance. However, a gain of altitude ∆h > 0 improves the value by a distance ∂Φ(h,x,t,l) ∂h ∆h > 0. Therefore, climbing a thermal can be interesting only if the gain in expected distance due to ∆h is greater than the loss due to ∆t, that is if ∂Φ(h,x,t,l) ∂h ∆h > -∂Φ (h,x,t,l) ∂t ∆t. This implies

dh dt ≈ ∆h ∆t > - ∂Φ(h,x,t,l) ∂t ∂Φ(h,x,t,l) ∂h ≥ 0
where dh dt represents the vertical speed of the plane in the thermal. It can be defined by dh dt = W + S min as in (1), where S min is the minimum sink the plane is capable of. The criterion for the climb-or-cruise decision λ is

W + S min > - ∂Φ(h,x,t,l) ∂t ∂Φ(h,x,t,l) ∂h = - ∂Φ( h,x, t) ∂ t ∂Φ( h,x, t) ∂ h = - Φt Φh := λ (7)
If the decision of climbing is taken, the plane will stay in the thermal until the criterion ( 7) is not respected anymore. Moreover, since the thermal strength is known for the thermal, the knowledge of the problem during all the climb is the same. Then, the value function and its derivatives are equal to those of the point where the plane leaves the thermal.

Boundary conditions

The values of the derivatives Φt and Φh used in the climbor-cruise criterion (7) as well as Φ x used to find the optimal speed (4) are recursively defined and thus have to be initialized. The only states where the value function can be expressed explicitly is at the boundaries, namely when the plane has landed (h = 0) or when the time is over (t = T , with T the maximal flight time). Therefore, it is necessary to start the computing of the optimal strategy from these ending conditions, and iterate back in time to the starting time (t = 0).

First, when the plane is on the ground (h = 0), the traveled distance x will not evolve and thus the value function is equal to Φ(0, x, t) = x:

• If the plane was at a distance x + dx, the value function would be Φ(0, x + dx, t) = x + dx. It means that the derivative with respect to x is equal to Φ(0,x+dx,t)-Φ(0,x,t) dx = dΦ(0,x,t) dx = ∂Φ(0,x,t) ∂x = 1. • If the plane was at an altitude 0 + dh, the value function would be Φ(dh, x, t) = x+L/D max dh, where L/D max is the best glide ratio the plane is capable of (the plane fly as much conservatively as possible).

Thus ∂Φ(0,x,t) ∂h = L/D max . Note that the derivative at altitude h = 0 is not continuous, and the one that is needed is the one for h ≥ 0.

• If the plane was at a time t + dt, the value function would be Φ(0, x, t + dt) = x because time has no importance when the plane has landed: ∂Φ(0,x,t) ∂t = 0.

Second, in the case the plane is still aloft (t = T ), the value function is equal to Φ(h, x, T ) = x:

• If the plane was at a distance x + dx, the value function would be Φ(h, x + dx, T ) = x + dx and ∂Φ(h,x,T ) ∂x = 1 as before. • If the plane was at an altitude h + dh, the value function would be Φ(h + dh, x, T ) = x because the time is over. Therefore ∂Φ(h,x,T ) ∂h = 0.

But since the optimal speed (4) is calculated using Φx Φh , the optimal speed would be infinity for all altitude which is not a good place to start (during one time step back all values of the derivatives will correspond to a point with h = ∞). Since the plane has to be on the ground at the end, it is possible to find the derivatives at the time t = T -∆t , because at this place the plane has full knowledge of the future conditions. In this latter case (t = T -∆t), the value function is

Φ(h, x, T -∆t) = x + v(h, t)∆t = x + v(h, t)(T -t) (8)
If no lift is present at the last step, then the optimal speed v(h, t) can be found by solving

S(v(h, t)) = -h ∆t = -h T -t (9)
because -h ∆t is the vertical speed the plane should have to be at the ground at time T . Finally, the partial derivatives are (see Appendix 6.2)

∂Φ(h, x, T -∆t) ∂t = S(v(h, t) S (v(h, t)) -v(h, t) ∂Φ(h, x, T -∆t) ∂h = -1 S (v(h, t)) ∂Φ(h, x, T -∆t) ∂x = 1 (10)
It can be noticed that derivatives at the boundaries are not dependent on x, as initially assumed, so the problem can be solved using only the two state variables h and t.

Iterative resolution

Remind the present optimization problem is solved recursively from time t = T back to time t = 0. First, consider the plane is at the ground at time t = T (initialisation). All three derivatives Φ x, Φh and Φt are known from time T -∆t for all possible h, see (10). Then, for all possible thermal strength W , the recursion process is:

(1) The optimal speed at time t is found for all altitudes by solving (4).

(2) The derivative values ∂Φ(h, ,t,l) ∂ are reported at their corresponding altitudes based on (5) and (3). Because the plane increases its altitude while flying (when iterating from the end to the start), the values of the derivatives will not be known at altitude h = 0. The boundary values of the derivatives when h = 0 are used instead (see above).

(3) For all altitudes, the climb-or-cruise decision is taken based on (7). If the decision is to climb, the derivative values at the current position are replaced by those of the point where the thermal will be left (criterion not respected, or maximal altitude met). The leaving point can be few altitude steps and time steps away in the future, where the derivatives have already been calculated. (4) The expected values ∂Φ(h, ,t) ∂ are found using (6).

Then, go to the next time step t new = t old -∆t, and return to the first step of the recursion process. Do this until reaching the time t = 0.

SIMULATION RESULTS

The optimal procedure has been implemented in Matlab. The derivatives ∂Φ(h, ,t) ∂h and ∂Φ(h, ,t) ∂t are stored in order to be usable for fly tests in simulation or in real conditions (the derivative with respect to time is always equal to 1). Since the MacCready setting λ(h, t) := -Φt Φh allows to find the optimal speed from (4) and the climb-or-cruise criterion from ( 7), the analysis of the result can only focus on this value. Moreover, MacCready setting is easy to interpret since it represents how aggressive or conservative the plane has to fly. The only parameters that need to be defined in order to solve the optimization problem are the polar, the couple (a, b) to define the inter thermal lift/sink, the maximal and minimal altitudes where thermaling is allowed (the plane needs to be visible by the pilot, and altitude needs to be saved to land properly) and the thermal distribution. Note that the computing time when using 60 different thermal strengths is only 4 times the one when using only 6 values. Moreover, the computing time with steps that gives relatively smooth curves is less than 1s with a standard laptop and thus could be implemented on an in-board computer.

For the plane's flight simulation results, the thermal map is randomly constructed based on the thermal repartition. The plane flies exactly in the same conditions as in the optimization, but from the beginning (t = 0) to the end. Simulation tests over 10000 flights show that the results from the optimal strategy seem consistent and rather promising. The average speed is more than 11m.s -1 including stops for climbs and premature landing cases, and only 7% of premature landings occurred. A typical flight is given in Fig. 2(a). The thermal strengths denote the climbing speed in the thermal (W + S min ). The plane always leaves the thermal at an altitude where λ = W + S min . In a second example, see Fig. 2(b), the maximal possible climbing speed in the optimization was 3m.s -1 , and the last thermal the plane met was of this strength. Therefore, it is not interesting to stop for another thermal until the end. This can be observed by the fact that the MacCreading setting stays constant and does not decrease at the end. Furthermore, for all tests, the plane always land less than 2s before or after t = T , if no land-outs occurred. This is due to the discretization imperfections. It can also be noticed that the MacCreading setting tends to always be higher close to the end because the probability of using all the possible time increases with the time if the plane did not land before.

Note that these results cannot be directly compared to Cochrane's approach because the problem is different: maximize the total distance covered in a limited time vs. minimize the time of flight for a given distance.

CONCLUSIONS AND FUTURE WORK

This project investigated the possibility of using vertical movements of the air mass as an energy source for an autonomous UAV glider. The thematic being wide, this study was considered in the more restricted framework of the GPS-triangle competition for manual gliders, thus defining clear objectives and constraints. Research was carried out in order to efficiently choose the sources of energy to be used, as well as the speed of flight to adopt. The problem has been simplified by different assumptions, in order to allow the optimal resolution of this optimization problem taking into account the uncertain nature of the occurrence of these sources. The procedure for finding the optimal answer, based on the concept of dynamic programming, was implemented under Matlab and tested on simulated glider flights.

Next step will be the implementation of the method in a low-scale glider for real tests. Different improvement areas will also be studied, like i) taking the centering time into account when the plane meets a thermal, ii) consider the wind that directly affects thermals and therefore the optimal strategy, or iii) build the thermal probability repartition P d (W ) in real time (instead of assuming it is known in advance).

APPENDICES

How to find (5)

If the decision to cruise is made, the knowledge on the problem is the same at the point (h, x, t) knowing the lift l, as at point ( h, x, t), and thus it is possible to write the value function iteration equation Φ(h, x, t, l) = max E(X|h, x, t, l) = max v Φ( h, x, t)

If the plane flies at the optimal speed v * to the next time step in the lift l, previous equation becomes Φ(h, x, t, l) = Φ( h, x, t) = Φ(h + (l -S(v The extension for the derivative with respect to x and t is trivial, that finally yields the general relation (5).

.

  (a) Trajectory with respect to the traveled distance of a typical flight simulation. A lap approximately corresponds to 2km. . (b) Other flight trajectory.

Fig. 2 .

 2 Fig. 2. Illustrative examples of flight tests in simulation.

6. 2

 2 How to find(10) On the first hand, the partial derivative of (8) with respect to time is∂Φ(h, x, T -∆t) ∂t = ∂v(h, t) ∂t (T -t) -v(h, t)Furthermore, the time derivative of S gives ∂S(v(h, t(h, t))(T -t) 2 Combining the latter relations yields the first relation in (10).On the other hand, the partial derivative with respect to altitude is∂Φ(h, x, T -∆t) ∂h = ∂v(h, t) ∂h (T -t)and the derivative of S with respect to h gives ∂S(, yields the second relation in (10).

  Thus the previous relation can be written as ∂Φ(h, x, t, l) ∂h = Φh 1 -Φ x Φh

									∂v ∂h	∆t + Φ x ∂v ∂h	∆t = Φh
	Its derivative with respect to h gives
	∂Φ(h, x, t, l) ∂h	=		∂Φ( h, x, t) ∂h	= Φh	∂ ∂h h	+ Φ x ∂ ∂h x	+ Φt	∂ ∂h t
	Using (3) and substituting (4) gives
	∂ ∂h h	= 1 -= 1 -	∂S(v) ∂v Φ x Φh ∂v ∂h	∂v ∂h ∆t	∆t = 1 -S (v)	∂v ∂h	∆t
	∂ ∂h x ∂ t ∂h	= = 0 ∂v ∂h	∆t

* ))∆t, x + ∆t v * , t + ∆t)


Salim El Ghouli
email: salim.elghouli@st.com

Jean Michel Sallese

André Juge

Patrick Scheer

Christophe Lallement

Transadmittance efficiency in presence of NQS transport in Asymmetric Double Gate FDSOI MOSFET

Keywords: NQS, HF, Transadmittance, Transconductance efficiency, gm over Id, Double-gate FETs, FDSOI, UTBB, Analog and RF, Low-Power, Low-Voltage I

The state of the art RF and millimeter-wave circuits design requires accurate prediction of the non-quasistatic (NQS) effects at high frequency for all levels of channel inversion. This work provides a practical insight to help high frequency performance assessment of UTBB FDSOI MOSFETs through a powerful frequency normalization scheme. Frequency dependence of small signal characteristics derived from experimental S-parameters are analyzed and reveal that the transconductance efficiency (gm/ID) concept, already adopted as a low frequency Analog figure-of-merit (FoM), can be generalized to high frequency, including under asymmetric operation. We report that the normalized frequency dependence of the generalized transadmittance efficiency (ym/ID) FoM only depends on the mobility and inversion coefficient (IC). In addition, this approach is also used to extract essential parameters such as the critical NQS frequency fNQS.

INTRODUCTION

RF CMOS based MOSFETs have been widely used in RF integrated circuits (RFIC) thanks to the geometry downscaling. Critical issues related to short-channel effects have been mitigated using innovative architectures coupled with new materials [START_REF] Raskin | FinFET versus UTBB SOI -a RF perspective[END_REF]. The new architectures allow for higher level of integration and higher frequency operation of the CMOS technology. For analog and RF applications, the ultra-thin Body and Box (UTBB) fully depleted (FD) Silicon-on-Insulator (SOI) technology has proven to be relevant, especially when power consumption matters [START_REF] Raskin | FinFET versus UTBB SOI -a RF perspective[END_REF]. UTBB FDSOI MOSFETs exhibit high analog and RF performances thanks to the reduced parasitic capacitances and resistances. Moreover, the independent back gate allows threshold voltage and channel mobility controls [START_REF] Noel | Multi-VT UTBB FDSOI Device Architectures for Low-Power CMOS Circuit[END_REF] [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF].

In this work, the UTBB FDSOI MOSFET transconductance efficiency (gm/ID) is generalized and studied at high frequency for different channel inversion levels. The main attention is paid to the impact of Non-Quasi-Static (NQS) effect on the proposed transadmittance efficiency Figure-of-Merit (FoM), and in particular how a proper choice of the normalization factors ends with simple representation of this fundamental figure of merit, with special focus on low-moderate inversion levels. The goal is to provide RF designers with a practical tool for high frequency performance assessment. The paper is organized as follows. The device architecture along with the high frequency characterization are described in Section II. Next in section III we discuss and validate experimentally the generalization of the transconductance efficiency concept in UTBB FDSOI MOSFETs derived at low frequency in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF]. The distributed gate effects are assessed using high frequency measurements of large gate devices in Section IV, and an onset frequency of the channel NQS effect is defined and compared with transit frequency in Section V. The back gate transadmittance efficiency is also analyzed using TCAD simulations in Section VI. Finally, discussion and conclusion are presented in Section VII and VIII, respectively.

II. DEVICE DESCRIPTION AND CHARACTERIZATION

Measured UTBB FDSOI MOSFETs were processed at STMicroelectronics. Si film, BOX and equivalent gate oxide thicknesses are 7 nm, 25 nm and 1.3 nm, respectively. More details on the process can be found in [START_REF] Planes | 28nm FDSOI technology platform for high-speed lowvoltage digital applications[END_REF]. DC and CV measurements are carried out to extract the device electrical parameters (e.g. threshold voltage VTH, oxide capacitance Cox, etc.). In particular, the Inversion Coefficient (IC) is derived following the method in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF]. Multi-finger MOSFETs are embedded into CPW (Coplanar-Waveguide) transmission lines for high frequency characterization. After a standard short, open, matched load, and through (SOLT) calibration for reference plan definition, 2-port S-parameters are measured using ground-signal-ground (GSG) configuration probes up to a frequency of 110 GHz, under several DC bias conditions. All modes of operation are covered: non-saturation, saturation, weak inversion (WI), moderate inversion (MI) and strong inversion (SI). Measurements are de-embedded using Open-Short method. The gate lengths of the measured devices range from 28 nm to 1 µm, while gate widths range from 200 nm to 5 µm. Both N and P MOSFET types are measured and analyzed in this work. Comparable results are obtained for both types with a more pronounced NQS effect in the case of P type due to a lower carrier mobility. However, for reasons of limited space, only NMOS results are shown here.

III. FREQUENCY DEPENDENCE OF THE TRANSADMITTANCE EFFICIENCY IN LONG CHANNEL UTBB FDSOI MOSFETS

In [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF], the invariance of the transconductance efficiency, a well-known feature of bulk MOSFETs, was shown to remain valid for UTBB FDSOI MOSFETs that are controlled by two independent gates. This finding was not expected and constituted a major result for low frequency design methodologies. The detrimental impact of the high inversion level on the transconductance efficiency is further amplified under NQS operation, as shown in reporting the measured transadmittance efficiency modulus versus IC and frequency (F) in Fig. 1. Here, IC is the inversion coefficient defined as [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF]:

𝐼𝐶 = 𝐼 𝑊 𝐿 • 𝐼 (2)
where ID, L, W and Ispec are respectively the DC drain current, the MOSFET length, width, and the 28 nm FDSOI technology current estimated for the front gate at Ispec = I1 = 0.7028 µA when assuming W = L [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF].

Various models of NQS effect have been reported for bulk [START_REF] Wang | Quasi-static and nonquasi-static compact mosfet models based on symmetric linearization of the bulk and inversion charges[END_REF][6] and for common double gate MOSFETs [7][8]. A practical large signal model has also been published in [START_REF] Scholten | A large signal non-quasi-static MOS model for RF circuit simulation[END_REF]. However, none of these models proposes an exact analytical solution of the NQS behavior. To the best of the author's knowledge, the first analytical model of NQS in bulk MOSFETs valid in all regions of operation was derived in [START_REF] Sallese | A novel approach to charge-based nonquasi-static model of the MOS transistor valid in all modes of operation[END_REF] which gave rise to a dimensionless powerful abacus where normalized transadmittances are represented in terms of normalized frequencies [START_REF] Porret | A compact non-quasi-static extension of a charge-based MOS model[END_REF].

In this respect, we propose to investigate the front and back gates high frequency NQS behavior, and verify if this invariant still holds for UTBB FDSOI MOSFETs. As for the gm/ID figure of merit [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF], we define the transadmittance efficiency as:

𝑇𝐸 = 𝑦 𝐼 = 𝑦 -𝑦 𝐼 ( 1
)
where ymi is the mutual transadmittance defined as the difference between the gate-to-drain and drain-to-gate transadmittances, i.e. ydg and ygd respectively. The transadmittance efficiency is a generalization of the transconductance efficiency defined at low frequency operation.

At high frequency, the front or back gate transconductance is an imaginary number with a modulus and phase. Since gmi represents the mutual interaction between the gate and drain, ymi keeps the same meaning even for high frequency (note that ymi reverts to gmi at low frequency). It should be noted that the capacitive coupling attributed to Cgd, which dominates yDG at high frequency, cancels in the definition of ymi, as yGD ≈ -jωCgd [START_REF] Tsividis | Operation and modeling of the MOS transistor[END_REF] in the low frequency limit. Then, gmi and ymi assess the gate control of the intrinsic MOSFET channel. At low frequency (quasi static operation) we have [START_REF] Tsividis | Operation and modeling of the MOS transistor[END_REF]:

𝑔 = 𝑟𝑒𝑎𝑙(𝑦 ) = |𝑦 | (3)
At any frequency, the ymi/ID FoM determines the gain and phase shift for a given DC current ID. The evolution of the transadmittance efficiency versus frequency provides a very valuable insight in terms of power performances and efficiency at high frequency in presence of NQS effects.

At low frequency, it is well known that the transconductance efficiency is the highest in weak inversion as evidenced in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF]. However, at high frequency as shown in Fig. 1, the NQS transport degrades the low frequency efficiency, in particular for weak and moderate inversion (IC < 10). This is also illustrated in Fig. 2, a side view of Fig. 1. The higher efficiency of the moderate and weak inversion modes (i.e. gm1/ID > 10) does not hold anymore beyond F = 4 GHz for L = 1 µm. Moreover, beyond 13 GHz and unlike in low frequency operation, strong inversion performs better with regards to the transadmittance efficiency. At high frequency, the lower cutoff frequency (lower current) in weak inversion results in a lower transconductance efficiency. Indeed, at high frequency, the transit time is large below the threshold, and the current lags behind the voltage excitation. Conversely, in SI, the transit time is small and, even though the DC gm/ID is lower, combining these effects give a transconductance efficiency that is greater than in weak inversion. Plotting the integral of each curve in Fig. 2 over the interval [0, FNQS] versus IC where FNQS is an onset frequency of NQS effect defined in (19) or (20), the maximum is found in SI (IC > 10) for L = 1 µm, while it is located in MI for L = 300 nm. This integral is a new FoM that we propose for the tradeoff between performance, power consumption, and speed. This FoM is more realistic than the previously defined fT.gm/ID in [START_REF] Shameli | Ultra-low power RFIC design using moderately inverted MOSFETs: an analytical/experimental study[END_REF] since it does not assume that the transconductance vanishes beyond fT and does not underestimate the performance in weak inversion. Indeed, as it will be discussed later, fT is not the onset of NQS operation in weak inversion. The real high frequency capability is better captured by a NQS onset frequency FNQS that will be defined in Section V. Moreover, in contradiction to the previously reported conclusions such as in [START_REF] Shameli | Ultra-low power RFIC design using moderately inverted MOSFETs: an analytical/experimental study[END_REF], the optimum is now located in SI for long channels, and no longer in moderate inversion. For shorter devices, the optimum is still located in MI, as already reported using the fT.gm/ID FoM. Fig. 3 shows the phase of the transadmittance efficiency versus frequency for various IC values. According to the approximation of the transadmittance efficiency expression in [START_REF] Scholten | A large signal non-quasi-static MOS model for RF circuit simulation[END_REF], in regard to the modulus, the phase is more sensitive to NQS effects and particularly in weak inversion. For instance, at lower frequency (F < 500 MHz) and for IC = 0.3, the phase delay is more than 31° while the modulus remains almost unaffected, i.e. 97 % of its DC value. Fig. 2 and Fig. 3 reveal that the effect of NQS transport on the transadmittance efficiency is gradual and depends on the channel inversion coefficient. Having a clear definition for the onset frequency of NQS operation is then needed to provide a more realistic insight into the high frequency performance of MOSFETs.

𝐹𝑂𝑀 = |𝑦 | 𝐼 𝑑𝐹 (4)
Exact small signal analytical expressions for NQS operation are presented in [START_REF] Sallese | A novel approach to charge-based nonquasi-static model of the MOS transistor valid in all modes of operation[END_REF] for single gate MOSFETs and simplified in [START_REF] Porret | A compact non-quasi-static extension of a charge-based MOS model[END_REF] to derive a compact model. The model reported in [START_REF] Porret | A compact non-quasi-static extension of a charge-based MOS model[END_REF] proposed a frequency normalization scheme using the first pole obtained from the first order approximation of the NQS model of bulk MOSFETs. The critical frequency Fcrit expression, which corresponds to the first pole of the frequency characteristic of the device, as a function of the inversion coefficient IC concept redefined in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF] for a double gate (DG) UTBB FDSOI MOSFET is:

F = 1 2𝜋 15 4 • 𝜇𝑈 𝐿 • 2 + 6 𝐼𝐶 + 3 𝐼𝐶 + 1 4 + 4 𝐼𝐶 + 1 4 1 + 2 𝐼𝐶 + 3 𝐼𝐶 + 1 4
(5)

≈ 15 4𝜋 × 𝜇𝑈 𝐿 • √𝐼𝐶 (in SI) (6)
where µ and UT are respectively the effective mobility and the thermal voltage. ( 6) is an approximation of (5) while neglecting smaller terms in strong inversion, with less than 1 % error on Fcrit for IC > 2. The effective mobility is taking into account the impact of the vertical electrical field on scattering. The velocity saturation effect is found to be negligible. Indeed, according to [START_REF] Roy | Compact Modeling of Anomalous High-Frequency Behavior of MOSFET's Small-Signal NQS Parameters in Presence of Velocity Saturation[END_REF], an anomalous frequency dependence of the gate-to-drain transadmittance yDG in saturation is expected. This would manifest by an increase of yDG with frequency, a trend which is not observed in our case.

It should be noted that Fcrit is proportional to the effective mobility and UT, inversely proportional to the square of gate length, and importantly, it accounts for the dependence of NQS on the inversion coefficient IC. This critical NQS frequency will be used for frequency normalization hereafter.

The transadmittance efficiency in saturation is derived from [START_REF] Porret | A compact non-quasi-static extension of a charge-based MOS model[END_REF] as :

𝑦 𝐼 = 𝑔 𝐼 𝜉 𝐹 𝐹 , 𝐼𝐶 (7)
where gm1 is the low frequency front gate transconductance, ID is the DC drain current, F is the operation frequency, Fcrit is given by relation ( 5), and ξm is a NQS small-signal function based on the first kind Bessel functions. Defining a critical time delay as:

𝜏 = 1 2𝜋 • 𝐹 (8)
An approximation of (7) using the first pole of ξm becomes:

𝑦 𝐼 = 𝑔 𝐼 1 (1 + 𝑗𝜔𝜏 ) (9)
where ω is the angular frequency (ω =2πF).

In Fig. 4, a normalized modulus of the transadmittance efficiency (|ym1|/ID)norm is shown for various values of IC versus the normalized frequency Fnorm1 defined as:

𝐹 = 𝐹 𝐹 (10)
In addition, a further normalization of |ym1|/ID is applied using the low frequency transconductance efficiency gm1/ID as follows:

|𝑦 | 𝐼 = |𝑦 | 𝐼 /(𝑔 /𝐼 ) (11)
All normalized moduli plotted in Fig. 4 from weak to strong inversion are almost superposed as a result of the frequency normalization process in [START_REF] Sallese | A novel approach to charge-based nonquasi-static model of the MOS transistor valid in all modes of operation[END_REF] (despite the very high frequency disturbance of the parasitic elements). The invariance is also evidenced for the phase in Fig. 5 up to a phase shift as large as 270°. It should be noted that a phase shift of 45° is reached at Fnorm1 = 1.

These results are of prime importance for modeling NQS effects in UTBB SOI FETs. For the first time, a simple but powerful analytical approach can predict the gate transconductance, modulus and phase, during highly non static operation of the intrinsic channel and in all the regions of operation. Additional normalizations in terms of the inversion coefficient and critical frequency were used to derive a kind of abacus, paving the way to a simple and efficient HF design methodology where device parameters, biasing conditions and frequency are treated on the same ground.

IV. IMPACT OF THE INTRINSIC GATE DEALY

In this Section, we investigate the breakdown of the normalized transadmittance efficiency invariance detailed in Section III. In particular, we will study the impact of the distributed effects along the gate on the critical frequency Fcrit used for normalization.

A closer look at the normalized curves depicted in Fig. 5 reveals that in SI regime the transadmittance efficiency decreases earlier than in WI while the frequency increases. This enhanced phase shift in SI can be explained either by a decrease in the mobility or by an additional delay caused by the gate distributed effect. Note that in our previous normalizations (Fig. 4 and Fig. 5), the gate delay generated by the 1 µm finger length (i.e. gate width) was ignored.

At high frequency, the gate can be considered as a transmission line because of the limited electromagnetic wave velocity along the gate finger length. An estimation of the time constant of the distributed gate finger is given by [START_REF] Hasani | Analysis and modeling of a MOS transistor with long gate finger for millimeter wave band applications[END_REF]:

𝜏 = 𝑅 𝐶 (12)
where Rgg and Cgg are gate resistance and average gate capacitance respectively. It should be noted that Cgg depends on the level of inversion, which could explain the lower critical frequency in SI case when gate distributed effect is important. Conversely, in weak inversion Cgg is lower and a higher critical frequency is expected. According to [START_REF] Razavi | Impact of distributed gate resistance on the performance of MOS devices[END_REF][17], we have:

𝑅 = 𝑅 𝜂 ( 13
)
where η is equal to 3 for a single-side connected gate and is equal to 12 for a double-side connected gate. Rg is the resistance of the MOSFET finger at low frequency (i.e. without the distributed effect) and extracted in SI. If we consider the additional pole due to the gate distributed effect, we can derive an equivalent global time constant as:

𝜏 = 𝜏 + 𝜏 (14)
Considering that delays add, the delay caused by both gate finger length (i.e. the width of the transistor) and NQS are added using the semi-empirical formula [START_REF] Roy | Compact Modeling of Anomalous High-Frequency Behavior of MOSFET's Small-Signal NQS Parameters in Presence of Velocity Saturation[END_REF]. Consequently, the time constant in ( 14) includes distributed effects along the channel (i.e. NQS) and through the gate. Thus, the new critical frequency is:

𝐹 = 1 2𝜋 • 𝜏 (15)
Finally, neglecting higher order terms in the denominator and using the same approach, the transadmittance efficiency expression given in (9) becomes:

𝑦 𝐼 ≈ 𝑔 𝐼 1 1 + 𝑗𝜔 • 𝜏 (16)
A new normalized frequency considering non static transport in the channel and along the gate is now given by:

𝐹 = 𝐹 𝐹 (17)
For L = 1 µm, the error on the frequency predicting a 45° phase shift while using the normalization Fnorm2 (not shown here) is less than 15 % (5 % enhancement versus Fig. 5) since gate distributed effects are negligible in regard to NQS effect. According to [START_REF] Wang | Quasi-static and nonquasi-static compact mosfet models based on symmetric linearization of the bulk and inversion charges[END_REF], this residual discrepancy could be explained by a mobility varying with IC, which would be estimated at about 15 % between WI and SI for a long channel (i.e. L = 1 µm). Conversely, the error on the phase shift prediction for a given frequency is less than 6 %.

For L = 300 nm, the normalization using either Fcrit or Fdis leads to similar conclusions as for the L = 1 µm case. Fig. 6 shows the phase shift for L = 300 nm MOSFET versus frequency. Based on the apparent critical frequency, the dominant distributed effect is the channel NQS. In Fig. 7, the frequency is normalized using Fcrit to account for the channel NQS effect only. Clearly, the gate delay is evidenced since a higher inversion coefficient is responsible for 'earlier' phase shift. This cannot be attributed to the channel where higher IC moves the phase shift to higher frequencies. In Fig. 8, the second normalization is introduced using Fdis instead of Fcrit to also account for the distributed gate delay. The normalization using Fdis is taking into account both delays: (1) the one caused by the distributed effect along the gate finger, and (2) the other caused by the NQS effect along the channel. Therefore, Fdis is representing same critical frequency for each value of IC. Normalization using Fdis in Fig. 8 wipes out the differences generated by IC in the critical frequencies through both the gate delay and NQS, which finally results in merging all the curves in a single one. As a result, the discrepancy of the frequency for a 45° phase shift is reduced to less than 14 % as shown in Fig. 9, and the error on the predicted phase is less than 6 %. Again, the remaining mismatch could be explained by a variation of the mobility with the gate voltage. This point will be investigated in future work. For shorter devices in RF circuits such as high gain amplifiers, the gate finger length is usually long (i.e. increased width for higher current drive or larger transconductance). In this case, the distributed effect along the gate finger adds an important delay that dominates the NQS effect in the shorter channel. This is evidenced on the phase delay versus the real frequency for the 28 nm channel length displayed in Fig 10 . In the short device, the inversion level (IC) seems to have minimal impact on the critical frequency, see Fig. 10, and the two counter effects almost compensate each other. The gate distributed delay slightly take over since the phase shift happens earlier when IC is increased, which should not be if the channel would still be the bottleneck in terms of delay. To follow up with the normalization scheme, we calculate the 'intrinsic channel normalization frequency' Fnorm1 from relation [START_REF] Sallese | A novel approach to charge-based nonquasi-static model of the MOS transistor valid in all modes of operation[END_REF] and plot the phase shift accordingly in Fig. 11. Whereas for 'long' channels, this first normalization almost condenses the bunch of curves in a single characteristics (see Fig. 5 and Fig. 7), here it spreads out. The gate distributed effect could explain this dispersion. Fig. 12 shows the phase shift for the 28 nm channel length device when using the normalization frequency Fnorm2 as given from ( 17) which now includes the first order non quasi static effect originating from the gate. The way the different characteristics are distributed in a narrow region is in favor of the predominance of the enhanced gate delay in short channel devices. Moreover, thanks to the mitigated channel NQS effect, the phase shift does not reach 45° within the frequency range from 100 MHz to 110 GHz for all levels of inversion. V. ONSET OF NQS EFFECTS Significant NQS effects pop up when the operation frequency reaches the critical frequency Fcrit defined in [START_REF] Wang | Quasi-static and nonquasi-static compact mosfet models based on symmetric linearization of the bulk and inversion charges[END_REF]. In particular, a phase shift greater than 45° is measured at this characteristic frequency, meaning that Fcrit overestimates the threshold of the NQS effect. To our knowledge, no clear definition is given for the emergence of NQS transport. Usually, the concept of transit frequency fT is used, i.e. the frequency where the modulus of the current gain HDG falls to unity (i.e. 0 dB), and the NQS frequency is often defined as being proportional to fT [START_REF] Tsividis | Operation and modeling of the MOS transistor[END_REF]. However, this definition relies on the low frequency transconductance gm, and therefore it is implicitly a quasi-static definition. Noting that the current gain HDG is defined for a 2port structure as:

𝐻 = 𝐼 𝐼 = 𝑦 𝑦 (18)
where ID and IG are the current phasors across port-2 (drain) and port-1 (front gate), respectively. VD is the voltage across port-2, which is impedance-matched in this case. In order to assess the relation between fT and the NQS critical frequency, the current gain is plot for several IC values with respect to the normalized Fnorm1 frequency. Fig. 13 shows the modulus of the current gain HDG as a function of the normalized frequency Fnorm1 for various IC and for a channel length of L = 1 µm. fT is proportional to Fcrit in MI and SI as can be expected from ( 6) and ( 2), and noting that fT is proportional to the square root of the drain current provided velocity saturation is ignored. However, for WI, the proportionality becomes less valid, as expected given that Fcrit becomes constant in WI while fT remains proportional to IC. If we use the SI approximation [START_REF] Zhu | Accurate RTA-Based Nonquasi-Static MOSFET Model for RF and Mixed-Signal Simulations[END_REF] for the frequency normalization for all the levels of inversion with L = 1 µm, all |HDG| plots are superimposed independently of IC as shown in inset of Fig. 13. This suggests that the transit frequency is equal to one third of the critical frequency in SI and that it remains proportional to the square root of the current all the way, from WI to SI.

Therefore, the proportionality between the onset of NQS frequency and fT holds only in SI for long channels. In WI and practically in MI, this rules of thumb (proportionality to fT) underestimates the critical NQS frequency. Depending on the criteria invoked, we propose to use two expressions for the onset of non-static transport:  10° phase shift of the transadmittance efficiency:

𝐹 °= 1 2𝜋 tan 𝜋 18 1 2𝜋 • 𝐹 + 𝑅 𝜂 • 𝐶 (19)
 3% decrease in the transadmittance efficiency modulus from the low frequency value:

𝐹 % = 1 2𝜋 0.25 1 2𝜋 • 𝐹 + 𝑅 𝜂 • 𝐶 (20)
Note that FNQS3% is practically equal to fT in SI and, higher in WI and MI.

VI. BACK GATE CONTROL

In a DG UTBB MOSFET, the independent back gate voltage controls the threshold voltage and consequently IC. The back gate also influences the apparent mobility of the carriers in the silicon film thanks to the volume inversion in low electrical field condition. Therefore, as the NQS characteristic frequency Fcrit depends on mobility and IC, the back gate should affect NQS effects. In this Section, the impact of the back gate voltage on |ym2|/ID for a fixed front gate voltage is assessed using TCAD simulations. Note that due to lack of dedicated back gate RF structures, TCAD simulations are used to validate the invariant representation of transadmittances when the RF signal is applied to the back gate. Using same normalization as in [START_REF] Jindal | Noise associated with distributed resistance of MOSFET gate structures in integrated circuits[END_REF], we report that the different curves merge as in Fig. 4 and Fig. 5 for the front gate transadmittance efficiency |ym1|/ID versus frequency while setting IC with the back gate voltage (not presented here). Interestingly, when the input signal is applied to the back gate, a back gate transadmittance efficiency |ym2|/ID can also be defined.

In order to evaluate the back gate performance at high frequency, as in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF], the transadmittance efficiency of the back gate is now assessed for a fixed front gate voltage VGS = 0, and for various inversion levels imposed by the back gate (various VbG). Fig. 14 shows the phase of the transadmittance efficiency of the back gate at VGS = 0 V. In this case, the channel is spread in the volume of the silicon film. Again, after applying the same normalization for the frequency and for the transadmittance to current ratio, a unique plot is obtained in Fig. 15. This supports that the normalization scheme used to predict non static transport in UTBB FDSOI MOSFETs is valid whatever gate is used.

VII. DISCUSSION

The developments leading to the definition of a NQS critical frequency Fcrit for a single gate bulk MOSFET developed in [START_REF] Sallese | A novel approach to charge-based nonquasi-static model of the MOS transistor valid in all modes of operation[END_REF] [START_REF] Porret | A compact non-quasi-static extension of a charge-based MOS model[END_REF] hold for a DG MOSFET, including with asymmetric operation. In particular, the normalization of the frequency leads to almost unique charts for long channel devices allowing fast and accurate high frequency performance assessment. Moreover, the invariance is also confirmed for short channel devices with a dominant gate distributed effects provided that the related time constant is taken into consideration. The transadmittance efficiency concept defined in this paper provides a deep insight into high frequency UTBB FDSOI MOSFET behavior and gives an estimation of the limiting MOSFET characteristic frequencies. The critical channel NQS frequency depends on the device length, free carrier mobility and inversion coefficient IC. The gate distributed effect, in turn, depends on the gate resistance and capacitance.

The frequency normalization using the distributed gate time constant generalizes the normalization scheme by merging the longitudinal (gate) distributed effect with the channel NQS effect. The invariant normalized transadmittance efficiency versus normalized frequency was confirmed extensively with both measured and simulated data. This novel approach provides a way to calculate the critical frequency beyond which the gain and the phase will be affected.

In Analog and RF design, oversimplified rules of thumb are commonly used, and bias the analysis. For instance, stating that operating frequency should be lower than fT/10 or fT/3 means that the designer will choose the shortest MOSFETs, which introduce detrimental mismatch and SCE, especially Drain Induced Barrier Lowering (DIBL). Such a rule of thumb can be overlooked using the transadmittance efficiency charts and critical frequency concepts developed in this paper based on the Inversion Coefficient exposed in [START_REF] Ghouli | Experimental gm/ID Invariance Assessment for Asymmetric Double-Gate FDSOI MOSFET[END_REF] since knowing Fcrit versus the channel length and DC operating point, the optimum geometry can be selected unambiguously.

VIII. CONCLUSION

In this work, we have shown that we can transpose the powerful normalizations derived for high frequency operation of bulk MOSFETs to the independent double gate FET in UTBB SOI technology. This ends up with very simple expressions for the phase and the magnitude for the front and back gate transadmittances versus frequency. In addition to the intrinsic channel delay, we have introduced the effect of the distributed gate thus taking into account lateral and longitudinal distributed effects in a single expression. This offers a practical insight to help high frequency performance assessments of UTBB FDSOI MOSFETs. Based on experimental results, we have shown that frequencies of interest where distributed effects occur can be predicted using the simple analytical expressions detailed in this work down to 28 nm channel length and up to 110 GHz. The normalization scheme is applicable to whatever the gates length and DC biasing, and supports the very general character of the NQS analysis proposed in this work.

Fig. 1 3

 1 Fig. 1 3-D plot of the measured transadmittance efficiency modulus over ID versus the frequency and inversion coefficient IC for L = 1 µm at VDS = 1 V.

Fig. 2

 2 Fig. 2 Measured transadmittance efficiency modulus versus the frequency for different IC for a NMOS (L = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V.

Fig. 3

 3 Fig. 3 Measured phase of the transadmittance over ID versus frequency at different IC for NMOS (L = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V.

7 Fig. 4

 74 Fig. 4 Normalized transadmittance modulus over ID versus normalized frequency Fnorm1 at different IC for NMOS (L = 1 µm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V (Measurements).

Fig. 5

 5 Fig. 5 Phase of the transadmittance over ID versus normalized frequency Fnorm1 at different IC for NMOS (L = 1 µm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V (Measurements).

7 Fig. 6

 76 Fig. 6 Phase of the transadmittance over ID versus frequency at different IC for NMOS (L = 300 nm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V (Measurements).

Fig. 7

 7 Fig. 7 Phase of the transadmittance over ID versus normalized frequency Fnorm1 with constant mobility at different IC for NMOS (L = 300 nm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V (Measurements).

Fig. 8

 8 Fig. 8 Phase of the transadmittance over ID versus normalized frequency Fnorm2 with constant mobility at different IC for NMOS (L = 300 nm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V (Measurements).

Fig. 9

 9 Fig.9Frequency error versus phase shift after normalization using Fdis (IC = 93.1 is taken as a reference)

Fig. 10

 10 Fig. 10 Measured phase of the transadmittance over ID versus frequency at different IC for NMOS (L = 28 nm & W = 1 µm) in saturation (VDS = 1 V) and back gate voltage VbG = 0 V.

Fig. 11

 11 Fig. 11 Phase of the transadmittance over ID versus the normalized frequency Fnorm1 estimated with a constant mobility at different IC for NMOS (L = 28 nm & W = 1 µm) in saturation (VDS = 1 V) and VbG = 0 V (Measurements).

Fig. 12

 12 Fig. 12 Phase of the transadmittance over ID versus normalized frequency Fnorm2 with constant mobility at different IC for NMOS (L = 28 nm & W = 1 µm) in saturation (VDS = 1 V) and VbG = 0 V (Measurements).

Fig. 13 Measured

 13 Fig. 13 Measured |HDG| versus normalized frequency Fnorm1 at different IC for NMOS (L = 1 µm & W = 1 µm), (VDS = 1 V) and back gate voltage VbG = 0 V. Inset shows |H21| versus F/FcritSI where FcritSI is given by (6).

Fig. 14

 14 Fig. 14 Simulated phase of the back transadmittance efficiency versus frequency at different VbG for NMOS (L = 1 µm) in saturation (VDS = 1 V) and front gate voltage VG = 0 V.

Fig. 15

 15 Fig. 15 Simulated phase of the back transadmittance efficiency versus normalized frequency Fnorm1 at different VbG for NMOS (L = 1 µm) in saturation (VDS = 1 V) and front gate voltage VG = 0 V.

Phase(y m2 /I D ) (deg)

V bG = 8 V IC = 13 V bG = 9 V IC = 24


Mattia Pesenti

Ziad Alkhoury

Maciej Bednarczyk

Hassan Omran

Bernard Bayle

Linear Parameter-Varying Identification of the EMG-Force Relationship of the Human Arm

Keywords: Human modeling, system identification, biomedical signal processing, electromyography

In this paper, we present a novel identification approach to model the EMG-Force relationship of the human arm, reduced to a single degree of freedom (1-DoF) for simplicity. Specifically, we exploit the Linear Parameter Varying (LPV) framework. The inputs of the model are the electromyographic (EMG) signals acquired on two muscles of the upper arm, biceps brachii and triceps brachii, and two muscles of the forearm, brachioradialis and flexor carpi radialis. The output of the model is the force produced at the hand actuating the elbow. Because of the position-dependency of the system, the elbow angle is used as scheduling signal for the LPV model. Accurate modeling of the human arm with this approach opens new possibilities in terms of robot control for physical Human-Robot Interaction and rehabilitation robotics.

I. INTRODUCTION

Skin electromyography (EMG) allows to acquire the electrical activity that is responsible of muscular contractions. This technique has become essential in a wide range of applications, as research on the human neuromotor mechanisms, clinical diagnostic methods [START_REF] Hogrel | Clinical applications of surface electromyography in neuromuscular disorders[END_REF] and neuro-prosthesis. Robotic rehabilitation devices show great potential when treating neurological injuries [START_REF] Bogue | Rehabilitation robots[END_REF][START_REF] Weber | The use of robots in stroke rehabilitation: A narrative review[END_REF], such as strokes or spinal cord injuries, as the treatment involves repetitive movements. Clinical studies show the increase of patient recovery when the robot provides active aid by physical interaction [START_REF] Kahn | Robot-assisted reaching exercise promotes arm movement recovery in chronic hemiparetic stroke: a randomized controlled pilot study[END_REF][START_REF] Basteris | Training modalities in robot-mediated upper limb rehabilitation in stroke: a framework for classification based on a systematic review[END_REF]. With the aim of improving the rehabilitation outcome, it would be possible to integrate an EMG-based model of the human limb into the robot control strategy.

EMG-driven models are developed considering as input the electromyographic signals acquired on the muscles of interest. The problem consists in determining the force developed by a set of muscles from the corresponding EMG signals. The EMG-Force relationship has been under research with two different approaches. The first one consists in phenomenological modeling by exploiting mathematical equations describing the neuromuscular system, such as the Hill model [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF]. This modeling approach is highly nonlinear and depends on several subject-specific parameters. Even though optimization techniques have been proven to improve the fitting of these parameters [START_REF] Buongiorno | A linear approach to optimize an EMG-driven neuromusculoskeletal model for movement intention detection in myo-control: A case study on shoulder and elbow joints[END_REF], the model structure remains of high complexity and nonlinearity. The second approach, used in this work, avoids the drawbacks of the first one and is based on black-box system identification. More precisely, the model parameters are identified directly from input/output datasets. Several works have considered the black-box approach to model the EMG-Force relationship. In [START_REF] Clancy | Identification of constant-posture EMGtorque relationship about the elbow using nonlinear dynamic models[END_REF], a nonlinear dynamic polynomial model was proposed to describe the EMG-Torque relationship around the elbow. Other works showed that parallel cascade Wiener models can perform well in the identification of the EMG-Force relationship considering both constant-posture trials [START_REF] Hashemi | EMG-force modeling using parallel cascade identification[END_REF] and free hand motions [START_REF] Hashemi | Enhanced dynamic EMG-force estimation through calibration and PCI modeling[END_REF]. Here, the elbow angle of the subject has been included as additional input in order to identify the dependency of the system on the geometric configuration of the arm. Indeed, the shape of the EMG-Force relationship scales depending on the elbow angle [START_REF] Liu | Influence of joint angle on EMG-torque model during constant-posture, quasi-constant-torque contractions[END_REF]. This modeling approach was extended also to the case of free hand motions [START_REF] Hashemi | Enhanced dynamic EMG-force estimation through calibration and PCI modeling[END_REF]: releasing the constant-posture constraint led to a dramatic decrease of the modeling accuracy.

The presented black-box models are still highly nonlinear and computationally complex, thus being unsuitable for the integration within the control system of rehabilitation robots. Moreover, the accuracy of the presented models is generally limited to the quasi-static case (isometric contractions), while the performance decreased considering dynamic contractions.

Linear Parameter Varying (LPV) modeling has been recently applied to study the human neuromusculoskeletal system. Indeed, the LPV approach was exploited to model the nonlinearity (due to a position dependency) of the joint impedance, as for the case of the wrist [START_REF] Van Eesbeek | LPV subspace identification of timevariant joint impedance[END_REF] and the ankle [START_REF] Golkar | Linear Parameter Varying Identification of Dynamic Joint Stiffness during Time-Varying Voluntary Contractions[END_REF]. The accuracy reported in these works serves as motivation to extend the LPV approach to the EMG-Force relationship of the 1-DoF human arm.

The present paper is structured as follows. In Section II, we introduce the class of LPV models. Then, in Section III, we describe the experimental setup. In Section IV, we present the main results, which are discussed in Section V. Finally, in Section VI, we conclude the paper.

II. PRELIMINARIES

The Linear Parameter Varying (LPV) representation provides a framework to describe certain classes of nonlinear systems. Many tools dedicated to Linear Time-Invariant (LTI) systems can be extended to the LPV framework. Indeed, LPV models can be considered as an extension of LTI models. The practical advantage of an LPV model is that its input/output relationship is linear, while the model itself is capable of describing nonlinear and time-varying systems.

A. LPV Models

A discrete-time, state-space LPV model is defined as follows

x(k + 1) = A (q(k)) x(k) + B (q(k)) u(k) (1a) y(k) = C (q(k)) x(k) + D (q(k)) u(k) (1b)
where x(k) ∈ IR n is the state vector, u(k) ∈ IR p is the input vector and y(k) ∈ IR m is the output vector, all at time instant k. The order of the LPV model (n) is the dimension of its state vector, while the integers p, m represent the input and output dimensions, respectively. The state matrix A (q(k)) ∈ IR n×n , input matrix B (q(k)) ∈ IR n×p , output matrix C (q(k)) ∈ IR m×n and feed-through matrix D (q(k)) ∈ IR m×p are all dependent on an external signal q(k), which is usually called scheduling signal. In this work, we consider the subclass of Affine LPV (ALPV) models [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF][START_REF] Alkhoury | Identifiability of affine linear parameter-varying models[END_REF]. In this representation, the model matrices depend in an affine way on the scheduling signal, namely

A (q(k)) = A 0 + A 1 q(k) (2)
and similarly for B(•) and C(•). Without loss of generality, the feed-through term is considered to be zero, i.e. D(•) = 0. We focus on the ALPV representation for two reasons: first, the realization theory for this subclass of LPV models is already mature [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF][START_REF] Alkhoury | Identifiability of affine linear parameter-varying models[END_REF]; in addition, the control theory [START_REF] Mohammadpour | Control of Linear Parameter Varying Systems with Applications[END_REF][START_REF] Sename | Robust control and linear parameter varying approaches, application to vehicle dynamics[END_REF] dedicated to ALPV models is well-developed. Moreover, many of the tools developed for the LTI case are extended to the LPV framework. This means that some LPV identification methods directly rely on identification tools developed for LTI models.

B. LPV Identification

LPV system identification can be divided into two main approaches based on the experiment design, namely the local and the global identification approaches [START_REF] Tóth | Modeling and identification of linear parametervarying systems[END_REF]. The global approach [START_REF] Van Overschee | Subspace identification for linear systems: TheoryImplementationApplications[END_REF][START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] requires performing only one experiment on the system, in which the input and the scheduling signals are excited in order to capture all the dynamics of the system. This is a natural way for obtaining datasets that describe the real system. In contrast, the local approach identifies multiple LTI models (usually called frozen LTI models) for multiple, constant scheduling signals. Then, the parameterdependent LPV model is fitted to the system matrices of the frozen LTI models using regressors formed from the scheduling signal q [START_REF] Lovera | Identification for gain scheduling: a balanced subspace approach[END_REF][START_REF] Zhu | A method of LPV model identification for control[END_REF]. The local approach can be teared down to the following three steps.

1) Identification of linear state-space models for several values of the scheduling parameter space. 2) Transformation of the obtained LTI models to a consistent state-space representation [START_REF] Sename | Robust control and linear parameter varying approaches, application to vehicle dynamics[END_REF]. 3) Interpolation of the state-space matrices of the local models with respect to the scheduling signal in order to obtain the LPV model. It is well-known that the local approach suffers from a systematic modeling error. However, this error is bounded and its bound depends on the speed of change of the scheduling signal [START_REF] Alkhoury | Comparing global input-output behavior of frozen-equivalent LPV state-space models[END_REF]. On the other hand, the local approach is more common for two reasons: it is more practical in terms of experimental design and it benefits of the well-developed theory of LTI systems identification.

Considering the LPV model here presented, the input of the system is a set of EMG signals, the output is the measured force, while the scheduling signal is a quantity proportional to the elbow angle of the subject (see §III).

III. EXPERIMENTAL SETUP

In our setup, input/output data are acquired while the user interacts with a collaborative robot, the KUKA LBR iiwa. EMG and force data are synchronously acquired and stored to be processed offline. Black-box identification methods are then used to identify a model mapping the muscular activation signals into the force developed by the arm.

The joint configuration of the robot is set in order to kinematically couple the human elbow angle with the elbow joint of the robot, as shown in Fig 1 . In this way, the user can rest their arm on the robot while standing next to it. Both the standing point and the arm resting point are marked for repeatability. This setup permits to achieve weight compensation of the arm without the use of any belts or braces. Moreover, the proportionality between the elbow angle of the user and the elbow joint of the robot allows to easily measure and vary the human elbow angle. In this way, the experimental setup allows to measure the interaction forces between the user and the robot in two modes. First, by fixing the robotic arm at a certain position (i.e. at a fixed angle), we can model the dynamics of the 1-DoF human arm in a constant-posture setting. This type of experiment is called local trial and it is used in the local approach for LPV identification. Second, we can apply a movement profile to the elbow joint of the robotic arm, which would lead to a time-varying elbow angle profile in the horizontal plane. This experiment is called hereafter global trial and results in datasets which are suitable for the global identification approach. During the experiments, the subjects are asked to stand still, maintaining a constant body posture, and focus on actuating only the elbow joint to generate force, while keeping the wrist at a fixed angle.

The EMG signals, which represent the input of the system, are acquired using the Delsys Trigno Wireless System (Delsys Inc) with a sampling rate of 2000 Hz. Four muscles of the arm are considered, namely biceps brachii (BIC), triceps brachii (TRI), flexor carpi radialis (FCR) and brachioradialis (BRD). The force, which is the output of the system, is acquired using a 6-DoF force/torque sensor (ATI Mini40, ATI Industrial Automation) mounted on the end effector of the robotic arm. For the model, we only consider the force component in the horizontal plane and perpendicular to the forearm of the subject. Real-time visual feedback of the force profile is provided to the user on a computer screen, in order to keep the applied forces within a maximum value of ± 20 N. Thanks to the aforementioned kinematic coupling, the coordinate of the elbow joint of the robot (q) is acquired to be used as the scheduling signal of the LPV model. The Fast Research Interface (FRI) [START_REF] Schreiber | The fast research interface for the kuka lightweight robot[END_REF] of the KUKA robot is exploited in order to acquire the force and scheduling signals at 100 Hz.

The study population consists in two male subjects with an average age of 35.5 years and an average Body Mass Index (BMI) of 24.5 kg/m 2 .

A. EMG Processing

The so-called muscular activation signal is derived from the raw EMG data for each muscle according to the following processing steps. The acquisition system applies a first bandpass filter between 20 and 450 Hz, in order to filter out offsets and baseline drifts as well as to avoid aliasing. Then, during the offline processing, the bandwidth of the signal is further reduced using a Butterworth filter that selects the frequencies between 20 and 350 Hz, implemented by means of the series of a fifth order high-pass and a second order low-pass. Full-Wave Rectification is later applied to the resulting signal by taking its absolute value. The last filtering step is the envelope extraction, which is performed using a second order Butterworth low-pass filter with a cut-off frequency of 1.775 Hz. Starting from state-of-the-art values, the parameters of the filtering stages have been optimized in order to get the highest modeling accuracy.

IV. IDENTIFICATION OF THE EMG-FORCE RELATIONSHIP

In this Section, the two LPV identification approaches exploited to obtain the model of the 1-DoF arm are described. The results of model validation have been computed by comparing the estimated force output with the measured signals in terms of goodness of FIT [%] and Variance Accounted For (VAF) [%].

A. Local Identification Approach

The basic idea of the local approach for LPV modeling is to consider the nonlinear, position-dependent system as a set of Linear Time-Invariant (LTI) models for different working points of the scheduling signal q. Hence, a successful description of the EMG-Force relationship at constant values of q by means of LTI models is essential to build the LPV model following this approach. Therefore, several identification experiments are performed at constant angular positions. Specifically, 8 local datasets are acquired at q = 10, 20, 25, 30, 40, 45, 50, 60 [ • ]. Considering the human elbow angle, q = -40 • corresponds to the full extension of the arm (i.e. 180 • ). The local LTI models are identified at each of these angular positions using numerical algorithms for state-space models based on sub-space identification (N4SID) [START_REF] Van Overschee | Subspace identification for linear systems: TheoryImplementationApplications[END_REF]. The identification is done in free form with prediction focus, thus minimizing the 1-step ahead prediction error.

1) Local LTI Models: The average estimation performance of the LTI models on the local datasets resulted in 87.75 % of FIT and 98.32 % of VAF. It is interesting to evaluate the accuracy of these local models as a function of the distance ∆q = q I -q V between the identification position q I and the validation position q V . More precisely, the LTI model identified at q I = 20 • is validated using the datasets acquired at q V = 10, 40, 60 [ • ]. As shown in Fig. 2, the performance of the LTI model drops as ∆q increases. In particular, the FIT is much more sensitive to the angular distance. The accuracy of the local linear models is acceptable for low displacements from the identification position, while it decreases dramatically for ∆q > 35 • . This corroborates the modeling hypothesis of a positiondependent system: the EMG-Force relationship does depend q = 60°, q = 40°F ig. 2. Force plots comparing the real output (solid black) with the estimated output of the LTI model identified at q = 20 • and validated at q = 10 • (top panel), q = 40 • (middle panel) and q = 60 • (bottom panel). In the plots, the green line corresponds to the LTI model identified with p = 4 inputs, while the orange line refers to LTI models identified from the signals of biceps and triceps only (p = 2). on the elbow angle.

2) Locally-identified LPV Model: The sub-space method used for the identification of the local LTI models may have produced LTI state-space models belonging to different bases. Hence, each LTI state-space model must be transformed to the same basis before the interpolation, in order for this to make sense. In this work, a Gramian-based balancing transformation of the state-space realizations is used [START_REF] Laub | Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms[END_REF][START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF]. After this, the coefficients of the ALPV model are estimated by means of the Levenberg-Marquardt nonlinear least-squares algorithm.

The LPV model obtained by means of the Local Identification Framework (denoted as LIF-LPV) can be validated with a twofold approach. On the one hand, the local dataset acquired at q = 40 • is used to evaluate the accuracy of the LPV model obtained by interpolation. On the other hand, it is more significant to perform model validation considering global trials, as they describe a more realistic scenario. In the former case, the average FIT of the LIF-LPV model is equal to 84.32 % and its average VAF is equal to 97.80 %. On the local trials, the validation performance of the LPV model is very similar to the estimation performance of the local LTI models, which consists in an average FIT equal to 85.83 % and an average VAF equal to 98.02 %. On the other hand, in the latter case, the validation performance of the LIF-LPV model decreases considerably. In particular, the average FIT is 74.10 % and the average VAF is 95.43 %, as reported in Table I. It can be noticed a different behavior of the two indicators: while the average VAF of the LIF-LPV model remains close to the performance obtained during local validation, the FIT drops significantly.

B. Global Identification Approach

The Global Identification Framework for LPV modeling is carried out exploiting the Prediction-Based Subspace Identification (PBSID) toolbox [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF] and the LPVcore toolbox [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF], both developed for MATLAB. The output of the PBSID toolbox is used to initialize the computation of the LPVcore toolbox. This allows to refine the estimation of the LPV model identified globally. This algorithm identifies the statespace LPV model by means of a gradient-based optimization technique. The LPV model obtained by means of the Global Identification Framework (denoted as GIF-LPV) is again a discrete-time, second order, ALPV model.

For each user, the acquisition of a global dataset is repeated four times. In this way, one of these datasets is used to identify the LPV model according to the global approach, while the other three are used to perform model validation.

The average FIT and VAF of the GIF-LPV model are 78.25 % and 95.66 %, respectively. As expected, the global identification approach for LPV modeling gives better performance with respect to the local one considering globaltrial validation. This is due to the different datasets used for the identification of the LPV model according to the two approaches. In particular, the GIF-LPV model is identified using a global dataset containing information about the arm dynamics during time variations of the elbow angle. On the other hand, this information is not contained in the localtrial datasets, as the arm is in a fixed position. Moreover, it should be noted that, while the performance of the LIF-LPV model decreased when validated on a global trial, the GIF-LPV model performs well also if validated on datasets obtained with local trials.

C. Performance Comparison

In this Section, the performance of the LPV model are compared to the results obtained by a Hammerstein model, which is identified as a reference term for the LPV approach. Moreover, it is interesting to evaluate the validation performance on global trials of the LTI approach, given its accuracy in the case of constant-posture trials (see §IV-A.1).

The Hammerstein model is identified using the System Identification Toolbox of MATLAB. On the other hand, the LTI model (hereafter, globally-identified LTI or gLTI for brevity) is identified using again the N4SID algorithm. For both models, the identification dataset is the same global trial used also in the case of the GIF-LPV model.

The performance of the globally-identified LTI model are presented first. The gLTI model outperforms the LPV model identified following the local approach (LIF-LPV), with an average FIT of 76.75 % and an average VAF of 94.75 %. Although this result may sound surprising, the reason behind it is twofold. The first one is the difference between the datasets used to identify the two models. In particular, the LIF-LPV model is identified using a batch of local-trial datasets, in which the scheduling signal is constant and the muscles contract in isometric conditions. On the other hand, the LTI model is in this case identified on a global trial, thus being able to describe the dynamics of the system due to the time variations of the scheduling signal and to the non-isometric muscular contractions. The second reason is what we call here the linearization effect actuated by the muscles of the forearm. These muscles (BRD and FCR) are not included in the models of the 1-DoF arm presented in the literature. Indeed, their effect is generally neglected as the wrist motion is externally constrained. On the other hand, they have a non-negligible role in the process of endpoint force production considering the 1-DoF arm. Moreover, their EMG-Force relationship is not influenced by the elbow angle. Consequently, their contribution helps reducing the overall nonlinearity of the system. This is confirmed by evaluating the performance of the LTI model identified using the activation signals of biceps and triceps only (i.e. with p = 2 inputs) instead of using the full set of input signals (i.e. p = 4). In the former case, the linear model is much more sensitive to the elbow joint coordinate and hence to ∆q. This is evident comparing the force signals in Fig. 2.

The results of the Hammerstein model are presented at last. This modeling approach scores an average FIT of 82.94 % and an average VAF of 97.09 % on global trials. As expected, the nonlinear model outperforms both the LTI and the LPV models, allowing the highest modeling accuracy. The validation results on global trials of all the presented models are summarized in Table I.

V. DISCUSSION

The first interesting result presented in this work is that local, linear models (LTI) are accurate enough to describe the EMG-Force relationship for constant-posture trials. This demonstrated the feasibility of modeling the human arm exploiting the Local Identification Framework for LPV models. Although feasible, the performance in terms of FIT decreased dramatically considering global validation trials instead of constant-posture trials. As argued above, this increase of the estimation error is due to the lack of information about the time-varying dynamics of the elbow angle in the local datasets, which also results in non-isometric muscular contractions. This is the main limitation of the local approach for LPV identification. Indeed, the behavior of the true LPV system cannot be fully determined interpolating a set of local models identified at several values of the scheduling signal, as the system is sensitive not only to q, but also to its time variations. Aside to this intrinsic limitation of the local identification approach, also the behavior of the muscles changes from local to global trials. Indeed, the constantposture constraint is released in the latter case. Therefore, the muscles work in isometric conditions in the case of the identification trials and in non-isometric conditions in the case of the validation trials. Hence, the more complex dynamics of the validation data cannot be described by the LIF-LPV model, since it is not contained in the dataset used for its identification.

As expected, the LPV model identified globally (GIF-LPV) outperforms both the LIF-LPV and the gLTI. Nevertheless, the Global Identification Framework is generally a less practical choice compared to the Local Identification Framework for LPV modeling. Moreover, the global identification of LPV models is still under ongoing research, affecting especially the development of algorithms and tools for the identification process itself. Consequently, the performance of this approach is expected to increase in the near future.

Another important result presented in this work is the importance of including the muscles of the forearm, namely the brachioradialis and the flexor carpi radialis, while releasing the external constraints on the wrist movement. Indeed, these muscles have an important role in force generation. This is mainly due to the unconscious, synchronous activity they manifest while biceps and triceps are consciously activated when the subjects interact with the robot. This may be partially related to the fact that the wrist was not constrained during the acquisitions. One could argue that the identified model actually describes the quasi-2-DoF arm, because of the activation of the muscles of the forearm actuating the wrist. Nevertheless, the aim of the identified model is to provide a sufficiently accurate estimation of the force produced by the human arm considering an experimental setup similar to a real-case scenario of HRI. As a consequence, it was chosen to let the wrist free from external constraints. As shown in this work, this aim has been achieved. The muscles of the forearm, in fact, induce a linearization effect of the EMG-Force relationship, thus allowing higher accuracy with a lower complexity modeling approach.

VI. CONCLUSION AND FUTURE DEVELOPMENTS

In this work, we have addressed the problem of modeling the EMG-Force relationship of the 1-DoF human arm with a low-complexity, control-oriented approach, by means of Linear Parameter Varying (LPV) models. These results show that an LPV model can be a good compromise in terms of accuracy and complexity with respect to state-of-theart alternatives. Releasing the external constraints off the wrist and including the forearm muscles appeared to be fundamental to achieve such accurate results. The reported performance metrics legitimate this approach, that can be further investigated with the aim of integrating the identified model into the control system of collaborative robots.

For the future, several improvements to the experimental setup could be done. First of all, the number of subjects of the study could be increased in order to have statistically significant results. Moreover, it could be useful to have a direct measurement of the elbow angle of the human subject. This could allow to investigate also different interaction scenarios between the robot and the human operator, as the kinematic coupling would be no more necessary. Finally, the model could be extended adding the degree of freedom due to the actuation of the wrist and/or the shoulder.

Fig. 1 .

 1 Fig. 1. Schematic representation of the acquisition setup: the user rests his arm on the KUKA robot while standing. EMG and force are synchronously acquired. The span of the scheduling signal q is shown with a dash-dotted line. A computer screen in front of the user provides the visual force feedback in real time.

Fig. 3 .

 3 Fig. 3. Validation of the locally-identified LPV model on a local trial (acquired at q = 40 • ). The measured force (solid black) is compared to the signal estimated by the LPV model (solid green).

Fig. 4 .

 4 Fig. 4. Validation force plot comparing the GIF-LPV model with the Hammerstein model and the globally-identified LTI model. The scheduling variable trajectory q(t) is shown in the bottom panel of the plot.

TABLE I VALIDATION

 I RESULTS COMPARING THE PERFORMANCE ON GLOBAL TRIALS OF LTI, LIF-LPV AND GIF-LPV. THE PERFORMANCE OF THE HAMMERSTEIN MODEL ARE REPORTED AS REFERENCE TERM.

			LTI	LIF-LPV	GIF-LPV	Hammerstein
	User 1	FIT [%] VAF [%]	72.61 93.12	67.77 94.63	75.37 94.88	80.12 96.19
	User 2	FIT [%] VAF [%]	80.89 96.38	80.42 96.22	81.12 96.44	85.76 97.98
	Average	FIT [%] VAF [%]	76.75 94.75	74.10 95.55	78.25 95.66	82.94 97.09

ACKNOWLEDGMENTS

This work has been published under the framework of the IdEx Unistra and benefits from a funding from the State managed by the French National Research Agency as a part of the Investments for the future program.