# Modelling and Control of a Micro-cantilever Array

S. Cogan, H. Hui, Michel Lenczner, E. Pillet, N. Ratier, Y. Yakoubi

► To cite this version:

S. Cogan, H. Hui, Michel Lenczner, E. Pillet, N. Ratier, et al., Modelling and Control of a Micro-cantilever Array. Micro, Nanosystems and Systems on Chips: Modeling, Control and Estimation, Wiley, pp.2-25, 2010, 9781118622520. hal-00676344

HAL Id: hal-00676344

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

Submitted on 20 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
# Modeling and Control of a Micro-cantilever Array

by Scott COGAN, Hui HUI, Michel LENCZNER, Emmanuel PILLET, Nicolas RATIER, and Y. YAKOUBI

Since its invention by G. Binnig [3], the Atomic Force Microscope (AFM) has opened up new possibilities for a number of operations at a nanoscale level, having an impact across various sciences and technologies. Today, the most popular application of it (AFM) is in the material sciences, biology and fundamental physics, see the reviews of F. Giessibl [14], D. Drakova [12] and R. Garcia and R. Perez [13] among others. The AFM is also used for the manipulation of an object or materials at the nanoscale, for example the parallel Lithography of Quantum Devices [17], [5], investigations into mechanical interactions at the molecular level in biology [6], [29], [23], manipulation of nano-objects [10], [24] and data storage [20], [11], [18], [28], [16], [26] to cite only few of them. A number of research laboratories are now developing large AFM Arrays which can achieve the same kind of task in parallel. The most advanced system is the Millipede from IBM [11] for data storage, but again, a number of new architectures are emerging, see [20], [11], [18], [17], [28], [16], [26], [4], [15], [25].

We are currently developing tools for modelling, identification and control of micro-cantilever arrays like those encountered in Atomic Force Microscope Arrays. In this chapter we report results in this direction. The thread of our approach is to provide light computational methods for complex systems. This concern modelling as well as control. Our mechanical structure model is based on a specific multi-scale technique. For control, we start with a general theory of optimal control applied to our simple cantilever array model and we provide an approximation of the control law which may be implemented on a semi-decentralized computing architecture. In particular it could be implemented under the form of a periodically distributed analog electronic circuit. Even if this implementation remains to be completed, we present in advance a general model of such periodically distributed electronic circuits. It will be applied to fast simulations of electronic circuits realizing our control approximation. The general model has been derived with a modified form of the multi-scale technique used for mechanical structures. In a near future, we intend to couple both multi-scale models so that to run light simulations for matrices of electro-mechanical systems. Associated to our light models we also develop a variety of identification tools. (A compléter par Emmanuel ou Scott).

## 1.1 Modelling and Identification of a Cantilever Array

We present a simplified model of mechanical behavior of large cantilever arrays with discoupled rows in the dynamic operating regime. Since the supporting bases are assumed to be elastic, cross-talk effect between cantilevers is taken into account. The mathematical derivation combines a thin plate asymptotic theory and the two-scale approximation theory, devoted to strongly heterogeneous periodic systems. The model is not standard, so we present some of its features. We explain how each eigenmode is decomposed into a products of a base mode with a cantilever mode. We explain the method used for its discretization, and report results of its numerical validation with full three-dimensional Finite Element simulations. A COMPLETER PAR EMMANUEL OU SCOTT

### 1.1.1 Geometry of the Problem

We consider a two-dimensional array of cantilevers. It is comprised of rectangle parallelepiped bases crossing the array in which rectangle parallelepiped cantilever are clamped. Bases are supposed to be connected in the $x_1$-direction only, so that
the system behaves as a set of discoupled rows. Each of them is clamped at its ends. Concerning the other ends, we report two cases, one for free cantilevers and one for cantilevers equipped with a rigid tip, as in Atomic Force Microscopes. The

<figure><img src="image_2.png" /><figcaption>Figure 1.1: FEM and Two-Scale Model Eigenvalues (a) and Absolute Errors (b)</figcaption></figure>

whole array is a periodic repetition of a same cell, in the two directions $x_1$ and $x_2$, see Figure 1.1 (a). We suppose that the number of columns and of rows of the array are sufficiently large, namely larger or equal to 10. Then, we introduce the small parameter $\epsilon^*$ equals to the inverse $1/N$ of the number of cantilevers in a row. We underline the fact that the technique presented in the rest of the paper can be extended to other geometries of cantilever arrays and even to other classes of microsystem arrays.

### 1.1.2 Two-Scale Approximation

Each point of the three-dimensional space, with coordinates $x = (x_1, x_2, x_3)$, is decomposed as $x = x^c + \epsilon y$, where $x^c$ represents the coordinates of the center of the cell to which $x$ belongs, $\epsilon = \begin{pmatrix} \epsilon^* & 0 & 0 \\ 0 & \epsilon^* & 0 \\ 0 & 0 & 1 \end{pmatrix}$, and $y = \epsilon^{-1}(x - x^c)$ is the dilated relative position of $x$ with respect to $x^c$. Points with coordinates $y$ vary in the so-called reference cell, see the two-dimensional view on Figure 1.1 (b), that is obtained through a translation and the $(x_1, x_2)$-dilatation $\epsilon^{-1}$ of any current cell in the array.

We consider the distributed field $u(x)$, of elastic deflections in the array, and we introduce its two-scale transform,

$$ \hat{u}^{\epsilon}(\tilde{x}, y) = u(x^c + \epsilon y), $$

for any $x = x^c + \epsilon y$ and $\tilde{x} = (x_1, x_2)$. By construction, the two-scale transform is constant, with respect to its first variable $\tilde{x}$, over each cell. Since it depends on the ratio $\epsilon^*$, then it may be approximated by the asymptotic field, denoted by $u^0$, obtained for large number of cells (in both $x_1$ and $x_2$-directions) or equivalently when $\epsilon^*$ approaches (mathematically) 0:

$$ \hat{u}^{\epsilon} = u^0 + O(\epsilon^*). $$

The approximation $u^0$ is called the two-scale approximation of $u$. We mention that, as a consequence of the asymptotic process, the partial function $\tilde{x} \mapsto u^0(\tilde{x}, \cdot)$ is continuous instead of being piecewise constant.

Now, we consider that the field of elastic deflections $u$ is a solution of the Love-Kirchhoff thin elastic plate equation in the whole mechanical structure, including bases and cantilevers. Furthermore, we assume that the ratio of cantilever thickness $h_C$ to base thickness $h_B$ is very small, namely

$$ \frac{h_C}{h_B} \approx \epsilon^{4/3}. \qquad (1.1) $$

This assumption is formulated so that the ratio of cantilever stiffness to base stiffness be very small, namely of the order of $\epsilon^{*4}$. The asymptotic analysis when $\epsilon^*$ vanishes shows that $u^0$ does not depend on the cell variable $y$ in bases and so depends only on the spatial variable $\tilde{x}$.

Next, we remark that $u^0(\tilde{x}, y)$ is a two-scale field, and therefore cannot be directly used as an approximation of the field $u(x)$ in the actual array of cantilevers. So, an inverse two-scale transform is to be applied to $u^0$. However, we remark that $\tilde{x} \mapsto u^0(\tilde{x}, y)$ is continuous, and so $u^0$ does not belong to the range of the two-scale transform operator and it has no
preimage. Hence we introduce an approximated inverse of the two-scale transform, $v(\tilde{x}, y) \mapsto \bar{v}(x)$, in the sense that for any sufficiently regular one-scale function $u(x)$ and two-scale function $v(\tilde{x}, y)$,

$$ \bar{u} = u + O(\varepsilon^*) \text{ and } \bar{v} = v + O(\varepsilon^*). $$

It turns out that $\bar{v}(x)$ is a mean over the cell including $x$ with respect to $\tilde{x} = (x_1, x_2)$ when $x$ belongs to a cantilever

$$ \bar{v}(x) = \langle v(., \varepsilon^{-1}(x-x^c)) \rangle_{\tilde{x}}, $$

and with respect to $x_2$ when $x$ belongs to a base

$$ \bar{v}(x) = \langle v(., \varepsilon^{-1}(x-x^c)) \rangle_{x_2}. $$

In total, we retain $\bar{u}^0$ as an approximation of $u$ in the actual physical system. Note that for the model in dynamics, the deflection $u(t, x)$ is a time-space function. In our analysis we do not introduce a two-scale transformation in time, so the time variable $t$ acts as a simple parameter.

### 1.1.3 Model Description

Now, we describe the model satisfied by the two-scale approximation $u^0(t, \tilde{x}, y)$ of $u(t, x)$. Remark that as the deflection in the Kirchhoff-Love model, $u$ is independent of $x_3$, thus $u^0$ is independent of $y_3$. For further simplicity, we neglect torsions effect i.e. the variations of $y_1 \mapsto u^0(t, \tilde{x}, y)$ in cantilevers. Cantilever motion is governed by a classical Euler-Bernoulli beam equation, in the microscopic variable $y_2$,

$$ m^C \partial_{tt} u^0 + r^C \partial_{y_2 \dots y_2}^4 u^0 = f^C $$

with $r^C = \varepsilon^{*4} E^C I^C$, where $m^C$ is a linear mass, $E^C$ the cantilever elastic modulus, $I^C$ the second moment of cantilever section, and $f^C$ a load per unit length in the cantilever. This model represents motion of an infinite number of cantilevers parameterized by all $\tilde{x} = (x_1, x_2)$.

Bases are also governed by an Euler-Bernoulli equation, in the macroscopic variable $x_1$, where part of loads comes from continuous distributions of cantilever shear forces,

$$ m^B \partial_{tt} u^0 + r^B \partial_{x_1 \dots x_1}^4 u^0 = -d^B \partial_{y_2 \dots y_2}^3 u^0 + f^B $$

with $r^B = E^B I^B$, where $m^B$, $E^B$, $I^B$, $d^B$ and $f^B$ are a linear mass, the base elastic modulus, the second moment of section of the base, a cantilever-base coupling coefficient and the load per unit length in the base.

In the model, cantilevers appear as clamped in bases. So at base-cantilever junctions,

$$ u_{\text{cantilever}}^0 = u_{\text{base}}^0 \quad \text{and} \quad (\partial_{y_2} u^0)_{\text{cantilever}} = 0, \qquad (1.2) $$

because $\partial_{y_2} u^0 = 0$ in bases. Other cantilever ends may be free, with equations,

$$ \partial_{y_2 y_2}^2 u^0 = \partial_{y_2 y_2 y_2}^3 u^0 = 0, \qquad (1.3) $$

or may be equipped with a rigid part (usually a tip in Atomic Force Microscopes), so their equation are

$$ \begin{aligned} & J^R \partial_{tt} \begin{pmatrix} u^0 \\ \partial_{y_2} u^0 \end{pmatrix} + \varepsilon r^C \begin{pmatrix} -\partial_{y_2 y_2 y_2}^3 u^0 \\ \partial_{y_2 y_2}^2 u^0 \end{pmatrix} \\ & = \begin{pmatrix} f_3^R \\ F_3^R + F_2^R \end{pmatrix} \end{aligned} $$

at junctions between elastic parts and rigid parts. Here, $J^R$ is a matrix of moments of the rigid part about the junction-plane, $f_3^R$ is a load in the $y_3$ direction, $F_3^R$ is a first moment of loads about the junction-plane, and $F_2^R$ the first moment of loads in the $y_2$ direction about the beam neutral plane. Finally, base ends are assumed to be clamped in a fixed support,

$$ u^0 = \partial_{x_1} u^0 = 0. \qquad (1.4) $$

