Semi-decentralized approximation of a LQR-based controller for a one-dimensional cantilever array

We apply the method of semi-decentralized approximation, introduced in Lenczner and Yakoubi [2009] and Yakoubi [2010], to the linear quadratic regulation of a one-dimensional array of cantilevers with regularly spaced actuators and sensors. It is based on two mathematical concepts, namely on functions of operators, and on the Cauchy integral formula. We evaluate its performances and the errors of approximation. We also propose its implementation in terms of an analog processor, namely a periodic network of resistors. The presented application is based on a two-scale model representing an array of cantilevers. We shortly explain its genesis before to state it in details, and to show validation results.


INTRODUCTION
In the past decade, a number of papers have been focused on semi-decentralized distributed optimal control for systems with distributed actuators and sensors. Most of them deal with infinite length systems, see Bamieh et al. [2002] and Paganini and Bamieh [1998] for systems governed by partial differential equations, and D'Andrea and Dullerud [2003] for discrete systems. In the papers Kader et al. [2000] and Kader et al. [2003] the authors have introduced an approximation of an optimal control to a finite length beam endowed with a periodic distribution of piezoelectric sensors and actuators. Even if it gives satisfactory results, it suffers from some limitations. In Lenczner and Yakoubi [2009] and Yakoubi [2010], a complete framework has been introduced so that to extend it, to cover a larger range of systems, and to increase its precision and robustness. The new method does not require that all operators involved are functions of a self-adjoint bounded operator Λ. They only need to be functions of Λ up to some change of variables. Regarding precision of our method, the Taylor series approximating a function of an operator has been replaced using the functional calculus followed by a quadrature rule for the contour integral.
In this paper, we apply this theory to approximate an optimal control of cantilever arrays, see Fig. 1.
The calculations have been carried out using a simple optimal control strategy, namely a Linear Quadratic Regulator (LQR), for the purpose of cancelling vibrations. Here, we are focused on the quality of our approximation method, so we study its precision and its cost. As in Kader et al. [2003], we also provide a realization of the semi-decentralized control scheme through a Periodic Network of Resistors (PNR). The latter implements a finite difference scheme for the partial differential operator Λ −1 in the Cauchy integral formula of the functional calculus. Finally, we notice that the entire approach can be extended to other linear optimal control problems, e.g. LQG or H ∞ controls. It will apply to any system including a cantilever array, for instance to parallel Atomic Force Microscopes (AFM) or to storage devices, like the millipede, see Eleftheriou et al. [2002].
The simplified model used in the control loop, was announced in Lenczner [2007], and its derivation is detailed in a submitted paper. It is rigorously justified thanks to an adaptation of the two-scale approximation 1 method introduced in Lenczner [1997], and to further results in Lenczner and Smith [2007]. Its main advantage is that it requires little computing effort and it is reasonably precise for large arrays. In this paper, we report validation results of the model in static and dynamic regimes.
The paper is organized as follows. In Section 2, we describe the array geometry, the two-scale approximation method, and our model, which we reformulate, in Section 3, in a way suitable for our semi-decentralized approximation.
The LQR control problem and the semi-decentralized approximation method are described in Section 4 and 5 respectively. In Section 6, we detail the implementation of the approximate optimal control with analog distributed electronic circuits. Numerical validations are reported in Section 7. Finally, two-scale model validation results are detailed in Appendix A.

A TWO-SCALE MODEL OF CANTILEVER ARRAYS
We consider a one-dimensional cantilever array comprised of an elastic base, and a number of clamped elastic cantilevers with free end, see Fig. 2. Assuming that the number of cantilevers is sufficiently large, a homogenized model was derived using a two-scale approximation method. This principle is exploited in the detailed paper Lenczner and Smith [2007] devoted to static regime. The corresponding model extended to dynamic regime is introduced in the letter Lenczner [2007]. Both papers were written in view of AFM application.

