S Augé

A Tamra

L Rigal

V Lobjois

B Ducommun

D Dubuc

K Grenier

Experimental Dosimetry Study of a Miniature RF Applicator Dedicated to the Evaluation of Severe RF Exposure Impact on a 3D Biological Model

Keywords: micro-dosimetry, RF applicator, microtechnologies, RF effects, biological cells, micro-tissues, SAR

The impact of RF radiation on the living constitutes a controversial topic, whereas wireless applications continue to grow drastically and our daily environment is surrounded of electromagnetic fields. In order to evaluate the potential impact of RF radiations even with severe RF exposures, it is important to conduct experiments with a very well controlled dosimetry and with calibrated RF exposure systems. In this paper, a miniature RF applicator, which enables to apply controlled electromagnetic fields radiation on micro-tissues, is introduced. An experimental method based on thermal measurements has been developed in such a micro-device to define the Specific Absorption Rate applied to micro-tissues depending on the applied power and the thermal increment in steady state. Due to miniature configuration of the applicator, a large range of SAR may be achieved in moderate thermal conditions.

I. INTRODUCTION

ElectroMagnetic (EM) fields are widely used for civil and military applications. In addition to the wireless standards and applications, such as the 5G transfer data generation, the connected objects, wireless power transfer and wireless implants, radars and EM-based weapon equipment for electronic warfare for example are increasing of interest and development. It includes not only low or medium RF exposure ranges (with Specific Absorption Rate -SAR-up to few W/kg, below or close to ICNIRP recommendations), but also severe ones with kW/kg values. The evaluation of potentially harmful effects on the living due to electromagnetic fields is consequently imperative. This is even more true that the interaction phenomena between EM fields and the living are far from being completely known.

Numerous studies related to RF exposure systems and RF effects have already been conducted in the past [START_REF] Wu | Safe for generations to come[END_REF][START_REF] Paffi | Considerations for Developing an RF Exposure System: A Review for in vitro Biological Experiments[END_REF][START_REF] Zhadobov | Exposure System and Dosimetry for In Vitro Studies of Biocompatibility of Pulse-Modulated RF Signals of Ultrahigh Field MRI[END_REF][START_REF] Perrin | Evaluation of the co-genotoxic effects of 1800 MHz GSM radiofrequency exposure and a chemical mutagen in cultured human cells[END_REF][START_REF] Schuderer | In vitro exposure systems for RF exposures at 900 MHz[END_REF][START_REF] Habauzit | Transcriptome Analysis Reveals the Contribution of Thermal and the Specific Effects in Cellular Response to Millimeter Wave Exposure[END_REF][START_REF] Merla | Microdosimetry for nanosecond pulsed electric field applications: a parametric study for a single cell[END_REF]. Some are revealing harmful effects, with even the appearance of different cancers on rats [START_REF] Lin | Cancer Occurrences in Laboratory Rats from Exposure to RF and Microwave Radiation[END_REF][START_REF] Wyde | Report of Partial findings from the National Toxicology Program Carcinogenesis Studies of Cell Phone Radiofrequency Radiation in Hsd: Sprague Dawley® SD rats (Whole Body Exposure)[END_REF][START_REF] Lin | Human exposure to RF, microwave and millimeterwaveelectromagnetic radiation[END_REF][11], whereas others do not present any EM effect on the living. To contribute to this controversial topic, it is important therefore important to use or develop [START_REF] Wu | Safe for generations to come[END_REF] precise and controlled RF exposure systems, which are well characterized in terms of SAR, EM field intensity, thermal increment, with a clear description of the applied RF signal, [START_REF] Paffi | Considerations for Developing an RF Exposure System: A Review for in vitro Biological Experiments[END_REF] with high range of SAR values, low, medium and up to severe level RF exposures.

In vitro evaluations of EM fields are traditionally performed with cell cultures grown as monolayer (2D) in Petri dishes. This cellular configuration presents however limitations, as far from real in vivo behavior. Moreover, due to ethic issues, there is a strong interest in developing RF exposure systems, which could avoid the use of animal testing while being biologically relevant. Employing micro-tissues, e.g mastered organized 3D aggregates of cells, is particularly attractive, as they present a strong cellular heterogeneity, which confers an intermediate level of this model toward in vivo investigations.

To our best knowledge, there is no RF applicator, which has been developed so far to radiate any controlled RF signals on micro-tissues from low to severe level RF exposures. This paper introduces such a structure and a method applied to characterize its SAR and its thermal increment depending on the applied RF power. This paper first introduces the advantages of using a 3D biological model, called micro-tissue. A description and a microwave characterization of the developed RF exposure applicator appropriate to radiate micro-tissues with controlled EM fields is then be given in the next section, as well as its technology. The section 4 presents then the experimental

CPW access Fluidic walls

Micro-tissue location dosimetry study developed to characterize the SAR and thermal increment within the miniature RF applicator. Finally, RF exposures are performed on micro-tissues and evaluated in term of genotoxic effect with DNA damages detection.

II. MATERIALS AND METHODS

A. Micro-tissue: a 3D biological model

Traditionally, 2D monolayer cell cultures are used for cellular-based dosimetric investigations. However, these biological models are not relevant of in vivo behavior.

On the other side, 3D multicellular models, also called micro-tissues, are constituted of an aggregate of multiple cells and exhibit important cellular heterogeneity, regionalization of growth and strong interactions between cells. The main interest of such a model is that it may reproduce in vitro the heterogeneity of a tissue at a lower complexity level, without requesting in vivo investigations and therefore animal testing.

Such 3D micro-tissues exhibit a rounded shape and present increasing diameters as they continue to proliferate and grow naturally while kept in culture medium with appropriate environmental conditions. Due to microfabrication constraints, 250 to 300 µm large micro-tissues have been chosen for our experiments.

B. Architecture description of the RF applicator, specific to radiate micro-tissues with controlled RF signals

As shown with the photography of Fig. 1, the miniature RF applicator is composed of a first coplanar section, which enables the positioning of a RF probe. While the coplanar grounds are merged on the opposite side of the device, the central conductor is ended by a capacitive patch in the middle of two fluidic walls. The micro-tissue under test is localized in the center part of the conductive open-ended transmission line section. The fluidic walls permit to maintain the micro-tissue in a culture medium, providing all nutrients necessary to the living cells. As the fluidic walls present a thickness of 500 µm, there is no liquid where the RF probe is applied. All dimensions are given in Fig. 2.

For the experiments, the fluidic channel is loaded with the culture medium of the micro-tissues, Dulbecco's Modified Eagle Medium (DMEM) in our case. To characterize the fabricated microdevice up to 40 GHz, the device under test is connected to a Vector Network Analyser using a coplanar RF probe and a coaxial cable, using a preliminary SOLT calibration step. The measured S11 module is then compared to the simulated one in Fig. 3, while the structure is loaded with DMEM. The good agreement between the two curves validates the fabrication of the micro-device. At 2.45 GHz, the device is highly reflective (|S11|=0.974 at 2.45 GHz). This is due to the low DMEM volume fraction combined with the moderate loss tangent of DMEM, exhibiting a 0.27 value at this frequency. This will have an impact on the SAR values, as only 5% of the incident power is dissipated within the DMEM liquid. Nevertheless, section IV will experimentally demonstrate that high (severe) SAR values may be reached.

C. RF applicator technology

The RF applicator fabrication involves three main steps based on standard microtechnology processes. First, the RF applicator is realized on top of a quartz substrate, which has been preliminary cleaned with different chemical baths. The RF metallization is then performed using a lift-off technique. The obtained waveguide is composed of a flash of titanium to enable the metal adhesion on the substrate, followed by a gold layer of 0.3 µm thick. To avoid any spreading and disassembling of the micro-tissue within the fluidic channel, a 5 µm thick isolation layer is deposited and patterned. 500 µm fluidic walls are finally realized on top, while keeping the micro-tissue in its liquid environment during the experiments. Moreover, the choices of all constitutive materials have been done to make the structure totally biocompatible.

III. EXPERIMENTAL DOSIMETRY METHODOLOGY

To experimentally predict the SAR level within the microstructure, a methodology has been selected, which consists in measuring the transient temperature immediately after RF exposure.

A. Experimental test setup

In the homogeneous host liquid, which is initially at thermal equilibrium, the SAR is related to the derivative of the temperature with respect to time, as long as no significant effect due to thermal conduction appears. This last requirement is obtained just after applying the RF exposure. The corresponding relationship is given by the following equation:

𝑆𝐴𝑅 = 𝑐 '( ') (1)
where the SAR is in W/kg, c is the specific heat capacity of the liquid in J/kg/K, dT is the change in temperature in K during dt in seconds.

To perform the dosimetry investigation, a test setup has been defined. The corresponding schematic is given in Fig. 4. The experimental test setup is composed of two parts, one dedicated to the RF signal generation and exposure, whereas the second one is related to the thermal measurements.

