Valentin Preda
email: valentin.preda@esa.int

Francesco Sanfedino

Samir Bennani

Fabrice Boquet

Daniel Alazard

Robust and adaptable dynamic response reshaping of flexible structures

Keywords:

This paper outlines a complete methodology for designing a control system that reshapes the dynamic response of the flexible structure to robustly match the dynamics of a given adaptable reference model. The procedure was experimentally verified on a setup developed at the European Space Agency, consisting of a cantilevered flexible plate actuated by two shakers. Angular displacements at the free tip of the plate were measured with sub-microradian resolution using a laser autocollimator. Following a comprehensive system identification phase, mathematical models of the uncertain plant were extracted. The models reliably fit the experimental data and were used to synthesize a low order and high bandwidth structured Linear Parameter Varying controller. The controller was designed by taking into account the limits of achievable performance and the closed loop effectively constrained the flexible structure to behave like it was made out of an adaptable material. The robust stability and worst case performances were assessed by means of a structured singular value analysis and excellent agreement between theoretical predictions and experimental results was observed.

Introduction

Background and motivation

The fairing size envelopes of current launchers impose substantial constraints on spacecraft design. To overcome these size restrictions, a significant number of modern telecommunication, Earth observation, and science missions rely on large deployable structures to meet their performance goals [START_REF] Puiga | RandoA review on large deployable structures for astrophysics missionsActa Astronaut[END_REF][START_REF] Gonaloc | JulianDeployable structures activities at the European space agency's structures sectionECSSMET[END_REF]. There is currently a strong interest in developing ultralight, compactly packable [START_REF] Royers | PellegrinoUltralight ladder-type coilable space structures[END_REF][START_REF] Aryan | PellegrinoUltralight structures for space solar power satellites3rd AIAA Spacecraft Structures Conference[END_REF] structures that can be deployed [START_REF] Rodrigues | AhnDeployable soft composite structuresSci[END_REF] or even manufactured [START_REF] Pratern | Snyder3D Printing in Zero G Technology Demonstration Mission: complete experimental results and summary of related material modeling effortsInt[END_REF] and assembled in space. However, as structures become larger and more flexible they become more susceptible to mechanical vibrations due to the low structural damping in the materials and the absence of other forms of damping in the vacuum of space. Moreover, both ESA and other space agencies are preparing missions that require extremely high pointing accuracy and an ultra-quiet environment for the vibration-sensitive instruments [START_REF] Dennehyo | SalazarSpacecraft Micro-vibration: A Survey of Problems, Experiences, Potential Solutions, and Some Lessons Learned[END_REF][START_REF] Komatsuh | UchidaMicrovibration in spacecraftMech[END_REF]. This is a very challenging problem since a multitude of critical on-board equipment such as reaction wheels, cryocoolers, solar array drives and antenna pointing mechanisms generate significant levels of mechanical vibrations during operation [START_REF] Liux | LiRecent advances in micro-vibration isolationMech[END_REF]. These vibrations can easily propagate and interact with the large number of lightly damped and closely spaced flexible modes of the spacecraft structure [START_REF]ECSS-E-HH-32-26ASpacecraft Mechanical Loads Analysis Handbook[END_REF]. With the advent of increasingly powerful imaging sensors and more lightweight materials, understanding and controlling the propagation of these microvibrations and their dynamic interaction is one of the most critically important problem areas in spacecraft systems engineering [START_REF] Dennehyo | SalazarSpacecraft Micro-vibration: A Survey of Problems, Experiences, Potential Solutions, and Some Lessons Learned[END_REF].

Over the years, a wide variety [START_REF] Liux | LiRecent advances in micro-vibration isolationMech[END_REF] of microvibration isolation systems have been developed or proposed in an effort to tackle this problem. Most of them focus on mechanical isolation of the noisy equipment and/or of the sensitive payload from the rest of the spacecraft structure [START_REF] Davisj | RodenHubble Space Telescope Reaction Wheel Assembly Vibration Isolation System[END_REF][START_REF] Davisd | HydeSecond-generation hybrid D-strutC.D.JohnsonProceedings of SPIE -The International Society for Optical Engineeringvol[END_REF][START_REF] Cameloa | BrennanDamping and isolation concepts for vibration suppression and pointing performance50th[END_REF][START_REF] Malyd | PargettVibration suppression for the gemini planet imagerSPIE Astronomical Telescopes+ Instrumentation[END_REF][START_REF] Leeg | HanHybrid isolation of micro vibrations induced by reaction wheelsJ[END_REF][START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF][START_REF]PredaRobust Microvibration Control and Worst-Case Analysis for High Pointing Stability Space MissionsPh[END_REF]. For some observation missions, active and adaptive optics [START_REF] Bastaitsd | PreumontSegment ed bimorph mirrors for adaptive optics: segment design and experimentAppl[END_REF][START_REF] Rauschs | WittrockUnimor ph deformable mirror for space telescopes: design and manufacturingOpt[END_REF][START_REF] Bowlerimaging | Extrasolar Giant Planets[END_REF] are a promising technology to perform wavefront correction and improve pointing stability at the payload level. Another complementary strategy is to actively control the spacecraft structure itself and limit its ability to propagate microvibrations. These so-called smart space structures have been investigated for a number of years with promising results [START_REF] Huangb | WuExperimental study and numerical simulation of active vibration control of a highly flexible beam using piezoelectric intelligent materialAero[END_REF][START_REF] Sparksj | JuangSurvey of experiments and experimental facilities for control of flexible structuresJ[END_REF][START_REF]BalasTrends in large space structure control theory: fondest hopes, wildest dreamsIEEE[END_REF][START_REF] Hydeactive | Vibration Isolation for Precision Space StructuresPh[END_REF][START_REF] Balasj | DoyleRobustness and performance trade-offs in control design for flexible structuresIEEE Trans[END_REF]. However, most of these control strategies focus solely on vibration isolation and operate under the assumption that performance requirements and flexible plant dynamics are static.

In this context, the goal of the paper goes beyond microvibration isolation. The main purpose is to design and experimentally validate a control system that actively modifies the dynamic response of a flexible structure to match the dynamics of an adaptable reference model. The objective is to facilitate the development of future lightweight, actively controlled adaptable space structures. Adaptability is a key point in this study, since space structures can undergo significant changes during operation. These alternations in the structural dynamics can be induced by thermal gradients or changes in the inertial properties due to equipment realignment or fuel consumption. Furthermore, the aim is to present how the proposed closed loop system can be driven to the limits of the achievable performance and how guaranteed performance and stability certificates can be established without relying on extensive Monte-Carlo validation campaigns.

Contributions and paper organization

In agreement with the previously stated goals, the aim of the paper is to present the complete controller design cycle for an active structure experiment developed at ESA. The paper introduces the following key contributions:

• a complete physical model of the uncertain plant is developed. This model captures the various subtle interactions between the flexible structure and the actuators and provides deep insight into the underlying dynamics and the sensitivity to changes in various physical parameters. This model is calibrated using a grey-box identification procedure based on non-smooth optimization techniques. The resulting model is able to explain the structural modes observed experimentally all the way up to about 100 Hz.

• a model reduction method suitable for the reduction of high order Linear Fractional Transformation (LFT) models of uncertain or parameter dependent mechanical systems is presented. This reduction method is critical to ensure that the resulting loworder plant models are simplified enough in terms of numerical complexity to allow the usage of modern analysis and controller synthesis tools.

• a detailed controller synthesis and analysis procedure for the design, implementation and experimental validation of a low order LPV controller for guaranteed adaptable model matching of a flexible structure. The resulting controller is discretized at 2 kHz and ensures that the dynamics of the plant tracks an adaptable reference model in a bandwidth of about 80 Hz. This paper is organized into three parts: system modeling, controller design and performance analysis. In the first part (section 2), the experimental setup is presented and mathematical models of the uncertain system dynamics are extracted. The first model relies purely on the input-output experimental data while the second model is derived from the physical equations of motion. A grey-box identification procedure is introduced to calibrate the physical model based on the experimental data. The models are used in the second part of the paper (section 3) to design and optimize a controller capable of reshaping the response of the setup and match the dynamics of an adaptable reference model. The design process is systematically outlined together with the constraints imposed on the feedback loop by the different requirements. Lastly, the third part of the paper (section 4) details the rigorous analysis procedure that was used to obtain robust performance and robust stability certificates prior to the actual controller implementation. The section also presents a series of experiments that show excellent agreement with the theoretical predictions and further demonstrate the applicability of the proposed methods. gain. The k-th element of the vector signal u is denoted as u{k}. 𝜎 (A) denotes the maximum singular value of the real matrix A.

Notations

A right-handed coordinate frame centered at the point  with basis vectors ⃗ x, ⃗ y, ⃗ z is denoted as

(  ; ⃗ x, ⃗ y, ⃗ z ) .

Problem formulation and system modeling

Experimental setup and control objectives

