# A multiscale model of cantilever arrays with a semi-decentralized approximation of an optimal control

Hui Hui, Youssef Yakoubi, Michel Lenczner, Nicolas Ratier

► To cite this version:

Hui Hui, Youssef Yakoubi, Michel Lenczner, Nicolas Ratier. A multiscale model of cantilever arrays with a semi-decentralized approximation of an optimal control. Hardware and Software Implementation and Control of Distributed MEMS (DMEMS), Jan 2010, Besançon, France. pp.56-63, 10.1109/dMEMS.2010.16. hal-00957605

HAL Id: hal-00957605

https://hal.science/hal-00957605

Submitted on 14 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multiscale model of cantilever arrays
with a semi-decentralized approximation of
an optimal control

H. Hui*‡, Y. Yakoubi†, M. Lenczner‡ and N. Ratier ‡

*School of Mechatronic Northwestern Polytechnical University,
127, Youyi Xilu, 710072 Xi’an Shaanxi, China

Email: hui.hui@femto-st.fr

†Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie,
Boîte courrier 187, 75252 Paris Cedex 05, France

Email: yyakoubi@ann.jussieu.fr

‡FEMTO-ST, Département Temps-Fréquences Université de Franche-Comté,
26, Chemin de l'Epitaphe, 25030 Besançon Cedex, France

Email: michel.lenczner@utbm.fr, nicolas.ratier@femto-st.fr

*Abstract*—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.

**Keywords**-Two-scale Model, Cantilevers Arrays, Distributed Control, Semi-Decentralized Approximation, Distributed Analog Electric Circuit

## 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_\infty$ 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 two-scale referential, which is a rectangle $\Omega = (0, L_B) \times (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 $\Gamma = \{(x, y) | x \in (0, L_B) \text{ and } y = 0\}$, and the rectangle $\Omega$ is filled by an infinite number of cantilevers. We describe the system motion by its

<figure><img src="image_2.png" /><figcaption>Fig. 1. Array of Cantilevers</figcaption></figure>

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)$, $\rho^B$, $R^B$, $R^C$ and $\ell_C^*$, the base governing equation states in $\Gamma$

$$ \rho^B \partial_{tt}^2 w + R^B \partial_{x\dots x}^4 w + \ell_C^* R^C \partial_{yyy}^3 w = u. \quad (1) $$

The base is assumed to be clamped, so the boundary conditions are