The RF power delivery part is composed of a RF generator followed by a power amplifier, which both enable to deliver a power level up to 30 dBm. The RF signal is then conveyed to the miniature RF applicator through a directional coupler, a cable and a RF probe. The coupler is in charge of monitoring the reflected power (5% of the incident one at 2.45 GHz) and permits (1) to monitor the dissipated power and (2) to detect the malfunction of the micro-device.

The thermal part includes a fiber optic temperature sensor connected to a signal conditioner. The tip of the sensor only exhibits a size of 120 µm wide, which is compatible with the 300 µm width of the fluidic channel. This fiber is placed in the vicinity of the isolation layer and permits to accurately measure the temperature modification within the fluidic channel, in place of the micro-tissue. Temperature versus time may then be measured with such an equipment.

B. SAR prediction methodology

The methodology is based on the analysis of temperature readings as a function of time, before and after RF exposure. An example of obtained curve is given in Fig. 5.

The SAR is extracted from the slope of the temperature versus time curve as soon as the RF exposure is applied. This slope is indicated with the orange dashed line in Fig. 5. SAR/c therefore corresponds to the "a" coefficient of the curve. Such a measurement also enables to evaluate the thermal increment in steady state (ΔT indicated in red in Fig. 5).

Moreover, the SAR depends on the power generated at the output of the power amplifier. Theoretically, a proportional relationship exists between these two quantities. Each evaluated SAR will consequently be normalized with respect to the injected power. Next section is dedicated to the experimental results.

IV. EXPERIMENTAL DOSIMETRY RESULTS

A RF signal at 2.45 GHz in CW mode has been used for the experiments. Beside the SAR, it is important to know if the experiment is in a thermal or an athermal regime. This is the reason why the thermal increment has been studied, which is also theoretically proportional to the output power of the amplifier. A 20 dBm source power level of the output power of the power amplifier has been tested. This enables a sufficiently sensitive measurement to determine both thermal increment and the slope of the time-dependent temperature curve.

Table 1 presents the experimental results obtained for the thermal increment and the SAR, both normalized with respect to the source power in mW. A non-significant thermal increment inferior to 0.1°C corresponds to SAR value of several hundred W/kg. Consequently, in athermal conditions and using the microdevice combined to medium power RF equipment, SAR values up to hundred, even thousand W/kg may be reached, which may be relevant for specific severe RF exposures studies (such as microwave based anti-personnel weapons for instance in military applications). With microfluidics, the surface area to volume ratio indeed increases, which translates to high SAR values generated with a moderate thermal increment.

Additionally, the "microscale" setup may also generate conditions in the range of W/kg and below, with thermal increase lower than one thousandth degree, which is highly non-significant.

V. BIOLOGICAL EXPERIMENTS After RF exposure, RF radiated micro-tissues are collected and analysed with traditional biological equipment. Depending on the employed techniques, further understanding of the interaction mechanisms of EM fields on living cells may be reached.

In this experiment, RF genotoxic effects are targeted. Four replicated sample types are considered:

-the RF exposed micro-tissues, -others, which are similarely prepared and placed apart to constitute the negative control, -the mock ones, which correspond to micro-tissues submitted to all the manipulations involved during the experiment without RF illumination, e.g. also placed in the micro-devices without application of EM fields, -and finally, the positive controls, i.e. micro-tissues submitted to a chemical treatment known as generating DNA damages. Fig. 6 presents the photography of a micro-tissue positioned in the center of the RF applicator. The micro-tissue is maintained in its culture medium, visible on the picture as the pink liquid within the fluidic channel. After manipulation or RF exposure, micro-tissues are fixed and cryo-sectioned, then analysed by immunofluorescence using the γH2AX label, a nuclear marker which is specific to the detection of DNA strand breaks (DNA damage). Fig. 7 presents examples of images for each type of micro-tissues. As expected, the positive control presents many positive nuclei, related to DNA damages, whereas the negative control, the mock and the RF exposed micro-tissues do not present significative labelling for an applied SAR of 59 kW/kg for 15 min at 2.45 GHz.

VI. CONCLUSIONS

To conclude, this paper presents a miniature RF applicator, which enables to apply controlled electromagnetic fields on micro-tissues. Such a biological model presents the major advantage of exhibiting a cellular complexity and heterogeneity close to the one present in living entities, such as animals and humans, without requesting animal testing. To perform a dosimetry study of the device, a dedicated test setup and a measurement method have been defined to obtain both thermal increment and SAR values for different applied powers. Due to the miniature size of the device, a large range of SAR values may be reached in athermal condition, which may present interest for specific applications, including severe RF exposures. Finally, such a device allied to biological postanalysis on micro-tissues enables investigations toward a further comprehension of the EM fields interaction mechanisms at the cellular level.

Micro-tissue in culture medium

Fig. 1 .

 1 Fig. 1. Photography of a RF applicator adapted to micro-tissue.

Fig. 2 .

 2 Fig. 2. Schematic of the RF applicator suited for micro-tissue exposure.

Fig. 3 .

 3 Fig. 3. Simulated and measured S parameter of the device loaded with cell culture medium.

Fig. 4 .

 4 Fig. 4. Test setup used to evaluate the SAR within the fluidic cavity.

Fig. 5 .

 5 Fig. 5. Example of a temperature modification measured within the fluidic cavity versus time.

Fig. 6 .

 6 Fig. 6. Photography of an RF applicator under test, loaded with a micro-tissue.

Fig. 7 .

 7 Fig. 7. Example of labeled images obtained after γH2AX immunolabelling for the 4 micro-tissues samples : negative control, positive control, mock and RF exposed. Scale bar : 100 µm.

Table 1 .

 1 Normalized thermal increment and SAR value.

ACKNOWLEDGMENT This work was supported in part by French projects, grant Plan Cancer Env 2012-RF3DDR, grant n° EST-16 RF-07 from ANSES, and grant n° ANR-15-ASTR-0017 from ANR and DGA, and in part by LAAS-CNRS micro and nano technologies platform, member of the French RENATECH network.


Davide Invernizzi
email: davide.invernizzi@polimi.it

Marco Lovera
email: marco.lovera@polimi.it

Luca Zaccarian
email: zaccarian@laas.fr

Hierarchical dynamic control for robust attitude tracking

In this paper robust attitude tracking for fully actuated rigid bodies is addressed. By exploiting the cascade structure of the underlying mathematical model, a hierarchical framework including a large number of dynamic feedback controllers is proposed. The closed loop results in error dynamics comprising an inner loop associated with the angular velocity error, and an outer loop associated with the attitude error. We then establish sufficient conditions for solving the attitude tracking problem from an (almost) global perspective by leveraging recent results on stabilization of nonlinear cascades and invariance principles for differential inclusions. The modular nature of the proposed approach allows one to conclude stronger stability properties than those available for existing dynamic control laws, such as PID loops, often employed in applications.

INTRODUCTION

The attitude control problem for systems that can be modeled as fully actuated rigid bodies by a control torque has been addressed by a huge body of literature. Applications enjoying this structure range from aerospace, robotics and underwater vehicles. Most of the works dealing with attitude control from a global perspective are concerned with the design of stabilizing controllers in an ideal setting, which is an already demanding challenge: due to the topology of SO(3), global stabilization by means of continuous state feedback is structurally impossible [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF][START_REF] Casau | Robust global exponential stabilization on the ndimensional sphere with applications to trajectory tracking for quadrotors[END_REF]. When accounting for parametric uncertainties and exogenous disturbances, the attitude control problem becomes even more involved.

Robust attitude control designs have begun to emerge in recent years by extending techniques developed for Euclidean spaces [START_REF] Camblor | A Comparison of Nonlinear PI and PID Inertia-Free Spacecraft Attitude Control Laws[END_REF][START_REF] Goodarzi | Geometric nonlinear PID control of a quadrotor UAV on se(3)[END_REF][START_REF] Maithripala | An intrinsic PID controller for mechanical systems on Lie groups[END_REF]. Exploiting the full actuation assumption, the control torque is designed to compensate for the gyroscopic torque and to enforce a decrease condition of suitable Lyapunov functions along the flow of the closed-loop system. However, the resulting control laws are typically complex and not easy to tune due to nontrivial conditions on the gains: one may not be able to enforce a desirable level of performance around the desired attitude without violating sufficient conditins guaranteeing global stability.

