Structure and controller design of a piezo-driven orientation stage for space antenna pointing

Mechanical Systems and Signal Processing 138 (2020) 106525

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Structure and controller design of a piezo-driven orientation stage for space antenna pointing Shubao Shao a, Yan Shao a, Siyang Song a, Minglong Xu a,⇑, Xiaofei Ma b a b

State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, PR China China Academy of Space Technology(Xi’an), Xi’an 710000, PR China

a r t i c l e

i n f o

Article history: Received 28 May 2019 Received in revised form 26 August 2019 Accepted 12 November 2019

Keywords: Piezo-driven Orientation stage Support vector regression Adaptive control

a b s t r a c t Space antenna pointing is key to establishing inter-satellite and earth-satellite communication links. Compared with the traditional electromagnetic driven mechanism for antenna pointing, the piezo-driven orientation stage (POS) is more compact in structure and has higher pointing accuracy. However, the main difficulties of designing a POS include the compromise between workspace and structural stiffness, as well as devising a controller that can compensate the hysteresis behavior of the piezo-driven unit and coupled motion between axes. This paper presents a 2-DOF POS in which the rhombic mechanism and 2D flexure hinge are used to enlarge the workspace and decouple the motion between axes. To improve the mechanical modeling efficiency, support vector regression (SVR) and finite element analysis (FEA) are employed in combination. Subsequently, the trade-off between the workspace and structural stiffness is achieved by implementing the structure optimization design. To compensate the hysteresis and remaining coupling motion between axes, multiple input multiple output (MIMO) adaptive control is used to improve the control accuracy, where the real-time nonlinear model of the controlled plant with hysteresis property is approached by the controlled auto regression and moving average (CARMA) model. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction In space applications, a huge boost in the capacity of communication satellites has been demanded to accommodate the expansion of human social and economic development. Satellites with multi-antenna and several small spot beams over the coverage service area tend to have much greater communication capacity due to the parallel communication of the antenna group, high antenna gain, and the increased frequency reuse. However, the small spot beams require ultra-high pointing accuracy [1,2] and independent pointing control of each antenna is required. Consequently, the 2-DOF orientation stage (the crucial mechanism for the antenna to achieve two-dimensional beam scanning) is subjected to tighter pointing accuracy and a more compact structure. Unfortunately, the conventional orientation stages that employ electromagnetic motors suffer from insufficient pointing accuracy and are bulky in structure. Hence, the use of electromagnetic orientation stages hinders the potential performance improvement and structural miniaturization of the antenna system for future space applications.

⇑ Corresponding author. E-mail address: [email protected] (M. Xu). https://doi.org/10.1016/j.ymssp.2019.106525 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

Compared with the electromagnetic pointing mechanism, the piezo-driven orientation stages (POS), which are driven by piezoelectric actuators (PEAs), can achieve much higher pointing accuracy and miniaturized structure. This is due to the advantages of PEAs which include extreme fine resolution, large force generation, and small volume, along with the merits of flexure-based mechanisms that include lack of backlash, wear-free property and compactness [29]. For antenna pointing applications, the wide range orientation is first generally achieved by the attitude control of satellite or the coarse adjustment mechanisms, then the fine orientation is accomplished by the fine pointing stages whose motion range are set as milliradian scale so as to acquire much higher pointing accuracy. In addition, the large structural stiffness of the fine pointing mechanisms are also needed to improve the structure’s fundamental frequency, thus maintaining the sufficient fast response capability. Although the amplifiers used to magnify the limited stroke output of PEA (and therefore expanding the workspace of the positioning system) can be obtained by appropriately designed flexible mechanisms, the maximization of the displacement amplification ratio, together with the stiffness design for proper dynamic response, has proved to be challenging. Some researchers focus on developing flexible joints that can achieve large displacement for quasi-static motions [30,31], and the flexible constant velocity universal joints that can accommodate high misalignment angles between input and output rotational axes to transfer rotary motion or energy between two angled shafts [32,33]. In above applications, the large stroke, high precision and tolerance ability are more critical concern compare to the dynamic response of compliant mechanism. However, the design of orientation stage for antenna pointing applications requires the trad-off between the flexibility and dynamic response of the flexure-based mechanisms. In 2011, Moon et al. [3] developed a 6-DOF micro-positioning stage featuring low parasitic motions achieved by properly arranging the flexure hinges and PEAs. The rotational range and the natural frequency were reported to be approximately 0.1mrad and 148.6 Hz, respectively. In 2018, Cai et. Al [4,5] proposed a 6-DOF micro-positioning stage with fundamental frequency of 528 Hz and rotational range of 0.1mrad. To improve the workspaces of the piezo-driven positioning stages, the flexible amplifiers of different types (such as level-type mechanism [6], bridge-type mechanism [7,8], and flexure rhombic mechanism [9,10]) were employed. Using the aforementioned amplifiers (with amplification ratios 5.5–27.5), the in-plane motion stages [11–13] and out-of-plane motion stages [6,14] with significant enlarged workspaces and compact constructions were proposed. Nevertheless, their resonant frequencies were dramatically reduced. It is clear that a trade-off between the motion range and resonant frequency is required to obtain the optimal design. Numerous studies document efforts to establish optimization models where the accurate theoretical derivation of the stiffness and the displacement amplification ratio of the positioning system composed of different flexure-hinges are required [11,12,14–16]. Unfortunately, this means that the methods are difficult to implement and are inefficient. In this paper, the stiffness and displacement amplification ratio of the flexible mechanism are modeled by employing the support vector regression (SVR) based on training data calculated from finite element analysis (FEA). Consequently, the modelling efficiency is improved while maintaining modelling precision. The high-precision control of the multi-DOF piezo-stages should deal with the coupling motions between axes and the hysteresis nonlinearity inherent in PEA. One method to reduce the coupling motions between axes is based on the decoupling structures [17,18], which present low coupling ratios between axes when driven by simple open-loop control only. An alternative method to address coupling motions between axes relies on complex control methods such as decoupling control [19] and independent modal space control [20], which can foreseeingly compensate the coupling disturbance from other axis. For the compensation of hysteresis nonlinearity of PEA, the feedback control may be considered with the inverse hysteresis model as the feedforward compensation [21]. The inverse hysteresis model can be derived from the Preisach model [22], Prandtl-Ishlinskii model [23], Duhem model [24], Bouc-Wen model [25], and geometry based hysteresis model [26]. However, these hysteresis models have issues such as complexity, parameters which are difficult to find, and the inverse model being difficult to obtain. In this paper, the decoupling structure is first designed to reduce the coupling ratio between axes, and then multiple input multiple output (MIMO) adaptive control is used to improve the control accuracy. The real-time nonlinear model of the controlled plant with hysteresis property is approached by the auto regression and moving average (CARMA) model. Decoupling control is intended to improve the individual accuracy of each DOF, and therefore may result in the loss of overall accuracy. In contrast, the MIMO adaptive control makes the compromise between two DOF to optimize the synthesized accuracy. The combination of the decoupling structure and MIMO adaptive control significantly improves the pointing accuracy.