Fig. 1 shows a photo and a schematic diagram of the experimental platform used in this study. The setup is mounted on an optical table and consists of autocollimator together with an aluminum cantilevered flexible plate with two proof-mass actuators (PMAs) and a mirror attached to it. The Newport CONEX-LDS autocollimator sends a laser beam that reflects of the mirror attached at the free tip of the plate and returns on a position sensor within the device. The angular deflection at the tip of the plate, corresponding to the angle between the outgoing and incoming ray, is subsequently computed down to sub 𝜇rad resolution with a sampling frequency of 2 kHz. This angular measurement y k (t) ∈ ℝ is then sent to a discrete controller implemented on a dSPACE MicroLabBox platform that computes the control voltages u k (t) ∈ ℝ 2 with a 2 kHz frequency. These signals are combined with the reference signals r(t) ∈ ℝ 2 to form u(t) ∈ ℝ 2 . The voltages u are then used to drive the two Wilcoxon F5B shakers through a set of two Kepco BOP-100 amplifiers. The purpose of the controller is to alter the dynamic response from the reference inputs r towards the angular measurements y k in order to robustly match the response of an adaptable reference model within a given control bandwidth. The aim of the setup is to provide a simplified model of a flexible space structure and demonstrate how the proposed design methodology can be employed for more complicated assemblies with more sensors, actuators and flexible modes. While the paper relies on this specific setup for experimental purposes, the control design process can be generalized to other space application scenarios. For example, the laser autocollimator can be replaced by a set of accelerometers or angular rate sensors while the proof-mass actuators can be replaced with piezoelectric patches. For an overview of actuators for space applications see Ref. [START_REF] Allegranzal | TosoActuators for space applications: state of the art and new technologiesProc[END_REF].

Model identification and mathematical modeling

The techniques presented in this paper rely on a good understanding of the main dynamics of the system in order to guarantee worst-case behavior. In essence, to design a controller that pushes the system to the limits of performance, it is critical to first develop a system model that includes the various perturbations and uncertainties acting on the plant. This is especially important in the space industry where systems need to be designed to work without maintenance for extended periods of time and withstand different structural changes induced by thermal deformations, gyroscopic effects or equipment realignment. In this section, two models of the ESA setup are introduced: an experimental black-box model and an analytical or symbolic model based upon the physical equations of motion. The experimental model is deduced purely based on the system response to various excitations. On its own, this empirical model can be used for controller design. However, this black-box representation provides no physical understanding of the dynamics of the system. As such, it is difficult to analyze and predict the changes in the overall dynamics as a result of variations in different physical parameters. To overcome such shortcomings, a second model was developed based upon the key physical principles and equations of motions. This model complemented the experimental one and provided deep insight into the sensitivity of the plant dynamics with respect to changes in different physical parameters. A procedure is also introduced to calibrate this analytical model and perform grey-box identification based on the experimental response.

Experimental model identification

Considering the plant schematic shown in Fig. 1b, the black-box model identification together with the uncertainty and noise characterization was performed in the following way: (a) In the first part, labeled with a ○ in Fig. 2a, the inputs u were kept at zero for 2 19 samples (≈ 4.37 min) in order to record in the y k deflection measurement the combined effect of sensor noise and environmental perturbations acting on the flexible plate. The estimate Φ nn (𝜔) of the power spectral density (PSD) of this signal was obtained at frequencies 𝜔 ∈ [7, 160] Hz using Welch's method with a Hann window of length 2 19 and 50% overlap. Fig. 2b shows the estimated peak amplitude spectral density1 (ASD) spectrum Φ nn (𝜔) 1/2 and cumulative root-mean square (CRMS), obtained during this phase across all experiments. The CRMS function provides a measure of the power in a signal, up to a given frequency 𝜔. Note that the ≈ 6.55 𝜇rad bias, occurring in the low frequencies below 0.1 Hz, is due to the static misalignment between the autocollimator and the mirror attached to the plate. Above this frequency, it can be seen that most of the power is concentrated around two key regions. Firstly, an increase of ≈ 0.53 𝜇rad occurs around 10.2 Hz and corresponds to deflections around the first bending mode of the plate due to unavoidable air and ground vibrations. This bending mode will be described in detail in the subsequent section. A second significant increase of ≈ 0.04 𝜇rad occurs around 50 Hz and is due to the electrical noise in the power supply of the actuator amplifiers.

For simplicity, the experimental model lumps all the different physical sources of perturbation observed in the laser measurement y k into a single output noise model. The measurement y k is therefore assumed to be the sum of a "noisefree" laser displacement response y and an overall noise signal n, i.e.

y k = y + n and n = W n d n with W n ∈ ℝℍ ∞ ; | | W n (j𝜔) | | ≥ Φ nn (𝜔) 1∕2 (1)
where d n is a zero mean unit variance white noise signal and the weighting filter W n is an upper bound on the ASD spectrum Φ nn (𝜔) 1/2 of the noise measurements.

(b) In the second and third phase (labeled b ○ and c ○ in Fig. 2a), each of the actuators was driven in turn by a finite energy stochastic signal. This random signal was obtained by passing zero mean unit variance white noise with a truncated normal distribution through a 12th order band-pass Butterworth filter with a pass-band of 5 Hz-200 Hz. The filter is scaled such that its  2 system norm (i.e. the variance of its output in response to unit white noise) is equal to the squared RMS value desired for each of the actuator signals. The length of each excitation sequence is equal to 2 21 samples for a total duration of ≈ 17.5 min. At the beginning of each sequence the input signal was slowly scaled up to its nominal value to avoid interactions with high frequency modes outside of the identification bandwidth. A similar scale down was performed at the end of each excitation sequence. Additionally a pause was inserted between these successive excitation phases to allow a decay of the flexible plate back to its resting state.

2. For each identification sequence, a Welch PSD estimate Φ y k y k (𝜔) of the output y k was computed using a Hann window of length 2 19 and 50% overlap. This provided a frequency resolution of about 8 mHz in the spectral estimate, averaged across 8 windowed intervals. Using the same method, the PSD estimates Φ u{1}u{1} (𝜔) and Φ u{2}u{2} (𝜔) were computed for each of the

G y k u{•} (𝜔) = Φ y k u{•} (𝜔) Φ u{•}u{•} (𝜔)
and

Γ y k u{•} (𝜔) = | | | Φ y k u{•} (𝜔) | | | 2 Φ u{•}u{•} (𝜔)Φ y k y k (𝜔) (2)
The coherence function quantifies at each frequency 𝜔, the fraction of the output PSD resulting from the input. The same function can also be seen as a measure of the causality between the inputs and the output response [START_REF][END_REF]. Nonlinearities, measurement noise or unwanted perturbations contribute to a reduction in the coherence spectrum. Hence, the function provides a strong indicator of the degree of uncertainty in the spectral estimate. 4. For each channel, a transfer function was fitted to the FRF results from the three experiments. The fit was performed with a vector fitting procedure [START_REF] Ozdemirs | GumussoyTransfer function estimation in system identification toolbox via vector fittingIFAC-PapersOnLine[END_REF] using MATLAB's tfest command. The normalized root mean squared error (NRMSE), measuring how well the response of the model fits the estimation data, averaged around 98%. After fitting each of the channels, the transfer functions were aggregated into a global system transfer matrix and a balanced order reduction was performed. An alternative identification method, explored in this study, was to directly estimate, based on both FRFs, the overall system transfer matrix using subspace methods and nonlinear least-squares fitting (see ssest command in MATLAB). However, this alternative method proved to be more time consuming and the fit quality was not significantly different. Finally, the resulting state space system G yu ∈ ℝℍ 1×2

∞ was put into the modal form

G yu ≔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ̇xG = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A 1 0 0 A 2 ⋱ ⋱ ⋱ 0 0 A n ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ A G x G + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ B 1 B 2 ⋮ B n ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏟⏟ B G u y = [ C 1 C 2 … C n ] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ C G x G
where

A i = [ Re(𝜆 i ) Im(𝜆 i ) -Im(𝜆 i ) Re(𝜆 i ) ] i = 1, … , n (3)
for complex conjugate eigenvalues 𝜆 i = Re(𝜆 i ) ± j Im(𝜆 i ) and A i = 𝜆 i for real eigenvalues. The state transformation to modal form was performed using MATLAB's canon command that computes a block-diagonal Schur factorization of the state matrix.

Fig. 3 shows the FRFs, minimum coherence spectrum and the resulting gains of the 14th order state space model G yu resulting from the previous identification procedure. It can be seen that the FRF estimates do not significantly vary across the three experiments and that the state space model reliably fits the response data in the bandwidth of interest. Furthermore, the coherence functions remains close to one except for frequencies near the anti-resonances and the 50 Hz power supply noise. This indicates that reliable FRF estimates were obtained using the experimental data.

The paper aims to demonstrate how the proposed control method can be applied even in the presence of significant model uncertainty. Linear Fractional Transformations (LFTs) are the one of the most widely employed means of representing uncertain or varying parameters and other nonlinearities [START_REF] Doylea | of LFTs, LMIs, and muConference on Decision and Control[END_REF] since any rational function can be expressed as an LFT [START_REF] Heckera | VargaSymbolic manipulation techniques for low order LFT-based parametric uncertainty modellingInt[END_REF]. Furthermore, interconnections of LFTs retain the LFT structure and therefore uncertainty at subsystem level can be aggregated into uncertainty at global system level. In the case of the plant considered in this paper, possible parameter variations and model inaccuracies are considered by augmenting the previously identified nominal model with an uncertainty structure. Three kinds of uncertainties are assumed to operate simultaneously:

