Ultrasound carotid plaque segmentation via image reconstruction-based self-supervised learning with limited training labels

Ran Zhou; Yanghan Ou; Xiaoyue Fang; M. Reza Azarpazhooh; Haitao Gan; Zhiwei Ye; J. David Spence; Xiangyang Xu; Aaron Fenster; Ran Zhou; Yanghan Ou; Xiaoyue Fang; M. Reza Azarpazhooh; Haitao Gan; Zhiwei Ye; J. David Spence; Xiangyang Xu; Aaron Fenster

doi:10.3934/mbe.2023074

Mathematical Biosciences and Engineering

2023, Volume 20, Issue 2: 1617-1636. doi: 10.3934/mbe.2023074

Previous Article Next Article

Research article Special Issues

Ultrasound carotid plaque segmentation via image reconstruction-based self-supervised learning with limited training labels

1.
School of Computer Science, Hubei University of Technology, Wuhan, China
2.
Robarts Research Institute, Western University, London, Canada
3.
Liyuan Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
Academic editor: Biswajeet Pradhan

Received: 09 August 2022 Revised: 09 October 2022 Accepted: 13 October 2022 Published: 03 November 2022

Carotid total plaque area (TPA) is an important contributing measurement to the evaluation of stroke risk. Deep learning provides an efficient method for ultrasound carotid plaque segmentation and TPA quantification. However, high performance of deep learning requires datasets with many labeled images for training, which is very labor-intensive. Thus, we propose an image reconstruction-based self-supervised learning algorithm (IR-SSL) for carotid plaque segmentation when few labeled images are available. IR-SSL consists of pre-trained and downstream segmentation tasks. The pre-trained task learns region-wise representations with local consistency by reconstructing plaque images from randomly partitioned and disordered images. The pre-trained model is then transferred to the segmentation network as the initial parameters in the downstream task. IR-SSL was implemented with two networks, UNet++ and U-Net, and evaluated on two independent datasets of 510 carotid ultrasound images from 144 subjects at SPARC (London, Canada) and 638 images from 479 subjects at Zhongnan hospital (Wuhan, China). Compared to the baseline networks, IR-SSL improved the segmentation performance when trained on few labeled images (n = 10, 30, 50 and 100 subjects). For 44 SPARC subjects, IR-SSL yielded Dice-similarity-coefficients (DSC) of 80.14–88.84%, and algorithm TPAs were strongly correlated ( $r = 0.962 - 0.993$ , $p$ < 0.001) with manual results. The models trained on the SPARC images but applied to the Zhongnan dataset without retraining achieved DSCs of 80.61–88.18% and strong correlation with manual segmentation ( $r = 0.852 - 0.978$ , $p$ < 0.001). These results suggest that IR-SSL could improve deep learning when trained on small labeled datasets, making it useful for monitoring carotid plaque progression/regression in clinical use and trials.

Keywords:

Citation: Ran Zhou, Yanghan Ou, Xiaoyue Fang, M. Reza Azarpazhooh, Haitao Gan, Zhiwei Ye, J. David Spence, Xiangyang Xu, Aaron Fenster. Ultrasound carotid plaque segmentation via image reconstruction-based self-supervised learning with limited training labels[J]. Mathematical Biosciences and Engineering, 2023, 20(2): 1617-1636. doi: 10.3934/mbe.2023074

Related Papers:

[1]	Qingyun Zhang, Xinhua Zhao, Liang Liu, Tengda Dai . Dynamics analysis of spatial parallel robot with rigid and flexible links. Mathematical Biosciences and Engineering, 2020, 17(6): 7101-7129. doi: 10.3934/mbe.2020365
[2]	Zhen Liu, Song Yang, Tao Ding, Ruimin Chai . Dynamic modeling and performance analysis of a lower-mobility parallel robot with a rotatable platform. Mathematical Biosciences and Engineering, 2023, 20(2): 3918-3943. doi: 10.3934/mbe.2023183
[3]	Jia Tan, ShiLong Chen, ZhengQiang Li . Robust tracking control of a flexible manipulator with limited control input based on backstepping and the Nussbaum function. Mathematical Biosciences and Engineering, 2023, 20(12): 20486-20509. doi: 10.3934/mbe.2023906
[4]	Yunxia Wei, Yuanfei Zhang, Bin Hang . Construction and management of smart campus: Anti-disturbance control of flexible manipulator based on PDE modeling. Mathematical Biosciences and Engineering, 2023, 20(8): 14327-14352. doi: 10.3934/mbe.2023641
[5]	Tongyu Wang, Yadong Chen . Event-triggered control of flexible manipulator constraint system modeled by PDE. Mathematical Biosciences and Engineering, 2023, 20(6): 10043-10062. doi: 10.3934/mbe.2023441
[6]	Dawei Ren, Xiaodong Zhang, Shaojuan Lei, Zehua Bi . Research on flexibility of production system based on hybrid modeling and simulation. Mathematical Biosciences and Engineering, 2021, 18(1): 933-949. doi: 10.3934/mbe.2021049
[7]	Xiaoyu Du, Yijun Zhou, Lingzhen Li, Cecilia Persson, Stephen J. Ferguson . The porous cantilever beam as a model for spinal implants: Experimental, analytical and finite element analysis of dynamic properties. Mathematical Biosciences and Engineering, 2023, 20(4): 6273-6293. doi: 10.3934/mbe.2023270
[8]	Dan Yang, Zhiqiang Xie, Chunting Zhang . Multi-flexible integrated scheduling algorithm for multi-flexible integrated scheduling problem with setup times. Mathematical Biosciences and Engineering, 2023, 20(6): 9781-9817. doi: 10.3934/mbe.2023429
[9]	Xinyu Shao, Zhen Liu, Baoping Jiang . Sliding-mode controller synthesis of robotic manipulator based on a new modified reaching law. Mathematical Biosciences and Engineering, 2022, 19(6): 6362-6378. doi: 10.3934/mbe.2022298
[10]	Zongchao Liu, Linhui Wu, Junwei Yang, Fangsen Cui, Pei Ho, Liping Wang, Jianghui Dong, Gongfa Chen . Thoracic aorta stent grafts design in terms of biomechanical investigations into flexibility. Mathematical Biosciences and Engineering, 2021, 18(1): 800-816. doi: 10.3934/mbe.2021042

Abstract

1. Introduction

Compared with serial manipulators, parallel manipulators offer the benefits of large stiffness, excellent dynamic performance, and small cumulative error ^[1]. With the rapid technological development in recent years, they are being widely used in motion simulation, aeronautical manufacturing, high-speed machining, medicine, and other fields ^[2,3,4,5,6]. However, it suffers from certain disadvantages such as a small workspace with many singularities and the difficulty in solving dynamics problems.

The singularity of the parallel manipulator can result in motion instability and the decline of stiffness and accuracy. Therefore, the design and application of parallel manipulators should avoid singularities as far as possible, particularly the second type. Numerous studies conducted on this subject have proposed certain methods to avoid singularities ^[7,8,9,10]. Specifically, the drive redundancy method ^[4,11] can effectively avoid singularities without increasing the degree of freedom (DOF) of the system.

To adapt to the working environment of high-speed motion in certain applications, the parallel manipulator frequently adopts a lightweight design. However, certain components may experience flexible deformation with the accumulation of the high-speed motion of the manipulator movement with time. Consequently, the dynamic performance and motion accuracy of the entire system is affected. In recent years, many studies have focused on the dynamic modeling of flexible parallel manipulators. Liang et al. studied a 2-DOF multi-drive parallel manipulator, established a rigid-flexible coupling dynamic model, and conducted control research ^[12,13]. Chen et al. studied a 3-RRR parallel manipulator, established a rigid-flexible coupling dynamic model, and conducted dynamic research ^[14]. Rosyid et al. studied a 3-PRR asymmetric planar parallel manipulator and established a rigid-flexible coupling dynamic model ^[15]. Sheng et al. established a rigid-flexible coupling dynamic model of a parallel manipulator with flexible intermediate links and verified it via modal analysis ^[16].

In addition, there are many scholars who have carried out many meaningful studies on the control strategy of flexible parallel robots ^{[17,18,19,20]}, which are of great significance for improving the accuracy of trajectory tracking and positioning of flexible parallel robots.

Owing to the modeling being based on a flexible multi-body system, the final result is generally a complex time-variant nonlinear partial differential equation. Thus, obtaining the analytical solution for complex flexible systems is challenging. Generally, the flexible body is discretized to reduce the complexity of the research, primarily through the following methods: 1) finite element method (FEM) ^[21,22], 2) assumed mode method (AMM) ^[23,24], 3) finite segment method (FSM) ^[25], and 4) lumped parameter method ^[26]; FEM and AMM have been used extensively. A study compared these methods, and the results were found to be similar ^[27]. However, FEM can be applied to the flexible members with complex shapes. More accurate system dynamics model can be obtained when the elements are reasonably divided; therefore, the method is suitable for the dynamic analysis of the system. In contrast, AMM is more suitable for the study of active control.

Shang et al. studied a 3-RRR parallel manipulator ^[28], established the kinematics and rigid body dynamics models, and studied the control strategy. However, this study did not consider the influence of elastic deformation of components on the dynamic performance of the entire system when the manipulator was in high-speed motion.

Compared with the previous work, this paper adopts the floating reference coordinate system and finite element method, combined with the augmented Lagrange multiplier method, to establish a dynamic model of the parallel robot composed of fully flexible Levers, and analyzes the dynamic performance of the simplified model in three driving modes. The modeling adopts a modular approach, which is conducive to generalizing to other models and is more conducive to programming, and the dynamic model of the fully flexible Lever system established is conducive to the optimization between the lightweight design and dynamic performance of the parallel robot, and also helps the subsequent control strategy research.

This study examined a 3-RRR flexible parallel manipulator. To facilitate modeling, a floating reference coordinate system was used to describe the deformed body. The methodology of the study conducted is as follows: In Section 2, The 3-RRR parallel mechanism was described. In Section 3, FEM was used to discretize a single flexible body and establish the dynamic equation of the module. The dynamic equation of a single flexible body was obtained by assembling the module. Further, the flexible bodies were assembled to obtain the dynamic equation of the open chain system. Subsequently, the constraint equation was established and combined with the dynamics equation of the open chain system to form a complete system dynamics equation. In Section 4, the driving torques in three different modes obtained from the inverse dynamics of the rigid body were used as feed-forward input, and the dynamic response was obtained using the ode23tb solver. In Section 5, the numerical simulation results were compared with Adams software simulation results to verify the correctness of the model. In Section 6, finally, certain conclusions were drawn from the results obtained.

2. Model overview