2. Antenna pointing driven by POS 2.1. Antenna pointing A conceptual schematic of antenna pointing is shown in Fig. 1(a). The POS mounted on the bottom of the antenna is used to accomplish the task of precise pointing of the antenna to the determined target to track and maintain the communication link between the satellite and ground station (or other satellite). Fig. 1(b) shows the driving principle of POS. A 2-DOF flexure hinge set in the middle of the stage allows only rotations about the x and y axes. When the force coupled in the yz plane (composed of F and -F) is applied on the stage, the rotational angle, Dhx , about the x-axis can be produced due to the moment equilibrium. Similarly, the rotational angle, Dhy , about the y-axis can be produced by applying the force coupled in the xz plane.

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

3

Fig. 1. Conceptual schematic of antenna pointing: (a) schematic and (b) driving principle.

2.2. Structure of POS The POS is intended to have rotational angles as large as possible and sufficient stiffness to meet the requirements for the desired orientation ranges, while guaranteeing the dynamic property for fast response. The mechanical design of the proposed POS is shown in Fig. 2(a)–(c) in which the CAD model is depicted from the side view, the oblique view, and top view, respectively. The POS is composed of three main parts: the stage, four identical flexible mechanisms (FM) with mounted PZT stacks, and the baseboard (see Fig. 2(a)). The four flexible mechanisms sit on the baseboard and support the stage. Two flexible mechanisms (FM-1 and FM-2, as shown in Fig. 2(b) and (c)) are arranged parallel to each other along the x-axis. The other two (FM-3 and FM-4, as shown in Fig. 2(b) and (c)) are arranged parallel along the y-axis to produce the in-plane force couples. Each flexible mechanism is composed of a rhombic amplifier (RA) used to magnify the stroke of the PZT stack and a 2D-flexure hinge (FH) that can transmit 2-DOF motion with low coupling ratio (see Fig. 2(d)). Note, the overall structure of

Fig. 2. CAD model of the proposed POS (a) side view, (b) oblique view, (c) top view, and (d) FM composed of a FH and a RA.

4

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

the POS is integral machined. Except for the installation of the PZT stack, no assembly is required and as such assembly error is avoided. 3. Structural design of POS 3.1. Analytic model One key consideration for the design of POS is the resonant frequency that limits the response time of POS. The dynamic model of the proposed POS for the frequency analysis can be depicted as in Fig. 3 where the four flexible mechanisms are considered as four individual flexible beams. The stage and the baseboard are treated as the rigid body and the clamped end, respectively. The four flexible beams have identical flexibility but are arranged in different orientations. The linear rela
T
T tion between the force F ¼ F x ; F y ; F z ; M x ; M y ; M z and the deformation X ¼ x; y; z; hx ; hy ; hz of the flexible beams (S1–S4) is established as KFM X ¼ F(see Fig. 3), where the stiffness matrix KFM is [27]:

2

KFM

kFM1

6 0 6 6 6 0 ¼6 6 0 6 6 4 kFM10 0

0

0

0

kFM7

kFM2 0

0 kFM3

kFM8 0

0 0

kFM9

0

kFM4

0

0

0

0

kFM5

0

0

0

0

0

3

7 7 7 7 7 0 7 7 7 0 5

0 0

ð1Þ

kFM6

In Eq. (1), kFM1  kFM6 are the stiffness of six degree of freedom (6-DOFs), and kFM1  kFM6 represent the coupled stiffness  T between axes. If we denote the 6-DOFs for the rigid body (stage) as q ¼ x; y; z; hx ; hy ; hz . The equation of motion (neglecting damping) for the rigid body is:

€ þ Kq ¼ Q Mq

ð2Þ


 where the mass matrix M has the form of M ¼ diag mx ; my ; mz ; J x ; J y ; J z . mx ; my and mz denote the mass of the rigid body. J x ; J y
 and J z represent the rotational inertia about the x-, y- and z-axes. K ¼ diag kx ; ky ; kz ; khx ; khy ; khz is the stiffness matrix whose   2 2 elements can be calculated as kx ¼ ky ¼ 2ðkFM1 þ kFM2 Þ, kZ ¼ 4kFM3 , khx ¼ khy ¼ 2 kFM4 þ kFM3 l , and khz ¼ 4kFM2 l . Q ¼ ½ Q 1 Q 2 ::: Q 6 T represent the force vector on the stage. According to vibration theory, the resonant frequencies of the mechanism can be obtained by solving the characteristic equation of [34]:

   1  kI  M K ¼ 0

ð3Þ

where vector k represents the eigenvector of Eq. (3), and I is the identity matrix. The natural frequencies of the mechanism are

fi ¼

1 pffiffiffiffi ki ði ¼ 1; 2; :::; 6Þ 2p

ð4Þ

where ki are the elements of the vector k. In addition to the fundamental frequency, the tilt ranges are also important to be considered and are expected to be maximal. The static model for determining the rotational range is shown in Fig. 4. Only the

Fig. 3. Dynamic model of the proposed POS and geometric dimensions of the FM.

5

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

Fig. 4. Static analysis model of (a) FM, and (b) POS.

model for the rotation about the x-axis is depicted. The one for the y-axis shares the same model due to the rotational symmetric structure of the POS. Based on Fig. 4(a), the static model of FM can be represented as: 





KFM X ¼ F

ð5Þ





