Synchronized Motion Proﬁles for Inverse-Dynamics-Based Online Control of Three Inextensible Segments of Trunk-Type Robot Actuators

: This study proposes a novel method for the positioning and spatial orientation control of three inextensible segments of trunk-type robots. The suggested algorithm imposes a soft constraint assumption for the end-effector’s endpoint and a mandatory constraint on its direction. Simultane-ously, the algorithm by-design enforces nonholonomic features on the robot segments in the form of arcs. An approximate robot spine curve is the key to the ﬁnal robot state conﬁguration based on the given conditions. The numeric simulation showed acceptable (less than 1 s) performance for single-core processing tasks. The parametric method ﬁnds the best proximate robot state solution and represents the gray box model in addition to existing learning or black-box inverse dynamics approaches. This study also shows that a multiple inverse kinematics answer constructs a single inverse dynamics solution that deﬁnes the robot actuators’ motion proﬁles, synchronized in time. Finally, this text presents rotational expressions and their outlines for controlling the manipulator’s tendons. data curation, M.M. and L.Z.; writing—original draft preparation, M.M.; writing—review and editing, R.U. and M.M.; visualization, M.M., B.K. and D.M.; supervision, G.D.; project administration, G.D. and R.U.; funding acquisition, G.D. and R.U. All authors have read and agreed to the published version of the manuscript.


Introduction
As robotics technology advances, there are more and more efforts to create and control continuum robots that excel in kinematic redundancy. These robots are agile and flexible because of their backbone-less structure, which the biological world has inspired [1,2]. Such robots are called elephant's trunk or snake-arm robots. This study focuses on three inextensible segments of trunk-type robot's position and orientation control. Most trunktype robots are currently used for robotically assisted surgery [3,4], specifically in minimally invasive surgery [5,6]. In the mentioned surgery area, trunk-type robot pose control is crucial, while the continuum robot has to move along a narrow path. This paper presents a control method for trunk-type robots that are more dedicated to tasks requiring robot end-effector spatial control while holding the end-effector position as close as possible to the target. For instance, the authors state that such a concept has the potential to apply continuum robots that are used in urban search and rescue tasks [7], bomb disposal [8], and inspection and repair of aero-engines [9,10]. Continuum robot flexibility and mobility features encourage scientists to look for new kinds of robot-type structure designs and control to adapt robots to other applications.
A continuum robot consists of flexible segments, which, in principle, have a bending motion [1,2]. This bending motion gives two degrees of freedom (DOF) for one section [11]. However, depending on the robot segment type, some segments can also have a linear motion, adding one more DOF. According to this continuum, the robot segment can either have two or three DOF, depending on its type [11,12]. Adding more degrees of of a trunk-type robot with three inextensible segments, with an inverse kinematics solution for the imposed end-effector state. This section also discusses tendon length computation, which is necessary for the inverse dynamics part. In Section 4, the simulation results represent the various bending situations of a three-segment continuum robot with end-effector position and orientation errors from the desired target. Section 5 includes robot segments' motion from the initial position to the target position figures and the electric motor motion profile diagrams.