1. Modal uncertainty. Variations in some structural parameters can lead to changes in the natural frequency and damping of some modes. To take into account these modal uncertainties, the blocks A i in (3) are replaced by

 i = [ (1 + r R i 𝛿 R i )Re(𝜆 i ) (1 + r I i 𝛿 I i )Im(𝜆 i ) -(1 + r I i 𝛿 I i )Im(𝜆 i ) (1 + r R i 𝛿 R i )Re(𝜆 i ) ]
for complex eigenvalues and

 i = [(1 + r R i 𝛿 R i )𝜆 i ] (4)
for real eigenvalues 𝜆 i . The parameters r R i , r I i are used to set the maximum percent of variation for the real and imaginary parts of each eigenvalue while 𝛿 R i , 𝛿 I i ∈ [-1, 1] are normalized real uncertainties. In this case, the new uncertain system matrix Â, replacing the nominal one in (3), is affine in 𝛿 R i , 𝛿 I i and can be expressed as the LFT:

 G = A G + W mL 𝚫 mod W mR = 𝚫 mod ⋆ [ 0 W mR W mL A G ]
with

𝚫 mod = diag ( 𝚫 mod 1 , … , 𝚫 mod n ) ⊂ ℝ n mod ×n mod and
𝚫 mod i i = 1, … , n = ⎧ ⎪ ⎨ ⎪ ⎩ [ 𝛿 R i I 2 𝛿 I i I 2 ] for complex eigenvalues 𝜆 i 𝛿 R i for real eigenvalues 𝜆 i ( 5
)
where the matrices W mL , W mR containing the scaling factors r R i , r I i can be computed by means of a singular value decomposition (see Ref. [START_REF] Heckera | VargaSymbolic manipulation techniques for low order LFT-based parametric uncertainty modellingInt[END_REF] for details). 2. Additive uncertainty. Outside the identification bandwidth or around the frequencies where the coherence spectrum ( 2) is low, the dynamics of the system is unknown or uncertain. These inaccuracies in the nominal system G sys were covered using the following additive uncertainty model:

Ĝyu = G yu + W add 𝚫 add where 𝚫 add = diag ( 𝛿 add 1 , 𝛿 add 2 ) and W add = [ r add 1 r add 2 ] (6)
The Linear Time Invariant (LTI) weights r add • ∈ ℝℍ ∞ are used to scale at different frequencies the magnitude of the additive

normalized LTI uncertainties 𝛿 add • with 𝜎 ( 𝛿 add • ) ≤ 1.
3. Multiplicative uncertainty. Neglected dynamics or gains fluctuations in the actuators and amplifiers were modeled as multiplicative uncertainties at the plant input. In this case, the new uncertain control signal û is equal to

û{•} = ( 1 + 𝛿 mul• r mul• ) u{•} or û = ( 𝚫 mul ⋆ [ 0 2 W mul I 2 I 2 ])
u with

𝚫 mul = diag ( 𝛿 mul 1 , 𝛿 mul 2 ) W mul = diag ( r mul 1 , r mul 2 ) (7)
where 𝛿 mul • are scalar normalized LTI uncertainties satisfying 𝜎

( 𝛿 mul • ) ≤ 1
and the magnitudes of the LTI weights r dsk • ∈ ℝℍ ∞ quantify for different frequencies the maximum percent of relative uncertainty.

The combined effect of all these uncertainties on the input-output behavior of the system can be studied by aggregating them into the following global uncertainty model:

Analytical modeling

In order to extract a symbolic representation of the dynamics, the plant is modeled as the interconnection of a finite element model (FEM) of the flexible cantilevered plate, two proof-mass actuators and the mirror load. Each of these subsystems is described below and afterwards assembled into a global structure. The numerical values and range of variation of the various system parameters that are used in this section are give in Table 1. These values were either measured directly, provided in the various datasheets or identified based on the experimental data using the techniques given in section 2.2.3.