Figure 1 shows the structure of the 3-RRR parallel manipulator. Specifically, A_i is the position of the active joint of the i^th branch chain, whose coordinates are marked as $({x_{Ai}}, {y_{Ai}})$ and angle as ${\alpha _i}$ ; B_i is the position of slave joint of the i^th branch chain, whose coordinates are marked as $({x_{Bi}}, {y_{Bi}})$ and angle as ${\beta _i}$ ; the position of the end-effector is marked as P and its coordinates as $({x_P}, {y_P})$ ; and the coordinates of the active joint are ${A_1}(0,289)$ , ${A_2}(250, - 144)$ and ${A_3}(- 250, - 144)$ . A redundant parallel manipulator is formed when A₁, A₂ and A₃ function as the driving joint, and a non-redundant parallel manipulator formed when A₁ and A₂ function as the driving joint.

Figure 1. Flexible deformation of 3-RRR parallel manipulator in motion.

DownLoad: Full-Size Img PowerPoint

3. Dynamics modeling of rigid-flexible coupling system

To improve their dynamic performance, parallel manipulators often adopt a lightweight design, which results in the deformation of certain components when the manipulator works at high speed. Consequently, the dynamic performance of the manipulator is affected. Therefore, certain components of the manipulator must be flexible and a rigid-flexible coupling dynamic model needs to be established.

To facilitate the modeling, first, the slave joint and end effector were disconnected, as shown in Figure 2, to form an independent flexible body. Subsequently, the dynamics equation of each independent flexible body was derived. Finally, various parts were assembled to form a complete system dynamics equation.

Figure 2. Schematic of 3-RRR parallel manipulator hinge uncoupling.

DownLoad: Full-Size Img PowerPoint

3.1. Dynamics equation of single flexible body j

For the convenience of description, all levers were numbered uniformly as $j(j = 1, 2, \cdots, 6)$ . As the section size of the lever was much smaller than the length of the lever, the Euler–Bernoulli beam theory was applied.

This study employed the floating reference coordinate system to describe the deformed body; that is, the motion of the flexible body was divided into two parts: elastic deformation in the body-connected coordinate system and large-scale motion of the entire body-connected coordinate system relative to the inertial reference system, as shown in . In this figure, $O - x - y$ is the global coordinate system of the system where the rack is located, ${O_j} - x' - y'$ is the flexible body-connected coordinate system, P is any point along the axis of the body-connected coordinate system, ${{\boldsymbol{r}}_{{O_j}}}$ is the radius vector of the origin of the body-connected coordinate system in the global coordinate system, ${{\boldsymbol{r}}_{O, {P_j}}} = {\left[{\begin{array}{*{20}{c}} {x'} & 0 \end{array}} \right]^{\text{T}}}$ is the radius vector in the body-connected coordinate system ${O_j} - x' - y'$ of any point P on the flexible body j before deformation, ${\delta _P} = {\left[{\begin{array}{*{20}{c}} {{u_j}} & {{v_j}} \end{array}} \right]^{\text{T}}}$ expresses the deformation displacement vector from P or P' in the body-connected coordinate system, and ${{\boldsymbol{r}}_P}$ is the position vector of the position P' of the deformed flexible body in the global coordinate system.

Figure 3. Deformation of any flexible body j in motion.

DownLoad: Full-Size Img PowerPoint

According to the closed vector relationship shown in Figure 3, the following equation is obtained:

${\boldsymbol{r}}_{Pi}^j = {{\boldsymbol{r}}_{{O_j}}} + {\boldsymbol{R}}({\varphi _j})({{\boldsymbol{r}}_{O, {P_j}}} + {{\boldsymbol{\delta }}_P}) ,$

(1)

where ${\boldsymbol{R}}({\varphi _j}) = \left[{\begin{array}{*{20}{c}} {\cos {\varphi _j}} & { - \sin {\varphi _j}} \\ {\sin {\varphi _j}} & {\cos {\varphi _j}} \end{array}} \right]$ is the rotation transformation matrix of the body-connected coordinate system relative to the global coordinate system, and ${\varphi _j}$ is the rotational angular displacement of the flexible body j.

3.1.1. Discretization of flexible body j

FEM was used to discretize the flexible body j. For this purpose, the beam element was selected, each containing the left and right nodes and each node containing three generalized coordinates ^[29], as shown in Figure 4.

Figure 4. Discretization of finite element of arbitrary flexible body j and one element i.

DownLoad: Full-Size Img PowerPoint

Here, $O_{1i}^j - x{'_i} - y{'_i}$ is the coordinate system of the module i, $w_{1, i}^j$ , $w_{2, i}^j$ and $w_{3, i}^j$ represent the axial elastic displacement, transverse elastic displacement, and elastic rotational angle of the node i, respectively. The generalized coordinate array of the flexible element i is expressed as:

${\boldsymbol{w}}_i^j = {\left[ {\begin{array}{*{20}{l}} {w_{1, i}^j}&{w_{2, i}^j}&{w_{3, i}^j}&{w_{1, i + 1}^j}&{w_{2, i + 1}^j}&{w_{3, i + 1}^j} \end{array}} \right]^{\text{T}}} .$

(2)

The deformation displacement vector of any point $P_i^j$ on the element i in the flexible body j in the element coordinate system is expressed as:

${\boldsymbol{\delta }}_{Pi}^j = \left[ {\begin{array}{*{20}{c}} {u_i^j(x{'_i}, t)} \\ {v_i^j(x{'_i}, t)} \end{array}} \right] = {\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{w}}_i^j ,$

(3)

where ${\boldsymbol{S}}_i^j(x{'_i})$ is the type function matrix of the element i, given by:

${\boldsymbol{S}}_i^j(x{'_i}) = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{S}}_{i, 1}^j(x{'_i})} \\ {{\boldsymbol{S}}_{i, 2}^j(x{'_i})} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{N}}_{i, 1}^j}&0&0&{{\boldsymbol{N}}_{i + 1, 1}^j}&0&0 \\ 0&{{\boldsymbol{N}}_{i, 2}^j}&{{\boldsymbol{N}}_{i, 3}^j}&0&{{\boldsymbol{N}}_{i + 1, 2}^j}&{{\boldsymbol{N}}_{i + 1, 3}^j} \end{array}} \right] ,$

(4)

where $N_{i, 1}^j = 1 - \frac{{x{'_i}}}{{l_i^j}}$ , $N_{i, 2}^j = 1 - 3{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)^2} + 2{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)^3}$ , $N_{i, 3}^j = l_i^j\left[{\frac{{x{'_i}}}{{l_i^j}} - 2{{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)}^2} + {{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)}^3}} \right]$ , $N_{i + 1, 1}^j = \frac{{x{'_i}}}{{l_i^j}}$ , $N_{i + 1, 2}^j = 3{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)^2} - 2{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)^3}$ , and $N_{i + 1, 3}^j = l_i^j\left[{{{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)}^3} - {{\left({\frac{{x{'_i}}}{{l_i^j}}} \right)}^2}} \right]$ . Specifically, $l_i^j$ is the length of the element i in the flexible body j without deformation.

Overall deformation displacement array of flexible body j is expressed as:

${{\boldsymbol{W}}^j} = \left[ {\begin{array}{*{20}{c}} {w_{1, 1}^j}&{w_{2, 1}^j}&{w_{3, 1}^j}& \cdots &{w_{1, {n_j} + 1}^j}&{w_{2, {n_j} + 1}^j}&{w_{3, {n_j} + 1}^j} \end{array}} \right]_{3({n_j} + 1) \times 1}^{{\text{ T}}} .$

(5)

After discretizing the flexible body j via FEM, the displacement vector of any point $P_i^j$ on its upper element i in the global coordinate system after deformation is

${\boldsymbol{r}}_{{p_i}}^j = {{\boldsymbol{r}}_{{O_j}}} + {{\boldsymbol{R}}^j}({\varphi _j})({\boldsymbol{r}}_{O, i}^j + {\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{B}}_i^j\underbrace {{{\boldsymbol{D}}^j}{{\boldsymbol{w}}^j}}_{{{\boldsymbol{W}}^j}}) ,$

(6)

where $\boldsymbol{r}_{\boldsymbol{O}, \boldsymbol{i}}^\boldsymbol{j} = {\left[{\begin{array}{*{20}{c}} {\sum\limits_{\boldsymbol{k} = 0}^{\boldsymbol{i} - 1} {\boldsymbol{l}_\boldsymbol{k}^\boldsymbol{j}} + \boldsymbol{x}{'_\boldsymbol{i}}} & 0 \end{array}} \right]^{\text{T}}}(\boldsymbol{l}_0^\boldsymbol{j} = 0)$ is the coordinate array of any point $P_i^j$ on the element i of the flexible body j in the body-connected coordinate system before the deformation, and ${\boldsymbol{B}}_i^j$ is the Boolean indicator matrix of the element i in the flexible body j.

$\boldsymbol{B}_\boldsymbol{i}^\boldsymbol{j} = {\left[ {\begin{array}{*{20}{c}} {{0_{3 \times 3}}}&{{0_{3 \times 3}}}& \cdots &{{\boldsymbol{I}_{3 \times 3}}}&{{0_{3 \times 3}}}& \cdots &{{0_{3 \times 3}}}&{{0_{3 \times 3}}} \\ {{0_{3 \times 3}}}&{{0_{3 \times 3}}}& \cdots &{{0_{3 \times 3}}}&{{\boldsymbol{I}_{3 \times 3}}}& \cdots &{{0_{3 \times 3}}}&{{0_{3 \times 3}}} \end{array}} \right]_{6 \times 3({\boldsymbol{n}_\boldsymbol{j}} + 1)}} .$

(7)

For the cantilever beam model, the axial and transverse elastic displacements and elastic rotation at the first node are zero according to its boundary conditions. Then, the independent overall deformation displacement array of the flexible body j is expressed as:

${{\boldsymbol{w}}^j} = \left[ {\begin{array}{*{20}{c}} {w_{1, 2}^j}&{w_{2, 2}^j}&{w_{3, 2}^j}& \cdots &{w_{1, n + 1}^j}&{w_{2, n + 1}^j}&{w_{3, n + 1}^j} \end{array}} \right] .$

(8)

Thus, the corresponding global generalized coordinate transformation matrix can be written as:

${{\boldsymbol{D}}^j} = {\left[ {\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{3 \times 3{n_j}}}} \\ {{{\boldsymbol{I}}_{3{n_j} \times 3{n_j}}}} \end{array}} \right]_{_{3({n_j} + 1) \times 3{n_j}}}} .$

(9)

The generalized coordinate array j of an arbitrary flexible body is defined as:

${{\boldsymbol{q}}_j} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_{{O_j}}^{\text{T}}}&{{\varphi _j}}&{{{\boldsymbol{w}}^{j{\text{T}}}}} \end{array}} \right]^{{\text{ T}}}} .$

(10)

To solve the first derivative of time for closed vector Eq (6), the following is considered:

${\boldsymbol{\dot r}}_{{p_i}}^j = {{\boldsymbol{\dot r}}_{{O_j}}} + \frac{{\partial {{\boldsymbol{R}}^j}({\varphi _j})}}{{\partial {\varphi _j}}}{\dot \varphi _j}({\boldsymbol{r}}_{O, i}^j + {\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}{{\boldsymbol{w}}^j}) + {\boldsymbol{R}}({\varphi _j}){\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}{{\boldsymbol{\dot w}}^j} .$

