Exponential Stability and Design of Sensor Feedback Amplifiers for Fast Stabilization of Magnetizable Piezoelectric Beam Equations

\IEEEPARstart

Piezoelectric ceramics, including lead-free materials like barium titanate, sodium potassium niobate, and sodium bismuth titanate, are renowned multifunctional smart materials that generate electric displacement in response to mechanical stress [26]. Their small size and high power density make them ideal for various industrial applications, including implantable biomedical devices [5, 25], wearable interfaces with PVDF sensors [7], biocompatible sensors [32], and ultrasound imagers and cleaners [24]. Their fast response, large mechanical force, and fine resolution contribute to their effectiveness [10].

Consider a piezoelectric beam, clamped on one side and free on the other, with sensors for tip velocity and current. he beam is of length $L$ and thickness $h$. Assume that the transverse oscillations of the beam are negligible, so the longitudinal vibrations, in the form of expansion and compression of the center line of the beam, are the only oscillations of note together with the electromagnetic effects. Existing mathematical models often oversimplify intrinsic mechanical or electromagnetic interactions, impacting boundary feedback stabilizability [1, 17, 20, 6, 27, 29].

While electrostatic and quasi-static approaches based on Maxwell’s equations are typically sufficient for non-magnetizable piezoelectric beams, assessing radiated electromagnetic power requires consideration of electromagnetic waves generated by mechanical fields [28]. Therefore, fully dynamic models of piezoelectric beams are essential. Denoting $v(x,t)$ as the longitudinal oscillations and $p(x,t)$ as the total charge accumulated at the electrodes, the equations of motion form following system of partial differential equations [17]

$\displaystyle\begin{array}[]{ll}\begin{array}[]{ll}\begin{bmatrix}\rho&0\\ 0&\mu\end{bmatrix}\begin{bmatrix}v_{tt}\\ p_{tt}\end{bmatrix}-\begin{bmatrix}\alpha&-\gamma\beta\\ -\gamma\beta&\beta\end{bmatrix}\begin{bmatrix}v_{xx}\\ p_{xx}\end{bmatrix}=\begin{bmatrix}0\\ 0\end{bmatrix},&\\ \end{array}\\ \begin{array}[]{ll}(v,p)(0,t)=0,&\\ \begin{bmatrix}\alpha&-\gamma\beta\\ -\gamma\beta&\beta\end{bmatrix}\begin{bmatrix}v_{x}\\ p_{x}\end{bmatrix}(L,t)=\begin{bmatrix}u_{1}(t)\\ u_{2}(t)\end{bmatrix},\quad t\in\mathbb{R}^{+}\end{array}\\ \left[v,p,v_{t},p_{t}\right]\left(x,0\right)=\left[v_{0},p_{0},v_{1},p_{1}\right]\left(x\right),~{}x\in\left[0,L\right]\end{array}$

(7)

where $\rho$, $\alpha$, $\beta$, $\gamma$, and $\mu$ are the mass density per unit volume, the elastic stiffness, the impermeability, the piezoelectric constant, and the magnetic permeability, respectively, and $u_{1}(t)$ and $u_{2}(t)$ are strain and voltage actuators, respectively.

Define $\alpha:=\alpha_{1}+\gamma^{2}\beta>0$ with $\alpha_{1}>0$, and

$\displaystyle\begin{array}[]{l}\zeta^{\mp}=\frac{1}{\sqrt{2}}\sqrt{\frac{\alpha\mu}{\alpha_{1}\beta}+\frac{\rho}{\alpha_{1}}\mp\sqrt{\left(\frac{\alpha\mu}{\alpha_{1}\beta}+\frac{\rho}{\alpha_{1}}\right)^{2}-\frac{4\rho\mu}{\beta\alpha_{1}}}}.\end{array}$

(9)

Note that $\sqrt{\frac{\rho}{\alpha}}$ and $\sqrt{\frac{\beta}{\mu}}$ are non-identical and represent significantly different wave propagation speeds in (7). The natural energy of the solutions is defined as