(  • ; ⃗ x a• , ⃗ y a• ,⃗ z a• ) [ 114.04 114.04 ] mg m 2 a • gain 1.5 N V -1 k • stiffness 26 N m -1 c • damping 10 N s m -1 r  •  •
attachment point location relative to outer casing CoM, i.e. ⃖⃖⃖⃖⃖⃖⃖⃖⃖ ⃗

 •  • [ 0 0 -1.35 ] ⊺ cm Mirror m  mass 24.3 g J  moment of inertia in CoM frame (  ; ⃗ x l , ⃗ y l ,⃗ z l ) [
(  ; ⃗ x o , ⃗ y o ,⃗ z o )
when the plate is at rest.

1. Flexible plate model. The flexible plate was modeled as a Kirchhoff thin plate subdivided into the assembly of interconnected four-node rectangular plate elements shown in Fig. 5a. Each rectangle element has a side length of 5 mm at rest. The inertial reference frame

(  ; ⃗ x o , ⃗ y o ,⃗ z o )
, visible in Figs. 1 and5a, is fixed to the node at the center of edge anchored to the base structure. Attached to each of the unclamped nodes is a frame

(  • ; ⃗ x  • , ⃗ y  • ,⃗ z  • )
coincident in the rest state with the inertial frame. Each node has three degrees of freedom

q  • = [ z  • 𝜽  • 𝝍  • ] ⊺
where z  • is the displacement along the ⃗ z  • axis, 𝜽  • the rotation angle around the ⃗ x  • axis and 𝝍  • the rotation angle around the ⃗ y  • axis. Rotations around the ⃗ z  • axis are not considered. After grouping the node displacements into an overall displacement vector q  , the linearized dynamics of the flexible plate becomes

M  q + C  q + K  q  = f  with C  = 𝛽 M M  + 𝛽 K K  and q  (t) ∈ ℝ 1512 (9)
where f  are the generalized forces acting on every node coordinate of the plate. The state q  has a dimension of 1512 since the plate is discretized into 9 × 57 nodes, each with 3 coordinates and 9 of the nodes are in clamped condition. The mass and stiffness matrices M  and K  depend on structural properties of the plate given in Table 1. In-depth details on the construction of the matrices and Kirchhoff plate modeling can be found in Refs. [START_REF]CookConcepts and Applications of Finite Element Analysis[END_REF][START_REF] Sanfedinod | BoquetFinite element based N-Port model for preliminary design of multibody systemsJ[END_REF]. Since damping can be a challenging phenomenon to model physically, the classical (or Rayleigh) damping matrix C  is defined in terms of two uncertain real scalars 𝛽 M and 𝛽 K . The model assumes that each of the two PMAs both as well as the mirror load connect to the flexible plate attachment nodes labeled with

(  • ; ⃗ x p• , ⃗ y p• ,⃗ z p• ) and (
 ; ⃗ x e , ⃗ y e ,⃗ z e )
in Fig. 5a. The plate's equations of motion given in ( 9) can be rewritten in terms of input forces and output accelerations at the attachment nodes as

[ q⊺  q⊺  1 q⊺  2 ] ⊺ = D  [ f ⊺  f ⊺  1 f ⊺  2 ] ⊺ with D  = 𝚫  ⋆ D ∈ ℝℍ 9×9 ∞ ( 10
)
where the uncertain block diagonal matrix 𝚫  isolates the uncertain part of the damping coefficients 𝛼 M and 𝛼 K .

2. Mirror load model. The mirror used to reflect the laser beam is fixed to the plate using an aluminum L-shape (see Fig. 1), that significantly increases the stiffness of the plate region directly underneath. For simplicity, the mirror, together with the L-shape and the whole width of the plate under it, are modeled as a single rigid body load of mass m  . This combined static load connects to the plate model at the frame

(  ; ⃗ x l , ⃗ y l ,⃗ z l )
. Multiplying these body frame accelerations by the dynamic model of the load results in generalized reaction forces f  in the load body frame. These forces were subsequently translated to reaction forces f  in the attachment node frame in order to couple the plate and load dynamics. The linearized dynamical equations can therefore be expressed as

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ [ q f  ] = [  (r  ) 0 0  (r  ) ⊺ ] [ q f  ] = 𝚫  L ⋆ Ť [ q f  ] f  = [ m  0 0 J  ] q = 𝚫  ⋆ D q with  (r  ) = S ⊺  [ I 3 [ r  ] × 3 I 3 ] S  S  = [ 0 3×2 I 3 0 3×1 ] ⊺ [ r  ] × = ⎡ ⎢ ⎢ ⎢ ⎣ 0 -r  {3} r  {2} r  {3} 0 -r  {1} -r  {2} r  {1} 0 ⎤ ⎥ ⎥ ⎥ ⎦ (11)
where r  ∈ ℝ 3 are the coordinates of the attachment point  relative to the mirror CoM , i.e. ⃖⃖⃖⃖⃗  expressed in the inertial frame

(  ; ⃗ x o , ⃗ y o ,⃗ z o )
. m  and J  are the mass and the 2 × 2 moment inertia tensor of the mirror in the body frame. The matrix J  is only 2 × 2 since the load body frame (

 ; ⃗ x l , ⃗ y l ,⃗ z l
) is always coincident with the attachment node frame 3. Mirror/laser misalignment model. The noise free laser measurement y was considered to be almost equal to the pitch angle 𝜽  of the mirror around its ⃗ x  axis. Because of small deformations in the mirror surface as well as the inevitable misalignment between the laser beam and the reflecting mirror, the measurement y is also influenced by the axial displacement z  and the roll angle 𝝍  of the mirror load. In order to take into account these effects, the following uncertain model was introduced [START_REF] Davisd | HydeSecond-generation hybrid D-strutC.D.JohnsonProceedings of SPIE -The International Society for Optical Engineeringvol[END_REF] where 𝛿 𝜓  , 𝛿 z  ∈ [-1, 1] are real parametric uncertainties and r 𝜓  , r z  ∈ ℝ are used to set the maximum degree of expected coupling. 4. Proof-mass actuator model. Fig. 5b illustrates the model of the first PMA. This actuator was modeled as a hollow cylinder of mass m c1 and enclosing a point mass m 1 . This small mass attaches to the outside casing with a viscoelastic connection of stiffness k 1 and damping c 1 . The body frame

(  1 ; ⃗ x a1 , ⃗ y a1 ,⃗ z a1
) is fixed to the center-of-mass of the outside shell. Additionally, this frame forms a rigid connection with the plate frame

(  1 ; ⃗ x p1 , ⃗ y p1 ,⃗ z p1
) corresponding to node directly underneath the actuator (see Fig. 5a for clarity). Since rotations around ⃖⃖⃖⃗ z p1 axis are ignored in flexible plate model and

(  1 ; ⃗ x p1 , ⃗ y p1 ,⃗ z p1
) remains coincident with

(  1 ; ⃗ x a1 , ⃗ y a1 ,⃗ z a1
) , the moment of inertia tensor of the shell in the body frame J c1 is only 2 × 2. The control signal u{1} generates a magnetic force of magnitude a 1 u{1} in the voice coil that is applied in opposite directions to both the shell and the small mass. In the case when the actuator is not clamped to a supporting structure, this force produces a displacement z 1 of the small mass and a displacement z c1 of the cylindrical shell casing in an inertial reference frame coincident with

(  1 ; ⃗ x a1 , ⃗ y a1 ,⃗ z a1
) at rest. In this case, the translational dynamics is described by the following set of equations:

[

m 1 0 0 m c1 ] [ z1 zc1
]

+ [ c 1 -c 1 -c 1 c 1 ] [ ̇z1 ̇zc1
]

+ [ k 1 -k 1 -k 1 k 1 ] [ z 1 z c1 ] = [ a 1 -a 1 ] u{1} + [ 0 1 ] f zc1 ( 13
)
where f zc1 is an external axial force applied to the enclosing shell. When the actuator is clamped to the beam, the previous set of equations need to be slightly adapted. Firstly, as performed for the mirror load in [START_REF] Davisj | RodenHubble Space Telescope Reaction Wheel Assembly Vibration Isolation System[END_REF], the static inertia of the shell casing was added to the plate node  1 . Secondly. the axial acceleration of the shell becomes equal to that of the supporting node, i.e. zc1 = q 1 {1}. In this case, the equations for the clamped PMA are

[ q 1 f  1 ] = [  (r  1  1 ) 0 0  (r  1  1 ) ⊺ ] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ T  1  1 =𝚫  1 ⋆ Ť 1  1 [ q 1 f  1 ] and ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ f  1 = [ m c1 0 0 J c1 ] q 1 + ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ f pma1 f pma1 = c 1 ̇ez1 + k 1 e z1 + a 1 u{1} e z1 = z 1 -q  1 {1} z1 = -m -1 1 f pma1 ( 14
)
where r  1  1 are the coordinates of the attachment point  1 relative to outer casing CoM  1 , the function  (•) constructs the truncated kinematic transport matrix as defined in [START_REF] Davisj | RodenHubble Space Telescope Reaction Wheel Assembly Vibration Isolation System[END_REF] and 𝚫  1 captures the uncertainties in the location of the center of mass. In terms of input-output behavior, the same PMA dynamics can also be expressed as the following transfer function:

f  1 = D  1 [ u{1} q 1 ]
with

D  1 = 𝚫  1 ⋆ D 1 ∈ ℝℍ 3×4 ∞ ( 15
)
where 𝚫  1 combines the block 𝚫  1 and the uncertainty in the other actuator parameters. The dynamical equations of the second PMA are almost completely similar to eqs. ( 13)-( 15), except for changes in each of the indices (for example  1 becomes  2 and u{1} becomes u{2}).

The subcomponent models defined in eqs. ( 9), ( 11) and ( 14) were subsequently combined into the following global model mapping actuator inputs u to laser displacement y: (16) where

q = [ q ⊺  z 1 z 2 ] ⊺
; q(t) ∈ ℝ n q ; n q = 1514 combines the node displacements q  of the flexible plane with those of the proof-masses z 1 , z 2 . The matrices , , , ,  are the uncertain mass, damping, stiffness, input and output matrices. An equivalent representation of the same dynamical model can be deduced by combining each of the subsystem LFTs from eqs. ( 10) to ( 12) and ( 15) into the following global LFT representation:

G sys = 𝚫 sys ⋆ I 2n q s ⋆ ∼ G sys with 𝚫 sys = diag ( 𝚫  1 , 𝚫  2 , 𝚫  1 , 𝚫  2 , 𝚫  L , 𝚫  , 𝚫  , 𝚫  ) (17)
Fig. 6 illustrates the internal structure of this overall LFT model and the interconnections between the various subcomponents. This block diagram LFT representation of G sys offers some advantages over the system level description given in [START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF]. Firstly, this method of assembly can be more versatile since the different component blocks can be easily interchanged with others from an existing library (see Ref. [START_REF] Alazardc | CumerSate Llite Dynamic Toolbox: Principle, User-Guide and Tutorials[END_REF] for an example). Secondly, since the block uncertainty is already isolated at the component level, a low order global uncertainty block 𝚫 sys can be easily constructed by just concatenating the individual uncertainty blocks. In the system level description from ( 16), care must be taken to ensure that the uncertain matrices are properly factorized to avoid unwanted repetitions of the uncertain parameters. However, the first representation also comes with some distinct advantages. In particular, the first representation was used to compute the nominal mass-normalized modal matrix

𝚽 = [ 𝝋 1 𝝋 2 … 𝝋 n q
] ∈ ℝ n q ×n q of the overall mechanical system. This matrix was calculated by assuming no damping and nominal values for the mass and stiffness matrices and satisfies

𝚽 ⊺ M𝚽 = I n q ; ( K -𝜔 2 i M ) 𝝋 i = 0 for i = 1, … , n q ; 𝚽 ⊺ K𝚽 = diag ( 𝜔 2 1 , … , 𝜔 2 n q ) ( 18
)
where 𝜔 i is the modal frequency associated with every mode shape 𝝋 i . The matrix Φ serves two primary purposes:

(a) Mode shape visualization. Each mode shapes 𝝋 i can be visualized in order to get insight into the various interactions between the subcomponents and the different ways the structure can vibrate at the resonant frequencies.

(b) System order reduction by modal truncation. It can be seen from ( 16) that the order of G sys is equal to 2n q = 3028.

Such a high order system would introduce significant numerical difficulties in any subsequent analysis. For this reason, the modal matrix 𝚽 was used to perform an order reduction. In this case, the generalized displacements q can be expressed in terms of the modal coordinates 𝜼(t) ∈ ℝ n q as q = Φ𝜂. From the modal orthogonality condition in (18) it follows that 𝚽 -1 = 𝚽 ⊺ M and therefore 𝜼 = 𝚽 -1 q = 𝚽 ⊺ Mq. Consider now a partitioning 𝚽 =

[ 𝚽 r 𝚽 t ] and 𝜼 = [ 𝜼 ⊺ r 𝜼 ⊺ t ] ⊺
where 𝚽 r ∈ ℝ n q ×n r contains the n r = 12 modal vectors to be retained and Φ t those that will be truncated. The original displacement vector q = Φ r 𝜂 r + Φ t 𝜂 t from ( 16) is therefore approximated with q ≈ Φ r 𝜂 r . Using again the modal orthogo- )

. In this case, [START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF] and Fig. 6a can be updated to reflect the newly reduced order system by replacing the integrator relationship on the state (

[ q ⊺ ̇q⊺ ] ⊺ such that
) 19
The system order is therefore reduced from 2n q = 3028 to 2n r = 24. The reduction procedure is equivalent to rewriting G sys from ( 16) as [START_REF] Bowlerimaging | Extrasolar Giant Planets[END_REF] However, the integrator reduction method proposed in ( 19) is slightly more advantageous, since it can be directly applied to the block diagram LFT representation from [START_REF]PredaRobust Microvibration Control and Worst-Case Analysis for High Pointing Stability Space MissionsPh[END_REF], also shown in Fig. 6. The truncation errors induced by this model reduction can be included into an additional additive uncertainty model as performed in Ref. [START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF].