(11)

To facilitate the subsequent calculation, it is transformed into the form of matrix

${\boldsymbol{\dot r}}_{{p_i}}^j = {\boldsymbol{H}}_i^{'j}{{\boldsymbol{\dot q}}_j} ,$

(12)

${\boldsymbol{H}}_i^{'j} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{I}}_{2 \times 2}}}&{\frac{{\partial {{\boldsymbol{R}}^j}({\varphi _j})}}{{\partial {\varphi _j}}}{{\dot \varphi }_j}({\boldsymbol{r}}_{O, i}^j + {\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}{{\boldsymbol{w}}^j})}&{{\boldsymbol{R}}({\varphi _j}){\boldsymbol{S}}_i^j(x{'_i}){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}{{{\boldsymbol{\dot w}}}^j}} \end{array}} \right]_{2 \times (3 + 3{n_j})}}$

(13)

3.1.2. Dynamics modeling of arbitrary element i

A beam element i as shown in Figure 5. According to the results of kinematic analysis, the kinetic energy of the element i can be expressed as:

Figure 5. Beam element i.

DownLoad: Full-Size Img PowerPoint

$T_i^j = \frac{1}{2}\int_{V_i^j} {\rho {\boldsymbol{\dot r}}_{{p_i}}^{j{\text{T}}}} {\boldsymbol{\dot r}}_{{p_i}}^jdV_i^j = \frac{1}{2}\int_{V_i^j} {\rho {\boldsymbol{\dot q}}_j^T} {\boldsymbol{H}'}_i^{j{\text{T}}}{\boldsymbol{H}'}_i^j{{\boldsymbol{\dot q}}_j}dV_i^j = \frac{1}{2}{\boldsymbol{\dot q}}_j^T{\boldsymbol{M}}_i^j{{\boldsymbol{\dot q}}_j} .$

(14)

The mass matrix corresponding to the volume ${\boldsymbol{M}}_i^j$ of the middle continuum of the element $V_i^j$ , with the material density of $\rho$ , is a symmetric positive definite matrix as follows:

${\boldsymbol{M}}_i^j = \int_{V_i^j} {\rho {\boldsymbol{H}'}_i^{j{\text{T}}}} {\boldsymbol{{ H}}'}_i^jdV_i^j ,$

(15)

${\boldsymbol{M}}_i^j = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{M}}_{i, rr}^j}&{{\boldsymbol{M}}_{i, r\varphi }^j}&{{\boldsymbol{M}}_{i, rf}^j} \\ {{\boldsymbol{M}}_{i, \varphi r}^j}&{{\boldsymbol{M}}_{i, \varphi \varphi }^j}&{{\boldsymbol{M}}_{i, \varphi f}^j} \\ {{\boldsymbol{M}}_{i, fr}^j}&{{\boldsymbol{M}}_{i, f\varphi }^j}&{{\boldsymbol{M}}_{i, ff}^j} \end{array}} \right] .$

(16)

To facilitate the subsequent calculation, several commonly used integrals are expressed as follows:

${\boldsymbol{G}}_{i, 1}^j = \int_{V_i^j} {\rho r_{O, i}^j} dV_i^j ,$

${\boldsymbol{G}}_{i, 2}^j = \int_{V_i^j} {\rho {\boldsymbol{S}}_i^j(x{'_i}} ){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}dV_i^j ,$

${\boldsymbol{\tilde k}} = \left[ {\begin{array}{*{20}{c}} 0&1 \\ { - 1}&0 \end{array}} \right] ,$

potential energy analysis of arbitrary element i is determined as follows: Since the 3-RRR workspace is on the $O - x - y$ plane, the influence of gravity can be ignored. Here, only the bending strain energy and tensile and compressible strain energies of the element are considered, while the shear strain energy and buckling strain energy are ignored.

The strain energy of arbitrary element i in the flexible body j is as follows:

${\boldsymbol{U}}_{p, i}^j = \frac{1}{2}\int_0^{l_i^j} {EI(} \frac{{{\partial ^2}v_i^j(x{'_i}, t)}}{{\partial x{'}_i^2}}{)^2}dx{'_i} + \frac{1}{2}\int_0^{l_i^j} {EA(} \frac{{\partial u_i^j(x{'_i}, t)}}{{\partial x{'_i}}}{)^2}dx{'_i} = \frac{1}{2}{\boldsymbol{q}}_j^T{\boldsymbol{K}}_i^j{{\boldsymbol{q}}_j} ,$

(17)

where E is the elastic modulus of the material, I is the cross-sectional moment of inertia of the element, A is the cross-sectional area of the beam element, and ${\boldsymbol{K}}_i^j$ is the stiffness matrix of the beam element,

${\boldsymbol{K}}_i^j = \left[ {\begin{array}{*{20}{l}} {{{\mathit{\pmb{0}}}_{2 \times 2}}}&{{{\mathit{\pmb{0}}}_{2 \times 1}}}&{{{\mathit{\pmb{0}}}_{2 \times 3{n_j}}}} \\ {{{\mathit{\pmb{0}}}_{1 \times 2}}}&0&{{{\mathit{\pmb{0}}}_{1 \times 3{n_j}}}} \\ {{{\mathit{\pmb{0}}}_{3{n_j} \times 2}}}&{{{\mathit{\pmb{0}}}_{3{n_j} \times 1}}}&{{\boldsymbol{K}}_{i, ff}^j} \end{array}} \right] ,$

(18)

where

$\begin{gathered} {\boldsymbol{K}}_{i, ff}^j = \int_0^{l_i^j} {EI{{\boldsymbol{D}}^{jT}}{\boldsymbol{B}}_i^{jT}{{\left( {\frac{{{\partial ^2}{\boldsymbol{G}}_{i, 2}^j(x{'_i})}}{{\partial x{'}_i^2}}} \right)}^{\text{T}}}\left( {\frac{{{\partial ^2}{\boldsymbol{G}}_{i, 2}^j(x{'_i})}}{{\partial x{'}_i^2}}} \right){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}dx{'_i}} + \\ \int_0^{l_i^j} {EA{{\boldsymbol{D}}^{j{\text{T}}}}{\boldsymbol{B}}_i^{j{\text{T}}}{{\left( {\frac{{\partial {\boldsymbol{G}}_{i, 1}^j(x{'_i})}}{{\partial x{'_i}}}} \right)}^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{G}}_{i, 1}^j(x{'_i})}}{{\partial x{'_i}}}} \right){\boldsymbol{B}}_i^j{{\boldsymbol{D}}^j}dx{'_i}} \\ \end{gathered} .$

Dynamics equation of arbitrary element i is determined as follows: Considering the kinetic and potential energies of the element into the Lagrange equation, the following is obtained:

$\frac{d}{{dt}}(\frac{{\partial \zeta _i^j}}{{\partial {{{\boldsymbol{\dot q}}}_j}}}) - \frac{{\partial \zeta _i^j}}{{\partial {{\boldsymbol{q}}_j}}} = {\boldsymbol{F}}_{e, i}^j ,$

(19)

where ${\boldsymbol{F}}_{e, i}^j$ is the generalized external force array of elements, and $\zeta _i^j$ is the Lagrange function.

$\zeta _i^j = T_i^j - U_{p, i}^j = \frac{1}{2}{\boldsymbol{\dot q}}_j^T{\boldsymbol{M}}_i^j{{\boldsymbol{\dot q}}_j} - \frac{1}{2}{\boldsymbol{\dot q}}_j^T{\boldsymbol{K}}_i^j{{\boldsymbol{\dot q}}_j} .$

Further, on rearranging the elements in Eq (19):

$\frac{{\partial \zeta _i^j}}{{\partial {{\boldsymbol{q}}_j}}} = \frac{{\partial (\frac{1}{2}{\boldsymbol{\dot q}}_j^T{\boldsymbol{M}}_i^j{{{\boldsymbol{\dot q}}}_j})}}{{\partial {{\boldsymbol{q}}_j}}} - {\boldsymbol{K}}_i^j{{\boldsymbol{q}}_j} = \frac{1}{2}{({\boldsymbol{I}} \otimes {{\boldsymbol{\dot q}}_j})^T}\frac{{\partial {\boldsymbol{M}}_i^j}}{{\partial {{\boldsymbol{q}}_j}}}{{\boldsymbol{\dot q}}_j} = {\boldsymbol{C}}_{i2}^j{{\boldsymbol{\dot q}}_j} ,$

(20)

$\frac{d}{{dt}}(\frac{{\partial \zeta _i^j}}{{\partial {{{\boldsymbol{\dot q}}}_j}}}) = \frac{d}{{dt}}({\boldsymbol{M}}_i^j{{\boldsymbol{\dot q}}_j}) = \frac{{d{\boldsymbol{M}}_i^j}}{{dt}}{{\boldsymbol{\dot q}}_j} + {\boldsymbol{M}}_i^j{{\boldsymbol{\ddot q}}_j} ,$

(21)

where

$\frac{{d{\boldsymbol{M}}_i^j}}{{dt}}{{\boldsymbol{\dot q}}_j} = \frac{{\partial {\boldsymbol{M}}_i^j}}{{\partial {\boldsymbol{q}}_j^{\text{T}}}}{\boldsymbol{\dot q}}_j^{\text{T}}{{\boldsymbol{\dot q}}_j} = \frac{{\partial {\boldsymbol{M}}_i^j}}{{\partial {\boldsymbol{q}}_j^{\text{T}}}}({{\boldsymbol{\dot q}}_j} \otimes {\boldsymbol{I}}){{\boldsymbol{\dot q}}_j} = {\boldsymbol{C}}_{i1}^j{{\boldsymbol{\dot q}}_j} .$

(22)

Through further rearrangement of Eq (19), the following is obtained:

${\boldsymbol{M}}_i^j{\ddot q_j} + {\boldsymbol{C}}_i^j{\dot q_j} + {\boldsymbol{K}}_i^j{q_j} = {\boldsymbol{F}}_{e, i}^j ,$

(23)

where

${\boldsymbol{C}}_i^j = {\boldsymbol{C}}_{i1}^j - {\boldsymbol{C}}_{i2}^j = \frac{{\partial {\boldsymbol{M}}_i^j}}{{\partial {\boldsymbol{q}}_j^T}}\left( {{{{\boldsymbol{\dot q}}}_j} \otimes {\boldsymbol{I}}} \right) - \frac{1}{2}{\left( {{\boldsymbol{I}} \otimes {{{\boldsymbol{\dot q}}}_j}} \right)^{\text{T}}}\frac{{\partial {\boldsymbol{M}}_i^j}}{{\partial {{\boldsymbol{q}}_j}}} ,$

(24)

where symbol $\otimes$ indicates the Kronecker product, and ${{\boldsymbol{I}}_{3n + 3}}$ indicates the unit diagonal matrix.

${\boldsymbol{C}}_i^j = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{C}}_{i, rr}^j}&{{\boldsymbol{C}}_{i, r\varphi }^j}&{{\boldsymbol{C}}_{i, rf}^j} \\ {{\boldsymbol{C}}_{i, \varphi r}^j}&{{\boldsymbol{C}}_{i, \varphi \varphi }^j}&{{\boldsymbol{C}}_{i, \varphi f}^j} \\ {{\boldsymbol{C}}_{i, fr}^j}&{{\boldsymbol{C}}_{i, f\varphi }^j}&{{\boldsymbol{C}}_{i, ff}^j} \end{array}} \right] .$