T Dv  and F ¼ ½ DF h DF v  are the displacement vector and force vector, respectively. These are related by " #   the stiffness matrix K ¼ k11 k12 . The electromechanical coupling equation of PZT implies that the force DF generated T

where X ¼ ½ Dh

FM





k21 by PZT stack can be given as:

h

k22

DF h ¼ kp ðd33 V  DhÞ

ð6Þ

where kp ,d33 and V are the stiffness, equivalent piezoelectric constant, and the applied driven voltage of the PZT stack, respectively. In addition, according to the static analysis model as depicted in Fig. 4(b), the static torque balance equation for the stage can be derived as:

2DF v l ¼ kh Dhx

ð7Þ

where l and Dhx are the length of arm of force and the tilt angle of the stage, respectively. kh is the rotational stiffness of pivot O and can be derived as kh ¼ 2kFM4 . Combing Eqs. (4) to (6) and noting that Dv ¼ lDhx due to the small rotational angle, the tilt angle of the stage can then be derived as:

N V M

Dhx ¼ 

where N ¼

ð8Þ

k 12 kp d33 l 

k 11 kp



k

, M ¼ k2h þ  12

kp

k 11 kp



 k22 . Eqs. (3) and (7) reveal that the determination of natural frequencies f i and tilt Dh 

relies on the mass matrix M, the stiffness matrix K, and KFM . For the thin square-plate stage with the edge length of a,
 and the mass of m, the corresponding mass matrix M can be readily derived as M ¼ diag m; m; m; ma2 =12; ma2 =12; ma2 =6 via the theory of rigid-body dynamics. In matrix M, the last three formulas represent the rotational inertias about three axes. 

However, the accurate theoretical derivation of stiffness matrix KFM and KFM is a challenging and inefficient method due to the relatively complex shape of FM. In this paper, the SVR based on training data calculated from FEA is used to determine the required stiffness matrix. Compared with other regression methods (such as the least squares and response surface methods), SVR has the distinct advantage of achieving multivariate regression with higher accuracy based on only small-scale training data. SVR is a powerful machine learning technology based on statistical learning theory. The training data derived from FEA can approximate well to the real values with appropriate finite element grid density. Hence, both the efficiency and accuracy of computation can be improved using the SVR and FEA based determination of the stiffness matrix.

6

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

3.2. Brief introduction of SVR SVR can achieve regression from multiple inputs to a single output, which takes the form [28]:

yðxÞ ¼ xT uðxÞ þ b

ð9Þ

where the nonlinear function uðxÞ maps the input space to the high-dimensional space. fxk ; yk gðk ¼ 1; 2; :::N Þ is given as the training set, where N is the sample size. x and b are parameters, which can be determined by solving the following optimization problem:

min J ðx; eÞ ¼ 12 xT x þ C x;b;n

N P k¼1

n2k

ð10Þ

s:t:yk ¼ x uðxk Þ þ b þ nk T

N P k¼1

In the objective function in Eq. (10), xT x represents the belief bound that influences the predictive power of the model. n2k is the experimental risk related to the training error. C is the regularization parameter, which balances the predictive

power and training error. The Lagrangian function of the optimization problem in Eq. (10) can be expressed as:

Lðx; b; n; aÞ ¼

N N X X   1 T x xþC n2k  ak xT uðxk Þ þ b þ nk  yk 2 k¼1 k¼1

ð11Þ

where ak are the Lagrangian multipliers. The optimal solutions of Eq. (11) meet the following conditions:

8 N P > @L > > ¼0)x¼ ak uðxk Þ > @x > > k¼1 > > @L > < ¼ 0 ) ak ¼ Cnk @n

ð12Þ

k

N > P > @L > > ¼0) ak ¼ 0 > @b > > k¼1 > > : @L ¼ 0 ) xT uðxk Þ þ b þ nk  yk ¼ 0 @x

Eliminating x and n leads to:



0

eT

  b

e

X þ IN =C

a

¼

 0 Y

ð13Þ

where e ¼ ½1; 1; :::; 1T , a ¼ ½a1 ; a2 ; :::; aN T , Y ¼ ½y1 ; y2 ; :::; yN T , IN is the N  N identity matrix, and the elements of matrix X are
 represented as Xij ¼ K xi ; xj . K ðÞ is the radial basis kernel function that takes the form:

 
 K xi ; xj ¼ exp k xi  xj k2 =r2

ð14Þ

where r is the kernel width parameter. After calculating b and a from Eq. (13), the SVR regression model can be obtained as:

yðxÞ ¼

N X

ak K ðxk ; xÞ þ b

ð15Þ

k¼1

3.3. Determination of stiffness based on SVR and FEA The determination of stiffness is necessary to obtain the regressive relation between the stiffness and all dimensional parameters of FM. However, the PZT stack with dimensions 10 mm  10 mm  18 mm is pre-selected. Consequently, several dimensional parameters are fixed to fit the PZT stack and maintain the basic shape of the structure. In this study, the fixed parameters are t 1 ¼ 3mm, lp ¼ 17:98mm, br ¼ 10:5mm, bf ¼ 5mm, l1 ¼ 1mm, and l2 ¼ 2mm. The leg’s thickness tr and the angle a of RA, the length lf , and the thickness t f of FH are the variable parameters. The ranges of these variable parameters are listed in Table 1. Then, regressive relation between the stiffness and these variable parameters can be represented as:

Table 1 Variable dimensional parameters of FM. Parameters

tr (mm)

ar (°)

lf (mm)

tf (mm)

Ranges Samples

[0.5, 3] 0.5, 1.125, 1.75, 2.375, 3

[0, 45] 0, 11.25, 22.5, 33.75, 45

[0.5, 3] 0.5, 1.125, 1.75, 2.375, 3

[0.5, 3] 0.5, 1.125, 1.75, 2.375, 3

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525


 kFMi ¼ f i t r ; a; lf ; tf ði ¼ 1; 2; :::; 10Þ

7

ð16Þ

where kFMi ði ¼ 1; 2; :::; 10Þ are the elements of the stiffness matrix K FM expressed according to Eq. (1). f i ðÞ ði ¼ 1; 2; :::; 10Þ are the maps from multi-input to single output. The SVR-based estimation of Eq. (8) requires a training set with N elements:

D¼