Grey-box identification

The reduced order analytical model G sys = 𝚫 sys ⋆ Ǧsys detailed in [START_REF] Rauschs | WittrockUnimor ph deformable mirror for space telescopes: design and manufacturingOpt[END_REF]

𝜖𝚫 ⋆ Ǧexp

) -(

𝚫 sys ⋆ Ǧsys )] W r ∞ ≤ 𝛾 ( 21
)
where the weighting function W e was fixed to a unit gain second order band-pass Butterworth filter with a pass-band of 5 Hz-80 Hz. The purpose of this weight is to put more emphasis on the model errors within the control bandwidth. The other weight W r can be used, in the general case, to optimize the model error for a specific class of input signals. For simplicity, this weight was fixed to W r = I 2 in this study and the value 𝜖 = 0.05 was selected to account for only 5% of the total experimental uncertainty 𝚫.

The optimization problem in [START_REF] Huangb | WuExperimental study and numerical simulation of active vibration control of a highly flexible beam using piezoelectric intelligent materialAero[END_REF] was treated as a robust synthesis problem involving the structured block 𝚫 sys . This class of problems is known to be NP-hard even when the controller to be designed is unstructured. Nevertheless, powerful heuristic methods, have been developed to help deal with such problems. In this work, the heuristic method of choice was the nonsmooth structured  ∞ design technique [START_REF]ApkarianNonsmooth mu-synthesisInt[END_REF]. Mixed-uncertainty structured synthesis methods such as the ones presented in

Ref. [START_REF] Da | ApkarianD.NollStructured robust control against mixed uncertaintyIEEE Trans[END_REF] can directly tackle the minimization problem expressed in [START_REF] Huangb | WuExperimental study and numerical simulation of active vibration control of a highly flexible beam using piezoelectric intelligent materialAero[END_REF]. However, in this work a slightly simplified version was optimized. In this case, the worst case error is not calculated across all 𝚫; 𝜎 (𝚫) ≤ 1 but rather across a smaller subset of random samples. This multi-model approach resulted in a faster tuning process but also required a subsequent 𝜇 analysis to verify if the  ∞ norm condition was satisfied across the larger set. If any uncertainty combination was found to invalidate the norm requirement using the 𝜇 analysis, then that particular uncertainty was added to the collection of random samples. The optimization process was then repeated until no other uncertainty could be found or a finite number of iterations was reached. The nominal identified nominal values of the various parameters are provided in Table 1.

Fig. 8 shows the result of this tuning procedure by comparing the gains u{•} → y of the tuned nominal analytical model 𝚫 sys ⋆ Ǧsys to the gains of the uncertain experimental model (𝜖𝚫) ⋆ Ǧexp . It can be seen that the mathematical model reliably fits the experimental model in the control bandwidth of interest of 5-80 Hz. The expected model discrepancy occurring at high frequency most likely occurs due to deviations from the ideal Kirchhoff plate model, nonlinear damping effects or actuator dynamics. However, the analytical model is sufficient to provide a great deal of physical insight into the system dynamics.

For example, the first six mode shapes 𝚽 i ; i = 1, …, 6 were calculated up to a frequency of 110 Hz using [START_REF] Bastaitsd | PreumontSegment ed bimorph mirrors for adaptive optics: segment design and experimentAppl[END_REF]. The mode shapes labeled 1 ○ to 6 ○ are illustrated in Fig. 9 and their respective natural frequency is indicated in Fig. 8. The newly fitted analytical model also enabled the possibility to perform detailed parametric sensitivity analysis such as the one shown in Fig. 8.

Here, several physical parameters were varied within a ±25% range around their nominal values and an upper bound on the transfer u{2} → y was computed using standard 𝜇 analysis. In this way, the impact of each parameter variation can be accurately predicted across the frequency range. For example, modifying the mass m  of the mirror load corresponds to a frequency shift of the first bending mode. Similarly, changing the damping coefficient c 2 of the second PMA, modifies both the damping of the corresponding actuator mode and the damping of the second bending mode. On the other hand, modifying the offset

r z = r  1  1 {3} = r  2  2 {3}
between the center of mass of the two PMAs and the plate attachment nodes causes a shift in both frequency and damping for all the modes above 40 Hz due to the change in the torsional moment of inertia.

Controller design and limits of performance

Control architecture and synthesis methodology

The purpose of the controller K is to produce an adequate control signal u k = Ky k that meets the following requirements: R1. Robust stability: the closed loop shall be stable across all the uncertainties 𝚫 modeled using [START_REF] Komatsuh | UchidaMicrovibration in spacecraftMech[END_REF].

R2. Robust performance: the control signals and the model error between the experimental plant and an adaptable reference model H must be minimized and guaranteed to be below a certain specified level for a restricted subset of uncertainties

𝜖𝚫 with 𝜖 ∈ [0, 1].
To ensure such requirements, the controller was optimized following  ∞ ∕𝜇 design practices, by first assembling the weighted interconnection shown in Fig. 10. This interconnection is composed of the following blocks:

1. Uncertain plant model G sys : For controller synthesis, either the experimental model [START_REF] Komatsuh | UchidaMicrovibration in spacecraftMech[END_REF] or the reduced physical model G sys = 𝚫 sys ⋆ Ǧsys from ( 19) can be used. The choice depends on the particular design objectives.

G exp = 𝚫 ⋆ [ G zw G zu G yw G yu ] from
For example, the analytical model provides the user with virtual access to any physical signal in the overall plant. This can be of considerable advantage since one can rely on these virtual measurements to monitor and possibly manipulate signals for which no sensor is available. On the other hand, the experimental model can be readily derived for a general flexible structure with a more complex shape. Since the aim of the paper is to outline a general control methodology for dynamic reshaping, the choice was made to use the experimental model G exp for control design and rely on the analytical model G sys for a more in-depth worst-case analysis. 3). The difference is that for H(𝛼), the pairs of complex eigenvalues 𝜆 i are parameterized in terms of a 𝛼 as 𝜆 i = f i (𝛼) ±jg i (𝛼). In the general case, the functions f i (𝛼), g i (𝛼) ∶ ℝ → ℝ can be any rational functions. However, for simplicity these functions were chosen as affine in 𝛼. In this case, H can be expressed as the following LFT: [START_REF] Sparksj | JuangSurvey of experiments and experimental facilities for control of flexible structuresJ[END_REF] where n H = 6 is the minimal number of repetitions of 𝛼 in the LFT description, m = 4 is the number of poles up to 70 Hz in G yu and the matrices B i and C i are the ones from [START_REF] Royers | PellegrinoUltralight ladder-type coilable space structures[END_REF] In this case w p was chosen as a 4th order bandstop filter with maximum tracking error ASD of 50 𝜇rad∕ √ Hz in the 7 Hz-70 Hz stopband and a maximum of 500 𝜇rad∕ √ Hz outside of it. To ensure that the requirements are not conflicting or too conservative, the values for the performance weights W u and W p were chosen after understanding some of the fundamental limits of performance imposed by the control architecture. Details about this analysis are provided in section 3.2.

Structured adaptable controller K(𝛼):

To facilitate the implementation of the control law but also provide sufficient adaptability to changes in the reference model H(𝛼), the following affine structure K(𝛼) is imposed on the controller:

where A K0 , A K1 , B K0 , B K0 , C K0 , C K0 , D K0 , D K0 are real matrices of appropriate dimension. The controller order n qk = 6, the number of repetitions n K = 6 of the scheduling parameter 𝛼 and the initial values for the controller matrices were chosen after following the iterative procedure detailed in Ref. [START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF]. To speed up the subsequent optimization, the matrices A K0 , A K1 are constrained to be tridiagonal. For ensuing experiments, the controller is implemented in discrete-time following a Tustin transformation with a sampling time T = 0.5 ms synchronized with the 2 kHz sampling frequency of the autocollimator sensor. In this case, each of the integrators 1 s in ( 23) is replaced with

T 2 1+z -1 1-z -1
, where z -1 is the unit delay. In compact notation this is equivalent to: [START_REF] Hydeactive | Vibration Isolation for Precision Space StructuresPh[END_REF] Using the previous component definitions, the closed-loop relationship between the various input and output signals shown in Fig. 10 and the controller K can be determined. In the open loop case, the following mapping exists between the input and output signals shown in Fig. 10: (25) Closing the loop with the controller K(𝛼) such that u k = K(𝛼)y k results in the following new mapping between the signals: [START_REF] Allegranzal | TosoActuators for space applications: state of the art and new technologiesProc[END_REF] where