Thus far, the mass, stiffness, and centrifugal force/Coriolis force matrices of the dynamics Eq (23) of any element i on the flexible body j have been obtained, and the dynamics equation modeling of the module i has been completed.

3.1.3. Dynamics equation of arbitrary flexible body j

As the dynamics Eq (23) of the element i is established in the generalized coordinate system of the flexible body j, the dynamics equation of the flexible body j can be obtained provided the dynamics equations of all the elements i in the flexible body j are accumulated.

${{\boldsymbol{M}}^j}{\ddot q_j} + {{\boldsymbol{C}}^j}{\dot q_j} + {{\boldsymbol{K}}^j}{q_j} = {{\boldsymbol{F}}^j} ,$

(25)

where

${{\boldsymbol{M}}_j} = \sum\limits_{i = 1}^{{n_j}} {{\boldsymbol{M}}_i^j} \text{,} {{\boldsymbol{C}}_j} = \sum\limits_{i = 1}^{{n_j}} {{\boldsymbol{C}}_i^j} \text{,} {{\boldsymbol{K}}_j} = \sum\limits_{i = 1}^{{n_j}} {{\boldsymbol{K}}_i^j} \text{,} {{\boldsymbol{F}}_j} = \sum\limits_{i = 1}^{{n_j}} {{\boldsymbol{F}}_{e, i}^j} .$

Here, ${n_j}$ is the number of discretized elements in the flexible body j. To facilitate discussion, the discretization number of all flexible bodies was considered same in the subsequent sections; that is ${n_j} = n(j = 1, \cdots, 6)$ .

3.2. Modeling of system dynamics equation

3.2.1. System dynamics equation without constraints

First, the system dynamics equation without constraints is established. The mechanism is composed of six flexible bodies. The generalized coordinate array of the system is written as

${{\boldsymbol{q}}^{(s)}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_{{O_1}}^{\text{T}}}&{{\varphi _1}}&{{{\boldsymbol{w}}^{1{\text{T}}}}}& \cdots &{{\boldsymbol{r}}_{{O_6}}^{\text{T}}}&{{\varphi _6}}&{{{\boldsymbol{w}}^{6{\text{T}}}}} \end{array}} \right] .$

(26)

By assembling the six flexible bodies, the system dynamics equations without constraints can be established as:

${{\boldsymbol{M}}^{(s)}}{{\boldsymbol{\ddot q}}^{(s)}} + {{\boldsymbol{C}}^{(s)}}{{\boldsymbol{\dot q}}^{(s)}} + {{\boldsymbol{K}}^{(s)}}{{\boldsymbol{q}}^{(s)}} = {\boldsymbol{F}}_e^{(s)} ,$

(27)

where

${{\boldsymbol{M}}^{(s)}} = diag(\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_1}}&{{{\boldsymbol{M}}_2}}& \cdots &{{{\boldsymbol{M}}_6}} \end{array}) , {{\boldsymbol{C}}^{(s)}} = diag(\begin{array}{*{20}{c}} {{{\boldsymbol{C}}_1}}&{{{\boldsymbol{C}}_2}}& \cdots &{{{\boldsymbol{C}}_6}} \end{array}) ,$

${{\boldsymbol{K}}^{(s)}} = diag(\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_1}}&{{{\boldsymbol{K}}_2}}& \cdots &{{{\boldsymbol{K}}_6}} \end{array}) , {\boldsymbol{F}}{{}_e^{(s)}} = diag{(\begin{array}{*{20}{c}} {{\boldsymbol{F}}_{e, 1}^{\text{T}}}&{{\boldsymbol{F}}_{e, 2}^{\text{T}}}& \cdots &{{\boldsymbol{F}}_{e, 6}^{\text{T}}} \end{array})^{\text{T}}} .$

3.2.2. Constraint equation of the system

To establish a complete system dynamics equation, certain constraints must be imposed on the joints and disconnected joints of the end effector based on the system dynamics equation without constraints. The constraints include the joint position and link constraints, totaling 18 constraint equations. The specific equations are as follows:

Position coordinate constraint equation:

$\text{Joint}\; A_{1} \; {{\boldsymbol{\varPsi }}_1}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{A_1}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_1}}}} \\ {{y_{{A_1}}}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(28)

$\text{Joint}\; A_{2 } \; {{\boldsymbol{\varPsi }}_2}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{A_2}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_2}}}} \\ {{y_{{A_2}}}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(29)

$\text{Joint}\; A_{3 } \quad {{\boldsymbol{\varPsi }}_3}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{A_3}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_3}}}} \\ {{y_{{A_3}}}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(30)

$\text{Joint}\; B_{1 } \; {{\boldsymbol{\varPsi }}_4}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{B_1}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_1}}}} \\ {{y_{{A_1}}}} \end{array}} \right] - {\boldsymbol{R}}({\varphi _1})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{1, x'}^{end}} \\ {w_{1, y'}^{end}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(31)

$\text{Joint}\; B_{2 } \; {{\boldsymbol{\varPsi }}_5}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{B_2}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_2}}}} \\ {{y_{{A_2}}}} \end{array}} \right] - {\boldsymbol{R}}({\varphi _3})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{3, x'}^{end}} \\ {w_{3, y'}^{end}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(32)

$\text{Joint}\; B_{3 } \; {{\boldsymbol{\varPsi }}_6}\mathop = \limits^\Delta {{\boldsymbol{r}}_{{B_3}}} - \left[ {\begin{array}{*{20}{c}} {{x_{{A_3}}}} \\ {{y_{{A_3}}}} \end{array}} \right] - {\boldsymbol{R}}({\varphi _5})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{5, x'}^{end}} \\ {w_{5, y'}^{end}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(33)

Link constraint equation:

$\text{Branch chain 1}\; {{\boldsymbol{\varPsi }}_7}\mathop = \limits^\Delta \left[ {\begin{array}{*{20}{c}} {{x_{{A_1}}}} \\ {{y_{{A_1}}}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _1})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{1, x'}^{end}} \\ {w_{1, y'}^{end}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _2})\left[ {\begin{array}{*{20}{c}} {{l_2} + w_{2, x'}^{end}} \\ {w_{2, y'}^{end}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{x_P}} \\ {{y_P}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} ,$

(34)

$\text{Branch chain 2}\; {{\boldsymbol{\varPsi }}_8}\mathop = \limits^\Delta \left[ {\begin{array}{*{20}{c}} {{x_{{A_2}}}} \\ {{y_{{A_2}}}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _3})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{3, x'}^{end}} \\ {w_{3, y'}^{end}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _4})\left[ {\begin{array}{*{20}{c}} {{l_2} + w_{4, x'}^{end}} \\ {w_{4, y'}^{end}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{x_P}} \\ {{y_P}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} , _{ }$

(35)

$\text{Branch chain 3}\; {{\boldsymbol{\varPsi }}_9}\mathop = \limits^\Delta \left[ {\begin{array}{*{20}{c}} {{x_{{A_3}}}} \\ {{y_{{A_3}}}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _5})\left[ {\begin{array}{*{20}{c}} {{l_1} + w_{5, x'}^{end}} \\ {w_{5, y'}^{end}} \end{array}} \right] + {\boldsymbol{R}}({\varphi _6})\left[ {\begin{array}{*{20}{c}} {{l_2} + w_{6, x'}^{end}} \\ {w_{6, y'}^{end}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{x_P}} \\ {{y_P}} \end{array}} \right] = {{\mathit{\pmb{0}}}_{2 \times 1}} ,$

(36)

where ${{\boldsymbol{q}}_P} = {\left[{\begin{array}{*{20}{c}} {{x_P}} & {{y_P}} \end{array}} \right]^{\text{T}}}$ is the coordinates of the end effector.

The constraint equation of the system is written as:

${\boldsymbol{\varPsi }}({{\boldsymbol{q}}^{(s)}}, t) = {{\mathit{\pmb{0}}}_{18 \times 1}} ,$

(37)

where ${\boldsymbol{\varPsi }}({{\boldsymbol{q}}^{(s)}}, t) = {\left[{\begin{array}{*{20}{c}} {{\boldsymbol{\psi }}_1^{\text{T}}} & {{\boldsymbol{\psi }}_2^{\text{T}}} & \cdots & {{\boldsymbol{\psi }}_9^{\text{T}}} \end{array}} \right]^{\text{T}}} \in {R^{18}}$ .

To facilitate the operation, the generalized variable array (26) is reorganized into the following form

${{\boldsymbol{\tilde q}}^{(s)}} = \left[ {\begin{array}{*{20}{c}} {\underbrace {{\boldsymbol{r}}_{{O_1}}^{\text{T}} \cdots {\boldsymbol{r}}_{{O_6}}^{\text{T}}}_{{{\boldsymbol{q}}_r}}}&{\underbrace {{\varphi _1} \cdots {\varphi _6}}_{{{\boldsymbol{q}}_\varphi }}}&{\underbrace {{{\boldsymbol{w}}^{1{\text{T}}}} \cdots {{\boldsymbol{w}}^{6{\text{T}}}}}_{{{\boldsymbol{q}}_w}}} \end{array}} \right]_{(18n + 18) \times 1}^{\text{T}} .$

(38)

The mapping between ${{\boldsymbol{q}}^{(s)}}$ and ${{\boldsymbol{\tilde q}}^{(s)}}$ can be expressed as:

${{\boldsymbol{q}}^{(s)}} = {\boldsymbol{\zeta }}{{\boldsymbol{\tilde q}}^{(s)}} ,$

(39)

where ${\boldsymbol{\zeta }}$ is the adjustment matrix between ${{\boldsymbol{q}}^{(s)}}$ and ${{\boldsymbol{\tilde q}}^{(s)}}$ .

${\boldsymbol{\zeta }} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\zeta }}_1^{\text{T}}}&{{\boldsymbol{\zeta }}_2^{\text{T}}}&{{\boldsymbol{\zeta }}_3^{\text{T}}}&{{\boldsymbol{\zeta }}_4^{\text{T}}}&{{\boldsymbol{\zeta }}_5^{\text{T}}}&{{\boldsymbol{\zeta }}_6^{\text{T}}} \end{array}} \right]^{\text{T}}} ,$

(40)