The loads $f^C$, $f^B$ and $f^R$ in the model are asymptotic loads which are generally not defined from the physical problem. In practical computations, they are replaced by the two-scale transforms $\hat{f}^C$, $\hat{f}^B$ and $\hat{f}^R$. To be complete, we mention that rows of cantilevers are discoupled, this is why $x_2$ plays only the role of a parameter.
### 1.1.4 Structure of Eigenmodes

There is an infinite number of eigenvalues $\lambda^A$ and eigenvectors $\phi^A(x_1, y_2)$ associated to the model. For convenience, we parameterize them by two independent indices $i$ and $j$, both varying in the infinite countable set $\mathbb{N}$. The first index $i$ refers to the infinite set of eigenvalues $\lambda_i^B$ and eigenvectors $\phi_i^B(x_1)$ of the Euler-Bernoulli beam equation associated to a base. The eigenvalues $(\lambda_i^B)_{i \in \mathbb{N}}$ constitutes a sequence of positive number increasing towards infinity. At each such eigenvalue corresponds another eigenvalue problem associated to cantilevers, which has also a countable infinity of solutions denoted by $\lambda_{ij}^C$ and $\phi_{ij}^C(y_2)$. The index $i$ of $\lambda_i^B$ being fixed, the sequence $(\lambda_{ij}^C)_{j \in \mathbb{N}}$ is a positive sequence increasing towards infinity. In the other side, for fixed $j$ and large $\lambda_i^B$, i.e. large $i$, the sequence $(\lambda_{ij}^C, \phi_{ij}^C)_{i \in \mathbb{N}}$ converges to an eigenelement of the clamped-free cantilever model. The eigenvalues $\lambda_{ij}^A$ of the model are proportional to $\lambda_{ij}^C$. Finally, each eigenvector $\phi_{ij}^A(x_1, y_2)$ is the product of a mode in a base by a mode in a cantilever $\phi_i^B(x_1)\phi_{ij}^C(y_2)$.

### 1.1.5 Model Validation

We report observations made on eigenmode computations. We consider a one-dimensional silicon array of $N$ cantilevers ($N = 10, 15$ or $20$), with base dimensions $500\mu m \times 16.7\mu m \times 10\mu m$, and cantilever dimensions $41.7\mu m \times 12.5\mu m \times 1.25\mu m$, see Figure 1.2 for the two possible geometries, with or without tips. We have carried out our numerical study on

<figure><img src="image_3.png" /><figcaption>Figure 1.2: Cantilever Array with tips (a) and without tips (b)</figcaption></figure>

both cases, but we limit the following comparisons to cantilevers without tips, because configuration including tips yields comparable results.

<figure><img src="image_4.png" /><figcaption>Figure 1.3: Eigenmode Density Distributions for Finite Element Model and for the Two-Scale Model</figcaption></figure>

We restrict our attention to a finite number $n^B$ of eigenvalues $\lambda_i^B$ in the base. Computing the eigenvalues $\lambda^A$, we observe that they are grouped in bunches of size $n^B$ accumulated around a clamped-free cantilever eigenvalues. A number of eigenvalues are isolated far from the bunches. It is remarkable that the eigenelements in a same bunch share a same cantilever mode shape, (close to a clamped-free cantilever mode) even if they correspond to different indices $j$. This is why, these modes will be called *cantilever modes*. Isolated eigenelements share also a common cantilever shape, which looks like a first clamped-free cantilever mode shape excepted that the clamped side is shifted far from zero. The induced global mode $\phi^A$ is then dominated by base deformations and therefore will be called *base modes*. Densities of square root of eigenvalues are reported in the sub-figures 2, 4 and 6¹ of Figure 1.3 for $n_B = 10, 15$ and $20$ respectively. These

<footnote n="1">¹Sub-figures are counted from top to down.</footnote>
figures show three bunches with size $n_B$ and isolated modes that remain unchanged.

We discuss the comparison with the modal structure of the three-dimensional linear elasticity system for the cantilever array discretized by a standard finite element procedure. The eigenvalues of the three-dimensional elasticity equations constitute also an increasing positive sequence that accumulate at infinity. As for the two-scale model, its density distribution exhibits a number of concentration points and also some isolated values. Here bunch sizes equal the number $N$ of cantilevers, see sub-figures 1, 3 and 5 in Figure 1.3 representing eigenmode distributions for $N = 10, 15$ and $20$. Extrapolating this observation shows that when the number of cantilevers increases to infinity bunch size increases proportionally. Since the two-scale model is an approximation in the sense of an infinitely large number of cantilevers, this explains why the two-scale model spectrum exhibit mode concentration with infinite number of elements. This remark provides guidelines for operating mode selection in the two-scale model. In order to determine an approximation of the spectrum for an $N$-cantilevers array, we suggest to operate a truncation in the mode list so that to retain a simple infinity of eigenvalues ($\lambda_{ij}^A$) $i=1,...,N$ and $j \in \mathbb{N}$. We stress the fact that $N$-eigenvalue bunches are generally not corresponding to a single column of the truncated matrix $\lambda_{ij}^A$. This comes from the base mode distribution in this list. When considered in increasing

<figure><img src="image_5.png" /><figcaption>Figure 1.4: (a) Superimposed Eigenmode Distributions of the simple Two-Scale Model with the full three-dimensional Finite Element Model (b) Errors in logarithmic scale</figcaption></figure>

order, base modes are located in consecutive lines of the matrix $\lambda^A$ but not necessary in a same column. We remark that a number of eigenvalues in the Finite Element model spectrum have not their counterpart in the two-scale model spectrum. The missing elements correspond to physical effects not taken into account in the Euler-Bernoulli models for bases and cantilevers.

The next step in the discussion is to compare the eigenmodes and especially those belonging to bunches of eigenvalues. To compare an eigenvector from the two-scale model with an eigenvector of the elasticity system, we use the *Modal Assurance Criterion*, see [2] which is equal to one when the shapes are identical and to zero when they are orthogonal, see Figure 1.5. We compare some eigenmodes which have MAC value is near to 1, see Figure1.6.

This test has been applied on transverse displacement only and a further selection has been developed so that to eliminate modes corresponding to physical effects not modeled by the Euler-Bernoulli models. Following this procedure, mode pairing is achieved successfully. In Figure 1.4 (a) paired eigenvalues have been represented and the corresponding relative errors are plotted on Figure 1.4 (b). Note that errors are far from being uniform among eigenvalues. In fact, the main error source resides in a poor precision of the Euler-Bernoulli model for representing base deformations in few particular cases. Indeed, a careful observation of Finite Element modes shows that base torsion is predominant for some modes. This is especially true for the first mode of the first cantilever mode bunch.

### 1.1.6 Model Identification

#### Global Sensitivity Analysis (GSA)

GSA has the objective of studying the effect of the sensitivity of parameters. We denote $m = [m_1 m_2 \cdots m_{n_p}]^T \in \mathfrak{M}$ as the vector of parameters which describe the model and $d = \{d_1, d_2, \cdots, d_{n_d}\} \in \mathfrak{D}$ describe observable data. The exact
<figure><img src="image_6.png" /><figcaption>Figure 1.5: MAC matrix between two-scale model modes and FEM modes</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 1.6: Eigenmode shapes of analytical mode and FEM mode</figcaption></figure>

relation between *m* and *d* is $d = g(m)$. In this model, the parameters are Young's modulus, Poisson ratio, volume mass, the thickness, length and width of base, cantilever and tip. All the parameters are used in GSA. The list of eigenmodes is $(\varphi_{ij}^A)_{i=1,...,10 \text{ and } j=1,2}$. The index *i* and *j* represent "base modes" and "cantilever modes" respectively, see Figure 1.7.

<figure><img src="image_8.png" /><figcaption>Figure 1.7: Eigenmodes of model (B=base mode, C=cantilever mode)</figcaption></figure>

For each parameter, the probability law is *uniform* and the variation between 0.8 and 1.2 times the nominal value. The correlation coefficients matrix is presented in Figure 1.8 (a), which horizontal base are the coefficients and vertical base are eigenmodes. Figure 1.8 (b) and (c) are figures of Singular Value Decomposition(SVD) matrix and Singular value.

<figure><img src="image_9.png" /><figcaption>Figure 1.8: (a)Correlation coefficients matrix, (b)Singular Value Decomposition matrix, (c)Singular values</figcaption></figure>

From Figure 1.8 (a), we can see that the parameters *hB* and *Lbeam* are influential. So, we only consider these two
parameters by following.

## Updating by Sensitivity

Parameter updating through sensitivity is an iterative procedure based on eigenelements sensitivities. The convergence algorithm is governed by the evolution of a cost function which returns the computation of the minimum of difference between experimental data and calculated data. According to previous analysis, we note that parameters *hB* and *Lbeam* are perturbed. We set *hB* to 1.3 and *Lbeam* to 0.8. After 9 iterations, the convergence is reached and the exact value of the reference parameters (all equal to 1) is returned, see Figure 1.9

<figure><img src="image_10.png" /><figcaption>Figure 1.9: Evolutions of (a)cost objective function (b)perturbed parameters</figcaption></figure>

## Inverse Identification

<figure><img src="image_11.png" /><figcaption>Figure 1.10: Evolutions of $\hat{\mu}_{n_s}^t$ for (a)hB and (b)observation 1; evolutions of $\hat{S}_{n_s}^t$ for (c)hB and (d)observation 1</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 1.11: Results of identification for the parameters (a)hB (b)Lbeam</figcaption></figure>

The inverse problem approach is to identify the characteristics parameters *m* in physical system and to find out effected observable data in the model. We adopt the formulations of inverse problems approach which we use in this paper have been developed by Tarantola since twenty years[1]. Here, we have $m = \{m_1, m_2\} = \{hB\_coef, Lbeam\_coef\}$. The observation data are the eigenvalues ($\lambda_{ij}^A$)$_{i=1,...,10 \text{ and } j=1,2}$ and $d^{obs} = \{d_1^{obs},...,d_{20}^{obs}\}$. If we are in the classical Bayesian inference, the marginal probability *a* *posteriori* of the model parameters $\sigma_M^E(m)$ represents the conditional probability of the observations *d* given any *m*. The second case appears when the experiences and the model uncertainties are Gaussian. In our case, it is Gaussian. It is often impossible to calculate *a* *posteriori* directly. So, we estimated the density by a Monte Carlo simulation. As proposed in[22], an algorithm of Metropolis-Hastings[21] is utilized.

With the utilization of MCMC algorithm, the convergence of average $\hat{\mu}_{n_s}^t$ and partial differences $\hat{S}_{n_s}^t$ is reached after 124 iterations. The densities are estimated with 500 samples. The vertical indicate the nominal values of the parameters from
which the observations of $d^{obs}$ have been simulated, see Figure 1.10 and Figure 1.11. The dispersion diagrams a posteriori

<figure><img src="image_13.png" /><figcaption>Figure 1.12: Dispersion diagrams between (a) observation 1 and parameters (c) parameters</figcaption></figure>

between $hB$, $Lbeam$ and observation 1 or between $hB$ and $Lbeam$ are also presented in Figure 1.12

The applications of model identification are still to be improved. The correlation coefficients matrix will be different with the one we presented in this paper when the parameters and the geometry of model have been changed. For instance, the number of cantilevers $N$ is changed 15 or 20. We also interest in the consequent results of cost function and inverse identification.

## 1.2 Semi-Decentralized Approximation of an Optimal Control applied to a Cantilever Array