n o xm ; ym ðm ¼ 1; 2; :::; NÞ

ð17Þ


   where xm ¼ tr ; ar ; lf ; t f m are the vectors of variable parameters, ym are the stiffnesses corresponding to xm . For the determination of the stiffness matrix K FM , if each variable parameter has 5 samples (as listed in Table 1) and each sample of each variable parameter combines with all other samples of the rest of the variable parameters to form the vector
  xm ¼ tr ; ar ; lf ; t f m ðm ¼ 1; 2; :::; nÞ, a large sized dataset is required (n ¼ 54 ¼ 625). To reduce the scale of the training set, FM is divided into two flexure parts (FH and RA) whose stiffness matrices are determined using SVR separately. Next, the stiffness matrix of FM is derived based on the mechanical theory. Fig. 5 shows the schematic diagram for determining the stiffness matrix of FM where the two parts (FH and RA) of FM are also considered as the flexure supporters with a similar linear relation between the force and deformation expressed as:

KRA X ¼ F

ð18Þ

KFH X ¼ F


ð19Þ T




T

where X ¼ x; y; z; hx ; hy ; hz and F ¼ F x ; F y ; F z ; M x ; M y ; M z are the displacement and force vectors, respectively. KFH and KRA are the stiffness matrices for FH and RA, and have the same form as KFM . If we represent the elements of flexibility matrices C FH and C RA (flexibility matrices C FH and C RA are the inverse of stiffness matrices KFH and KRA , respectively) as

Fig. 5. Flexible mechanism is divided into two cascaded parts with their individual 6-DOF stiffnesses.

8

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

cFHi ði ¼ 1; 2; :::; 10Þ and cRAi ði ¼ 1; 2; :::; 10Þ, respectively. Consequently, the elements cFMi ði ¼ 1; 2; :::; 10Þ of flexibility matrix CFM (flexibility matrices CFM is the inverse of stiffness matrices KFM ) can be derived as:

8 > < cFMðiÞ ¼ cFHðiÞ þ cRAðiÞ þ cRAðiþ6Þ  lFH ði ¼ 1; 2Þ cFMðiÞ ¼ cFHðiÞ þ cRAðiÞ ði ¼ 3; 4; :::; 8Þ > : cFMðiÞ ¼ cFHðiÞ þ cRAðiÞ þ cRAði5Þ  lFH ði ¼ 9; 10Þ

ð20Þ

where lFH is the total length of FH, the elements cFHi ði ¼ 1; 2; :::; 10Þ and cRAi ði ¼ 1; 2; :::; 10Þ are determined by SVR-based estimation and expressed as:

cFHðiÞ ¼ f FHðiÞ ðtr ; ar Þ ði ¼ 1; 2; :::; 10Þ

ð21Þ


 cRAðiÞ ¼ f RAðiÞ lf ; tf ði ¼ 1; 2; :::; 10Þ

ð22Þ 

According to Eqs. (21) and (22), the vector, xm , for determining the flexibility elements of each part has only 2 variable
   parameters (i.e., xm ¼ ðt r ; ar Þ or xm ¼ lf ; tf ). With the same 5 samples of each variable parameter (as listed in Table 1) 

and the same method of combining parameters mentioned previously to form vector xm , the size of training dataset for each part is n ¼ 52 ¼ 25. Hence, the total number of the training data for two parts is only 2n ¼ 50. In Eqs. (18) and (19), each column of flexibility matrices KFH and KRA represents the force and moment required for each of the displacements taken separately. The elements of the stiffness matrix can be found from the force vector induced by applying each displacement separately. For example, consider a displacement vector X ¼ ð1; 0; 0; 0; 0; 0ÞT is applied on the
T loading surface of FH or RA. Consequently, the force vector F ¼ F x ; F y ; F z ; M x ; M y ; M z is induced. The first column elements
T of KFH or KRA are found to be 1=F x ; 1=F y ; 1=F z ; 1=M x ; 1=M y ; 1=M z . Similarly, the remaining column elements of the stiffness matrix can be found. The above procedure is implemented using FEA simulation during which the variable dimensional parameters are tuned one by one with samples (as listed in Table 1), the FEA models of FH and RA are established by employing the solid45 element and the materiel constant of steel, the solution of FEA takes the ramped loading of single step. 

Finally, the training sets with size of n = 25 for each part are obtained. Note, the stiffness matrix KFM of Eq. (5) is obtained  1

by calculating the inverses of the flexibility matrix KFM determined from FEA due to the convenience of force loading and data processing. To verify the SVR-based modeling method, nine tests with different dimensional parameters were implemented. Fig. 6 shows nine test data sets for validation of stiffness matrix KRA of RA where the stiffnesses along x and y-axes (kRA11 and kRA22 ) were determined from both SVR and FEA. Note, the numbers1-9 on the horizontal axis correspond to the nine differ

ent dimensional parameter sets xm ¼ ðtr ; ar Þ, as indicted by the red arrows, in addition, the selected hyper-parameters p,e and c of SVR are also shown. It can be seen that the maximum full scale errors between SVR and FEA-based calculations are approximately 2.5% for kRA11 and 5% for kRA22 . Fig. 7 shows nine test data sets for validation of stiffness matrix KFH of FH where the stiffnesses along x and y-axes (kFH11 and kFH22 ) were determined from both SVR and FEA. Note, the numbers1-9 on the horizontal axis correspond to the nine
  different dimensional parameter sets xm ¼ t f ; lf , as indicted by the red arrows, in addition, the selected hyper-

Fig. 6. Stiffnesses determined from SVR and FEA for KRA

9

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

Fig. 7. Stiffnesses determined from SVR and FEA for KFH .

parameters p,e and c of SVR are also shown. The maximum full scale errors between SVR and FEA-based calculations of kFH11 and kFH22 are found to be approximately 3.5% and 7.6%, respectively.  1



Fig. 8 shows nine test data sets for validation of flexibility matrix KFM where the flexibility along x-axis, c FM11 , and the 

coupling flexibility from x-axis to y-axis, c FM12 , were determined from both SVR and FEA. Note, the numbers1-9 on the hor

izontal axis correspond to the nine different dimensional parameter sets xm ¼ ðt r ; ar Þ, as indicted by the red arrows, in addition, the selected hyper-parameters p,e and c of SVR are also shown. The maximum full scale errors between SVR and FEA