Related Work
The trunk-type tendon robot's control differs from traditional industrial robots, like Cartesian, cylindrical, polar, and SCARA robots, and from robotic arms or parallel delta robots. The main reason trunk-type tendon control differs from traditional industrial robot control lies in the robot segment construction and motion constraints. While conventional industrial robots utilize rigid links, which move in spatial space with inextensible joints [27], the trunk-type robot does not have rigid links [1,2] but uses a connecting disk ( Figure 1). As mentioned before, a trunk-type tendon robot segment spine is usually flexible, and the connecting disk between adjacent segments neither directly causes bending the section nor is used for segment motion. Joints, in this case, serve as connecting pieces in series to get one continuum robot spine ( Figure 1). The only way to control this trunk-type robot is to bend robot segments. Each of the robot segments' motion is like the arc of a circle, in which the radius is continuously changing while the section is being bent (Figure 2a,b). Knowing this and having a clear view of the robot segment structure sheds light on trunk-type robot design approaches. Some methods implement robot control by using forward kinematics solutions. One example is an operator that directly controls robot actuators via an open-loop user interface [16]. The following example is a robot actuator torque closed-loop control system [28]. Here, the idea is that a robot is used mainly for grasping the object. The robot segment bends around the object and stops when robot actuators produce the same torque as the opposing load torque. The only way to control this trunk-type robot is to bend robot segments. Each of the robot segments' motion is like the arc of a circle, in which the radius is continuously changing while the section is being bent (Figure 2a,b). Knowing this and having a clear view of the robot segment structure sheds light on trunk-type robot design approaches. Some methods implement robot control by using forward kinematics solutions. One example is an operator that directly controls robot actuators via an open-loop user interface [16]. The following example is a robot actuator torque closed-loop control system [28]. Here, the idea is that a robot is used mainly for grasping the object. The robot segment bends around the object and stops when robot actuators produce the same torque as the opposing load torque.
However, there are more trunk-type robot control approaches with inverse kinematics solutions. One proposal is to analyze the trunk-type robot kinematic characteristics with a computer program when the robot's trajectory end-moving platform is assigned [29]. Then, from model parameters using the spatial backbone modal method, the robot's kinematic model is established. Nevertheless, the latter method is only valid for hyperredundant trunk-type robot control. There are also trunk-type robot control techniques where researchers delve into motion planning. One example is the motion planning method for solving the inverse kinematic problems of endoscopic operations of continuum manipulators [30]. The proposed method suggests robot control in a predefined complex environment. Another strategy has the operator control the surgical continuum robot on an arbitrary path in real time and not via a predefined path [31]. The latter paper proposes two-path generation algorithms. One of those algorithms represents an optimization Appl. Sci. 2021, 11, 2946 4 of 28 method with sequential quadratic programming. The other algorithm uses differential kinematics with a PID (Proportional Integral Derivative) control algorithm. The robot inverse kinematic model from the joint space to the wire-length space is the basis for these algorithms' operation. However, there are more trunk-type robot control approaches with inverse kinematics solutions. One proposal is to analyze the trunk-type robot kinematic characteristics with a computer program when the robot's trajectory end-moving platform is assigned [29]. Then, from model parameters using the spatial backbone modal method, the robot's kinematic model is established. Nevertheless, the latter method is only valid for hyperredundant trunk-type robot control. There are also trunk-type robot control techniques where researchers delve into motion planning. One example is the motion planning method for solving the inverse kinematic problems of endoscopic operations of continuum manipulators [30]. The proposed method suggests robot control in a predefined complex environment. Another strategy has the operator control the surgical continuum robot on an arbitrary path in real time and not via a predefined path [31]. The latter paper proposes two-path generation algorithms. One of those algorithms represents an optimization method with sequential quadratic programming. The other algorithm uses differential kinematics with a PID (Proportional Integral Derivative) control algorithm. The robot inverse kinematic model from the joint space to the wire-length space is the basis for these algorithms' operation.
There is a robot inverse dynamics approach with Euler-Lagrange equations, specifically for trunk-type robots with spherical piezoelectric actuators [32]. In [32], the researchers compared the Euler-Lagrange dynamics method with the analytical potential method. The authors conclude that the latter method is more accurate and efficient than the former.
One of the inverse kinematics methods implements robot control using the modified Denavit-Hartenberg procedure (modified D-H table) and Jacobian matrices using this D-H table [33]. This method relies on the idea that the trunk-type robot segment "is bending with constant curvature". Consequently, parameters that define each of the robot segments are an outcome of circle arc parameters (the arc length is ; the curvature is ) and one angle , which defines segment orientation around the z-axis (Figure 3a). To improve the D-H method, researchers propose adaptive fuzzy-based fault-tolerant control, as in [6]. The authors claim that their proposed solution significantly reduces position error and increases the robot control reliability. There is a robot inverse dynamics approach with Euler-Lagrange equations, specifically for trunk-type robots with spherical piezoelectric actuators [32]. In [32], the researchers compared the Euler-Lagrange dynamics method with the analytical potential method. The authors conclude that the latter method is more accurate and efficient than the former.
One of the inverse kinematics methods implements robot control using the modified Denavit-Hartenberg procedure (modified D-H table) and Jacobian matrices using this D-H table [33]. This method relies on the idea that the trunk-type robot segment "is bending with constant curvature". Consequently, parameters that define each of the robot segments are an outcome of circle arc parameters (the arc length is L; the curvature is k) and one angle ϕ, which defines segment orientation around the z-axis (Figure 3a). To improve the D-H method, researchers propose adaptive fuzzy-based fault-tolerant control, as in [6]. The authors claim that their proposed solution significantly reduces position error and increases the robot control reliability.  The circle arc assumption on the robot segment influences most control metho trunk-type robots [34][35][36]. The main distinctions between these methods are the c approach, the number of sections used, and priorities between end-effector positio straint or end-effector orientation constraint. The latter prioritization is crucial b control algorithms might not necessarily enforce both restrictions due to nonholo and dynamics constraints. In this text, the represented method suits a tendon-driven The circle arc assumption on the robot segment influences most control methods for trunk-type robots [34][35][36]. The main distinctions between these methods are the control approach, the number of sections used, and priorities between end-effector position constraint or end-effector orientation constraint. The latter prioritization is crucial because control algorithms might not necessarily enforce both restrictions due to nonholonomic and dynamics constraints. In this text, the represented method suits a tendon-driven robot orientation with three inextensible segments, and is predefined with determined vector and robot-end positioning control, predefined with point coordinates.

Approximate Inverse Kinematics Solution for Imposed End-Effector State
The proposed control method for the position and orientation of a trunk-type robot with three inextensible segments consists of four main steps: 1.
Trunk-type robot bent spine curve approximation in the 3D Cartesian coordinate system.

2.
The bending angles θ n calculation for each robot segment.

3.
Segments' endpoints' x, y, and z coordinates in the Cartesian coordinate system, and sections' endpoint orientation computation.

4.
Robot metal wire (tendon) lengths' calculation according to robot spine curvature and segments' orientations.
Sections 3.1-3.4 explain the above steps in more detail. In the first part, the primary purpose is to get a curve in 3D space such that its length would be equal to the robot spine length, which is predefined and fixed. Here, the soft constraint condition is the goal position, and the mandatory constraint condition is the orientation. The approximate curve endpoint's orientation must be the same as the preset orientation (as a unit vector). However, the curve endpoint must be as close as possible to the predefined robot's endeffector point. Then, further calculations will involve estimation of the approximate target position and mandatory target orientation.
In the second and third steps, the top hard constraint condition is that the robot segment shape has to consist of a regular circle arc form. The approximate robot spine will be the prerequisite for the next curve composed of arcs connected in series in these sections. Figure 3a depicts bending angles θ n , segments endpoints P n , and segments orientation s orienn unit vectors found in this approach.
In the final step, the second and third steps' expressions lead to the estimation of the final segment cables' lengths. Computed robot tendon lengths are necessary for the robot segment position and orientation control (Figure 3b).