S(𝛼) = [ I + G yu K(𝛼)
] -1 denotes the Output Sensitivity Function and the dependency on 𝛼 of K, H and S was omitted for clarity. After closing the uncertainty channels with the scaled uncertainty 𝜖𝚫 ; 𝜖 ∈ [0

⋆ P ⋆ Ǩ ⋆ ( 𝛼I n K ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ M(𝛼) ∞ ≤ 𝛾 (27)
Any controller for which 𝛾 < 1 satisfies the robust performance and stability requirements. More precisely, it represents a robust and structured Linear Parameter Varying (LPV) synthesis problem. However, in this work the time-varying aspect of the parameter 𝛼 is not explicitly taken into account and the constraints on the L 2 system gain in [START_REF][END_REF] are only enforced for fixed values of 𝛼. Therefore, the optimization is comparable in complexity to the grey-box fitting method proposed in [START_REF][END_REF] and the same set of nonsmooth  ∞ tools were used to address the problem. In order to improve readability of the paper, additional details about all the various steps and heuristics needed to solve [START_REF][END_REF] were suppressed. A complete explanation of the process used to obtain the structured LPV controller Ǩ ∈ 𝕂 using non-smooth  ∞ tools can be found in Refs. [START_REF] Predaj | FalcozRobust microvibration mitigation and pointing performance analysis for high stability spacecraftInt[END_REF][START_REF]PredaRobust Microvibration Control and Worst-Case Analysis for High Pointing Stability Space MissionsPh[END_REF][START_REF] Predaj | FalcozA H-infinity/mu solution for microvibration mitigation in satellites: a case studyJ[END_REF]. In the case of the optimization given in ( 27), a controller Ǩ was found to achieve a performance level of 𝛾 = 0.989. Fig. 12 shows the gains and pole map of this controller

K(𝛼) = Ǩ ⋆ ( 𝛼I n K ) for different values of the scheduling parameter 𝛼 ∈ [-1, 1].
It is important to understand the shape of the resulting controller and also the trade-offs involved in the overall control design. A detailed analysis of the limits of closed-loop performance is therefore provided in the following subsection.

Limits of performance and trade-off analysis

During the requirements specification phase, it is important to have an intuition about what can be achieved using feedback, subject to the various constraints imposed by available bandwidth, actuator limitations, environmental noise and model uncertainty. This sort of understanding about the limits of performance can simplify the process of selecting the desired performance requirements and the corresponding weights prior to the controller tuning. Additionally, before any experimentation, it is important to understand the mechanisms of action of the controller and cross-check the results returned by automatic synthesis tools. The goal of this section is to show how this trade-off analysis was performed prior to any control synthesis and how the results can be used to explain the gains of the synthesized controller shown in Fig. 12.

For any stabilizing controller K(𝛼), the sensitivity S(𝛼) = [ I + G yu K(𝛼) ] -1 is subject to the following types of constraints for all fixed values of 𝛼 ∈ [-1, 1]: 1. Analytical constraints. Bode's Sensitivity Integral also known as the waterbed effect can be seen as a conservation law valid for all stabilizing controllers [START_REF]SteinRespect the unstableIEEE Control Syst[END_REF]. Briefly, this constraint states that if the controller pushes on the sensitivity function S at one frequency it automatically increases its value at some other frequencies by the same amount. More precisely, consider an arbitrary MIMO plant G yu and a controller K compatible in size such that each of the entries in the open loop transfer matrix G yu K are rational functions with at least two more poles than zeros. In this case, if the closed loop is stable, then the singular values 𝜎 i of S with i = 1, …, n 𝜎 satisfy [START_REF] Freudenbergd | LoozeFrequency Domain Properties of Scalar and Multivariable Feedback Systemsvol[END_REF]:

n 𝜎 ∑ i=1 ∫ ∞ 0 log 𝜎 i [ S (j𝜔)
]

d𝜔 = 𝜋 N p ∑ i=1 Re(p i ) (28)
where p i represents the N p right-hand plane poles. The assumption on the excess number of poles is not restrictive for the type of plants G yu considered in this study, namely flexible structures. This is because both actuators and sensors have a strong roll-off above a certain frequency due to the amplifiers. In the case of the setup considered in this paper, both the ⊺ of the closed loop transfer matrix M from [START_REF] Allegranzal | TosoActuators for space applications: state of the art and new technologiesProc[END_REF]. A brief examination of this matrix reveals that the controller K has an influence on all of the sub-blocks. However, the algebraic dependency on K is different in every case and sometimes of an inverse nature. Therefore it is not possible to make the gain in any sub-block arbitrarily small without amplifying the gain of another block. Furthermore, each performance requirement places a different constraint on the sensitivity S and is important that these requirements are not conflicting. Consider that the optimization [START_REF][END_REF] found a controller that achieves a certain peak gain of 𝛾 at system level. This implies that the gain between any pairs of disturbance and performance signals is also kept below 𝛾. )W r (j𝜔)) ≤ 𝛾

; ∀𝜔 ∈ ℝ (30)
Using the submultiplicative property of the induced matrix 2-norm, the previous inequality holds if

𝜎 ( H(j𝜔) -S(j𝜔)G yu (j𝜔) ) ≤ ( 𝛾 𝜎 ( W p (j𝜔) ) 𝜎 ( W r (j𝜔) ) = γ(j𝜔) ) ; ∀𝜔 ∈ ℝ (31)
For clarity, the dependency on the frequency 𝜔 will be subsequently omitted. Considering now that,

| | 𝜎 (B) -𝜎 (A) | | ≤ 𝜎 (A -B) ; ∀A, B ∈ ℂ n×m , the previous inequality becomes 𝜎 (H) -γ ≤ 𝜎 ( SG yu ) ≤ 𝜎 (H) + γ (32)
Finally, as S is a scalar function, it follows from the homogeneity of the matrix 2-norm that 𝜎

( SG yu ) = |S| 𝜎 ( G yu
) and therefore (32) reduces to

𝜎 (H) - γ 𝜎 ( G yu ) ≤ |S| ≤ 𝜎 (H) + γ 𝜎 ( G yu ) (33)
It is clear that in frequency regions where 𝜎 (H) ≫ γ the bounds become very tight. In the case of the control effort performance, a peak gain 𝛾, i.e. 𝜎

𝜎 ( G yu ( I -KSG yu )) ≤ 𝜎 ( G yu ) 𝜎 ( I -KSG yu ) ⟺ |S| ≤ 𝜎 ( I -KSG yu ) ≤ γ (35)
For brevity, the conditions for the robust performance and stability requirements (i.e. the ones involving the w Δ and z Δ channels) as well as the nominal performance conditions involving the sensor noise channel d n are omitted. 33) tightly constraints the sensitivity function around the regions where significant control action is applied, i.e. regions where the gain of S is small. Additionally, the waterbed effect can be clearly observed in the growth of the sensitivity around regions where the same function is forced to a low value. Since the control effort requirements limits the maximum allowed increase in the sensitivity, it therefore also limits to some extent the maximum bandwidth over which significant attenuation of S can take place. In general, care must be taken to make sure that all distortions away from unity of the sensitivity function happen within the allowed bandwidth to avoid unwanted interactions with unmodeled dynamics [START_REF]SteinRespect the unstableIEEE Control Syst[END_REF]. A more in-depth worst-case analysis is provided in section 4.

Performance impact of actuator and sensor placement

The actuator and sensor placement plays a critical role on the tracking performance of the closed loop system. Intuitively, the actuators need to have sufficiently authority over the plant dynamics and mode shapes in the bandwidth of interest in order to reshape the dynamics and track the reference model. Within the same bandwidth, the modes of interest need to be observable from the sensors measurements to perform feedback control. The dynamical equations of the plant from ( 16) as well as the modal matrix Φ from (18) can be recomputed for different actuator and sensor layouts. Subsequently, the new mode shapes of the flexible structure can be visualized, in a manner similar to the one in Fig. 9, to understand the new structural dynamics in a given bandwidth. Furthermore, the open loop gains from the actuators towards the sensor measurements, shown in Fig. 8, can be recomputed to assess the controllability and observability of the new layout. These gains can then be optimized in order to increase control authority and tracking performance by tuning the location of the various actuators and sensors. This actuator and sensor layout optimization was beyond the scope of this paper but would be an interesting topic for future research efforts.

Performance and stability analysis

Worst case analysis

Before implementing the control law, the robust stability of the closed loop interconnection 𝚫 ⋆ M shown in Fig. 10 ; 𝜎 (𝚫) ≤ 1 given in [START_REF] Bastaitsd | PreumontSegment ed bimorph mirrors for adaptive optics: segment design and experimentAppl[END_REF]. Fig. 14b illustrates side views of the bounds (upper and lower) of the same function but also bounds on 𝜇 𝜹 computed only with respect different subsets of uncertainty. The minimums in stability margin, i.e. the peaks of 𝜇 𝜹 occurs mostly around the regions where control action was applied to dampen the modes. However, even in the combined uncertainty case, 𝜇 𝜹 remains below 0.56 and therefore the loop can tolerate an increase in the uncertainty of 1/0.56 = 78% while maintaining stability. The impact of uncertainties on different performance indicators was also assessed, prior to experimentation, using structured singular value calculations. Fig. 15 illustrates the upper bounds on the peak gain for different performance signals, across all frequencies 𝜔 and scheduling parameter values 𝛼. The first performance transfer that was studied was

