# Semi-Decentralized Approximation of a LQR-Based Controller for a One-Dimensional Cantilever Array

H. Hui, Y. Yakoubi, Michel Lenczner, N. Ratier

► To cite this version:

H. Hui, Y. Yakoubi, Michel Lenczner, N. Ratier. Semi-Decentralized Approximation of a LQR-Based Controller for a One-Dimensional Cantilever Array. 18th IFAC World Congress, Aug 2011, Milano, Italy. pp.1422-1428, 10.3182/20110828-6-IT-1002.03451. [hal-00677113](https://hal.science/hal-00677113v1)

HAL Id: hal-00677113

https://hal.science/hal-00677113v1

Submitted on 28 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 CC BY 4.0 - Attribution - International License
Semi-decentralized approximation of a
LQR-based controller for a one-dimensional
cantilever array

H. Hui*,** Y. Yakoubi ** M. Lenczner ** N. Ratier **

* School of Mechatronic Northwestern Polytechnical University, 127 Youyi Xilu, 710072 Xi'an Shaanxi, China (e-mail: hui.hui@femto-st.fr).

** FEMTO-ST, Département Temps-Fréquences Université de Franche-Comté, 26, Chemin de l'Epitaphe, 25030 Besançon Cedex, France, (e-mail: yyakoubi@ann.jussieu.fr, michel.lenczner@utbm.fr, nicolas.ratier@femto-st.fr)

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

**Keywords:** Cantilevers Arrays; Distributed Control; Riccati Equation; Semi-Decentralized Control; Functional Calculus; Cauchy Integral Formula; Two-scale Model; Distributed Analog Electric Circuit; Network of Resistors.

# 1. 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 opti-
mal control strategy, namely a Linear Quadratic Regulator
(LQR), for the purpose of cancelling vibrations. Here, we

<figure><img src="image_2.png" /><figcaption>Fig. 1. One-dimensional arrays of AFM. Courtesy of André Meister and of Thomas Overstolz, CSEM Neuchatel Switzerland.</figcaption></figure>

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 $\Lambda^{-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_\infty$ 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 an-
nounced 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<fnref n="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.

## 2. 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.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Array of Cantilevers</figcaption></figure>

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 $\Omega = (0, L_B) \times (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 $\Gamma = \{(x, y) | x \in (0, L_B) \text{ and } y = 0\}$, and the rectangle $\Omega$ is filled by an infinite number of

<figure><img src="image_4.png" /><figcaption>Fig. 3. Two-scale transform and inverse two-scale transform in two-scale domain</figcaption></figure>

cantilevers. 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$ the 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 + R^B \partial_{xx}^2 w \partial_{xx}^2 v) |_{\Gamma} dx \\ & + \ell_C^* \int_{\Omega} \rho^C \partial_{tt}^2 w + R^C \partial_{yy}^2 w \partial_{yy}^2 v dy dx \quad (5) \\ & = \int_0^{L_B} (u \cdot v) |_{\Gamma} dx + \ell_C^* \int_{\Omega} u \cdot v dy dx, \end{aligned} $$

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 base-cantilever junction.

## 3. 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$

<footnote n="1">¹ The approximation is in the sense of small ratio of a cell size to the whole array size.</footnote>
$\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,\dots,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\dots x}^4 w_n + K_{m,n}^C w_n = \sum_{n,m=1}^{N} \tilde{B}_{m,n} u_n \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,\dots,N} \in L^2(\Gamma)^N$ and for the cost functional

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

Notice that this functional is appropriate to vibration suppression.

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