In this work we appeal to ideas coming from reduction theory [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF] to propose a hierarchical control design based on an inner-outer loop paradigm that exploits the cascade structure of the dynamical attitude control model. In this way, we split the challenges of the control design in two simpler sub-problems: A) stabilization of the attitude error kinematics, which evolves in a boundaryless compact manifold and B) stabilization of the angular velocity error dynamics, which evolves in a Euclidean space. Although inner-outer loop paradigms already exist for attitude control [START_REF] Rudin | A robust attitude controller and its application to quadrotor helicopters[END_REF][START_REF] Forbes | Passivity-based attitude control on the special orthogonal group of rigid-body rotations[END_REF], the proposed approach is different in several ways including the fact that it encompasses a more general class of dynamic stabilizing laws. We start by developing a modular cascade control law based on static and dynamic feedback for attitude and angular velocity tracking, respectively. By requiring boundedness of the outer loop output (a desirable property), we are able to embed the closed-loop solutions in the funnel generated by a differential inclusion enjoying a cascaded property, without a priori restricting the class of allowable inner loop controllers. Recent results about Input-to-State Stability (ISS) for multistable systems [START_REF] Angeli | Stability robustness in the presence of exponentially unstable isolated equilibria[END_REF] allow us to point out basic properties that each component of the control law should satisfy to guarantee closed-loop stability. In particular, we show that the widely adopted attitude stabilizer corresponding to the gradient of the modified trace function on SO(3) [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF]] fits within our framework. Finally, we suggest a rather general class of dynamic controllers for angular velocity tracking, combining passivity concepts and invariance principles for differential inclusions [START_REF] Seuret | A nonsmooth hybrid invariance principle applied to robust event-triggered design[END_REF]. For instance, our control design makes it possible to prove almost global tracking, in the presence of constant disturbances and regardless of the gain values, for a control law behaving like a P/PI cascaded linear controller near the desired attitude motion.

NOTATIONS AND MATHEMATICAL PRELIMINARIES

R (R >0 , R ≥0 ) denotes the set of (positive, nonnegative) real numbers, R n denotes the n-dimensional Euclidean space and R m×n the set of m × n real matrices. The i-th vector of the canonical basis in R n is denoted as e i (vector with 1 in the i-th position and zeros elsewhere) and the identity matrix in R n×n is

I n := [e 1 • • • e i • • • e n ]. Given x ∈ R n , y ∈ R m , we denote (x, y) := [ x y ] ∈ R n+m . For x = (x 1 , . . . , x n ) ∈ R n , x := x 2 1 + . . . + x 2 n is its Euclidean norm, x ∞ := max i∈{1,...,n} (|x i |) is its Euclidean norm while for a matrix A ∈ R n×n , tr(A) is its trace, A F := tr(A A) is the Frobe- nius norm and skew(A) := A-A 2 is the skew-symmetric part of A. The set SO(3) := {R ∈ R 3×3 : R R = I 3 , det(R) = 1}
denotes the third-order Special Orthogonal group while S n := {q ∈ R n+1 : q = 1} denotes the n-dimensional unit sphere. The tangent space to SO(3) at any R ∈ SO(3) is denoted as T R SO(3). The normalized distance with respect to I 3 , induced by the Frobenius norm, is denoted as R SO(3) := 1 8 R -

I 3 F = 1 4 tr(I 3 -R) ∈ [0, 1]. Given ω ∈ R 3 , the S map S(•) : R 3 → so(3) := {Ω ∈ R 3×3 : Ω = -Ω } is such that S(ω)y = -S(y)ω = S(y) ω = ω × y, ∀y ∈ R 3
where × represents the cross product in R 3 . The inverse of the S map is the denoted as S(•) -1 : so(3) → R 3 . We use standard comparison functions from [START_REF] Khalil | Nonlinear systems, 3rd Edition[END_REF]: s) is decreasing for each r ≥ 0 and converges to zero as s → +∞.