based calculations of c FM11 and c FM12 are found to be approximately 3.8% and 6.2%, respectively. Compared with the direct determination of KFM through SVR-based method, the above method can reduce the size of required training set from N = 625 to N = 50, and the accuracy is also maintained. The reason why is that the complex structure of FM with four geometrical variables is divided into two flexure parts with simpler structure and only two geometrical variables. 3.4. Structure optimization design For the space antenna pointing application, large orientation range with fast dynamic response time is required. Consequently, the optimization model of POS for the structure design can be written as:


 max J t r ; ar ; tf ; lr ¼ Dh s:t:xn;1 P 450Hz

ð23Þ

 1

Fig. 8. Flexibilities determined from SVR and FEA for KFM

10

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525 



0:5 mm 6 t r 6 2 mm; 0 6 ar 6 45

0:5 mm 6 t f 6 2 mm; 1 mm 6 lr 6 3 mm where Dh and xn;1 represent the tilt range and fundamental frequency of POS. Additionally, the design parameters are tr ; ar ; t f and lf . The determination of tilt range Dh and fundamental frequency xn;1 relies on the SVR-based calculation of the related stiffnesses. Hence, it is difficult to determine the optimal design parameters manually. Herein, the particle swarm optimization (PSO) algorithm (an iteration-based evolutionary computation technique that simulates the phenomenon of bird predation [10]) is used to determine the best solution of the optimization problem expressed in Eq. (14). The optimal parameters found by PSO are listed in Table 2. The corresponding tilt range and fundamental frequency of POS are ± 1.3mrad and 496 Hz, respectively. The first three modes of POS determined by FEA are shown in Fig. 9 where the first mode is the rotation about z-axis with resonant frequency of 472 Hz. The other two modes are the translation along the x-axis with frequency of 770 Hz and the translation along the y-axis with frequency of 771 Hz. The second and third frequencies are close to each other due to the rotationally symmetric construction of POS. With the antenna mounting on the POS, the equivalent inertia mass is added by the antenna, hence the modes of vibration and resonant frequency change dramatically. Fig. 10 shows first three modes of the assemble structure, mounting an antenna with mass of 970 g on POS. We can see that the first mode is the rotation about z-axis with resonant frequency of 35 Hz, the other two modes are the rotations about x-axis and y-axis with frequency of about 87 Hz. Although the significant changes in dynamic property appear when the antenna is taken into account, the antenna only contributes to the equivalent inertia mass, so it’s reasonable to improve the fast response ability of antenna pointing by improving the stiffness of POS itself with sufficient orientation ranges. Table 2 Variable dimensional parameters of FM. Parameters

tr (mm)

ar (°)

lf (mm)

t f (mm)

values

1.16

10.36

2.56

0.5

Fig. 9. First three mode shapes of POS

Fig. 10. First three mode shapes of the assembled structure (Antenna + POS)

11

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

4. Design of control schema 4.1. Dynamic model The differential equation for dynamic motion of the POS can be written as:

M€h þ Ch_ þ Kh ¼ f 

T

ð24Þ 

T

where h ¼ hx ; hy and f ¼ f x ; f y are the generalized displacement and force vectors. M, C, and K are 2  2 matrices that represent mass, damping, and stiffness, respectively. Accordingly, the dynamic motion of POS can be described as a system with an input vector and output vector that are related by the transfer function matrix:





h x ð sÞ

¼

h y ð sÞ

Gxx ðsÞ

Gxy ðsÞ

Gyx ðsÞ Gyy ðsÞ

"

f x ð sÞ

# ð25Þ

f y ð sÞ

 T  T where hx ðsÞ; hy ðsÞ ¼ hðsÞ and f x ðsÞ; f y ðsÞ ¼ F ðsÞ are the Laplace transforms of output tilt and input force vectors, respec

 Gxx ðsÞ Gxy ðsÞ ¼ GðsÞ is the transfer function matrix. The theoretical calculation of the parameters of the transfer tively. Gyx ðsÞ Gyy ðsÞ function matrix usually produces large deviation from the actual model. Herein, system identification is employed. In the control of POS, the control signals generated by the controller are amplified by the power amplifier, and then the driving 

voltages of the PZT stacks are produced. Hence, the transfer function matrix GðsÞ that relates the control signal V ðsÞ and output tilt hðsÞ is needed for the controller design: 

hðsÞ ¼ GðsÞV ðsÞ 

where GðsÞ ¼

" Gxx ðsÞ 

ð26Þ 

Gxy ðsÞ 

#

  T and V ðsÞ ¼ ux ðsÞ; uy ðsÞ . Note, the transfer function matrix GðsÞ linearly maps to the transfer

Gyx ðsÞ Gyy ðsÞ function matrix GðsÞ, this is because the force f generated by the PZT stacks can be derived from the driving voltage V accord

 ing to the relation f ¼ ke kp d33 = ke þ kp V, where ke is the stiffness of the structure that PZT stack is against. Using a sine wave control signal with fixed amplitude and linear scanning frequencies from 10 Hz to 1KHz, the tilts of POS were measured by two laser displacement sensors (configured for differential measurement). Then, using the system identification 

tool in MATLAB, the transfer functions of GðsÞ were calculated as:

8 4:788109 > > Gxx ðsÞ ¼ s2 þ71:38sþ1:9910 7 > > > >  > 7 > 6:83810 < G ðsÞ ¼ xy

s2 þ79:78sþ1:94107

ð27Þ

 > 5:681107 > Gyx ðsÞ ¼ s2 þ88:72sþ1:9710 > 7 > > > > > : G ðsÞ ¼ 4:585109 yy s2 þ70:67sþ1:951107

4.2. Controller design 

Despite using 2D flexure hinges, the remaining coupling effect between axes were found in experiments (i.e., Gxy ðsÞ–0 

and Gyx ðsÞ–0). Additionally, the inherent hysteresis property of the PZT stack would introduce errors to the dynamic model

Fig. 11. Schematic diagram of MVSTC.

12

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