Fig. 2. Array of Cantilevers
The two-scale model derivation steps are illustrated in Fig. 3. First, (a) the two-scale transform (also called the unfolding operator) and the two-scale approximation are successively applied to map a thin plate model in bending from the physical domain to a two-scale domain comprised of a reference cell and the macroscopic domains. Then, (b) the displacement variation in the width direction of cantilevers is neglected. In (c), base displacements in the reference cell are explicitly calculated and eliminated to yield the model in the so-called two-scale domain where the optimal control is implemented. Finally, (d) an inverse two-scale transform technique is applied to map the solutions in two-scale domain back to the physical domain.
The approximate homogenized model is expressed in the minimal two-scale domain which is a rectangle Ω = (0, L B ) × (0, L * C ), see Fig. 3. The parameters L B and L * C represent respectively the base length in the macroscale x 1 −direction and the scaled cantilever length in the microscale y 2 −direction. For the sake of simplicity, in the following we denote x 1 and y 2 by x and y. The base is modeled by the line Γ = {(x, y) | x ∈ (0, L B ) and y = 0}, and the rectangle Ω is filled by an infinite number of 1 The approximation is in the sense of small ratio of a cell size to the whole array size. The bending displacement, the mass per unit length, the bending coefficient of base and of cantilevers, and the scaled cantilever width being denoted by w(t, x, 0), ρ B , R B , R C and ℓ * C , the base governing equation states in Γ (1) The base is assumed to be clamped, so the boundary conditions are w = ∂ x w = 0, (2) at its ends. Each cantilever is oriented in the y-direction, and its motion is governed by the Euler-Bernoulli equation distributed along the y-direction. It is subjected to a control force u(t, x, y) taken as distributed along each whole cantilever. It can be replaced by any other realistic force distribution. Denoting by w(t, x, y) and ρ C the bending displacements and the mass per unit length, the governing equation in (x, y) ∈ Ω is endowed with the boundary conditions ∂ y w = 0 at y = 0, ∂ 2 yy w = ∂ 3 yyy w = 0 at y = L * C , representing an end clamped in the base, and a free end. The weak formulation associated to (1-4) states as for any regular function v, satisfying in particular the conditions: v = ∂ x v = 0 at both end of the base and ∂ y v = 0 at y = 0 at base-cantilever junction.

MODEL REFORMULATION
To simplify the model, but keeping its distributed feature, we discretize in the y-direction projecting on a basis K n (y) = y 0 yT ′ n (y)dy, where T n (y) is the basis of Chebyshev polynomial. We define the approximations of the displacement and of the control by w(t, x, y) ≈ N n=1 w n (t, x)K n (y) and u(t, x, y) ≈ N n=1 u n (t, x)K n (y), where w n (t, x) and u n (t, x) are the polynomial coefficients in the approximation of w and u respectively.
We also choose v ≈ N m=1 v m (t, x)K m (y), so we find that (w n (t, x)) n=1,2,··· ,N are the solutions to a set of equations posed on Γ, The boundary conditions are w(t, 0, 0) = ∂ x w(t, 0, 0) = 0, and w(t, L B , 0) = ∂ x w(t, L B , 0) = 0. In (6), we use the notations for the matrices M, K B , K C and B, The LQR problem is set for control variables (u n ) n=1,··· ,N ∈ L 2 (Γ) N and for the cost functional Notice that this functional is appropriate to vibration suppression.

FORMULATION OF THE LQR PROBLEM
Now, we write the above LQR problem in an abstract setting, see Curtain and Zwart [1995], even if we do not detail the functional framework. We set z T = (w n ∂ t w n ) n=1,2,··· ,N the state variable, u T = (u n ) n=1,2,··· ,N the control variable, Consequently, the LQR problem, consisting in minimizing the functional under the constraint (6), may be written under its usual form as with the minimized cost functional (7). Here, A is the infinitesimal generator of a continuous semigroup on the separable Hilbert space We also know that (A, B) is stabilizable and that (A, C) is detectable, in the sense that there exist G ∈ L (Z, U ) and F ∈ L (Y, Z) such that A − BG and that A − F C are the infinitesimal generators of two uniformly exponentially stable continuous semigroups. It follows that for each z 0 ∈ Z, the LQR problem (8) admits a unique solution u * = −Kz (9) where K = S −1 B * P z, and P ∈ L (Z) is the unique selfadjoint nonnegative solution of the operational Riccati equation

SEMI-DECENTRALIZED APPROXIMATION
This section is devoted to formulate the approximation method. The mathematical derivation has been introduced in a paper Lenczner and Yakoubi [2009] and in a thesis Yakoubi [2010]. We denote by Λ, the mapping: The spectrum σ (Λ) is discrete and made up of real eigenvalues λ k . They are solutions to the eigenvalue problem Λφ k = λ k φ k with ||φ k || L 2 (Γ) = 1. In the sequel, I σ = (σ min , σ max ) refers to an open interval that includes the complete spectrum. For a given real valued function g, continuous on I σ , g(Λ) is the linear self-adjoint operator

