A Multiscale Model of Cantilever Arrays with a Semi-decentralized Approximation of an Optimal Control

In this paper, we present a two-scale model including an optimal active control for a one dimensional cantilever array with regularly spaced actuators and sensors. With the purpose of implementing the control in real time, we propose an approximation that may be realized by an analog distributed electronic circuit. More precisely, our analog processor is made by Periodic Network of Resistances (PNR). The control approximation method is based on two general concepts, namely functions of operators and on the Dunford-Schwartz representation formula. We conducted careful validations of the control approximation method as well as of its effect in the complete control loop.


I. 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 are dealing with infinite length systems, see [1] and [10] for systems governed by partial differential equations, and [3] for discrete systems. In the papers [4] and [5] 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 was giving satisfactory results, it was suffering from some limitations. In [9] and [12], it has been extended so that to cover a larger range of systems and to increase its precision and robustness. Indeed, the new method does not require that each operator of the state equation and of the cost functional be functions of a same operator but they must be only functions of a same operator up to some change of variable operators. Regarding precision, the Taylor series approximating a function of an operator has been replaced by the use of the Dunford-Schwartz representation formula followed by a quadrature rule for the contour integral.
Here we apply our new method to a recently developed and validated two-scale model of cantilever arrays, submitted in the paper [8]. It is rigorously justified thanks to an adaptation of the two-scale approximation method introduced in [6] and detailed in [7]. Its main advantage is that in the same time it requires little computing effort and it is reasonably precise.
This paper presents results from a full implementation of the new semi-decentralized optimal control strategy on the two-scale model of cantilever arrays. The calculations have been carried out using a simple optimal control strategy, namely a Linear Quadratic Regulator (LQR), with purpose vibrations cancellation. We are focused on studying the quality of our approximation method, studying its precision and cost but not on practical applications and their results. As in [5], we also provide a realization of the semi-decentralized control scheme through a Periodic Network of Resistances (PNR), implementing a finite difference scheme for the partial differential operator in the Dunford-Schwartz formula. Finally, we quote that the entire approach can be extended to other linear optimal control problems, i.e. LQG or H ∞ controls as well as to more physical actuating and sensing principles.