$$ A = \begin{pmatrix} 0_{N \times N} & I_{N \times N} \\ -(M^{-1}(K^B \partial_{x\dots 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} I_{N \times N} \partial_{xx}^2 & 0_{N \times N} \\ 0_{N \times N} & 0_{N \times N} \end{pmatrix}$ the observation operator, $S = I_{N \times N}$ the weight operator and the functional $J = \int_0^{+\infty} \|Cz\|_Y^2 + (Su, u)_U dt$. 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) $$

with the minimized cost functional (7). Here, $A$ is the infinitesimal generator of a continuous semigroup on the separable Hilbert space $Z = (H_0^2(\Gamma))^N \times (L^2(\Gamma))^N$ with dense domain $D(A) = (H^4(\Gamma) \cap H_0^2(\Gamma))^N \times (H_0^2(\Gamma))^N$. It is known that the control operator $B \in \mathcal{L}(U, Z)$, the

observation operator $C \in \mathcal{L}(Z, Y)$, and $S \in \mathcal{L}(U, U)$, where $Y = L^2(\Gamma)^{2N}$ and $U = L^2(\Gamma)^N$. We also know that $(A, B)$ is stabilizable and that $(A, C)$ is detectable, in the sense that there exist $G \in \mathcal{L}(Z, U)$ and $F \in \mathcal{L}(Y, Z)$ such that $A - BG$ and that $A - FC$ are the infinitesimal generators of two uniformly exponentially stable continuous semigroups. It follows that for each $z_0 \in Z$, the LQR problem (8) admits a unique solution

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

where $K = S^{-1}B^*Pz$, and $P \in \mathcal{L}(Z)$ is the unique self-adjoint nonnegative solution of the operational Riccati equation

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

for all $z \in D(A)$. The adjoint $A^*$ of the unbounded operator $A$ is defined from $D(A^*) \subset Z$ to $Z$ by the equality $(A^*z, z')_Z = (z, Az')_Z$ for all $z \in D(A^*)$ and $z' \in D(A)$. The adjoint $B^* \in \mathcal{L}(Z, U)$ of the bounded operator $B$ is defined by $(B^*z, u)_U = (z, Bu)_Z$, the adjoint $C^* \in \mathcal{L}(Y, Z)$ being defined similarly.

## 5. 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 $\Lambda$, the mapping: $\Lambda: f \to w$, where $w$ is the unique solution of $\partial_{x\dots 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\phi_k \, dx$.

### 5.1 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} I_{N \times N} \Lambda^{\frac{1}{2}} & 0_{N \times N} \\ 0_{N \times N} & I_{N \times N} \end{pmatrix} \in \mathcal{L}(L^2(\Gamma)^{2N}, Z)$, $\Phi_U = I_{N \times N} \in \mathcal{L}(L^2(\Gamma)^N, U)$ and $\Phi_Y = \begin{pmatrix} I_{N \times N} \partial_{xx}^2 \Lambda^{\frac{1}{2}} & 0_{N \times N} \\ 0_{N \times N} & I_{N \times N} \end{pmatrix} \in \mathcal{L}(L^2(\Gamma)^{2N}, Y)$, 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

$$ \begin{aligned} a(\lambda) &= \begin{pmatrix} 0_{N \times N} & I_{N \times N} \\ \tilde{M} & 0_{N \times N} \end{pmatrix}, & b(\lambda) &= \begin{pmatrix} 0_{N \times N} \\ (M^{-1}\tilde{B})_{N \times N} \end{pmatrix}, \\ c(\lambda) &= \begin{pmatrix} I_{N \times N} & 0_{N \times N} \\ 0_{N \times N} & 0_{N \times N} \end{pmatrix}, & s(\lambda) &= I_{N \times N}, \end{aligned} $$

where $\tilde{M} = -(M^{-1}(K^B \lambda^{-1/2} + K^C \lambda^{1/2}))_{N \times N}$. Endowing $Z, U$ and $Y$ with the inner products $(z, z')_Z = (\Phi_Z^{-1}z, \Phi_Z^{-1}z')_{(L^2(\Gamma))^2}$, $(u, u')_U = (\Phi_U^{-1}u, \Phi_U^{-1}u')_{(L^2(\Gamma))^N}$,
and $(y, y')_Y = (\Phi_Y^{-1}y, \Phi_Y^{-1}y')_{(L^2(\Gamma))^{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

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

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. \quad (11)$$

**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 $\lambda \in I_\sigma$ which is checked numerically.

**Remark 2.** In this example, $\Phi_U$ and $\Phi_Z$ are some matrices of functions of $\Lambda$, and so is $K$,

$$K = k(\Lambda). \quad (12)$$

Thus, the approximation is developed directly on $k(\Lambda)$, but we emphasize that in more generic situations it is pursued on $q(\Lambda)$.

**Remark 3.** Introducing the isomorphisms $\Phi_Z$, $\Phi_Y$, and $\Phi_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 $\partial_{xx}^2$. This is taken into account in $\Phi_Y$ in a manner in which $\Phi_Y^{-1}C\Phi_Z$ is a function of $\Lambda$ only.

**Remark 4.** We indicate how the isomorphisms $\Phi_Z$, $\Phi_Y$, and $\Phi_U$ have been chosen. The choice of $\Phi_Z$ comes directly from the expression of the inner product $(z, z')_Z = (\Phi_Z^{-1}z, \Phi_Z^{-1}z')_{(L^2(\Gamma))^{2N}}$ and from

$$(z_n, z'_n)_{H_0^2(\Gamma)} = \left( (\Delta^2)^{\frac{1}{2}} z_n, (\Delta^2)^{\frac{1}{2}} z'_n \right)_{L^2(\Gamma)}$$

with $n = 1, .., N$. For $\Phi_Y$, we start from $C = \Phi_Y c(\Lambda) \Phi_Z^{-1}$ and from the relation

$$(y, y')_Y = (\Phi_Y^{-1}y, \Phi_Y^{-1}y')_{(L^2(\Gamma))^{2N}}$$

which implies that $\partial_{xx}^2 = (\Phi_Y)_{i,i} c_{i,i}(\Lambda) \Lambda^{-\frac{1}{2}}$ and $0 = (\Phi_Y)_{j,j} c_{j,j}(\Lambda)$ with $i = 1, .., N$ and $j = N + 1, .., 2N$. The expression of $\Phi_Y$ follows. Choosing $\Phi_U$ is straightforward.

## 5.2 Approximation of Functions of $\Lambda$

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(\Lambda)$ is not known, the spectrum $\sigma(\Lambda)$ cannot be easily determined, so firstly, the function is approximated by a highly accurate rational approximation. We notice that $k(\lambda)$ may be a very singular function, see Fig. 4, so for each component $k_{ij}(\lambda)$, we introduce a rational approximation componentwise, based on the logarithm of $\lambda$,

$$k_R(\lambda) = \frac{\sum_{m=0}^{R^N} d_m (\ln \lambda)^m}{\sum_{m'=0}^{R^D} d_{m'}'(\ln \lambda)^{m'}}, \quad (13)$$

<figure><img src="image_5.png" /><figcaption>Fig. 4. One component of the function $k(\lambda)$</figcaption></figure>

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}(\lambda)$ which is simple to discretize, and which yields an accurate approximation. To do so, we use the Cauchy integral formula,

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

because it involves only the resolvent $(\zeta I - \Lambda)^{-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_\sigma$ but no poles. It 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$, see Fig. 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. The contour in the Cauchy integral formula</figcaption></figure>

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 \mathcal{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{array}{l} \zeta_1 v_1^\zeta - \zeta_2 v_2^\zeta - \Lambda v_1^\zeta = \mathrm{Re}(-i\zeta' k_R(\zeta)) z, \\ \zeta_2 v_1^\zeta + \zeta_1 v_2^\zeta - \Lambda v_2^\zeta = \mathrm{Im}(-i\zeta' k_R(\zeta)) z. \end{array} \right. \quad (14)$$

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 (15)$$

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 (14) corresponding to the M quadrature nodes $\zeta(\theta_l)$. The matrices $k_R(\zeta(\theta_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.

## 6. SPATIAL DISCRETIZATION

The final step in the approximation consists in a spatial discretization and synthesis of Equation (14). 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^\zeta, v_2^\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} = \mathrm{Re}(-i\zeta' k_R(\zeta)) z_h, \quad (16) $$

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

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 (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 $\alpha = \mathrm{Re}(-i\zeta' k_R(\zeta)) z_h$, $\beta = \mathrm{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 (18) $$

$$ 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 (19) $$

### 6.1 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 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, \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 Fig. 7 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} - \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)} = \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} u_1^{(k-2)} - 4u_1^{(k-1)} + 6u_1^{(k)} - 4u_1^{(k+1)} + u_1^{(k+2)} = \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$. 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 Fig. 7 is controlled by the voltage $v_1^{(k)}$ through the equation $i_1^{(k)} = gv_1^{(k)}$. The four boundary cells are represented on Fig. 8. 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.

### 6.2 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

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

a proper choice of resistances, Fig. 9 solves (18),

$$ v_1^{(k)} = \frac{R_1 + R_2}{R_1} \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, \quad (20) $$

where $\frac{1}{R_u} = \frac{1}{R_a} + \frac{1}{R_b} + \frac{1}{R_c} + \frac{1}{R_d}$.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Analog computation of $\Lambda_h v_1$.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>Fig. 8. Four boundary cells.</figcaption></figure>

## 6.3 Analog computation of equation (19)

In a very similar way, the analog computation of Equation 19 can be made by an array of classical difference summing amplifiers of Fig. 10.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Analog computation of the k-th equation (19).</figcaption></figure>

Analysis of the circuit of Fig. 10 leads to (21). With a proper choice of resistances, Fig. 10 solves (19),

$$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 (21)$$

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

# 7. 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\mu m \times 16.7\mu m \times 10\mu m$, and one cantilever dimensions $25\mu m \times 10\mu m \times 1.25\mu m$. The model parameters of base and cantilever are: 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.0233 kg/m$, $\rho^C = 0.00291 kg/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

$$E_R = \frac{\|k_R - k\|_{L^2(I_\sigma)}}{\|k\|_{L^2(I_\sigma)}},$$

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,

$$E_{R,M} = \frac{\|k_{R,M} - k\|_{L^2(I_\sigma)}}{\|k\|_{L^2(I_\sigma)}},$$

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 $\Lambda^{-1}$.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Relative error between the exact solution $k$ and $k_{R,M}$</figcaption></figure>

The approximation error,

$$E_{R,M,h} = \frac{\|k_{N,M,h}z_h - k(\Lambda)z\|_{\ell^2(I_\sigma)}}{\|k(\Lambda)z\|_{\ell^2(I_\sigma)}},$$

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.

<table><caption>Table 1. Relative errors with respect to $M$ and $h$</caption><thead><tr><th>M \ h</th><th>L<sub>B</sub>/10</th><th>L<sub>B</sub>/20</th><th>L<sub>B</sub>/30</th><th>L<sub>B</sub>/40</th><th>L<sub>B</sub>/50</th></tr></thead><tbody><tr><td>10</td><td>1.45e-1</td><td>1.17e-1</td><td>1.08e-1</td><td>1.06e-1</td><td>1.05e-1</td></tr><tr><td>20</td><td>8.38e-2</td><td>4.08e-2</td><td>2.08e-2</td><td>1.30e-2</td><td>9.09e-3</td></tr><tr><td>50</td><td>8.47e-2</td><td>4.06e-2</td><td>2.10e-2</td><td>1.28e-2</td><td>8.57e-3</td></tr></tbody></table>

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.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Displacement evolution of the first cantilever mode</figcaption></figure>

## 8. CONCLUSION

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.

## ACKNOWLEDGEMENTS

This work is partially supported by the European Territorial Cooperation Programme INTERREG IV A France-Switzerland 2007-2013.

## REFERENCES

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

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

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

E. Eleftheriou, T. Antonakopoulos, G. K. Binnig, G. Cherubini, M. Despont, A. Dholakia, U. Durig, M. A. Lantz, H. Pozidis, H. E. Rothuizen, and P. Vettiger. Milipede - mems-based scanning-probe data-storage system. *IEEE Transactions on Magnetics*, 39(2):938–945, 2002.

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.

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

M. Lenczner. Homogeneisation d'un circuit electrique. *C. R. Acad. Sci. Paris*, Serie II b, t. 324(9):537–542, 1997.

M. Lenczner. A multiscale model for atomic force microscope array mechanical behavior. *Applied Physics Letters*, 90:091908, 2007.

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.

M. Lenczner and Y. Yakoubi. Semi-decentralized approximation of optimal control for partial differential equations in bounded domains. *Comptes Rendus Mécanique*, 337:245–250, 2009.

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

N. Ratier. Towards 2d electronic circuits in the spatial domain. *Proceedings of the 13th WSEAS international conference on Circuits*, pages 212–218, 2009.

Y. 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é, 2010.

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