1 Introduction
Many of the objects we handle every day are highly deformable and with prevalent plastic behavior. However, humans manipulate those objects naturally, with high dexterity, and without any particular issue. In facts, manipulating these deformable objects has a wide variety of uses in domestic facilities and healthcare, such as robotic surgery [17], assistive dressing, garment sorting or folding clothing [25], [31]. Moreover, it is involved in manufacturing, aerospace [35], automotive, and electromechanical industries generally [18],[7].
On the contrary, the manipulation of deformable objects is still a challenging activity for robots. This is the main reason why many assembly procedures involving such deformable objects are still performed manually. One of the main reasons why robots have such limitations in deformable object manipulation is due to their complex behaviors, unpredictable initial configuration, and limited capability in measuring their state.
A thoughtful survey can be found in [32] that focuses on deformable object manipulation by robots in industrial and domestic applications. Actually, the dynamics of deformable objects is complex and nonlinear. Therefore, the state estimation of a deformable object is challenging, and the forward prediction is expensive. Robust and effective methods to manipulate deformable objects and predict their behavior remain extremely difficult to build, despite the several applications and attempts made by the robotics community [32].
In this work, the numerical integration of Deformable Linear Objects (DLOs), such as strings, cables, electric wires, catheters, ropes, and so on, is addressed, since efficient solutions to this problem will enable the implementation of robotized solution in many relevant subfields of the large and diverse industrial manufacturing. DLO manipulation is fundamental in automotive manufacturing, e.g. for wiring harness preparation and electrical cable installation inside the vehicle structure, as well as in medical surgery, e.g. in suturing in which a flexible wire in a straight configuration needs to go to a knot [26], and have a vital role in other fields such as architecture and power distribution. The relevant work on DLO manipulation can be split into three overlapping categories: modeling, simulation, and planning. Various manipulation tasks such as knot tying [16], rope untangling [22], string insertion [39], and shape manipulation [30] can be executed on DLOs.
A model of DLO deformation/flexibility is needed in order to predict its behavior successfully. This model implies that both the geometry and the mechanical behavior of the concerned parts can be represented accurately. Once the model is properly defined, a computationally efficient way is needed to evaluate it over time or to solve queries of motion planning. Virtual prototyping is used to reduce development costs and to boost consistency. Also, it enables the early detection of possible problems. It allows for studying the efficiency of assembly. The most distinctive characteristic of DLOs is the variations in shape that they undergo following the influence of forces and environmental constraints. Such shape variations may be the purpose of the planning itself, for example if they are caused by obstacles along the path [19]. Overall, DLO modeling is a major and complex challenge, with a broad range of applications. The PDE model can be established [10, 12, 11, 15, 14], and the associated structurepreserving method can be developed for the dynamic problems of the DLOs. Various research in the literature focused on modeling and manipulating DLOs for various purposes and many models and strategies were created. There are several various approaches for physically modeling DLOs such as Massspring [23], Multibody [34], Elastic rod [21], Dynamic spline [38], Finite element [5], and other models. The methodology, advantages and disadvantages of each model method are discussed in detail in [24]. The dynamic spline is one of these methods which provide a good theoretical basis, continuous model, higher authenticity. A geometrically reliable model of DLOs is generated and adopted in [36] to execute numerical simulations on the motion of the object under gravity and during environmental interaction. For interaction simulation, a linearized splinebased DLO model named quasi dynamic splines is created in [37].
Owing to the inherent tradeoff between precision and realtime capability, it is hard to deal with DLOs due to their complex model, resulting then in timeconsuming integration processes to predict their behavior. To solve numerically the differential equations, many types of integrators can be found in literature, such as RungeKutta, Euler, Taylor type integrators and many others. Despite their general applicability, the main problems associated with the use of these integrators are that they are often inefficient, and prone to instability. In [29] a model of DLO is implemented using a multivariate dynamic spline. The integration process is achieved using the traditional type RungeKutta method. This integration process is very time consuming and provide unstable results with many stuck situations. Also, it is not practical for a longtime experiment. Therefore, it is important to find another integration method to reduce simulation time.
The dynamics of a mechanical system, such as a DLO, can be easily represented by means of Hamiltonian equations. It is also known that a symplectic transformation is the solution of a Hamiltonian system [20]. The numerical solutions obtained by several numerical techniques, such as the RungeKutta method and the primitive Euler scheme, are not energypreserving map, resulting in spuriously damped or excited behaviors.
On the other hand, the symplectic integrator proposed in [2] exploits the energypreserving symplectic transformation to avoid spuriously damped or excited solutions. The symplectic integrator scheme has been widely applied to the calculations of the longterm evolution of chaotic Hamiltonian systems [33], and [8]. In the actions of a Hamiltonian system, the symplectic integrator does not generate a secular truncation error. The number of force function evaluations of the fourthorder symplectic integrator is smaller than the ones of the RungeKutta integrator of the same order. The energy conservation and the longtime stability for the symplectic scheme are investigated and verified in [2], [13], [11], and [9]. There are many merits for symplectic integrator which are discussed and outlined in deeper detail in [3], [4], and [27]. For these advantages, the use of the symplectic integrator is investigated in this paper for DLO simulation instead of the conventional integrators. To this end, a comparison between the results obtained from this symplectic integrator and other types of integrators will be executed.
The main objective of this paper is to compare the properties and results obtained by the symplectic integrator with the ones obtained with other numerical methods. This activity is justified by the fact that, at the best of our knowledge, no data are available in literature about the performances of numerical integration algorithms for DLOs. For the DLO modeling, multivariate dynamic splines are used in this paper. A comparison between the symplectic, RungeKutta, and Zhai methods is performed in two cases. In the first case, the effect of gravitational and internal forces only are considered, while in the second case the application of an external force in addition to the previously mentioned forces is considered.
The remainder of this article is arranged as follows. Sec. 2 outlines the key features of the mathematical model of the DLO. The conversion from Lagrangian into Hamiltonian and the symplectic integrator are discussed in Sec 3. Moreover, Sec. 4 introduces the preliminary simulation results and comparisons. Finally, conclusions and future work are draft out in Sec. 5.
2 Dynamic Model of the DLO
A thirdorder spline basis can effectively represent the shape of the DLO. It is a function of a free coordinate . This coordinate represents the position on the DLO which is equal to zero at the first endpoint. Also, equals at the opposite endpoint, provided that is the length of the DLO. This can be written mathematically as:
(1) 
where is the fourthdimensional configuration functional space of the DLO. It includes the three linear coordinates of the DLO position at point , and the DLO’s axial twisting . Also, is the spline polynomial basis employed to describe the shape of the DLO. Moreover, are properly chosen coefficients, typically called control points, used to correctly interpolate the shape of the DLO through the basis functions. is the number of control points.
For a variety of reasons, this DLO mathematical model is highly successful. Firstly, the spatial derivatives calculation is straightforward, i.e.
(2) 
where is the derivative of . In fact, it can be defined by the same coefficients and easytocompute derivatives of . Second, the spline basis proprieties guarantee that the DLO model curvature, which represents the DLO’s physical behavior [1], is minimized. Third, this model enables the representation of a generalized nonlinear function with smoothness characteristics as a linear combination of the nonlinear function basis , which relies on the free variable only, by the linear coefficients .
Referring to the system’s Lagrange equations, the DLO’s dynamic model can be expressed as a function of the control points as:
(3) 
where is the total kinetic energy, is the external force that is applied to the control point, and is the total potential energy acting on the DLO. This potential energy is generated due to the influence of gravity, stretching, torsion, and bending effects on the DLO.
2.1 Kinetic Energy of the DLO
Owing to the translational motion of the control points and the rotational motion of the cross sections, the total kinetic energy is generated. In order to express this total kinetic energy, a function of the control points can be used as :
(4) 
where is the DLO’s generalized density matrix, is the displacement element and can be computed from , denotes the linear density, and
denotes the polar moment of inertia.
According to the procedure explained in [28], we can drive:
(5) 
Then, by expanding this definition to the overall system, the total inertial forces of the DLO can be written as , in which denotes the inertia matrix of the DLO and
is a vector that represents the control point accelerations.
2.2 Potential Energy of the DLO
The total potential energy is made up of the gravitational effect and the strain energy of the DLO due to stretching, torsion and bending. Although the derived energy of gravity is quite simple, the concept of strain energy plays a crucial role in the modeling of DLO. By introducing the strain vector that includes to represent the stretching term, to denotes the torsional term, and to denotes the bending term such that:
(6)  
where and are the strain’s linear and the torsional component, respectively. The strain energy can therefore be written as:
(7) 
where
is the DLO’s plastic strain, which enables to consider the plasticity of the material. is called the residual strain that equals . Also, is the element stiffness matrix, is the crosssection diameter of the DLO. Moreover, and denote the Young’s modulus and the shear modulus of the DLO’s material, respectively.
It is possible to write the rightside term of Eq. (3) as:
(8) 
where represents the elastic forces because of the DLO deflection.
2.3 Dynamic Model of the DLO
To write the overall dynamic model of the DLO, the Eqs. (3), (5), and (8) are extended to the whole control points
(9) 
where refers to the vector that contains all external forces, including gravity, while represents the vector that contains all elastic forces.
At every simulation step, the system can be solved by utilizing a simple LU decomposition. Moreover, at each time step, the accelerations are integrated to get the control points velocities and positions. Several integration methods can be used, but some problems may arise, including numerical instability, large time required for the integration process and the unsuitability for longtime predictions. These problems can be mitigated by using the symplectic integrator, as will be shown in the following. The next section will discuss this symplectic integrator and its requirements in more details.
3 Symplectic Integrator
To solve numerically the differential equations, we decided to use the symplectic integrator because of its advantages with respect to other methods as will be shown in the following. These subsections will discuss in detail calculating the Hamiltonian from Lagrangian and the construction of the symplectic integrator.
3.1 Conversion from a Lagrangian to a Hamiltonian
Symplectic integrator is widely used in nonlinear dynamics. It is designed as the numerical solution for Hamilton’s equations, which provided as:
(10) 
where denotes the momentum coordinates, refers to the position coordinates, and is the Hamiltonian which can be found from:
(11) 
where and denote the kinetic and potential energy, respectively. The set of position and momentum coordinates are named canonical coordinates. In our case, the first step to use the symplectic integrator is to make a conversion from Lagrangian into Hamiltonian [6]. This conversion is achieved according to the following procedure.
Provided a Lagrangian as a function of the generalized coordinates and generalized velocities , where , the Hamiltonian can be calculated according to the following steps:

