## 2 Design and Modeling
### Modeling
As shown in Figure 2, it can be found that the interaction forces between excitation magnets and center magnet are attractive force, and the magnetic moment vectors \(\mathbf{\mu}_{i}\) of excitation magnets \(E_{i}\) can be given by  
\[\mathbf{\mu}_{i}=(-1)^{N}\big{(}M_{E}V_{E}\cos\beta_{i}\mathbf{\hat{e}}_{x}+M _{E}V_{E}\sin\beta_{i}\mathbf{\hat{e}}_{y}\big{)} \tag{6}\]  
where \(N=0\) represents the circular parallel and in phase arrangement with the attractive interaction forces between the excitation magnets and center magnet, \(N=1\) represents the circularly parallel andout-of-phase arrangement with the repulsive interaction forces between the excitation magnets and center magnet, and \(N=i-1\) represents the arrangement with alternative parallelity where attractive and repulsive forces occurs alternatively. \(M_{E}\) is the magnitude of the magnetization vector and it can be calculated as \(M_{E}=B_{r}/u_{0}\). \(V_{E}\) is the volume of the magnet \(E_{i}\), and \(\beta_{i}\) can be given as \(\beta_{i}=\beta+2(i-1)\)\(\pi/n\), where \(n\) is the number of the excitation magnets. The magnetic fields produced by the magnet \(C\) upon magnets \(E_{i}\) can be obtained as  
\[\textbf{B}_{i}=-\frac{\mu_{0}}{4\pi}\
abla\frac{\textbf{\mu}_{C}\cdot\textbf{ r}_{i}}{\left\|\textbf{r}_{i}\right\|_{2}^{3}} \tag{7}\]  
where \(\left\|.\right\|_{2}\) and \(\
abla\) denote Euclidean norm and vector gradient operator, respectively. The overall magnetic potential energy between the excitation magnets and center magnet can be given as  
\[U_{in}=-\sum_{i=1}^{n}\textbf{B}_{i}\cdot\textbf{\mu}_{i} \tag{8}\]  
\(F_{ix}\) is the \(x\)-direction interaction force between the center magnet and the excitation magnets. \(F_{ix}\) can be calculates as  
\[F_{ix}=M_{i}\left\{ \begin{array}{l}\frac{3(R+d-R\cos\beta_{i})\cos\beta_{i}+( R\sin\beta_{i}-w)\sin\beta_{i}}{\left[(R+d-R\cos\beta_{i})^{2}+(R\sin\beta_{i}-w) ^{2}\right]^{5/2}}\\ -\frac{5(R+d-R\cos\beta_{i})^{2}[(R+d-R\cos\beta_{i})\cos\beta_{i}+(R\sin\beta _{i}-w)\sin\beta]}{\left[(R+d-R\cos\beta_{i})^{2}+(R\sin\beta_{i}-w)^{2} \right]^{7/2}}\end{array} \right\} \tag{9}\]  
where \(M_{i}=3M_{E}V_{E}M_{C}V_{C}\mu_{0}(-1)^{j}/(4\pi)\), \(j=1\) and \(2\) denote the repulsive and attractive interaction force between the center magnet and the excitation magnet, respectively. The magnetic force applied to the buckled beam can be regarded as a concentrate force \(F_{in}\), given as  
\[F_{in}=\sum_{i=1}^{n}F_{ix} \tag{10}\]  
The potential energy of the magnetic field between the tip magnet and the fixed magnet on the frame can be written as  
Figure 2: Schematic of dynamic analysis of the MBBP-REH with the parallelity arrangement of the excitation magnets. \[U_{m}=\frac{3u_{0}}{2\pi}\Bigg{(}\frac{\upmu_{T}\upmu_{F}}{\|\upmu_{TF}\|_{2}^{3}} \Bigg{)} \tag{11}\]  
where \(\upmu_{T}\) and \(\upmu_{F}\) are the magnetic moment vectors of the tip magnet and the fixed magnet, and \(\uptau_{TF}\) is the center distance from \(\upmu_{T}\) to \(\upmu_{F}\). Thus, the Lagrange function of the system is defined as the difference between the kinetic energy and the potential energy  
\[L=T-U-U_{m} \tag{12}\]  
Taking into account of the primary mechanical damping, electrical damping, and input concentrate force, the virtual work done by the non-conservative force can be given by  
\[\delta W=-c_{1}\int_{0}^{l_{1}}\dot{w}\delta wdx-c_{2}\dot{u}\delta u+Q\delta V +F_{in}\delta w \tag{13}\]  
where \(c_{1}\) and \(c_{2}\) are the damping coefficients to approximate the energy loss of viscous damping and frictions, and \(Q\) is the electric charge output of the piezoelectric patches. The time rate change of electric charge \(Q\) is the electric current transiting through the resistive loading \(R_{L}\), and \(F_{in}\) is the interaction force between the center magnet and excitation magnets, which is regarded as the input force for the buckled beam system.