Determination of Approximate Robot Spine State Configuration
The solution of the three inextensible segments for the continuum robot starts with a consideration of initial conditions. First, the definition of continuum robot spine full-length l r is necessary. l r is the sum of all three segments' lengths L n (Figure 3a). The three segments' lengths. L 1 , L 2 , and L 3 , explicitly define the full robot spine length. Point G in 3D space represents the target point that the robot end-effector has to reach. Furthermore,ô r . defines the desired robot end-effector orientation at the point G. Figure 4 depicts all the necessary initial conditions and state.
The robot spine's length is a fixed-length estimate and should match l r because the robot segments are inextensible. Knowing this and having information about the preset target point G and orientationô r parametric curve equations are necessary. First, the authors propose an empiric tuning parameter λ a and its expression (1) that will play in final spine state expressions: The solution of the three inextensible segments for the continuum robot starts with a consideration of initial conditions. First, the definition of continuum robot spine fulllength is necessary. is the sum of all three segments' lengths ( Figure 3a). The three segments' lengths. , , and , explicitly define the full robot spine length. Point in 3D space represents the target point that the robot end-effector has to reach. Furthermore, defines the desired robot end-effector orientation at the point . Figure 4 depicts all the necessary initial conditions and state. The robot spine's length is a fixed-length estimate and should match because the robot segments are inextensible. Knowing this and having information about the preset target point and orientation parametric curve equations are necessary. First, the authors propose an empiric tuning parameter and its expression (1) that will play in final spine state expressions: Here, = −31 + 120 + + + (10 + ). The target point is = ( , , ) (in the Cartesian coordinate system). = ( , , ) is the unit vector, which defines the robot's third segment endpoint's direction. is a theoretic speed constant ( = 1). The trunk-type robot spine's full length is , which equals the sum of all three segments' lengths. The closed-form expressions for coordinates of point are contingent upon tuning the coefficient expression in Equation (1). Point expresses the estimated robot spine's endpoint. The latter point is the closest solution to the target point . Point coordinates' equations are necessary for creating estimated robot spine curve parametric equations as functions of time: Here, k r = −31 + 120v + o 2 rx + o 2 ry + o rz (10 + o rz ). The target point is G = x g , y g , z g (in the Cartesian coordinate system).ô r = o rx , o ry , o rz is the unit vector, which defines the robot's third segment endpoint's direction. v is a theoretic speed constant (v = 1). The trunk-type robot spine's full length is l r , which equals the sum of all three segments' lengths. The closed-form expressions for coordinates of point P e3 are contingent upon tuning the coefficient λ a expression in Equation (1). Point P e3 expresses the estimated robot spine's endpoint. The latter point is the closest solution to the target point G. Point P e3 coordinates' equations are necessary for creating estimated robot spine curve parametric equations as functions of time: x e3 = λ a l r o rx + 10 l r x g 12 λ a + 10 l r , y e3 = λ a l r o ry + 10 l r y g 12 λ a + 10 l r , z e3 = l r λ a + λ a o rz + 10 z g 2 (6 λ a + 5 l r ) , Finally, Equations (6)-(8) define robot spine curve parametric equations as a function of time, which will define the spine's curvature: P ec (t) = (x ec (t),ŷ ec (t),ẑ ec (t)), where t represents the time variable, which expresses the spine curve's trajectory in 3D space (in a 3D Cartesian coordinate system), because curve speed is a constant and is preset to v = 1. Equations (10)- (12) represent the robot segment endpoints on a curve as follows: where the n-th segment length is L n (in this case, all robot segments are of the same length). The calculated points from Equations (10)- (12) are shown in Figure 5a-c. PointP ec3 , retrieved using Equation (12), is equal to point P e3 , which originates from Equations (2)-(4). PointP ec3 is the robot curve endpoint near target point G. Figure 5c represents the robot spine curve's solution.
where represents the time variable, which expresses the spine curve's trajectory in 3D space (in a 3D Cartesian coordinate system), because curve speed is a constant and is preset to = 1. Equations (10)- (12) represent the robot segment endpoints on a curve as follows: where the n-th segment length is (in this case, all robot segments are of the same length). The calculated points from Equations (10)- (12) are shown in Figure 5a-c. Point , retrieved using Equation (12), is equal to point , which originates from Equations (2)-(4). Point is the robot curve endpoint near target point . Figure 5c represents the robot spine curve's solution.

Segment's Bending Angles Calculation
Computation of an approximate spine curve in Section 3.1 allows for further estimation of each segment's bending angle. First of all, spine curve parametric equations ( ), ( ), ̂ ( ), as in Equations (6)- (8), have first and second function derivatives. Appendix A gives detailed expressions of the computation formulas for , , , , , and . According to [37], Equations (A1)-(A6) are necessary for the evaluation of arc curvature . Equation (13) provides an expression of arc curvature for a parametrically defined space curve in three dimensions (Cartesian coordinates): Inserting Equation (13) into the expression of arc = returns average bending angles for any segment based on curvature representations:

Segment's Bending Angles θ n Calculation
Computation of an approximate spine curve in Section 3.1 allows for further estimation of each segment's bending angle. First of all, spine curve parametric equationŝ x ec (t),ŷ ec (t),ẑ ec (t), as in Equations (6)  z ec . According to [37], Equations (A1)-(A6) are necessary for the evaluation of arc Appl. Sci. 2021, 11, 2946 8 of 28 curvature k. Equation (13) provides an expression of arc curvature k for a parametrically defined space curve in three dimensions (Cartesian coordinates): Inserting Equation (13) into the expression of arc θ n = kL n returns average bending angles θ n for any segment based on curvature representations: The average bending angle θ n of the n-th segment is one of the main parameters represented in Figure 3a. The following section will present the approach for finding the remaining parameters of the robot segments.

Calculation of Arc Segments' Endpoints and Orientations
All three segments' bending angles θ n (Equation (14)) and endpoints (Equations (10)-(12)) on the curve allow for finding actual segments' endpoint coordinates in the Cartesian system. The further derivation will consider robot segments' shape in an arc form, similar to the ellipse arc seen in the spine curve ( Figure 5c). Arc geometry in the x-z plane 2D workspace ( Figure 6) allows for obtaining the z coordinates of all three robot segments endpoints P n .  As seen in Figure 6, the triangle ∆ will express the coordinate , defining a segment endpoint's z coordinate. Before that, the arc chord length requires determination, for which the formula is Based on circle arc, is equal to Then, the chord length can be expressed as follows: Figure 6. Robot segment in the 2D workspace x-z plane (the segment form is in a regular circle arc shape). d n -chord length (n -segment number); r n -arc radius (n -segment number); L -circle arc length; O -arc start point; P n -arc endpoint (n -segment number); C n -circle center point (n -segment number); A n -middle point of the OP n straight line (n -segment number); α n -angle between p n and c n vectors (n -segment number); θ n -arc angle (n -segment number).
As seen in Figure 6, the triangle ∆OB n P n will express the coordinate z n , defining a segment endpoint's z coordinate. Before that, the arc chord length requires determination, for which the formula is Based on circle arc, r n is equal to Then, the chord length can be expressed as follows: Triangle ∆OP n C n ( Figure 6) serves for determination of the angle α n that lies between vectors p n and c n : From triangle ∆OB n P n , Substituting d n and α n from Equations (17) and (18) into Equation (20), its simplification produces z n Equations (15)-(21) help to estimate coordinates z n of robot 1, 2, 3 segment endpoints P n . On the other hand, point P n coordinates x n , y n can be found similarly to coordinates z n . In this case, all coordinate x n , y n computations are in the first robot segment coordinate system, around point O, (Figure 6). Calculations begin with the estimation of the first segment endpoint's P 1 coordinates, x 1 , y 1 .
First of all, vector p ec1 (pointP ec1 has coordinate estimates based on Equation (10)) has to be normalized so that its magnitude is equal to chord length d n (Equation (17)). According to [38], the general formula for vector normalization and scaling with a given length is where the output vector length is l v , the output vector coordinate is u p (p stands for x, y, or z), the input vector is v, and the input vector coordinate is v p (p stands for x, y, or z). According to the generic normalization and scaling formula (Equation (22)), where the first robot segment chord length is d 1 . Figure 7, based on Equations (23)-(25), represents a new vector, p Nec1 . The magnitude of the calculated new vector p Nec1 is equal to chord length d 1 . However, the direction in the workspace is the same as the vector's p ec1 direction:

‖ ‖
where the first robot segment chord length is . Figure 7, based on Equations (23)-(25), represents a new vector, . The magnitude of the calculated new vector is equal to chord length . However, the direction in the workspace is the same as the vector's direction: Vector p Nec1 gets new values for its coordinates x 1Norm and y 1Norm . Length-invariant direction adaptation of the scaled vector p Nec1 is still necessary. It should rotate in the x-y plane to a new vector p 1 , which defines the final first robot segment endpoint ( Figure 8). To do so, triangle ∆OB 1 P 1 returns the magnitude of b 1 , as in Figure 6: Substituting Equations (17) and (18) into Equation (28) and simplifying ‖ ‖ yields Normalization of and scaling it to the magnitude ‖ ‖ produces new , values: Combining coordinates , , enables the expression of the first segment endpoint . Then, the segment is in a standard arc form, as follows: The first segment endpoint's coordinates have already been obtained. Point and bending angle fully express the robot's first segment. The segment rotation angle around the z-axis in Figure 3a is unnecessary. The approach of finding the parameters of the second and third segments will use the quaternion technique. This method requires knowledge of the orientation of the previous segment. In this case, the retrieval of the second segment's parameters requires knowledge of the first segment's orientation. Initially, for the calculation of the first segment arc center point (Figure 8a), there is a need to define the circle arc radius with the formula Substituting Equations (17) and (18) into Equation (28) and simplifying b 1 yields Normalization of p Nec1 and scaling it to the magnitude b 1 produces new x 1 , y 1 values: Combining coordinates x 1 , y 1 , z 1 enables the expression of the first segment endpoint P 1 . Then, the segment is in a standard arc form, as follows: The first segment endpoint's P 1 coordinates have already been obtained. Point P 1 and bending angle θ 1 fully express the robot's first segment. The segment rotation angle ϕ 1 around the z-axis in Figure 3a is unnecessary. The approach of finding the parameters of the second and third segments will use the quaternion technique. This method requires knowledge of the orientation of the previous segment. In this case, the retrieval of the second segment's parameters requires knowledge of the first segment's orientation. Initially, for the calculation of the first segment arc center point C 1 (Figure 8a), there is a need to define the circle arc radius with the formula Then, point C 1 coordinates estimation yields where z c1 = 0, because the first segment arc's center must lie on the x-y axis plane. In Equations (34)- (36), the coordinates all denote Quaternions will allow for estimation of the first segment's orientation. The calculation starts with finding an axis around which a unit vector defines the first segment orientation rotation. To get this, arc center point C 1 needs to be rotated around the z-axis 90 • counterclockwise ( Figure 9). The definition of the quaternion's axis vector is as follows: By using vector normalization and scaling, as in Equation (22), the vector c r1 magnitude yields a unit vectorĉ r1ĉ r1 = (x urc1 , y urc1 , z urc1 ).
Quaternion q expresses rotation around the unit vector axis → OC um1 through the angle θ 1 , with q = cos θ n 2 + x urcn sin θ n 2 i + y urcn sin θ n 2 j + z urcn sin θ n 2 k = q r + q in i + q jn j + q kn k. (40) Quaternions will allow for estimation of the first segment's orientation. The calculation starts with finding an axis around which a unit vector defines the first segment orientation rotation. To get this, arc center point needs to be rotated around the z-axis 90° counter-clockwise (Figure 9). The definition of the quaternion's axis vector is as follows: By using vector normalization and scaling, as in Equation (22), the vector magnitude yields a unit vector ̂ ̂ = ( , , ).
For any unit vector that will express first segment orientation and rotation, a rotation matrix exists. This rotation matrix rotates any vector by an angle counter-clockwise (right-hand rule) by using For any unit vector that will express first segment orientation and rotation, a rotation matrix R n exists. This rotation matrix rotates any vector by an angle θ n counter-clockwise (right-hand rule) by using 2s q in q jn − q kn q rn 2s q in q kn + q jn q rn 2s q in q jn + q kn q rn 1 − 2s q 2 in + q 2 kn 2s q jn q kn − q in q rn 2s q in q kn − q jn q rn 2s q jn q kn + q in q rn where the vector's scalar/real part is s, that is s = ||q || −2 . The quaternion will only rotate the unit vector; therefore, s = 1. For clockwise rotation around the same axis by an angle θ n , the quaternion rotation matrix is as follows: 2s q in q jn + q kn q rn 2s q in q kn − q jn q rn 2s q in q jn − q kn q rn 1 − 2s q 2 in + q 2 kn 2s q jn q kn + q in q rn 2s q in q kn + q jn q rn 2s q jn q kn − q in q rn The quaternion formula for vector rotation is where the output vector is v r2 and the input vector is v r1 . If v r1 = (0, 0, 1) then the first segment's orientation vector isŝ The sum of vectorŝ orien1 coordinates and point P 1 coordinates is equal to S orien1 . The vector s orienN1 represents the first segment orientation vector ( Figure 10). s orienN1 continues from point P 1 to point S orienN1 , while the vector s orienN1 magnitude ( Figure 10) is scaled up for illustration purposes. segment's orientation vector is = . (44) The sum of vector coordinates and point coordinates is equal to . The vector represents the first segment orientation vector ( Figure 10). continues from point to point , while the vector magnitude ( Figure  10) is scaled up for illustration purposes. To calculate the second segment endpoint's coordinates and orientation vector's coordinates, the first segment calculations from Equation (22) to Equation (44) are still valid. However, the second segment vector (from point to point , Figure 5c), from the second segment coordinate system, needs to be transformed into a first segment coordinate system. Therefore, vector reduces vector to get , which starts from point = (0,0,0): Vector needs to be rotated clockwise with the first segment quaternion rotation matrix in Equation (42) and using the quaternion rotation formula in Equation (43): Vector is a result of the translation to the first segment coordinate system. To calculate the second segment endpoint's P 2 coordinates and orientation vector's s orien2 coordinates, the first segment calculations from Equation (22) to Equation (44) are still valid. However, the second segment vector p ec12 (from pointP ec1 to pointP ec2 , Figure 5c), from the second segment coordinate system, needs to be transformed into a first segment coordinate system. Therefore, vector p ec1 reduces vector p ec2 to get p sec 12 , which starts from point O = (0, 0, 0): Vector p sec 12 needs to be rotated clockwise with the first segment quaternion rotation matrix R alt1 in Equation (42) and using the quaternion rotation formula in Equation (43): Vector p ec2 T is a result of the translation to the first segment coordinate system. Figure 11 depicts the p ec2 T vector. Now the second segment chord length d 2 allows for scaling vector p ec2 T : The estimated vector lies in the first segment coordinate system. Getting the x-y coordinates requires applying Equations (22)-(44) similarly for the second section. The estimated second segment's endpoint and its endpoint's orientation vector exist in the first segment coordinate system (Figure 12). The estimated vector p Nec2 T lies in the first segment coordinate system. Getting the x-y coordinates requires applying Equations (22)-(44) similarly for the second section. The estimated second segment's endpoint and its endpoint's orientation vector exist in the first segment coordinate system ( Figure 12). The estimated vector lies in the first segment coordinate system. Getting the x-y coordinates requires applying Equations (22)-(44) similarly for the second section. The estimated second segment's endpoint and its endpoint's orientation vector exist in the first segment coordinate system ( Figure 12). Translation of the second robot segment endpoint , arc center point , and orientation vector back to the second segment coordination system requires their rotation according to the first segment counter-clockwise quaternion rotation matrix in Equation (41). After that, , and add to the first segment endpoint . Translation of the second robot segment endpoint P 2 T , arc center point C 2 T , and orientation vector s orien2 T back to the second segment coordination system requires their rotation according to the first segment counter-clockwise quaternion rotation matrix in Equation (41). After that, P 2 T , C 2 T and s orienN2 T add to the first segment endpoint P 1 . The second robot arc-like segment's position and orientation are evaluated in the second segment coordinate system ( Figure 13). The second robot arc-like segment's position and orientation are evaluated in the second segment coordinate system ( Figure 13). These second segment calculations apply similarly to the third robot segment. The only difference in this situation is that, in the beginning, vector needs to reduce vector to get ( ), which starts from point = (0,0,0). Then, needs to be rotated clockwise with the first segment quaternion rotation matrix and then Figure 13. Second segment endpoint P 2 and orientation vector s orienN2 in second segment's coordination system.
These second segment calculations apply similarly to the third robot segment. The only difference in this situation is that, in the beginning, vector p ec2 needs to reduce vector p ec3 to get p ec23 (p sec 23 ), which starts from point O = (0, 0, 0). Then, p sec 23 needs to be rotated clockwise with the first segment quaternion rotation matrix R alt1 and then with the second segment quaternion rotation matrix R alt2 . For this purpose, similar calculations using Equations (22)-(44) are necessary.
The transposition of the third robot segment endpoint P 3 T , arc center point C 3 T , and orientation vector s orien3 T (Figure 14) back to the third segment coordination system is then necessary. These points need to be rotated first with the second segment counter-clockwise quaternion rotation matrix and then with the first segment counter-clockwise quaternion rotation matrix. The third segment endpoint P 3 , arc center C 3 point, and orientation vector s orienN3 add up to the second segment endpoint P 2 . The third robot segment or robot end-effector position and orientation endpoint vector lie in the third segment coordinate system, as shown in Figure 15.

