A control algorithm for a centrifuge motion simulator

Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

Contents lists available at ScienceDirect

Robotics and Computer-Integrated Manufacturing journal homepage: www.elsevier.com/locate/rcim

A control algorithm for a centrifuge motion simulator Vladimir M. Kvrgic n, Jelena Z. Vidakovic, Maja M. Lutovac, Goran Z. Ferenc, Vojkan B. Cvijanovic Lola Institute, Kneza Viseslava 70a, 11030 Belgrade, Serbia

art ic l e i nf o

a b s t r a c t

Article history: Received 13 September 2013 Received in revised form 17 January 2014 Accepted 29 January 2014 Available online 20 February 2014

Pilots of modern combat aircraft are exposed to the devastating effects of high acceleration forces. The pilots' ability to perform tasks under these extreme flight conditions must be examined. A centrifuge motion simulator for pilot training is designed as a 3DoF manipulator with rotational axes. Through rotations about these axes, acceleration forces that act on aircraft pilots are simulated. Because of the possibilities of present actuators, it is notably difficult to produce a centrifuge that can realise all of the desired changes of the acceleration forces completely accurately. For this reason, it is necessary to make a compromise in the centrifuge's design with regard to the motor choices and link designs. A new control algorithm that contains a new algorithm for the inverse dynamics of the robots (based on the recursive Newton–Euler algorithm) and that accounts for the possible motor actions has been developed in this study. This algorithm first calculates the successive actuator torques of the links, which are required for the given motion during each interpolation period. Next, the algorithm checks whether actuators can achieve these torques, and if they cannot, it calculates the maximal successive link angular accelerations that motors can achieve. Based on this, control unit sends appropriate control inputs. As a result, the quality of the motion control is improved, and a precise calculation of the forces and the moments that act on the centrifuge links (which is necessary to calculate the link strengths) is performed. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Centrifuge flight simulator Kinematics Dynamics Control algorithm Robotics

