Motion planning and tracking control for coupled flexible beam structures

Control Engineering Practice 84 (2019) 389–398

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Motion planning and tracking control for coupled flexible beam structures Andreas Kater, Thomas Meurer ∗ Chair of Automatic Control, Kiel University, Kaiserstraße 2, 24143 Kiel, Germany

ARTICLE

INFO

Keywords: Distributed parameter system Flatness-based feedforward control Motion planning Coupled structure Euler–Bernoulli beam Experiment Flexible structure

ABSTRACT A systematic flatness-based motion planning approach for a structure consisting of coupled flexible bending beams is introduced. The beams are equipped with in-domain piezoelectric macro-fiber composite (MFC) actuators and are coupled by an elastic string, modeled as spring, which connects the tip masses of adjacent beams. This setup provides a challenging benchmark experiment to evaluate control concepts. Hamilton’s principle leads to the equations of motion involving spatially varying parameters due to the different material characteristics of the carrier structure and the attached MFC-actuators. The flat output is constructed by exploiting the spectral system representation and enables a differential parameterization of the system states and inputs. This parameterization transfers the motion planning problem into a trajectory assignment problem for the flat output. The resulting flatness-based feedforward control law is analyzed and evaluated at an experimental setup. To increase the performance and the robustness w.r.t. model uncertainties and external disturbances the flatness-based feedforward control is amended by a feedback controller with observer in the spirit of the two-degrees-of-freedom control approach.

1. Introduction Adaptive elasto-mechanical structures, also called smart structures, are typically characterized by a spatially extended carrier structure with embedded actuators and/or sensors (Banks, Smith & Wang, 1996; Clark, Saunders, & Gibbs, 1998; Janocha, 2007; Stanewsky, 2001), e.g., to achieve transient shape control. This ability offers various opportunities for a large field of applications including, e.g., adaptive optics (Preumont, 2011), spoilers in automotive engineering (Mohd Jani, Leary, Subic, & Gibson, 2014) and adaptive wings (Barbarino, Bilgen, Ajaj, Friswell, & Inman, 2011; Yousefi-Koma & Zimcik, 2003) in space and aerospace engineering. Furthermore, with a properly developed feedforward and feedback control strategy smart structures are able to reach a desired dynamic state like resonant behavior or in contrary to reduce unintended oscillations in transient rest-to-rest motion (Banks, Smith et al., 1996; Meurer, 2013). The focus of this contribution is the transient shape control of an elastically coupled set of bending beams (depicted in Fig. 1) with embedded piezoelectric macro fiber composite (MFC) actuators. Bending beams are often used in literature as a benchmark example or test setup for shape control based on infinite-dimensional control approaches (Kugi & Thull, 2005; Kugi, Thull, & Kuhnen, 2006), 𝐻2 as well as 𝐻∞ control (Banks, Demetriou, Michael & Smith, 1996; Kugi, 2001; Morris & Vest, 2016), fuzzy control (De Abreu & Ribeiro, 2002), or neural network control techniques (Abdeljaber, Avcia, & Inman, 2016; Qiu, Zhang, & Ye, 2012). In the previous work (Kater & Meurer, 2015) of the authors a twodegrees-of-freedom control approach was presented for a simplified

setup. The present contribution introduces multiple extensions by dealing with different damping effects, in-domain measurements requiring the design an implementation of a state observer to evaluate the state feedback control, and a more complex configuration resulting on the one hand in the need for a higher order system description to capture the dynamics. On the other hand, the systematic model order reduction is included to decrease the computational burden for the controller and observer implementation. The considered flexible beam structure of Fig. 1 is motivated by the desire to create a rather challenging benchmark example to evaluate control strategies for flexible structures. While the control design for single beam structures has been studied extensively this coupled setup is significantly complicated by the elastic connection between the beams. This results in the clustering of eigenvalues with small real part, i.e., weakly damped modes. Moreover the coupled beam structure can be considered as an 1D approximation of more complex 2D or 3D structures that comprise beam elements which are connected by thin films or membranes. The mathematical modeling, presented in Section 2, leads to coupled partial differential equations (PDEs) for the beam’s deflections. There are in principle two approaches to dynamically analyze and to develop control algorithms for distributed parameter systems: the ‘‘early lumping’’ and the ‘‘late lumping’’ approach. ‘‘Early lumping’’ approaches directly apply approximation techniques like the finite difference, the finite volumes or the finite elements method to the

∗ Corresponding author. E-mail addresses: [email protected] (A. Kater), [email protected] (T. Meurer).

https://doi.org/10.1016/j.conengprac.2018.12.012 Received 5 September 2018; Received in revised form 16 December 2018; Accepted 18 December 2018 Available online xxxx 0967-0661/© 2018 Elsevier Ltd. All rights reserved.

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

2. Closed-loop tracking control is addressed using a two-degrees-offreedom control concept by amending the feedforward control with a linear quadratic regulator and a Kalman filter to estimate the tracking error. 3. Systematic model order reduction is included to decrease the computational burden in view of the implementation. 4. Weighting matrices to formulate the quadratic minimization problems for the control and the observer are determined by exploiting the Hankel singular values used for model order reduction. 5. Experimental results illustrate the resulting tracking performance even for very weakly damped coupled structures. 2. Mathematical modeling In the following the equations of motion for the transversal deflection of the considered coupled beam structure (Fig. 1) are derived under the assumptions of the Euler–Bernoulli beam theory. For simplification of notation the direction tag of the independent spatial coordinate is discarded, i.e., 𝑧 = 𝑧1 . The deflection 𝑤𝑗 (𝑧, 𝑡) of each beam is aligned with the 𝑧3 -axis, where 𝑗 denotes the related beam. For modeling and control design the structure is assumed to consist of 𝑟 ∈ N bending beams with identical dimensions. All beams are clamped at the position 𝑧 = 0 and free at 𝑧 = l𝑗 = l . Any beam is equipped with a tip mass 𝑚𝑗 at its free end. Neighboring beams are connected with each other by a string (see Fig. 1), which is modeled as a spring with the spring constant 𝑘𝑖,𝑗 with {𝑖, 𝑗} ∈ {{1, 2}, … , {𝑟 − 1, 𝑟}}. The beams are equipped with 𝑝𝑗 piezoelectric MFC-patch pairs, which are symmetrically located at a longitudinal position 𝑧𝑗,a with 𝑘 ∈ {1, … , 𝑝𝑗 } on the front and the 𝑘 back of beam 𝑗. Note that every parameter marked with the index a𝑘 refers to the 𝑘th MFC-patch pair. Applying an asymmetric voltage to each patch pair in the form

Fig. 1. Schematics of the coupled beam structure. The red dots mark the points considered for the shape control (subsequently called points of interest).

distributed parameter system to obtain a finite-dimensional representation in the form of ordinary differential equations (ODEs). Contrary the ‘‘late lumping’’ approach takes into account the spatial and time dependent description of the system. Due to this, the system dynamics is preserved throughout analysis and design but this comes at a higher effort. The motion planning objective is to obtain a smooth rest-to-rest motion of the transversal deflection of the structure. Therefore, a twodegrees-of-freedom control approach consisting of a feedforward control solving the motion planning problem, and a feedback control to achieve robustness against model uncertainties and external disturbances, is utilized. For the feedforward control design and the motion planning the concept of differential flatness is exploited, which is based on a parameterization of the system states and inputs by a flat output and its derivatives. In recent years the flatness-based concept (Fliess, Lévine, Martin, & Rouchon, 1995) was extended to linear and nonlinear distributed parameter systems (Knüppel & Woittennek, 2015; Meurer, 2013; Meurer & Kugi, 2009; Rudolph, 2003; Schörkhuber, Meurer, & Jüngel, 2013). The properties of a Riesz spectral systems are considered in Meurer (2011, 2013) to obtain a flatness-based design systematic taking into account spatially distributed inputs. This approach, which has been successfully applied in recent works considering single bending beams (Schröck, Meurer, & Kugi, 2010b, 2013), is extended to coupled structures. This contribution concentrates on the case of spatially varying parameters of the beam structure so that an early lumping approach is developed based on a Galerkin approximation in Section 3. The resulting feedforward control is amended by a state feedback controller in form of the linear quadratic regulator (LQR) to deal with disturbances and model uncertainties. The necessary state information is provided by a Kalman filter utilizing the in-domain measurements in terms of strain gauges. Herein the connection is shown between the Hankel singular values used for model order reduction and the choice of the weighting matrices in the LQR design and the process noise covariance matrix in the setup of the Kalman filter. Section 4 provides the experimental evaluation of the in-domain measurement as well as the performance of the two-degrees-of-freedom control approach. Finally, Section 5 concludes the paper with some final remarks. The main contributions of the paper can be summarized as follows:

f 𝑈𝑗,a (𝑡) = 𝑈0 + 𝑈𝑗,𝑘 (𝑡), 𝑘

b 𝑈𝑗,a (𝑡) = 𝑈0 − 𝑈𝑗,𝑘 (𝑡) 𝑘

(1)

generates a bending moment around the 𝑧2 -axis (Banks, Smith et al., 1996). Here, the superscripts refer to the location of the MFC-patch of the 𝑘th pair on the front (f) or on the back (b) of a beam. The voltage range of a MFC-patch is given by 𝑈𝑗f,b ∈ [−500, 1500] V (Smart Materials GmbH, 2014). To obtain a symmetric input range of 𝑈𝑗,𝑘 ∈ [−1000, 1000] V the constant supply voltage 𝑈0 = 500 V is applied. Remark 1. In general MFC-patches show nonlinear behavior due to hysteresis and creeping. These effects can be neglected if a compensation algorithm is used (see Kugi et al. (2006) and Schröck, Meurer, and Kugi (2010a)) implying linear properties of each MFC-patch. The equations of motion are derived by applying the variational formulation of the extended Hamilton’s principle taking into account nonconservative damping effects, i.e., 𝑡1 (

∫ 𝑡0

) 𝛿𝑊kin (𝑡) − 𝛿𝑊pot (𝑡) + 𝛿𝑊nc (𝑡) d𝑡 = 0

(2)

with the variational operator 𝛿 and 𝑊kin , 𝑊pot , 𝑊nc the kinetic energy, the potential energy, and the virtual work of the nonconservative forces (see, e.g., Nowacki (1975) and Reddy (1984)). The kinetic energy is stored in the mass of the 𝑟 beams and their tip masses, i.e., 𝑟 ( ∑ 𝑚𝑗 ( )2 𝜕𝑡 𝑤𝑗 (l , 𝑡) 𝑊kin (𝑡) = 𝑊𝑗,kin (𝑡) + 2 𝑗=1 (3) ) 𝐽𝑗 ( )2 + 𝜕𝑡 𝜕𝑧 𝑤𝑗 (l , 𝑡) , 2 where 𝐽𝑗 denotes the inertia of the tip mass. The kinetic energy of a single beam (see Schröck et al. (2010b)) is given by

1. Model-based motion planning and feedforward control are systematically developed and evaluated for a complex flexible benchmark structure consisting of elastically coupled beams with piezoelectric in-domain actuation.

𝑊𝑗,kin (𝑡) = 390

l ( )2 1 𝜇 (𝑧) 𝜕𝑡 𝑤𝑗 d𝑧 2 ∫0 𝑗

(4)

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

with

for 𝑧 ∈ (0, l ) and 𝑡 > 0. The boundary conditions are 𝑝

𝜇𝑗 (𝑧) = 𝜌𝑗 𝐴𝑗 +

𝑗 ∑

𝑤𝑗 = 0,

2𝜌𝑗,a𝑘 𝐴𝑗,a𝑘 𝛺𝑗,a𝑘 (𝑧).

𝑘=1

for 𝑧 = 0, 𝑡 > 0 and

[ ] 𝑚𝑗 𝜕𝑡2 𝑤𝑗 + k𝑗 (𝒘) = 𝜕𝑧 𝛬𝑗 (𝑧)𝜕𝑧2 𝑤𝑗 − 𝛾𝑗sv (𝑧)𝜕𝑡 𝜕𝑧 𝑤𝑗

Herein, 𝐴𝑗 = w𝑗 h𝑗 , 𝐴𝑗,a𝑘 = w𝑗,a𝑘 h𝑗,a𝑘 describe the cross section area defined by the width w𝑗 , w𝑗,a𝑘 as well as the thickness h𝑗 , h𝑗,a𝑘 (see also Fig. 1) and 𝜌𝑗 , 𝜌𝑗,a𝑘 characterize the density of the carrier and actuator, respectively. The spatial characteristics of the MFC-actuators are denoted by 𝛺𝑗,a𝑘 (𝑧) = ℎ(𝑧 − 𝑧𝑗,a ) − ℎ(𝑧 − 𝑧𝑗,a + l𝑗,a𝑘 ) 𝑘

𝑘

𝐽𝑗 𝜕𝑡2 𝜕𝑧 𝑤𝑗 = −𝛬𝑗 (𝑧)𝜕𝑧2 𝑤𝑗

𝑊pot (𝑡) =

(5)

with 𝒘(l ) = [𝑤1 (l ), … , 𝑤𝑟 (l )]𝑇 . The initial condition of the system at 𝑡 = 0 is given by

𝑊𝑗,pot (𝑡)

𝑗=1 𝑟−1 ∑ 𝑘𝑗,𝑗+1 ( )2 𝑤𝑗 (l , 𝑡) − 𝑤𝑗+1 (l , 𝑡) . + 2 𝑗=1

𝑤𝑗 = 𝑤𝑗,0 ,

(6)

l ( )2 1 𝛬 (𝑧) 𝜕𝑧2 𝑤𝑗 d𝑧 2 ∫0 𝑗 𝑝

+

𝑗 ∑

𝑘=1

(7)

l

∫0

(11e)

𝜕𝑡 𝑤𝑗 = 𝑤𝑗,1 .

Due to the spatial characteristics 𝛺𝑗,a𝑘 (𝑧) of the MFC-patches defined in (5) the system (11) has to be interpreted in a distributional sense leading to the introduction of the weak formulation (see also Banks, Smith et al. (1996), Schröck et al. (2013) and Wloka (1987)). The in-domain measurement is achieved by strain gauges, which provide the longitudinal strain averaged over the covered area. In terms of the Euler–Bernoulli theory the longitudinal strain is given by

Following Schröck et al. (2010b) the potential energy stored in a single beam is given by 𝑊𝑗,pot (𝑡) =

(11c)

at the free but coupled end 𝑧 = l . The effect of the coupling springs is taken into account by ) ( ⎧𝑘 𝑤 − 𝑤 , 𝑗=1 2 ⎪ 1,2 (1 ) ⎪𝑘𝑗,𝑗+1 𝑤𝑗 − 𝑤𝑗+1 (11d) 𝑗 ∈ {2, … , 𝑟 − 1} k𝑗 (𝒘) = ⎨ ) ( ⎪−𝑘𝑗−1,𝑗 𝑤𝑗−1 − 𝑤𝑗 , ( ) ⎪−𝑘 𝑗 = 𝑟, ⎩ 𝑟−1,𝑟 𝑤𝑟−1 − 𝑤𝑟 ,

with ℎ(⋅) the Heaviside function, l𝑗,a𝑘 the length and 𝑧𝑗,a the longitudinal 𝑘 position of the mounted patch. Analogously, the potential energy is given by the energy stored in the beams and the connecting springs 𝑟 ∑

(11b)

𝜕𝑧 𝑤𝑗 = 0,

𝛤𝑗,a𝑘 (𝑧)𝑈𝑗,𝑘 (𝑡)𝜕𝑧2 𝑤𝑗 d𝑧

𝜖𝑗,1 (𝑧, 𝑡) =

h𝑗 2

𝜕𝑧2 𝑤𝑗 (𝑧, 𝑡).

(12)

with 𝑝

𝛬𝑗 (𝑧) = 𝐸𝑗 𝐼𝑗 +

𝛤𝑗,a𝑘 (𝑧) =

𝑗 ∑

+ h𝑗,a𝑘 )

2𝛽𝑗,a𝑘 ,11 𝑒𝑠

Remark 2. In the following, the inertia of the tip mass 𝐽𝑗 is neglected since it does not show significant effect either in the static or in the transient behavior of the experimental setup.

(8)

2𝐸𝑗,a𝑘 𝐼𝑗,a𝑘 𝛺𝑗,a𝑘 (𝑧),