$\displaystyle\begin{array}[]{l}E(t)=\frac{1}{2}\int^{L}_{0}\left[\rho\left|v_{t}\right|^{2}+\mu\left|p_{t}\right|^{2}+\alpha_{1}|v_{x}|^{2}\right.\\ \hskip 128.0374pt\left.+\beta\left|\gamma v_{x}-p_{x}\right|^{2}\right]dx.\end{array}$

(12)

The exact observability result for the model (7) with $u_{1}(t),~{}u_{2}(t)=0$ and with only one sensor measurement is not possible [17]. However, by employing two sensor measurements, $v_{t}(L,t)$ (tip velocity) and $p_{t}(L,t)$ (total current accumulated at electrodes), a suboptimal observation time is achieved [23]. This result is later refined to include the optimal observation time [21].

When the observed sensor signals $v_{t}(L,t)$ and $p_{t}(L,t)$ are amplified and fed back to (7), a closed-loop system is formed. The amplification range depends on the sensor limits. While using only one sensor feedback can make the system energy dissipative, it is insufficient for exponential stabilization in $\mathcal{H}$. This is because the closed-loop system with this control design imposes strict conditions on stabilization results [17]. Exponential stability is jeopardized for a large class of material parameters and is only achievable for a small subset [18].

Let $\xi_{1},\xi_{2}>0$ be the feedback amplifiers for the two sensor measurements of the closed-loop system $v_{t}(L,t)$ and $p_{t}(L,t)$, respectively. The strain and voltage actuators $u_{1}(t)$ and $u_{2}(t)$ are chosen to be proportional to the sensor feedback amplifiers,

$$\begin{array}[]{ll}\begin{bmatrix}u_{1}(t)\\ u_{2}(t)\end{bmatrix}=-\begin{bmatrix}\xi_{1}&0\\ 0&\xi_{2}\end{bmatrix}\begin{bmatrix}\dot{v}(L,t)\\ \dot{p}(L,t)\end{bmatrix}.\end{array}$$

(15)

The exponential stability of the closed-loop system (7), (15) has been studied in [23], where the proof relies on decomposing the system into conservative and dissipative components. However, this method does not explicitly describe the decay rate or allow optimization of the feedback amplifiers.

The primary objective of this paper is to establish the existence of an exponential decay rate for the system described by (7) and (15) through the meticulous construction of a Lyapunov function. Given the impracticality of traditional spectral analysis due to the strong coupling in the wave system, we adopt a multiplier approach combined with an optimization argument to define a safe range of intervals for each feedback amplifier $\xi_{1}$ and $\xi_{2}$. This ensures that the decay rate provided by the Lyapunov approach is achievable for any type of initial conditions $\mathcal{H}$.

Our methodology provides a maximal decay rate for the system, serving as an upper bound to the optimal decay rate since the actual decay rate could potentially be even better. Existing literature on optimal actuator designs commonly offers different approaches for infinite-dimensional systems [8, 13, 15, 22], and finite-dimensional systems [4, 9]. Our approach is particularly valuable for model reductions by Finite Differences or Finite Elements [21] for (7). Notably, the Lyapunov-based exponential stability proof [2], uniformly as the discretization parameter $h\to 0$ with the recently proposed order-reduced Finite Differences [12], closely resembles the approach employed here.

The decay rate $-\sigma$ in Theorem 2 provides an upper bound for the exponential decay rate of the system. In this section, we provide an optimization process for feedback amplifiers $\xi_{1}$ and $\xi_{2}$ that ensure the maximal value of $\sigma$, given by

$$\sigma_{max}(\delta)=\frac{1}{4\eta L}~{}{\rm achieved~{}at}~{}\delta=\frac{1}{2\eta L}.$$

(29)