1. Introduction Modern thrust-vectored jet aircraft have the capability of developing multi-axis accelerations, especially during the performance of “super- manoeuvres” [1–3]. These “agile” aircraft are capable of unconventional flight with high angles of attack, high agile motions with thrust-vectored propulsion in all 3 aircraft axes, rotations around those axes and accelerations of up to 9g (g is Earth's acceleration), with acceleration rates (jerk) of up to 9g/s [4,5]. The human consequences of this agile flight environment are still not completely known [3]. The connexion between the jerk and the movements of the human body is shown in [6–10]. Hence, the destructive effects of the high acceleration forces and the rapid changes of these forces on the pilot's physiology and the ability to perform tasks under these flight conditions must be tested. A human centrifuge is used for the reliable generation of high G onset rates and high levels of sustained G, to test the reactions and the tolerances of the pilots. Here, acceleration force G ¼a/g,

n Corresponding author. Tel.: þ 381 112 541 303, mobile: þ 381 63 288 349; fax: þ381 112 544 096. E-mail address: [email protected] (V.M. Kvrgic).

http://dx.doi.org/10.1016/j.rcim.2014.01.002 0736-5845 & 2014 Elsevier Ltd. All rights reserved.

a ¼ ða2n þa2t þ g 2 Þ1=2 is the magnitude of acceleration acting on the pilot, an is normal, and at is the tangential acceleration. The centrifuge (Fig. 1) has the form of a three degree-of-freedom (3DoF) manipulator with rotational axes, where the pilot's head (or chest for some of the training) is considered to be the end-effector [11,12]. The arm rotation around the vertical (planetary) axis is the main motion that achieves the desired acceleration force. Centrifuge flight simulators must achieve velocity, acceleration and jerk of the pilot through suitable rotations of the centrifuge arm about this axis. The arm carries a gimballed gondola system, with two rotational axes providing pitch and roll capabilities. The roll axis lies in the plane of the arm rotation, perpendicular to the main rotational axis, i.e., in the x-axis direction. The pitch (y) axis is perpendicular to the roll axis. A similar centrifuge, driven by three hydraulic actuators, is described in [13]. In [14,15], another realisation of the centrifuge is described. Its second axis is perpendicular to the first axis and is along the horizontal line when the centrifuge is in a neutral position. The task of the roll and pitch axes is to direct the acceleration force into the desired direction. It is considered that the pilot's head (chest) is placed in the intersection of the gondola's roll and pitch axes. In this way, the centrifuge produces the transverse Gx, lateral Gy and longitudinal Gz acceleration forces ^ x , pitch ω ^ y and yaw ω ^ z angular velocities to simulate and the roll ω the aircraft's acceleration forces and angular velocities.

400

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

Nomenclature a

magnitude of acceleration acting on the simulator pilot an, at, g normal (radial), tangential and Earth's acceleration G simulator pilot's acceleration force; G0 is this force for q1 ¼0 Gx, Gy, Gz transverse, lateral and longitudinal acceleration forces Gn, Gt, Gv normal (radial), tangential and vertical acceleration forces ω^ x , ω^ y , ω^ z simulator pilot's roll, pitch and yaw angular velocities q1, q2, q3 arm (ψ), roll ring (ϕ) and gondola rotation angle (θpitch) q_ i , q€ i angular velocities and angular accelerations of the link i a1 arm length n Tm , n Dm , n pm homogenous transformation matrix, orientation matrix and position vector Rotðx; qÞ rotation transformation matrix of rotating q about x Transðx0 ; a1 Þ translation transformation matrix of translating a1 along x0 Δt interpolation period Δx change of x in one interpolation period n increase in the acceleration force G nd desired acceleration rate of change Δn increase/decrease of n in one interpolation period at the starting or ending stage of the acceleration change xprev value of x in the previous interpolation period as, ae initial and desired acceleration of the simulator pilot cac linear acceleration change part of the total acceleration change abs(x), sign(x) absolute value and sign of x Na number of interpolation periods in the starting or ending stage of the G load change Nc number of interpolation periods for n¼ const q_ 1r , q_ 1 max rated and maximal angular speed of link 1 P1r, M1r, n1r rated power, rated torque and rated number of revolutions of the axis 1 motor n1m, n1max number of revolutions and maximal number of revolutions through field weakening of the axis 1 motor

The presented dynamic environment simulator gimballed centrifuge is aimed not only at improving þGz tolerance but also at the combined Gy/Gz and Gx/Gz exposure. Multi-axis sustained accelerations can either enhance or reduce the þ Gz tolerance of the pilot, depending on the direction of the net gravitoinertial force. Gy acceleration in conjunction with Gz acceleration can enhance G tolerance. Gx acceleration in addition to Gz acceleration can reduce the G tolerance [3]. Although the centrifuge is capable of generating acceleration forces of up to 15g for materials testing purposes, forces that are less than or equal to 9g are used for pilot training. Because of the possibilities of present actuators, it is very difficult to produce a centrifuge that can realise all of the given changes of the acceleration forces completely accurately. Actuators that have the desired power have very large weights. For this reason, it is necessary to make compromises in the centrifuge design regarding the powers and weights of the chosen motors and the strengths, masses, and mass distribution (mass centre positions and moments of inertia) of the links. The gondola must be large enough to carry the pilot, his seat and the equipment.

f1p k1, I1

overload capability of the axis 1 motor gear ratio and efficiency of the gearbox of the axis 1 moment of inertia of the rotor and the gear box elements of the axis 1 brought down on that rotor M i max , q_ i max maximal torque and maximal speed for that torque, i¼ 2, 3 Mia maximal motor torque related to the motor speed, i¼ 2, 3 _i ωi , ω angular velocity and angular acceleration of the link i (i¼1,2,3) vi , v_ i linear velocities and linear accelerations of the link i (i¼1,2,3) v_ cm linear accelerations of the link i (i¼ 1,2,3) center of i mass and the external load (i¼4) center of mass rcm position of the link i (i¼1,2,3) center of mass and the i external load (i ¼4) center of mass with respect to the link i coordinates expressed in the base coordinates cm r^ i position of the link i (i¼ 1,2,3) and the external load (i¼4) center of mass with respect to the link i coordinates expressed in the link i coordinates mi mass of the link i (i¼1,2,3) and the external load (i¼4) Icm moment of inertia matrix of link i (i¼1,2,3) and the i external load (i¼ 4) about the centre of mass of link i expressed in the base link coordinates ^Icm moment of inertia matrix of link i (i¼1,2,3) and the i external load (i¼ 4) about the centre of mass of link i expressed in link i coordinates Fi , Mi total force and the total moment exerted on link i (i¼1,2,3) f i , mi force vector and the moment vector exerted on link i by link i  1 with respect to the base coordinate frame (i¼1,2,3) ^i f^ i , m force and moment exerted on link i by link i  1 in link i  1 coordinates (i ¼1,2,3) ^ zi , m ^ xyi torque of the joint i actuator and the moment acting m on the bearing i (i¼1,2,3) f^ ai , f^ ri axial and radial force of the bearing i (i ¼1,2,3) I ti total moment of inertia of the link i and the external load, reduced to the axis zi-1 (i¼1,2,3) j r xi , j r yi , j r zi x, y and z coordinates of the link i (i¼ 1,2,3) and the external load (i ¼4) center of mass with respect to the link j coordinates expressed in the base coordinates

η1

Easy access for the medical staff in the case of undesired health problems caused by the inertial loads must be provided as well. On the other hand, it has to be as light as possible. The same requirements apply for the roll ring and the centrifuge arm. Within the design of high-performance machines that imply rapid movements and high values of velocities, such as the human centrifuge, the dynamics of the manipulator plays an important role in achieving such high-speed performance [16,17]. A dynamic model can be used for a computer simulation of a robotic system. Through an overview of the behaviour of the model under various operating conditions, it is possible to predict how a real system will behave. The control algorithm for a centrifuge motion simulator that is presented in this paper calculates all of the centrifuge kinematic and dynamic parameters in each interpolation period, to predict its dynamic behaviour. This method includes a new algorithm for inverse dynamics of robots that calculates first the successive actuator torques and the angular accelerations of the links that are needed for the given motion. Then, it checks whether the actuators can achieve these torques and accelerations, and if they

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

401

of q€ 1 , q€ 2 and q€ 3 , and the control algorithm of the centrifuge movement are given in Section 4. Verification of the results from using the proposed control algorithm is presented in Section 5.

2. Kinematics and dynamics of the centrifuge

Fig. 1. Centrifuge with 3 degrees of freedom.

cannot, it calculates the maximal successive link angular accelerations that the motors can achieve. Instead of sending unachievable commands to actuators, the control unit sends commands that give the maximum possible values for the angular accelerations and angular speeds. This strategy improves the quality of the motion control and enables a more precise calculation of the forces and moments that act on the centrifuge links, which is necessary for the link strength calculations that are performed during the centrifuge design. The robot dynamic model can be derived by different methods, such as the Newton–Euler method, the Lagrange equation and the virtual work principle [18–20]. The method based on the Lagrange formulation is conceptually simple and systematic [18,21]. The method based on the Newton–Euler formulation yields a model in a recursive form. It is composed of a forward computation of the velocities and accelerations of each link, followed by a backward computation of the forces and moments in each joint [16–18,21]. Algorithm is computationally more efficient because it exploits the typically open structure of the manipulator kinematic chain [21,22]. On the other hand, the Newton–Euler procedure is difficult to use in a control application of the closed structure because the expense of the calculation is very high [23–25]. The new algorithm for the inverse dynamics of robots is based on the Newton–Euler equations of motion because they incorporate all of the forces that act on the individual links of a robot. This is essentially for sizing the links and bearings during the design stage. The developed control algorithm includes a simulation of the inertial loads that were exerted on the pilot during flight. This approach allows for virtual centrifuge prototyping and analysis of the interaction between the engine control, the simulation of inertial forces acting on the pilot and the structure of the centrifuge. 3D simulation for human centrifuge, used for testing and verification of presented control algorithms, is described in [26]. Although the centrifuge links are flexible to some degree, the modelling and dynamic analysis of robots with flexible links [27] has not been used. During kinematic and dynamic modelling of the centrifuge, the small elastic deformation of the centrifuge links is neglected. As a result, the obtained kinematic and dynamic model is not completely accurate, but the calculated forces and moments do not differ much from their actual values. A kinematics and dynamics analysis of the centrifuge is given in Section 2. The calculation of the acceleration forces that act on the pilot in the gondola and required link angles for second and third axis of the centrifuge is given in Section 3. Calculation of the angular acceleration of the centrifuge arm q€ 1 , the smoothing of its profile, the calculation of the desired and maximal possible values

This section defines the coordinate frames for the centrifuge links (Fig. 2) and the matrices that determine their relations. The centrifuge links and their coordinate frames are denoted by using the Denavit–Hartenberg convention (D–H). The base is denoted by 0, the arm by 1, the roll ring by 2 and the gondola by 3. The arm rotation angle is denoted by q1 ¼ ψ, the roll ring rotation angle by q2 ¼ ϕ and the gondola rotation angle (pitch) by q3 ¼ θ. The centrifuge that was developed as a research result presented in this paper has the following features: arm length a1 ¼ 8 m, roll axis rotation range of 7 1801 and pitch axis rotation range 7 3601. The allowed gondola payload is 250 kg. The centrifuge base coordinates are denoted by x0y0z0, the arm coordinates by x1y1z1 (link 1), the roll ring coordinates by x2y2z2 (link 2), the gondola coordinates by x3y3z3 (link 3) and the pilot coordinates by xyz. Here, x3 ¼x, y3 ¼y and z3 ¼z. The D–H parameters for the 3-axis centrifuge components are given in Table 1. The 4  4 homogenous matrix that transforms the coordinates of a point from frame xnynzn to frame xmymzm is denoted by " # n n Dm pm n Tm ¼ 0 0 0 1 and from x0y0z0 to xmymzm by Tm. n Dm is a 3  3 orientation matrix, and n pm is a 3  1 position vector. This matrix, which describes the relation between one link and the next, is called i  1Ai ¼ A(i  1,i). By using the convenient shorthand notation, sin ðq1 Þ ¼ s1 , cos ðq1 Þ ¼ c1 , sin ðq2 Þ ¼ s2 , cos ðq2 Þ ¼ c2 , sin ðq3 Þ ¼ s3 and cos ðq3 Þ ¼ c3 , the following homogenous matrices for the relation between the centrifuge link coordinate frames are defined to derive the kinematic equations for the machine, as follows: T1 ¼ 0 Α1 ¼ Rotðz0 ; q1 ÞTransðx0 ; a1 ÞRotðx00 ; 90 3 Þ 2 3 s1 a1 c 1 c1 0 6s 7 6 1 0  c 1 a1 s 1 7 ¼6 7 40 1 0 0 5 0

0

0

ð1Þ

1

Fig. 2. Coordinate frames of the 3-axis centrifuge links.

Table 1 D–H parameters for the 3-axis centrifuge links. Link

Variable [deg]

a [mm]

d [mm]

α [deg]

1 2 3

q1 q2 þ 90 q3 þ 90

a1 0 0

0 0 0

90 90  90

402

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

2 1

 s2

6 c 2 4 0

Α2 ¼ Rotðz1 ; q2 ÞRotðz1 ; 90 3 ÞRotðx01 ; 90 3 Þ ¼ 6 6

0 2 2

c2

0 1

s2 0

07 7 7 05

0

0

1

0

 c3

0

 s3

1

0

0

0

 s3

6 c 3 4 0

Α3 ¼ Rotðz2 ; q3 ÞRotðz2 ; 90 3 ÞRotðx02 ;  90 3 Þ ¼ 6 6

0

0

3

0

ð2Þ

0

3

07 7 7 05 1 ð3Þ

The forward kinematics that is related to the robot geometry is used to calculate the position and orientation of the links and endeffector (in this case, the pilot's head/chest) with respect to the centrifuge variables q1, q2 and q3. It is determined from the following matrix: 2 3 c1 s2 s1 c 1 c 2 a1 c 1 6 s s c1 s1 c2 a1 s1 7 1 2 6 7 T 2 ¼ 0 Α1 1 Α2 ¼ 6 ð4Þ 7 4 c2 0 s2 0 5 0 2

0 c1 s2 s3 þ s1 c3

6 s s s c c 1 3 6 1 2 3 T 3 ¼ 0 Α1 1 Α2 2 Α3 ¼ 6 4  c2 s3 0

0

1  c1 c2

c1 s2 c3  s1 s3

 s1 c2

s1 s2 c3 þ c1 s3

s2

 c2 c3

0

0

a1 c1

3

a1 s1 7 7 7 0 5 1

Fig. 3. The transverse, lateral and longitudinal acceleration force components Gx, Gy and Gz, which act on the pilot in the simulator.

ð5Þ On the basis of these transformational matrices, the equations of forward kinematics that relate to the velocities and accelerations of the links and end-effector-pilot's head/chest are developed. From those equations, inverse dynamics algorithms, which involve the determination of actuator torques that are required to generate a desired trajectory of the robot-centrifuge and other torques and forces in the centrifuge axes, are then developed. These algorithms are given in Appendix A.

3. Acceleration force components and roll and pitch link angles 3.1. Calculation of the simulator pilot acceleration force components Fig. 3 shows the transverse Gx, lateral Gy and longitudinal Gz acceleration force G components that act on the pilot's head (chest) in the simulator. The three main axes of the coordinate frame attached to the human body are: the x axis, which extends from the face to the back, the y axis, which extends from the left to the right side, and the z axis, which extends from the head to the pelvis. Fig. 4 shows coordinate frames, angles, angular velocities and acceleration forces of the centrifuge. The linear acceleration experienced by the simulator pilot at the intersection point of the roll and pitch axes is h iT h iT v_ 1 ¼ v_ 2 ¼ v_ 3 ¼ v_ 1x v_ 1y v_ 1z ¼ a1 s1 q_ 1 þ c1 q21 c1 q_ 1  s1 q21 0 ð6Þ Based on Eq. (6) for q1 ¼0 and adding the gravitational acceleration g, the orthogonal components Gn, Gt and Gv for the normal (radial), tangential and vertical acceleration force G components, respectively, that are experienced by the simulator pilot are the following: 3 2 2 3 3 2 3 2 an Gx0 Gn a1 q21 =g 7 6 6 G 7 16 a 7 6 7 ð7Þ 4 y0 5 ¼ 4 t 5 ¼ 4  a1 q_ 1 =g 5 ¼ 4  Gt 5 g g Gz0  Gv 1

Fig. 4. Coordinate frames, angles, angular velocities and acceleration forces of the centrifuge.

The link angles q2 ¼ ϕ and q3 ¼ θ and the angular velocity q_ 1 of the arm define the orthogonal components Gx, Gy and Gz of the resultant vector G that are experienced by the simulator pilot. Based on Eqs. (5) and (7), the resultant vector G is h iT h iT G ¼ Gx Gy Gz ¼ D3 1 Gx0 Gy0 Gz0 ð8Þ Gx ¼ s3 ðGx0 s2 þc2 Þ  Gy0 c3

ð9Þ

Gy ¼  Gx0 c2 þ s2

ð10Þ

Gz ¼ c3 ðGx0 s2 þ c2 Þ þ Gy0 s3

ð11Þ

Angles q1 ¼ ψ, q2 ¼ ϕ and q3 ¼ θ and their derivatives define the roll, ^ x, ω ^ y and ω ^ z , which are pitch and yaw angular velocities of ω experienced by the simulator pilot; they are given in the following equation (for q1 ¼0): ^x ω ^y ω ^ z  T ¼ D  1 ω3 ^ ¼½ω ω 3

ð12Þ

ω^ x ¼ q_ 2 c3  q_ 1 c2 s3

ð13Þ

ω^ y ¼  q_ 3  q_ 1 s2

ð14Þ

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

ω^ z ¼  q_ 2 s3  q_ 1 c2 c3

ð15Þ

3.2. Calculation of the centrifuge roll and pitch angles The centrifuge roll angle is calculated by Eq. (10), which uses the given lateral force Gy, in the following way: q2 ¼ ϕ ¼ atan2ðGx0 þGy ð1 G2y þ G2x0 Þ1=2 ; 1  G2y Þ

ð16Þ

If Gy o 0 and G2y 4 1, then q2 ¼ q2 þ π . The roll angle can be calculated only if G2x0 þ 1 Z G2y . Otherwise, it is not possible to achieve the given lateral force Gy. For Gy ¼0, Eq. (16) yields: q2 ¼ ϕ ¼ atan2ðGx0 ; 1Þ

ð17Þ

Eqs. (9) and (11) show that it is not possible to achieve both of the given Gx and Gz forces, even when they are in the allowed ranges. As a result, the centrifuge pitch angle is calculated by Eq. (9), using the given transverse force Gx, or by Eq. (11), using the given longitudinal Gz force. Eq. (9) yields: q3 ¼ θ ¼ atan2ðGy0 b þ Gx ðb þ G2y0  G2x Þ1=2 ; b  G2x Þ 2

2

ð18Þ

2

where b ¼ Gx0 s2 þ c2 . If b þ G2y0 o G2x , then it is not possible to achieve the given transverse force Gx. For Gx ¼ 0, Eq. (17) yields: q3 ¼ θ ¼ atan2ðGy0 ; bÞ

ð19Þ

Eq. (11) yields the following: q3 ¼ θ ¼ atan2ðGy0 d  Gz ðb þ G2y0  G2z Þ1=2 ; G2z  G2y0 Þ 2

If

Gz o 0 and G2z 4G2y0 , then 2 If b þ G2y0 o G2z , then it

ð20Þ

is q3 ¼ q3  π . is not possible to achieve the given

longitudinal Gz force. Basic pilot training implies that Gz ¼G (Gx ¼0 and Gy ¼ 0). Consequently, the roll and pitch angles are given by Eqs. (17) and (19).

4. The control algorithm of the centrifuge movement 4.1. Calculation of the angular acceleration q€ 1 Eq. (7) gives the resulting force that is experienced by the simulator pilot at the intersection point of the roll and pitch axes (for q1 ¼ 0) as a function of the angular velocity and acceleration of the centrifuge arm, which is 1 1 2 G ¼ ðG2x0 þ G2y0 þG2z0 Þ1=2 ¼ ða2n þ a2t þ g 2 Þ1=2 ¼ ½a21 ðq_ 41 þ q€ 1 Þ þg 2 1=2 g g ð21Þ According to the requirement that the increase in the acceleration force G should be constant and equal to n, the following is valid: dG d 1 2 _4 2 1=2 €2 ¼ dt ¼ n, g ½a1 ðq1 þ q1 Þ þ g  dt

which yields dð½a21 ðq_ 41 þ q€ 21 Þ þ g2 1=2 Þ ¼ngdt

If we assign the resulting acceleration with a ¼G g, then the previous equation will be: 2

da ¼ dð½a21 ðq_ 41 þ q€ 1 Þ þ g 2 1=2 Þ ¼ ngdt

ð22Þ

The previous differential equation does not have a solution in the general case. In each interpolation period, the robot controller determines the angular velocities of each motor link. An interpolation period of Δt¼ 0.005 s is adopted here. During this period, the servo system of the controller compares (every 0.001 s) the given and achieved motor rotor positions and corrects rotor angular velocities. Based on these observations, an approximated solution from Eq. (22) using a discretisation technique is obtained in the

403

following manner. This approach allows us to solve this differential equation for each interpolation period Δt, which simplifies the solution. For the given rate of change of acceleration Δa=Δt ¼ ng, the acceleration a will first be calculated on the basis of this acceleration in the previous interpolation period, aprev, in the following way: a ¼ aprev þ Δa;

Δa ¼ ngΔt

ð23Þ

If we assign the angular velocity of the centrifuge arm in the previous interpolation period with q_ 1prev , we obtain: q_ 1 ¼ q_ 1prev þ q€ 1 Δt

ð24Þ

If we substitute q_ 1 calculated in this manner into the equation 2 a2 ¼ a21 q_ 41 þ a21 q€ 1 þ g 2 and neglect the terms with Δt3 and Δt4, the following equation for calculating the centrifuge arm acceleration is obtained: q€ 1 ¼

 2q_ 31prev Δt þ ½ð1 þ 6q_ 21prev Δt 2 Þða2  g 2 Þ=a21  2q_ 61prev Δt 2  q_ 41prev 1=2 1 þ 6q_ 1prev Δt 2 2

ð25Þ The previous equation is valid for the movement that has a positive acceleration onset. For the movement that has a negative acceleration onset, the discriminant ð1 þ 6q_ 21prev Δt 2 Þða2 g 2 Þ= a21  2q_ 61prev Δt 2  q_ 41prev is mostly negative, which means that this equation cannot be used directly. In that case, a simple solution is used, in which the values of q€ 1 for the positive acceleration onset n of the same magnitude are reversed. In [28], Eq. (22) is solved, for every interpolation period, using principle given in Eq. (23) as well. Solution is obtained in the form of Jacobi elliptic integrals. 4.2. Smoothing the angular acceleration profile of q€ 1 While a normal acceleration an ¼ a1 q_ 21 and gravity g act all of the time during the centrifuge movement, the tangential acceleration at ¼ a1 q€ 1 does not act during the movement that has a constant load G. The transition from a varying to a constant acceleration causes a sudden increase in q_ 1 , while the transition from a constant to a varying acceleration causes a sudden decrease in q_ 1 . The abrupt change in q_ 1 causes an abrupt change in q2 (which can be seen from Eqs. (16) and (17)) and a very large value of q€ 1 , which causes an abrupt change in q3; these relationships can be seen from Eqs. (18) and (19). To prevent the abrupt change in q_ 1 , it is adopted that before and after the given linear change in G, i.e., acceleration a, there is a period of smoothing of the acceleration curve. The change in the acceleration has three different stages. Within the first period, smoothing is performed in such a way that the acceleration onset n changes linearly from 0 to a given value (snap). Thereafter, the total acceleration change is denoted by ae  as, where as is the initial value of the acceleration, and ae is the desired acceleration value. The linear acceleration change part of the total acceleration change is denoted by cac. This part varies depending on the size of the absolute value of ae as. Here, cac ¼0.12 for abs(ae  as)o1.59 g, and cac ¼0.8 for abs(ae  as) 410.59 g. If 1.59 go abs(ae  as)o 10.59 g, then cac changes linearly from 0.12 to 0.8. Within the smoothing algorithm, it is necessary to determine the following parameters: Na, which is the number of interpolation periods in the starting or ending stage of the G load change; Nc, which is the number of interpolation periods for n ¼ const; nd, which is the desired acceleration rate of

change; and Δn, which is the increase/decrease of n in one interpolation period at the starting or ending stage of the acceleration change. If we assign sign(ae  as)¼ 1 if a increases and sign (ae  as)¼  1 if a decreases, then these parameters will be

404

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

calculated in the following way: N c ¼ cac ðae  as Þ=ðnd g ΔtÞ

ð26Þ

N a ¼ signðae  as Þ½ðae  as Þ=ðnd g ΔtÞ  signðae  as ÞNc   1

ð27Þ

Δn ¼ nd =Na

ð28Þ

n¼

signðae  as Þðae  as Þ  ðN a þN c þ1Þnd g Δt ð2N a þN c Þg Δt

ð29Þ

Even though the acceleration profile is smoothed in transitions from a varying to a constant acceleration a (i.e., G and the angular velocity q_ 1 of the centrifuge arm), in the first interpolation period after the transition, a small tangential acceleration appears. This appearance affects the acceleration load G, which has already been achieved. To nullify these effects, it is necessary to make a small correction in the angular velocity q_ 1 , which gives a new angular acceleration q€ 1 ¼ ðq_ 1  q_ 1prev Þ=Δt of axis 1. The motor of this axis can achieve a desired angular velocity q_ 1 in one interpolation period. However, with regard that angle q3 , depends on q€ 1 , Eqs. (18) and (19), there could be an abrupt pitch rotation of the cockpit. To avoid this rotation, an algorithm is developed that yields a smooth approach to the given G load. 4.3. Calculation of the desired and maximal possible values of q€ 1 , q€ 2 and q€ 3 Inverse kinematics is used to determine the positions, velocities and accelerations of the robot axes by using the desired motion of the end-effector. Based on these kinematic parameters, inverse dynamics gives torques and forces of the robot's axes. This usual way of determination of the robot's axes kinematics and dynamics parameters is sufficient when their actuators can achieve desired motion. On the contrary, in some interpolation periods the control unit will send to the actuators inputs-velocities (positions) which they cannot achieve. Control inputs will be limited to the maximal possible values, but in these interpolations periods the differences between the given and the achieved positions (servo errors) of some links will be too great. Large allowed servo errors reduce the safety during the robot's work. To improve the quality of the centrifuge control algorithm, it is very useful to give, in each interpolation period, realistic (feasible) link positions and to allow only acceptably small servo errors. This will increase the safety of the motion and not lessen the smoothing of the link motions. On the other hand, the computation of the axes forces and moments based on desired but unachievable velocities and accelerations of the robot links would give bigger values of these forces and moments than existed. This will cause over-dimensioning the axes bearings and links that will unnecessary increase the robot mass. For the given reasons, in this paper is suggested a new algorithm for the calculation of the maximal possible values of the axes accelerations q€ 1 , q€ 2 and q€ 3 in each interpolation period. It will limit maximal given axes velocities q_ 1 , q_ 2 and q_ 3 as well. In this way, the real forces and moments that act on each link will be calculated. Because of the possibilities of present actuators, it is very difficult for the motor actuating axis 1 of the centrifuge to achieve complete accuracy given the linear increase in the acceleration force of 9g/s from the lower baseline level of 1.41g to the upper baseline level of, for example, 9g or 15g, and for the motors of axes 2 and 3 to provide that acceleration components Gx and Gy have zero values during this whole acceleration increase. To accomplish such a demanding request, all of the motors should have been more powerful than the motors that are usually used on the

centrifuges. These motors are very heavy and, thus, are unacceptable. It is very important that the roll ring and the gondola, with their motors, are as light as possible. The centrifuge arm motion greatly affects the inertial forces that act on links 2 and 3. The motion of link 2 greatly affects the inertial forces and moments that act not only on link 2 but also on link 3, and vice versa. Conversely, the motions of links 2 and 3 have a very small influence on the inertial forces and the moments that act on the link 1-centrifuge arm because of its size and mass. For this reason, for each interpolation period, the maximal value of q€ 1 can be examined separately, while the maximal values of q€ 2 and q€ 3 must be examined together. 4.3.1. Calculation of the maximal possible value of q€ 1 For axis 1, the DC motor is chosen with a rated torque of M1r ¼ 41,200 Nm, a rated power of P1r ¼1610 kW, a rated number of revolutions of n1r ¼374 min  1, and a maximal number of revolutions through field weakening of n1max ¼680 min  1. If the motor number of revolution is assigned to be n1m, then for n1m rn1r, it will be M1r ¼ const and P1 ¼P1r n1m/n1r, and for n1m 4n1r, it will be M1r ¼P1r n1r/n1m and P1 ¼P1r ¼const. If we assign the short time maximal torque (until 5 s) to be M1max, then the overload capability will be f1p ¼ M1max/M1r ¼2. The chosen gearbox has a gear ratio of k1 ¼ n1m π =ð30q_ 1r Þ ¼ 16:5 and an efficiency of η1 ¼0.94. Now, the rated angular speed of link 1 is q_ 1r ¼ 2:3736 s  1 and maximal angular speed of link 1 q_ 1 max ¼ 4:3157 s  1 . The moment of inertia of the rotor of the motor with the gear box elements brought down on that rotor is I1 ¼400 kgm2.

First, the maximal possible value q€ 1 max of q€ 1 will be calculated. If the angular velocity q_ 1 A ½  q_ 1r ; q_ 1r , then we have the following: if

q€ 1 signðae  as Þ Z 0;

then q€ 1 max ¼ b1 η1

ð30Þ

and if

q€ 1 signðae  as Þ o 0;

then q€ 1 max ¼ b1 =η1

ð31Þ

If the angular velocity q_ 1 A ½  q_ 1max ;  q_ 1r  or q_ 1 A ½q_ 1r ; q_ 1 max , then we have the following: if

q€ 1 signðae  as Þ Z 0;

then q€ 1max ¼ b1 n1r η1 π =ð30q_ 1 Þ

ð32Þ

then q€ 1max ¼ b1 n1r π =ð30q_ 1 η1 Þ

ð33Þ

and if

q€ 1 signðae  as Þ o 0;

If q€ 1 4 q€ 1max , then q€ 1 ¼ q€ 1max , and if q€ 1 o  q€ 1max , then q€ 1 ¼  q€ 1max . Next, we calculate q_ 1 by Eq. (24). If q_ 1 4 q_ 1 max , then q_ 1 ¼ q_ 1 max , and if q_ 1 o  q_ 1max , then q_ 1 ¼  q_ 1max . Now, we have q€ 1 ¼ ðq_ 1  q_ 1prev Þ=Δt and q1 ¼ q1prev þ q_ 1 Δt, with q_ 1prev ¼ q_ 1 , and q1prev ¼ q1 . In Eqs. (29)–(32), we have b1 ¼ M 1r f 1p k1 =I t1 . The total moment of inertia It1 of all three links and the external load reduced to the axis z0 is given by Eq. (A31). 4.3.2. Calculation of the desired and maximal possible values of q€ 2 and q€ 3 Axis 2 and 3 are each actuated by two torque motors. Each of the four motors has a maximal torque of Mimax ¼10,900 Nm, a maximal speed for that torque of q_ i max ¼ 3:979 s  1 , a rated torque of Mir ¼5760 Nm and a maximal speed at that torque of q_ ir ¼ 7:1209 s  1 , i¼2, 3. Their torque decreases linearly from the maximal to the rated value as a function of the speed. Calculation of the desired and maximal possible values of q€ 2 and q€ 3 will be performed in the following six steps. 1. The maximal motor speeds are limited in the following way: if q_ i 4 q_ ir , then q_ i ¼ q_ ir , and if q_ i o  q_ ir , then q_ i ¼  q_ ir , i¼2, 3.

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

2. The maximal values of the motor torques Mia that are related to the motor speed are calculated in the following way: For

q_ i A ½  q_ imax ; q_ imax ; we have M ia ¼ M imax

ð34Þ

for q_ i 4 q_ imax , we have M ia ¼ M imax  ðq_ ir  q_ imax ÞðM imax  M ir Þ= ðq_ ir  q_ imax Þ, and for q_ i o  q_ imax , we have M ia ¼ M imax þ ðq_ ir þ ωimax ÞðMimax  Mir Þ=ðq_ ir  q_ imax Þ, i¼2, 3 3. The accelerations q€ i are limited in the following way. If q€ i 4M ia =I ti for the accelerated motion in the positive direction and the decelerated motion in the negative direction ðq€ i 40Þ, then we have q€ i ¼ M ia =I ti . If q€ i o  M ia =I ti for the decelerated movement in the positive direction and the accelerated movement in the negative direction ðq€ i o 0Þ, then we have q€ i ¼  M ia =I ti . Next, we again calculate q_ i ¼ q_ iprev þ q€ i Δt and _i restrict q_ i in the same way as in Step 1. The new value of ω and qi will be obtained now from q€ i ¼ ðq_ i  q_ iprev Þ=Δt and qi ¼ qiprev þ q_ i Δt, i ¼ 2, 3. The total moment of inertia It2 of the links 2 and 3 and the external load reduced to the axis z1 is given by Eq. (A33). The total moment of inertia It3 of link 3 and the external load reduced to the axis z2 is given by Eq. (A36). 4. A new algorithm is developed for the inverse dynamics of the robots based on the recursive Newton–Euler algorithm, which ^ z2 and m ^ z3 as functions of the calculates the actuator torquesm angular accelerations q€ 2 and q€ 3 and some coefficients on which we do not have influence. This algorithm is applied in _ cm the following way. If in Eq. (A8) for v_ cm 2 , Eq. (A10) for v3 and _vcm 4 , Eq. (A19) for M2 and Eq. (A20) for M3 and M4 , the terms q_ 1 , q_ 2 and q_ 3 are replaced with q_ 1 ¼ q_ 1prev þ q€ 1 Δt, q_ 2 ¼ q_ 2prev þ q€ 2 Δt and q_ 3 ¼ q_ 3rev þ q€ 3 Δt, the actuator torques ^ z3 , given in Eqs. (A26) and (A24), will be obtained in ^ z2 and m m the form of the following two equations: ^ z2 ¼ b2 þ ða21 q€ 1 þ a22 q€ 2 þ a23 q€ 3 ÞΔt m þ ða211 q€ 1 q€ 1 þ a212 q€ 1 q€ 2 þ a213 q€ 1 q€ 3 þ a222 q€ 2 q€ 2 þ a223 q€ 2 q€ 3 þ a233 q€ 3 q€ 3 ÞΔt Δt

ð35Þ

^ z3 ¼ b3 þ ða31 q€ 1 þ a32 q€ 2 þ a33 q€ 3 ÞΔt m þ ða311 q€ 1 q€ 1 þ a312 q€ 1 q€ 2 þ a313 q€ 1 q€ 3 þ a322 q€ 2 q€ 2 þ a323 q€ 2 q€ 3 þ a333 q€ 3 q€ 3 ÞΔt Δt

ð36Þ

In the previous equations, the terms with ΔtΔt can be neglected. Because the motions of links 2 and 3 have a very small influence on the inertial forces and the moments that act on link 1, we could, for each interpolation period, consider that the terms c2 ¼ b2 þ a21 q€ 1 Δt and c3 ¼ b3 þ a31 q€ 1 Δt are known coefficients. Coefficients c22 ¼ a22 Δt, c23 ¼ a23 Δt, c32 ¼ a32 Δt and c33 ¼ a33 Δt can also be introduced. In that ^ z2 and m ^ z3 , given by way, we can obtain the actuator torques m Eqs. (A21), (A24), (A26) and (A28), written in the following, much simpler way: ^ z2 ¼ c2 þ c22 q€ 2 þc23 q€ 3 m

ð37Þ

^ z3 ¼ c3 þ c32 q€ 2 þc33 q€ 3 m

ð38Þ

The coefficients c2, c22 and c23 in Eq. (37) are functions of the variables q€ 1 , q_ 2prev , q1 , q2 and q3 and the constants Δt, a1 , r^ x2 , r^ y2 , r^ z2 , r^ x3 , r^ y3 , r^ z3 , r^ x4 , r^ y4 , r^ z4 , ^I x2 , ^I y2 , ^I z2 , ^I x3 , ^I y3 , ^I z3 , ^I x4 , ^I y4 , ^I z4 , m2, m3 and m4. The coefficients c3, c32 and c33 in Eq. (38) are functions of the variables q€ 1 , q_ 2prev , q_ 3prev , q1 , q2 and q3 and the constants Δt, a1 , r^ x3 , r^ y3 , r^ z3 , r^ x4 , r^ y4 , r^ z4 , ^I x3 , ^I y3 , ^I z3 , ^I x4 , ^I y4 , ^I z4 , m3 and m4. The performed simulations have shown that the ^ z2 and m ^ z3 obtained by Eqs. (37) and numerical result for m (38) differ by less than 1% from the results that are obtained from Eqs. (A21), (A24), (A26) and (A28). This difference is

405

negligibly small. Eqs. (37) and (38) allow a very simple calculation of q€ 2 and q€ 3 , while it is almost impossible to perform that calculation with Eqs. (A21)–(A28). If necessary, the actuator torque for the first axis could be easily added to ^ z1 ¼ c1 þc11 q€ 1 þ c12 q€ 2 þ c13 q€ 3 , and Eqs. (37) and the equation m (38) could be extended with the terms c12 q€ 1 and c13 q€ 1 , respectively. In that way, we would obtain three equations with three unknowns, which we could easily calculate. 5. In this step, check is performed whether the values of required torques of second and third axis for given angular accelerations ^ z2 4 M 2a for the are greater than maximum possible. If m accelerated motion in the positive direction and the decelerated motion in the negative direction ðq€ 2 4 0Þ, then this relationship ^ z2 ¼ M 2a . If m ^ z2 o M 2a for the decelerated motion in will be m the positive direction and the accelerated motion in the negative ^ z2 ¼  M 2a . direction ðq€ 2 o 0Þ, it will be m 6. The maximal possible values of q€ 2 and q€ 3 that motors can achieve will be calculated based on Eqs. (37) and (38) and the ^ z2 and m ^ z3 , as determined in maximal values of the torques m Step 5; this calculation will be performed in the following way: ^ z2  c2 Þc33  ðm ^ z3  c3 Þc23 ðm ; c22 c33  c23 c32 ^ ^ ðmz3  c3 Þc22  ðmz2  c2 Þc32 q€ 3 ¼ c22 c33  c23 c32 q€ 2 ¼

ð39Þ

Next, q_ i and qi are calculated again in the following way: q_ i ¼ q_ iprev þ q€ i Δt; qi ¼ qiprev þ q_ i Δt; q_ iprev ¼ q_ i ; qiprev ¼ qi ; i ¼ 2; 3 4.4. Control algorithm The algorithm for the calculation of q€ 1 , q_ 1 , q1 , q€ 2 , q_ 2 , q2 , q€ 3 , q_ 3 and q3 is based on the equations that are presented in Sections 4.1–4.3. First, Nc, Na, Δn and n are calculated with Eqs. (26)–(29), respectively. Then, the centrifuge kinematic parameters are calculated in three phases. For the positive acceleration onset in the first phase, the onset of the centrifuge inertial force (jerk ¼ Δa/Δt) increases linearly from zero to the programmed value. In the second phase, this force increases linearly (Δa/Δt ¼const¼ n), and in the third stage, the onset of this force decreases linearly from the programmed value to zero. In the first and third phases, smoothing of the acceleration force profile is performed. For the deceleration motion, the smoothing is performed vice versa. These three phases of the centrifuge control algorithm are shown below. First phase. for i¼1 to Na {n ¼ n þ Δn, a ¼ aprev þ Δa ¼ aprev þ signðae  as Þ ng Δt; _ 1 by Eq. (25) for the positive Step 1. Calculation of ω acceleration onset or the reading of the array q€ 1 ¼  arrayðq€ 1 ½2Na þ N c  iÞ for the negative acceleration onset; Step 2. Calculation of the desired and maximal possible _ 1 and the modified values of q_ 1 and q1 (Section value of ω 4.3.1); Step 3. Calculation of q2 by Eqs. (16) or (17) and q3 by Eqs. (18) or (19), respectively, and q_ i ¼ ðqi  qiprev Þ=Δt, q€ i ¼ ðq_ i  q_ iprev Þ=Δt, i ¼2, 3; Step 4. Calculation of the desired and maximal possible _ 2 and ω _ 3 and the corrected values for q_ 2 , q2 , q_ 3 values of ω and q3 (subsection 4.3.2); } Second phase. for i¼ 1 to Nc {a ¼ aprev þ Δa ¼ aprev þsignðae  as Þng Δt; Calculating q€ 1 by Eq. (25) for the positive acceleration onset or reading the array q€ 1 ¼  seriesðq€ 1 ½N a þN c iÞ for the negative

406

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

acceleration onset; The remainder of this phase is the same as in the first phase.} Third phase. for i¼ 1 to Na {n ¼ n  Δn, a ¼ aprev þ Δa ¼ aprev þ signðae  as Þ ng Δt; Calculation of q€ 1 by Eq. (25) for the positive acceleration onset or reading of the array q€ 1 ¼  seriesðq€ 1 ½Na  iÞ for the negative acceleration onset; The remainder of this phase is the same as in the first phase. Steps 2, 3 and 4 are the same as in the first phase.}

5. Results: verification for the proposed control algorithm The proposed control algorithm for the centrifuge motion was tested on the off-line programming system of the robot controller developed at the Lola Institute. Gz, Gx and Gy acceleration force profiles are provided by this programming system. Lola-Industrial Robot Language (L-IRL), based on [29], is here extended with the move instruction and the additional parameters that are required for the centrifuge motion programming. The centrifuge movement control algorithm, presented above, is also added.

Fig. 5. Kinematics and dynamics parameters of the centrifuge motion of the program test 1: (а1) G, (a2) Gx, (a3) Gy, (a4) Gz, (b1) q_ 1 , (b2) q_ 2 , (b3) q_ 3 [1/s], (c1) q€ 1 , (c2) q€ 2 , ^ z1 , (g2) m ^ z2 , (g3) m ^ z3 [Nm], (h1) m ^ xy1 , (h2) m ^ xy2 , (h3) m ^ xy3 [Nm], (i1) P1, (i2) P2 (c3) q€ 3 [1/s2], (d1) q2 ¼ϕ, (d2) q3 ¼ θ [1], (e1) f^ a1 , (e2) f^ a2 , (e3) f^ a3 [N], (f1) f^ r1 , (f2) f^ r2 , (f3) f^ r3 [N], (g1) m and (i3) P3 [kW].

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

Fig. 5 shows kinematics and dynamics parameters of an example of the centrifuge motion program. The masses, the mass centre coordinates and the matrices of inertia of the external load and the centrifuge links are given in Appendix A of this paper. For that example, Fig. 5а1, a2, a3 and a4 shows G, Gx, Gy and Gz, respectively. It is shown that the centrifuge can achieve the mean growth of the acceleration load of 9 g/s in the period of constant acceleration change (cac) in the range between 1.41 and 15 g. The greatest absolute difference between the programmed and realised Gz change is 0.216 g. During the programmed motion, the deviation of Gx was between 0.38 and  0.5 g. There was no deviation in Gy. Fig. 5b1, b2 and b3 shows q_ 1 , q_ 2 and q_ 3 . It can be observed that the presented algorithm limits the values for these velocities: q_ 1 is not larger than 4.28 s  1, q_ 2 is between 1.47 and  1.35 s  1, and q_ 3 is between 3.98 and 3.98 s  1. Fig. 5c1, c2 and c3 shows q€ 1 , q€ 2 and q€ 3 .

Fig. 6. (a) Gz, (b) Gx and (c) q3 ¼ θ from the example given by the program test 1, when the roll and pitch motors are used as well as in the scenarios that have the field and torque weakening.

407

These values are also limited: q€ 1 is between 3.27 and  3.11 s  2, q€ 2 is between 9.88 and  9.4, and q€ 3 is between 52.1 and  52.1 s  2. Fig. 5d1 and d2 shows the angles q2 ¼ ϕ and q3 ¼ θ, which are obtained with the algorithm presented in Section 4.4. Here, q2 is between 0 and 86.21, and q3 is between 43.6 and  57.31. Fig. 5e , e and e shows the axial forces f^ , f^ and f^ and 1

2

3

a1

a2

a3

Fig. 5f1, f2 and f3 shows the radial forces f^ r1 , f^ r2 and f^ r3 of the axes bearings. ^ z1 , m ^ z2 and m ^ z3 of the Fig. 5g1, g2 and g3 shows the torques m actuators. It can be observed that the presented algorithm limits ^ z1 is between 1,228,670 and the values of these torques: m ^ z2 is between 19,664 and  21,777 Nm, and  1,107,000 Nm, m ^ z3 is between 21,756 and  21,696 Nm. Fig. 5h1, h2 and h3 shows m ^ xy1 , m ^ xy2 and m ^ xy3 that act on the bearings. the torques m Fig. 5i1, i2 and i3 shows the powers P1, P2 and P3 of the motors. It can be observed that the presented algorithm limits the values for these velocities; P1 is between 2900 and -2559 kW, P2 is between 10 and  16.417 kW and P3 is between 28.646 and 22.115 kW. ^ z2 and m ^ z3 (Fig. 5g2 and g3) are presented for The diagrams for m the case when the algorithm for the calculation of the maximal possible values of q€ 2 and q€ 3 (Section 4.3.2) has been applied. If this algorithm is not used, the simulation program shows that, to completely accurately achieve a given change in Gz, the required torques that the actuators would have to produce are 26.287 Nm for axis 2 and 24.723 Nm for axis 3. Fig. 5 shows the kinematic and dynamic parameters that were obtained for the given example, in the case of using roll and pitch motors only in the range of the maximal torques. It can be observed that in that case, the chosen motors fulfil all of the requested demands. In contrast, Fig. 6(a) shows Gz, (b) shows Gx, and (c) shows q3 when these motors are used as well in the ranges in which the magnetic fields and torques weaken, until the maximal speeds; in this case, q_ i 4 q_ imax or q_ i o  q_ imax , i ¼2, 3 (in Eq. (34)). It can be observed that in that case, the chosen motors do not satisfy the requested characteristics. During the programmed motion in that case, the deviation of Gx is between 5.55 and  8.82 g, while the greatest difference between the programmed

Fig. 7. A case in which a positive or negative Gz is explicitly given: (a) G, (b) Gz, (c) Gx, (d) q_ 1 , (e) q2 and (f) q3.

408

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

Fig. 8. A case in which a positive or negative Gx is explicitly given: (a) G, (b) Gz, (c) Gx, (d) q_ 1 , (e) q2 and (f) q3.

Fig. 9. A case in which a positive or negative Gy is explicitly given: (a) G, (b) Gz, (c) Gy, (d) q_ 1 , (e) q2 and (f) q3.

and the realised Gz change is 3.529 g. In this way, the implementation of the presented control algorithm shows that using motors for the third axis in the scenario in which the field and torque weakening is included gives unacceptable results. Fig. 7 shows a case in which the positive or negative Gz is explicitly given. Here, G, Gz, Gx, q_ 1 , q2 and q3 are shown in Fig. 7a, b, c, d, e and f, respectively. Fig. 8 shows a case in which a positive or negative Gx is explicitly given. Here, G, Gz, Gx, q_ 1 , q2 and q3 are shown on Fig. 8a, b, c, d, e and f, respectively. Fig. 9 shows a case in which a positive or negative Gy is explicitly given. Here, G, Gz, Gy, q_ 1 , q2 and q3 are shown on Fig. 9a, b, c, d, e and f, respectively.

6. Conclusions In this paper, a new centrifuge motion simulator for pilot training has been presented. This device is modelled as a manipulator that has three rotational axes. The centrifuge produces the desired transverse Gx, the lateral Gy and the longitudinal Gz ^ x , pitch, ω ^ y and yaw, ω ^ z angular acceleration forces and the roll, ω velocities, which simulate the acceleration forces and angular velocities that are accomplished by state-of-the-art aircraft. Because it is difficult to design a centrifuge that would realise all of the tasks completely accurately, it was necessary to develop a

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

new algorithm for the inverse dynamics of the robots; this new algorithm calculates first, in each interpolation period, the successive actuator toques of the links that are required for the given motion. Next, the algorithm checks whether the actuators can achieve those torques, and if they cannot, it calculates the maximal successive link angular accelerations that motors can achieve. The algorithm does not allow giving motion commands to the centrifuge motors that the motors cannot realise. This approach has improved the quality of the motion control and enabled the precise calculation of forces and moments that act on the centrifuge links, which was necessary for the link strength calculations. Using the presented algorithm has enabled us to design the centrifuge links and select its motors in such a way that the centrifuge motion that is realised is notably close to the given motion. Implementation of the presented control algorithm shows that using motors for the third centrifuge axis in the scenario with the field and torque weakening yields unacceptable results.

_ i þ 1  pni þ 1 þ ωi þ 1  ðωi þ 1  pni þ 1 Þ, The linear accelerations v_ i þ 1 ¼ v_ i þ ω i¼ 1, 2, 3 of the centrifuge links are given by Eq. (6). If we assign the position of the link i center of mass (CM) with respect to the link i coordinates expressed in the link i coordinates h iT cm with r^ i ¼ r^ xi r^ yi r^ zi this vector in the centrifuge base

r coordinates is rcm i ¼ ½ xi

0

The forces and moments that act on the centrifuge axes are determined here using the recursive Newton–Euler algorithm. Because the forces and moments that are exerted on link i by link i  1 in link i 1 coordinates do not depend on the position of link 1, the algorithm can calculate them for the angle q1 ¼0. This step reduces the computational cost of the algorithm significantly. For this case, the centrifuge kinematics and dynamics parameters will be as follows. Centrifuge kinematics-velocities and accelerations The centrifuge links angular velocities ωi þ 1 ¼ ωi þ zi q_ i þ 1 , i¼1, 2, 3 (Fig. 2) are: h iT ω1 ¼ 0 0 q_ 1 ; ðA1Þ 



ω2 ¼ ω1 þ 0  1 0 T q_ 2 ; 

ðA2Þ

0

2 3 3 2 3 r z2  r y2  a1  r x2 6 0 7€ 6 7€ 6 7 ¼ 4 a1 þ r x2 5q1 þ 4 5q2 þ 4  r y2 5q_ 1 q_ 1 r x2 0 0 2 3 2 3  r x2 0 6 7 6 7 þ 24 r z2 5q_ 1 q_ 2 þ 4 0 5q_ 2 q_ 2 2

Acknowledgements

Appendix A

r zi T ¼ Di r^ cm (Fig. 10). Here is the i

r yi

external load (weight of the pilot, his seat and equipment) that affects the link 3 marked with 4. Further we can calculate the linear acceleration of the centrifuge links CM with the following equation:h iT cm v_ cm v_ cm v_ cm _ i  rcm v_ cm ¼ v_ i þ ω xi yi zi i ¼ i þ ωi  ðωi  ri Þ. Based on this are: 2 3 2 3  r y1 a1  r x1 6 7€ 6 7  r y1 5q_ 1 q_ 1 v_ cm ðA8Þ 1 ¼ 4 a1 þ r x1 5q1 þ 4

v_ cm 2

This research is supported by the Ministry of Education, Science and Technological Development of Serbia under the project “Development of devices for pilot training and dynamic simulation of modern fighter planes flights: 3DoF centrifuge and 4DoF spatial disorientation trainer” (2011–2014), no. 35023.

409

2

 r y3;4

ðA9Þ

 r z2

0 3

2

 r z3;4

3

2

 s2 r y3;4

3

6 0 7€ 6 7 6 7€ v_ cm 5q2 þ 4 s2 r x3;4  c2 r z3;4 5q€ 3 3;4 ¼ 4 a1 þ r x3;4 5q1 þ 4 c r r 0 2 y3;4 x3;4 2 3 2 3 2 3  r x3;4 a1 þ r x3;4 0 6 0 7_ _ 6 7_ _ 6 7  4 r y3;4 5q1 q1 þ 4 5q2 q2 þ 24 r z3;4 5q_ 1 q_ 2  r z3;4 0 0 2 3 2 3  c2 r y3;4 c2 r z3;4 6 7_ _ 6 7 0 þ 24 s2 r y3;4 5q_ 1 q_ 3 þ 24 5q2 q3  s2 r y3;4 0 2 3 s2 ðc2 r z3;4  s2 r x3;4 Þ 6  c2 r 7_ _ 2 þ4 2 y3;4  s2 r x3;4 5q 3 q3 c2 ðs2 r x3;4 r z3;4 Þ

ðA10Þ

For the centrifuge links shown in this paper, the mass centre coordinates  are: T  T cm cm 0:782 0:002 , r^ 2 ¼ 0: 0:  0:004 , r^ 1  ¼  8:344    T T cm cm r^ 3 ¼ 0: 0: 0:064 and r^ 4 ¼  0:5 0:005 0:55 . Here, we have: h iT ðA11Þ rcm ¼ r^ x1  r^ z1 r^ y1 1



ω3 ¼ ω2 þ c2 0 s2 T q_ 3

ðA3Þ

where zi is the unit vector that represents the zi axis. The centrifuge links linear velocities vi þ 1 ¼ vi þ ωi þ 1  pni þ 1 , i¼1, 2, 3 are:  T ðA4Þ v1 ¼ v2 ¼ v3 ¼ a1 0 1 0 q_ 1    T T where pni þ 1 ¼ pi þ 1  pi , pn1 ¼ a1 0 0 , pn2 ¼ pn3 ¼ 0 0 0 . The centrifuge links angular accelerations _ i þ zi q€ i þ 1 þ ωi  zi q_ i þ 1 , i ¼1, 2, 3 are: _ iþ1 ¼ ω ω h iT _ 1 ¼ 0 0 q€ 1 ; ω ðA5Þ

h c2 r^ z2 s2 r^ x2 rcm 2 ¼

 r^ y2

c2 r^ x2 þs2 r^ z2

3 s2 ðs3 r^ x3;4 þ c3 r^ z3;4 Þ  c2 r^ y3;4 6 7 c3 r^ x3;4 þ s3 r^ z3;4 r3;4 ¼ 4 5 ^ ^ ^  c2 ðs3 r x3;4 þ c3 r z3;4 Þ  s2 r y3;4

iT

ðA12Þ

2

ðA13Þ

In this paper's indexing notation, the subscript 3,4 denotes axis 3 or the external load. Centrifuge dynamics

_ 3 ¼ω _ 2 þ c2 0 s2 T q€ 3 þ 0 c2 0 T q_ 1 q_ 3 þ  c2 0 c2 T q_ 2 q_ 3 ω

The 3  3 matrix on the inertia of link i about the centre of mass of link i (expressed in the base link coordinates) is assigned to Icm i , and the 3  3 matrix on the inertia of link i about the centre of cm mass of link i (expressed in link i coordinates) is assigned to I^ .

ðA7Þ

cm In matrix ^Ii , only the axial moments of the inertia ^I xi , ^I yi and ^I zi









_ 2 ¼ω _ 1 þ 0  1 0 T q€ 1 þ 1 0 0 T q_ 1 q_ 2 ω 







ðA6Þ 



i

410

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

Fig. 10. Coordinate frame and the centre of mass of (a) link 3, (b) link 2 and (c) link 3.

Fig. 11. Forces and moments exerted on (a) link 3, (b) link 2, and (c) link 1.

and m4 ¼ 250 kg. In this way, for axis 1, we have:

are accounted for; in other words, the centrifugal moments of the is a function of the inertia of link i are neglected. Matrix Icm i cm position and orientation of link i, whereas I^ is constant. Matrix

I z1 ¼ ^I y1

Icm i

For axis 2, we have:

i

is calculated with: 2

Icm i ¼

cm Di ^Ii DTi

I xi

6 ¼ 4  I xyi  I xzi

 I xyi I yi  I yzi

 I xzi

2

3

 I yzi 7 5; where I zi

cm I^i

^I xi 6 6 ¼4 0 0

0 ^I yi 0

3

0 7 07 5; ^I zi

2

ðA15Þ

I x2

6 Icm ¼ 4  I xy2 2  I xz2

 I xy2 I y2  I yz2

3

2

^I x2 s2 þ ^I z2 c2 2 2 6 7  Iyz2 5 ¼ 6 0 4 I z2 s c ð^I  ^I Þ  I xz2

2 2

z2

x2

0 ^Iy2

s2 c2 ð^I z2  ^Ix2 Þ

0

^I x2 c2 þ ^I z2 s2 2 2

0

3 7 7 5

ðA16Þ i ¼ 1; 2; 3

ðA14Þ

The calculated axial moments of inertia about the link centres of mass in the link coordinate systems and the masses of the centrifuge links are: ^I x1 ¼ 91; 905, ^I y1 ¼ 219; 978, ^I z1 ¼ 243; 913, ^I x2 ¼ 3; 243, ^I y2 ¼ 1; 365, ^I z2 ¼ 2; 010, ^I x3 ¼ 666, ^I y3 ¼ 217, ^I z3 ¼ 650, ^I x4 ¼ 47, ^I y4 ¼ 61, ^I z4 ¼ 46 kg m2 , m1 ¼45,500, m2 ¼1139, m3 ¼566

and for axis 3 and external load (4), we have: I x3;4 ¼ s22 ð^I x3;4 s23 þ ^I z3;4 c23 Þ þ ^I y3;4 c22 ; I y3;4 ¼  ^I x3;4 c23 þ ^I z3;4 s23 ; I z3;4 ¼ c22 ð^I x3;4 s23 þ ^I z3;4 c23 Þ þ ^I y3;4 s22 ; I xy3;4 ¼ s2 s3 c3 ð^I x3;4  ^I z3;4 Þ I xz3;4 ¼ c2 ðs2 ð^I x3 s3 s3 þ ^I z3;4 c3 c3 Þ  ^I y3 s2 Þ; I yz3;4 ¼ c2 s3 c3 ð  ^I x3;4 þ ^I z3;4 Þ ðA17Þ

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

The total force Fi and the total moment Mi exerted on link i are obtained from the Newton–Euler equation: h iT h cm iT cm cm Fi ¼ F xi F yi F zi ¼ mi v_ xi v_ yi v_ zi  g and h Mi ¼ M xi

M yi

M zi

iT

 I xz2

3

2

 ð  I xy2 Þ

3

2

 ð  I yz2 Þ

3

6 7 6 7 6 7 M2 ¼ 4 I yz2 5q€ 1 þ 4  I y2 5q€ 2 þ 4 ð  I xz2 Þ 5q_ 1 q_ 1  ð  I Þ I z2 0 yz2 2 3 2 3 ð  I yz2 Þ I x2 þ I y2  I z2 6 7_ _ 6 7 0 þ4 5q1 q2  4 ð  I yz2 Þ 5q_ 2 q_ 2

 I xz3;4

3

2

ðA19Þ

I xy3;4

3

2

c2 I x3;4  s2 I xz3;4

3 I xz3;4 s1  I yz3;4 6 7_2 6 7 0 þ 4  I xz3;4 5q_ 21 þ 4 5q 2 I 0 xy3;4 2 3  s2 ðI x3;4  I y3;4 þ I z3;4 Þ 6 2ðs I 7_ _ 2 xy3;4  c2 I yz3;4 Þ 5q þ4 2 q3 2

According to Eq. (A21) this force and moment in link 2 coordinates are:  f y3

^ 3 ¼ ½  s2 mx3 þ c2 mz3 m

c2 f x3 þs2 f z3 T ;

 my3

c2 mx3 þ s2 mz3 T

2

 c2 ðI y3;4  I x3;4 þ I z3;4 Þ

s2 c2 I xy3 ðc22  s22 ÞI yz3;4

The total force and total moment that exert on link 2 are: F2 ¼ f 2  f 3 ,M2 ¼ m2  m3  r2  F2 . Based on this arrangement, the force and moment that exert on link 2 by link 1 with respect to the base coordinate frame are (Fig. 11b):

¼ M2 þ ½ r y2 F z2 r z2 F y2

r z2 F x2  r x2 F z2

r x2 F y2 r y2 F x2 T þ m3 ðA25Þ

The total force and total moment that exert on link 1 are: F1 ¼ f 1  f 2 andM1 ¼ m1  m2  ðpn1 þ r1 Þ  F1  pn1  f 2 . In this way, the force and moment that exert on link 1 with respect to the base coordinate frame are (Fig. 11c): f 1 ¼ F1 þ f 2 ; m1 ¼ M1 þ m2 þðpn1 þ r1 Þ  F1 þpn1  f 2 ¼ M1 2 3 r y1 F z1  r z1 F y1 6  ða1 þ r x1 ÞF z1 þr z1 F x1  a1 f 7 þ4 z2 5 þ m2 ða1 þ r x1 ÞF y1  r y1 F x1 þ a1 f y2

3

6 7 þ 4 s2 c2 ½ðI x3;4  I z3;4 Þ þ ðc22  s22 ÞI xz3;4 5q_ 33 2 c2 I xy3;4  c2 s2 I yz3;4

ðA24Þ

According to Eq. (A21), this force and moment in link 1 coordinates are: h iT h iT ^ 2 ¼ mx2 mz2  my2 f^ 2 ¼ f x2 f z2 f x2  f y2 ; m ðA26Þ

3

3 2 3 2s2 Iyz3;4 Ix3;4 þ Iy3;4  Iz3;4 6 7_ _ 6 7 0 þ4 5q1 q2 þ 4 c2 ðIx3;4 þ Iy3;4  Iz3;4 Þ  2s2 Ixz3;4 5q_ 1 q_ 3  2c2 I yz3;4  2Ixz3;4

3

þ m4

m2 ¼ M2 þ r2  F2 þ m3

2

I yz3;4

iT

f 2 ¼ F2 þ f 3 ;

6 7 6 7 6 7 M3;4 ¼ 4  I yz3;4 5q€ 1 þ 4 I y3;4 5q€ 2 þ 4 c2 I xy3;4  s2 I yz3;4 5q€ 3 I yz3;4 I z3;4 s2 I z3;4  c2 I xz3;4

2

r x3 F y3 r y3 F x3

ð  I xy2 Þ

2ð  I xz2 Þ 2

r z3 F x3  r x3 F z3

ðA23Þ

f^ 3 ¼ ½  s2 f x3 þ c2 f z3

_ i þ ωi  ðIcm ¼ Icm ω ωi Þ i i

The total forces F1, F2 and F3 can be obtained from Eqs. (A8)–(A10), and the total moments are:  T ðA18Þ M1 ¼ 0 0 I z1 q€ 1 2

h þ r y3 F z3 r z3 F y3

411

ðA20Þ

h iT and the We will assign the force vector with f i ¼ f xi f yi f zi h iT moment vector with mi ¼ mxi myi mzi exerted on link i by link i -1 with respect to the base coordinate frame. From robot dynamics, it is known that the total force and the total moment exerted on link i are: Fi ¼ f i  f i þ 1 and Mi ¼ mi mi þ 1  ðpni þrcm i Þ f i þ rcm i  f i þ 1 . Solving for f i and mi gives f i ¼ Fi þ f i þ 1 and n mi ¼ Mi þmi þ 1 þ ðpni þrcm The forces and i Þ  Fi þ pi  f i þ 1 . moments exerted on link i by link i  1 in link i  1 coordinates are: h iT h iT ^ xi m ^ yi m ^ zi ¼ DT mi ^i¼ m f^ i ¼ f^ xi f^ yi f^ zi ¼ DTi 1 f i and m i1

ðA27Þ

Based on Eq. (A21), this force and moment with respect to the base coordinate frame are: h iT h iT ^ x1 m ^ y1 m ^ z1 ¼ mx1 my1 mz1 ; ^1¼ m m h f^ 1 ¼ f^ x1

f^ y1

f^ z1

iT

h ¼ f x1

f y1

f z1

iT

ðA28Þ

^ zi , and the moment acting on The torque of the joint i actuator is m the bearing i is ^ 2xi þ m ^ 2yi Þ1=2 ^ xyi ¼ ðm m

ðA29Þ

The axial force of the bearing i is f^ ai ¼ f^ zi , and radial force of the bearing i is: 2 2 f^ ri ¼ ðf^ xi þ f^ yi Þ1=2

ðA30Þ

ðA21Þ The effect of the external forces and moments on link 3 will be represented by f 4 and m4 (Fig. 11a). In this way, F4 ¼ f 4 and M4 ¼ m4  r4  F4 . The forces and the moments of the external load that exert on link 3 with respect to the base coordinate frame are: f 4 ¼ F4 ; m4 ¼ M4 þr4  F4

ðA22Þ

The total force and the total moment that exert on link 3 are: F3 ¼ f 3  f 4 and M3 ¼ m3 m4  r3  F3 . Based on this arrangement, the total force and the total moment that exert on link 3 by link 2 with respect to the base coordinate frame are: f 3 ¼ F3 þ f 4 ; m3 ¼ M3 þ r3  F3 þ m4 ¼ M3

Total moments of inertia for axes z0, z1 and z2 The total moments of inertia are calculated also for q1 ¼0. The total moment of inertia of all three links and the external load, reduced to the axis z0, is I t1 ¼ I z1 þ I z2 þ I z3 þ I z4 þ ½ð0 r x1 Þ2 þ ð0 r y1 Þ2 þ ð0 r z1 Þ2 m1 þ½ð0 r x2 Þ2 þð0 r y2 Þ2 þ ð0 r z2 Þ2 m2 þ ½ð0 r x3 Þ2 þ ð0 r y3 Þ2 þ ð0 r z3 Þ2 m3 þ½ð0 r x4 Þ2 þ ð0 r y4 Þ2 þð0 r z4 Þ2 m4 þ I 1 k1

ðA31Þ

The values of Iz1, Iz2, Iz3,4 are given in Eqs. (A15)–(A17), respectively; k1 represents the gear ratio, and I1 is the moment of inertia of the motor and gear box of axis 1. The components of the

412

V.M. Kvrgic et al. / Robotics and Computer-Integrated Manufacturing 30 (2014) 399–412

vectors0 r1 , 0 r2 , 0 r3 and 0 r4 are:   0 r1 ¼ 0 r x1 0 r y1 0 r z1 T ¼ D1 ð a1 0  0

r2 ¼

0

r3 ¼

0

r4 ¼

 T r x2 0 r y2 0 r z2  ¼ D1 a1

0

0

0

 T r x3 0 r y3 0 r z3  ¼ D1 a1

0

0

0

 T r x4 0 r y4 0 r z4  ¼ D1 a1

0

0

 

h

þ r^ 1 Þ ¼ r x1 þ a1

r y1

r z1

T

h þ D2 r^ 2 ¼ r x2 þ a1

r y2

r z2

T

h þ D3 r^ 3 ¼ r x3 þ a1

r y3

r z3

T

h þ D3 r^ 4 ¼ r x4 þ a1

r y4

r z4

0

0

T

iT iT iT iT

ðA32Þ In previous equations, the matrices D1 , D2 and D3 were given in Eqs. (1), (4) and (5), and the vectors r1 , r2 and r3;4 were given in Eqs. (A11)–(A13), respectively. The total moment of inertia of links 2 and 3 and the external load reduced to axis z1 is I t2 ¼ 1 I z2 þ 1 I z3 þ 1 I z4 þ ½ðr^ x2 Þ2 þ ðr^ z2 Þ2 m2 þ ½ð1 r x3 Þ2 þ ð1 r y3 Þ2 þ ð1 r z3 Þ2 m3 þ ½ð1 r x4 Þ2 þ ð1 r y4 Þ2 þ ð1 r z4 Þ2 m4 ðA33Þ In 1 cm I2 1 cm I3

(A33), the value 1 I z2 ¼ ^I y2 is obtained from cm ¼ D2 I^2 1 D2 1 , the value 1 I z3 ¼ ^I x3 c23 þ ^I z3 s23 is obtained from cm ¼ 1 D3 I^ 1 D  1 , and the value 1 I z4 ¼ ^I x4 c2 þ ^I z4 s2 is obtained

Eq.

1

3

3 cm

3

3

from 1 Icm ¼ 1 D3 I^4 1 D3 1 . Eq. (A33) applies the following: 4 2 3 2 3 s2 s3 r^ x3  c2 r^ y3 þ s2 c3 r^ z3 61 1 1 7 1 6 1 7 ^ ^ ^ 7 r3 ¼ 6 4 r x3 r y3 r z3 5 ¼ D3 r^ 3 ¼ 4  s3 c3 r x3  s2 r y3  c2 c3 r z3 5 c3 r^ x3 s3 r^ z3 ðA34Þ 2

1

3

2 3 s2 s3 r^ x4  c2 r^ y4 þ s2 c3 r^ z4 61 1 1 7 1 6 7 ^ ^ ^ 7 r4 ¼ 6 4 r x4 r y4 r z4 5 ¼ D3 r^ 4 ¼ 4  s3 c3 r x4  s2 r y4  c2 c3 r z4 5 c3 r^ x4 s3 r^ z4 ðA35Þ 1

1

In previous equations, we have D3 ¼ D2 D3 . The matrices 1 D2 2

and 2 D3 are given in Eqs. (2) and (3), respectively. The total moment of inertia of link 3 and the external load reduced to the axis z2 is: I t3 ¼ ^I y3 þ ^I y4 þ ½ðr^ x3 Þ2 þ ðr^ z3 Þ2 m3 þ ½ðr^ x4 Þ2 þ ðr^ z4 Þ2 m4 In 2 cm I3 2 cm I4

ðA36Þ

(A36), the value 2 I z3 ¼ ^I y3 is obtained from cm ¼ D3 I^3 2 D3 1 , and the value 2 I z4 ¼ ^I y4 is obtained from cm ¼ 2 D3 I^ 2 D  1 .

Eq.

2

4

3

References [1] Albery WB. Human centrifuges: the old and the new. In: Proceedings of the 36th annual symposium of the SAFE association, Phoenix AZ; 14–16 September 1998.

[2] Albery WB. Current and future trends in human centrifuge development. SAFE J 1999;29(2). [3] Albery WB. Acceleration in other axes affects þGz tolerance: dynamic centrifuge simulation of agile flight. Aviat Space Environ Med 2004;75(1):1–6. [4] Schot SH. Jerk: the time rate of change acceleration. Am J Phys 1978;46 (11):1090–4. [5] Gallardo-Alvarado J. Jerk analysis of a six-degrees-of-freedom three-legged parallel manipulator. Robotics Comput-Integr Manuf 2012;28:220–6. [6] Morasso P. Spatial control arm movements. Exp Brain Res 1981;42:223–7. [7] Flash T, Hogan N. The coordination of arm movements: an experimentally confirmed mathematical model. J Neurosci 1985;5:1688–703. [8] Uno Y, Kawato M, Suzuki R. Formation and control of optimal trajectory in human multi joint arm movements. Biol Cybern 1989;61:89–101. [9] Viviani P, Schneider R. A developmental study of the relationship between geometry and kinematics in drawing movements. J Exp Psychol: Hum Percept Perform 1991;17(1):198–218. [10] Viviani P, Flash T. Minimum-jerk, two-third spower, law and isochrony: converging approaches to movement planning. J Exp Psychol: Hum Percept Perform 1995;21(1):32–53. [11] Crosbie, RJ. Dynamic flight simulator control system. United States Patent, Number 4,751,662; Jun. 14, 1988. [12] Crosbie, RJ, Kiefer, DA. Controlling the human centrifuge as a force and motion platform for the dynamic flight simulator technologies. In: Proceedings of the AIAA flight simulation conference, St Louis, MO, AIAA Paper; 1985. p. 37–45. [13] Tsai M.-H., Shih M.-C. G-load tracking control of a centrifuge driven by servo hydraulic systems. Proc IMechE G: J Aerosp Eng, 2009; Vol. 223pp. 669–682. [14] Chen YC, Reperger DW. A study of the kinematics dynamics and control algorithms for a centrifuge motion simulator. Mechatronics 1996;6(7):829–52. [15] Chen YC, Reperger DW, Roberts R. Study of the kinematics, dynamics and control algorithms for a centrifuge motion simulator. In: Proceedings of the American control conference, vol. 3; 1995. p. 1901–1905. [16] Tsai L-W. Robot analysis: the mechanics of serial and parallel manipulators. John Wiley and Sons, Inc.; 1999. [17] Nof SY. Handbook of industrial robotics. John Wiley and Sons, Inc.; 1985. [18] Wu J, Wang JS, You Z. An overview of dynamic parameter identification of robots. Robotics Comput-Integr Manuf 2010;26(5):414–9. [19] Wu J, Wang J, Wang L, Li T. Dynamics and control of a planar 3-DOF parallel manipulator with actuation redundancy. Mech Mach Theory 2009;44 (4):835–49. [20] Wu J, Chen X, Li T, Wang L. Optimal design of a 2-DOF parallel manipulator with actuation redundancy considering kinematics and natural frequency. Robotics Comput Integr Manuf 2013;29(1):80–5. [21] Siciliano B, Sciavicco L, Villani L, Oriolo G. Robotics: modelling, planning and control. Springer-Verlag London Limited; 2009. [22] Djuric A, Saidi RA, ElMaraghy W. Dynamics solution of n-DOF global machinery model. Robot Comput-Integr Manuf 2012;28:621–30. [23] Tsai L-W. Robot analysis the mechanics of serial and parallel manipulator. New York: Wiley; 1999. [24] Staicu S, Liu XJ, Li J. Explicit dynamics equations of the constrained robotic systems. Nonlinear Dyn 2009;58:217–35, http://dx.doi.org/10.1007/s11071009-9473-4. [25] Gherman B, Pisla D, Vaida C, Plitea N. Development of inverse dynamic model for a surgical hybrid parallel robot with equivalent lumped masses. Robotics Comput-Integr Manuf 2012;28:402–15. [26] Lutovac M, Kvrgic V, Ferenc G, Dimic Z, Vidakovic J. 3D simulator for human centrifuge motion testing and verification. In: Proceedings of the 2nd Mediterranean conference on embedded computing MECO; 2013. p. 160–163. [27] Renato Vidoni R, Gasparetto A, Giovagnoni M. Design and implementation of an ERLS-based 3-D dynamic formulation for flexible-link robots. Robotics Comput-Integr Manuf 2013;29:273–82. [28] Vidakovic J, Ferenc G, Lutovac M, Kvrgic V. Development and implementation of an algorithm for calculating angular velocity of main arm of human centrifuge source. In: Proceedings of the 15th international power electronics and motion conference, EPE-ECCE 2012 Europe Congress; 2012. pp. DS2a.17. [29] DIN 66312 Teil 1, Industrieroboter Programmiersprachen, Entwurf; 1992.