II. 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 Figure 1. Assuming that the number of cantilevers is sufficiently large, a homogenized model was derived using a two-scale approximation method. This is reported in the detailed paper [7] devoted to static regime. The corresponding model extended to dynamic regime is introduced in the letter [6]. The modelling papers were written in view of Atomic Force Microscopy application.
After a number of simplifications, the approximate homogenized model expressed in the twoscale referential, which is a rectangle Ω = (0, L B )× (0, L * C ). The parameters L B and L * C represent respectively the base length in the macroscale x−direction and the scaled cantilever length in the microscale y−direction. The base is modelled by the line Γ = {(x, y) | x ∈ (0, L B ) and y = 0}, and the rectangle Ω is filled by an infinite number of cantilevers. We describe the system motion by its Fig. 1. Array of Cantilevers bending displacement only. So, the base is governed by an Euler-Bernoulli beam equation with two kinds of distributed forces, one exerted by the attached cantilevers and the other, denoted by u(t, x, 0), originating from an actuator distribution. 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 Γ The base is assumed to be clamped, so the boundary conditions are at its ends. Each cantilever is oriented in the ydirection, 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 bending displacements and the mass per unit length, the governing equation in (x, y) ∈ Ω is endowed with the boundary conditions representing an end clamped in the base, and a free end. The weak formulation associated to (1)(2)(3)(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 the junction.

III. 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 and u n (t, x) are the polynomial coefficients in the approximation of w and u respectively. We also The boundary conditions are w(t, (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,2,··· ,N ∈ L 2 (Γ) N and for the cost functional The choice of the functional is related to vibration stabilization of the microcantilever array.

IV. CLASSICAL FORMULATION OF THE LQR PROBLEM
Now, we write the above LQR problem in a classical abstract setting, see [2], 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, observation operator, and S = I N ×N the weight operator. 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). We know that (A, B) is stabilizable and that (A, C) is detectable, in the sense that the system is controllable and observable. It follows that for each z 0 , the LQR problem (8) admits a unique solution where K = S −1 B * P, and P is the unique selfadjoint nonnegative solution to the operational Riccati equation This Section is devoted to formulate the approximation method. The mathematical derivation has been introduced in a paper [9]. 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

A. 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 (9), the optimal controller K admits the factorization where and where for all λ ∈ σ, p(λ) is the unique self-adjoint nonnegative matrix solving the algebraic Riccati equation

B. Approximation of the Functions of Λ
We build the approximation in two steps. Firstly, we use a rational approximation k R (Λ) of k(Λ), then it is approximated by another function k R,M which is simple to discretize, and yields an accurate approximation. To do so, we use the Dunford-Schwartz formula, see [13], representing a function of an operator, because it involves only the operator (ζI − Λ) −1 which may be simply and accurately approximated. Since the function k(Λ) is not known, the spectrum σ (Λ) cannot be easily determined, so we approximate k(λ) by a highly accurate rational approximation k R (Λ), then the Dunford-Schwartz formula is applied to k R (Λ) with a path tracing out ellipses including I σ but no poles.
Since the interval I σ is bounded, for each function k ij (λ) have a rational approximation over I σ , we write under a global formulation, where d m , d m are matrices of coefficients and R = R N , R D is the couple comprised of the matrices R N of numerator polynomial degrees and the matrices R D of denominator polynomial degrees. The path C, in the Dunford-Schwartz formula, 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, we introduce the 2Ndimensional vector field Decomposing v ζ into its real part v ζ 1 and its imag- Thus, combining the rational approximation k R and the quadrature formula yields an approximate realization k R,M (Λ) of k (Λ) , 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 M systems like (13) corresponding to the M quadrature nodes ζ(θ l ). The matrices k R (ζ(θ l )) could be computed "offline" once and for all, and stored in memory, so their determination would not penalize a rapid realtime 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.

VI. SPATIAL DISCRETIZATION
The final step in the approximation consists in a spatial discretization and synthesis of Equation (13). 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, 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 (15-16) by a distributed electronic circuit that could be integrated in the mechanical structure. For this purpose, the system is rewritten under the manageable form (17-18). For the sake of simplification, we use the notations

A. Analog computation of Λ h v 1 and Λ h v 2
The analog computation of Λ h v 1 and Λ h v 2 are made by Periodic Network of Resistances (PNR) circuits [11]. 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 (Figure 2), 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 Figure 3 with its two adjacent cells on each side. We call ρ 1 the resistance value between the potentials u 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 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 (k) 1 . The values of the resistances inside a cell depend only on the circuit topology and are easily expressed as a function of ρ 1 or ρ 2 .  Figure 4. The imposed values of the voltages correspond to the clamping boundary condition. 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.

B. Analog computation of equation (17)
The analog computation of Equation (17) can be made by an array of classical non inverting summing amplifiers of Figure 5. Notice that there is no current exchange between these circuits and PNR inputs and outputs, see buffers in Figure 3. Analysis of the circuit of Figure 5 leads to (19). With a proper choice of resistances, Figure 5 solves

C. Analog computation of equation (18)
In a very similar way, the analog computation of Equation 18 can made by an array of classical difference summing amplifiers of Figure 6. Analysis of the circuit of Figure 6 leads to (20). With a proper choice of resistances, Figure 6 solves

VII. NUMERICAL SIMULATION
In this Section, we validate the approximation method, established in Section V, by a numerical simulation. We consider a silicon array comprised of an elastic base clamped of 10 elastic cantilevers, with base dimensions L B × l B × h B = 500μm × 16.7μm × 10μm, and one cantilever dimensions The model parameters of base and cantilever: 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 , can be in the order of 10 −8 . This choice has no effect on real-time computation time.
Numerical integrations have been performed with a standard trapezoidal quadrature rule. The relative , between the exact functions and final approximations are shown in Figure 7, for the number of nodes M varying from 5 to 20. It may be easily tuned without changing spatial complexity associated with the finite difference discretization of Λ −1 .
We present Figure 8 to illustrate the displacement We also present the ratio of the computation time of solving the whole system for varying number of nodes M to the reference computation time of solving the whole system for M = 20, see Figure  9. In this paper, we have presented a semidecentralized approximation of a linear optimal control operator applied to a two-scale model of microcantilever arrays. This model is discretized in y-direction projecting on the transformed basis of Chebyshev polynomials. A semi-decentralized optimal controller is implemented by a set of distributed electronic circuits. From numerical simulations we have shown that the computation time for control operator realization is almost linearly increasing with respect to the number of quadrature nodes. A simulation of the displacement evolution has confirmed that the approached optimal control is effective on this model. Furthermore, this method can be extended to other optimal control theories.