$$ w = \partial_x w = 0, \quad (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 $\rho^C$ bending displacements and the mass per unit length, the governing equation in $(x, y) \in \Omega$ is

$$ \rho^C \partial_{tt}^2 w + R^C \partial_{y\dots y}^4 w = u, \quad (3) $$

endowed with the boundary conditions

$$ \left\{ \begin{array}{ll} \partial_y w = 0 & \text{at } y = 0, \\ \partial_{yy}^2 w = \partial_{yyy}^3 w = 0 & \text{at } y = L_C^*, \end{array} \right. \quad (4) $$

representing an end clamped in the base, and a free end. The weak formulation associated to (1-4) states as

$$ \begin{aligned} & \int_0^{L_B} (\rho^B \partial_{tt}^2 w \, v + R^B \partial_{xx}^2 w \partial_{xx}^2 v) |_{\Gamma} \, dx \\ & + \ell_C^* \int_{\Omega} \rho^C \partial_{tt}^2 w \, v + R^C \partial_{yy}^2 w \partial_{yy}^2 v \, dy dx \\ = & \int_0^{L_B} (u \, v) |_{\Gamma} \, dx + \ell_C^* \int_{\Omega} u \, v \, dy dx, \end{aligned} \quad (5) $$

for any regular function $v$, satisfying in particular the conditions: $v = \partial_x v = 0$ at both end of the base and $\partial_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) = \int_0^y y T'_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) \approx \sum_{n=1}^{N} w_n(t, x)K_n(y)$ and $u(t, x, y) \approx \sum_{n=1}^{N} 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 \approx \sum_{m=1}^{N} 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 $\Gamma$,

$$ \sum_{n,m=1}^{N} M_{m,n} \partial_{tt}^2 w_n + K_{m,n}^B \partial_{x...x}^4 w_n + K_{m,n}^C w_n = \sum_{n,m=1}^{N} \tilde{B}_{m,n} u_n \quad \text{in } [0, \infty) \times \Gamma. \quad (6) $$

The boundary conditions are $w(t, 0, 0) = \partial_x w(t, 0, 0) = 0$, and $w(t, L_B, 0) = \partial_x w(t, L_B, 0) = 0$. In (6), we use the notations for the matrices $M, K^B, K^C$ and $\tilde{B}$,

$$ \begin{aligned} M_{m,n} &= \rho^B(K_m K_n)|_{\Gamma} + \ell_C^* \rho^C \int_0^{L_C^*} K_m K_n \, dy, \\ K_{m,n}^B &= R^B(K_m K_n)|_{\Gamma}, \\ K_{m,n}^C &= \ell_C^* R^C \int_0^{L_C^*} \partial_{yy}^2 K_m \partial_{yy}^2 K_n \, dy, \\ \tilde{B}_{m,n} &= (K_m K_n)|_{\Gamma} + \ell_C^* \int_0^{L_C^*} K_m K_n \, dy. \end{aligned} $$

The LQR problem is set for control variables $(u_n)_{n=1,2,...,N} \in L^2(\Gamma)^N$ and for the cost functional

$$ \mathcal{J} = \int_{0}^{+\infty} \sum_{n=1}^{N} \left\| \partial_{xx}^2 w_n(t, x) \right\|_{L^2(\Gamma)}^2 \quad (7) \\ + \left\| u_n(t, x) \right\|_{L^2(\Gamma)}^2 \, dt. $$

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 \ \partial_t w_n)_{n=1,2,...,N}$ the state variable, $u^T = (u_n)_{n=1,2,...,N}$ the control variable,

$$ A = \begin{pmatrix} 0_{N \times N} & I_{N \times N} \\ -(M^{-1}(K^B \partial_{x...x}^4 + K^C))_{N \times N} & 0_{N \times N} \end{pmatrix} $$