By differentiating the Lagrangian with respect to the generalized velocities , the momenta are determined:
(12) 
By inverting the expressions in the former step, the velocities are formulated in terms of the momenta .
(13) where

Using the Lagrangian relation ( ), we can conclude that:
(14) The velocities are then substituted from Equ. (13),
(15) 
The Hamiltonian is determined by employing the typical definition of as the Legendre transformation of :
(16) Substitution for from Equ. (13) and the Lagrangian from Equ. (15) into Equ. (16), will lead to:
(17) The last equation is the require one as it is equivalent to the Hamiltonian equation that stated in Equ. (11).
3.2 Symplectic Integrator of the fourthorder
Forest [3] and Neri [27] introduced the generalized form of the fourthorder symplectic integrator as:
(18) 
where is the time stepsize, and are the initial values, and and are the numerical solution after . Moreover, and are numerical coefficients which can be determined uniquely from:
(19) 
The reader can refer to the reference [4] for the derivation and the values of the numerical coefficients.
The number of force function evaluations , that is the most timeexpensive operation in the integration process, is three rather than four since . On the other hand, four force function evaluations are required in the conventional fourthorder RungeKutta integrator. Thus, the CPU time can be decreased by approximately percent with the symplectic integrator with respect to the fourthorder RungeKutta method.
Other than the energypreserving property, a remarkable characteristic of the symplectic integrator is that the aggregation of the truncation error in the total energy has not a secular component and the positional errors rise in proportional to the first order of time [20]. On the other hand, conventional integrators generate a secular energy error. Moreover, the error due to position discretization increase with the square of time. The symplectic integrator is thus well suited to the longtime study of a dynamic system. Generally, in the actions of a Hamiltonian system, the symplectic integrator does not generate a secular truncation error [20].
4 Simulation Results
Numerical simulations have been carried out in MATLAB for the assessment of the mathematical framework discussed in the preceding sections. The issue that is aimed in this section is determining the DLO state evolution, i.e. the control points evolution over time, while influenced by a internal elastic and inertial forces as well as external forces like gravitational force or contacts forces generated by obstacles. In the case of contacts, we can reshape the system as a constrained dynamic system in which the Lagrange approach for constraints is adopted. Just add some constraints to the solution of the dynamic system. This will not change the problem, but just increase the number of dynamic equations that existed in the model. In [36], the solution of dynamic equations representing the DLO dynamics including contacts is reported. The same approach can be applied also to the symplectic integrator by changing the dynamic equations accordingly. The simulations have been performed using the dynamic model in Eq. (9). It is worth mentioning that all simulations of this work have been performed on Ubuntu 18.04.5 LTS operating system, processor intel core i53210M CPU@2.50GHz x 4, RAM 8GiB.
In our simulations, DLO’s material is assumed to be aluminum, and its length is 2 meters. The DLO has a circular crosssection with a diameter D that equals 2 millimeters. The number of control points is selected as = 9. Also, the number of sample points along the DLO is chosen to be =101. These values are chosen to be sufficiently low but provide a good interpolation capability. During these simulations, the DLO is supposed to be placed initially on a straight position along the axis, and the two extreme endpoints are constrained to move along the axis and attached by springs with stiffness to their initial position. In other words, one endpoint of the DLO is located at , while the opposite endpoint is located at . Two simulation scenarios will be used in this section. In the first scenario, the DLO is affected by the internal forces and the gravitational force only as shown in Fig. 1a. While in the second scenario, the DLO is affected by the previously mentioned forces in addition to an external sinusoidal force which is applied to the DLO center in the Ydirection as illustrated in Fig. 1b.
Fig. 2, and 3 illustrate the computed solutions accomplished utilizing the splinebased model using symplectic integrator where the DLO is affected by the gravity and the internal forces only without any other external forces. Fig. 2 presents the component of each control point, while Fig. 3 shows the component of each control point. This simulation is achieved for 10 seconds using 2 milliseconds as a time stepsize . This time stepsize is the largest value that preserves the stability of the integration method. It is worth mentioning that, during the whole period, the component for each control point is equal to zero as the DLO is positioned along the axis and the external force is not applied yet. In Fig. 3, due to the symmetry of the DLO, each two control points’ trajectories are located above each other in one curve except the middle point. A video that shows this numerical simulation exists in the link ^{1}^{1}1https://www.dropbox.com/s/4tb3ohyod5v71yj/videoforspline.avi?dl=0. Fig. 4 shows the sequence of trajectories that illustrate the DLO motion starting from the initial position according to the first scenario, i.e. the DLO is affected by the internal forces and the gravitational force only without any other external forces.
The simulation is performed again in the case of applying an external force on a DLO point. We assumed that the external force takes the shape of a sinusoidal trajectory. This external force is applied to the center of the DLO in the Yaxis direction. Fig. 5, 6, and 7 show the , , and components, respectively. The solution discussed here assumes the knowledge of the external forces (the input). The computation will give the control points (the output) that describe the DLO motion during the simulation interval. A video that illustrates this numerical simulation, in case of applying an external force, is found in the link ^{2}^{2}2https://www.dropbox.com/s/b1tc8v1k627o33r/videoforspline.avi?dl=0. These results prove the ability of the splinebased dynamic model and the symplectic integrator to represent the state evolution of the DLO accurately. Fig. 8 shows the sequence of trajectories that illustrate the DLO motion starting from the initial position according to the second scenario, i.e. the DLO is affected by the internal forces, the gravitational force, and an external sinusoidal force applied to the DLO center in the Ydirection.
The DLO model simulation time of the symplectic integrator is compared with the ones of both the classical RungeKutta integrator and Zhai integrator. Zhai integrator is a new simple fast explicit time integration method, the reader can refer to [40] for the related details and advantages of Zhai integrator. In MATLAB/Simulink, the three integrators are implemented by using the same splinebased model defined in Sec. 2. In the case the DLO is affected by internal forces and gravity only, Table 1 reports the comparison between the three integrators’ simulation time in seconds for a certain period (10 sec) with 2 milliseconds as a time stepsize . The results show that the symplectic integrator is faster than both the RungeKutta and Zhai integration methods. It is worth mentioning that the Zhai method needs optimization for two parameters [40] . However, we made a manual tunning to these parameters. This may be the reason for the notable slowness of Zhai compared to the RungeKutta method.
Table 2 presents the comparison between the three integrators’ computation time in seconds for a simulation of 10 seconds in the case where the DLO is affected by the external force in addition to the previously mentioned forces. The results show that the symplectic integrator is able to solve the problem while, on the other hand, the RungeKutta and Zhai integrators have been stuck and could not compute any solution. It is worth mentioning that, in this work, we did not investigate the reason why other methods are unstable.
The simulation implemented in MATLAB, using the fastest integrator, takes an average time of 470 seconds for the simulation of 10 seconds like the one reported in Fig. 2 and 3. This long execution time occurs due to two reasons. The first reason is the specifications of the processing unit or the PC. The second reason is the implementation on MATLAB which is not really efficient. To avoid this reason, another simulation is repeated using C++ language. In these simulations and in the related comparisons, the Zhai integrator is excluded as it is the slowest and it looks not really efficient in solving the DLO dynamics. Hence, the current comparison will be focused on the symplectic and the RungeKutta integration methods and their implementation in C++. Tables 3 and 4 introduce the comparison between the two integrators’ execution time in the case of internal forces only and sinusoidal external force, respectively. From Tables 3 and 4, it can be concluded that the symplectic integrator gives the fastest execution time. Also, the C++ execution time is reduced to be onethird of the execution time of MATLAB results.
To check the effect of changes in the DLO model granularity, and on the execution time, simulations on both MATLAB and C++ are performed. Tables 5 and 6 present MATLAB results of the execution time by changing and for both cases of internal forces only and sinusoidal external force, respectively. Each cell in the two tables introduces the execution time in seconds for a simulation period of 10 seconds. Likewise, Table 7, and 8 illustrate C++ results of execution time with changing and . From these tables, we conclude that C++ is faster about three times more than MATLAB implementations. So, C++ implementation with symplectic integrator will be better for the practical implementation of the simulation system. Moreover, reducing and/or will require less execution time, but it will reduce interpolation precision in representing the DLO. It is worth mentioning that the time step size is chosen to be unique and equal to 0.8 milliseconds for all of these comparisons. As and are increased, more calculations are required and hence the integration time increases. The same step size is selected for all comparisons for fairness.
The model with the highest resolution ( and ) is considered as the reference model. Each model resolution considered in this work is compared to the reference model and the error is computed. For each resolution, the DLO position is measured at 10 fixed points along the DLO. Moreover, the error over time is computed by considering 20 equidistant time intervals. Fig. 9 shows the evolution of this error over the time and the DLO length. It is noticed that the maximum error occurred at the center of the DLO and the error decreases over time. The evolution of the error with different values for and is presented in Fig. 10. As increases, the error decreases and vice versa, while has no noticeable effect on the error. Moreover, as the number of control points increases, the error decreases over time and with respect to the position along the DLO as indicated in Fig. 11 and 12, respectively.
Simulation time  Symplectic  RungeKutta  Zhai 