𝑘=1 𝐴𝑗,a𝑘 𝑎11 (h 𝑗,a𝑘 ,1 𝑗

𝛺𝑗,a𝑘 (𝑧).

(9) 3. Flatness of a coupled beam structure

These substitutions are derived under the assumption of an asymmetric voltage supply, the linear piezoelectric constitutive law and the underlying geometry. Herein, 𝐸𝑗 and 𝐸𝑗,a𝑘 denote the elasticity modulus of the carrier structure and the MFC-patches and 𝐼𝑗 = w𝑗 h𝑗3 ∕12 as well as 𝐼𝑗,a𝑘 = w𝑗,a𝑘 [(h𝑗 ∕2 + h𝑗,a𝑘 )3 − (h𝑗 ∕2)3 ] describe the area moment of inertia of the carrier and MFC-actuators. The parameters 𝑎11 , 𝛽𝑗,a𝑘 ,11 𝑗,a𝑘 ,1 and 𝑒𝑠 represent a piezoelectric parameter, the dielectric constant of the MFC-patch pair and the distance between neighboring electrodes on an MFC-patch. Damping is introduced by taking into account the variational form of the virtual work of nonconservative forces, i.e., 𝛿𝑊nc = −

𝑗=1

+

l

𝑟 ∑

∫

( 𝛾𝑗v (𝑧)𝜕𝑡 𝑤𝑗 (𝑧, 𝑡)𝛿𝑤𝑗 (𝑧, 𝑡)

0

)

𝛾𝑗sv (𝑧)𝜕𝑡 𝜕𝑧 𝑤𝑗 (𝑧, 𝑡)𝛿𝜕𝑧 𝑤𝑗 (𝑧, 𝑡)

A semi-analytic flatness-based approach is presented to solve the motion planning problem. For this purpose the weak formulation of the equations of motion (11) is used, which taking into account a test function 𝜙𝑗 (𝑧) ∈ 𝐻02 (0, l ) with 𝐻02 (0, l ) = {𝜙 ∈ 𝐻 2 (0, l )|𝜙(0) = 𝜕𝑧 𝜙(0) = 0}. The weak formulation also enables us to obtain a Galerkin approximation of the distributed parameter system. 3.1. Weak formulation and Galerkin approximation The weak form is deduced by the multiplication of (11) with a test function 𝜙𝑗 (𝑧), chosen from the set of eigenfunctions of the undamped system, followed by an integration by parts considering the boundary conditions (11b) and (11c). This yields

(10)

⟨𝜇𝑗 𝜕𝑡2 𝑤𝑗 , 𝜙𝑗 ⟩ + ⟨𝛬𝑗 𝜕𝑧2 𝑤𝑗 , 𝜕𝑧2 𝜙𝑗 ⟩ + ⟨𝛾𝑗v 𝜕𝑡 𝑤𝑗 , 𝜙𝑗 ⟩

d𝑧.

+ ⟨𝛾𝑗sv 𝜕𝑡 𝜕𝑧 𝑤𝑗 , 𝜕𝑧 𝜙𝑗 ⟩ + 𝑚𝑗 𝜕𝑡2 𝑤𝑗 (l )𝜙𝑗 (l )

Herein, two distinct damping models are considered. Viscous damping v + ∑𝑝𝑗 2𝛾 v 𝛺 is introduced by 𝛾𝑗v (𝑧) = 𝛾𝑗,c 𝑗,a𝑘 𝑗,a𝑘 (𝑧). In addition, structural 𝑘=1 sv + ∑𝑝𝑗 2𝛾 sv 𝛺 v v damping is defined as 𝛾𝑗sv (𝑧) = 𝛾𝑗,c 𝑗,a𝑘 𝑗,a𝑘 (𝑧). Herein 𝛾𝑗,c , 𝛾𝑗,a𝑘 , 𝑘=1 sv and 𝛾 sv denote the damping coefficients of 𝛾𝑗,c the carrier and the MFC𝑗,a𝑘 actuators, respectively. It has to be mentioned, that the combination of viscous and structural damping is required to match the behavior observed in experiments. The evaluation of (2) with (3)–(10) and the successive partial integration combined with the application of the fundamental lemma of variational calculus leads to the equations of motion of the coupled structure [ ] 0 = 𝜇𝑗 (𝑧)𝜕𝑡2 𝑤𝑗 + 𝛾𝑗v (𝑧)𝜕𝑡 𝑤𝑗 − 𝜕𝑧 𝛾𝑗sv (𝑧)𝜕𝑡 𝜕𝑧 𝑤𝑗 𝑝

𝑗 [ ] ∑ + 𝜕𝑧2 𝛬𝑗 (𝑧)𝜕𝑧2 𝑤𝑗 + 𝜕𝑧2 𝛤𝑗,a𝑘 (𝑧)𝑈𝑗,𝑘 ,

𝑝

+ k𝑗 (𝒘)𝜙𝑗 (l ) =

𝑗 ∑

(13)

− ⟨𝛤𝑗,a𝑘 𝑈𝑗,𝑘 , 𝜕𝑧2 𝜙𝑗 ⟩, 𝑘=1 2

where ⟨⋅, ⋅⟩ denotes the 𝐿 (0, l ) inner product. This formulation allows us to determine the finite dimensional Galerkin approximation of (11) by defining a specific set of linear independent test functions {{𝜙𝑗,𝑖 }𝑟𝑗=1 }𝑁 𝑖=1 with 𝑁 ∈ N by assuming the expansion 𝑤𝑗 (𝑧, 𝑡) =

𝑁 ∑

𝑞𝑗,𝑖 (𝑡)𝜙𝑗,𝑖 (𝑧),

(14)

𝑖=1

which results in a separation of the spatial and temporal variable. Substitution of (14) into (13) for 𝑗 = 1, 2, … , 𝑟 and evaluating the arising integrals leads to

(11a)

̈ + 𝐶 𝒒(𝑡) ̇ + 𝐾𝒒(𝑡) = 𝐿𝒖(𝑡), 𝑀 𝒒(𝑡)

𝑘=1

391

(15)

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

with 𝒒(𝑡) = [𝑞1,1 (𝑡) … 𝑞1,𝑁 (𝑡) … 𝑞𝑟,1 (𝑡) … 𝑞𝑟,𝑁 (𝑡)]𝑇 ∈ R𝑟𝑁 and 𝒖(𝑡) = [𝑈1,1 (𝑡) … 𝑈1,𝑝1 (𝑡) … 𝑈𝑟,1 (𝑡) … 𝑈𝑟,𝑝𝑟 (𝑡)]𝑇 ∈ R𝑛𝑢 denoting the vector of generalized coordinates and inputs, respectively. The elements of the mass matrix 𝑀, the damping matrix 𝐶, the stiffness matrix 𝐾, and the input matrix 𝐿 are summarized in the Appendix. By introducing the state vector 𝒙(𝑡) = [𝒒 𝑇 (𝑡) 𝒒̇ 𝑇 (𝑡)]𝑇 ∈ R2𝑟𝑁 system (15) can be rewritten as ̇ 𝒙(𝑡) = 𝐴𝒙(𝑡) + 𝐵𝒖(𝑡), 𝑡 > 0, 𝒙(0) = 𝒙0 , with [ ] [ ] 0 𝐼 0 𝐴= , 𝐵 = , −𝑀 −1 𝐾 −𝑀 −1 𝐶 𝑀 −1 𝐿

where the output matrix is chosen as 𝐻 = 𝐼 to later enable the state observer design. Both symmetric and positive definite matrices are partitioned into pairs of block matrices [ ] [ ] c,p c,12 o,p o,12 c = , o = , (23) c,21 c,v o,21 o,v whereby the second index denotes the connection of the Gramians to the position or the velocity contributions. Under the constraint that both Gramians are regular the system is controllable and observable, i.e., there exists a matrix , that transforms ̃ c,p =  ̃ o,p and  ̃ c,v =  ̃ o,v the system into a balanced form where 

(16) (17)

holds (Reis & Stykel, 2008). The construction of  requires at first the decomposition of the sub-Gramians by the Cholesky factorization

where 𝐼 denotes the identity matrix. In addition, the system is characterized by the two output vectors (18)

𝑇 c,p = c,p c,p ,

where 𝒚 pi ∈ R𝑛pi denotes the deflection at the points of interest of the structure 𝒛pi ∈ R𝑛pi , for which the motion is controlled (see also Fig. 1), and 𝒚 pm ∈ R𝑛pm denotes the measured strain (12) of the points of measurement 𝒛pm ∈ R𝑛pm . Considering (14) the output matrices are given by

𝑇 o,p o,p ,

𝒚 pi (𝑡) = 𝐻pi 𝒙(𝑡),

⎡𝒉𝑇 1,pi|pm

⎢ 𝐻pi|pm = ⎢ ⎢ ⎣

0 0

𝒚 pm (𝑡) = 𝐻pm 𝒙(𝑡),

0 ⋱ 0

0 0 𝒉𝑇𝑛

pi|pm ,pi|pm

o,p =

𝑇 o,p c,p = 𝑈p 𝛴p 𝑉p𝑇 ,

(20)

(25)

with (21)

−1

p = c,p 𝑉p 𝛴p 2 ,

3.2. Model order reduction and spectral representation

−1

v = c,v 𝑉v 𝛴v 2 ,

The dimension of the state vector dim(𝒙) = 2𝑟𝑁 rises with the number of beams 𝑟 and the number of used basis functions 𝑁 defining, together with the choice of test function, the precision of the Galerkin approximation. Hence, the computational effort especially for a implementation of a state feedback controller with observer on a real time system accumulates with larger number of beams. Applying model order reduction (MOR) techniques enables us to reduce the dimension while retaining a sufficiently accurate behavior of the resulting reduced order system compared to the full order dynamics. For mechanical systems there are four major approaches for MOR: modal truncation, balanced truncation, proper orthogonal decomposition (POD), and methods based on Krylov subspaces. The paper Antoulas (2005) provides a summary with detailed information about the basics and common variations of these methods. Due to the advantages of preserving the system’s stability and providing an a priori computable error bound of the reduction error, a particular form of the balanced truncation is used. The main idea, see, e.g., Antoulas (2005), of the general concept is to identify and truncate the states of the system, which are simultaneously difficult to reach and to observe. For the simultaneous identification of these properties the controllability and the observability Gramians c , o , respectively, ̃c =  ̃o are used and transformed into a balanced form, in which  holds. The balanced form relies on the determination of the Hankel singular values, which simultaneous quantify the controllability and the observability of each state. This quantification can also be used for the adjustment of a state controller and observer as will be demonstrated later. Addressing the special structure of the input and system matrix (17) arising from the second order system (15), the so-called second order balanced truncation is used to preserve the physical meaning of the representation (Chahlaoui, Lemonnier, Vandendorpe, & Van Dooren, 2006; Reis & Stykel, 2008). The controllability and the observability Gramians can be calculated by solving the associated Lyapunov equations 𝐴c + c 𝐴𝑇 = −𝐵𝐵 𝑇 , 𝐴𝑇 o + o 𝐴 = −𝐻 𝑇 𝐻,

𝑇 o,v c,v = 𝑈v 𝛴v 𝑉v𝑇 .

This decomposition provides directly the Hankel singular values 𝜁p|v,𝑖 in terms of the elements of the diagonal matrix 𝛴p|v . With these preparations the transformation matrix and its inverse are given by [ ] ] [  −1 0  0 ,  −1 = p = p (26) 0 v 0 v−1

(19)

with

[ ] 𝒉𝑇𝑓 ,pi = 𝜙𝑗,1 (𝑧𝑓 ,pi ) … 𝜙𝑗,𝑁 (𝑧pi,𝑓 ) , ] h𝑗 [ 2 ) … 𝜕𝑧2 𝜙𝑗,𝑁 (𝑧pm,𝑔 ) . 𝜕 𝜙 (𝑧 𝒉𝑇𝑔,pm = 2 𝑧 𝑗,1 pm,𝑔

(24)

𝑇 o,v = o,v o,v .

The associated lower triangular matrices are multiplied and a singular value decomposition is applied

𝟎𝑇 ⎤

⎥ ⋮ ⎥, 𝟎𝑇 ⎥ ⎦

𝑇 c,v = c,v c,v ,

−1

𝑇 p−1 = 𝛴p 2 𝑈p𝑇 o,p ,

(27)

−1

𝑇 v−1 = 𝛴v 2 𝑈v𝑇 o,v .

The reduction is achieved by considering only the 𝑟𝑁r position and the 𝑟𝑁r velocity states with the highest associated Hankel singular values 𝜁p|v,1 to 𝜁p|v,𝑟𝑁r and neglecting all other states. This can be done by a partition of the transformation sub-matrices [ ] −1 p|v,1 ] [ −1 , (28) p|v = p|v,1 p|v,2 , p|v = −1 p|v,2 −1 =  −1 ∈ R𝑟𝑁r ×𝑟𝑁 are chosen. where ̃p|v = p|v,1 ∈ R𝑟𝑁×𝑟𝑁r and ̃p|v p|v,1 −1 ̃ ̃ ̃ ̃ Introducing  = diag{p , v } and  = diag{̃p−1 , ̃v−1 } allows us to determine the reduced order system

𝒙̇ r (𝑡) = 𝐴r 𝒙r (𝑡) + 𝐵r 𝒖(𝑡), 𝒚 pi (𝑡) = 𝐻r,pi 𝒙r (𝑡),

𝑡 > 0,

𝒙r (0) = 𝒙r,0 ,

𝒚 pm (𝑡) = 𝐻r,pm 𝒙r (𝑡),

𝑡 ≥ 0,

(29)

and the reduced order state is given by 𝒙r = ̃ −1 𝒙 ∈ R2𝑟𝑁r with 𝐴r = ̃ −1 𝐴̃ ,

𝐵r = ̃ −1 𝐵,

𝐻r,pi|pm = 𝐻pi|pm ̃ .

(30)

The applied damping model leads to pairs of complex conjugate eigenvalues of the full system 𝝀 ∈ C2𝑟𝑁 with ℜ(𝝀) < 𝟎. As stated in Reis and Stykel (2008, Corollary 3.2) the applied second order balance truncation preserves the stability in the case of a symmetric system (15) with a positive definite mass matrix 𝑀, damping matrix 𝐶, and stiffness matrix 𝐾 as it is given here. This statement is confirmed by numerical computations. The eigenvalues of the reduced order system 𝜆r,𝑙 ∈ C with 𝑙 ∈ L = {1, … , 2𝑟𝑁r } fulfill ℜ(𝜆r,𝑙 ) < 0 for the system parameters listed in Tables 1 and 2. The spectral representation of the system is obtained by the similarity transformation 𝒙r (𝑡) = 𝑉 𝒛(𝑡) with 𝑉 = [𝒗r,1 , … , 𝒗r,2𝑟𝑁r ] the matrix of eigenvectors of 𝐴r . With this (29) reduces to 𝑧̇ 𝑙 (𝑡) = 𝜆r,𝑙 𝑧𝑙 (𝑡) + 𝒃𝑇𝛥,𝑙 𝒖(𝑡),

𝑙∈L

(31)

where 𝒃𝑇𝛥,𝑙 denotes the 𝑙th row of the transformed input matrix 𝐵𝛥 = 𝑉 −1 𝐵r .

(22) 392

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

steady state the input parameterization (36) reduces to 𝒖st = 𝝃 st , the required relation between 𝝃 ∗st and 𝒚 ∗st is given by

3.3. Construction of a flat output Applying the Laplace transformation to the spectral state representation (31) yields 𝑧̂ 𝑙 (𝑠) =

1 1 1 ̂ =− ̂ 𝒃𝑇 𝒖(𝑠) 𝒃𝑇 𝒖(𝑠), 𝑠 − 𝜆𝑙 𝛥,𝑙 𝜆𝑙 1 − 𝑠 𝛥,𝑙 𝜆

−1 𝒙r,st = −𝐴−1 r 𝐵r 𝒖st = −𝐴r 𝐵r 𝝃 st , 𝒚 pi,st = 𝐻r,pi 𝒙r,st { ( )−1 −𝐻r,pi 𝐴−1 𝒚 , 𝑛pi = 𝑛𝑢 r 𝐵r ⇒ 𝝃 st = ( )† pi,st −𝐻r,pi 𝐴−1 𝐵 𝒚 , 𝑛pi ≠ 𝑛𝑢 , r r pi,st

(32)

𝑙

̂ with 𝒛(𝑠) unity yields

̂ 𝒛(𝑡) and 𝒖(𝑠)

𝒖(𝑡). Extending the term 1∕(1 −

∏ (1 − 𝜆𝑠 ) 1 𝑖∈L,𝑖≠𝑙 𝑖 ̂ 𝑧̂ 𝑙 (𝑠) = − 𝒃𝑇 𝒖(𝑠) ∏ 𝑠 𝜆𝑙 ) 𝛥,𝑙 𝑖∈L (1 −

𝑠 𝜆𝑙

𝑛pi ×𝑛𝑢 where 𝐻r,pi 𝐴−1 . Inserting the desired values for 𝒚 ∗st at r 𝐵r ∈ R 𝑡 = 0 and 𝑡 = 𝑇 directly provides the start and end value of the flat output. Substitution of (38) into the input parameterization (36) defines the feedforward control law and the state parameterization defines the associated transient state trajectory.

) by

(33)

𝜆𝑖

The substitution ) ∏( 𝑠 ̂ ̂ ̂ (𝑠)𝝃(𝑠) ̂ = 𝒖(𝑠) 1− 𝝃(𝑠) = 𝐷 𝑢 𝜆 𝑖 𝑖∈L leads to 𝑧̂ 𝑙 (𝑠) = − =−

( ) 1 ∏ 𝑠 𝑇 ̂ 1− 𝒃 𝝃(𝑠) 𝜆𝑙 𝑖∈L,𝑖≠𝑙 𝜆𝑖 𝛥,𝑙 1 ̂ ̂ 𝐷 (𝑠)𝒃𝑇𝛥,𝑙 𝝃(𝑠). 𝜆𝑙 𝑧

Remark 3. The presented flatness-based approach to motion planning and feedforward control can be directly extended to incorporate an arbitrary number of modes corresponding to infinite-dimensional continuum model (11). For this, the index set L has to be replaced by a countable infinite set, typically N referring to the infinite number of modes. The determined parametrizations (34) or (35), respectively, then represent infinite products or differential operators of infinite order, respectively. Motion planning and the convergence of these expressions can be addressed similar to (38) by taking into account smooth functions of a certain Gevrey order, see, e.g., Meurer (2011, 2013) and Schröck et al. (2013).

(34a)

(34b)

Hence, (34) can be interpreted as an input and state parameterization ̂ ̂ (𝑠) and 𝐷 ̂ (𝑠) as in terms of 𝝃(𝑠) in the operator domain. Rewriting 𝐷 𝑢 𝑧 polynomials yields ̂ (𝑠) = 𝐷 𝑢

|L| ∑

𝑢 𝑐𝑖 𝑠𝑖 ,

̂ (𝑠) = 𝐷 𝑧

|L|−1 ∑

𝑧 𝑐𝑖 𝑠𝑖

Remark 4 (Advantages of Flatness-Based Motion Planning). Contrary to many other techniques including input shaping or optimal control the presented flatness approach allows us (i) to compute the correct input trajectory that is required to achieve rest-to-rest motion along a predefined path; (ii) to utilize the state parameterization, which provides us with the time evolution of the nominal (desired) state variables, for the design of the tracking controller (see Fig. 2); (iii) to obtain the state and input parameterization by only algebraic computations without the need to integrate ODEs or to solve an optimization problem; (iv) to a priori compute and analyze the desired state evolution and thus to iteratively reassign the flat output trajectories to fulfill user-defined conditions concerning the transient system behavior.

(35)

𝑖=0

𝑖=0

with |L| = max L = 2𝑟𝑁r . Recalling that 𝑠 corresponds to time differentiation, so that (34a) with (35) leads to 𝒖(𝑡) =

|L| ∑

𝑢

𝑐𝑖 𝝃 (𝑖) (𝑡), 𝑧𝑙 (𝑡) = −

𝑖=0

|L|−1 1 ∑ 𝑧 𝑇 (𝑖) 𝑐 𝒃 𝝃 (𝑡). 𝜆𝑙 𝑖=0 𝑖 𝛥,𝑙

(39)

(36)

Similarly 𝒙r (𝑡) = 𝑉 𝒛(𝑡) can be computed in terms of 𝝃(𝑡) and its derivatives, so that 𝝃(𝑡) denotes a flat output for (29). 3.4. Motion planning and feedforward control

This systematic approach enables the precalculation of the state and input trajectories. Hence, additional requirements of the experimental setup such as input saturations, restricted deflections, etc. can be taken into account by modifying 𝝃 ∗ (𝑡).

The objective of motion planning is to achieve a desired rest-torest motion of the deflection profiles of the beam structure within a finite time span 𝑡 ∈ [0, 𝑇 ] starting at an initial steady state deflection. Therefore, the input and state parameterization (36) is used to compute ̂ (𝑠) and a corresponding feedforward control. The series expansions of 𝐷 𝑧 ̂ (𝑠) require that the flat output fulfills 𝝃(𝑡) ∈ 𝐶 |L| (R). This requirement 𝐷 𝑢 can be satisfied by the usage of Fourier series, polynomials or splines. Subsequently, a locally non-analytic smooth function 𝛩𝑔,𝑇 (𝑡) is used defined by ⎧0 , ⎪ ⎪1 , 𝛩𝑔,𝑇 (𝑡) = ⎨ 𝑡 ⎪ ∫0 𝜃𝑔,𝑇 (𝜏)𝑑𝜏 , ⎪ 𝑇 ⎩ ∫0 𝜃𝑔,𝑇 (𝜏)𝑑𝜏

3.5. Observer design and feedback tracking control In order to deal with unknown disturbances, model uncertainties, and inaccuracies in the hysteresis compensation for the MFC-actuators the feedforward control is amended by a state feedback control strategy in the so-called two-degrees-of-freedom control approach. Fig. 2 displays the corresponding block diagram. Due to the fact, that the measurements are restricted to a few spatially discrete positions on the interconnected beam structure an observer has to be used to estimate the state vector of the system.

𝑡≤0 𝑡≥𝑇

(37)

𝑡 ∈ (0, 𝑇 ).

3.5.1. Kalman filter The usage of the strain gauges comes with the cost, that the provided electrical signal is noisy. To reduce the effect of noise to address the implementation on a real time control device the discrete linear Kalman filter is used. The discrete system description of (29) is given by

Herein, 𝜃𝑔,𝑇 (𝑡) = 0 for 𝑡 ∉ (0, 𝑇 ) and 𝜃𝑔,𝑇 (𝑡) = exp(−([1 − 𝑇𝑡 ] 𝑇𝑡 )−𝑔 ) for 𝑡 ∈ (0, 𝑇 ) . With this, the desired flat output trajectory is build by ( ) 𝝃 ∗ (𝑡) = 𝝃 ∗0 + 𝝃 ∗𝑇 − 𝝃 ∗0 𝛩𝑔,𝑇 (𝑡) (38) with 𝝃 ∗0 and 𝝃 ∗𝑇 defining the start and end value of the flat output, respectively. In the next step, these values have to be expressed in term of the points of interest (18) as the desired initial and final deflection profile 𝒚 ∗pi (0) and 𝒚 ∗pi (𝑇 ), respectively. The specification of single points of the structure is sufficient to accomplish a desired profile since a restto-rest motion is the objective, so that in the start and in the end all transients have to be settled and the profile is given by the specified points w.r.t. geometric constraints. Taken into account, that in the

𝒙𝑘 = 𝐴d 𝒙𝑘−1 + 𝐵d 𝒖𝑘−1 + w𝑘−1 , 𝑘 > 1, 𝒙r (0) = 𝒙r,0 , 𝒚 pi,𝑘 = 𝐻r,pi 𝒙𝑘 ,

𝒚 pm,𝑘 = 𝐻r,pm 𝒙𝑘 + v𝑘 ,

(40)

with 𝒙𝑘 = 𝒙r (𝑘𝑡d ), 𝒖𝑘 = 𝒖(𝑘𝑡d ) and 𝒚 pi|pm,𝑘 = 𝒚 pi|pm (𝑘𝑡d ) where 𝑡d denotes the sample time. In addition, the system is extended by process noise w𝑘−1 ∼  (𝟎, 𝑄o ) and measurement noise v𝑘 ∼  (𝟎, 𝑅o ). For these, white noise signals are assumed with 𝟎 expected values and covariance 393

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

Fig. 2. Block diagram of the two-degrees-of-freedom control scheme, with 𝛴 ∗ denoting the flatness-based trajectory generation and feedforward control as well as 𝛴tr representing the tracking controller.

matrices 𝑄o and 𝑅o . The matrices in (40) follow from Åström and Wittenmark (2011) as

Fig. 3. Experimental setup of the interconnected flexible beam structure.

𝑘𝑡d

𝐴d = 𝑒𝐴r 𝑡d , 𝐵d =

∫

𝑒𝐴r (𝑘𝑡d −𝜏) 𝐵r d𝜏, 𝐻d,pi|pm = 𝐻r,pi|pm .

(41) 3.5.3. Choice of weighting matrices In the practical implementation the covariance matrices of the Kalman filter are used as tuning factors. For the measurement noise due to strain gauges (type HBM 1-Ly48-6/350), the covariance is chosen as 𝑅o = 20𝐼. The choice for the process noise covariance matrix is based on the Hankel singular values 𝜁p|v,𝑖 of (25), which provide information about the observability of each state. Due to the partitioning into position and velocity related states the matrix is set in the same manner

(𝑘−1)𝑡d

The Kalman filter can be split into in the extrapolation phase ̂ 𝑘−1 + 𝐵d 𝒖𝑘−1 , 𝑃𝑘− = 𝐴d 𝑃𝑘−1 𝐴𝑇d + 𝑄o 𝒙̂ − 𝑘 = 𝐴d 𝒙

(42a)

and the update phase ̂− 𝒙̂ 𝑘 = 𝒙̂ − 𝑘 + 𝐿𝑘 (𝒚 pm,𝑘 − 𝐻d,pm 𝒙 𝑘 ), 𝑃𝑘 = (𝐼 − 𝐿𝑘 𝐻d,pm )𝑃𝑘− , 𝐿𝑘 =

𝑇 𝑇 𝑃𝑘− 𝐻d,pm (𝐻d,pm 𝑃𝑘− 𝐻d,pm

(42b) −1

+ 𝑅o ) ,

𝑄o = 𝑞o,𝑡 diag{𝒒 o,p , 𝒒 o,v }, [ ] 𝒒 𝑇o,p|v = 𝑞o,p|v 𝜁p|v,1 … 𝜁p|v,𝑁r

where 𝑃𝑘 denotes the error covariance matrix with the expected start value of 𝑃0 = 𝐸(𝒙̃ o,0 𝒙̃ 𝑇o,0 ) and the observer error 𝒙̃ o,𝑘 = 𝒙𝑘 − 𝒙̂ 𝑘 (Gelb, 1974; Lewis, Xie, & Popa, 2008). In case of stable observer error dynamics the error covariance approaches to a steady state lim𝑘→∞ 𝑃𝑘− = 𝑃, that is given by the solution of an algebraic Riccati equation. Hence, the Kalman gain becomes constant so that 𝑇 𝑇 𝐿 = 𝑃𝐻d,pm (𝐻d,pm 𝑃𝐻d,pm + 𝑅o )−1 and the update phase simplifies to ̂ 𝑘−1 ). 𝒙̂ 𝑘 = 𝒙̂ − 𝑘 + 𝐿(𝒚 pm,𝑘−1 − 𝐻d,pm 𝒙

∑𝑁 ∑𝑁 with 𝑞o,𝑡 = 50, 𝑞o,p = 5 × 103 ∕ 𝑖=1r 𝜁p,𝑖 , and 𝑞o,v = 10∕ 𝑖=1r 𝜁v,𝑖 . The numeric factors are determined during hardware-in-the-loop tests. The usage of the Hankel singular values is advantageous, due to the reduction of complexity of the possible combinations of the elements of 𝑄o . Similar to the observer design the second order representation allows a straightforward choice of the elements of the weighting matrices in form

(42c)

3.5.2. Linear quadratic regulator To ensure that the desired trajectory 𝒙∗𝑘 is reached a linear quadratic regulator (LQR) is used to minimize the tracking error 𝒙̃ c,𝑘 = 𝒙∗𝑘 − 𝒙̂ 𝑘 . Therefore, the quadratic cost function (Bryson, 2002) for an infinite horizon is used ∞ ( ) ∑ 𝒙̃ 𝑇c,𝑘 𝑄c 𝒙̃ c,𝑘 + 𝒖𝑇c,𝑘 𝑅c 𝒖c,𝑘 . (43) 𝐽 (𝒙̃ c,𝑘 , 𝒖c,𝑘 ) =

𝑄c 𝑇 𝒒 c,p|v

= 𝑞c,𝑡 diag{𝒒 c,p , 𝒒 c,v }, 𝑅c = 0.01𝐼, [ ] = 𝑞c,p|v 𝜁p|v,1 … 𝜁p|v,𝑁r

(47)

∑𝑁 ∑𝑁 with 𝑞c,𝑡 = 2 × 102 , 𝑞c,p = 10∕ 𝑖=1r 𝜁p,𝑖 , and 𝑞c,v = 1 × 10−3 ∕ 𝑖=1r 𝜁v,𝑖 , whereby the numeric factors are similarly determined by hardware-inthe-loop tests. The required Hankel singular values provided by the second order balance truncation are listed in Table 3.

𝑘=1

The positive definite matrix 𝑅c and positive semi-definite matrix 𝑄c can be interpreted as weighting matrices for the input as well as for the error state. The cost function has to be minimized under the constraint of the open loop dynamics. Recalling the separation theorem the control design is independent of the observer. Hence, the dynamics can be deduced as 𝒙∗𝑘 = 𝐴d 𝒙∗𝑘−1 + 𝐵d 𝒖∗𝑘−1 and 𝒖𝑘 = 𝒖∗𝑘 + 𝒖c,𝑘 with 𝒖c,𝑘 = 𝐾 𝒙̃ c,𝑘 , where 𝐾 denotes the state feedback gain. Following standard arguments, the control law simplifies for 𝑘 → ∞ to 𝒖c,𝑘 = 𝐾 𝒙̃ c,𝑘 , with 𝐾 = (𝐵d𝑇 𝑆𝐵d + 𝑅c )−1 𝐵d 𝑆𝐴d

(46)

4. Experimental results The considered experimental setup shown in Fig. 3 consists of 𝑟 = 3 coupled beams. The beams are composed of two layers of carbon fibers bond by epoxy glue. The geometric extensions of the beams and the embedded MFC-actuators (type M-8525-P1 (Smart Materials GmbH, 2014) characterized by an electrode spacing of 𝑒𝑠 = 5 × 10−4 m) are summarized in Table 1. In contrast to the other beams the second beam is equipped with two MFC-patch pairs.

(44)

and 𝑆 = lim𝑘→∞ 𝑆𝑘 is given by the solution of an algebraic Riccati equation. To deal with remaining static errors, that are, e.g. caused by residual parts of hysteresis (see Remark 1) the controller is extended by an integral part ∑ 𝒖c,𝑘 = 𝐾 𝒙̃ c,𝑘 + 𝑘I 𝐻d,pi 𝒙̃ c,𝑘 (45)

4.1. Parameter identification Table 2 displays the identified material parameters, which differ between beams, due to the manufacturing process. Their determination follows two steps to sequentially determine disjoint subsets of the in overall 37 system parameters (cf. Table 2). In the first step, the springs are detached so that the beams can be considered decoupled. Based on step responses one the hand the

𝑘

weighting the tracking error at the points of interest by the gain value 𝑘I . 394

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398 Table 2 Set of identified system parameters with 𝑟 = 3 beams.

Table 1 Dimensions, points of interest and measurements positions in mm.

𝐿2 –norm of the difference between measured and simulated deflection 𝒚 pm (𝑡) at 𝒛pm ∈ R6 is minimized during transients to identify 𝜌𝑗 , 𝛾𝑗v , 𝛾𝑗sv , v , and 𝛾 sv for each beam 𝑗 ∈ 1, 2, 3 and each actuator a . 𝑚𝑗 , 𝜌𝑗,a𝑘 , 𝛾𝑗,a 𝑘 𝑗,a𝑘 𝑘 On the other hand the steady stated difference is minimized to obtain 𝐸𝑗 , 𝐸𝑗,a𝑘 and 𝑎11 ∕𝛽 . This separation is motivated by the fact, that 𝑗,a𝑘 ,1 𝑗,a𝑘 ,11 the first parameter set has a major impact on the frequency, amplitude and phase of the resulting oscillations while the second parameter set implies the steady state value. In the second step, the springs are attached and the steady state difference between measured and simulated deflection at 𝑧pm,𝑗 , 𝑗 = 1, … , 6 is minimized by adjusting the coupling spring constants 𝑘1,2 and 𝑘2,3 .

Table 3 Hankel singular values in ×10−2 .

4.2. Galerkin approximation and model order reduction The test functions {{𝜙𝑗,𝑖 }3𝑗 }𝑁 𝑖 with 𝑁 = 10 utilized for the evaluation of the second order system matrices (see the Appendix) are chosen as the eigenfunctions of a uniform undamped Euler–Bernoulli beam in clamped-free configuration. The full order system is of dimension 2𝑟𝑁 = 60. The usage of balanced truncation with 𝑁r = 5 leads to the dimension 2𝑟𝑁r = 30 of the reduced order system, which covers 82% of the full order system’s energy considering the associated Hankel singular values listed in Table 3. This choice constitutes a compromise between high accuracy regarding the behavior of the experimental system, and computational effort, required by the implementation on a real time system. The resulting reduced order system satisfy the requirement of a high accuracy regarding the behavior of the experimental system, which can be seen later on in the experimental results, as well as the requirement of a reduced computational effort, that is required to implement the presented approaches on a real time system. For the implementation and recording of measurements a dSPACE MicroLabBox (dSPACE GmbH Product Information, 2017) with sample time 𝑡d = 2 × 10−4 s is used.

In addition Fig. 4 provides a comparison between the direct laser measurements 𝒚 pi,L (𝑡), which take place at the tip of each beam 𝑧pi,𝑖 with 𝑖 ∈ {1, 3, 4} (see Fig. 1 and Table 1), and the respective observer output 𝒚̂ pi (𝑡) driven by the measurements of the strain gauges located at 𝑧pm,𝑚 ≠ 𝑧pi,𝑖 ∀𝑖, 𝑚. As it can be seen, the observed output coincides with the laser measurements, which indicates a good performance of the observer and the underlying MOR with 𝑁r = 5. This is crucial for the implementation of the presented state feedback controller.

4.3. Step response

4.4. Feedforward control

To illustrate the dynamic of the interconnected structure Fig. 4 shows the response of the transversal deflection of the points of interest located near the beam tip’s 𝑦pi,𝑖 = 𝑤(𝑧pi,𝑖 , 𝑡) with 𝑖 ∈ {1, 3, 4} for a voltage step [ 𝒖(𝑡) = 817

−770

−807

]𝑇 781 V ℎ(𝑡 − 1 s).

To analyze the performance of the presented feedforward control a rest-to-rest motion connecting 𝒚 0 = 𝟎 m and 𝒚 𝑇 = 𝒚 pi (𝑇 ) is considered. The corresponding steady state is the same as in the step response scenario, so that in the steady state the same voltages emerge as before. The considered transition time of the feedforward control defined by (36) and (38) is chosen as 𝑇 = 0.9 s with a 𝛼 = 1.9. The resulting voltage trajectory 𝒖∗ (𝑡) is depicted by the blue line in the second row of Fig. 5. In the first row the associated desired reference trajectory 𝒚 ∗ (𝑡), computed a priori by evaluating the state parameterization (36) at the points of interest 𝒛pi ∈ R𝑛pi , is displayed as green line as well as the

(48)

These voltage values correspond to a steady state deflection of 𝒚 pi (𝑇 ) = [−8, 1, 4, −6]𝑇 × 10−3 m. The applied voltage step leads to rather large oscillations of all three beams with amplitudes between [4, 6] × 10−3 m. Due to the beam configuration, the second beam endures the highest relative overshoot of 150% w.r.t. the steady state values. The settling time of the oscillations is about 18 s. 395

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

Fig. 4. Step response for the voltage input (48): Transversal deflection of the beam tips 𝒚 pi (𝑡) = 𝒘(𝑧pi,𝑖 , 𝑡) with 𝑖 ∈ {1, 3, 4}. For comparison purposes laser measurements 𝒚 pi,L (𝑡) are opposed with the respective results obtained from the presented observer (42) driven by the strain gauge measurements 𝒚̂ pi (𝑡).

Fig. 5. Experimental results of the flexible beam structure. The first row presents the controlled outputs (points of interest) of the open loop 𝒚 pi,ff (𝑡) and the closed loop 𝒚 pi,tr (𝑡) behavior as well as the desired reference trajectory 𝒚 ∗pi (𝑡). The second row shows the associated input signals.

actual resulting trajectory 𝒚 pi,ff (𝑡) plotted as blue line. Compared to the step response the feedforward control results in a far smaller oscillation of the transversal deflection on all four points of interest after 𝑡 ≈ 2.5 s, where the steady voltage is reached. The largest relative overshooting can be seen by the second point of interest 𝑖 = 2 with an amplitude of 10%, which is, however, rapidly decaying.

of the remaining oscillation to 2 s. Furthermore, stationary accuracy is accomplished. Fig. 6 presents the experimental results of the flexible structure under disturbances. At 𝑡 = 1 s a pulse hammer applied to the second beam excites the complete structure. The first row displays the observed beam deflections of the open loop 𝒚 pi,ff (𝑡) in blue as well as the deflections in the closed loop system 𝒚 pi,tr (𝑡), plotted as dashed red line, at the points of interest. In the second row the associated closed loop input signal 𝒖tr (𝑡) is depicted as dashed red line. Compared to the uncontrolled structure, which displays an oscillation characterized by amplitudes larger than 0.03 m and a decay time of 25 s, the extended LQ controller displays a good performance by damping the oscillation and ensuring the desired deflection within 4.2 s. In summary, the experimental results show that the motion planning problem of the flexible interconnected structure can be solved by the presented model-based two-degrees-of-freedom control approach.

Remark 5 (Online Implementation of Feedforward Control). The proposed flatness-based motion planning and feedforward control is capable to almost completely achieve the desired rest-to-rest motion without the need to perform intensive computations or optimization. The approach is real time capable, is implemented in the real time processing unit and allows us to adjust the motion planning online. 4.5. Two-degrees-of-freedom control The application of the two-degrees-of-freedom control approach results in a further improvement of the performance as depicted by the measured deflections 𝒚 pi,tr (𝑡) by means of the red dash-dotted lines in the first row of Fig. 5. The tracking controller reduces the decay time

5. Conclusions In this contribution a flatness-based feedforward control approach is introduced providing a solution of the motion planning problem for 396

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398

Fig. 6. Experimental results of the flexible beam structure under disturbances introduced by a pulse hammer at 𝑡 = 1 s. The first row presents the controlled outputs (points of interest) of the open loop 𝒚 pi,ff (𝑡) and the closed loop 𝒚 pi,tr (𝑡) behavior. The second row shows acting input signal of the controller 𝒖tr (𝑡).

a flexible interconnected beam structure. The obtained control law is applied to the in-domain actuation in terms of MFC-actuators. The structure’s deflection is locally measured by strain gauges mounted on the beam surfaces. Hamilton’s principle is used to determine the equations of motion taking into account the different material characteristics of the specific parts of the structure. The flatness-based parameterization of the reduced order states and system inputs provides the desired rest-to-rest motion of the deflection profiles. Model order reduction is addressed by the usage of the second order balanced truncation method. The flat output is constructed by exploiting the spectral system representation, enabling a differential parameterization of the system states and inputs. To increase the performance and robustness to model uncertainties and external disturbances the feedforward control is extended by a state feedback controller with additional integral part using a Kalman filter to estimate the state vector. The resulting two-degrees-of-freedom control strategy is evaluated at an experimental setup and displays a very good performance.

where the second index characterizes the component of the structure, c denotes the carrier structure and a denotes the actuator. The elements are determined according to l

[𝑀𝑗,c ]𝑖,𝑚 = 𝜌𝑗 𝐴𝑗

𝜙𝑗,𝑖 𝜙𝑗,𝑚 d𝑧,

∫0

[𝑀𝑗,a ]𝑖,𝑚 = 𝜌𝑗,a 𝐴𝑗,a

𝑧𝑗,a +la 𝑘

∫𝑧

𝜙𝑗,𝑖 𝜙𝑗,𝑚 d𝑧,

𝑗,a𝑘

[𝑀𝑗,𝑚𝑗 ]𝑖,𝑚 = 𝑚𝑗 𝜙𝑗,𝑖 (l )𝜙𝑗,𝑚 (l ) l

[𝐶𝑗,c ]𝑖,𝑚 = [𝐶𝑗,a ]𝑖,𝑚 =

sv v 𝛾𝑗,c 𝜕𝑧 𝜙𝑗,𝑖 𝜕𝑧 𝜙𝑗,𝑚 + 𝛾𝑗,c 𝜙𝑗,𝑖 𝜙𝑗,𝑚 d𝑧,

∫0

𝑧𝑗,a +la 𝑘

∫𝑧

sv v 𝛾𝑗,a 𝜕𝑧 𝜙𝑗,𝑖 𝜕𝑧 𝜙𝑗,𝑚 + 𝛾𝑗,a 𝜙𝑗,𝑖 𝜙𝑗,𝑚 d𝑧,

𝑗,a𝑘

l

[𝐾𝑗,c ]𝑖,𝑚 = 𝐸𝑗 𝐼𝑗

∫0

𝜕𝑧2 𝜙𝑗,𝑖 𝜕𝑧2 𝜙𝑗,𝑚 d𝑧,

[𝐾𝑗,a ]𝑖,𝑚 = 𝐸𝑗,a𝑘 𝐼𝑗,a𝑘 Appendix. Second order system matrices

𝑧𝑗,a +la 𝑘

∫𝑧

𝜕𝑧2 𝜙𝑗,𝑖 𝜕𝑧2 𝜙𝑗,𝑚 𝑑𝑧,

𝑗,a𝑘

[𝐾sp,pre,𝑗 ]𝑖,𝑚 = 𝑘𝑗−1,𝑗 𝜙𝑗−1,𝑖 (l𝑗−1 )𝜙𝑗,𝑚 (l𝑗 ), Subsequently the structures and elements of the mass, damping, stiffness and input matrix used in (15) are summarized. The structure of the matrices is given by

[𝐾sp,mid,𝑗 ]𝑖,𝑚 = (𝑘𝑗−1,𝑗 + 𝑘𝑗,𝑗+1 )𝜙𝑗,𝑖 (l𝑗 )𝜙𝑗,𝑚 (l𝑗 ), [𝐾sp,post,𝑗 ]𝑖,𝑚 = 𝑘𝑗,𝑗+1 𝜙𝑗+1,𝑖 (l𝑗+1 )𝜙𝑗,𝑚 (l𝑗 ). The spring constant of the first and the last beam is 𝑘0,1 = 𝑘𝑟,𝑟+1 = 0. The elements of the input vectors follow as

𝑀 = diag{𝑀𝑗,c + 𝑀𝑗,a + 𝑀𝑗,𝑚𝑗 }, 𝐶 = diag{𝐶𝑗,c + 𝐶𝑗,a },

[𝒍𝑗,a𝑘 ]𝑖 = −𝛤𝑗,a𝑘

𝐾 = diag{𝐾𝑗,c + 𝐾𝑗,a }+ ⎡ 𝐾sp,mid,1 ⎢ −𝐾sp,pre,2 +⎢ 0 ⎢ 0 ⎢ 0 ⎣ ⎡𝒍1,1 𝐿=⎢ 0 ⎢ ⎣ 0

0 ⋱ 0

−𝐾sp,post,1

0

𝐾sp,mid,2

−𝐾sp,post,2

⋱ ⋱ −𝐾sp,pre,𝑗−1 𝐾sp,mid,𝑗−1 0

−𝐾sp,pre,𝑟

0

⎤ ⎥ ⎥, ⋱ −𝐾sp,post,𝑗−1 ⎥ ⎥ 𝐾sp,mid,𝑟 ⎦

𝑧𝑗,a +la 𝑘

∫𝑧

𝜕𝑧2 𝜙𝑗,𝑖 d𝑧.

𝑗,a𝑘

0

References Abdeljaber, O., Avcia, O., & Inman, D. J. (2016). Active vibration control of flexible cantilever plates using piezoelectric materials and artificial neural networks. Journal of Sound and Vibration, 363, 33–53.

0 ⎤ 0 ⎥, ⎥ 𝒍𝑟,𝑝𝑟 ⎦

Antoulas, A. C. (2005). Approximation of large-scale dynamical systems. In Advances in Design and Control, Vol. 53 (pp. 1689–1699). SIAM. 397

A. Kater and T. Meurer

Control Engineering Practice 84 (2019) 389–398 Meurer, T. (2011). Flatness-based trajectory planning for diffusionreaction systems in a parallelepipedon - a spectral approach. Automatica, 47(5), 935–949. Meurer, T. (2013). Control of higher-dimensional PDEs. Heidelberg: Springer. Meurer, T., & Kugi, A. (2009). Trajectory planning for boundary controlled parabolic PDEs with varying parameters on higher-dimensional spatial domains. IEEE Transactions on Automatic Control, 54(8), 1854–1868. Mohd Jani, J., Leary, M., Subic, A., & Gibson, M. A. (2014). A review of shape memory alloy research, applications and opportunities. Materials & Design, 56, 1078–1113. Morris, K. A., & Vest, A. (2016). Design of damping for optimal energy dissipation of vibrations. In Proc. 55th IEEE Conference on Decision and Control (pp. 532–536). Nowacki, W. (1975). Dynamic problems of thermoelasticity. Warszawa: Noordhoff Int. Publ., PWN–Polish Scientific Publ. Preumont, A. (2011). Vibration control of active structures, an introduction, (3rd ed.). Berlin: Springer. Qiu, Z., Zhang, X., & Ye, C. (2012). Vibration suppression of a flexible piezoelectric beam using bp neural network controller. 25(4), 417–428. Reddy, J. (1984). Energy and variational methods in applied mechanics. New York: Wiley– Interscience. Reis, T., & Stykel, T. (2008). Balanced truncation model reduction of second-order systems. Mathematical and Computer Modelling of Dynamical Systems, 14(5), 391–406. Rudolph, J. (2003). Flatness based control of distributed parameter systems. In Berichte aus der Steuerungs– und Regelungstechnik, Aachen: Shaker–Verlag. Schörkhuber, B., Meurer, T., & Jüngel, A. (2013). Flatness of semilinear parabolic PDEs A generalized Cauchy-Kowalevski approach. IEEE Transactions on Automatic Control, 58(9), 2277–2291. Schröck, J., Meurer, T., & Kugi, A. (2010a). Control of a flexible beam actuated by macro-fiber composite patches: II. Modeling and feedforward trajectory control. Smart Materials and Structures, 20(1), 015015. Schröck, J., Meurer, T., & Kugi, A. (2010b). Control of a flexible beam actuated by macro-fiber composite patches: I. Modeling and feedforward trajectory control. Smart Materials and Structures, 20(1), 015015. Schröck, J., Meurer, T., & Kugi, A. (2013). Motion planning for piezo-actuated flexible structures: Modeling, design, and experiment. IEEE Transactions on Control Systems Technology, 21(3), 807–819. Smart Materials GmbH (2014). Datasheet. Dresden, Germany. http://www.smartmaterial.com. Stanewsky, E. (2001). Adaptive wing and flow control technology. Progress in Aerospace Sciences, 37(7), 583–667. Wloka, J. (1987). Partial differential equations. Cambridge University Press. Yousefi-Koma, A., & Zimcik, D. G. (2003). Applications of smart structures to aircraft for performance enhancement. Canadian Aeronautics and Space Journal, 49(4), 163–172.

Åström, K. J., & Wittenmark, B. (2011). Computer-controlled systems: theory and design (3rd ed.). Dover books on electrical engineering, (p. 555). Dover Publications. Banks, H. T., Demetriou, M. A., & Smith, R. C. (1996). An 𝐻 ∞ MinMax periodic control in a two-dimensional structural acoustic model with piezoceramic actuators. IEEE Transactions on Automatic Control, 41(7), 943–959. Banks, H., Smith, R., & Wang, Y. (1996). Smart material structures: modeling, estimation and control. Chichester: John Wiley & Sons. Barbarino, S., Bilgen, O., Ajaj, R. M., Friswell, M. I., & Inman, D. J. (2011). A review of morphing aircraft. Journal of Intelligent Material Systems and Structures, 22(9), 823– 877. Bryson, A. E. (2002). Applied linear optimal control: examples and algorithms. In Applied mechanics reviews. Cambridge: Cambridge University Press. Chahlaoui, Y., Lemonnier, D., Vandendorpe, A., & Van Dooren, P. (2006). Second-order balanced truncation. Linear Algebra and its Applications, 415(2–3), 373–384. Clark, R. L., Saunders, W. R., & Gibbs, G. P. (1998). Adaptive structures: dynamics and control (p. 488). Wiley. De Abreu, G. L. C. M., & Ribeiro, J. F. (2002). A self-organizing fuzzy logic controller for the active control of flexible structures using piezoelectric actuators. Applied Soft Computing, 1(4), 271–283. dSPACE GmbH - Product Information (2017). MicroLabBox. pp. 1–12. www.dspace.com. Fliess, M., Lévine, J., Martin, P., & Rouchon, P. (1995). Flatness and defect of non–linear systems: introductory theory and examples. International Journal of Control, 61, 1327– 1361. Gelb, A. (1974). Applied optimal estimation. In Applied optimal estimation, MIT Press. Janocha, H. (2007). Adaptronics and smart structures (p. 554). Berlin, Heidelberg: SpringerVerlag. Kater, A., & Meurer, T. (2015). Flachheitsbasierte Bewegungsplanung für gekoppelte elastische Balken. At-Automatisierungstechnik, 63(9), 684–699. Knüppel, T., & Woittennek, F. (2015). Control design for quasi-linear hyperbolic systems with an application to the heavy rope. IEEE Transactions on Automatic Control, 60(1), 5–18. Kugi, A. (2001). Non-linear control based on physical models electrical, mechanical and hydraulic systems. In Lecture notes in control and information sciences, London: SpringerVerlag. Kugi, A., & Thull, D. (2005). Infinite-dimensional decoupling control of the tip position and the tip angle of a composite piezoelectric beam with tip mass. In Meurer, T. and Graichen, K. and Gilles, E.D. (editors), Control and observer design for nonlinear finite and infinite dimensional systems, Vol. 322. (pp. 351–368). Berlin Heidelberg: Springer. Kugi, A., Thull, D., & Kuhnen, K. (2006). An infinite-dimensional control concept for piezoelectric structures with complex hysteresis. Structural Control and Health Monitoring, 13(6), 1099–1119. Lewis, F. L., Xie, L., & Popa, D. (2008). Optimal and robust estimation: with an introduction to stochastic control theory (2nd ed.). In Automation and Control Engineering, CRC Press.

398