We apply a recently developed general theory of optimal control approximation to the cantilever array model. The theory applies to the field of finite length distributed systems where actuators and sensors are regularly spaced. It yields approximations implementable on semi-decentralized architectures. Our result is limited to the Linear Quadratic Regulator, but its extension to other optimal control theories for linear distributed systems like LQG or $H_\infty$ controls is in progress. We focus on illustrating the method more than on providing a mathematically rigorous treatment. In the sequel, we begin with transforming the two-scale model of cantilever arrays into an appropriate form. Then, all construction steps of the approximate Linear Quadratic Regulator are fully presented. Finally, we report numerical simulation results.

### 1.2.1 General Notations

The norm and the inner product of an Hilbert space $E$ are denoted by $\|\cdot\|_E$ and $(\cdot, \cdot)_E$. For a second Hilbert spaces $F$, $\mathcal{L}(E, F)$ denotes the space of continuous linear operators defined from $E$ to $F$. In addition, $\mathcal{L}(E, E)$ is denoted by $\mathcal{L}(E)$. One says that $\Phi \in \mathcal{L}(E, F)$ is an isomorphism from $E$ to $F$ if $\Phi$ is a one-to-one continuous mapping with a continuous inverse.

### 1.2.2 Reformulation of the Two-Scale Model of Cantilever Arrays

We reformulate the two-scale model presented in Section 1.1 in a set of notations which is more usual in control theory of infinite dimensional systems. We adopt the configuration of the cantilevers without tip, see Figure 1.13. The model

<figure><img src="image_14.png" /><figcaption>Figure 1.13: Array of Cantilevers</figcaption></figure>

expressed in the two-scale referential appears as posed in 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 direction $x$ and the scaled cantilever length in the microscale variable $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 cantilevers. We recall that the system motion is described by its bending displacement only. The base being governed by an Euler-Bernoulli beam equation, here we consider two kind of distributed forces, one exerted by the attached cantilevers and the other, denoted by $u_B(t,x)$, originating from an actuator distribution. The bending displacement, the mass per unit length, the bending coefficient and the width being denoted by $w_B(t,x)$, $\rho^B$, $R^B$ and $\ell_C$, the base governing equation states

$$ \rho^B \partial_{tt}^2 w_B + R^B \partial_{x...x}^4 w_B = -\ell_C R^C \partial_{yyy}^3 w_C - \partial_{xx}^2 u_B. \quad (1.5) $$

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

$$ w_B = \partial_x w_B = 0 \quad (1.6) $$

at both ends. Cantilevers are oriented in the y-direction, and we recall that their motions are governed by an infinite number of Euler-Bernoulli equations distributed along the x-direction. Here, each cantilever is subjected to a control force $u_C(t,x)$ taken, for simplicity, constant along cantilevers. This choice does not affect the method presented hereafter, so it can be replaced by any other realistic force distribution. Denoting by $w_C(t,x,y)$, $\rho^C$ and $R^C$ cantilever bending displacements, mass per unit length, and bending coefficient, the governing equation in $(x,y) \in \Omega$ is

$$ \rho^C \partial_{tt}^2 w_C + R^C \partial_{y...y}^4 w_C = u_C, $$

endowed with the boundary conditions