0.1  6.66  9.43  34.85 
0.2  11.39  16.42  67.57 
0.3  16.14  22.43  96.69 
0.4  20.86  28.41  120.45 
0.5  25.60  34.62  145.42 
0.6  30.42  41.29  166.22 
0.7  35.17  48.62  186.46 
0.8  39.92  56.07  206.70 
0.9  44.61  62.59  226.81 
1  49.30  68.89  247.39 
2  97.18  140.00  478.94 
3  144.21  210.41  685.46 
4  191.20  278.22  904.89 
5  238.24  347.53  1127.37 
6  285.31  418.16  1348.17 
7  332.29  487.72  1542.90 
8  379.23  555.20  1733.71 
9  426.17  624.36  1970.40 
10  473.18  697.04  2204.93 
Simulation time  Symplectic  RungeKutta  Zhai  

0.1  6.83 



0.2  11.75  
0.3  16.70  
0.4  21.64  
0.5  26.60  
0.6  31.56  
0.7  36.55  
0.8  41.48  
0.9  46.44  
1  51.40  
2  101.14  
3  150.73  
4  200.30  
5  249.89  
6  298.60  
7  346.63  
8  394.63  
9  442.50  
10  490.59 
Simulation time  Symplectic  RungeKutta 

0.1  2.04  5.88 
0.2  3.85  10.44 
0.3  5.65  13.90 
0.4  7.55  16.75 
0.5  9.38  19.93 
0.6  11.17  23.09 
0.7  12.95  26.43 
0.8  14.73  29.82 
0.9  16.56  33.00 
1  18.35  36.39 
2  36.57  69.57 
3  54.67  102.85 
4  72.88  136.00 
5  90.99  169.06 
6  109.11  201.95 
7  127.30  235.07 
8  145.86  268.21 
9  164.16  301.24 
10  182.46  334.40 
Simulation time  Symplectic  RungeKutta 