where ${{\boldsymbol{\zeta }}_j} = \left[{\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{2 \times 2(j - 1)}}} & {{{\boldsymbol{I}}_{2 \times 2}}} & {{{\mathit{\pmb{0}}}_{2 \times (18 + 6n - 2j)}}} \\ {{{\mathit{\pmb{0}}}_{1 \times (j + 11)}}} & {{{\boldsymbol{I}}_{1 \times 1}}} & {{{\mathit{\pmb{0}}}_{1 \times (6 + 6n - j)}}} \\ {{{\mathit{\pmb{0}}}_{n \times (nj - n + 18)}}} & {{{\boldsymbol{I}}_{n \times n}}} & {{{\mathit{\pmb{0}}}_{n \times n(6 - j)}}} \end{array}} \right]$ $(j = 1, \cdots, 6)$

Solving the first derivative of time for Eq (37), the following is obtained:

${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t){{\boldsymbol{\dot q}}^{(s)}} + {{\boldsymbol{\psi }}_t} = {{\mathit{\pmb{0}}}_{18 \times 1}} ,$

(41)

where ${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t) = \frac{{\partial {\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t)}}{{\partial {{\boldsymbol{q}}^{(s){\text{T}}}}}}$ is the system constraint Jacobian matrix. As the constraint equation is time independent, therefore ${{\boldsymbol{\psi }}_t} = \frac{{\partial {\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t)}}{{\partial t}} = {{\mathit{\pmb{0}}}_{18 \times 1}}$ . Thus, it can be concluded that

${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t){{\boldsymbol{\dot q}}^{(s)}} = {0_{18 \times 1}} .$

(42)

According to Eq (26), the first derivative of the generalized coordinates of the system can be expressed as:

${{\boldsymbol{\dot q}}^{(s)}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\dot r}}_{{O_1}}^T}&{{{\dot \varphi }_1}}&{{{{\boldsymbol{\dot w}}}^{1{\text{T}}}}}& \cdots &{{\boldsymbol{\dot r}}_{{O_6}}^T}&{{{\dot \varphi }_6}}&{{{{\boldsymbol{\dot w}}}^{6{\text{T}}}}} \end{array}} \right]^{\text{T}}} .$

(43)

The second derivative of the generalized coordinates of the system can be expressed as:

${{\boldsymbol{\ddot q}}^{(s)}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\ddot r}}_{{O_1}}^T}&{{{\ddot \varphi }_1}}&{{{{\boldsymbol{\ddot w}}}^{1{\text{T}}}}}& \cdots &{{\boldsymbol{\ddot r}}_{{O_6}}^T}&{{{\ddot \varphi }_6}}&{{{{\boldsymbol{\ddot w}}}^{6{\text{T}}}}} \end{array}} \right]^{\text{T}}} .$

(44)

Subsequently, while determining the second derivative of time for Eq (37), the following equation is obtained:

${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t){{\boldsymbol{\ddot q}}^{(s)}} = \eta ,$

(45)

where,

$\eta = - \frac{\partial }{{\partial {{\boldsymbol{q}}^{(s){\text{T}}}}}}\left( {\frac{{\partial {\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t)}}{{\partial {{\boldsymbol{q}}^{(s){\text{T}}}}}}{{{\boldsymbol{\dot q}}}^{(s)}}} \right){{\boldsymbol{\dot q}}^{(s)}} - 2\frac{\partial }{{\partial t}}\left( {\frac{{\partial {\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t)}}{{\partial {{\boldsymbol{q}}^{(s){\text{T}}}}}}} \right){{\boldsymbol{\dot q}}^{(s)}} - \frac{{{\partial ^2}{\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t)}}{{\partial {t^2}}} .$

(46)

3.2.3. Complete system dynamics equation

By combining the open-chain system Eq (27) with the constraint Eq (37) and introducing Lagrange multipliers, a complete system dynamics equation can be obtained as follows:

$\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{M}}^{(s)}}{{{\boldsymbol{\ddot q}}}^{(s)}} + {{\boldsymbol{C}}^{(s)}}{{{\boldsymbol{\dot q}}}^{(s)}} + {{\boldsymbol{K}}^{(s)}}{{\boldsymbol{q}}^{(s)}} + {{\boldsymbol{J}}_q}^{\text{T}}({{\boldsymbol{q}}^{(s)}}, t)\lambda = {\boldsymbol{F}}_e^{(s)}} \\ {{\boldsymbol{\psi }}({{\boldsymbol{q}}^{(s)}}, t) = {{\mathit{\pmb{0}}}_{18 \times 1}}} \end{array}} \right. ,$

(47)

where $\lambda \in {\mathbb{R}^{18}}$ is the Lagrange multiplier array.

$\begin{gathered} {\boldsymbol{F}}_e^{(s)} = \left[ {\begin{array}{*{20}{c}} {\underbrace {\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{1 \times 2}}}&{{\tau _1}}&{{{\mathit{\pmb{0}}}_{1 \times 3{n_1}}}} \end{array}}_{{\text{Flexible body 1}}}}&{\underbrace {{{\mathit{\pmb{0}}}_{1 \times (3{n_2} + 3)}}}_{{\text{Flexible body 2}}}}&{\begin{array}{*{20}{c}} {\underbrace {\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{1 \times 2}}}&{{\tau _3}}&{{{\mathit{\pmb{0}}}_{1 \times 3{n_3}}}} \end{array}}_{{\text{Flexible body 3}}}}&{\underbrace {{{\mathit{\pmb{0}}}_{1 \times (3{n_4} + 3)}}}_{{\text{Flexible body 4}}}} \end{array}} \end{array}} \right. \\ {\left. {\begin{array}{*{20}{c}} {\underbrace {\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{1 \times 2}}}&{{\tau _5}}&{{{\mathit{\pmb{0}}}_{1 \times 3{n_5}}}} \end{array}}_{{\text{Flexible body 5}}}}&{\underbrace {{{\mathit{\pmb{0}}}_{1 \times (3{n_6} + 3)}}}_{{\text{Flexible body 6}}}} \end{array}} \right]^{\text{T}}} \\ \end{gathered} .$

(48)

The generalized independent coordinates of the system can be marked as:

${\boldsymbol{q}} = \left[ {\begin{array}{*{20}{c}} {\underbrace {\begin{array}{*{20}{c}} {{x_P}}&{{y_P}} \end{array}}_{{{\boldsymbol{q}}_P}}}&{\underbrace {\begin{array}{*{20}{c}} {{{\boldsymbol{w}}^{1{\text{T}}}}}& \cdots &{{{\boldsymbol{w}}^{{\text{6T}}}}} \end{array}}_{{{\boldsymbol{q}}_w}}} \end{array}} \right]_{(18n + 2) \times 1}^{\text{T}} .$

(49)

To ease the calculation, the constraint Eq (37) is expressed as

${{\boldsymbol{\psi }}_1} \cdots {{\boldsymbol{\psi }}_6}\mathop = \limits^\Delta {{\boldsymbol{\psi }}_a}({{\boldsymbol{q}}_r}, {{\boldsymbol{q}}_\varphi }, {{\boldsymbol{q}}_w}) = {\mathit{\pmb{0}}} ,$

(50)

${{\boldsymbol{\psi }}_7} \cdots {{\boldsymbol{\psi }}_9}\mathop = \limits^\Delta {{\boldsymbol{\psi }}_b}({{\boldsymbol{q}}_P}, {{\boldsymbol{q}}_\varphi }, {{\boldsymbol{q}}_w}) = {\mathit{\pmb{0}}} .$

(51)

By solving the first derivative with respect to time at both ends of Eq (50), the following equation is obtained:

${{\boldsymbol{\dot \psi }}_a}\mathop = \limits^\Delta \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_r^{\text{T}}}}{{\boldsymbol{\dot q}}_r} + \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}}{{\boldsymbol{\dot q}}_\varphi } + \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}}{{\boldsymbol{\dot q}}_w} = 0 ,$

(52)

where $\partial {{\boldsymbol{\psi }}_a}/\partial {\boldsymbol{q}}_r^{\text{T}}$ is the 12^th order reversible square matrix. By rearranging Eq (52), the following equation is obtained:

${{\boldsymbol{\dot q}}_r} = - {(\frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_r^{\text{T}}}})^{ - 1}}\frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}}{{\boldsymbol{\dot q}}_\varphi } - {(\frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_r^{\text{T}}}})^{ - 1}}\frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}}{{\boldsymbol{\dot q}}_w} .$

(53)

Further, by solving the first derivative with respect to time at both ends of Eq (51), the following equation is obtained:

${{\boldsymbol{\dot \psi }}_b}\mathop = \limits^\Delta \frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}}{{\boldsymbol{\dot q}}_\varphi } + \frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}}{{\boldsymbol{\dot q}}_w} - {\boldsymbol{Z}}{{\boldsymbol{\dot q}}_P} = 0 ,$

(54)

where ${\boldsymbol{Z}} = {\left[{\begin{array}{*{20}{c}} {{{\boldsymbol{I}}_2}} & {{{\boldsymbol{I}}_2}} & {{{\boldsymbol{I}}_2}} \end{array}} \right]^{\text{T}}}$ , ${{\boldsymbol{I}}_2}$ is the identity matrix, and $\partial {{\boldsymbol{\psi }}_b}/\partial {\boldsymbol{q}}_\varphi ^{\text{T}}$ is a 6^th order reversible square matrix. Eq (54) can be rearranged as:

${{\boldsymbol{\dot q}}_\varphi } = \left[ {\begin{array}{*{20}{c}} {{{\left( {\frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}}} \right)}^{ - 1}}{\boldsymbol{Z}}}&{ - {{\left( {\frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}}} \right)}^{ - 1}}\frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] .$

(55)

For simplification, it is expressed as follows:

${{\boldsymbol{J}}_{b\varphi }} = \frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_\varphi ^T}} \text{,} {{\boldsymbol{J}}_{bw}} = \frac{{\partial {{\boldsymbol{\psi }}_b}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}} \text{,} {{\boldsymbol{J}}_{ar}} = \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_r^{\text{T}}}} \text{,} {{\boldsymbol{J}}_{a\varphi }} = \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_\varphi ^{\text{T}}}} \text{,} {{\boldsymbol{J}}_{aw}} = \frac{{\partial {{\boldsymbol{\psi }}_a}}}{{\partial {\boldsymbol{q}}_w^{\text{T}}}} .$

(56)

Substituting Eq (55) into Eq (53) and canceling ${{\boldsymbol{\dot q}}_\varphi }$ ,

${{\boldsymbol{\dot q}}_r} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varPhi }}_1}}&{{{\boldsymbol{\varPhi }}_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] = {\boldsymbol{\varPhi }}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] ,$

(57)

where,

${{\boldsymbol{\varPhi }}_1} = - {\left( {{{\boldsymbol{J}}_{ar}}} \right)^{ - 1}}{{\boldsymbol{J}}_{a\varphi }}{\left( {{{\boldsymbol{J}}_{b\varphi }}} \right)^{ - 1}}{\boldsymbol{Z}} ,$

(58)

${{\boldsymbol{\varPhi }}_2} = {\left( {{{\boldsymbol{J}}_{ar}}} \right)^{ - 1}}{{\boldsymbol{J}}_{a\varphi }}{\left( {{{\boldsymbol{J}}_{b\varphi }}} \right)^{ - 1}}{{\boldsymbol{J}}_{bw}} - {\left( {{{\boldsymbol{J}}_{ar}}} \right)^{ - 1}}{{\boldsymbol{J}}_{aw}} ,$