the state operator, $B = \begin{pmatrix} 0_{N \times N} \\ (M^{-1}\tilde{B})_{N \times N} \end{pmatrix}$ the control operator, $C = \begin{pmatrix} \partial_{xx}^2 I_{N \times N} & 0_{N \times N} \\ 0_{N \times N} & 0_{N \times N} \end{pmatrix}$ the observation operator, and $S = I_{N \times 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

$$ \partial_t z(t, x) = Az(t) + Bu(t) \quad (8) \\ \text{for } t > 0 \text{ and } z(0) = z_0, $$

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

$$ u^* = -Kz, \quad (9) $$

where $K = S^{-1}B^*P$, and $P$ is the unique self-adjoint nonnegative solution to the operational Riccati equation

$$ A^*P + PA - PBS^{-1}B^*P + C^*C = 0. \quad (10) $$

## V. SEMI-DECENTRALIZED APPROXIMATION

This Section is devoted to formulate the approximation method. The mathematical derivation has been introduced in a paper [9]. We denote by $\Lambda$, the mapping: $\Lambda: f \to w$, where $w$ is the unique solution of $\partial_{x...x}^4 w = f$ in $\Gamma$ with the boundary conditions $w = \partial_x w = 0$ for $x = \{0, L_B\}$. The spectrum $\sigma(\Lambda)$ is discrete and made up of real eigenvalues $\lambda_k$. They are solutions to the eigenvalue problem $\Lambda\phi_k = \lambda_k\phi_k$ with $\|\phi_k\|_{L^2(\Gamma)} = 1$. In the sequel, $I_\sigma = (\sigma_{\min}, \sigma_{\max})$ refers to an open interval that includes the complete spectrum. For a given real valued function $g$, continuous on $I_\sigma$, $g(\Lambda)$ is the linear self-adjoint operator on space $L^2(\Gamma)$ defined by $g(\Lambda)z = \sum_{k=1}^{\infty} g(\lambda_k)z_k\phi_k$, where $z_k = \int_{\Gamma} z_k\phi_k \, dx$.

### A. Factorization of $K$ by a Matrix of Functions of $\Lambda$

In this part, we introduce the factorization of the controller $K$ under the form of a product of a matrix of functions of $\Lambda$. To do so, we introduce
the change of variable operators $\Phi_Z = \begin{pmatrix} \Lambda^{\frac{1}{2}} & 0 \\ 0 & I \end{pmatrix}$,

$\Phi_U = I$ and $\Phi_Y = \begin{pmatrix} \partial_x^2 \Lambda^{\frac{1}{2}} & 0 \\ 0 & I \end{pmatrix}$, from which we introduce the matrices of functions of $\Lambda$, $a(\Lambda) = \Phi_Z^{-1} A \Phi_Z$, $b(\Lambda) = \Phi_Z^{-1} B \Phi_U$, $c(\Lambda) = \Phi_Y^{-1} C \Phi_Z$ and $s(\Lambda) = \Phi_U^{-1} S \Phi_U$, simple to implement on a semi-decentralized architecture. A straightforward calculation yield

$$a(\lambda) = \begin{pmatrix} 0 & I \\ \tilde{M} & 0 \end{pmatrix}, b(\lambda) = \begin{pmatrix} 0 \\ M^{-1} \tilde{B} \end{pmatrix},$$

$$c(\lambda) = \begin{pmatrix} I & 0 \\ 0 & 0 \end{pmatrix}, \text{ and } s(\lambda) = I,$$

where $\tilde{M} = -M^{-1}(K^B \lambda^{-1/2} + K^C \lambda^{1/2})$. From (9), the optimal controller $K$ admits the factorization

$$K = k(\Lambda) = \Phi_U q(\Lambda) \Phi_Z^{-1}, \quad (11)$$

where $q(\lambda) = s^{-1}(\lambda)b^T(\lambda)p(\lambda)$, and where for all $\lambda \in \sigma$, $p(\lambda)$ is the unique self-adjoint nonnegative matrix solving the algebraic Riccati equation

$$a^T(\lambda)p + pa(\lambda) - pb(\lambda)s^{-1}(\lambda)b^T(\lambda)p + c^T(\lambda)c(\lambda) = 0.$$

## B. Approximation of the Functions of $\Lambda$

We build the approximation in two steps. Firstly, we use a rational approximation $k_R(\Lambda)$ of $k(\Lambda)$, 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 $(\zeta I - \Lambda)^{-1}$ which may be simply and accurately approximated. Since the function $k(\Lambda)$ is not known, the spectrum $\sigma(\Lambda)$ cannot be easily determined, so we approximate $k(\lambda)$ by a highly accurate rational approximation $k_R(\Lambda)$, then the Dunford-Schwartz formula is applied to $k_R(\Lambda)$ with a path tracing out ellipses including $I_\sigma$ but no poles. Since the interval $I_\sigma$ is bounded, for each function $k_{ij}(\lambda)$ have a rational approximation over $I_\sigma$, we write under a global formulation,

$$k_R(\lambda) = \frac{\sum_{m=0}^{R^N} d_m \lambda^m}{\sum_{m'=0}^{R^D} d'_{m'} \lambda^{m'}}, \quad (12)$$

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,

$$k_R(\Lambda) = \frac{1}{2i\pi} \int_C k_R(\zeta) (\zeta I - \Lambda)^{-1} d\zeta,$$

is chosen to be an ellipse parameterized by $\zeta(\theta) = \zeta_1(\theta) + i\zeta_2(\theta)$, with $\theta \in [0, 2\pi]$. The parametrization is used as a change of variable, so the integral can be approximated by a quadrature formula involving $M$ nodes $(\theta_l)_{l=1,...,M} \in [0, 2\pi]$, and $M$ weights $(\omega_l)_{l=1,...,M}$, $I_M(g) = \sum_{l=1}^{M} g(\theta_l) \omega_l$.

In the following equations, we state that the matrices $k_R(\zeta)$ associated to the rational approximation of the couple $(R^N, R^D)$. So, for each $z \in L^2(\Gamma)^{2N}$ and $\zeta \in C$, we introduce the $2N$-dimensional vector field

$$v^{\zeta} = -i\zeta' k_R(\zeta) (\zeta I - \Lambda)^{-1} z.$$

Decomposing $v^\zeta$ into its real part $v_1^\zeta$ and its imaginary part $v_2^\zeta$, the couple $(v_1^\zeta, v_2^\zeta)$ is solution of the system

$$\left\{ \begin{aligned} \zeta_1 v_1^\zeta - \zeta_2 v_2^\zeta - \Lambda v_1^\zeta \\ = \operatorname{Re}(-i\zeta' k_R(\zeta)) z, \\ \zeta_2 v_1^\zeta + \zeta_1 v_2^\zeta - \Lambda v_2^\zeta \\ = \operatorname{Im}(-i\zeta' k_R(\zeta)) z. \end{aligned} \right. \quad (13)$$

Thus, combining the rational approximation $k_R$ and the quadrature formula yields an approximate realization $k_{R,M}(\Lambda)$ of $k(\Lambda)$,

$$k_{R,M}(\Lambda) z = \frac{1}{2\pi} \sum_{l=1}^{M} v_1^{\zeta(\theta_l)} \omega_l. \quad (14)$$

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}(\Lambda)z$, requires solving $M$ systems like (13) corresponding to the $M$ quadrature nodes $\zeta(\theta_l)$. The matrices $k_R(\zeta(\theta_l))$ could be computed "offline" 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.

## VI. SPATIAL DISCRETIZATION

The final step in the approximation consists in a spatial discretization and synthesis of Equation (13). The interval $\Gamma$ is meshed with regularly spaced nodes separated by a distance $h$, we introduce $\Lambda_h^{-1}$ the finite difference discretization of $\Lambda^{-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 $\zeta$ we introduce $(v_{1,h}^{\zeta}, v_{2,h}^{\zeta})$, a discrete approximation of $(v_1^{\zeta}, v_2^{\zeta})$, solution of the discrete set of equations,

$$ \zeta_1 v_{1,h}^{\zeta} - \zeta_2 v_{2,h}^{\zeta} - \Lambda_h v_{1,h}^{\zeta} = \operatorname{Re} (-i\zeta' k_R(\zeta)) z_h, \quad (15) $$

$$ \zeta_2 v_{1,h}^{\zeta} + \zeta_1 v_{2,h}^{\zeta} - \Lambda_h v_{2,h}^{\zeta} = \operatorname{Im} (-i\zeta' k_R(\zeta)) z_h. \quad (16) $$

Finally, an approximate optimal control, intended to be implemented in a set of spatially distributed actuators, could be estimated from the nodal values,

$$ k_{R,M,h} z_h = \frac{1}{2\pi} \sum_{l=1}^{M} v_{1,h}^{\zeta_l} \omega_l, $$

estimated at mesh nodes in the following. We shall propose a synthesis 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 $\alpha = \operatorname{Re}(-i\zeta' k_R(\zeta)) z_h$, $\beta = \operatorname{Im}(-i\zeta' k_R(\zeta)) z_h$, $v_1 = v_{1,h}^{\zeta}$, and $v_2 = v_{2,h}^{\zeta}$.

$$ v_1 = \frac{\zeta_1}{\zeta_1^2 + \zeta_2^2} (\alpha + \Lambda_h v_1) + \frac{\zeta_2}{\zeta_1^2 + \zeta_2^2} (\beta + \Lambda_h v_2), \quad (17) $$

$$ v_2 = \frac{\zeta_1}{\zeta_1^2 + \zeta_2^2} (\beta + \Lambda_h v_2) - \frac{\zeta_2}{\zeta_1^2 + \zeta_2^2} (\alpha + \Lambda_h v_1). \quad (18) $$

### A. Analog computation of $\Lambda_h v_1$ and $\Lambda_h v_2$

The analog computation of $\Lambda_h v_1$ and $\Lambda_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, \dots, 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(\Lambda_h v_1)$ and so $\Lambda_h v_1$. The computation of $\Lambda_h v_2$ is done in the same way. The interior cell $k$ which computes $(\Lambda_h v_1)_k$ is represented on Figure 3 with its two adjacent cells on each side. We call $\rho_1$ the resistance value between the potentials $u_1^{(k)}$ and $u_1^{(k \pm 2)}$, and $\rho_2$ the resistance value between the potentials $u_1^{(k)}$ and $u_1^{(k \pm 1)}$. By applying the Kirchhoff Current Law (KCL) at node $u_1^{(k)}$, rearranging some terms and dividing by $h^4$, the equation of the cell $k$ can be written under the form:

$$ \frac{1}{h^4} \left[ -\frac{1}{\rho_1} u_1^{(k-2)} - \frac{1}{\rho_2} u_1^{(k-1)} + 2u_1^{(k)} \left( \frac{1}{\rho_1} + \frac{1}{\rho_2} \right) - \frac{1}{\rho_2} u_1^{(k+1)} - \frac{1}{\rho_1} u_1^{(k+2)} \right] = \frac{1}{h^4} i_1^{(k)}. $$

If one choose the negative potential $\rho_1 = -h^4\rho_0$ and the positive potential $\rho_2 = h^4\rho_0/4$, then the potential at node $u_1^{(k)}$ is expressed as a function of its neighbor voltages as

$$ \frac{1}{h^4} \left[ u_1^{(k-2)} - 4u_1^{(k-1)} + 6u_1^{(k)} - 4u_1^{(k+1)} + u_1^{(k+2)} \right] = \rho_0 i_1^{(k)}, $$

which is the stencil of the differential operation $\Lambda^{-1}$. Consequently, the whole electronic circuit composed of $N-1$ cells computes the finite differences approximation $u_1 = \Lambda_h i_1 = g(\Lambda_h v_1)$. The numerical value of $\rho_0$ only changes the magnitude of the voltages $u_1^{(k)}$. The values of the resistances inside a cell depend only on the circuit topology and are easily expressed as a function of $\rho_1$ or $\rho_2$.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Analog computation of $\Lambda_h v_1$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Five adjacent interior cells.</figcaption></figure>

Here the resistances of the cells can be taken as $r_1 = r_3 = r_4 = r_6 = \rho_1/4$ and $r_2 = r_5 = \rho_2/2$.

The VCCS (Voltage Controlled Current Source) $i_1^{(k)}$ of Figure 3 is controlled by the voltage $v_1^{(k)}$ through the equation $i_1^{(k)} = gv_1^{(k)}$. The four boundary cells are represented on 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).

<figure><img src="image_5.png" /><figcaption>Fig. 5. Analog computation of the k-th equation (17).</figcaption></figure>

With a proper choice of resistances, Figure 5 solves (17),

$$v_1^{(k)} = \frac{R_1 + R_2}{R_1} \left[ \begin{aligned} & \frac{R_u}{R_a} \alpha + \frac{R_u}{R_b} g(\Lambda_h v_1)_k \\ & + \frac{R_u}{R_c} \beta + \frac{R_u}{R_d} g(\Lambda_h v_2)_k \end{aligned} \right], \quad (19)$$

where $\frac{1}{R_u} = \frac{1}{R_a} + \frac{1}{R_b} + \frac{1}{R_c} + \frac{1}{R_d}$.

## C. Analog computation of equation (18)

In a very similar way, the analog computation of Equation 18 can be made by an array of classical
<figure><img src="image_6.png" /><figcaption>Fig. 4. Four boundary cells.</figcaption></figure>

difference summing amplifiers of Figure 6.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Analog computation of the k-th equation (18).</figcaption></figure>

Analysis of the circuit of Figure 6 leads to (20). With a proper choice of resistances, Figure 6 solves (18),

$$v_2^{(k)} = \frac{R_v}{R_w} \frac{R'_2}{R'_a} \beta + \frac{R_v}{R_w} \frac{R'_2}{R'_b} g(\Lambda_h v_2)_k - \frac{R'_2}{R'_c} \alpha - \frac{R'_2}{R'_d} g(\Lambda_h v_1)_k, \quad (20)$$

where $\frac{1}{R_v} = \frac{1}{R'_a} + \frac{1}{R'_b} + \frac{1}{R'_c}$ and $\frac{1}{R_w} = \frac{1}{R'_c} + \frac{1}{R'_d} + \frac{1}{R'_e}$.

## 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 \times l_B \times h_B = 500\mu m \times 16.7\mu m \times 10\mu m$, and one cantilever dimensions $L_C \times l_C \times h_C = 41.7\mu m \times 12.5\mu m \times 1.25\mu m$. The model parameters of base and cantilever: the bending coefficient $R^B = 1.09 \times 10^{-5} N/m$, $R^C = 2.13 \times 10^{-4} N/m$ the mass per unit length $\rho^B = 0.0233kg/m$, $\rho^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$, $e = \frac{||k_R - k||_{L^2(I_\sigma)}}{||k||_{L^2(I_\sigma)}}$, 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 errors, $E = \frac{||k_{N,M} - k||_{L^2(I_\sigma)}}{||k||_{L^2(I_\sigma)}}$, 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 $\Lambda^{-1}$.

We present Figure 8 to illustrate the displacement

<figure><img src="image_8.png" /><figcaption>Fig. 7. The relative error between the exact solution and the final approximation</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Displacement evolution of first cantilever mode with approximation of optimal control</figcaption></figure>

evolution $w(t, x, y)$ of the first cantilever mode at the coordinate $(x, y) = (L_B/2, L_C/2)$ with different number of nodes $M$. We choose the displacement evolution for $M = 20$ as a reference.

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.

<figure><img src="image_10.png" /><figcaption>Fig. 9. The ratio of computation time</figcaption></figure>

## VIII. CONCLUSION

In this paper, we have presented a semi-decentralized 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.

## REFERENCES

[1] B. Bamieh, F. Paganini, and M.Dahleh. *Distributed control of spatially invariant systems*. *IEEE Transactions on Automatic Control*, 47(7):1091–1107, 2002.

[2] R. F. Curtain and H. Zwart. An introduction to infinite-dimensional linear systems theory. Texts in Applied Mathematics. Springer-Verlag, 1995.

[3] R. D'Andrea and G. E. Dullerud. *Distributed control design for spatially interconnected systems*. *IEEE Trans. Automat. Control*, 48(9):1478–1495, 2003.

[4] M. Kader, M. Lenczner, and Z. Mrcarica. *Approximation of an optimal control law using a distributed electronic circuit: application to vibration control*. 328(7):547–53, 2000.

[5] M. Kader, M. Lenczner, and Z. Mrcarica. *Distributed optimal control of vibrations: a high frequency approximation approach*. 12(3):437–446, 2003.

[6] M. Lenczner. *A Multiscale Model for Atomic Force Microscope Array Mechanical behavior*. *Applied Physics Letters*, 90:091908, 2007.

[7] M. Lenczner and R. C. Smith. *A Two-scale model for Atomic Force Microscopes Arrays in static operating regime*. *Maths. and. Compt. Modelling*, 46:776–805, 2007.

[8] Michel Lenczner, Emmanuel Pillet, Scott Cogan, and Hui Hui. *A Multiscale Model of Cantilever Arrays and its Updating*. 7:125–135, 2009.

[9] Michel Lenczner and Youssef Yakoubi. *Semi-decentralized Approximation of Optimal Control for Partial Differential Equations in Bounded Domains*. *Comptes Rendus Mécanique*, 337:245–250, 2009.

[10] F. Paganini and B. Bamieh. *Decentralization properties of optimal distributed controllers*. 2(9):1877 – 1882, 1998.

[11] Nicolas Ratier. *Towards 2D electronic circuits in the spatial domain*. *Proceedings of the 13th WSEAS international conference on Circuits*, pages 212–218, 2009.

[12] Youssef Yakoubi. Deux Méthodes d'Approximation pour un Contrôle Optimal Semi-Décentralisé pour des Systèmes Distribués. PhD thesis, Université de Franche-Comté, Besançon, FRANCE, 2010.

[13] K. Yosida. Functional analysis. Classics in Mathematics. Springer-Verlag, reprint of the sixth edition edition, 1980.