Robot Cables' (Tendons') Length Calculation According to the Robot Spine Curvature
The key to controlling the shape of a trunk-type robot is cabling (tendons). One trunktype robot segment can have three or four (or more) tendons dedicated to control. This section will discuss the robot segments' options with only three cables for the control of each segment. Each segment will have equations for three robot cable lengths (nsegment number, i-cable number).

Robot Cables' (Tendons') Length Calculation According to the Robot Spine Curvature
The key to controlling the shape of a trunk-type robot is cabling (tendons). One trunktype robot segment can have three or four (or more) tendons dedicated to control. This section will discuss the robot segments' options with only three cables for the control of

Robot Cables' (Tendons') Length Calculation According to the Robot Spine Curvature
The key to controlling the shape of a trunk-type robot is cabling (tendons). One trunktype robot segment can have three or four (or more) tendons dedicated to control. This section will discuss the robot segments' options with only three cables for the control of each segment. Each segment will have equations for three robot cable lengths l ni (n-segment number, i-cable number).
Some of the cable length formulas originate from [33]. The authors adjusted these equations and used them in the robot control approach. The main reason for choosing these equations was that [33] used robot segment spine curvature and orientation information about robot segment construction. In [33], the equations were based on assumptions that their cables are in straight lines between the section connecting disks and not in the arc-like segment spine. However, the authors of [33] assessed the number of cable guide connecting disks in the section (Figure 16), which is essential for accurate robot cable length calculation.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 17 of 29 guide connecting disks in the section (Figure 16), which is essential for accurate robot cable length calculation. Figure 16. A trunk-type robot with three segments, each having three tendons.
In [33], the equations for the first segment cable length's calculation are − cos where the robot first segment cable lengths are , , ; the first segment spine curvature is ; the first segment circle arc radius is ; the distance from the robot spine center point to the robot cable center point on the robot connecting disk plane of the first segment is ; is the angle between a vector and the first segment coordinate system x-axis on the x-y plane (the angle around the z-axis). The angle range is 0,2π from the x-axis to vector on the x-y plane in a counter-clockwise direction. comes from this equation: where the n-th segment length is and the number of spaces between the connecting disks of the n-th segment is . The first robot cable holes' positions on the connecting disks are necessary to determine the cable lengths for the second and third segments of the robot. Cables conclude three groups: for the first robot segment control , , ; for the second robot segment , , , and for the third robot segment , , . According to this, nine cable holes have to be on the robot's first segment connecting disks, six cable holes have to be on the robot second segment connecting disks, and three cable In [33], the equations for the first segment cable length's calculation are where the robot first segment cable lengths are l 11 , l 12 , l 13 ; the first segment spine curvature is k 1 ; the first segment circle arc radius is r 1 ; the distance from the robot spine center point to the robot cable center point on the robot connecting disk plane of the first segment is D 1 ; ϕ 1 is the angle between a vectorŝ orien1 and the first segment coordinate system x-axis on the x-y plane (the angle around the z-axis). The angle ϕ 1 range is [0, 2π] from the x-axis to vectorŝ orien1 on the x-y plane in a counter-clockwise direction. q 1 comes from this equation: q n = 2m n sin k n L n 2m n 1 kn =r n → 2m n sin L n 2m n r n , where the n-th segment length is L n and the number of spaces between the connecting disks of the n-th segment is m n . The first robot cable holes' positions on the connecting disks are necessary to determine the cable lengths for the second and third segments of the robot. Cables conclude three groups: for the first robot segment control l 11 , l 12 , l 13 ; for the second robot segment l 21 , l 22 , l 23 , and for the third robot segment l 31 , l 32 , l 33 . According to this, nine cable holes have to be on the robot's first segment connecting disks, six cable holes have to be on the robot second segment connecting disks, and three cable holes have to be on the robot third segment disks. Punches used for robot cables from the same cable group have to be spaced 120 • along all connecting disks. Cable hole positions were first defined in the first segment connecting disks and later on the second and third segment cable guide disks. The starting robot segment cable hole on the first segment guide disk is on the y-axis for first segment control cable l 11 (Figure 17a). The holes for cables l 12 , l 13 are placed correspondingly. Disk holes for cables l 21 , l 22 , l 23 on the first segment are simply the cable l 11 , l 12 , l 13 holes rotated 40 • clockwise around point O. Furthermore, the cable l 31 , l 32 , l 33 holes were rotated the same way as l 21 , l 22 , l 23 but by 80 • . Figure 17a- Third, position error results from multiple connected arcs that approximately replace the estimated bent robot spine curve. In most simulations, the spine curve endpoint is closer to the desired goal point, but robot configuration requires arcs as section shapes. Therefore, a minor drift of the spine endpoint is acceptable to impose the segment shape nonholonomic constraint. During this recalculation process, the first task is to get average bending angles for each arc segment. We assume that the angles resulting from curvature integration (Equation (14)) must match similar arc sections' curvature.
Additionally, the lengths of both arc and elliptical curves are the same. However, the chord lengths of arc and elliptical angles are different. The latter explains why position errors between estimated robot spine curves and calculated robot segments arcs endpoints occur in all cases. Only when → 0 does this position error between the arc and elliptical curve vanish to a minimum.     Segment length L n (cm) 40 40 40 Distance from robot spine to the tendon on connecting disk D n (cm) 5 5 5 In this paper, the simulations' goal is to evaluate the developed method's efficiency and accuracy. In this case, robot position and orientation errors will express the accuracy of the method. Position error is essentially the scalar Euclidean distance for the robot pose [39]: The authors verified the study routine's scalability when inverse dynamics problems affect the planning of all robot actuators' synchronous motion. For this purpose, the authors transformed mathematical scripts into a C program and visually inspected the potential benefits of applying the routines for solving inverse dynamics tasks. A single state (a single target endpoint) spine estimation, (Figures 18-21) lasted 200-300 ms on a single core, where the processor had 1.80 GHz, and the RAM was 4 GB. Therefore, the robot control method presented herein is acceptable for both online and offline scenarios. Currently, experimental-practical validation of the three-segment trunk-type tendon robot is ongoing. The authors will reveal research results in the future.