of POS. Accordingly, accurate two-dimensional tracking pointing of POS should address the above two problems. Herein, multivariate self-tuning control (MVSTC) is designed to improve the tracking accuracy of POS. Fig. 11 shows the schematic h  iT  h  iT T diagram of the closed-loop control system in which hx ; hy , ux ; uy , and hx ; hy are the desired pointing angles, the control signal, and the actual pointing angles, respectively. In MVSTC, the controller of the plant adaptively adjusts its paramh  iT  T eters online according to the adaptive law whose updates to the rule are driven by the sampling values of hx ; hy , ux ; uy ,  T and hx ; hy at each sampling instance (such that the dynamic modeling errors caused by the hysteresis behavior of PZT stack would be compensated). In addition, MVSTC uses the composite error of two tilt angles to calculate its adaptive law in order that the coupling effect between axes is compensated. MVSTC requires the difference-equation of the plant, the discrete forms of Eq. (27) can be obtained as Eq. (28) through the relation z ¼ eT s S , where T s ¼ 0:0001s is the sample interval.

8
 23:49z1 þ23:43z2 > Gxx z1 ¼ 11:798z > 1 þ0:991z2 > > > 
>  > < Gxy z1 ¼ 0:3355z11þ0:3347z2 2 11:801z

þ0:992z

ð28Þ


 0:2766z1 þ0:2778z2 > > Gyx z1 ¼ 11:798z > 1 þ0:991z2 > > > >
 : 1 22:5z1 þ22:45z2 Gyy z ¼ 11:802z1 þ0:993z2

The difference-equation of the dynamic model for POS can be expressed as:




 A z1 YðkÞ ¼ z1 B z1 UðkÞ þ C z1 eðkÞ

ð29Þ

8




1  1:983 0 3:599 0 > 1 > A z þ z2 ¼ I þ z > < 0 1:984 0 3:598




 > 23:43 0:3355 23:49 0:3347 1 > > þ z : B z1 ¼ 0:2766 22:50 0:2778 22:45

ð30Þ


 where z1 represents the delay (i.e., z1 UðkÞ ¼ Uðk  1Þ), eðkÞ is the noise with its coefficient matrix C z1 ¼ I,
T
T YðkÞ ¼ hx ðkÞ; hy ðkÞ , and UðkÞ ¼ ux ðkÞ; uy ðkÞ represent the two-dimensional vectors of pointing angle and control input,

 respectively. A z1 and B z1 are the corresponding coefficient matrices that can be derived from Eq. (28) as:

The controller is designed to derive the best control input to ensure that the actual output is close enough to the reference trajectory in each sampling step. The cost function J (reflecting the control objective) is defined as follows:

n o T T J ¼ E ½Yðk þ 1Þ  Yr ðkÞ ½Yðk þ 1Þ  Yr ðkÞ þ ½KUðkÞ ½KUðkÞ

ð31Þ

where Efg is the mathematical expectation operation, Yðk þ dÞ and Yr ðkÞ are the actual output and the given reference signal, respectively. K ¼ lI is the weight coefficient matrix of the input vector UðkÞ. The first term is linked to the objective of minimizing the synthesized error between the actual output and the reference signal. The second term reflects the consideration given to the size of UðkÞ when the objective function J is made to be as small as possible.


 is expressed as A z1 ¼ I  z1 G z1 , in which To obtain the predicted output Yðk þ 1Þ, A z1




 3:599 0 1:983 0 þ z1 . Multiplying both sides of Eq. (29) by z1 gives: G z1 ¼ 0 1:984 0 3:598


 I  z1 G Yðk þ 1Þ ¼ BUðkÞ þ Ceðk þ 1Þ

ð32Þ

Yðk þ 1Þ ¼ GYðkÞ þ BUðkÞ þ Ceðk þ 1Þ

ð33Þ




 where G ¼ G z1 , B ¼ B z1 , and C ¼ C z1 . The predicted output Yðk þ 1Þ can then be derived as:

Substituting Eq. (33) into Eq. (31) and noting that YðkÞ and UðkÞ are irrelevant to eðk þ 1Þ, the cost function J can be derived as:

n o T T T J ¼ E ½GYðkÞ þ BU ðkÞ  Yr ðt Þ ½GYðkÞ þ BUðkÞ  Yr ðt Þ þ ½Ceðk þ 1Þ ½Ceðk þ 1Þ þ ½KUðkÞ ½KUðkÞ

ð34Þ

To ensure that J is minimized, the partial derivative with respect to input UðkÞ should satisfy the following constraint:

@J ¼ 2BT0 ½GYðkÞ þ BUðkÞ  Yr ðkÞ þ 2KT0 KUðkÞ ¼ 0 @UðkÞ 

Thus, the solution of the optimal control input UðkÞ of each sampling step can be derived as:

ð35Þ

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

1  1  UðkÞ ¼ B þ BT0 KT0 K ½Yr ðkÞ  GYðkÞ

13

ð36Þ



It is clear that UðkÞ is dependent on the parameters of the dynamic system expressed in Eq. (29) (identified based on the assumption of a linear relationship between the input and output). However, the hysteresis behavior inherent in the PZT stack would introduce model error and then influence the accuracy of trajectory tracking control. Herein, the adaptive control which adjusts the parameters of Eq. (36) according to Yr ðkÞ,Yðk þ 1Þ, and UðkÞ is used to compensate the model error caused by hysteresis.

   1 Suppose rðk þ 1Þ ¼ Yðk þ 1Þ  Yr ðkÞ þ BT0 KT0 K UðkÞ where Yr ðkÞ,Yðk þ 1Þ, and UðkÞ are known in each k + 1 sampling instant. Then combining with Eq. (33) gives:

rðk þ 1Þ ¼ GYðkÞ þ HUðkÞ  Yr ðkÞ þ Ceðk þ 1Þ

ð37Þ

 1 where H ¼ B þ BT0 KT0 K. Note, G and H are exactly the parameters of Eq. (36), such that the parameter identification of Eq. (37) gives the best solution to the controller parameters of Eq. (36) in each sampling instant k + 1.