[ d r {1} d n ] ⊺ → p
corresponding to the maximum model tracking error using sensor noise and reference signals in the first actuator. In this case, it has been observed that the worst-case gains channel are mainly sensitive to modal uncertainties 𝚫 mod around the flexible mode frequencies. This drop in performance can be explained by the fact the uncertainties act only on the plant model G exp and not on the adaptable reference model H(𝛼) in the closed-loop interconnection (see Fig. 10). Therefore, in the worst-case, because of the significant modal uncertainty considered for the plant model, big differences can occur around the resonance frequencies between the gains of G exp and those of the reference model (see Fig. 12a for the gains of H(𝛼)). Simply put, the peaks in uncertain plant G exp can occur at frequencies that are very far from those of the desired reference model. Therefore, the controller needs to work much harder in this situation in order to reshape the plant dynamics to match those of the reference model. However, due to the hard constraints on the actuators this is not physically possible in all scenarios and a drop in the tracking performance . The gains were computed for the controller scheduled at 𝛼 ∈{-1, 0, 1} using the analytical model G sys and the interconnection from Fig. 10.

naturally occurs. More details about these performance limitations are discussed in section 3.2.

The worst-case gains of a second performance channel It must be mentioned that the methodology outlined in this paper, relies on the assumption that the system undergoes small deflections and that the dynamics remain fairly linear. This assumption covers a wide variety of space applications since most spacecraft tend to avoid by design large deflections and nonlinearities in the structural dynamics. However, if this condition doesn't hold, the specific type of nonlinearity needs to be considered during the tuning and validation procedures. Most likely, the design procedure would be based on a couple of iteration between a linear controller design and a Monte-Carlo analysis on the full nonlinear plant model.

Experimental results

Satisfied with robust stability and robust performance characteristics, the controller was subsequently evaluated in a series of experiments. using the experimental (section 2.2.1) and analytical plant models (sections 2.2.2 and 2.2.3). While both models reliably anticipate the closed-loop behavior, in the neighborhood of 35 Hz the analytical model correctly predicts the clear separation of the two PMA modes around this frequency. This is a natural consequence of the fact that the experimental model only fitted a single mode at that frequency (see Fig. 3). On the other hand, the analytical model physically models the two modes as distinct but close in frequency (see the corresponding mode shapes in Fig. 9). Predictions about the closed-loop ASD and CRMS of the control signals were also made in Fig. 17b. The figure illustrates these prediction on top of the experimental results for the closed loop scheduled at 𝛼 ∈{-1, 0, 1}. As visible in these plots, the control signals can be accurately predicted before the experiments are performed.

Interestingly, the analytical model can be used to perform a closed loop examination of certain signals without relying on additional sensors. The model therefore acts as a window through which the typically black-box dynamics can be observed. For example, Fig. 18 displays the calculated nominal ASD spectrum of the two proof-mass displacements, i.e. the singular values of This result is consistent with the results presented in Fig. 16.

Conclusions

This paper presented a complete experimental modeling, analytical modeling and control design methodology to actively alter the response of a flexible structure. The framework ensures that the closed loop plant dynamics robustly matches that of a given adaptable reference model even in the presence of significant model uncertainty and sensor noise. The controller synthesis procedure includes a step by step description of the design process including details about the interplay between the different requirements and the limits of performance. The extended robust performance assessment performed prior to the various closed loop experiments was a necessary step to ensure the safety and reliability of the proposed control law. The work also demonstrated that accurate predictions can be made about the closed loop behavior, provided that sufficient effort is placed in extracting a physical model of the plant and characterizing the different disturbances acting on it. The paper outlines the process of deriving such an uncertain model either in a purely data-driven fashion or as an analytical model calibrated from experimental data. These results provide the confidence to extend the proposed methodology to more complex scenarios involving other flexible space structures equipped with more sensors and actuators.