0.1  1.98  5.45 
0.2  3.98  9.85 
0.3  5.92  13.63 
0.4  7.93  16.88 
0.5  9.99  20.12 
0.6  11.91  23.42 
0.7  13.97  26.81 
0.8  15.99  30.52 
0.9  17.90  34.53 
1  19.90  38.17 
2  40.10  76.95 
3  60.16  115.53 
4  80.41  154.16 
5  101.47  193.16 
6  123.51  231.97 
7  146.46  270.46 
8  167.34  309.58 
9  189.53  349.07 
10  210.46  388.04 
51  81  101  121  151  251  

5  488.28  688.22  825.99  960.52  1288.93  2094.12 
7  588.53  811.56  998.6  1201.32  1419.81  2357.17 
9  713.17  967.77  1299.68  1421.18  1634.21  2787.05 
11  770.00  1129.01  1344.34  1691.19  2033.16  3171.8 
13  908.49  1332.11  1566.67  1893.44  2304.9  3909.13 
19  1331.77  1714.65  2348.26  2479.57  3104.75  5090.28 
51  81  101  121  151  251  

5  489.98  677.93  888.11  1004.88  1254.07  2015.19 
7  538.59  881.56  1007.06  1214.88  1554.83  2379.88 
9  643.93  970.64  1124.15  1405.01  1675.14  2970.31 
11  717.9  1101.62  1385.67  1692.81  2006.02  3140.75 
13  874.72  1291.78  1568.77  1886.01  2313.28  3903.67 
19  1386.62  1826.59  2102.17  2667.17  2990.85  5035.5 
51  81  101  121  151  251  