function α : R ≥0 → R ≥0 is of class K if it is continuous, zero at zero and strictly increasing. It is of class K ∞ if it is also unbounded. β : R ≥0 × R ≥0 → R ≥0 is of class KL if r → β (r, s) is of class K for each s ≥ 0 and s → β (r,

RIGID BODY ATTITUDE CONTROL PROBLEM

Consider an inertial reference frame and a body fixed frame the origin of which is located at the center of mass of a rigid body. The configuration of the rigid body is described by the orientation of the body fixed frame with respect to the inertial frame, and it is globally and uniquely represented by a rotation matrix

R := [ b 1 b 2 b 3 ] ∈ SO(3) where b i ∈ S 2 , i ∈ {1, 2, 3}
are three orthogonal vectors with components resolved in the inertial frame forming a right-handed triad. The attitude motion of the rigid body is described by the following well-known dynamics:

Ṙ = RS(ω) (1) J ω = -S(ω)Jω + τ c + τ e (2)
where J = J ∈ R 3×3 >0 is the inertia matrix expressed in the body frame, ω ∈ R 3 is the body angular velocity, τ c ∈ R 3 is the control torque exerted by the actuators and τ e ∈ R 3 is the disturbance torque accounting for unknown exogenous effects. Our goal is to introduce dynamic controllers guaranteeing attitude tracking for a fully actuated rigid body in the presence of disturbances. For such a problem, we assume the following standard assumption for the desired attitude trajectory that should be tracked. Assumption 1. (Smoothness and boundedness of the desired trajectory). The desired trajectory t

→ (R d (t), ω d (t)) ∈ SO(3)× R 3 satisfies Ṙd (t) = R d (t)S(ω d (t)) ∀t ≥ 0 and t → ω d (t) is continuously differentiable and uniformly bounded.
The state feedback dynamic attitude tracking problem for a fully actuated rigid body can be then formalized as follows. Problem 1. Consider the attitude dynamics in equations (1)-(2). Given a desired trajectory t → (R d (t), ω d (t)) ∈ SO(3) × R 3 satisfying Assumption 1, design a state-feedback dynamic controller delivering a control torque τ c ∈ R 3 such that the trajectory t → (R d (t), ω d (t)) is asymptotically tracked.

It is worthwhile to mention that the attitude kinematics model (1) evolves on SO(3), a boundaryless compact manifold (not diffeomorphic to a Euclidean Space), while on the other hand, the angular velocity dynamics (2) evolves on a linear manifold R 3 and, in practical applications, is affected by exogenous disturbances and parameter uncertainties. The motivation for the solution proposed here is to parametrize hierarchical dynamic control architectures allowing for the use of PI-like robust control solutions independently addressing these uncertainty problems for each one of the two subsystems in (1)-(2).

CONTROL LAW AND STABILITY ANALYSIS

We propose here a hierarchical control architecture solving Problem 1 based on an inner-outer loop paradigm. Our framework leaves the door open for different selections of stabilizers of the inner and outer loops. Some examples are given in Section 5.

Control law and closed-loop dynamics

Let the attitude and angular velocity error coordinates be defined, respectively, as

R e := R T d R ∈ SO(3) (3
)

ω e := ω v -ω ∈ R 3 , (4) and consider the following control law ẋc = γ c (x c , ω e , ω v (R e ,t)) (5) τ c := S(ω)Jω + J ωv (R e , ω e ,t) + γ ω (x c , ω e , ω v ) (6) where ω v (R e ,t) := γ R (R e ) + R e ω d (t), (7)
whose derivative ωv , exploiting Ṙe in ( 9), proven in Proposition 1, can be expressed as

ωv (R e , ω e ,t) := γR (R e , ω e ) -S(γ(R e ) -ω e )R e ω d (t) + R e ωd (t) (8) and x c ∈ R n c is the controller state and γ R (•) : SO(3) → R 3 , γ ω (•, •, •) : R n c × R 3 × R 3 → R 3 , γ c (•, •, •) : R n c × R 3 × R 3 → R 3
are continuous stabilizers to be defined.

The inner-outer loop architecture of controller ( 5)-( 6) is represented in Figure 1: the inner loop ( 6) is in charge of assigning a suitable torque τ c to track the angular velocity reference ω v provided by the outer loop (7). The controller dynamics in ( 5) is included in our architecture to allow for suitable dynamic approaches, such as PID loops, often necessary in practical applications, where unknown external disturbances or unmodeled dynamics affect the dynamics of ω and may negatively affect closed-loop stability and performance.

Outer loop ( 7) The following proposition well characterizes the dynamics of the closed-loop system. Proposition 1. Consider the plant dynamics (1)-( 2) and the dynamic controller ( 5)-( 6). Using the tracking errors in (3)-( 4), the closed-loop dynamics reads Ṙe = R e S(γ R (R e )) + R e S(ω e ) (9)

R Inner loop (5)-(6) ω R d , ω d ω v , R e , R d , ω d , ωd τ c
J ωe = -γ ω (x c , ω e , ω v (R e ,t)) -τ e (10) ẋc = γ c (x c , ω e , ω v (R e ,t)). ( 11
)
Proof. Taking the time derivative of (3) along (1), and using the properties of ω d and R d in Assumption 1, we get Ṙe

= Ṙ d R + R d Ṙ = -S(ω d )R d R + R d RS(ω) = -S(ω d )R e + R e S(ω) = R e (S(ω) -R e S(ω d )R e ) (12)
Recall now that S -1 (RS(x)R ) = Rx ∀x ∈ R 3 , which implies S(Rx) = RS(x)R , and add and subtract ω v from the previous equation, to get, based on ( 7),

Ṙe = R e (S(ω v ) + S(ω) -S(ω v ) -R e S(ω d )R e ) = R e (S(γ R (R e )) + R e S(ω d )R e -S(ω e ) -R e S(ω d )R e )
= R e S(γ R (R e )) -R e S(ω e ) = R e S(γ R (R e )) + R e S(ω e ) .

(13) Finally, computing J ωe = J( ωvω) and substituting equations ( 2), ( 7) and ( 6) into (4), we get (10).

• Remark 2. As shown in Figure 2, the error dynamics ( 9)-( 11) presents a nontrivial feedback structure perturbed by the (timevarying) desired velocity t → ω d (t). It is emphasized that one could transform this feedback structure into an autonomous cascaded form by selecting γ ω (x c , ω e , ω v ) = γω (x c , ω e ) and γ c (x c , ω e , ω v ) = γc (x c , ω e ). While mathematically elegant, such a selection can result in limited degrees of freedom in the design of stabilizers that do not cancel non-harmful nonlinearities (for example only partially canceling the term -S(ω)Jω as done in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in thrust-vectoring UAVs[END_REF]). Such stabilizers are also used in our design of Section 5. Our general formulation results in the feedback interconnection of Figure 2: while being more challenging from the point of view of stability analysis, this non cascaded solution gives more flexibility to the control designer. Equations ( 10)-( 11) Equation ( 9)

ω d ω e R e
Fig. 2. Closed-loop error dynamics: feedback interconnection.

Stability analysis

In this section we address the stability analysis of the feedback interconnection structure in Figure 2 of the closed-loop dynamics ( 9)-( 11) characterized in the previous section. The main goal is to guarantee asymptotic stability of the equilibrium 1 (R e , ω e , x c ) = (I 3 , 0, 0), which captures the tracking requirements in Problem 1. Due to the presence of the external perturbation t → ω d (t) in the error system a possible strategy for the study of asymptotic stability of the time-varying error dynamics ( 9)-( 11) is to ensure that γ R in ( 7) is uniformly bounded over SO(3) and exploit the uniform boundedness of ω d in Assumption 1, to conclude that there exists a uniform bound ω M > 0 on t → ω v (t) along all solutions to (9)-( 11), and then embed these solutions in the larger funnel of solutions of the following differential inclusion:

Ṙe = R e S(γ R (R e )) + R e S(ω e ) (14) ωe ẋc ∈ F ω (ω e , x c ), ( 15
)
1 It is assumed hereafter that one can write the perturbed closed-loop error dynamics as ( 9)-( 11) for a given perturbation τ e with an equilibrium point having x c = 0, possibly after introducing a suitable change of coordinates (see for instance Example 1).

where

F ω (ω e , x c ) := co ω v ≤ω M -J -1 γ ω (x c , ω e , ω v ) γ c (x c , ω e , ω v ) , (16)
where co denotes the convex hull and where, to the end of studying stability, we select the input τ e = 0. It is instructive to observe that embedding the solution set of ( 9)-( 11) in the larger funnel of solutions to ( 14), ( 15) allows us to restore the cascaded structure obtained in Remark 2 when restricting the degrees of freedom of the control architecture. Remark 3. A basic requirement for the well posedness (in terms of existence and sequential compactness of the set of solutions) of the differential inclusion ( 14), ( 15) is that the map γ R be a continuous function and that F ω be outer semicontinuous and locally bounded relative to

R 3 × R n c , F ω (ω e , x c
), in addition, to being nonempty and convex for each (ω e , x c ) ∈ R 3 × R n c . These properties are automatically satisfied if γ c and γ ω in the original control law ( 5), ( 6) are continuous functions.

To analyze asymptotic stability of the equilibrium (R e , ω e , x c ) = (I 3 , 0, 0), we use the following result, which is a slightly modified version of [Angeli 2004, Theorem 2], adapted to the closed-loop system ( 9)-( 11) evolving on SO(3

) × R 3 × R n c .
This lemma is instrumental to prove our main result. Lemma 4. Let x ∈ R n and R ∈ SO(3) and consider the cascade

interconnection ẋ = f (x), Ṙ = Q(R) + G(R, x) where f : R n × R ≥0 → R n , Q : SO(3) → T SO(3) and G : SO(3) × R n → T SO(3) are smooth vector fields with f (0) = 0, Q(I 3 ) = 0 and G(I 3 , 0) = 0. Suppose that the equilibrium R = I 3 is almost global Input-to-State Stable (aISS) for Ṙ = Q(R) + G(R, d)
with respect to d, namely, that R = I 3 is locally asympotically stable for d = 0 and that there exists γ ∈ K such that for each essentially bounded and measurable d : R → R n , there exists a zero volume set

B d ⊂ SO(3) such that, for all R(t 0 ) ∈ SO(3) \ B d , lim sup t→+∞ R(t) SO(3) ≤ γ( d ∞ ). ( 17
)
If the equilibrium x = 0 is GAS for ẋ = f (x), then (x, R) = (I 3 , 0) is almost globally asymptotically stable. Remark 5. The notion of almost global ISS, originally introduced in [START_REF] Angeli | An almost global notion of input-to-state stability[END_REF], is exploited here to deal with the fact that, due to topological obstructions, there exists no continuous function γ R (•) globally asymptotically stabilizing the desired attitude. In the context of nonlinear observer design [START_REF] Vasconcelos | Combination of Lyapunov and density functions for stability of rotational motion[END_REF], this notion has been exploited to show that the widely adopted stabilizer γ R (R e ) := -k R 2 S -1 (skew(R e )), k R ∈ R >0 , makes the equilibrium R e = I 3 aISS for (9) with respect to ω e for a sufficiently large gain.

The proof of the Lemma 4 follows exactly the same steps as the proof given in [START_REF] Angeli | An almost global notion of input-to-state stability[END_REF] and is therefore omitted. Note that herein we assume the GAS property (instead of almost GAS) for the perturbing subsystem (a more desirable property for our purposes), which is not hard to obtain as both the angular velocity and the controller dynamics evolve in Euclidean spaces.

In the spirit of a modular control design [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust vtol uavs[END_REF], we now present some basic properties that the functions γ R (•), γ ω (•, •, •) and γ c (•, •, •) must satisfy in order to guarantee desirable stabilizing properties for the closed-loop system. Property 1. (Outer loop). The attitude stabilizer R e → γ R (R e ) is at least C 2 , uniformly bounded, and such that

• R e = I 3 is a locally asymptotically stable equilibrium point for Ṙe = R e S(γ R (R e )); • subsystem ( 9) is aISS with respect to ω e as in (17). Remark 6. We emphasize that the uniform boundedness assumption of γ R comes for free when γ R is the gradient of a smooth function on SO(3). Property 2. (Inner loop). Given a uniform bound γR > 0 on function γ R , and a uniform bound ωd > 0 from Assumption 1, define ω M := γR + ωd . The angular velocity error stabilizer (x c , ω e , ω v ) → γ ω (x c , ω e , ω v ) and the controller (x c , ω e , ω v ) → γ c (x c , ω e , ω v ) are continuous functions such that (ω e , x c ) = (0, 0) is a GAS equilibrium point for ( 15), ( 16).

Based on the above properties and Lemma 4, the main result of the paper is given by the following theorem. Theorem 1. Consider the closed-loop system described by ( 1)-(2) controlled by ( 5)-( 6). If γ R (•), γ ω (•, •, •) and γ c (•, •, •) satisfy Properties 1 and 2, then for any desired trajectory t → (R d (t), ω d (t)) ∈ SO(3) × R3 satisfying Assumption 1, the control law ( 7)-( 6) solves Problem 1. In particular, the point (R e , ω e , x c ) = (I 3 , 0, 0) is almost globally asymptotically stable.

Proof. From the definition of ω M in Property 2 and expression (7), we have that along all the closed-loop solutions it holds that |ω v (R e ,t)| ≤ γR + ωd = ω M . Hence, the solutions of the error dynamics ( 9)-( 11) are contained in the solution funnel of the differential inclusion ( 14)-( 15). The cascaded structure of ( 14)-( 15) comprises an upper subsystem (15), for which the equilibrium point (ω e , x c ) = (0, 0) is GAS due to Property 2, and a lower subsystem ( 14), for which the stabilizer γ R (•) guarantees local asymptotic stability of R e = I 3 and aISS due to Property 1. The proof is then completed by applying Lemma 4 to these two subsystems.

•

Before concluding this section, we emphasize that an extension of our main result (Theorem 1) to global asymptotic stability could be achieved by employing suitable hybrid attitude stabilizers [START_REF] Casau | Robust global exponential stabilization on the ndimensional sphere with applications to trajectory tracking for quadrotors[END_REF] guaranteeing ISS-like properties.

SAMPLE STABILIZERS DESIGN

In this section we present a sample selection of the basic components of our control design, namely functions γ R (•), γ ω (•, •, •) and γ c (•, •, •). Then, we prove that they satisfy Properties 1 and 2, thereby guaranteeing the applicability of Theorem 1.

For the outer loop controller (7), we employ a widely adopted solution for attitude stabilization [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF]]. As for the inner loop controller, we propose a structure involving an important linear dynamic component, and then use passivity arguments to satisfy Property 2. Our general selection includes robust control solutions such as PI control loops. As an important outcome of these results (and of Theorem 1) our hierarchical construction allows proving almost global tracking results in the presence of constant disturbances τ e in (2), without imposing conditions on the PI controller gains [START_REF] Maithripala | An intrinsic PID controller for mechanical systems on Lie groups[END_REF]Berg 2015, Goodarzi et al. 2013].