For the partitioned matrices M = [M 11 M 12 M 21 M 22 ] 12 N 21 N 22 ]M ⋆ N 11 M 21 N ⋆ M 22 ]

 =11122122122122112122 of appropriate dimensions, the Redheffer star product of M and N is M ⋆ N = [ 12 (I -N 11 M 22 ) -1 N 12 N 21 (I -M 22 N 11 ) -1 N where the existence of inverses is assumed. If M or N don't have an explicit 2 × 2 structure, then the star product reduces to a linear fractional transformation (LFT). In this case the operator ⋆ is associative, the lower LFT of M and N is M ⋆ N = M 11 + M 12 N(I -M 22 N) -1 N 21 and the upper LFT of M and N is N ⋆ M = M 22 + N 21 N(I -M 11 N) -1 M 12 . The set ℝℍ n×m ∞ represents the set of finite gain transfer matrices with n outputs and m inputs. In the case of a SISO transfer function, this set reduces to ℝℍ ∞ . For G ∈ ℝℍ n×m ∞ , the value ‖G‖ ∞ represents the L 2 system

Fig. 1 .

 1 Fig. 1. (a) Top-down photo of the system. (b) Schematic diagram of the setup.

Fig. 2 .

 2 Fig. 2. (a) Actuator voltages and measured angular deflection during the first identification experiment. Labeled time regions: a ○ environmental noise phase; b ○ first actuator phase; © second actuator phase. (b) Maximum sensor amplitude spectral density and cumulative root mean square (dotted line) under the influence of environmental perturbations and zero actuator signals. Labeled points: 1 ○ first bending mode; 2 ○ power supply electric noise.

Fig. 3 .

 3 Fig. 3. Gains of the estimated FRFs and nominal fitted plant model G yu together with the minimum coherence spectrum across all three identification experiments: (a) u{1} → y channel; (b) u{2} → y channel.

3 .

 3 The estimates G y k u{•} (𝜔) of the frequency response functions (FRF) from inputs u{•} to output y k together with the corresponding coherence functions Γ y k u{•} (𝜔) were computed as:

( 8 )

 8 Fig. 4a highlights this LFT structure of the global uncertain model while fig. b shows the separate effects of each type of uncertainty on the gains of the transfer function from the input u{2} to the output y.

Fig. 4 .

 4 Fig. 4. (a) Internal structure of the uncertain experimental model of the plant. (b) The effects of different sets of uncertainties on the gains of the channel u{2} → y.

Fig. 5 .

 5 Fig. 5. (a) Finite element model of the flexible plate together with corresponding coordinate frames and illustration of the four node plate element. (b) Diagram of the first proof-mass actuator model in unclamped configuration illustrating the displacements z 1 and z c1 relative to the rest state (dashed lines).

(

  ; ⃗ x e , ⃗ y e ,⃗ z e ) corresponding to the central node at the end of the flexible part of the plate. Coincident with this attachment frame is a body frame (  ; ⃗ x l , ⃗ y l ,⃗ z l ) fixed at the center-of-mass of the load. The load was coupled to the plate dynamics by first translating plate accelerations q = [ ; ⃗ x e , ⃗ y e ,⃗ z e ) to accelerations q = [ z θ ψ ] ⊺ of the load body frame

(O 3 I 3 ]

 33 ; ⃗ x e , ⃗ y e ,⃗ z e ) and rotations around the ⃗ z e are not considered in the flexible plate model. The blocks 𝚫  L and 𝚫  capture the uncertainties on the geometric and inertial parameters of the mirror load. The selection matrix S  is needed to truncate the kinematic transport matrix [ since each plate node only has three coordinates (z  • , 𝜽  • and 𝜑  • ) instead of the usual six for a general rigid body.

Fig. 6 .

 6 Fig. 6. (a) Block diagram of the uncertain plant together with the various subcomponents written in LFT form. Note: the grid on the D flexible plate block corresponds to the FEM node grid and the colored circles indicate the various attachment nodes. (b) Equivalent global LFT form.

r

  M𝚽 r = I n r it follows that 𝚽 + r = 𝚽 ⊺ r M is a generalized left inverse of Φ r i.e. 𝚽 +r 𝚽 r = I n r and therefore

Fig. 7 .

 7 Fig. 7. (a) Interconnection used for grey-box identification. (b) Solid lines: upper bounds computed using 𝜇-analysis on the gains of the u{2} → y channel for ±25% variation around the nominal values of various physical parameters. Dashed lines: gains for different values of the corresponding uncertain parameter.

2 . 3 .

 23 Disturbance weights W n and W r : The measurement noise weight W n was introduced in (1) and is used to model the upper bound on the expected ASD spectrum of the experimental closed-loop noise measurements. Similarly, the filter W r specifies the upper bound on the ASD of the closed-loop reference signals r. The reference inputs r used in the closed loop experiments were chosen to be the same as the signals u used during the experimental identification phase (see section 2.2.1 and Fig. 2a). Therefore, W r was chosen as a scaled 4th order band-pass Butterworth filter with a pass-band of 5 Hz-200 Hz and the same amplitude as the one used to color the white noise open loop identification signals. Adaptable reference model H(𝛼): This model represents the target dynamic response for the close-loop plant and is parameterized in terms of a normalized scheduling parameter 𝛼 ∈ ℝ ; |𝛼| ≤ 1. In this paper, H(𝛼) ∈ ℝℍ 1×2 ∞ is based on the nominal open loop plant G yu expressed in canonical modal from in (

4 .

 4 Fig. 11 shows the gains and pole maps of H for different values of 𝛼 compared to the open loop plant G yu . It can be seen that for 𝛼 = -1, the damping coefficients of the PMA modes around 35 Hz and the second bending mode around 63 Hz are increased by over an order of magnitude compared to the open loop. When 𝛼 = 1, the damping of the first open loop bending mode around 10 Hz is increased by two orders of magnitude. At the same time, the natural frequency of this flexible mode is raised up to 14 Hz. The torsional mode around 45 Hz is kept at the open loop values.The reason why H(𝛼) was chosen to have this dependency on 𝛼 was to demonstrate that the controller can distinctly and selectively change the structural behavior of the plant at different frequencies. This particular value for H was also selected based on the limits of performance analysis detailed in section 3.2. Performance weights W u and W p : The purpose of the weight W u = w -1 u I 2 ; w u ∈ ℝℍ ∞ is to impose a desired closed loop upper bound of |w u (j𝜔)| on the worst-case ASD of the actuator signals u at different frequencies 𝜔. In this way, the specification on the maximum actuator RMS can be guaranteed across any frequency band since the RMS is equal to the square root of the area under the PSD curve. For the specific actuator used in the study, the RMS of the input voltage must stay below 0.1 V to avoid damage to the coil. Therefore w u was fixed to w u = 5.0 • 10 -3 to ensure a bound of 5 × 10 -3 V∕ √ Hz on the ASD and a maximum RMS of 0.1 V up to 200 Hz. Likewise, the weight W p = w -1 p ; w p ∈ ℝℍ ∞ is used to specify the desired upper bound |w p (j𝜔)| on the ASD of the error p between the output y h of the reference model H and the output y of the plant model G exp .

Fig. 8 .

 8 Fig. 8. Comparison between the gains of the experimental and analytical models: (a) u{1} → y channel; (b) u{2} → y channel. (Note: the region shaded in corresponds

Fig. 9 .

 9 Fig. 9. Nominal mode shapes of the undamped analytical model up to 110 Hz numbered in Fig. 8: 1 ○ 1st bending mode; 2 ○ 1st actuator mode (in-phase proof-mass displacements); 3 ○ 2nd actuator mode (out-of-phase proof-mass displacements); 4 ○ 1st torsional mode; 5 ○ 2nd bending mode; 6 ○ 2nd torsional mode.

Fig. 10 .

 10 Fig. 10. (a) System architecture used for controller synthesis and worst-case analysis. (b) Equivalent standard form of the interconnection.

Fig. 11 . 29 ) 2 .

 11292 Fig. 11. (a) Comparison between the gains of the adaptable reference model H(𝛼) for different 𝛼 ∈ [-1, 1] and the nominal open loop system G yu . (b) Comparison between the poles of H(𝛼) and G yu (Note: due to pole symmetry, only the upper complex half-plane is shown).

  In the case of the pair d r → e p this is related to the model tracking error in response to the reference signals. A peak gain of 𝛾 for this transfer implies that 𝜎(W p (j𝜔)(H(j𝜔) -G yu (j𝜔)(I -K(j𝜔)S(j𝜔)G yu (j𝜔)) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ = S(j𝜔)G yu (j𝜔)

Fig. 12 .

 12 Fig. 12. (a) Gains and (b) pole maps of the controller K(𝛼) for different values of the scheduling parameter 𝛼 ∈ [-1, 1].

Fig. 13 .

 13 Fig. 13. Closed loop sensitivity functions S(𝛼) for different values of the controller scheduling parameter: (a) 𝛼 = -1; (b) 𝛼 = 0; (c) 𝛼 = 1. Shaded regions correspond to the algebraic constraints imposed by the requirements on the model tracking (33) and control effort (35).

Fig. 14 .

 14 Fig. 14. Robust stability plots: (a) upper bound on 𝜇 𝜹 across a dense grid of frequencies and scheduling parameters. (b) side views of the bounds (solid line for upper and dashed lines for lower) on 𝜇 𝜹 computed only with respect to different subsets of uncertainty.

  was evaluated. This assessment was done by calculating the bounds on the function 𝜇 𝛿 ( M w→z (𝛼, j𝜔) ) across a dense grid of frequencies 𝜔 and for multiple values of the scheduling parameter 𝛼 ∈ [-1, 1]. The function 𝜇 𝛿 (M) represents the structured singular value [39] defined as zero if no uncertainty 𝚫 makes I -M𝚫 singular and 𝜇 𝜹 (M) = 1∕min 𝚫 { 𝜎 (𝚫) , det (I -M𝚫) = 0 } otherwise. This function therefore offers precise information about the magnitude of uncertainty need to destabilize the loop at any frequency. Furthermore, efficient algorithms [40] have been developed to bound this function. Fig. 14a displays the upper bound on this function for the set of uncertainties 𝚫 = diag ( 𝚫 mod , 𝚫 add , 𝚫 mul )

Fig. 15 .Fig. 16 .Fig. 17 .

 151617 Fig. 15. Upper bounds on the gains of different performance channels with respect to multiple and uncertainty sets: (a) the model tracking error channel [ d r {1} d n ] ⊺ → p;

Fig. 18 .z 1 z 2 ]

 182 Fig. 18. Predicted nominal closed loop ASD spectrum of the two proof-mass displacements, i.e. the gains of the channel d r → [ z 1 z 2 ]

[ d r {1} d n ] ⊺ → u{1} are shown in Fig. 15 .

 15 The gains correspond to the maximum control effort ASD in the presence of sensor noise and reference signals in the first actuator. It can be observed that the channel maintains values close to nominal ones even in the presence of significant model uncertainty.

  Fig. 16 shows the gains of the experimentally identified closed loop frequency response functions r{•} → y for five different experiments with the controller scheduled at 𝛼 ∈{-1, -0.5, 0, 0.5, 1}. The reference signals r used during the experiments are the same as the open loop identification signals u shown in Fig. 2a and discussed in section 2.2.1. Fig. 16 also includes the coherence plots and the predicted closed loop uncertain response. The predictions were made using the experimental model arranged in the architecture shown in Fig. 10. To compute the different predictions, the model uses a reduced uncertainty subset scaled to 10%, i.e. 𝚫 ; 𝜎 (𝚫) ≤ 0.1. Excellent agreement between experiment and theoretical prediction can be observed. For the same reduced uncertainty set and 𝛼 = -1, Fig. 17a illustrates the different closed loop predictions made

Fig. 19 .

 19 Fig. 19. Experimental time domain performance in open and closed loop for a 10.2 Hz sine excitation. In the first phase, the sine input in one of the actuators excites the first bending mode. In the second phase, the excitation is stopped and the controller (scheduled at 𝛼 ∈{-0.5, 1}) is activated.

2 ]

 2 using the interconnection from Fig.10and replacing the plant model G exp with the analytical model G sys . The gains are calculated in both open loop and in closed loop with 𝛼 ∈{-1, 0, 1}. The cross-coupling between the two actuators as well as the interactions with the flexible modes of the supporting plate are thus revealed. For example, when 𝛼 = -1 and the controller focuses more on dampening the first bending mode around 10 Hz, the mass displacements internal to the actuators are also significantly reduced. Fig.19reveals the time domain experimental performance results in open and closed loop. In the first half of the experiment, the first actuator is driven by a constant sine input at 10.2 Hz around the first bending mode frequency of the plate. In the second phase, the excitation is stopped and the controller is activated with 𝛼 ∈{-0.5, 1}. Compared to the open loop, the controller significantly increases the damping of this mode and the structure settles back to equilibrium almost 100 times faster for 𝛼 = 1.

Table 1

 1 Nominal parameter values of the overall system.

	Subsystem	Parameter	Description

moment of inertia of the outer casing in CoM frame

  contains k real uncertain parameters 𝛿 p • in the block diagonal block 𝚫 sys . In order to calibrate the nominal values of these various physical parameters, a grey-box identification procedure was used. The proposed method relies on the interconnection shown in Fig.7abetween G sys and the uncertain experimental model G exp = 𝜖𝚫 ⋆ Ǧexp from[START_REF] Komatsuh | UchidaMicrovibration in spacecraftMech[END_REF], where the parameter 𝜖 ≤ 1 was introduced to scale the set of normalized experimental uncertainty. The idea is to search for a set of parameter values for the uncertainty block 𝚫 sys that minimizes the worst-case weighted additive error between the two models across the reduced subset of experimental uncertainties 𝜖𝚫.

					This takes the
	form of the following optimization procedure:
	mininimize 𝚫 sys	𝛾 s.t.	sup 𝚫 ; 𝜎(𝚫)≤1	W e	[(

Throughout this paper, the two-sided ASDs and PSDs were calculated.