5  155.28  232.5  284.52  339.25  416.16  673.25 
7  197.56  295.17  363.86  435.05  529.7  861.73 
9  245.61  369.37  455.4  523.42  652.07  1065.55 
11  286.957  436.69  534.16  623.1  765.75  1271.19 
13  319.786  489.52  601.8  713.24  887.17  1437.82 
19  450.95  690.64  840.81  990.89  1236.14  2021.42 
51  81  101  121  151  251  

5  154.18  231.83  283.43  335.97  414.29  680.41 
7  195.2  295.25  361.43  427.65  530.54  863.73 
9  238.25  358.64  442.07  521.7  645.43  1060.26 
11  278.72  426.39  523.94  617.38  764.78  1253.72 
13  319.12  488.51  600.33  709.03  889.24  1436.29 
19  442.82  677.01  835.82  988.3  1228.88  2023.9 
5 Conclusions
This paper reports a comparison among integration methods to solve the dynamic model of Deformable Linear Objects. The adopted DLO model is based on the multivariate dynamic spline and employs the symplectic integrator. This symplectic integrator received particular attention due to its advantages over alternative methods in solving the dynamic equations of Hamiltonian systems. Simulations on MATLAB and using C++ code are performed to compare the performance of the symplectic integrator against the RungeKutta and Zhai integrators. The results prove that the symplectic integrator is the fastest and the most stable integrator among the considered ones. Moreover, the proposed comparison shows that C++ is about three times faster than the corresponding MATLAB implementation. Hence, a splinebased model with a symplectic integrator looks to be the best choice for practical implementation. Finally, comparisons are performed to check the effects of changing the model granularity. These results show that reducing the control points and sample points will reduce the execution time, but it will reduce the interpolation precision in representing the DLO behavior.
In future work, practical experiments will be used to validate the splinebased model and the results obtained with the symplectic integrator. Moreover, this model will be exploited to optimize the DLO manipulation process. Also, an extension to the method will be applied to multibranch DLOs, like the preassembled wiring harness utilized in the aerospace and automotive industries.
acknowledgment
This work was supported by the European Commissions Horizon 2020
Framework Programme with the project REMODEL  Robotic technologies
for the manipulation of complex deformable linear objects  under
grant agreement No 870133.
References
 [1] Carl De Boor. A practical guide to splines; rev. ed., ser. applied mathematical sciences. 2001.
 [2] Kang Feng. On difference schemes and symplectic geometry. In Proceedings of the 5th international symposium on differential geometry and differential equations, 1984.
 [3] Etienne Forest. Canonical integrators as tracking codes (or how to integrate perturbation theory with tracking). Technical report, Lawrence Berkeley Lab., 1987.
 [4] Etienne Forest and Ronald D Ruth. Fourthorder symplectic integration. Physica D: Nonlinear Phenomena, 43(1):105–117, 1990.
 [5] Leopoldo Greco and Massimo Cuomo. Bspline interpolation of kirchhofflove space rods. Computer Methods in Applied Mechanics and Engineering, 256:251–269, 2013.
 [6] Patrick Hamill. A student’s guide to Lagrangians and Hamiltonians. Cambridge University Press, 2014.
 [7] Tomas Hermansson, Robert Bohlin, Johan S Carlson, and Rikard Söderberg. Automatic assembly path planning for wiring harness installations. Journal of manufacturing systems, 32(3):417–422, 2013.
 [8] Jialin Hong, Chuying Huang, and Xu Wang. Symplectic runge–kutta methods for hamiltonian systems driven by gaussian rough paths. Applied Numerical Mathematics, 129:120–136, 2018.
 [9] Weipeng Hu, Zichen Deng, Songmei Han, and Wenrong Zhang. Generalized multisymplectic integrators for a class of hamiltonian nonlinear wave pdes. Journal of Computational Physics, 235:394–406, 2013.
 [10] Weipeng Hu, Zhen Wang, Yunping Zhao, and Zichen Deng. Symmetry breaking of infinitedimensional dynamic system. Applied Mathematics Letters, 103:106207, 2020.
 [11] Weipeng Hu, Mengbo Xu, Jiangrui Song, Qiang Gao, and Zichen Deng. Coupling dynamic behaviors of flexible stretching hubbeam system. Mechanical Systems and Signal Processing, 151:107389, 2021.
 [12] Weipeng Hu, Juan Ye, and Zichen Deng. Internal resonance of a flexible beam in a spatial tethered system. Journal of Sound and Vibration, 475:115286, 2020.
 [13] Weipeng Hu, Tingting Yin, Wei Zheng, and Zichen Deng. Symplectic analysis on orbitattitude coupling dynamic problem of spatial rigid rod. Journal of Vibration and Control, 26(1718):1614–1624, 2020.
 [14] Weipeng Hu, Lingjun Yu, and Zichen Deng. Minimum control energy of spatial beam with assumed attitude adjustment target. Acta Mechanica Solida Sinica, 33(1):51–60, 2020.
 [15] Weipeng Hu, Chuanzeng Zhang, and Zichen Deng. Vibration and elastic wave propagation in spatial flexible damping panel attached to four special springs. Communications in Nonlinear Science and Numerical Simulation, 84:105199, 2020.
 [16] Sandy H Huang, Jia Pan, George Mulcaire, and Pieter Abbeel. Leveraging appearance priors in nonrigid registration, with application to manipulation of deformable objects. In 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 878–885. IEEE, 2015.
 [17] Jagadeesan Jayender, Rajnikant V Patel, and Suwas Nikumb. Robotassisted active catheter insertion: Algorithms and experiments. The International Journal of Robotics Research, 28(9):1101–1117, 2009.
 [18] Xin Jiang, Kyongmo Koo, Kohei Kikuchi, Atsushi Konno, and Masaru Uchiyama. Robotized assembly of a wire harness in a car production line. Advanced Robotics, 25(34):473–489, 2011.
 [19] P Jiménez. Survey on modelbased manipulation planning of deformable objects. Robotics and computerintegrated manufacturing, 28(2):154–163, 2012.
 [20] Hiroshi Kinoshita, Haruo Yoshida, and Hiroshi Nakai. Symplectic integrators and their application to dynamical astronomy. Celestial Mechanics and Dynamical Astronomy, 50:59–71, 1991.
 [21] Joachim Linn and Klaus Dreßler. Discrete cosserat rod models based on the difference geometry of framed curves for interactive simulation of flexible cables. In Math for the Digital Factory, pages 289–319. Springer, 2017.
 [22] Wen Hao Lui and Ashutosh Saxena. Tangled: Learning to untangle ropes with rgbd perception. In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 837–844. IEEE, 2013.
 [23] Naijing Lv, Jianhua Liu, Xiaoyu Ding, Jiashun Liu, Haili Lin, and Jiangtao Ma. Physically based realtime interactive assembly simulation of cable harness. Journal of Manufacturing Systems, 43:385–399, 2017.
 [24] Naijing Lv, Jianhua Liu, Huanxiong Xia, Jiangtao Ma, and Xiaodong Yang. A review of techniques for modeling flexible cables. ComputerAided Design, 122:102826, 2020.
 [25] Rene J Moreno Masey, John O Gray, Tony J Dodd, and Darwin G Caldwell. Guidelines for the design of lowcost robots for the food industry. Industrial Robot: An International Journal, 2010.
 [26] Mark Moll and Lydia E Kavraki. Path planning for deformable linear objects. IEEE Transactions on Robotics, 22(4):625–636, 2006.
 [27] Filippo Neri. Lie algebras and canonical integration. Dept. of Physics, University of Maryland, 1987.
 [28] Olivier Nocent and Yannick Remion. Continuous deformation energy for dynamic material splines subject to finite displacements. In Computer Animation and Simulation 2001, pages 87–97. Springer, 2001.
 [29] Gianluca Palli. Modelbased manipulation of deformable linear objects by multivariate dynamic splines. In 2020 IEEE Conference on Industrial Cyberphysical Systems (ICPS), volume 1, pages 520–525. IEEE, 2020.
 [30] Matthias Rambow, Thomas Schauß, Martin Buss, and Sandra Hirche. Autonomous manipulation of deformable objects based on teleoperated demonstrations. In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 2809–2814. IEEE, 2012.
 [31] Arnau Ramisa, Guillem Alenya, Francesc MorenoNoguer, and Carme Torras. Using depth and appearance features for informed robot grasping of highly wrinkled clothes. In 2012 IEEE International Conference on Robotics and Automation, pages 1703–1708. IEEE, 2012.
 [32] Jose Sanchez, JuanAntonio Corrales, BelhassenChedli Bouzgarrou, and Youcef Mezouar. Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey. The International Journal of Robotics Research, 37(7):688–716, 2018.
 [33] JM9728120655 SanzSerna. Rungekutta schemes for hamiltonian systems. BIT Numerical Mathematics, 28(4):877–883, 1988.
 [34] Martin Servin and Claude Lacoursiere. Rigid body cable for virtual environments. IEEE Transactions on Visualization and Computer Graphics, 14(4):783–796, 2008.
 [35] Ankit Shah, Lotta Blumberg, and Julie Shah. Planning for manipulation of interlinked deformable linear objects with applications to aircraft assembly. IEEE Transactions on Automation Science and Engineering, 15(4):1823–1838, 2018.
 [36] Adrien Theetten, Laurent Grisoni, Claude Andriot, and Brian Barsky. Geometrically exact dynamic splines. ComputerAided Design, 40(1):35–48, 2008.
 [37] Adrien Theetten, Laurent Grisoni, Christian Duriez, and Xavier Merlhiot. Quasidynamic splines. In Proceedings of the 2007 ACM symposium on Solid and physical modeling, pages 409–414, 2007.
 [38] Pier Paolo Valentini and Ettore Pennestrì. Modeling elastic beams using dynamic splines. Multibody system dynamics, 25(3):271–284, 2011.
 [39] Weifu Wang, Dmitry Berenson, and Devin Balkcom. An online method for tighttolerance insertion tasks for string and rope. In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 2488–2495. IEEE, 2015.
 [40] Wanming Zhai. Numerical method and computer simulation for analysis of vehicle–track coupled dynamics. In Vehicle–track coupled dynamics, pages 203–229. Springer, 2020.