Since $\delta$ and $\sigma$ are functions of $\epsilon$ and the sensor feedback amplifiers $\xi_{1}$ and $\xi_{2}$ in (15), the following results provide intervals for the feedback amplifiers and $\epsilon$ that ensures the decay rate (29) is attained.

The amplifier $\xi_{1}$ maximizing $h_{1}$ and the amplifier $\xi_{2}$ minimizing $h_{2}$ are always in the respective intervals, i.e. $\xi_{1}=\frac{\rho}{{2}\eta}\in{(}c_{1}^{-},c_{1}^{+}{)}$ and $\xi_{2}=\frac{\mu}{{2}\eta}\in{(}c_{2}^{-},c_{2}^{+}{)}$. Moreover, as $\epsilon$ approaches to the upper (lower) bound ${(}c_{1}^{-},c_{1}^{+}{)}$ gets smaller (larger) and ${(}c_{2}^{-},c_{2}^{+}{)}$ gets larger (smaller), see Fig. 2. Feedback amplifiers chosen within the intervals ensure that $f_{1}(\xi_{1})\geq\frac{1}{{2}\eta}$ and $f_{2}(\xi_{2})\geq\frac{1}{{2}\eta}$ are satisfied, see Figs. 3 and 4.

In this section, we present numerical simulations primarily aimed at demonstrating the following:

• •

  The exponential stability of the system (7) with (15) is showcased using feedback amplifiers $\xi_{1}\in(c_{1}^{-},c_{1}^{+})$ and $\xi_{2}\in(c_{2}^{-},c_{2}^{+})$, as obtained in Theorem 3.

• •

  Comparison of the reduced model exponential decay rates for a selection of feedback amplifiers.

Moreover, we address a common challenge encountered in standard model reductions, such as Finite Differences (FD) and Finite Elements (FE), which introduce spurious high-frequency vibrational modes to the system. These spurious modes significantly impact the decay rate of the system. In fact, reduced models lack the uniform observability property as the mesh parameter approaches zero, necessitating filtering techniques such as the direct Fourier filtering [3, 11] or the indirect filtering [30].

To mitigate the adverse effects of approximation methods on the decay rate, we consider a recently introduced technique known as Order-Reduced Finite Differences (ORFD) [12]. The ORFD method has been observed to retain the decay rate of the infinite-dimensional system. For rigorous stability analyses of all three approximations, interested readers can refer to [2].

Consider a given natural number $N\in\mathbb{N}$, which denotes the number of nodes in the spatial semi-discretization. Let’s introduce the mesh size as $h:=\frac{1}{N+1}$, and discretize the interval $[0,L]$ as $0=x_{0}<x_{1}<\dots<x_{j}=j*h<\ldots<x_{N}<x_{N+1}=L.$ Let $(v_{j},p_{j})=(v_{j},p_{j})(t)\approx(v,p)(x_{j},t)$ be the approximation of the solution $(v,p)(x,t)$ of (7)-(15) at the point space $x_{j}=j\cdot h$ for any $j=0,1,...,N,N+1$, and let $\vec{v}=[v_{1},v_{2},...,{v_{N+1}}]^{T}$ and $\vec{p}=[p_{1},p_{2},...,{p_{N+1}}]^{T}$.

In the ORFD method, in addition to the standard nodes $\left\{x_{j}\right\}_{j=0}^{N},$ we also consider the in-between middle nodes of each subinterval, denoted by $\left\{x_{j+\frac{1}{2}}:=\frac{x_{j+1}+x_{j}}{2}\right\}_{j=0}^{N}.$ Define the average, $v_{j+\frac{1}{2}}:=\frac{v_{j+1}+u_{j}}{2}$ and difference operators

$\displaystyle\begin{array}[]{ll}\delta_{x}v_{j+\frac{1}{2}}:=\frac{v_{j+1}-v_{j}}{h},&\delta_{x}^{2}v_{j}:=\frac{v_{j+1}-2v_{j}+v_{j-1}}{h^{2}},\end{array}$

(43)