$$ \left\{ \begin{array}{ll} w_C = w_B \text{ and } \partial_y w_C = 0 & \text{at } y = 0 \\ \partial_{yy}^2 w_C = \partial_{yyy}^3 w_C = 0 & \text{at } y = L_C \end{array} \right. \quad (1.7) $$

representing an end clamped in the base, and a free end. Finally, both equations are supplemented with initial conditions on displacements and velocities,

$$ \begin{aligned} w_B &= w_{B,0}, \quad \partial_t w_B = w_{B,1}, \\ w_C &= w_{C,0}, \quad \partial_t w_C = w_{C,1}. \end{aligned} $$

The LQR problem is set for control variables $(u_B, u_C) \in U = H^2 \cap H_0^1(\Gamma) \times L^2(\Gamma)$ and for the cost functional

$$ \begin{aligned} \mathcal{J} (w_{B,0}, w_{B,1}, w_{C,0}, w_{C,1}; u_B, u_C) = & \int_0^\infty \|w_B\|_{H_0^2(\Gamma)}^2 + \|\partial_{yy}^2 w_C\|_{L^2(\Omega)}^2 \\ & + \|u_B\|_{H^2 \cap H_0^1(\Gamma)}^2 + \|u_C\|_{L^2(\Gamma)}^2 dt. \end{aligned} \quad (1.8) $$

### 1.2.3 Model Reformulation

The first step, in applying the method, consists in transforming the control problem into a control problem into another problem with internal distributed control and observation. To do so, we are lead to make additional assumptions yielding model simplifications. We set $\bar{w}_C = w_C - w_B$, solution of an Euler-Bernoulli equation in cantilevers with homogeneous boundary conditions

$$ \left\{ \begin{array}{ll} \rho^C \partial_{tt}^2 \bar{w}_C + R^C \partial_{y...y}^4 \bar{w}_C = u_C - \rho^C \partial_{tt}^2 w_B & \text{in } \Omega, \\ \bar{w}_C = \partial_y \bar{w}_C = 0 & \text{at } y = 0, \\ \partial_{yy}^2 \bar{w}_C = \partial_{yyy}^3 \bar{w}_C = 0 & \text{at } y = L_C. \end{array} \right. \quad (1.9) $$

We introduce the basis of normalized eigenfunction $(\psi_k)_k$, solution of the corresponding eigenvalue problem

$$ \left\{ \begin{array}{l} \partial_{y...y}^4 \psi = \lambda^C \psi \text{ in } (0, L_C), \\ \psi(0) = \partial_y \psi(0) = 0, \\ \partial_{yy}^2 \psi(L_C) = \partial_{yyy}^3 \psi(L_C) = 0, \\ \|\psi_k\|_{L^2(0,L_C)} = 1. \end{array} \right. \quad (1.10) $$

It is well known that, in most practical applications, a very small number of cantilever modes is sufficient to properly describe the system. For the sake of simplicity, we take into account only the first one, keeping in mind that the method can handle more than one mode. Therefore, we adopt the approximation

$$ \bar{w}_C(t,x,y) \approx \bar{w}_C^1(t,x) \psi_1(y), $$

where $\bar{w}_C^1$ is the coefficient of the first mode $\psi_1$ in the modal decomposition of $\bar{w}_C$. Introducing the mean $\bar{\psi}_1 = \int_0^{L_C} \psi_1 dy$, and $u_C^1 = \int_0^{L_C} u_C \psi_1 dy$, we find that $\bar{w}_C^1$ is solution of

$$ \rho^C \partial_{tt}^2 \bar{w}_C^1 + R^C \lambda_1^C \bar{w}_C^1 = u_C^1 - \rho^C \bar{\psi}_1 \partial_{tt}^2 w_B. $$
In order to avoid the term $\partial_{tt}^2 w_B$, we introduce $\tilde{w}_C = \tilde{w}_C^1 + \tilde{\psi}_1 w_B$, so as to make $\tilde{w}_C$ be solution of

$$ \rho^C \partial_{tt}^2 \tilde{w}_C + R^C \lambda_1^C \tilde{w}_C - R^C \lambda_1^C \tilde{\psi}_1 w_B = u_C^1. \quad (1.11) $$

Since,

$$ \partial_y^3 w_C = \partial_y^3 (\tilde{w}_C + w_B) = \partial_y^3 [\tilde{w}_C^1 \psi_1 + w_B \psi_1] = \partial_y^3 \psi_1 \tilde{w}_C, $$

we set $c_1 = \partial_y^3 \psi_1(0)$, and obtain that the couple $(w_B, \tilde{w}_C)$ is solution of the system of equations posed on $\Gamma$,

$$ \left\{ \begin{array}{ll} \rho^B \partial_{tt}^2 w_B + R^B \partial_{x_{xx}}^4 w_B + \ell_C R^C c_1 \tilde{w}_C = - \partial_{x_{xx}}^2 u_B & \text{in } \Gamma, \\ \rho^C \partial_{tt}^2 \tilde{w}_C + R^C \lambda_1^C \tilde{w}_C - R^C \lambda_1^C \tilde{\psi}_1 w_B = u_C^1 & \text{in } \Gamma, \end{array} \right. \quad (1.12) $$

with the boundary conditions (1.6). The cost functional is simplified accordingly,

$$ \mathcal{J} \approx \int_0^\infty \|\partial_{xx}^2 w_B (t,x)\|_{L^2(\Gamma)}^2 + \|\lambda_1^C \tilde{w}_C (t,x)\|_{L^2(\Gamma)}^2 \\ + \|\partial_{xx}^2 u_B\|_{L^2(\Gamma)}^2 + \|u_C^1\|_{L^2(\Gamma)}^2 \, dt. \quad (1.13) $$

### 1.2.4 Classical Formulation of the LQR Problem

Now, we write the above LQR problem in a classical abstract setting, see [8]. We set $z^T = (w_B \quad \tilde{w}_C \quad \partial_t w_B \quad \partial_t \tilde{w}_C)$

the state variable, $u^T = (u_B \quad u_C^1)$ the control variable, $A = \begin{pmatrix} 0 & 0 & I & 0 \\ 0 & 0 & 0 & I \\ -R^B \partial_{x_{xx}}^4 / \rho^B & -\ell_C R^C c_1 / \rho^B & 0 & 0 \\ R^C \lambda_1^C \tilde{\psi}_1 / \rho^C & -R^C \lambda_1^C / \rho^C & 0 & 0 \end{pmatrix}$ the state operator,

$B = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ -\partial_{xx}^2 & 0 \\ 0 & I/\rho^C \end{pmatrix}$ the control operator, $C = \begin{pmatrix} I & 0 & 0 & 0 \\ 0 & \lambda_1^C I & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$ the observation operator, $S = I$ the weight operator

and the functional $J(z_0, u) = \int_0^\infty \|Cz\|_Y^2 + (Su, u)_U dt$. Consequently, the LQR problem, consisting in minimizing the functional under the constraint (1.12), may be written under its usual form as

$$ \frac{dz}{dt}(t) = Az(t) + Bu(t) \quad \text{for } t > 0 \text{ and } z(0) = z_0, \quad (1.14) \\ \min_{u \in U} J(z_0, u). $$

Here, $A$ is the infinitesimal generator of a continuous semigroup on the separable Hilbert space $Z = H_0^2(\Gamma) \times L^2(\Gamma)^3$ with dense domain $D(A) = H^4(\Gamma) \cap H_0^2(\Gamma) \times L^2(\Gamma) \times H_0^2(\Gamma) \times L^2(\Gamma)$. 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 = Z$. 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 (1.14) admits a unique solution

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

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

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.

### 1.2.5 Semi-Decentralized Approximation

This section is devoted to formulate, step by step, the method of approximation.

#### Matrices of Functions of a Self-Adjoint Operator

Since the approximation method of $P$ is based on the concept of matrices of functions of a self-adjoint operator, this section is devoted to their definition. We discuss only the simplest case of compact operators which avoid spectral theory
technicalities, because it is enough for the present example, see [9] for the general theory. From now on, we denote by X the separable Hilbert space L²(Γ) and by Λ the self-adjoint operator (∂_x^4..._x)^-1 with domain

$$D(\Lambda) = H^4(\Gamma) \cap H_0^2(\Gamma) \text{ in } X.$$

As Λ is self-adjoint and compact, its spectrum σ(Λ) is discrete, bounded and made up of real eigenvalues λ_k. They are solutions of the eigenvalue problem Λφ_k = λ_kφ_k with ||φ_k||_X = 1. In the sequel, I_σ = (σ_min, σ_max) refers to an open interval that includes σ(Λ).

For a given real valued function f, continuous on I_σ, f(Λ) is the linear self-adjoint operator on X defined by

$$f(\Lambda)z = \sum_{k=1}^{\infty} f(\lambda_k) z_k \phi_k$$

where z_k = (z, φ_k)_X, with domain

$$D(f(\Lambda)) = \{z \in X \mid \sum_{k=1}^{\infty} |f(\lambda_k)z_k|^2 < \infty\}.$$

Then, if f is a n₁ × n₂ matrix of real valued functions fᵢⱼ, continuous on I_σ, f(Λ) is a matrix of linear operators fᵢⱼ(Λ) with domain

$$D(f(\Lambda)) = \{z \in X^{n_2} \mid \sum_{k=1}^{\infty} \sum_{j=1}^{n_2} |f_{ij}(\lambda_k)(z_j)_k|^2 < \infty \quad \forall i = 1..n_1\}.$$

## Factorization by a Matrix of Functions of Λ

The second step in the semi-decentralized control approximation method is the factorization of K under the form of a product of a function of Λ with operators admitting a natural semi-decentralized approximation. To do so, we introduce three isomorphisms Φ_Z ∈ L(X^4, Z), Φ_U ∈ L(X^2, U), and Φ_Y ∈ L(X^4, Y) mapping a power of X into Z, U, and Y respectively, so that

$$a(\Lambda) = \Phi_Z^{-1}A\Phi_Z, \quad b(\Lambda) = \Phi_Z^{-1}B\Phi_U,$$

$$c(\Lambda) = \Phi_Y^{-1}C\Phi_Z, \text{ and } s(\Lambda) = \Phi_U^{-1}S\Phi_U$$

be some matrices of functions of Λ. In the present example, we propose

$$\Phi_Z = \begin{pmatrix} \Lambda_2^1 & 0 & 0 & 0 \\ 0 & I & 0 & 0 \\ 0 & 0 & I & 0 \\ 0 & 0 & 0 & I \end{pmatrix}, \quad \Phi_U = \begin{pmatrix} (-\partial_{xx}^2)^{-1} & 0 \\ 0 & I \end{pmatrix}, \quad \text{and} \quad \Phi_Y = \Phi_Z.$$

This choice yields

$$a(\lambda) = \begin{pmatrix} 0 & 0 & \lambda^{-1/2} & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{R^B}{\rho^B} \lambda^{-1/2} & -\frac{\ell_C R^C c_1}{\rho^B} & 0 & 0 \\ \frac{R^C \lambda_1^C \bar{\psi}_1}{\rho^C} \lambda^{1/2} & -\frac{R^C \lambda_1^C}{\rho^C} & 0 & 0 \end{pmatrix},$$

$$b(\lambda) = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ \frac{1}{\rho^B} & 0 \\ 0 & \frac{1}{\rho^C} \end{pmatrix}, \quad c(\lambda) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \lambda_1^C & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad s(\lambda) = 1.$$

Endowing Z, U and Y with the inner products (z, z')_Z = (Φ_Z^{-1}z, Φ_Z^{-1}z')_{X^4}, (u, u')_U = (Φ_U^{-1}u, Φ_U^{-1}u')_{X^2}, and (y, y')_Y = (Φ_Y^{-1}y, Φ_Y^{-1}y')_{X^4}, we find the subsequent factorization of the controller K in (1.15) 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(λ) = s⁻¹(λ)bᵀ(λ)p(λ), and where for all λ ∈ σ, p(λ) is the unique self-adjoint nonnegative matrix solving the algebraic Riccati equation

$$\begin{aligned} a^T(\lambda)p + pa(\lambda) - pb(\lambda)s^{-1}(\lambda)b^T(\lambda)p \\ + c^T(\lambda)c(\lambda) = 0. \end{aligned} \quad (1.17)$$
**Sketch of the proof** The algebraic Riccati equation can be found after replacing A, B, C and S by their decomposition in the Riccatti equation (1.16).

In the sequel, we require that the algebraic Riccati equation (1.17) admits a unique solution for all $\lambda \in I_\sigma$ which is checked numerically.

**Remark 2** In this example, $\Phi_Z$ is some matrix of function of $\Lambda$, and so is $\Phi_U^{-1}K$,

$$k(\Lambda) = \Phi_U^{-1}K. \qquad (1.18)$$

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 control operator B is composed with the operator $-\partial_{xx}^2$. This is taken into account in $\Phi_U$ in a manner in which $\Phi_Z^{-1}B\Phi_U$ 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$ and $\Phi_Y$ comes directly from the expression of the inner product $(z, z')_Z = (\Phi_Z^{-1}z, \Phi_Z^{-1}z')_{X^4}$ and from $(z_1, z'_1)_{H_0^2(\Gamma)} = ((\Delta^2)^{\frac{1}{2}} z_1, (\Delta^2)^{\frac{1}{2}} z'_1)_{L^2(\Gamma)}$. For $\Phi_Y$, we start from $B = \Phi_Z b(\Lambda) \Phi_U^{-1}$ and from the relation $(u, u')_Y = (\Phi_U^{-1}u, \Phi_U^{-1}u')_{X^2}$ which implies that $-\partial_{xx}^2/\rho^B = b_{3,1}(\Phi_U)_{1,1}$ and $I/\rho^C = b_{4,2}(\Lambda)(\Phi_U)_{2,2}$. The expression of $\Phi_U$ follows.

## Approximation of the Functions of $\Lambda$

The third step in the method consists in an approximation of a general function of $\Lambda$ by a simpler function of $\Lambda$ easily discretized and implemented in a semi-decentralized architecture. The strategy must be general, and in the same time the approximation must be accurate. A simple choice would be to adopt a polynomial or a rational approximation, but their discretization yields very high errors due to the powers of $\Lambda$. This can be avoided when using the Dunford-Schwartz formula, see [27], representing a function of an operator, because it involves only the operator $(\zeta I - \Lambda)^{-1}$ which may be simply, and accurately approximated. However, this formula requires the function be holomorphic inside an open vicinity of $\sigma$. Since the function is generally not known, this set cannot be easily determined, so we prefer to proceed within two steps. First, the function is approximated through a highly accurate rational approximation, then the Dunford-Schwartz formula is applied to the rational approximation, with a path tracing out an ellipse including $I_\sigma$ but no poles.

Since the interval $I_\sigma$ is bounded, each function $k_{ij}(\lambda)$ have a rational approximation over $I_\sigma$, that we write under a global formulation,

$$k_N(\lambda) = \frac{\sum_{m=0}^{N^N} d_m \lambda^m}{\sum_{m'=0}^{N^D} d_{m'} \lambda^{m'}}, \qquad (1.19)$$

where $d_m, d_{m'}$ are matrices of coefficients and $N = (N^N, N^D)$ is the couple comprised of the matrix $N^N$ of numerator polynomial degrees and the matrix $N^D$ of denominator polynomial degrees. The path $\mathcal{C}$, in the Dunford-Schwartz formula,

$$k_N(\Lambda) = \frac{1}{2i\pi} \int_{\mathcal{C}} k_N(\zeta) (\zeta I - \Lambda)^{-1} d\zeta,$$

is chosen to be an ellipse parameterized by

$$\zeta(\theta) = \zeta_1(\theta) + i\zeta_2(\theta), \text{ with } \theta \in [0, 2\pi].$$

The parametrization is used as a change of variable, so the integral is rewritten on the form $I(g) = \int_0^{2\pi} g(\theta) d\theta$, and may be approximated by a quadrature formula involving $M$ nodes $(\theta_l)_{l=1,...,M} \in [0, 2\pi]$, and $M$ weights $(w_l)_{l=1,...,M}$,

$$I_M(g) = \sum_{l=1}^{M} g(\theta_l) w_l.$$

For each $z \in X^4$ and $\zeta \in \mathcal{C}$, we introduce the four-dimensional vector field

$$v^{\zeta} = -i\zeta' k_N(\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 = \operatorname{Re}(-i \zeta' k_N(\zeta)) z, \\ \zeta_2 v_1^\zeta + \zeta_1 v_2^\zeta - \Lambda v_2^\zeta = \operatorname{Im}(-i \zeta' k_N(\zeta)) z. \end{array} \right. \qquad (1.20)$$
Thus, combining the rational approximation $k_N$ and the quadrature formula yields an approximate realization $k_{N,M}(\Lambda)$ of $k(\Lambda)$,

$$k_{N,M}(\Lambda)z = \frac{1}{2\pi} \sum_{l=1}^{M} v_{1}^{\zeta(\theta_l)} w_l. \qquad (1.21)$$

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_{N,M}(\Lambda)z$, requires solving $M$ systems like (1.20) corresponding to the $M$ nodes $\zeta(\theta_l)$. The matrices $k_N(\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.

## Spatial Discretization

The final step consists in a spatial discretization of Equation (1.20), it does not represent a specific novelty, so we do not discuss it through numerical simulations. For the sake of simplicity, the interval $\Gamma$ being meshed with regularly spaced nodes separated by a distance $h$, we introduce $\Lambda_h^{-1}$ the finite difference discretization of $\Lambda^{-1}$,

$$\Lambda_h^{-1} = \frac{1}{h^4} \begin{pmatrix} h^4 & 0 & 0 & 0 & 0 & \dots & 0 \\ -\frac{3}{2}h^3 & 2h^3 & -\frac{1}{2}h^3 & 0 & 0 & \dots & 0 \\ 1 & -4 & 6 & -4 & 1 & \dots & 0 \\ 0 & \dots & \dots & \dots & \dots & \dots & 0 \\ \vdots & \dots & \dots & \dots & \dots & \dots & \vdots \\ 0 & \dots & 0 & 1 & -4 & 6 & -4 & 1 \\ 0 & \dots & 0 & 0 & 0 & -\frac{1}{2}h^3 & 2h^3 & -\frac{3}{2}h^3 \\ 0 & \dots & 0 & 0 & 0 & 0 & 0 & \hbar^4 \end{pmatrix}$$

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,

$$\left\{ \begin{aligned} \zeta_1 v_{1,h}^{\zeta} - \zeta_2 v_{2,h}^{\zeta} - \Lambda_h v_{1,h}^{\zeta} &= \mathrm{Re}(-i\zeta' k_N(\zeta)) z_h, \\ \zeta_2 v_{1,h}^{\zeta} + \zeta_1 v_{2,h}^{\zeta} - \Lambda_h v_{2,h}^{\zeta} &= \mathrm{Im}(-i\zeta' k_N(\zeta)) z_h. \end{aligned} \right.$$

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

$$\Phi_{U,h} k_{N,M,h} z_h = \Phi_{U,h} \frac{1}{2\pi} \sum_{l=1}^{M} v_{1,h}^{\zeta_l} w_l,$$

estimated at mesh nodes, where $\Phi_{U,h}$ is the discretization of $\Phi_U$ which requires the discretization of $-\partial_{xx}^2$ which can be done as for $\Lambda$ by using a finite difference method.

## 1.2.6 Numerical Validation

To build a rational interpolation $k_N$ of the form (1.19) over $I_\sigma$, we mesh the interval with $L+1$ distinct nodes $\lambda_0, \dots, \lambda_L$. Then all $p(\lambda_n)$ solutions of the algebraic Riccati equation are accurately computed with a standard solver. Computing the rational approximation start by imposing $L+1$ conditions

$$k_N(\lambda_n) = k(\lambda_n),$$

or equivalently that

$$\sum_{m=0}^{N^N} d_m \lambda_n^m - k(\lambda_m) \sum_{m'=0}^{N^D} d_{m'}^{m'} \lambda_n^{m'} = 0,$$

for $n=0, \dots, L+1$. Then, when $L$ is large enough, the resulting system with $N^N+N^D+2$ unknowns, $[d, d'] = [d_0, \dots, d_{N^N}, d'_0, \dots, d'_{N^D}]$ is overdetermined, so it is solved in the mean square sense.

In a numerical experiment, we have set all coefficients $R^B, \rho^B, \ell_C, R^C, \rho^C$ and $L_C$ to one, and $L_B = 4.73$. Thus, all eigenvalues of $\Lambda$ turns to be included in $(0, 1)$, the first cantilever eigenvalue turns to be equal to $\lambda_1^C = 12.36, \bar{\psi}_1 = -0.78$ and $c_1 = 9.68$. Moreover, we have chosen $L = 100$ nodes logarithmically distributed along $I_\sigma = (10^{-2}, 1)$. We remark that the shapes of all spectral functions $k_{ij}$ involved in $K$, represented in Figure 1.14, exhibit a singular behavior at the
<figure><img src="image_15.png" /><figcaption>Figure 1.14: Shapes of the Spectral Functions k</figcaption></figure>

<table><caption>Table 1.1: Errors in Rational Approximations</caption><thead><tr><th>(i, j)</th><th>N<sub>ij</sub></th><th>e<sub>ij</sub> × 10<sup>-7</sup></th></tr></thead><tbody><tr><td>(1, 1)</td><td>(7, 19)</td><td>4.78</td></tr><tr><td>(1, 2)</td><td>(7, 20)</td><td>0.69</td></tr><tr><td>(1, 3)</td><td>(13, 8)</td><td>3.83</td></tr><tr><td>(1, 4)</td><td>(7, 19)</td><td>1.19</td></tr><tr><td>(2, 1)</td><td>(8, 20)</td><td>1.81</td></tr><tr><td>(2, 2)</td><td>(7, 19)</td><td>1.19</td></tr><tr><td>(2, 3)</td><td>(20, 10)</td><td>0.89</td></tr><tr><td>(2, 4)</td><td>(19, 7)</td><td>0.53</td></tr></tbody></table>
origin. This shows that this example is by no means trivial. In Table 1.1, we report polynomial degrees $N = (N^N, N^D)$ and relative errors

$$e_{ij} = \frac{\|k_{ij,N} - k_{ij}\|_{L^2(I_\sigma)}}{\|k_{ij}\|_{L^2(I_\sigma)}}$$

between the exact $k$ and its rational approximation $k_N$. The degrees $N^N$ and $N^D$ can be chosen sufficiently large so that errors are sufficiently small, since this has no effect on on-line control computation time. Numerical integrations

<figure><img src="image_16.png" /><figcaption>Figure 1.15: Errors between $k$ and $k_{N,M}$</figcaption></figure>

have been performed with a standard trapezoidal quadrature rule. Relative errors, between the exact functions and final approximations,

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

are reported in Figure 1.15, in logarithmic scale, for $M$ varying from 10 to $10^3$. The results are satisfactory. Accuracy is proportional to the number of nodes. So it may be easily tuned without changing spatial complexity governed by the operator $\Lambda$.

## 1.3 Simulation of Large-Scale Periodic Circuits by a Homogenization Method

This section focuses on the simulation of spatially periodic circuits that may come, for instance, from realization of our control approximations. The periodic unit cell is limited to linear and static components but its number can be very large. Our theory allows one to simulate arrays of electronic circuits which are far away from the possibility of regular circuit simulators like Spice. It is an adaptation of the two-scale approach used in Section 1.1 and has been introduced in [19]. The resulting model consists in a partial differential equations (PDE), related to a macroscopic electric potential, coupled with local circuit equations. In the following, we present the general framework illustrated through a simple example. The numerical resolution of the PDE can be done with usual computational tools. Soving this PDE and postprocessing its solution leads to an approximation of all voltages and currents. Theoretically, more the number of cells is large, more the model is accurate. The method is illustrated on a basic circuit to allow hand calculations, which are mostly matrix multiplications.

### 1.3.1 Linear Static Periodic Circuits

We consider the class of periodic circuits in $d$ space dimensions. An example of such circuit in two space dimensions is shown in Figure 1.16. The circuit cell is detailed on Figure 1.17. Some voltage or current sources, whose value may be zero, are placed on the boundary to realize specific boundary conditions. We assume that the number of cells is large in all the $d$ directions. Mathematically, it is easier to formalize the problem by considering that the whole circuit occupies a unit square $\Omega = (0,1)^d$ and that the period lengths, in all directions, are equal to an identical small parameter $\epsilon$ (cf. Figure 1.16).

We limit ourselves to the study of circuits whose cell is linear and static. Precisely, the components of a cell are limited to the Spice elements $R, V, I, E, F, G, H$. All ports of any multiport component $E, F, G, H$ must belong to a same cell.
<figure><img src="image_17.png" /><figcaption>Figure 1.16: Circuit example.</figcaption></figure>

The expanded cell is arbitrarily defined in a unit cell $Y = (-1/2, +1/2)^d$ (see Figure 1.17). We map any discrete node $n$ onto a continuous coordinates $(y_1, \dots, y_d)$. The vector $y(n) \in \mathbb{R}^d$ is the coordinate vector of a node $n$. For example, the coordinates of the nodes in Figure 1.17 are $y(1, \dots, 6)$,

$$y(1, \dots, 6) = \begin{pmatrix} -1/2 & 0 & 1/2 & 0 & 0 & 1/4 \\ 0 & 0 & 0 & 1/2 & -1/2 & -1/4 \end{pmatrix}.$$

In particular, the coordinates of the node $n=3$ is the vector $(1/2, 0)^T$.

<figure><img src="image_18.png" /><figcaption>Figure 1.17: Expanded cell of the circuit.</figcaption></figure>

The maps of voltages and currents from the whole circuit (global network) to the cell circuit (local network) are defined as follows. First, we denote by

$\mathcal{E}$ = the branch set of the whole circuit,

$\mathcal{N}$ = the node set of the whole circuit,

$E$ = the branch set of the unit cell circuit,

$N$ = the node set of the unit cell circuit,

and we define three indices

* the global index $\mathcal{I}$ references all the branches of the whole circuit,

* the multi-integer $\mu = (\mu_1, \dots, \mu_d) \in \{1, \dots, m\}^d$ enumerates all the cells $Y_\mu^\epsilon$ in the circuit $\Omega$.

* the local index $j \in \{1, \dots, |E|\}$ enumerates all the branches of the unit cell $Y$.

Each branch voltage or current can then be referenced by the index $\mathcal{I}$ or by the couple $(\mu, j)$. This is a one-to-one correspondence denoted by $\mathcal{I} \sim (\mu, j)$. Using this correspondence, for each vector $\mathbf{u} \in \mathbb{R}^{|\mathcal{I}|}$, one may define a unique tensor $U_{\mu j}$ with $(\mu, j) \in \{1, \dots, m\}^d \times \{1, \dots, |E|\}$ by $U_{\mu j} = \mathbf{u}_{\mathcal{I}}$ for $(\mu, j) \sim \mathcal{I}$.
### 1.3.2 Circuit Equations

The electrical state of a circuit can be characterized [7] by the vectors $(\phi, v, i)$ where,

$\phi \in \mathbb{R}^{|\mathcal{N}|}$ = the nodal voltages (or electric potentials),

$v \in \mathbb{R}^{|\mathcal{E}|}$ = the branch voltages,

$i \in \mathbb{R}^{|\mathcal{E}|}$ = the branch currents.

We can formulate the circuit equations under the form

$$v = \mathcal{A}^T \phi, \qquad (1.22)$$

$$\mathcal{R}i + \mathcal{M}v = u_s, \qquad (1.23)$$

$$i^T w = 0, \qquad (1.24)$$

$$\text{for all } w = \mathcal{A}^T \psi \text{ with } \psi \in \Psi. \qquad (1.25)$$

where $u_s \in \mathbb{R}^{|\mathcal{E}|}$ represents voltage and current sources merged in single vector completed by some zeros. Equation (1.22) is the Kirchhoff's Voltage Law. Equation (1.23) represents the constitutive equations and Equations (1.24, 1.25) correspond to the Tellegen theorem. Here $\Psi$ is the set of admissible potentials for the circuit problem, that is to say

$$\Psi = \{ \psi \in \mathbb{R}^{|\mathcal{N}|} \text{ such that } \psi_{\mathcal{J}} = 0 \text{ for all ground nodes } n_{\mathcal{J}} \}.$$

As the matrices $\mathcal{M} \in \mathbb{R}^{|\mathcal{E}|} \times \mathbb{R}^{|\mathcal{E}|}$, $\mathcal{R} \in \mathbb{R}^{|\mathcal{E}|} \times \mathbb{R}^{|\mathcal{E}|}$ and the vector $u_s \in \mathbb{R}^{|\mathcal{E}|}$ are exclusively deduced from the branch equations of the circuit, they can be expressed in terms of two reduced matrices $M \in \mathbb{R}^{|E|} \times \mathbb{R}^{|E|}$ and $R \in \mathbb{R}^{|E|} \times \mathbb{R}^{|E|}$ and a reduced vector $u_s \in \mathbb{R}^{|E|}$. The reduced matrices and vector are simply derived from the constitutive equations of the unit cell, which are in the example,

$$-v_1 + r i_1 = 0,$$

$$-v_2 + r i_2 = 0,$$

$$-v_3 + r i_3 = 0,$$

$$-v_4 + r i_4 = 0,$$

$$i_5 = i_s.$$

The transpose $\mathcal{A}^T \in \mathbb{R}^{|\mathcal{E}|} \times \mathbb{R}^{|\mathcal{N}|}$ of the incidence matrix can also be expressed in terms of a reduced matrix noted by $A^T$ (with a little abuse of notation). Notice that we cannot find a reduced matrix for the incidence matrix itself. We introduce the local (complete) incidence matrix $A \in \mathbb{R}^{|\mathcal{N}|} \times \mathbb{R}^{|E|}$,

$$A_{i,j} = \begin{cases} +1 & \text{if branch } j \text{ leaves node } i, \\ -1 & \text{if branch } j \text{ enters node } i, \\ 0 & \text{if branch } j \text{ does not touch node } i. \end{cases}$$

The solution of the simplified model introduced in this section realizes an approximation of the solution of (1.22–1.25) for small values of $\varepsilon$ ($\varepsilon \ll 1$). It is derived as a limit of the latter when the cell length $\varepsilon$ diminishes towards zero.

### 1.3.3 Direct Two-scale Transform $T_E$

The general idea of the two-scale transform rests on gathering the voltages or currents both denoted by $\hat{u}$ of a same branch $j$ of all cells. Indeed, the voltages or currents of all branches $j$ are defined by a function $x \mapsto \hat{u}_j(x)$, which depends on the parameter $\varepsilon$ and whose limit when $\varepsilon \to 0$ will be calculated.

Let us first denote by $\chi_{Y_\mu^\varepsilon}(x)$ the characteristic function of the cell $Y_\mu^\varepsilon$ equals to 1 when $x \in Y_\mu^\varepsilon$ and 0 otherwise. As an example, the characteristic function $\chi_{Y_{32}^\varepsilon}$ of the cell $\mu = (3, 2)$ is represented in Figure 1.18.

The two-scale transform $\hat{u}$ of the branch vector $u \in \mathbb{R}^{|\mathcal{E}|}$ is the vector of $Y_\mu^\varepsilon$-piecewise constant functions $\hat{u} \in \mathbb{P}^0(\Omega)^{|E|}$,

$$\hat{u}_j(x) = \sum_{\mu \in \{1, \dots, m\}^d} \chi_{Y_\mu^\varepsilon}(x) U_{\mu j}. \qquad (1.26)$$

where $U_{\mu j} = u_j$ with $(\mu, j) \sim \mathcal{J}$. For example, $\hat{v}_j(x)$ is the voltage $V_{\mu j}$ of the branch referred by the local index $j$ of the cell $\mu$ to which $x$ belongs. By construction, the function $x \mapsto \hat{v}_j(x)$ is constant on all cells. Figure 1.19 illustrates this concept by representing a component of an arbitrary vector $\hat{v}$. It indicates that the voltage $\hat{v}_2(x)$ related to the branch $j=2$ (cf. Figure 1.17) of the cell $(\mu_1, \mu_2) = (1, 4)$ is equal to 2V.

We denote by $T_E$ the linear map $u \mapsto \hat{u}$ from $\mathbb{R}^{|\mathcal{E}|}$ to $\mathbb{P}^0(\Omega)^{|E|} \subset L^2(\Omega)^{|E|}$, with $\mathbb{P}^0(\Omega)$ the set of piecewise constant functions over the cells. Our model is derived for the limit, when $\varepsilon \to 0$, of all vectors involved in the circuit equations. The actual circuit voltages and currents are then computed by inverting the two-scale transform with the physical value of $\varepsilon$. Next Section is devoted the construction of $T_E^{-1}$.
<figure><img src="image_19.png" /><figcaption>Figure 1.18: Characteristic Function $\chi_{Y_{32}^\varepsilon}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 1.19: One component $\hat{v}_2(x)$ of a two-scale transform.</figcaption></figure>

### 1.3.4 Inverse Two-scale Transform $T_E^{-1}$

The calculation of the inverse two-scale transform $T_E^{-1}$ is done by computing the adjoint $T_E^*$ and then proving two identity properties between these transforms. The notation for inner products and norms in $\mathbb{R}^{|\mathcal{E}|}$ and $L^2(\Omega)^{|E|}$ are collected in the following table,

<table><thead><tr><th></th><th><b>u, v</b> &isin; &#8477;<sup>|&mathcal{E}|</sup></th><th><i>u</i>, <i>v</i> &isin; L<sup>2</sup>(&Omega;)<sup>|E|</sup></th></tr></thead><tbody><tr><td>Inner product</td><td>[<b>u</b>, <b>v</b>]</td><td>(<i>u</i>, <i>v</i>)</td></tr><tr><td>Norm</td><td>|<b>v</b>| = [<b>v</b>, <b>v</b>]<sup>1/2</sup></td><td>||<i>u</i>|| = (<i>u</i>, <i>u</i>)<sup>1/2</sup></td></tr></tbody></table>

and defined by

$$[\mathbf{u}, \mathbf{v}] = \varepsilon^d \mathbf{u}^T \cdot \mathbf{v},$$

$$(u, v) = \sum_{j=1}^{|E|} \int_{\Omega} u_j(x) v_j(x) dx.$$

For all $u \in L^2(\Omega)^{|E|}$, the adjoint $T_E^* u$ is defined through the equality

$$[T_E^* u, v] = (u, T_E v) \quad \text{for all } v \in \mathbb{R}^{|\mathcal{E}|}. \qquad (1.27)$$
The calculation of $T_E^*$ from (1.27) is given in Appendix and leads to

$$ (T_E^* u)_\mathcal{F} = \varepsilon^{-d} \int_{Y_\mu^\varepsilon} u_j(x) dx. \qquad (1.28) $$

Moreover, Appendix proves that $T_E^*T_E = I_\mathcal{F}$ on $\mathbb{R}^{|\mathcal{F}|}$ and $T_E T_E^* = I_E$ on $\mathbb{P}^0(\Omega)^{|E|}$. As $T_E$ is one-to-one from $\mathbb{R}^{|\mathcal{F}|}$ to $\mathbb{P}^0(\Omega)^{|E|}$, these two identities show that

$$ T_E^{-1} = T_E^*. $$

### 1.3.5 Two-scale transform $T_N$

The two-scale transform $\hat{\varphi}$ of the nodal vector $\varphi \in \mathbb{R}^{|\mathcal{N}|}$ is the vector of $Y_\mu^\varepsilon$-piecewise constant functions $\hat{\varphi} \in \mathbb{P}^0(\Omega)^{|\mathcal{N}|}$ defined by

$$ \hat{\varphi}_j(x) = \sum_{\mu \in \{1, \dots, m\}^d} \chi_{Y_\mu^\varepsilon}(x) \Phi_{\mu j}, \qquad (1.29) $$

where $\Phi_{\mu j} = \varphi_\mathcal{F}$ with $(\mu, j) \sim \mathcal{F}$. We denote by $T_N$ the linear map $\varphi \mapsto \hat{\varphi}$ from $\mathbb{R}^{|\mathcal{N}|}$ to $\mathbb{P}(\Omega)^{|\mathcal{N}|} \subset L^2(\Omega)^{|\mathcal{N}|}$. As the nodes located on the cell boundary belong to two adjacent cells $T_N$ is not one to one.

### 1.3.6 Behavior of “Spread” Analog Circuits

We start by illustrating the scaling of currents and voltages in a one-dimensional circuit. A circuit spread out over a large region may have some pathes linking oposite sides. In view of deriving a partial differential equation for the electric potential, we assume that voltages are increments of the order $\varepsilon$ along such pathes. Flowing current result of numerous $(1/\varepsilon)$ additive sources, so it has a magnetude of 1 as soon as the sources are of the order of $\varepsilon$ (Figure 1.20).

<figure><img src="image_21.png" /><figcaption>Figure 1.20: Illustration of the magnitude orders of $\varphi$ and $v$.</figcaption></figure>

A branch which does not belongs to any crossing path is necessarily part of a path to the ground, so its voltage magnitude is 1. We choose its magnitude current be of the order $\varepsilon$ as it may be a source for a crossing path (Figure 1.21).

The periodicity of the circuit implies that each node $n$ located on the boundary of the unit cell has its counterpart $n'$ on the opposite side. We assume that each such couple is linked by at least a crossing path. We introduce the set $E_C \subset E$ constituted of all the branches of at least one path linking each couple $(n, n')$. Of course, a link between $(n, n')$ which includes a ground node is not considered as a path. The complementary set $E - E_C$ is denoted by $E_{NC}$ (non-crossing pathes).

In the case where many crossing pathes are linking $n$ and $n'$, the designer is free to decide which are included in $E_C$ and which are not, with regard to the above discussion about current and voltage magnitudes. The subset $E_C$ is partitioned in its $n_c$ connected components $E_C = \cup_{p=1}^{n_c} E_{Cp}$. In the following, the main result on the circuit equations will be derived for the connected components of $E_C$ and not for $E_C$ itself.

The subsets $N_C$ and $N_{NC}$ of $N$ are defined as the set of nodes involved in at least one of the branch of $E_C$ and $E_{NC}$ respectively. Since the branches of $E_C$ and $E_{NC}$ have common nodes, these two subsets $N_C$ and $N_{NC}$ are not a partition
<figure><img src="image_22.png" /><figcaption>Figure 1.21: Illustration of the magnitude order of i.</figcaption></figure>

of $N$, $N_C \cap N_{NC} \neq \emptyset$. $N_{Cp}$ is the set of nodes involved in the branches of $E_{Cp}$. The subsets $N_{Cp}$ constitute a partition of $N_C$, $N_C = \cup_{p=1}^{n_c} N_{Cp}$, as the subsets $E_{Cp}$ constitute a partition of $E_C$ and they are disconnected. The definitions and some properties of these sets are depicted on Figure 1.22.

<figure><img src="image_23.png" /><figcaption>Figure 1.22: Depiction of node and branch sets.</figcaption></figure>

Finally, $N_{C \times 1}$ is a set of $n_c$ nodes constituted of one arbitrary node of each connected component $N_{Cp}$ and $N_0$ is the set containing only the cell ground node.

### 1.3.7 Cell Equations (Problem Micro)

The model formulation is decomposed in four parts. We first formulate the linear relation between mean electric potentials $\phi_C^0$ along crossing paths and the other fields as branch currents and voltages. This relation is strictly local in each cell. In the next Section, the linear relation is simply rewritten introducing linear operators. They are then used for coefficients of the boundary value problem on $\phi_C^0$. Finally, actual voltages and currents are computed thanks to the inverse two-scale transform.

The previous assumptions about voltage and current magnitudes is formulated using the scaling matrices $S_v$, $S_c$ and $S_s$ applied to the two-scale transforms,

$$ \hat{\mathbf{i}}^{\epsilon} = S_c \hat{\mathbf{i}}, $$

$$ \hat{\mathbf{v}}^{\epsilon} = S_v \hat{\mathbf{v}}, $$

$$ \hat{\mathbf{u}}_s^{\epsilon} = S_s \hat{\mathbf{u}}_s. $$

with the $|E| \times |E|$ scaling matrices defined as

$$ S_v = \epsilon^{-1} I_{E_C} + I_{E_{NC}}, \qquad (1.30) $$

$$ S_c = I_{E_C} + \epsilon^{-1} I_{E_{NC}}, \qquad (1.31) $$

$$ S_s = \Pi_c S_c + \Pi_v S_v. \qquad (1.32) $$

Here the $|E| \times |E|$ matrices $I_{E_C}$ and $I_{E_{NC}}$ are the projectors on the subspaces of $\mathbb{R}^{|E|}$ generated by vectors with non vanishing
values on $E_C$ and $E_{NC}$,

$$ (I_{E_C})_{jk} = \begin{cases} \delta_{jk} & \text{if } e_j \in E_C, \\ 0 & \text{otherwise,} \end{cases} $$

$$ (I_{E_{NC}})_{jk} = \begin{cases} \delta_{jk} & \text{if } e_j \in E_{NC}, \\ 0 & \text{otherwise.} \end{cases} $$

Moreover, each branch equation in (1.23) is homogeneous to a current or to a voltage, this leads to a partition of E into two subsets. The $|E| \times |E|$ matrices $\Pi_c$ and $\Pi_v$ (for currents and voltages respectively) are defined as the projectors on these two subsets.

The transform $\hat{\varphi}$ is not scaled, it is rewritten as $\hat{\varphi}^{\varepsilon}$ for the sake of notation uniformity,

$$ \hat{\varphi}^{\varepsilon} = \hat{\varphi}. $$

The reduced matrices $M$ and $R$ of $\mathcal{M}$ and $\mathcal{R}$ are scaled in a consistent manner,

$$ M^{\varepsilon} = S_s M S_v^{-1}, $$

$$ R^{\varepsilon} = S_s R S_c^{-1}. $$

The scaled reduced matrices $M^{\varepsilon}$ and $R^{\varepsilon}$ are assumed to converge towards some limit $M^0$ and $R^0$ when $\varepsilon \to 0$. If the norms $||\hat{\mathbf{i}}^{\varepsilon}||$, $||\hat{\mathbf{v}}^{\varepsilon}||$, $||\hat{\varphi}^{\varepsilon}||$ and $||\hat{\mathbf{u}}_s^{\varepsilon}||$ are bounded then $(\hat{\mathbf{i}}^{\varepsilon}, \hat{\mathbf{v}}^{\varepsilon}, \hat{\varphi}^{\varepsilon}, \hat{\mathbf{u}}_s^{\varepsilon})$ is weakly converging when $\varepsilon \to 0$ towards a limit $(i^0, v^0, \varphi^0, u_s^0)$ in $L^2(\Omega)$ [27].

We prove that the vector of electric potential $\varphi^0(x)$ is a constant $\varphi_{C_p}^0(x)$ in each connected component of cell crossing pathes. So, we split it according to $\varphi^0 = I^0 \varphi_C^0 + \varphi_{NC}^0$, $I^0$ being defined at (1.39), with $\varphi_C^0 = (\varphi_{C_p}^0)_{p=1,...,n_c}$ and $\varphi_{NC}^0(x)$ being the electric potentials at nodes not in crossing pathes. Electric potential variations within connected components of crossing pathes are recovered thanks to the corrector $\varphi_C^1$,

$$ \varphi_C = \varphi^0 + \varepsilon \varphi_C^1. $$

We are ready to state the cell equations. We begin with assuming that $\varphi_C^0$ is known.

For given $\varphi_C^0 \in \Psi^H$, $\Psi^H$ defined in (1.47), and $u_s^0 \in L^2(\Omega)^{|E|}$ there exist $\varphi_C^1 \in L^2(\Omega; \mathbb{R}_{per}^{|N|})$ such that $\varphi_{NC}^0 \in L^2(\Omega; \mathbb{R}^{|N|})$, $i^0 \in L^2(\Omega)^{|E|}$ and $v \in L^2(\Omega)^{|E|}$ are solution of the algebraic cell circuit equations at each location of $\Omega$,

$$ v = I_{EC} A^T \varphi_C^1 + I_{ENC} A^T \varphi_{NC}^0, \quad (1.33) $$

$$ R^0 i^0 + M^0 v = u_s^0 - M^0 (\tau \nabla \varphi_C^0 + I_{ENC} A^T I^0 \varphi_C^0), \quad (1.34) $$

$$ I^{0T} w = 0, \quad (1.35) $$

$$ \text{for all } w = I_{EC} A^T \psi_C^1 + I_{ENC} A^T \psi_{NC}^0 \text{ with } (\psi_C^1, \psi_{NC}^0) \in \Psi^m. \quad (1.36) $$

Moreover, the vector $v^0 \in L^2(\Omega)^{|E|}$ is expressed by

$$ v^0 = v + \tau \nabla \varphi_C^0 + I_{ENC} A^T I^0 \varphi_C^0. $$

We assume that the solution is unique. This assumption is generally satisfied once the global circuit equations has a unique solution.

The admissible nodal voltage set is

$$ \Psi^m = \{ (\psi_C^1, \psi_{NC}^0) \in L^2(\Omega; \mathbb{R}_{per}^{|N|}) \times L^2(\Omega; \mathbb{R}^{|N|}) \text{ such that } I_{(N-N_C) \cup N_{C \times 1}} \psi_C^1 = 0 \text{ and } I_{N_C \cup N_0} \psi_{NC}^0 = 0 \}. $$

The set $\mathbb{R}_{per}^{|N|}$ is defined as,

$$ \mathbb{R}_{per}^{|N|} = \{ \phi \in \mathbb{R}^{|N|} \text{ such that } \phi_j = \phi_{j'} \text{ for all couple } (n_j, n_{j'}) \text{ of opposite nodes} \}. $$

The $|E| \times d \times n_c$ tensor $\tau$ is defined by

$$ \tau_{lkp} = \begin{cases} \sum_{j:n_j \in N_{C_p}} y_k(n_j) A_{jl} & \text{for } e_l \in E_{C_p}, \\ 0 & \text{otherwise.} \end{cases} \quad (1.37) $$

We recall that $y(n) \in \mathbb{R}^d$ is the coordinate vector of a node $n$, and we use the tensor product notation,

$$ (\tau \theta)_l = \sum_k \sum_p \tau_{lkp} \theta_{kp}, \quad (1.38) $$
where the summation is on the two last indices of $\tau$.

The $|N| \times n_c$ matrix $I^0$ is defined by

$$I_{jp}^{0} = \begin{cases} 1 & \text{if } n_j \in N_{C_p}, \\ 0 & \text{otherwise,} \end{cases} \qquad (1.39)$$

$N_{C_p}$ is the set of nodes involved in the branches of $E_{C_p}$.

### 1.3.8 Reformulation of the Problem Micro

In the next section, we state that $\varphi_C^0$ is the solution of a partial differential equation, so once $\varphi_C^0$ is known, $i^0$ and $v^0$ can be computed. The equations (1.33–1.36) being linear, there exists some matrices $L_x$, $H_x$ and a third order tensor $P_x$ such that $i^0$, $\varphi_{NC}^0$ and $v$ that can be expressed as function of $\varphi_C^0$, its gradient $\nabla\varphi_C^0$ and the vector source $u_s^0$,

$$i^0 = \mathcal{L}_i \varphi_C^0 + \mathcal{P}_i \nabla \varphi_C^0 + \mathcal{H}_i u_s^0, \qquad (1.40)$$

$$\varphi_{NC}^{0} = \mathcal{L}_{\varphi} \varphi_{C}^{0} + \mathcal{P}_{\varphi} \nabla \varphi_{C}^{0} + \mathcal{H}_{\varphi} u_{s}^{0}, \qquad (1.41)$$

$$v = \mathcal{L}_v \varphi_C^0 + \mathcal{P}_v \nabla \varphi_C^0 + \mathcal{H}_v u_s^0. \qquad (1.42)$$

The computation of the vector $v^0$ is then unchanged,

$$v^0 = v + \tau \nabla \varphi_C^0 + I_{E_{NC}} A^T I^0 \varphi_C^0. \qquad (1.43)$$

The terms $L_\alpha$, $P_\alpha$ and $H_\alpha$ are of course independent of $\varphi_C^0$, $\nabla\varphi_C^0$ and $u_s^0$ and therefore can be used to express the coefficients in the equation of $\varphi_C^0$.

### 1.3.9 Homogenized Circuit Equations (Problem Macro)

In this Section, we state the equation satisfied by $\varphi_C^0$.

The vector $\varphi_C^0 \in \Psi^H$ is solution of the $n_c$ partial differential equations, so-called homogenized equations, with its boundary conditions,

$$A^H (\mathcal{P}_i \nabla \varphi_C^0 + \mathcal{L}_i \varphi_C^0) = -A^H \mathcal{H}_i u_s^0, \qquad (1.44)$$

$$\varphi_{C_p}^0 = 0 \text{ on } \Gamma_{0p}, \qquad (1.45)$$

$$(\mathcal{P}_i \nabla \varphi_C^0 + \mathcal{L}_i \varphi_C^0) n_\tau = 0 \text{ on } \Gamma - \Gamma_{0p}. \qquad (1.46)$$

$\Gamma_{0p}$ is the part of the boundary $\Gamma$ of $\Omega$ where the $p^{th}$ connected component is grounded. The operator $A^H$ is defined by

$$A^H = -\partial_{\tau^*} + I^{0T} A I_{E_{NC}},$$

where $\partial_{\tau^*} i = \tau^* \nabla i$ with $\tau_{pkl}^* = \tau_{lkp}$ and the use of convention (1.38). The derivative $\partial_\tau \varphi_C^0$ and the normal $n_\tau$ are defined by

$$\partial_\tau \varphi_C^0 = \tau \nabla \varphi_C^0,$$

$$(n_\tau)_{lp} = \sum_{k=1}^{d} \tau_{lkp} n_k,$$

$\nabla$ being the gradient $(\partial_{x_k})_{k=1..d}$ and $n = (n_k)_{k=1..d}$ being the outward unit normal vector to the boundary $\Gamma$ of $\Omega$. Remark that the coefficients $A^H$ and the derivatives $\nabla_\tau$ depend on node coordinates inherited from the expression (1.37) of $\tau$. Finally, the admissible set of macroscopic potential is

$$\Psi^H = \{ \psi \in L^2(\Omega)^{n_c} \text{ such that } \partial_\tau \psi \in L^2(\Omega)^{|E|} \text{ and } \psi_k(x) = 0 \text{ on } \Gamma_{0k} \}. \qquad (1.47)$$

Let us turn to the example depicted on Figure 1.16, $\varphi_{C1}^0$ has only one component $\varphi_{C1}^0$ ($n_c = 1$), we assume that $r = \varepsilon r^0$ and $i_s = \varepsilon i_s^0$. $\varphi_C^0$ is then solution of the partial differential equation

$$\frac{\partial^2 \varphi_{C1}^0}{\partial x_1^2} + \frac{\partial^2 \varphi_{C1}^0}{\partial x_2^2} = -2r^0 i_s^0 \text{ in } \Omega,$$

$$\varphi_{C1}^0 = 0 \text{ on } \Gamma_{0,1},$$

$$\nabla \varphi_{C1}^0 \cdot n_\tau = 0 \text{ on } \Gamma - \Gamma_{0,1}.$$
Once the solution $\varphi_{C1}^0$ is computed, the two-scale limits $(v^0, i^0)$ are expressed by (1.40)–(1.43). For the example, the two-scale current and voltage are given by

$$ 
\begin{aligned}
v^0 &= -\frac{1}{2} \left( \frac{\partial \varphi_{C1}^0}{\partial x_1}, \frac{\partial \varphi_{C1}^0}{\partial x_1}, \frac{\partial \varphi_{C1}^0}{\partial x_2}, \frac{\partial \varphi_{C1}^0}{\partial x_2}, 2\varphi_{C1}^0 \right)^T, \\
i^0 &= -\frac{1}{2} \left( \frac{1}{r^0} \frac{\partial \varphi_{C1}^0}{\partial x_1}, \frac{1}{r^0} \frac{\partial \varphi_{C1}^0}{\partial x_1}, \frac{1}{r^0} \frac{\partial \varphi_{C1}^0}{\partial x_2}, \frac{1}{r^0} \frac{\partial \varphi_{C1}^0}{\partial x_2}, -2i_s^0 \right)^T.
\end{aligned}
$$

### 1.3.10 Computation of Actual Voltages and Currents

Actual voltages and currents may be then recovered through the inverse two-scale transform (1.28) and inverse scalings (1.30–1.32),

$$ 
\begin{aligned}
\mathbf{v} &\approx T_E^{-1} S_v^{-1} v^0, \\
\mathbf{i} &\approx T_E^{-1} S_i^{-1} i^0.
\end{aligned}
$$

The approximation of the node voltages $\varphi$ is realized on a different manner. From the solution $\varphi_C^0$ of (1.44–1.46) the approximation of the node voltages is built as follows (neglecting the first order correction),

$$ \varphi_{\mathcal{I}} \approx \varphi_C^0(x_{\mu}^{\varepsilon} + \varepsilon y(n_j)) \text{ for } \mathcal{I} \sim (\mu, j) \text{ and } n_j \in N_C, $$

where $x_{\mu}^{\varepsilon} \in [0, 1]^d$ is the coordinate vector of the center of the cell $\mu$ (Figure 1.16) and we recall that $y(n_j) \in [-1/2, 1/2]^d$ is the coordinate vector of the node $n_j$ (Figure 1.17). For the voltage at the nodes which belong to a non crossing path, the approximation is more simple because it does not need to refer to the spatial location in the cell

$$ \varphi_{\mathcal{I}} \approx \varphi_C^0(x_{\mu}^{\varepsilon}) \text{ for } \mathcal{I} \sim (\mu, j) \text{ and } n_j \in N - N_C - N_0. $$

A comparison of the solutions computed by the two-scale model with those obtained by a direct circuit simulation is has been done for the example circuit with 10 × 10 cells. In the following simulation results, all the circuit boundaries are connected to ground and the component values are $r^0 = 10 \text{ k}\Omega$ and $i_s^0 = 1 \text{ mA}$. The comparison focus only on node voltages magnitude, their location on the domain being shown on Figure 1.23.

<figure><img src="image_24.png" /><figcaption>Figure 1.23: Location of node voltages.</figcaption></figure>

The computation of $\varphi_C^0$ has been done using a Finite Element Method (FEM) on a regular mesh of 40 × 40 squares (see Figure ??). We must emphasize that the mesh size is chosen to obtain an accurate numerical solution of the PDE and is not related at all to the number of the cells of the circuit. The maximum amplitude is 1.4742 V. A direct simulation of the periodic circuit has been made with Spice. The maximum amplitude is 1.4723V. In Figure ?? the continuous FEM solution $\varphi_{\mathcal{I}} = \varphi_C^0(x_{\mu}^{\varepsilon} + \varepsilon y(n_j))$ is represented by the mesh while all voltage nodes computed by Spice located by bullets. These results show a good qualitative agreement between the two methods even if the number of cells is not large (10 in one direction, which corresponds to $\varepsilon = 0.1$).
<figure><img src="image_25.png" /><figcaption>Figure 1.24: Two-scale model versus Spice computation.</figcaption></figure>

### 1.3.11 Appendices

#### Basic properties of some integrals on cells

$$ \int_{Y_{\mu}^{\epsilon}} dx' = \epsilon^{d} $$

$$ \int_{Y_{\mu}^{\epsilon}} \chi_{Y_{\lambda}^{\epsilon}}(x') dx' = \epsilon^{d} \delta_{\mu\lambda} $$

#### Derivation of the expression of $T_E^*$

$$ (T_E v, u) = \int_{\Omega} (T_E v) \cdot u(x) dx $$

$$ = \sum_{j=1}^{|E|} \int_{Y_{\mu}^{\epsilon}} (T_E v)_j(x) \cdot u_j(x) dx $$

$$ = \sum_{\mu \in \{1, \dots, m\}^d} \sum_{j=1}^{|E|} \epsilon^{-d} \int_{Y_{\mu}^{\epsilon}} \chi_{Y_{\mu}^{\epsilon}}(x) u_j(x) dx V_{\mu j} $$

$$ = \epsilon^d \sum_{\mu \in \{1, \dots, m\}^d} \sum_{j=1}^{|E|} \epsilon^{-d} \int_{Y_{\mu}^{\epsilon}} u_j(x) dx V_{\mu j} $$

$$ [T_E^* u, v] = \epsilon^d (T_E^* u)^T \cdot v $$

$$ = \epsilon^d \sum_{\mu \in \{1, \dots, m\}^d} \sum_{j=1}^{|E|} (T_E^* u)_{\mu j} V_{\mu j} $$

$$ \Rightarrow (T_E^* u)_{\mu j} = \epsilon^{-d} \int_{Y_{\mu}^{\epsilon}} u_j(x) dx $$
**Proof that $T_E^*T_E = I_E$ on $\mathbb{R}^{\mathcal{S}}$**

Let $u \in \mathbb{R}^{\mathcal{S}}$ and $\mathcal{S} \sim (\mu, j)$,

$$ (T_E^* T_E u)_{\mathcal{S}} = T_E^* \left( \sum_{\lambda \in \{1, \dots, m\}^d} U_{\lambda j} \chi_{Y_\lambda^\varepsilon}(x) \right) \\ = \varepsilon^{-d} \int_{Y_\mu^\varepsilon} \sum_{\lambda \in \{1, \dots, m\}^d} \chi_{Y_\lambda^\varepsilon}(x) dx U_{\lambda j} \\ = \varepsilon^{-d} \varepsilon^d U_{\mu j} \\ = u_{\mathcal{S}} $$

**Proof that $T_E T_E^* = I_E$ on $\mathbb{P}^0(\Omega)^{|E|}$**

Let $u \in \mathbb{P}^0(\Omega; \mathbb{R}^{|E|})$,

$$ (T_E T_E^* u)_j(x) = (T_E(\varepsilon^{-d} \int_{Y_\mu^\varepsilon} u_j(x') dx'))_j(x) \\ = \sum_{\mu \in \{1, \dots, m\}^d} \varepsilon^{-d} \int_{Y_\mu^\varepsilon} u_j(x') dx' \chi_{Y_\mu^\varepsilon}(x) \\ = \sum_{\mu \in \{1, \dots, m\}^d} \sum_{\lambda \in \{1, \dots, m\}^d} \varepsilon^{-d} \int_{Y_\mu^\varepsilon} \chi_{Y_\lambda^\varepsilon}(x') dx' U_{\lambda j} \chi_{Y_\mu^\varepsilon}(x) \\ = \sum_{\mu \in \{1, \dots, m\}^d} \sum_{\lambda \in \{1, \dots, m\}^d} \delta_{\mu \lambda} U_{\lambda j} \chi_{Y_\mu^\varepsilon}(x) \\ = \sum_{\mu \in \{1, \dots, m\}^d} U_{\mu j} \chi_{Y_\mu^\varepsilon}(x) \\ = u_j(x) $$
Bibliography

[1] Tarantola A. *Inverse Problem Theory and Methods for Model Parameter Estimation*. SIAM, 2005.

[2] R. J. Allemang. The modal assurance criterion - twenty years of use and abuse. *S V Sound and Vibration*, 37(8):14-23, 2003.

[3] G. Binnig, C.F. Quate, and C. Gerber. Atomic force microscope. *Physical Review Letters*, 56(9):930-3, 1986/03/03.

[4] David Bullen, Sung-Wook Chung, Xuefeng Wang, Jun Zou, Chad A. Mirkin, and Chang Liu. Parallel dip-pen nanolithography with arrays of individually addressable cantilevers. *Applied Physics Letters*, 84(5):789-791, 2004.

[5] David Bullen and Chang Liu. Electrostatically actuated dip pen nanolithography probe arrays. *Sensors and Actuators, A: Physical*, 125(2):504-511, 2006.

[6] M. Carrion-Vazquez, A.F. Oberhauser, T.E. Fisher, P.E. Marszalek, Hongbin Li, and J.M. Fernandez. Mechanical design of proteins-studied by single-molecule force spectroscopy and protein engineering. *Progress in Biophysics; Molecular Biology*, 74(1-2):63-91, 2000.

[7] L.O. Chua, A.D. Desoer, and S.K. Kuh. *Linear and Nonlinear Circuits*. Series in Electrical Engineering. McGraw-Hill, 1987.

[8] R. F. Curtain and H. Zwart. *An introduction to infinite-dimensional linear systems theory*, volume 21 of *Texts in Applied Mathematics*. Springer-Verlag, New York, 1995.

[9] R. Dautray and J.-L. Lions. *Mathematical analysis and numerical methods for science and technology*. Springer-Verlag, Berlin, 1990.

[10] S. Decossas, L. Patrone, A.M. Bonnot, F. Comin, M. Derivaz, A. Barski, and J. Chevrier. Nan manipulation by atomic force microscopy of carbon nanotubes on a nanostructured surface. *Surface Science*, 543(1-3):57-62, 2003.

[11] M. Despont, J. Brugger, U. Drechsler, U. Durig, W. Haberle, M. Lutwyche, H. Rothuizen, R. Stutz, R. Widmer, G. Binnig, H. Rohrer, and P. Vettiger. Vlsi-nems chip for parallel afm data storage. *Sensors and Actuators, A: Physical*, 80(2):100-107, 2000. Thermomechanical writing; Electrical stability; Time-multiplexed method;

[12] D. Drakova. Theoretical modelling of scanning tunnelling microscopy, scanning tunnelling spectroscopy and atomic force microscopy. *Reports on Progress in Physics*, 64(2):205-90, 2001.

[13] R. Garcia and R. Perez. Dynamic atomic force microscopy methods. *Surface Science Reports*, 47(6-8):197-301, 2002.

[14] F.J. Giessibl. Advances in atomic force microscopy. *Reviews of Modern Physics*, 75(3):949-83, 2003.

[15] John-Bruce D. Green and Gil U. Lee. Atomic force microscopy with patterned cantilevers and tip arrays: Force measurements with chemical arrays. *Langmuir*, 16(8):4009-4015, 2000.

[16] Gen-Wen Hsieh, Ching-Hsiang Tsai, We-Chih Lin, Chao-Chiun Liang, and Yu-Wen Lee. Bond-and-transfer scanning probe array for high-density data storage. *IEEE Transactions on Magnetics*, 41(2):989-991, 2005.

[17] Kuniyuki Kakushima, Toshiyuki Watanabe, Kouji Shimamoto, Takushi Gouda, Manabu Ataka, Hidenori Mimura, Yoshimasa Isono, Gen Hashiguchi, Yutaka Mihara, and Hiroyuki Fujita. Atomic force microscope cantilever array for parallel lithography of quantum devices. *Japanese Journal of Applied Physics, Part 1: Regular Papers and Short Notes and Review Papers*, 43(6 B):4041-4044, 2004.

[18] Young-Sik Kim, Hyo-Jin Nam, Seong-Moon Cho, Jae-Wan Hong, Dong-Chun Kim, and Jong U. Bu. Pzt cantilever array integrated with piezoresistor sensor for high speed parallel operation of afm. *Sensors and Actuators, A: Physical*, 103(1-2):122-129, 2003.
[19] M. Lenczner. Homogenization of linear spatially periodic electronic circuits. *Netw. Heterog. Media*, 1(3):467–494 (electronic), 2006.

[20] M. Lutwyche, C. Andreoli, G. Binnig, J. Brugger, U. Drechsler, W. Haberle, H. Rohrer, H. Rothuizen, P. Vettiger, G. Yaralioglu, and C. Quate. 5x5 2d afm cantilever arrays a first step towards a terabit storage device. *Sensors and Actuators, A: Physical*, 73(1-2):89–94, 1999.

[21] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. *JCP J. Chem. Phys.*, 1(6):1087–1092, 1953.

[22] K. Mosegaard and A. Tarantola. *International Handbook of Earthquake and Engineering Seismology (part A), chapitre Probabilistic Approach to Inverse Problems*. Academic Press, 2002.

[23] Matthias Rief and Helmut Grubmller. Single molecule force spectroscopy in biology using the atomic force microscope. *European journal of chemical physics and physical chemistry*, 3(3):255–61, 2002.

[24] Metin Sitti. Atomic force microscope probe based controlled pushing for nanotribological characterization. *IEEE/ASME Transactions on Mechatronics*, 9(2):343–348, 2004.

[25] Zunxian Yang, Xinxin Li, Yuelin Wang, Haifei Bao, and Min Liu. Micro cantilever probe array integrated with piezoresistive sensor. *Microelectronics Journal*, 35(5):479–483, 2004.

[26] Zunxian Yang, Ying Yu, Xinxin Li, and Haifei Bao. Nano-mechanical electro-thermal probe array used for high-density storage based on nems technology. *Microelectronics Reliability*, 46(5-6):805–810, 2006.

[27] K. Yosida. *Functional analysis*. Classics in Mathematics. McGraw–Hill, Berlin, 1995. Reprint of the sixth (1980) edition.

[28] Xiaomei Yu, Dacheng Zhang, Ting Li, Xiaobao Wang, Yong Ruan, and Xianfeng Du. Fabrication and analysis of micromachined cantilever array. *Pan Tao Ti Hsueh Pao/Chinese Journal of Semiconductors*, 24(8):861–865, 2003.

[29] J. Zlatanova, S.M. Lindsay, and S.H. Leuba. Single molecule force spectroscopy in biology using the atomic force microscope. *Progress in Biophysics; Molecular Biology*, 74(1-2):37–61, 2000.