Electric Motor Speed Control Profiles
Frequently, an electric motor controls tendon lengths [1,2,5], and such a drive represents a typical actuator. More specifically, electric motors can control trunk-type robot metal wires (tendons) by using Equations (48)-(57) for all three robot segments. Such open-loop (prediction/planning) motion profiles will automatically impose final rigid body motion that practically solves the inverse dynamic problem. That is why this section further elaborates on simulation results dedicated to motor speed control profiles.
First, this study presents angular speed profiles for any generic electric motor instal- Here, the robot third segment endpoint found with the proposed method is P 3 . The target point that the robot end-effector has to reach is G. According to the authors' decision, the dot product of the third segment orientationŝ orienN3 and predefined orientationô r , at the target point G, will determine the orientation error. Lipschutz et al. [40] give the mentioned dot product between the two vectors: where the first vector is a and the second vector b is the angle θ between them. Redesigning the dot product via Equation (59) produces the orientation error expression: Table 2 exposes orientation errors between vectorsô r and s orienN3 and position errors between goal point G and bent robot third segment endpoint P 3 . From Table 2, it is clear that the algorithm imposes both constraint conditions: soft for segment goal point and hard for robot end orientation. However, substantial position errors, defined as the distance between desired and achieved robot endpoints, can be seen in some situations, Table 2. For example, the fifth simulation results contain a 27.0281 cm position error. Three main concerns affect the bent robot position accuracy. First, some goal points that have to reach the robot end effector's target position exceeded the robot's configuration space. Second, the position error is noticeable because the robot position is soft-constrained, while the primary condition is orientation. Figures 18-21 show that reaching the desired locations is possible in most simulations but not with the desired robot end orientations. This is the main reason for the substantial position errors in some cases. Third, position error results from multiple connected arcs that approximately replace the estimated bent robot spine curve. In most simulations, the spine curve endpoint is closer to the desired goal point, but robot configuration requires arcs as section shapes. Therefore, a minor drift of the spine endpoint is acceptable to impose the segment shape nonholonomic constraint. During this recalculation process, the first task is to get average bending angles θ n for each arc segment. We assume that the angles θ n resulting from curvature integration (Equation (14)) must match similar arc sections' curvature.
Additionally, the lengths of both arc and elliptical curves are the same. However, the chord lengths of arc and elliptical angles are different. The latter explains why position errors between estimated robot spine curves and calculated robot segments arcs endpoints occur in all cases. Only when θ n → 0 does this position error between the arc and elliptical curve vanish to a minimum.
The authors verified the study routine's scalability when inverse dynamics problems affect the planning of all robot actuators' synchronous motion. For this purpose, the authors transformed mathematical scripts into a C program and visually inspected the potential benefits of applying the routines for solving inverse dynamics tasks. A single state (a single target endpoint) spine estimation, (Figures 18-21) lasted 200-300 ms on a single core, where the processor had 1.80 GHz, and the RAM was 4 GB. Therefore, the robot control method presented herein is acceptable for both online and offline scenarios. Currently, experimental-practical validation of the three-segment trunk-type tendon robot is ongoing. The authors will reveal research results in the future.