Outer loop stabilizer

The stabilization of the attitude kinematics in a global sense is a challenging problem due to the non-Euclidean nature of SO(3) [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF]]. Among the large number of (nonlinear) stabilizers that have been developed in the literature, we select one of the most common solutions, for which proving the aISS of Property 1 is straightforward, thanks to the tools developed in [START_REF] Angeli | Stability robustness in the presence of exponentially unstable isolated equilibria[END_REF].

Proposition 1. Given γ R (R e ) := -1 2 S -1 (skew(K R R e )), where K R ∈ R 3×3 is a symmetric matrix with distinct eigenvalues satisfying tr(K R )I 3 -K R ∈ R 3×3
>0 , Property 1 is satisfied.2 Proof. The proof stems from the results of [Angeli and Praly 2011, Proposition 2] combined with those of [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF]. As shown in the latter, γ R (R e ) makes the equilibrium R e = I 3 locally asymptotically stable for Ṙe = R e γ R (R e ) while the other equilibria R e ∈ Ω R := {R ∈ SO(3) \ I 3 : γ R (R) = 0} are isolated and exponentially unstable, i.e., at least one eigenvalue of the tangent map of γ R (R e ) restricted to R e ∈ Ω R has positive real part. Thus, the sufficient conditions A1-A2 of [Angeli and Praly 2011, Proposition 2] are matched. Finally, aISS holds upon verification that the solutions to ( 9) are ultimately bounded according to [Angeli and Praly 2011, Definition 2], which is trivial for our case where the state evolves on a compact manifold. By [Angeli and Praly 2011, Proposition 2], subsystem (9) with γ R (•) defined as in Proposition 1 is aISS with respect to ω e and thereby satisfies Property 1.

•

Inner loop stabilizer

The inner loop stabilizer comprises the stabilizer γ ω (•, •, •) for the angular velocity dynamics and the function γ c (•, •, •), which shapes the controller dynamics. Consider the selection

γ c (x c , ω e , ω v ) := A c x c + B c ω e (18) γ ω (x c , ω e , ω v ) := C c x c + D c ω e + K ω ω e + S(ω e )J(ω v -ω e ) (19)
where

A c ∈ R n c ×n c , B c ∈ R n c ×3 , C c ∈ R 3×n c , D c ∈ R 3×3 and K ω ∈ R 3×3 .
The control torque obtained by substituting ( 19) in (5), i.e., τ c = S(ω v )Jω + J ωv + C c x c + (D c + K ω )ω e , does not cancel the gyroscopic term -S(ω)Jω: this may be helpful in practical implementations to reduce the risk of actuators saturation, e.g., when the desired angular velocity ω d is much smaller than the initial angular velocity. It is emphasized that the dynamics in ( 18)-( 19) comprises some linear dynamics. The following proposition proves Property 2 whenever that linear dynamics is passive (a special case of this being PI controllers). Proposition 2. Consider the inner loop controller ( 5),( 6) with the selections ( 18)-( 19). Suppose that G ω (s

) := C c (sI n c - A c ) -1 B c + D c is a n c × n c positive real transfer function, where (A c , B c
) and (C c , A c ) are controllable and observable pairs, respectively. 3 Then, for any K ω ∈ R 3×3 >0 and any ω M > 0, Property 2 is satisfied.

Proof. When using selection ( 18)-( 19), the differential inclusion in (15) can be written as ωe ẋc ∈

ω v ≤ω M J -1 (S(ω e ) J(ω v -ω e )-C c x c -(D c +K ω )ω e ) A c x c +B c ω e =: F ω , (20)
indeed F ω is already convex because ω v enters linearly in the dynamics. To prove our proposition we rely on an invariance principle for differential inclusions similar to [Invernizzi et al. 2020, Proposition 3], based on the results in [START_REF] Seuret | A nonsmooth hybrid invariance principle applied to robust event-triggered design[END_REF]. Consider the following Lyapunov candidate V ω := 1 2 ω e Jω e + 1 2 x c Px c , where P ∈ R n c ×n c >0 . Next, define Vω (ω e , x c ) := max

f ∈F ω (ω e ,x c )
∇V (ω e , x c ), f . From (20), we obtain:

Vω (ω e , x c ) = ω e (-C c x c -(D c + K ω )ω e ) + 1 2 x c PA c + A c P x c + x c PB c ω e (21)
where we exploited the property ω e S(ω e ) J(ω vω e ) = -ω e (ω e × J(ω vω e )) = 0 ∀ω v . Since G ω (s) is positive real, from the positive real lemma [Khalil 2002, Lemma 6.2], there exist 21) can be written as:

P ∈ R n c ×n c >0 , L, W such that PA c + A c P = -L L, PB c = C c -L W and W W = D c + D c . Hence, equation (
Vω (ω e , x c ) = -ω e K ω ω e -ω e C c x c - 1 2 ω e W W ω e - 1 2 x c L Lx c + x c C c ω e -x c L W ω e (22) = -ω e K ω ω e - 1 2 (Lx c +W ω e ) (Lx c +W ω e ) ≤ 0, ∀(ω e , x c ) ∈ R 3 × R n c
. Combining this result and the fact that V ω is radially unbounded, compact sets of the form

Ω c := (ω e , x c ) ∈ R 3 × R n c : V ω (ω e , x c ) ≤ c , c ∈ R >0
, are forward invariant. Finally, since the viability condition of [Goebel et al. 2012, Pag. 124] is satisfied, solutions starting in Ω c are complete. Consider now any solution t → φ (t) := (ω e (t), x c (t)) starting in Ω c at t = t 0 ≥ 0 such that V ω (φ (t 0 )) = a = 0. Then, if φ (t 0 ) ∈ (ω e , x c ) ∈ R 3 × R n c : ω e = 0 , the function V ω (φ (t)) has to decrease in time by continuity and by the negative upper bound in (22). Instead, if φ (t 0 ) ∈ G :=

(ω e , x c ) ∈ R 3 × R n c : ω e = 0 \ {(0, 0)}, then, according to the closed-loop dynamics restricted to G,

ωe ẋc ∈ F ω (ω e , x c )| G = ω v ≤ω M -J -1 C c x c A c x c = -J -1 C c x c A c x c (23)
φ (t) will exit the set G for some small t ≥ t 0 . Indeed, if the solution were to evolve in G, due to linearity of the dynamics (23), x c (t) = exp(A c (tt 0 ))x c (t 0 ), and remaining in G would also impose J ωe (t) = C c exp(A c (tt 0 ))x c (t 0 ) = 0 for all t ≥ 0. Since (C c , A c ) is observable, this establishes a contradiction, thus proving all the items of the invariance principle and establishing Property 2.

Example 1. (PI-like controller.) Consider the dynamic part of the controller corresponding to an integral action, i.e., A c = 0,

B c = I 3 , C c = diag(k I 1 , k I 2 , k I 3 ) ∈ R 3×3 >0 , D c = 0.
Due to the decoupled structure for each axis, checking positive realness (as required by Proposition 2) is trivial because G ω (s) is diagonal and each element G ω ii (s) = k I i /s is positive real. The overall control law corresponds to a PI-like loop, comprising the proportional term K ω . This solution is particularly interesting as it allows rejecting a constant unknown torque τ e ∈ R 3 in (2). This can be shown by referring to the perturbed closed-loop dynamics J ωe = S(ω e ) J(ω vω e ) -C c x c -K ω ω e + τ e and by applying the change of coordinates xc = x c -C -1 c τ e : the closed-loop dynamics in the new coordinates has the same structure as equation (20). Hence, we can conclude that (ω e , x c ) = (0,C -1 c τ e ) is GAS from the proof of Proposition 2.

Remark 7. More complex inner loop controllers achieving desirable performance and disturbance rejection capabilities could be considered. For instance, PID control, typically employed in aircraft applications and multirotor UAV control [PX4-Community 2018], fits within the class of controllers outlined in Proposition 2 . Disturbance rejection of harmonic disturbances with known frequency but unknown phase and amplitude could be also considered, by including an internal model along the lines of [START_REF] Camblor | A Comparison of Nonlinear PI and PID Inertia-Free Spacecraft Attitude Control Laws[END_REF].

NUMERICAL RESULTS

We present a simulation example to show that the control law ( 5)-( 6) with the inner and outer loop stabilizers selected as in Propositions 1 and 2 solves the robust tracking problem in the presence of constant disturbances τ e in (2). We employ a PIlike controller for the inner loop corresponding to the selection

A c = 0, B c = I 3 , C c = diag(k I 1 , k I 2 , k I 3 ) ∈ R 3×3 >0 , D c = 0, K ω = diag(k ω 1 , k ω 2 , k ω 3 ) ∈ R 3×3
>0 , from which we expect asymptotic tracking in the considered scenario, following Example 1.

The desired attitude motion is characterized by a polynomial spin up maneuver with a steady-state angular velocity ω d ss = (0.5, 0.5, 0.5)rad/s (see Figure 3). The initial state that we select corresponds to a largely misplaced configuration R e (0) = I 3 + sin(π)S(e 1 ) + (1cos(π))S(e 1 ) 2 and to a significant angular velocity error e ω (0) := ω(0) -R e (0)ω d (0) = (3, 3, 3)rad/s with respect to the target velocity. The rigid body that we use in the numerical simulation has inertia matrix J = diag(1, 2, 3)kgm 2 and is affected by a piece-wise constant disturbance torque τ e = (1, 1, 1)Nm for t < 15s and τ e = (3, 3, 3)Nm for t ≥ 15s. The gains of the controller are tuned to have a predefined behavior of the closed-loop system in the proximity of the desired attitude motion. Specifically, the gains of the inner loop stabilizers G c i (s) = k ω i + k I i s are chosen to have the linearized inner loop error dynamics4 , namely, J ωe i = -k I i x ik ω i ω e i , ẋI i = ω e i , behaving like first order systems with bandwidth 10rad/s (k

ω 3 = 2k ω 2 = k ω 1 = k I 3 = 2k I 2 = 3k I 1 = 3.33).
The outer loop gains are tuned to have a 1rad/s bandwidth for the linearized outer loop error5 dynamics (K R = diag(1, 1.001, 0.999)). The attitude tracking performance of the proposed controller τ ) with bound M = 3Nm on each component.6 Finally, we also considered for comparison a PD-like controller τ pd c := S(R e ω d )JR e ω d + JR e ωd + γ R (R e ) + K ω (R e ω dω) (magenta dotted line) borrowed from the literature, tuned to have a similar transient response (K R = diag(25, 12.5, 0), K ω = diag(10, 20, 30)). From Figures 45, one sees that the feedback linearizing controller achieves a faster tracking at the expense of a larger actuation effort than the other PI-like controllers. The PD-like controller cannot compensate for the disturbance effects as shown by the non-null steady state error in Figure 4 and furthermore, it requires much more actuation effort than the other controllers while showing a similar transient (Figure 5). All the PI-like controllers promptly react to the abrupt change of the disturbance torque occurring at t = 15s. It is worth mentioning that the saturated controller is capable of asymptotically tracking the desired motion but has the worst transient performance.