5. Experiments The prototype of the developed POS is shown in Fig. 12(b). The structure of POS is made from steel and is manufactured using the wire electrical discharge machining (WEDM) technique. Four identical PZT stacks with the maximum free stroke of 18 lm under the driving voltage of 120 V are mounted on four FMs. Note, the overall structure of POS is integral machined. Hence, no assembly process is required and assembly error is reduced. Fig. 12(a) shows the experimental setup for testing response time, resolution, and trajectory tracking control. The dSPACE DS1103 system acts simultaneously as the control unit and the data acquisition unit. The control input produced by dSPACE is amplified by the high voltage power amplifier (HVPA) with an amplification ratio of 10. The titles about x- or y-axis are measured by laser displacement sensors system in which two sensors (A and B) are set to differential mode. The sampling time interval was set as t = 0.0001 s and the tilt ranges of hx and hy were validated to be approximately ± 1.21mrad and ± 1.19mrad, respectively.

5.1. Response time Fast response is required for antenna pointing. Hence, tests of step responses under MVSTC were implemented. The test results of rotation about the x- and y-axes are shown in Figs. 13 and 14, respectively. The results show that the 3% settling time of the step responses of rotations about x- and y-axes are approximately 26 ms and 21 ms, respectively.

Fig. 12. Experimental setup for tests.

14

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

Fig. 13. Step response of rotation about x-axis under MVSTC.

Fig. 14. Step response of rotation about y-axis under MVSTC.

5.2. Resolution The resolutions of hx and hy were tested by tracking the multi-step trajectory under MVSTC. The corresponding test results are shown in Figs. 15 and 16, respectively. The minimum resolutions of hx and hy can be resolved to be approximately 15lrad. 5.3. Trajectory tracking control The 2-DOF tracking performance of POS was tested by tracking a circular trajectory centered at point (0,0) with a radius of 600lrad and an inclined elliptical trajectory centered at point (0,0), which can be expressed as Eqs. (38) and (39), respectively.



hx ¼ 600sinð2pt Þ

hy ¼ 600cosð2pt Þ



hx ¼ 600sinð2pt Þ

hy ¼ 600cosð2pt þ p=5Þ

ð38Þ

ð39Þ

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

Fig. 15. Resolution of rotation about x-axis under MVSTC.

Fig. 16. Resolution of rotation about y-axis under MVSTC.

Fig. 17. Circular trajectory tracking under MVSTC.

15

16

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

The test data over a four second period are plotted and the discrepancies between the actual trajectories and their references for two trajectories are shown in Figs. 17 and 18, respectively. According to the tracking error depicted in Figs. 17(b) and 18(b), the maximal tracking errors of both circular and inclined elliptical trajectories are less then 10lrad. Additionally, a proportional relationship between the tracking error and the amplitude of trajectory is observed.

Fig. 18. Elliptical trajectory tracking under MVSTC.

Fig. 19. Independent PID control schematic of POS

Fig. 20. Circular trajectory tracking under independent PID control

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

17

Fig. 21. Elliptical trajectory tracking under independent PID control

For comparing study, the 2-DOF tracking of the same circular and elliptical trajectories expressed as Eqs. (38) and (39) by using independent PID method is also carried out. In the independent PID control, the two dofs are under independent PID control where the cross-coupling motions are considered as disturbances, and the controller parameters for each axis are Px = 0.051, Ix = 20.56 and Dx = 0.00001 for rotation about x-axis, and Py = 0.052, Iy = 20.85 and Dy = 0.00001 for rotation about y-axis. Fig. 19 shows the independent PID control schematic of POS, where Gxx ðsÞ and Gyy ðsÞ are the transfer functions of rotation about x- and y-axes, respectively, Gxy ðsÞ and Gyx ðsÞ represent the cross-coupling transfer functions. The test results of tracking of circular and elliptical trajectories are shown in Figs. 20 and 21, respectively. According to the tracking error depicted in Figs. 20(b) and 21(b), the maximal tracking errors of both circular and inclined elliptical trajectories are close to 20lrad, which shows the degradation of the pointing accuracy compared to MVSTC. 6. Conclusions This paper presents a 2-DOF POS for space antenna pointing, in which the rhombic amplifier and 2D flexure hinge are used to enlarge the workspace and decouple the coupling motion between axes. Based on the FEA and SVR modelling method, the trade off between pointing range and structure stiffness is achieved using an iterative optimization design. In addition, to compensate the hysteresis and remaining coupling between axes, MIMO adaptive control, in which the real-time nonliear model of the controlled plant with hysteresis property is approached by the CARMA model, is used to improve the control accuracy. The validation tests show that the proposed POS can achieve approximately ± 1.2mrad tilt range for hx and hy , as well as the fundamental frequency of 472 Hz. Under the MIMO adaptive control, the trajactory tracking error is less than 10lrad and the minimum resolution can be resolved to be approximately 15lrad. Acknowledgement This work is supported by the National Natural Science Foundation of China (Grant Nos. 11902244 and 11872050). References [1] Mickaël Ménez, Xavier Deplancq, Baptiste Palacin and Hugo Gonzalez Perez. Antenna pointing error effects on global capacity for broadband ACM systems, 28th AIAA International Communications Satellite Systems Conference. 30 August-2 September 2010, Anaheim, California. [2] H. Gonzalez-Perez, C. Laporte and B. Palacin. CNES and multibeam antennas, a decade of joint venture, The 8th European Conference on Antennas and Propagation, (2014) 175-177. [3] Jun-Hee Moon, Heui Jae Pahk, Bong-Gu Lee, Design, modeling, and testing of a novel 6-DOF micropositioning stage with low profile and low parasitic motion, Int. J. Adv. Manuf. Technol. 55 (2011) 163–176. [4] K.H. Cai, Y.L. Tian, F.J. Wang, D.W. Zhang, X.P. Liu, B. Shirinzadeh, Design and control of a 6-degree-of-freedom precision positioning system, Rob. Comput. Integr. Manuf. 44 (2017) 77–96. [5] K.H. Cai, Y.L. Tian, X.P. Liu, S. Fatikow, F.J. Wang, L.Y. Cui, B. Shirinzadeh, Modeling and controller design of a 6-DOF precision positioning system, Mech. Syst. Sig. Process. 104 (2018) 536–555. [6] J.L. Qu, W.H. Chen, J.B. Zhang, W.J. Chen, A piezo-driven 2-DOF compliant micropositioning stage with remote center of motion, Sens. Actuat. A 239 (2016) 114–126. [7] H.W. Ma, S.M. Yao, L.Q. Wang, Z. Zhong, Analysis of the displacement amplification ratio of bridge-type flexure hinge, Sens. Actuat. A 132 (2006) 730– 736. [8] K.Q. Qi, Y. Xiang, C. Fang, Y. Zhang, C.S. Song, Analysis of the displacement amplification ratio of bridge-type mechanism, Mech. Mach. Theory 87 (2015) 45–56. [9] J.L. Chen, C.L. Zhang, M.L. Xu, Y.Y. Zi, X.N. Zhang, Rhombic micro-displacement amplifier for piezoelectric actuator and its linear and hybrid model, Mech. Syst. Sig. Process. 50–51 (2015) 580–593.