It is worth noting that by considering odd-number derivatives at the in-between nodes within the uniform discretization of $[0,L]$, higher-order approximations are achieved [12].

The ORFD approximation of equations (7) with (15) is

$$\left(\bm{C}_{1}\otimes\bm{M}\right)\begin{bmatrix}{\ddot{\vec{v}}}\\ {\ddot{\vec{p}}}\\ \end{bmatrix}+\left(\bm{C}_{2}\otimes\bm{A}_{h}\right)\begin{bmatrix}{{\vec{v}}}\\ {{\vec{p}}}\\ \end{bmatrix}\\ +\left(\bm{C}_{3}\otimes\mathcal{B}\right)\begin{bmatrix}{\dot{\vec{v}}}\\ {\dot{\vec{p}}}\\ \end{bmatrix}=\vec{0},$$

(44)

where $C_{1}$ and $C_{2}$ are matrices for material parameters, and $C_{3}$ is the matrix for the state feedback amplifiers, defined as

$$\bm{C}_{1}=\begin{bmatrix}\rho&0\\ 0&\mu\end{bmatrix},\quad\bm{C}_{2}=\begin{bmatrix}\alpha&-\gamma\beta\\ -\gamma\beta&\beta\end{bmatrix},\quad\bm{C}_{3}=\begin{bmatrix}\xi_{1}&0\\ 0&\xi_{2}\end{bmatrix},$$

the $(N+1)\times(N+1)$ mass matrix $\bm{M}$ is defined by

$$\bm{M}=\frac{1}{4}\begin{bmatrix}2&1&0&\dots&\dots&\dots&0\\ 1&2&1&0&\dots&\dots&0\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\\ 0&\dots&\dots&0&1&2&1\\ 0&\dots&\dots&\dots&0&1&1\\ \end{bmatrix},$$

the $(N+1)\times(N+1)$ central difference matrix $\bm{A}_{h}$ is

$$\bm{A}_{h}=\frac{1}{h^{2}}\begin{bmatrix}2&-1&0&\dots&\dots&\dots&0\\ -1&2&-1&0&\dots&\dots&0\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\\ 0&\dots&\dots&0&-1&2&-1\\ 0&\dots&\dots&\dots&0&-1&1\\ \end{bmatrix},$$

and the $(N+1)\times(N+1)$ boundary node matrix $\mathcal{B}$ is a zero matrix except for the $(N+1)\times(N+1)$-th entry, which is $\frac{1}{h},$ due to the boundary damping being injected at the last node. Here, $\otimes$ denotes the matrix Kronecker product. It’s worth noting that equation (44) can be rewritten in the first-order form as

$$\frac{d}{dt}\begin{bmatrix}{{\vec{v}}}&{{\vec{p}}}&{\dot{\vec{v}}}&{\dot{\vec{p}}}\end{bmatrix}^{T}=\bm{\mathcal{A}}\begin{bmatrix}{{\vec{v}}}&{{\vec{p}}}&{\dot{\vec{v}}}&{\dot{\vec{p}}}\end{bmatrix}^{T},$$

(45)

$$\bm{\mathcal{A}}=\begin{bmatrix}\bm{0}_{2N+2}&I_{2N+2}\\ -\left(\bm{C}_{1}^{-1}\bm{C}_{2}\right)\otimes\left(\bm{M}^{-1}\bm{A}_{h}\right)&-\left(\bm{C}_{1}^{-1}\bm{C}_{3}\right)\otimes\left(\bm{M}^{-1}\mathcal{B}\right)\end{bmatrix}$$

(46)

The discretized energy corresponding to (45) is