CONCLUSIONS

In this paper we revisited the attitude tracking control problem and proposed an inner-outer loop design achieving almost global tracking results. By means of the proposed approach one can exploit dynamic controllers for the inner loop, such as PIDlike loops, to achieve desirable performance and disturbance rejection capabilities. Future work is oriented to obtaining global tracking results and to showing the potential of the proposed approach in handling model uncertainties and saturation effects.

Fig. 1 .

 1 Fig. 1. Proposed inner-outer loop dynamic control architecture.

  prop c := S(ω v )Jω + J ωv +C c x c + K ω ω e is illustrated in Figure 4 (black dotted line) where it is compared to the results obtained with a feedback linearizing version of the control law (red solid line), namely, τ fl c := S(ω)Jω + J ωv + C c x c + K ω ω e ) and with a saturated version (blue dash-dotted line) τ sat c := σ M (τ prop c

Fig. 3 .Fig. 4 .

 34 Fig. 3. Spin up motion: desired angular velocity ω d .

Fig. 5 .

 5 Fig. 5. Control input: infinity norm of the control torque τ c .

K R having distinct eigenvalues guarantees isolated solutions to γ R (R e ) = 0, a condition required to use the arguments in[START_REF] Angeli | Stability robustness in the presence of exponentially unstable isolated equilibria[END_REF] to prove aISS. The case K R = k R I

, k R ∈ R >0 , has been analyzed in[START_REF] Vasconcelos | Combination of Lyapunov and density functions for stability of rotational motion[END_REF], whose results can be used to extend the proof of Proposition 1 also for this gain selection.3 Note that assuming controllability and observability of (A c , B c ) and (C c , A c ) is without loss of generality, as a matter of fact any unobservable or uncontrollable subcomponents would not cause any effect on the closed loop, besides possibly perturbing the initial transient.

Note that we assumed ω d = const ≈ 0.

Considering small error angles θ e ∈ R 3 such that R e ≈ I 3 + S(θ e ), the linearized outer loop error dynamics is given by θe = -KR θ e , where KR := 1 2 diag k R 2 + k R 3 , k R 1 + k R 3 , k R 1 + k R 2 .Slightly different gains are used to respect the conditions of Proposition 1. See Footnote 2 about the use of equal gains in the attitude stabilizer γ R (•).

6 σ M (x) := (sat M (x 1 ), sat M (x 2 ), sat M (x 3 )) for x := (x 1 , x 2 , x 3 ), where sat M (•) := min(max(•, -M), M) for M ∈ R >0


Ruben Begnis
email: ruben.begnis@studenti.unitn.it

Giulio Panzani
email: giulio.panzani@polimi.it

Mirko Brentari

Luca

Luca Zaccarian
email: zaccarian@laas.fr

An LMI-based approach for the control of semi-active magnetorheological suspensions

Keywords: Semi-active suspensions, vehicle dynamics, magnetic suspension, ride comfort, passive suspensions

come

INTRODUCTION

Semi-active suspensions correspond to shock absorbers whose dynamical nature can be manipulated by way of some suitable control inputs with the constraint that no energy can be injected into the system (as opposed to active suspensions that can inject energy, see [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF]). These suspensions have become popular in the automotive field because they are cost-effective solutions that can be employed to considerably increase the driving comfort.

The first developed semi-active suspensions technology consisted in the so-called linear adaptive configuration, wherein the viscous damping coefficient c of an electrohydraulic (EH) shock absorber can be adjusted in real-time between a bounded (and positive) set of allowable ranges [c min , c max ], by manipulating the size of the orifices within the suspension. While the passive EH configuration is fully linear (because viscous damping is a linear phenomenon), the input c of the linear adaptive solution enters nonlinearly in the vehicle vertical dynamics, thereby posing interesting control design challenges for improved passenger comfort and/or road handling. Different approaches have been proposed for the control of such devices; among them, switching strategies [START_REF] Savaresi | Mixed sky-hook and add: Approaching the filtering limits of a semi-active suspension[END_REF], optimal control [START_REF] Poussot-Vassal | Optimal skyhook control for semi-active suspensions[END_REF], MPC [START_REF] Canale | Semiactive suspension control using "fast" model-predictive techniques[END_REF] and LPV approaches [START_REF] Poussot-Vassal | A new semi-active suspension control strategy through lpv technique[END_REF] can be mentioned.

A different semi-active suspensions technology is based on the magnetorheological (MR) effects of certain fluids [START_REF] Wereley | Magnetorheology: advances and applications[END_REF], usually oil. When subjected to a magnetic field, the oil exhibits an increased viscosity, to the point of becoming a viscoelastic solid: such an effect is employed to change the equivalent damping force exerted by the suspension, without the need of changing its geometry. The control input in this case is once again nonlinear and corresponds to the yield stress of the fluid, accurately adjusted by varying the magnetic field intensity with a simple electromagnet. The peculiarity of the magnetorheological technology lies in the highly non-linear characteristic of its damping force. Typically, it appears that the non-linear characteristic resembles the overall effect of a constant viscous term, plus a saturation-like term, whose amplitude is proportional to the control input.

Even in the passive configuration (namely when the control input is held constant), the magnetorheological suspension provides a nonlinear behaviour: as a consequence, controlling these devices is an even more complicated task. There are two main streams in the scientific literature: on the one side state-of-the-art solutions developed for the EH technology have been reused in the MR case, see e.g. [START_REF] Yao | Mr damper and its application for semi-active control of vehicle suspension system[END_REF]; [START_REF] Corno | Design and validation of a full body control semi-active suspension strategy for a supercar[END_REF]. On the other side, researchers also developed custom control solutions explicitly addressing the nonlinearity of the MR technology, see e.g. [START_REF] Du | Semi-active h∞ control of vehicle suspension with magneto-rheological dampers[END_REF]; [START_REF] Turnip | Control of a semi-active mr-damper suspension system: A new polynomial model[END_REF]; [START_REF] Zapateiro | Semiactive control methodologies for suspension control with magnetorheological dampers[END_REF].

In this paper we set in the latter mainstream, by proposing a model-based controller design for a semi-active MR suspension, to improve the passenger comfort. With respect to the available literature, the MR suspension is here modelled by exploiting the saturation nonlinear elements. As a result, the control problem is addressed using the Linear Matrix Inequalities (LMIs) framework stemming from regional sector properties of the saturation function, as summarized in [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF]. The proposed control strategy is preliminarily validated against a passive configuration and a Skyhook strategy, a consolidated switching approach for the control of EH semi-active suspensions.

The reminder of the paper is as follows: in Section 2 the vehicle and the MR suspension model are introduced, along with the road generation and the performance cost function. The benchmark strategies are presented in Section 3 and the proposed LMI-based solution is given in Section 4. The simulations in Section 5 show the effectiveness of the proposed controller and its characteristics. Concluding remarks are finally given in Section 6.

SYSTEM MODELLING AND PERFORMANCE METRIC

Quarter-car model

In order to test the proposed control strategy and compare it with known state-of-the-art solutions, a quarter car model is used, as typically proposed in the consolidated scientific literature (see e.g. [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF]).

The quarter car model parameters are tuned to match those of an SUV-like vehicle and are listed in Table 1. Indeed such vehicles benefit from the use of MR suspensions for the vertical dynamics control, given the increased ground clearance for the off-road usage, and the significant forces required at low speed to damp the high sprung mass.

The differential equations that describe the quarter car model follow:

m s zs = -kz -f d (z) m u zu = -k t (z u -z r ) + kz + f d ( ż), (1)
where z := z s -z u . In equations (1), the suspension is modelled as a linear spring, coupled with the damper whose force depends only on the suspension stroke speed ż, as further elaborated in the following.

MR semi-active suspension model

Modelling a semi-active magnetorheological suspension is a known topic in the scientific literature. There are two distinguishing features of such devices:

(1) A highly nonlinear force-speed relationship, made up by a pre-yield and a post-yield region: in the former the force increases drastically with the stroke speed whereas in the second one the force exerted by the device bends. Thanks to the semi-active capability, the pre-yield region can be increased or decreased. (2) A hysteretic behaviour highly dependent on the fluid properties, that can be more or less pronounced and which is usually neglected when dealing with control oriented models, see for example In the present work, the MR damper is modelled by focusing on the nonlinear force-speed relationship, neglecting the possible hysteretic behaviour of the device. This choice is consistent with the classical scientific literature approach, as discussed, and with the objective of this paper, namely highlighting the benefits of designing a force-speed shape aware control algorithm, as compared to the state-of-the-art solutions that are typically agnostic with respect to the device characteristics.

The upper plot of Figure 1 shows an example of the typical force-speed curves of a magnetorheological suspension, for different values of the applied current. It can be seen that when the pre-yield region is null, the suspension characteristic resembles a linear one (red dash-dotted line); for higher current values, the pre-yield region expands with a significant increase of the damping force. The mathematical model used in this work to describe the force-speed nonlinearity f d ( ż) of the MR damper is

f d = c min ( żs -żu ) + fd (u MR , żs -żu ).
(2) which is composed by two elements: a linear one, representing the minimum damping (the red dashed-dotted line in the upper plot of Figure 1) and the nonlinear term fd accounting for the nonlinearity of the suspension. In particular, the nonlinear term can be described by the saturation-like function shown in the bottom plot of Figure 1, modelled as:

fd (u MR , żs -żu ) = sat uMR (k 0 ( żs -żu )), u MR ∈ [0, fmax ].
(3)

where sat ū(s) := max{min{s, ū}, -ū} represents the symmetric scalar saturation function having limits ±ū, and the constant k 0 corresponds to the slope of the saturation function near the origin. Notice that in equation ( 3) the control variable of the semi-active suspension has been explicitly defined, and corresponds to the saturation level u MR .

The parameters that completely define equations ( 2) and

(3) are listed in Table 1: they have been chosen according to the results presented in [START_REF] Goldasz | Parametric study on the performance of automotive mr shock absorbers[END_REF], to represent an average automotive MR damper, that suites also for an SUV vehicle.

Performance characterization and road definition

The performances of the passive suspension and the semiactive control strategies are compared in terms of passenger perceived comfort. As widely accepted in the scientific literature, a means to quantitatively evaluate this aspect is the index J corresponding to the rms of the output acceleration y = zs :

J(z s ) = 1 T T 0 |z s (τ )| 2 dτ 1 2 , ( 4
)
where T is the duration of the simulation (or experiment).

The driving comfort is evaluated by simulating the quarter car model (1) for a road profile z r generated according to the ISO-8608 standard, as presented in [START_REF] Agostinacchio | The vibrations induced by surface irregularities in road pavements -a matlab R approach[END_REF]. All the simulations and the results of the next sections are achieved by considering a road profile of class C-D (representing a light off-road condition, see [START_REF] Agostinacchio | The vibrations induced by surface irregularities in road pavements -a matlab R approach[END_REF]), with a velocity of v = 90 km h . Figure 2 shows an example of the considered road profile, as a function of the longitudinal travelled distance. In order to avoid the variability due to the specific realization of the random Gaussian noise, the considered road has been selected to be sufficiently long (that is, 2.5km long) so that the functional J directly benefits of an averaging effect.

Overall system model

Combining all the elements of the simulation framework discussed above, it is possible to write a compact expression for the considered system. By defining the state variable vector x = [z s z żs ż] , considering as external input the excitation input z r coming from the road, and as the performance output the vertical acceleration y = zs , the overall system can be described by the following statespace equations:

Σ MR : ẋ = A(c min )x -B fd (u MR , ż) + Ez r y = zs = C(c min )x -D fd (u MR , ż) (5)
with

A(c min ) B E C(c min ) D =         0 0 1 0 0 0 0 0 0 1 0 0 0 - k ms 0 - c min ms 1 ms 0 kt mu - k ms - kt + k mu 0 - c min (ms + mu) mums ms + mu mums - kt mu 0 - k ms 0 - c min ms 1 ms         ,
The proposed model is generically nonlinear, due to the saturation that defines fd (u MR , ż); when u MR = 0 the model turns into a linear one, with minimum suspension damping, equal to c min .

PASSIVE AND SEMI-ACTIVE BENCHMARKS

The control strategy proposed in Section 4 is tested against two classical benchmarks when dealing with semi-active suspensions: the fully passive and the Skyhook semi-active control strategy.

Passive benchmark

The basic comparison of any semi-active control strategy is the corresponding passive suspension one. It is well-known that the choice of the suspension damping, in the passive case, undergoes a performances trade-off and the best comfort performances are found for intermediate (neither too low, nor too high) damping values. Thus, a sensitivity analysis of the comfort index with respect to different, constant, values of the control variable u MR has been performed: as a matter of fact, the best performances are obtained for the passive case (J ∼ = 1.6 m s 2 ) corresponding to setting u MR = 0. This result is consistent with the selected value of the linear damping c min which is very close to the optimal value that can be obtained for a linear suspension; within this perspective it is clear that introducing an additional MR damping force in the system (by using values of u MR > 0) worsens the comfort performances in the same way as it would happen for the linear case.

Skyhook benchmark

The Skyhook semi-active controller is a classical comfortoriented control strategy developed for semi-active suspensions, [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF]. In particular, the so called two-state Skyhook is a switching control strategy that switches the suspension damping between its minimum and maximum values, according to the following law:

c = c min , if żs ż ≤ 0 c max , if żs ż > 0. (6)
Despite being conceived and designed for linear semi-active suspensions (typically realised with the Electro-Hydraulic technology) this control strategy has been directly applied also to the MR suspensions, as recalled in the introduction.

In the present work, we consider this strategy as another benchmark to be compared with the LMI-based solutions of the next sections.

The MR Skyhook control strategy mimics the EH one, but rather than switching between two damping factors c min and c max , the control input u MR in (3) switches directly between 0 and fmax :

u MR = 0, if żs ż ≤ 0 fmax , if żs ż > 0 (7)
When simulating the quarter-car model with the control law (7) the resulting performance index is J ∼ = 2.35 m s 2 , even higher than the passive benchmark. On the one side this result is consistent with the Skyhook formulation which is designed to limit the chassis velocity, rather than its accelerations. However, such a performance worsening is amplified in the MR case, mainly due to the nonlinearities of the force-speed maps of these dampers.

The obtained result motivates the development of the MR control strategy proposed in the next section, which explicitly includes the shape of the force-speed maps in the controller design.

LMI-BASED CONTROLLER

The results of the previous section show that the twostates skyhook controller merely applied to the MR suspension is not effective at improving the passive performance. Motivated by this fact, in this section we design a specific controller for the MR suspension by exploiting Linear Matrix Inequality (LMI) techniques certifying stability and performance of systems subject to saturation nonlinearities [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF]. Since those techniques mostly focus on designing a saturated state feedback K in dynamical systems of type ẋ = Ax + B sat ū(K x), for a fixed value of ū, we perform below some manipulations of the MR dynamics (1), ( 2) to partially match this standard structure.

Approximate representation of fd

To the end of providing a representation compatible with existing LMI-based techniques for feedback control design, we consider the MR suspension output in (3), which enters linearly in the dynamics (5) (recall that c min is constant) and we provide the following approximate representation 8) it becomes then convenient to design the bounded equivalent input v MR as sat fmax 2 (K 2 x), where K 2 is a four element row vector representing a state feedback gain to be optimized.