(59)

${{\boldsymbol{\dot q}}_\varphi } = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varGamma }}_1}}&{{{\boldsymbol{\varGamma }}_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] = {\boldsymbol{\varGamma }}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] ,$

(60)

further,

${{\boldsymbol{\varGamma }}_1} = {\left( {{{\boldsymbol{J}}_{b\varphi }}} \right)^{ - 1}}{\boldsymbol{Z}} ,$

(61)

${{\boldsymbol{\varGamma }}_2} = - {\left( {{{\boldsymbol{J}}_{b\varphi }}} \right)^{ - 1}}{{\boldsymbol{J}}_{bw}} ,$

(62)

${{\boldsymbol{\dot q}}_w} = {\boldsymbol{U}}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] ,$

(63)

where, ${\boldsymbol{U}} = \left[{\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{18n \times 2}}} & {{{\boldsymbol{I}}_{18n \times 18n}}} \end{array}} \right]$ .

By combining Eqs. (58), (61) and (64) with the generalized coordinate variable array (38), the following is obtained:

${{\boldsymbol{\dot {\tilde q}}}^{(s)}} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_r}} \\ {{{{\boldsymbol{\dot q}}}_\varphi }} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] = {\boldsymbol{\varTheta }}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot q}}}_P}} \\ {{{{\boldsymbol{\dot q}}}_w}} \end{array}} \right] = {\boldsymbol{\varTheta \dot q}} ,$

(64)

where,

${\boldsymbol{\varTheta }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varPhi }}_{12 \times (18n + 2)}}} \\ {{{\boldsymbol{\varGamma }}_{6 \times (18n + 2)}}} \\ {{{\boldsymbol{U}}_{18n \times (18n + 2)}}} \end{array}} \right] .$

(65)

This is the first-order mapping matrix between the generalized independent coordinate array and generalized coordinate variable array.

Further, solving the first-order derivative of time for Eq (65), the following equation is obtained:

${{\boldsymbol{\ddot {\tilde q}}}^{(s)}} = {\boldsymbol{\dot \varTheta \dot q}} + {\boldsymbol{\varTheta \ddot q}} ,$

(66)

where,

${\boldsymbol{\dot \varTheta }} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\dot \varPhi }}} \\ {{\boldsymbol{\dot \varGamma }}} \\ {{\boldsymbol{\dot U}}} \end{array}} \right] .$

(67)

Each element in Eq (68) can be expressed as:

${\boldsymbol{\dot \varPhi }} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot \varPhi }}}_{\boldsymbol{1}}}}&{{{{\boldsymbol{\dot \varPhi }}}_{\boldsymbol{2}}}} \end{array}} \right] ,$

(68)

${\boldsymbol{\dot \varGamma }} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot \varGamma }}}_1}}&{{{{\boldsymbol{\dot \varGamma }}}_2}} \end{array}} \right] ,$

(69)

${\boldsymbol{\dot U}} = {{\mathit{\pmb{0}}}_{18n \times (18n + 2)}} .$

(70)

The mapping relationship between the independent generalized coordinate variable array (49) and generalized coordinate variable array (38) is

${{\boldsymbol{\tilde q}}^{(s)}} = {{\boldsymbol{\varTheta }}^{'}}{\boldsymbol{q}} ,$

(71)

where,

${{\boldsymbol{\varTheta }}^{'}} = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\pmb{0}}}_{18 \times 2}}}&{{{\mathit{\pmb{0}}}_{18 \times 18n}}} \\ {{{\mathit{\pmb{0}}}_{18n \times 2}}}&{{{\boldsymbol{I}}_{18n \times 18n}}} \end{array}} \right] .$

(72)

To facilitate the calculation, the following settings are made

$\begin{gathered} {{{\boldsymbol{\dot J}}}_{ar}} = [\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{{\boldsymbol{\psi }}_a}}}{{\partial {{({q_{r1}})}^2}}}{{\dot q}_{r1}}}&{\frac{{{\partial ^2}{{\boldsymbol{\psi }}_a}}}{{\partial {{({q_{r2}})}^2}}}{{\dot q}_{r2}}}& \cdots &{\frac{{{\partial ^2}{{\boldsymbol{\psi }}_a}}}{{\partial {{({q_{r12}})}^2}}}{{\dot q}_{r12}}} \end{array}], \hfill \\ {{{\boldsymbol{\dot J}}}_{b\phi }} = [\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{{\boldsymbol{\psi }}_b}}}{{\partial {{({q_{\phi 1}})}^2}}}{{\dot q}_{\phi 1}}}&{\frac{{{\partial ^2}{{\boldsymbol{\psi }}_b}}}{{\partial {{({q_{\phi 2}})}^2}}}{{\dot q}_{\phi 2}}}& \cdots &{\frac{{{\partial ^2}{{\boldsymbol{\psi }}_b}}}{{\partial {{({q_{\phi 6}})}^2}}}{{\dot q}_{\phi 6}}} \end{array}], \hfill \\ {{{\boldsymbol{\dot J}}}_{bw}} = [\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_b}}}{{\partial {{({q_{w1}})}^2}}}{{\dot q}_{w1}}}&{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_b}}}{{\partial {{({q_{w2}})}^2}}}{{\dot q}_{w2}}}& \cdots &{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_b}}}{{\partial {{({q_{w18n}})}^2}}}{{\dot q}_{w18n}}} \end{array}], \hfill \\ {{{\boldsymbol{\dot J}}}_{a\phi }} = [\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{\phi 1}})}^2}}}{{\dot q}_{\phi 1}}}&{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{\phi 2}})}^2}}}{{\dot q}_{\phi 2}}}& \cdots &{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{\phi 6}})}^2}}}{{\dot q}_{\phi 6}}} \end{array}], \hfill \\ {{{\boldsymbol{\dot J}}}_{aw}} = [\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{w1}})}^2}}}{{\dot q}_{w1}}}&{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{w2}})}^2}}}{{\dot q}_{w2}}}& \cdots &{\frac{{{\partial ^2}{{\boldsymbol{\varPsi }}_a}}}{{\partial {{({q_{w18n}})}^2}}}{{\dot q}_{w18n}}} \end{array}]. \hfill \\ \end{gathered}$

(73)

Thus far, the mapping relationship between the generalized coordinate variable array and system independent generalized coordinate array (including the first-order form and second-order form) is established.