18

S. Shao et al. / Mechanical Systems and Signal Processing 138 (2020) 106525

[10] S.B. Shao, M.L. Xu, S.W. Zhang, S.L. Xie, Stroke maximizing and high efficient hysteresis hybrid modeling for a rhombic piezoelectric actuator, Mech. Syst. Sig. Process. 75 (2016) 631–647. [11] K.B. Choi, J.J. Leea, S. Hata, A piezo-driven compliant stage with double mechanical amplification mechanisms arranged in parallel, Sens. Actuat. A 161 (2010) 173–181. [12] L.J. Lai, Z.N. Zhu, Design, modeling and testing of a novel flexure-based displacement amplification mechanism amplification mechanism, Sens. Actuat. A 266 (2017) 122–129. [13] J. Guo, S.K. Chee, T.K. Yano, T.S. Higuchi, Micro-vibration stage using piezo actuators, Sens. Actuat. A 194 (2013) 119–127. [14] H.J. Lee, H.C. Kim, H.Y. Kim, D.G. Gweon, Optimal design and experiment of a three-axis out-of-plane nano positioning stage using a new compact bridge-type displacement amplifier, Rev. Sci. Instrum. 84 (2013) 115103. [15] M.X. Ling, L.L. Howell, J.Y. Cao, A pseudo-static model for dynamic analysis on frequency domain of distributed compliant mechanisms, J. Mech. Rob. 10 (2018). [16] M.X. Ling, J.Y. Cao, Z. Jiang, M.H. Zeng, Q.S. Li, Optimal design of a piezo-actuated 2-DOF millimeter-range monolithic flexure mechanism with a pseudo-static model, Mech. Syst. Sig. Process. 115 (2019) 120–131. [17] W.L. Zhu, Z.W. Zhu, P. Guo, B.F. Ju, A novel hybrid actuation mechanism based XY nanopositioning stage with totally decoupled kinematics, Mech. Syst. Sig. Process. 99 (2018) 747–759. [18] S.B. Shao, Z. Tian, S.Y. Song, M.L. Xu, Two-degrees-of-freedom piezo-driven fast steering mirror with cross-axis decoupling capability, Rev. Sci. Instrum. 89 (2018) 055003. [19] Y.Q. Xie, Y.H. Tan, R.L. Dong, Nonlinear modeling and decoupling control of XY micropositioning stages with piezoelectric actuators, IEEE/ASME Trans. Mechatron. 18 (3) (2013). [20] W. Zhu, F.F. Yang, X.T. Rui, Robust independent modal space control of a coupled nano-positioning piezo-stage, Mech. Syst. Sig. Process. 106 (2018) 466–478. [21] G.Y. Gu, L.M. Zhu, C.Y. Su, H. Ding, S. Fatikow, Modeling and control of piezo-actuated nanopositioning stages-a survey, IEEE Trans. Autom. Sci. Eng. 13 (1) (2016). [22] Y.H. Yu, Z.H. Xiao, N.G. Naganathon, R.V. Dukkipati, Dynamic Preisach modelling of hysteresis for the piezoceramic actuator system, Mech. Mach. Theory 37 (2002) 75–89. [23] F. Stefanski, B. Minorowicz, J. Persson, A. Plummer, C. Bowen, Non-linear control of a hydraulic piezo-valve using a generalised Prandtl-Ishlinskii hysteresis model, Mech. Syst. Sig. Process. 82 (2017) 412–431. [24] C.J. Lin, P.T. Lin, Tracking control of a biaxial piezo-actuated positioning stage using generalized Duhem model, Comput. Math. Appl. 64 (5) (2012) 766– 787. [25] W. Zhu, D.H. Wang, Non-symmetrical Bouc-Wen model for piezoelectric ceramic actuators, Sens. Actuat. A 181 (2012) 51–60. [26] A. Milecki, M. Pelic, Application of geometry-based hysteresis modelling in compensation of hysteresis of piezo bender actuator, Mech. Syst. Sig. Process. 78 (2016) 4–17. [27] Y. Koseki, T. Tanikawa, T. Arai and N. Koyachi. Kinematic analysis of translational 3-DOF micro parallel mechanism using matrix Method, in: Proceedings of the 2000 IEEE/RSJ International Conference on Intelligent Robots and Systems, (2000),pp. 786–792. [28] X.F. Mao, Y.J. Wang, X.D. Liu, Y.G. Guo, A hybrid feedforward-feedback fysteresis compensator in piezoelectric actuators based on least-squares support vector machine, IEEE Trans. Ind. Electron. 65, NO (JULY 2018.) 7. [29] L.L. Howell, Compliant Mechanisms, John Wiley & Sons, 2001. [30] B.P. Trease, Y.M. Moon, S. Kota, Design of large-displacement compliant joints, J. Mech. Des. 127 (4) (2005) 788–798. [31] E.G. Merriam, J.E. Jones, S.P. Magleby, L.L. Howell, Monolithic 2 DOF fully compliant space pointing mechanism, Mech. Sci. 4 (2) (2013) 381–390. [32] D.F. Machekposhti, N. Tolou, J.L. Herder, A review on compliant joints and rigid-body constant velocity universal joints toward the design of compliant homokinetic couplings, J. Mech. Des. 137 (3) (2015) 032301. [33] D.F. Machekposhti, N. Tolou, J.L. Herder, A fully compliant homokinetic coupling, J. Mech. Des. 140 (1) (2018) 012301. [34] T. William, Thomson and Marie Dillon Dahleh. Theory of vibration with applications (In Chinese), Tsinghua University Press, 2005.