fd = v MR + sat fmax 2 (k 0 ż), v MR ∈ -
In particular, as shown in Figure 3, the force fd generated by the actual suspension (3) driven by the input

u MR = fmax 2 + sgn( ż) sat fmax 2 ([k 1 k 2 k 3 k 4 ] x), = fmax 2 + sgn( ż) sat fmax 2 (K 2 x) vMR , (9)
perfectly coincides with the force exerted by the double saturation model ( 8) when |v MR | ≥ fmax 2 , namely whenever |K 2 x| ≥ fmax 2 . It turns out that the important feature ensured by the actual model (3) as compared to the approximated one (8) is to remove points from the graph of Figure 3 corresponding to injection of energy (namely the points in the second and fourth quadrant) so that the predicted performance J by the approximate model is not very different from the one of the actual model (this emerges from our numerical simulations).

LMI-based Controller design

Based on the approximated MR model (8) derived in the previous section and on the saturated state-feedback selection for v MR given in (9), we may write system (5), (8) as:

ẋ = A(c min )x -B sat fmax 2 (K 1 x) -B sat fmax 2 (K 2 x) + Ez r z = C(c min )x -D sat fmax 2 (K 1 x) -D sat fmax 2 (K 2 x),
(10) where matrix K 1 is fixed and provides the second contribution in (8), while matrix K 2 is to be designed and corresponds to the selection of v MR in (8). In particular,