$\left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{q}}^{(s)}} = {\boldsymbol{\zeta }}{{{\boldsymbol{\tilde q}}}^{(s)}} = {\boldsymbol{\zeta }}{{\boldsymbol{\varTheta }}^{'}}{\boldsymbol{q}}} \\ {{{{\boldsymbol{\dot q}}}^{(s)}} = {\boldsymbol{\zeta }}{{{\boldsymbol{\dot {\tilde q}}}}^{(s)}} = {\boldsymbol{\zeta \varTheta \dot q}}} \\ {{{{\boldsymbol{\ddot q}}}^{(s)}} = {\boldsymbol{\zeta }}{{{\boldsymbol{\ddot {\tilde q}}}}^{(s)}} = {\boldsymbol{\zeta \dot \varTheta \dot q}} + {\boldsymbol{\zeta \varTheta \ddot q}}} \end{array}} \right. .$

(74)

Substituting Eq (75) into constraint Eq (43) and combining the relationship between ${{\boldsymbol{q}}^{(s)}}$ and ${{\boldsymbol{\tilde q}}^{(s)}}$ , the following relation is obtained:

${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t){\boldsymbol{\xi \varTheta \dot q}} \equiv {\mathit{\pmb{0}}} .$

(75)

As q is a generalized independent coordinate array of the system, the following relationship is always true for any q in the system workspace:

${{\boldsymbol{J}}_q}({{\boldsymbol{q}}^{(s)}}, t){\boldsymbol{\xi \varTheta }} \equiv {\mathit{\pmb{0}}} \Leftrightarrow {({\boldsymbol{\xi \varTheta }})^{\text{T}}}{\boldsymbol{J}}_q^{\text{T}}({{\boldsymbol{q}}^{(s)}}, t) \equiv {\mathit{\pmb{0}}} .$

(76)

To eliminate the generalized constraining force ${\boldsymbol{J}}_q^{\text{T}}({{\boldsymbol{q}}^{(s)}}, t)\lambda$ that is difficult to obtain in the closed-chain system equation, the two ends of the closed-chain system dynamics Eq (48) are multiplied by ${({\boldsymbol{\xi \varTheta }})^{\text{T}}}$ , and ${({\boldsymbol{\xi \varTheta }})^{\text{T}}}{\boldsymbol{J}}_q^{\text{T}}({{\boldsymbol{q}}^{(s)}}, t)\lambda$ is canceled. Thus, the following equation is obtained:

${({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{{\boldsymbol{M}}^{(s)}}{{\boldsymbol{\ddot q}}^{(s)}} + {({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{{\boldsymbol{C}}^{(s)}}{{\boldsymbol{\dot q}}^{(s)}} + {({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{{\boldsymbol{K}}^{(s)}}{{\boldsymbol{q}}^{(s)}} = {({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{\boldsymbol{F}}_e^{(s)} .$

(77)

The mapping relationship (75) is plugged into Eq (78) and rearranged to obtain the complete dynamics equation of the system represented by the generalized independent coordinates of the system,

${\boldsymbol{M\ddot q}} + {\boldsymbol{C\dot q}} + {\boldsymbol{Kq}} = {\boldsymbol{F}} ,$

(78)

where,

${\boldsymbol{M}} = {({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{{\boldsymbol{M}}^{(s)}}{\boldsymbol{\zeta \varTheta }} ,$

${\boldsymbol{C}} = {\left( {{\boldsymbol{\zeta \varTheta }}} \right)^{\text{T}}}{{\boldsymbol{M}}^{(s)}}{\boldsymbol{\zeta \dot \varTheta }} + {\left( {{\boldsymbol{\zeta \varTheta }}} \right)^{\text{T}}}{{\boldsymbol{C}}^{(s)}}{\boldsymbol{\zeta \varTheta }} ,$

${\boldsymbol{K}} = {\left( {{\boldsymbol{\zeta \varTheta }}} \right)^{\text{T}}}{{\boldsymbol{K}}^{(s)}}{\boldsymbol{\zeta \varTheta }}' ,$

${\boldsymbol{F}} = {({\boldsymbol{\zeta \varTheta }})^{\text{T}}}{\boldsymbol{F}}_e^{(s)} .$

This completes the dynamics modeling.

4. Numerical simulation and analysis

To establish a generic dynamic model and facilitate the rigid-flexible coupling dynamic analysis, 1–6 levers are all considered flexible. The study found that the elastic deformation of the driving lever is more evident than that of the follower lever. Thus, to reduce the difficulty of numerical solution, only 1, 3 and 5 levers were considered flexible, and each lever was divided into three sections. Whereas, 2, 4 and 6 levers were regarded as rigid bodies. When the manipulator moves at high speed, its motion essentially includes two parts: large-scale rigid and small-scale flexible motions, which affect each other. The ode23tb solver was used to solve the system dynamics Eq (79) established in the previous section.

In this simulation, the planning path was that the center point P of the terminal console moved from the point (0, 0) to the point (−14.42 mm, 36.15 mm) with a spiral trajectory within 0–0.5 s. The trajectory, velocity, and acceleration of point P are shown in Figures 6 and 7. The eccentric spiral was selected owing to the increasing speed and acceleration of the spiral trajectory, which is convenient for dynamic analysis. In addition, the influence of the asymmetric trajectory of eccentric spiral on the system was studied.

Figure 6. Trajectory of Point P.

DownLoad: Full-Size Img PowerPoint

Figure 7. (a) P velocity of end-effector; (b) P acceleration of end-effector.

DownLoad: Full-Size Img PowerPoint

The parameters of a group of 3-RRR parallel manipulators are listed in Table 1. It was assumed that all the levers and terminal consoles were composed of aluminum with uniform texture, and all joints were light ideal bearings (i.e., massless, frictionless, clearanceless, etc.).

Table 1. Parameters of 3-RRR parallel manipulator (unit: mm).

l₁	l₂	Section A(mm²)	Lever thickness h	$\rho$ (g/mm³)
244	244	40	5	0.0027

| Show Table

DownLoad: CSV

When the end-effector moves in a planned trajectory, the driving torque has a unique solution for the non-redundant drive mode with two driving joints. However, for redundant drive mode with three driving joints, the distribution of driving torque is multi-solution. This study selected three types of driving torque distribution modes: (1) under non-redundant drive (A1, A2 are driving joints); (2) under redundant drive (A1, A2 and A3 are driving joints), the distribution of driving torque with minimum Euclidean norm of internal force of slave joints; (3) under redundant drive (A1, A2, A3 are driving joints), the distribution of driving torque with minimum European norm of driving joint power. The driving torque calculated via the inverse rigid body dynamic model (as shown in Figure 8) was used as the feed-forward input of the rigid-flexible coupling dynamic model to obtain the positive motion solution, thereby obtaining the dynamic response of the system.

Figure 8. Three modes of feed-forward input torque of inverse rigid body dynamics: (a) non-redundant drive mode; (b) minimum Euclidean norm mode of internal force of redundant driving joint; (c) minimum Euclidean norm mode of active joint power.

DownLoad: Full-Size Img PowerPoint

Figures 9 and 10 show that among the three modes, the trajectory of the end-effector is very close to the expected trajectory at the beginning, and the error increases gradually with the passage of time. In particular, the deviation of Mode 1 from the expected trajectory is relatively large at the final stage; this is more clearly observable in Figure 11. The redundant drive adopted in Modes 2 and 3 was found to effectively improve the dynamic performance of the parallel manipulator.

Figure 9. Trajectory of end-effector under the three driving torque distribution modes.

DownLoad: Full-Size Img PowerPoint

Figure 10. Motion error of end-effector under the three driving torque distribution modes: (a) x-direction error; (b) y-direction error.

DownLoad: Full-Size Img PowerPoint

Figure 11. Angular error of three driving torque distribution modes. (a) angular error for Lever 1; (b) angular error for Lever 3; (c) angular error for Lever 5; (d) angular error for Lever 2; (e) angular error for Lever 4; (f) angular error for Lever 6.

DownLoad: Full-Size Img PowerPoint

Figures 12 and 13 show that the oscillation amplitude of angular velocity and angular acceleration of Levers 1, 3 and 5 is much larger than those of Levers 2, 4 and 6. This is because the former are flexible, while the latter are treated as rigid, and the oscillations generated by the latter are caused by the flexible body. In addition, the oscillation amplitude of Mode 1 is much larger than that of Modes 2 and 3, indicating that the dynamic performance of redundant drive is significantly better than that of the non-redundant drive, and the stability of the system driven by redundant drive is significantly improved.

Figure 12. Angular velocity under three driving torque distribution modes. (a) Lever 1; (b) Lever 3; (c) Lever 5; (d) Lever 2; (e) Lever 4; (f) Lever 6.

DownLoad: Full-Size Img PowerPoint

Figure 13. Angular acceleration under three driving torque distribution modes. (a) Lever 1; (b) Lever 3; (c) Lever 5; (d) Lever 2; (e) Lever 4; (f) Lever 6.

DownLoad: Full-Size Img PowerPoint

The elastic deformation at the end of the three flexible bodies can be clearly observed in Figures 14–16. The axial elastic deformation under the three modes is considerably small, i.e., three orders of magnitude smaller than the transverse deformation. Compared with the other two modes, Mode 2 exhibited smaller elastic deformation, more stable change, and better performance. The elastic deformation of Mode 1 in Figure 16 showed a significant decrease mainly owing to the non-redundant drive of Mode 1 and follower Lever 5.

Figure 14. Elastic parameters of Lever 1 under three driving torque distribution modes: (a) axial deformation; (b) transverse deformation; (c) elastic angle.

DownLoad: Full-Size Img PowerPoint

Figure 15. Elastic parameters of Lever 3 under three driving torque distribution modes: (a) axial deformation; (b) transverse deformation; (c) elastic angle.

DownLoad: Full-Size Img PowerPoint

Figure 16. Elastic parameters of Lever 5 under three driving torque distribution modes: (a) axial deformation; (b) transverse deformation; (c) elastic angle.

DownLoad: Full-Size Img PowerPoint

5. Simulation verification

The simulation verification was based on the parametric model of the rigid-flexible coupling system of the 3-RRR parallel manipulator established on Adams. The modeling used a method similar to the finite segment discretization, wherein the three driving levers were discretized (Levers 1, 3 and 5), and each lever was divided into eight rigid segments linked by massless torsion springs (total of 24 segments). The three follower levers (Levers 2, 4 and 6) were not divided; the details are shown in the Figure 17. To facilitate the observation of the section of the lever, the parameters of Lever 1 were hidden and amplified, as shown in Figure 17(b). The driving torque under the previous Mode 2 was used as the feed-forward input of the simulation model, and certain simulation results were compared with the numerical model simulation results of the previous Mode 2, as shown in the Figures 18–21.

Figure 17. Adams simulation model and trajectory.

DownLoad: Full-Size Img PowerPoint

Figure 18. (a) Displacement of the end effector in x-direction; (b) Displacement of the end effector in y-direction.

DownLoad: Full-Size Img PowerPoint

Figure 19. (a) Velocity of end-effector in x-direction; (b) Velocity of end-effector in y-direction.

DownLoad: Full-Size Img PowerPoint

Figure 20. Angular velocity of three flexible driving levers: (a) Lever 1; (b) Lever 3; (c) Lever 5.

DownLoad: Full-Size Img PowerPoint

Figure 21. Angular velocity of three rigid follower levers: (a) Lever 2; (b) Lever 4; (c) Lever 6.

DownLoad: Full-Size Img PowerPoint

To maintain generality, certain dynamic response parameters of the end-effector and each lever in Mode 2 were selected and compared with the simulation results. It was found that the trajectory and angular velocity of the end-effector in Figures 18 and 19 were nearly identical, and the velocity of the three follower levers in Figure 21 was consistent, with small deviations. Further, the velocity of the three driving levers in Figure 20 was consistent in general, with certain deviations. This is mainly because the discretization technology used in the two methods is different. The numerical simulation uses FEM while the software simulation uses FSM. Thus, it was confirmed that the rigid-flexible coupling dynamic model established in this study is accurate.

6. Conclusions

This study determined the configuration of 3-RRR parallel manipulator and performed a flexible treatment to establish a rigid-flexible coupling dynamic model. The modular modeling method is conducive for extension to other models and soft programming. The dynamic response of three different driving torque distribution modes was obtained via numerical simulation. Finally, FSM was used to establish the model in Adams, and the established model was verified. The conclusions of this study are as follows:

(1) Establishing dynamics model of flexible link system

The rod was finitely discretized by the element method, and the dynamic equation of the module was established using the floating reference coordinate system. Further, the module was linked to establish the dynamic equation of the flexible body. All the flexible bodies were combined to obtain the dynamic equation of the open chain system, and a complete system dynamics equation was established. The modular modeling method was applied in the entire process. First, the independent flexible body dynamics equation was established. Second, based on the constraint equation, each segment was grouped together to form a complete equation, which is conducive to the promotion of other modeling and software programming systems, as well as lightweight design and system dynamic performance optimization.

(2) Numerical simulation and analysis

The driving torque in the three modes obtained by inverse rigid body dynamics was used as the feed-forward input of numerical simulation to obtain the dynamic response using the ode23tb solver. The analysis revealed that the redundant drive had better dynamic performance than the non-redundant drive and exhibited a suitable inhibitory effect on vibration. In addition, the dynamic performance of Mode 2 was smoother, with smaller vibration, compared to that of Mode 3. From the flexible deformation perspective, the axial deformation was much smaller than the transverse deformation, and they differed by three orders of magnitude. Moreover, the deformation of the driven flexible body was smaller than that of the active flexible body. Furthermore, the deformation of driving lever in Mode 2 was smaller than that in Modes 1 and 3, and the change was more stable with the performance being better than the other two modes.

(3) Finally, FSM was used to establish the simulation model in Adams.

The driving torque of Mode 2 is selected as the feed-forward input, and certain dynamic response results were compared with the numerical simulation results to verify the correctness of the established rigid-flexible coupling dynamic model.

This paper mainly carries out the dynamic modeling and dynamic analysis of flexible Lever parallel robots.it will be carried out around the control strategy of parallel robots in the follow-up research, which improve its trajectory tracking and positioning accuracy, the potential application fields of the parallel robot including fast packaging, fast sorting, high-speed motion simulation, etc.

Acknowledgments

The authors would like to appreciate the editors and the reviewers for their significant contributions to the improvement of the paper. Thanks to Xi 'an Jiaotong University for their help in this study. In this paper, the MATLAB calculations and ADAMS dynamics simulations are completed in Xi'an Jiaotong University.

Conflict of interest

The authors declare that they have no conflict of interest.

References

[1]	T. Vos, S. S. Lim, C. Abbafati, K. M. Abbas, M. Abbasi, M. Abbasifard, et al., Global burden of 369 diseases and injuries in 204 countries and territories, 1990–2019: A systematic analysis for the global burden of disease study 2019, Lancet, 369 (2020), 1204–1222. https://doi.org/10.1016/S0140-6736(20)30925-9 doi: 10.1016/S0140-6736(20)30925-9
[2]	J. D. Spence, Technology insight: Ultrasound measurement of carotid plaque-patient management, genetic research, and therapy evaluation, Nat. Clin. Pract. Neurol., 2 (2006), 611–619. https://doi.org/10.1038/ncpneuro0324 doi: 10.1038/ncpneuro0324
[3]	M. L. Bots, A. W. Hoes, P. J. Koudstaal, A. Hofman, D. E. Grobbee, Common carotid intima-media thickness and risk of stroke and myocardial infarction: The rotterdam study, Circulation, 96 (1997), 1432–143. https://doi.org/10.1161/01.cir.96.5.1432 doi: 10.1161/01.cir.96.5.1432
[4]	J. D. Spence, Carotid ultrasound phenotypes are biologically distinct, Arterioscler., Thromb., Vasc. Biol., 35 (2015), 1910–1913. https://doi.org/10.1161/ATVBAHA.115.306209 doi: 10.1161/ATVBAHA.115.306209
[5]	A. V. Finn, F. D. Kolodgie, R. Virmani, Correlation between carotid intimal/medial thickness and atherosclerosis: A point of view from pathology, Arterioscler., Thromb., Vasc. Biol., 30 (2010), 177–181. https://doi.org/10.1161/ATVBAHA.108.173609 doi: 10.1161/ATVBAHA.108.173609
[6]	R. L. Pollex, R. A. Hegele, Genetic determinants of carotid ultrasound traits, Curr. Atheroscler. Rep., 8 (2006), 206–215. https://doi.org/10.1007/s11883-006-0075-z doi: 10.1007/s11883-006-0075-z
[7]	E. B. Mathiesen, S. H. Johnsen, T. Wilsgaard, K. H. Bønaa, M. L. Løchen, I. Njølstad, Carotid plaque area and intima-media thickness in prediction of first-ever ischemic stroke: A 10-year follow-up of 6584 men and women: the tromsø study, Stroke, 42 (2011), 972–978. https://doi.org/10.1161/strokeaha.110.589754 doi: 10.1161/strokeaha.110.589754
[8]	C. Loizou, S. Petroudi, M. Pantziaris, A. Nicolaides, C. Pattichis, An integrated system for the segmentation of atherosclerotic carotid plaque ultrasound video, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 61 (2014), 86–101. https://doi.org/10.1109/tuffc.2014.6689778 doi: 10.1109/tuffc.2014.6689778
[9]	J. Cheng, H. Li, F. Xiao, A. Fenster, X. Zhang, X. He, et al., Fully automatic plaque segmentation in 3-D carotid ultrasound images, Ultrasound Med. Biol., 39 (2013), 2431–2446. https://doi.org/10.1016/j.ultrasmedbio.2013.07.007 doi: 10.1016/j.ultrasmedbio.2013.07.007
[10]	F. Destrempes, J. Meunier, M. F. Giroux, G. Soulez, G. Cloutier, Segmentation of plaques in sequences of ultrasonic b-mode images of carotid arteries based on motion estimation and a bayesian model, IEEE Trans. Biomed. Eng., 58 (2011), 2202–2211. https://doi.org/10.1109/tbme.2011.2127476 doi: 10.1109/tbme.2011.2127476
[11]	S. Delsanto, F. Molinari, P. Giustetto, W. Liboni, S. Badalamenti, J. S. Suri, Characterization of a completely user-independent algorithm for carotid artery segmentation in 2-D ultrasound images, IEEE Trans. Instrum. Meas., 56 (2007), 1265–1274. https://doi.org/10.1109/TIM.2007.900433 doi: 10.1109/TIM.2007.900433
[12]	R. M. Menchón-Lara, M. C. Bastida-Jumilla, J. Morales-Sánchez, J. L. Sancho-Gómez, Automatic detection of the intima-media thickness in ultrasound images of the common carotid artery using neural networks, Med. Biol. Eng. Comput., 52 (2014), 169–181. https://doi.org/10.1007/s11517-013-1128-4 doi: 10.1007/s11517-013-1128-4
[13]	C. Azzopardi, Y. A. Hicks, K. P. Camilleri, Automatic carotid ultrasound segmentation using deep convolutional neural networks and phase congruency maps, in 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017), (2017), 624–628. https://doi.org/10.1109/ISBI.2017.7950598
[14]	S. Savaş, N. Topaloğlu, O. Kazcı, P. N. Koşar, Classification of carotid artery intima media thickness ultrasound images with deep learning, J. Med. Syst., 43 (2019), 273. https://doi.org/10.1007/s10916-019-1406-2 doi: 10.1007/s10916-019-1406-2
[15]	C. Qian, E. Su, X. Yang, Segmentation of the common carotid intima-media complex in ultrasound images using 2-D continuous max-flow and stacked sparse auto-encoder, Ultrasound Med. Biol., 46 (2020), 3104–3124. https://doi.org/10.1016/j.ultrasmedbio.2020.07.021 doi: 10.1016/j.ultrasmedbio.2020.07.021
[16]	M. Jiang, Y. Zhao, B. Chiu, Segmentation of common and internal carotid arteries from 3D ultrasound images based on adaptive triple loss, Med. Phys., 48 (2021), 5096–5114. https://doi.org/10.1002/mp.15127 doi: 10.1002/mp.15127
[17]	R. Zhou, A. Fenster, Y. Xia, J. D. Spence, M. Ding, Deep learning-based carotid media-adventitia and lumen-intima boundary segmentation from three-dimensional ultrasound images, Med. Phys., 46 (2019), 3180–3193. https://doi.org/10.1002/mp.13581 doi: 10.1002/mp.13581
[18]	R. Zhou, F. Guo, M. R. Azarpazhooh, J. D. Spence, E. Ukwatta, M. Ding, A. Fenster, A voxel-based fully convolution network and continuous max-flow for carotid vessel-wall-volume segmentation from 3D ultrasound images, IEEE Trans. Med. Imaging, 39 (2020), 2844–2855. https://doi.org/10.1109/tmi.2020.2975231 doi: 10.1109/tmi.2020.2975231
[19]	R. Zhou, F. Guo, M. R. Azarpazhooh, S. Hashemi, X. Cheng, J. D. Spence, et al., Deep learning-based measurement of total plaque area in b-mode ultrasound images, IEEE J Biomed. Health Inform., 25 (2021), 2967–2977. https://doi.org/10.1109/jbhi.2021.3060163 doi: 10.1109/jbhi.2021.3060163
[20]	C. Doersch, A. Zisserman, Multi-task self-supervised visual learning, in 2017 IEEE International Conference on Computer Vision(ICCV), (2017), 2070–2079. https://doi.org/10.1109/ICCV.2017.226
[21]	Z. Ren, Y. J. Lee, Cross-domain self-supervised multi-task feature learning using synthetic imagery, in 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, (2018). http://dx.doi.org/10.1109/CVPR.2018.00086
[22]	X. Zheng, Y. Wang, G. Wang, J. Liu, Fast and robust segmentation of white blood cell images by self-supervised learning, Micron, 107 (2018), 55–71. https://doi.org/10.1016/j.micron.2018.01.010 doi: 10.1016/j.micron.2018.01.010
[23]	L. Chen, P. Bentley, K. Mori, K. Misawa, M. Fujiwara, D. Rueckert, Self-supervised learning for medical image analysis using image context restoration, Med. Image Anal., 58 (2019), 101539. https://doi.org/10.1016/j.media.2019.101539 doi: 10.1016/j.media.2019.101539
[24]	X. Zhuang, Y. Li, Y. Hu, K. Ma, Y. Yang, Y. Zheng, Self-supervised feature learning for 3D medical images by playing a rubik's cube, in International Conference on Medical Image Computing and Computer-Assisted Intervention, (2019), 420–428. https://doi.org/10.1007/978-3-030-32251-9_46
[25]	Q. Lu, Y. Li, C. Ye, White matter tract segmentation with self-supervised learning, in International Conference on Medical Image Computing and Computer-Assisted Intervention, (2020), 270–279. https://doi.org/10.1007/978-3-030-59728-3_27
[26]	J. D. Spence, M. Eliasziw, M. Dicicco, D. G. Hackam, R. Galil, T. Lohmann, Carotid plaque area: A tool for targeting and evaluating vascular preventive therapy, Stroke, 33 (2002), 2916–2922. https://doi.org/10.1161/01.str.0000042207.16156.b9 doi: 10.1161/01.str.0000042207.16156.b9
[27]	S. H. Johnsen, E. B. Mathiesen, O. Joakimsen, E. Stensland, T. Wilsgaard, M.L. Løchen, et al., Carotid atherosclerosis is a stronger predictor of myocardial infarction in women than in men: A 6-year follow-up study of 6226 persons: the tromsø study, Stroke, 38 (2007), 2873–2880. https://doi.org/10.1161/strokeaha.107.487264 doi: 10.1161/strokeaha.107.487264
[28]	M. Noroozi, P. Favaro, Unsupervised learning of visual representations by solving jigsaw puzzles, in European Conference on Computer Vision, (2016), 69–84. https://doi.org/10.1007/978-3-319-46466-4_5
[29]	O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, in International Conference on Medical Image Computing and Computer-Assisted Intervention, (2015), 234–241. https://doi.org/10.1007/978-3-319-24574-4_28
[30]	Z. Zhou, M. M. R. Siddiquee, N. Tajbakhsh, J. Liang, Unet++: Redesigning skip connections to exploit multiscale features in image segmentation, IEEE Trans. Med. Imaging, 39 (2020), 1856–1867. https://doi.org/10.1109/tmi.2019.2959609 doi: 10.1109/tmi.2019.2959609
[31]	R. Zhou, M. R. Azarpazhooh, J. D. Spence, S. Hashemi, W. Ma, X. Cheng, et al., Deep learning-based carotid plaque segmentation from b-mode ultrasound images, Ultrasound in Medicine & Biology, 47 (2021), 2723–2733. https://doi.org/10.1016/j.ultrasmedbio.2021.05.023 doi: 10.1016/j.ultrasmedbio.2021.05.023
[32]	K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, Proceedings of the IEEE conference on computer vision and pattern recognition, 2016,770–778. https://doi.org/10.1109/CVPR.2016.90
[33]	T. Heimann, B. V. Ginneken, M. A Styner, Y. Arzhaeva, V. Aurich, C. Bauer, et al., Comparison and evaluation of methods for liver segmentation from CT datasets, IEEE transactions on medical imaging, 28 (2009), 1251–1265. https://doi.org/10.1109/TMI.2009.2013851 doi: 10.1109/TMI.2009.2013851

This article has been cited by:

1.	Guoning Si, Wenkai Li, Hanjing Lu, Zhuo Zhang, Xuping Zhang, Vibration Modeling and Analysis of a Flexible 3-PRR Planar Parallel Manipulator Based on Transfer Matrix Method for Multibody System, 2023, 11, 2075-1702, 505, 10.3390/machines11050505
2.	Yunxia Wei, Yuanfei Zhang, Bin Hang, Construction and management of smart campus: Anti-disturbance control of flexible manipulator based on PDE modeling, 2023, 20, 1551-0018, 14327, 10.3934/mbe.2023641
3.	Decio de M. Rinaldi, Tarcisio A. Hess-Coelho, Modular modeling of parallel mechanisms with actuation redundancy, 2025, 47, 1678-5878, 10.1007/s40430-025-05671-1

Reader Comments

Your name:*

Email:*
© 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)