Electric Motor Speed Control Profiles
Frequently, an electric motor controls tendon lengths [1,2,5], and such a drive represents a typical actuator. More specifically, electric motors can control trunk-type robot metal wires (tendons) by using Equations (48)-(57) for all three robot segments. Such open-loop (prediction/planning) motion profiles will automatically impose final rigid body motion that practically solves the inverse dynamic problem. That is why this section further elaborates on simulation results dedicated to motor speed control profiles.
First, this study presents angular speed profiles for any generic electric motor installation with or without a gearbox. The only customization parameter value is the radius of a motor rotor or a gear. The parameter helps translate rotational mechanical motion to longitudinal motion, which directly controls a tendon. Accordingly, the speed motion profiles and the kit radii provide an opportunity to plan torque profiles for electric motors and their gearboxes, if any, in the future.
Next, an expression for tendon lengths' rate of change is necessary to arrive at the angular speed contours. Based on Equations (48)-(57), the approximation of the first-order derivative of cable lengths, as a function of time (Figure 22), is as follows: where a cable length l nm belongs to the robot's n-th segment and its m-th cable, the observation's discrete-time index is i, and ∆t is the sampling interval in seconds. Calculating the robot cable length longitudinal speed motor angular speed is then done as follows: where the motor rotor or gear radius is r ms . Motor shaft angular speed also involves the representation of rotational speed: where motor or gear rotational speed is n nm , which belongs to the n-th segment and its m-th cable. Combining and simplifying Equations (62) and (63) produces the robot motor shaft rotational speed as follows: where a cable length belongs to the robot's n-th segment and its m-th cable, the observation's discrete-time index is i, and ∆ is the sampling interval in seconds. Calculating the robot cable length longitudinal speed motor angular speed is then done as follows:   The solution of inverse dynamics requires the robot's end-effector initial state and target state set. The fifth robot simulation (Figure 20a) presents the target G and robot end-effector orientation unit vectorô r . There are two scenarios discussed that differ in their initial state. The first scenario corresponds to a vertically straight form of the manipulator. The second case will circumscribe the robot as initially bent. Figures 23 and 24 show both situations and their corresponding discrete surfaces that practically represent the inverse dynamics' numeric solution. Equations (61)-(64) allow for estimations of speed control profiles in both scenarios (Figures 23 and 24). Simulations involved a radius value of 1 cm for all actuators. The total modeling time for the robot to reach the target position from its initial state was 10 s.  As in Figure 25, the profile information brings future opportunities for online and offline robot motion planning and prediction when obstacles are present and different trajectory surface configurations come into effect. Simultaneously, these profiles are synchronous, i.e., their open-loop execution will result in final experimental trajectories close to the theoretical solutions proposed by this study.

Conclusions
This study proposes an approach for finding synchronized actuator profiles that (by design) solve inverse dynamics problems for a trunk-type robot manipulator consisting of three segments. The intention of the method is the spatial control of the three inextensible segments' positions and orientations. The technique imposed the soft constraint assumption for the segment's goal point and the mandatory constraint condition for the end-effector's alignment. The method's principle is to obtain an approximate curve of the spatial spine that meets the listed position and orientation conditions. Then, a robot is "built" around the estimated backbone by using circular arcs. The calculated final robot configuration determines all segments' spatial positions and orientations. The obtained robot configuration also presents all lengths of the corresponding tendons so that their   As in Figure 25, the profile information brings future opportunities for online and offline robot motion planning and prediction when obstacles are present and different trajectory surface configurations come into effect. Simultaneously, these profiles are synchronous, i.e., their open-loop execution will result in final experimental trajectories close to the theoretical solutions proposed by this study.

Conclusions
This study proposes an approach for finding synchronized actuator profiles that (by design) solve inverse dynamics problems for a trunk-type robot manipulator consisting of three segments. The intention of the method is the spatial control of the three inextensible segments' positions and orientations. The technique imposed the soft constraint assumption for the segment's goal point and the mandatory constraint condition for the end-effector's alignment. The method's principle is to obtain an approximate curve of the spatial spine that meets the listed position and orientation conditions. Then, a robot is "built" around the estimated backbone by using circular arcs. The calculated final robot configuration determines all segments' spatial positions and orientations. The obtained robot configuration also presents all lengths of the corresponding tendons so that their motion implicitly is synchronous. The average time of the code execution varied between 200 and 300 ms.
This study shows that, depending on the specified target point and orientation, the end-effector orientation simulation error ranged from 0 • to 10 • . Since the robot's position was not the highest priority in this study, the maximum recorded position error was 27.028 cm and the lowest was 2.242 cm. The total length of the robot's spine during the simulations was 120 cm. When estimating the error, it is necessary to consider that the authors provide the error statistics regardless of the robot workspace. However, the described method suggested finding the best robot configuration under defined conditions with closed-form expressions, avoiding neural networks [41] or robot motion learning [42] techniques. However, as this work clarifies, the accuracy of the end-effector's state is the main trade-off.
The approach also presents an approximate discrete inverse dynamics solution relevant for online and offline trunk-type robot motion planning and navigation tasks. This text discusses the rotational profile diagrams for two scenarios. Since the provided motor speed profiles are synchronized, the manipulator's motion would result in trajectories close to those shown in Figures 23 and 24. The experimental research on a three-segment tendon robot with similar technical specifications to the robot in the simulations is ongoing. The results of this research will appear in the future.