$$\begin{split}E_{h}(t):=&\frac{h}{2}\left\langle(\bm{C}_{1}\otimes\bm{M})\begin{bmatrix}{\dot{\vec{v}}}\\ {\dot{\vec{p}}}\end{bmatrix},\begin{bmatrix}{\dot{\vec{v}}}\\ {\dot{\vec{p}}}\end{bmatrix}\right\rangle\\ &\qquad\qquad+\frac{h}{2}\left\langle(\bm{C}_{1}\otimes\bm{A}_{h})\begin{bmatrix}{{\vec{v}}}\\ {{\vec{p}}}\end{bmatrix},\begin{bmatrix}{{\vec{v}}}\\ {{\vec{p}}}\end{bmatrix}\right\rangle.\end{split}$$

(47)

To show the importance of our analysis for the choices of sensor feedback amplifiers via the optimization process outlined in the previous section, sample numerical simulations are presented for the material constants in Table 1.

Note that $\sigma_{max}(\delta)$ in (29) depends solely on the material parameters. In particular, $\sigma_{\rm max}=\frac{1}{4\eta L}\approx 102.04$ for the material parameters in Table 1. Letting $\epsilon=1$, which is a valid choice for any system according to equations (41), the intervals for optimal sensor feedback amplifiers in (30), corresponding to $\sigma_{\rm max},$ are as follows:

$$\begin{split}(c_{1}^{-},c_{1}^{+})&\approx(7.17\times 10^{5},~{}4.18\times 10^{6}),\\ (c_{2}^{-},c_{2}^{+})&\approx(1.02\times 10^{-4},~{}9.78\times 10^{9}).\end{split}$$

(48)

Here note that both the maximal decay rate $-\sigma_{\rm max}$ and the corresponding intervals for feedback amplifiers in (48) are independent of the choice of initial conditions.

For simulations, we take $N=80$ nodes, and we choose triangular (hat-type) initial conditions that are composed of high and low-frequency vibrational nodes. The simulations are performed for $T=0.1$ sec to show the rapid exponential decay of solutions. Indeed, the solutions of the system (44) with the choices of optimal feedback amplifiers $\xi_{1}=10^{6}\in(c_{1}^{-},c_{1}^{+})$ and $\xi_{2}=10^{9}\in(c_{2}^{-},c_{2}^{+})$, are observed drive all states to zero state rapidly, e.g. see Fig. 5, and the total energy follows the exponentially decay in a higher rate than claimed $\sigma_{\rm max}$, Fig. 7.

Our analysis indicates that even a relatively small perturbation in the values of the feedback amplifiers $\xi_{1}$ can lead to suboptimal performance. Indeed, when $\xi_{1}=10^{4}$ and does not fall within the range $(c_{1}^{-},c_{1}^{+})$, the system still continues to exhibit low and high-frequency longitudinal oscillations for $T=0.1$ seconds, as demonstrated in Figure 6. Furthermore, the energy of the system decays in an unstable manner, with a decay rate significantly larger than $-\sigma_{\rm max}$, as depicted in Figure 7.

In summary, using the Lyapunov approach, we derived the exponential decay rate $``-\sigma"$ for system (7), with the design of state feedback amplifiers $\xi_{1}$ and $\xi_{2}$ as described in (15). This decay rate ensures fast stabilization, though improvements are possible. The rate depends on material parameters and feedback amplifiers. Our numerical results confirm the robustness of this decay rate, with deviations affecting it from the theoretical value. Interactive simulations are available through the Wolfram Demonstrations Project [31].

A full numerical analysis and proof of exponential stability for the semi-discretized model (44), as $h\to 0$, and determining optimal designs of feedback amplifiers remain beyond this paper’s scope and are ongoing research. However, the Lyapunov approach here serves as the basis for results in [2], where the discretized model achieves the same decay rate as the PDE model.

Future work can extend these concepts to multi-layer magnetizable piezoelectric beams, which have numerous practical applications [19]. The design of feedback amplifiers may become more delicate when there are more than two state feedback amplifiers, as highlighted in [19]. There is also potential to explore general hyperbolic PDE systems with multiple boundary dampers. Relevant studies address systems with random inputs [14] and the optimal design and placement of actuators in higher-dimensional systems [16, 15], both of which are valuable for future research.