Factorization of K by a Matrix of Functions of Λ
In this part, we introduce the factorization of the controller K under the form of a product of a matrix of functions of Λ. To do so, we introduce the where M = −(M −1 (K B λ −1/2 + K C λ 1/2 )) N ×N . Endowing Z, U and Y with the inner products (z, and (y, 2N , we find the subsequent factorization of the controller K in (9) which plays a central role in the approximation. Proposition 1. The controller K admits the factorization , and where for all λ ∈ σ, p(λ) is the unique self-adjoint nonnegative matrix solving the algebraic Riccati equation Proof. The algebraic Riccati equation can be found after replacing A, B, C and S by their decomposition in the Riccatti equation (10).
In the sequel, we require that the algebraic Riccati equation (11) admits a unique solution for all λ ∈ I σ which is checked numerically. Remark 2. In this example, Φ U and Φ Z are some matrices of functions of Λ, and so is K, K = k(Λ).
(12) Thus, the approximation is developed directly on k(Λ), but we emphasize that in more generic situations it is pursued on q(Λ). Remark 3. Introducing the isomorphisms Φ Z , Φ Y , and Φ U allows to consider a broad class of problems where the operators A, B, C and S are not strictly functions of a same operator. In this particular application, the observation operator C is composed with the operator ∂ 2 xx . This is taken into account in Φ Y in a manner in which Φ −1 Y CΦ Z is a function of Λ only. Remark 4. We indicate how the isomorphisms Φ Z , Φ Y , and Φ U have been chosen. The choice of Φ Z comes directly from the expression of the inner product (z, ., N and j = N + 1, .., 2N . The expression of Φ Y follows. Choosing Φ U is straightforward.

Approximation of Functions of Λ
Our approximation method is based on the Cauchy integral formula of the functional calculus, see Yosida [1980] representing a function of an operator. We build the approximation in two steps. Since the function k(Λ) is not known, the spectrum σ (Λ) cannot be easily determined, so firstly, the function is approximated by a highly accurate rational approximation. We notice that k(λ) may be a very singular function, see Fig. 4, so for each component k ij (λ), we introduce a rational approximation componentwise, based on the logarithm of λ, where d m , d ′ m ′ are two coefficient matrices, and R = R N , R D is a couple of matrices of polynomial degrees. Then, we approximate it by another function k R,M (λ) which is simple to discretize, and which yields an accurate approximation. To do so, we use the Cauchy integral formula, because it involves only the resolvent (ζI − Λ) −1 , which may be simply and accurately approximated. We apply it to the rational approximation with a path C tracing out an ellipse including I σ but no poles. It is chosen to be an ellipse parameterized by ζ(θ) = ζ 1 (θ) + iζ 2 (θ), with θ ∈  In the following equations, we state that the matrices k R (ζ) associated to the rational approximation of the couple R N , R D . So, for each z ∈ L 2 (Γ) 2N and ζ ∈ C, Decomposing v ζ into its real part v ζ 1 and its imaginary part v ζ 2 , the couple (v ζ 1 , v ζ 2 ) is solution of the system Thus, combining the rational approximation k R and the quadrature formula yields an approximate realization This formula is central in the method, so it is the center of our attention in the simulations. A fundamental remark is that, a "real-time" realization, k R,M (Λ) z, requires solving  (14) corresponding to the M quadrature nodes ζ(θ l ). The matrices k R (ζ(θ l )) could be computed "off-line" once and for all, and stored in memory, so their determination would not penalize a rapid real-time computation. In total, the ultimate parameter responsible of accuracy in a real-time computation, apart from spatial discretization discussed in next Section, is M the number of quadrature points.

SPATIAL DISCRETIZATION
The final step in the approximation consists in a spatial discretization and synthesis of Equation (14). The interval Γ is meshed with regularly spaced nodes separated by a distance h, we introduce Λ −1 h the finite difference discretization of Λ −1 , associated with the clamping boundary condition. In practice, the discretization length h is chosen small compared to the distance between cantilevers. Then, z h denoting the vector of nodal values of z, for each ζ we introduce (v ζ 1,h , v ζ 2,h ), a discrete approximation of (v ζ 1 , v ζ 2 ), solution of the discrete set of equations, (17) Finally, an approximate optimal control, intended to be implemented in a set of spatially distributed actuators, could be estimated from the nodal values, estimated at mesh nodes in the following. We shall propose a synthesization of (16-17) by a distributed electronic circuit that could be integrated in the mechanical structure. For this purpose, the system is rewritten under the manageable form (18)(19). For the sake of simplification, we use the notations The analog computation of Λ h v 1 and Λ h v 2 are made by Periodic Network of Resistances (PNR) circuits Ratier [2009]. These electronic circuits have been developed to solve a large class of PDEs by analog computation. More exactly, PNR circuits compute the finite difference solution of a PDE. PNR circuits are gathering of cells (Fig. 6), the interior cells are indexed by k = 1, . . . , N − 1, while the boundary cells correspond to k = −1, 0, N and N + 1. We will show that the circuits solve the equations Au 1 = i 1 . If the current sources i 1 are replaced by voltage controlled current sources defined by i 1 = gv 1 (with g is a real number), the voltage outputs of the circuits u 1 solve g(Λ h v 1 ) and so Λ h v 1 . The computation of Λ h v 2 is done in the same way. The interior cell k which computes (Λ h v 1 ) k is represented on Fig. 7 with its two adjacent cells on each side. We call ρ 1 the resistance value between the potentials u . By applying the Kirchhoff Current Law (KCL) at node u (k) 1 , rearranging some terms and dividing by h 4 , the equation of the cell k can be written under the form: 1 If one choose the negative potential ρ 1 = −h 4 ρ 0 and the positive potential ρ 2 = h 4 ρ 0 /4, then the potential at node u (k) 1 is expressed as a function of its neighbor voltages as 1 which is the stencil of the differential operation Λ −1 . Consequently, the whole electronic circuit composed of N − 1 cells computes the finite differences approximation u 1 = Λ h i 1 = g (Λ h v 1 ). The numerical value of ρ 0 only changes the magnitude of the voltages u Remark that the terminals denoted by a cross are not connected, so the resistances which are linked by one side at them can be removed without changing the behavior of the circuits. They are saved to show clearly the real difference between interior cells and boundary cells.

Analog computation of equation (18)
The analog computation of Equation (18) can be made by an array of classical non inverting summing amplifiers of Fig. 9. Notice that there is no current exchange between these circuits and PNR inputs and outputs, see buffers in Fig. 7. Analysis of the circuit of Fig. 9 leads to (20). With Fig. 9. Analog computation of the k-th equation (18). a proper choice of resistances, Fig. 9 solves (18), where 1

Analog computation of equation (19)
In a very similar way, the analog computation of Equation  19 can made by an array of classical difference summing amplifiers of Fig. 10. Fig. 10. Analog computation of the k-th equation (19).
Analysis of the circuit of Fig. 10 leads to (21). With a proper choice of resistances, Fig. 10 solves (19),

NUMERICAL SIMULATION
In this section, we validate the approximation method, established in Section 5, by a numerical simulation. We consider a silicon array of 10 cantilevers, with base dimensions 500µm × 16.7µm × 10µm, and one cantilever dimensions 25µm × 10µm × 1.25µm,. The model parameters of base and cantilever are: the bending coefficient R B = 1.09 × 10 −5 N/m, R C = 2.13 × 10 −4 N/m the mass per unit length ρ B = 0.0233kg/m, ρ C = 0.00291kg/m. In the rational approximation, the numerator polynomial degrees R N and the denominator polynomial degrees R D can be chosen sufficiently large (namely R N = R D = 20) so that the relative errors between the exact solution k and its rational approximation k R , can be in the order of 10 −8 . The large R N and R D has no effect on the real-time computation.
Numerical integrations have been performed with a standard trapezoidal quadrature rule. The relative errors, between the exact functions and final approximations are shown in Fig. 11, for M the number of quadrature nodes varying from 5 to 20. We found that the relation between M and E R,M is almost linear. It may be easily tuned without changing spatial complexity associated with the finite difference discretization of Λ −1 . in the two-scale domain with respect to M and h the spatial mesh size in the finite difference scheme is represented in Table 1. For different spatial mesh size, the error E R,M,h is well controlled with a relatively small number M = 20 of quadrature points. Fig. 12, represents evolution of cantilever displacement w at the center of the fifth cantilever for different M number of quadrature nodes when only one array mode is excited. Notice that the reference curve has been computed with M = 20 quadrature nodes. We observe that the displacement evolution for M = 6 is already close to the reference. Our semi-decentralized approximation method has been applied to a linear quadratic regulation of a cantilever array. The system was represented through a two-scale model which validation has been carefully carried out and presented. We have proposed a possible implementation of the semi-decentralized controller as a set of distributed electronic circuits. The method has been validated, and all sources of errors have been quantified. We arrive to the conclusion that the main limitation comes from the spatial mesh size h which need to be quite small to reach a good resolution. Conversely, the number M of quadrature nodes is not needed to be large. This may be interpreted in terms of analog circuit implementation by saying that a large number of resistors is needed in the circuit, and a relatively small number of global analog computations is required to get accurate results. Further applications are now possible, for instance to more complex systems, as two-dimensional arrays, and to more sophisticated optimal control laws involving Riccati equations or inequalities.