K 1 K 2 = 0 0 0 k 0 k 1 k 2 k 3 k 4 ,
where we may exploit the degree of freedom of designing its second row K 2 . In the rest of this section, to simplify the notation, we will write A (respectively, C) in place of A(c min ) (respectively C(c min )), and sat in place of sat fmax because all the saturation limits correspond to the same value fmax 2 . We provide below two control design techniques, ensuring local asymptotic stability properties and suitable optimized disturbance rejection properties measured in terms of L 2 and H 2 performance.

When focusing on regional L 2 performance, we solve the following optimization problem, for a fixed value of s, which allows adjusting the size z r 2 of the maximum disturbance taken into account in the regional L 2 gain characterization.

γ 2 L (s) := min Q=Q >0,X,U1,U2,Y1,Y2,γ γ 2 subject to:

He

        (A -BK 1 )Q -BX BU 1 BU 2 E 0 K 1 Q + Y 1 -U 1 0 0 0 X + Y 2 0 -U 2 0 0 0 0 0 - I 2 0 (C -DK 1 )Q -DX DU 1 DU 2 0 -γ 2 I 2         < 0, Q Y 1 Y 1 f 2 max /(4s 2 ) ≥ 0, Q Y 2 Y 2 f 2 max /(4s 2 ) ≥ 0,
The following proposition establishes desirable properties when implementing K 2 = XQ -1 in (10), stemming from a solution to the optimization (11). Proposition 1. Given a scalar s > 0, consider any solution to the optimization (11). Then selecting K 2 = XQ -1 , the closed loop ( 10) is such that for any input z r having L 2 norm z r 2 ≤ s, it holds that the performance output satisfies z 2 ≤ γ L (s) z r 2 .

The second LMI-based design that we propose for the gain K 2 is based on a performance index of the H 2 type (see, e.g., (Boyd et al., 1994, §10.1.1 and page 137)), which is motivated by the fact that there is no feedthrough term between z r and z. In particular, we propose below to minimize the H 2 norm of the linear transfer function seen by the disturbance z r when the amplitude is small enough to not activate the deadzone functions. This norm can be minimized while, at the same time, ensuring asymptotic stability of the origin for the saturated closed loop, by solving the following optimization γ 2 H (s) := min Q=Q >0,X,U1,U2,Y1,Y2,γ γ 2 subject to: ( 12)

He (A -BK 1 )Q -BX BU 1 BU 2 K 1 Q + Y 1 -U 1 0 X + Y 2 0 -U 2 < 0, He (A -BK 1 )Q -BX 0 (C -DK 1 )Q -DX - 1 2 < 0, γ 2 E E Q > 0, Q Y 1 Y 1 f 2 max /(4s 2 ) ≥ 0, Q Y 2 Y 2 f 2 max /(4s 2 ) ≥ 0,
where the optimized variable γ is an upper bound on the maximum energy experienced at the output z 2 as compared to the size of the initial condition x 0 = ηE, selected in the subspace spanned by the input matrix E, namely z 2 ≤ γ|η|. In optimization ( 12), the positive parameter s establishes an indirect bound on the size of the ellipsoidal set of initial conditions that is guaranteed to be contained in the basin of attraction of the origin.

The following proposition establishes desirable properties when implementing K 2 = XQ -1 in (10), stemming from a solution to the optimization (12). Proposition 2. Given a scalar s > 0, consider any solution to the optimization (12). Then selecting K 2 = XQ -1 , the closed loop ( 10) is (locally) asymptotically stable to the origin with the set x Q -1 x ≤ s 2 being contained in the basin of attraction. Moreover, any (small signal) response remaining below the saturation limits for all times and starting in x 0 = ηE for some η ∈ R, satisfies

z 2 ≤ γ H (s)|η|
The proofs of Propositions 1 and 2, are based on the use of the regional sector conditions discussed in [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF] and are here omitted for space constraints.

SIMULATION RESULTS WITH LMI-BASED CONTROLLERS

Using the numerical parameters in Table 1, both the LMI optimization problems in ( 11) and ( 12) illustrated in Section 4 have been solved in Matlab, using the solver MOSEK over a range of possible values for the variable s between 0.1 and 1.5, which provides interesting trends for the quantity γ.

For each one of these solutions (corresponding to a specific selection of the four gains in K 2 ), the performance criterion J in (4) has been evaluated by simulation, by implementing two scenarios, corresponding to the regions represented in Figure 3:

• Double Saturation (DS) Model: represented by the bottom curve of Figure 3, corresponds to directly implementing the closed loop (10), namely the closed loop for which Propositions 1 and 2 provide stability and performance guarantees; • Variable Saturation (VS) Model: more faithful to reality, corresponds to only allowing pairs ( fd , ż) to take values in the actual region spanned by the feasible input-output pairs of the real suspension. This set of simulations has been performed by running model ( 5) with input fd obtained from (3) when the suspension input u MR is selected as in (9).

A second interesting result that can be observed from Figure 4 is that the green-black and red-blue curve pairs are very close to each other, thereby showing that the DS and VS saturation models considered in our validation Fig. 4. Evolution of the index J as a function of s for the L 2 -and the H 2 -based syntheses evaluated with the VS and the DS models. The horizontal purple line is the value of J obtained from the skyhook method of Section 3. provide essentially equivalent results, despite the presence of infeasible areas (those circled in red in Figure 3) when using the DS model. Such an evidence confirms that the approximation in Section 4.1 is reasonable and the controller design is effective even with this approximation.

To conclude, the simulation tests confirm that the proposed LMI-based techniques stemming from the approximated model of Section 4 are effective in providing a useful gain tuning technique for the semi-active MR suspension.

CONCLUSIONS

In this paper we addressed the control synthesis for a semiactive suspension designed based on magnetorheological effects. We first provided a performance characterization based on a prototype simulation scenario, then showed that the classical Skyhook control technique does not provide significant improvements as compared to the passive scenario (when the suspension input is kept constant), and we finally presented an approximated model allowing for optimized tuning of a state feedback gain, based on a semidefinite program written in terms of Linear Matrix Inequalities for efficient tuning. Simulation results confirmed the effectiveness of the proposed technique.

  [START_REF] Du | Semi-active h∞ control of vehicle suspension with magneto-rheological dampers[END_REF]; Do et al. (2012); Turnip et al. (2008); Zapateiro et al. (2011); Pepe et al. (2019).

Fig. 1 .

 1 Fig. 1. Typical force-speed curve associated to a MR suspension (upper plot) and the corresponding nonlinearity (lower plot).

Fig. 2 .

 2 Fig. 2. Example of a generated road profile.

  Fig. 3. Variable saturation model (upper plot) and double saturation model (lower plot) of the MR suspension. The red areas in the lower graph are regions where the double saturation model is not accurate. u MR ∈ [0, f max ] in (3). With model (8) it becomes then convenient to design the bounded equivalent input v MR as sat fmax

Figure 4

 4 Figure 4 reports the curves obtained from the evaluation of J in the various cases. A first result stemming from this figure is that the L 2 -based optimization does not provide interesting results because it gives large values of J. The H 2 -based scheme, instead, provides a minimum cost J * = 1.48 m s 2 obtained for s * = 0.1, which is desirably lower than the minimum cost achieved with the Skyhook linear controller (around 2.35 m s 2 ) and the passive solution (1.6 m s 2 ), leading to an overall comfort performance increase of ≈ 7%. The gain K 2 resulting from this optimal point corresponds to K 2 = [-18901 -45920 22704 -36338] . (13)