Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 57, 1, pp. 99-113, Warsaw 2019 DOI: 10.15632/jtam-pl.57.1.99 AN INVERSE KINEMATIC MODEL OF THE HUMAN TRAINING CENTRIFUGE MOTION SIMULATOR Rafał Lewkowicz Military Institute of Aviation Medicine, Warsaw, Poland e-mail: rlewkowicz@wiml.waw.pl Grzegorz Kowaleczko Air Force Institute of Technology, Warsaw, Poland e-mail: grzegorz.kowaleczko@itwl.pl The paper presents an inverse kinematic model for a centrifuge motion simulator used to verify newly defined absolute accelerationprofiles.Themodelling is concernedwith a human training centrifuge with three degrees of freedom. The values of kinematic parameters have beenobtained for this three-jointedmanipulator.Validationof thedevelopedmodel hasbeen performed by comparing the results obtained from the centrifugemotion simulator with the results of numerical simulations. The simulation revealed that the inverse kinematic model enabled calculation of the angular displacement, velocity and acceleration of the links that are needed for the given linear acceleration of the simulator cabin. Keywords: inverse kinematic, simulator, centrifuge, motion system 1. Introduction High-performance aircraft pilots as well as civilian aerobatic pilots are exposed to high linear accelerations duringflight (Newman, 2015). In order toproperlyprepare forwork in this environ- ment, pilots are evaluated and trained to increase their acceleration tolerance level (Wojtkowiak, 1991). This training is carried out in a specially designed centrifuge motion simulator or, so- called, human training centrifuge (HTC) (Dančuo et al., 2012b; Truszczyński and Kowalczuk, 2012). From the standpoint of classical mechanics, the task of the HTC is to achieve accelera- tions through rotations around three axes that simulate the load the pilot is exposed to in a real flight. This type of simulatormakes it possible to create high and prolonged linear accelerations. Moreover, the centrifuge provides a safe ground-based platform to train pilots, especially in the field of increasing pilot’s awareness about unwanted effects of accelerations, such as G-induced loss of consciousness or spatial disorientation. In addition, the HTC is an important tool for researchers to understand the changes taking place in human physiology during accelerative stress. In this paper, a dynamic flight simulator has been considered, namely the HTC shown in Fig. 1, manufactured by the AMST-Systemtechnik GmbH (Austria), located at the Military Institute of Aviation Medicine (Warsaw, Poland). The gondola/cabin of the centrifuge is assembled on an eight-meter-long arm and allowed longitudinal accelerations (in the direction frompilot’s head to foot) to be achieved in the range from −3g to +16g (g is the Earth’s gravitational acceleration) with the maximal onset rate of acceleration (G-onset rate) at n = 14.5g/s. Additionally, the gyroscopic suspension of the cabin allowed it to achieve transversal and lateral accelerations in the range of ±10g and ±6g, respectively. 100 R. Lewkowicz, G. Kowaleczko Fig. 1. Dynamic flight simulator – human training centrifuge HTC-07 2. Problem formulation Within the centrifuge simulator, there is a pre-programmed standard open loopmode that con- tains predefined nonlinear profiles of the absolute acceleration in the centre of the cabin. In these profiles, which are independent for each of the three axes of the pilot’s head-fixed coor- dinate system, a positive or negative acceleration that is constant for given periods of time is determined. When this acceleration profile is created, the angular accelerations in the indivi- dual links of the simulator motion system may be exceeded (Table 1). There is also a problem concerning hypogravity (< 1g), which cannot be obtained in this simulator. In that case, the simulator software reports an error, which is then eliminated using a trial and error method. This problem makes it difficult to definemore complex acceleration profiles that should not be tested on a real device. The solution to this problem should be sought in the inverse kinematics of the HTC’s motion system. The inverse kinematics of the centrifuge motion simulator will be based on calculating the angular position, velocity and acceleration of each motion system link. This approach will make it possible to indicate maximum values of angular acceleration that are necessary to achieve a given linear acceleration of the simulator cabin. In this way, how the position of individual links of the motion system should change over time in order to achieve the desiredmovement of the cabin will be determined. 3. Physical model of the HTC simulator A centrifuge motion simulator is modelled as a three-joint manipulator (Fig. 2) with rotational axes, where the pilot’s head is considered to be the end-effector (Crosbie, 1988). The model consists of three links: arm, ring and cabin (Fig. 2). The arm rotation around the vertical axis is the main motion that achieves the desired acceleration force. The arm carries a gimballed cabin systemwith two rotational axes providing pitch and roll capabilities. The task of the roll and pitch axes is to direct the acceleration force into the desired direction. The pilot’s head is placed in the intersection of the cabin roll and pitch axes. The arm rotation angle is denoted by ψA, the roll ring rotation angle by ϕR and the pitch cabin rotation angle by θC. The centrifuge has the following parameters: arm length dA =8m, roll and pitch axis rotation range of±360◦. Other parameters of the simulatormotion system kinematics are shown in Table 1. An inverse kinematic model of the human training centrifuge motion simulator 101 Fig. 2. A physical model of the centrifuge motion simulator HTC-07 Table 1.Motion capabilities of the HTC-07 simulator (AMST-Systemtechnik GmbH, 2011) Parameter zc-axis yc-axis xc-axis Maximum acceleration rate [g/s] 14.5 6 10 Maximum angular acceleration [rad/s2] 2.82 8 5 4. Kinematics of the centrifuge motion simulator The forward kinematics that is related to the simulator motion system geometry is used to calculate the linear acceleration components of the end-effector (the pilot’s head), whe- re GH = [GxH ,GyH,GzH ] T with respect to the centrifuge variables ψA, φR and θC (the- se angles are determined in the next Section). Thus, for a given joint coordinate vector q= [ψ̈A, ψ̇A,θC, θ̈C,φR, φ̈R] T, the forward kinematics equation must be solved as follows GH = f(q) (4.1) where f is a nonlinear, continuous and differentiable function. The simulator kinematic model can be derived by different methods, such as the Lagrange equation (Siciliano et al., 2009; Wu et al., 2010), the Newton-Euler method (Grotjahn et al., 2004; Tsai, 1999) and the virtual work principle (Wu et al., 2009, 2013; Zhao and Gao, 2009). The method based on the Lagrange formulation is conceptually simple and systematic. Themethod based on theNewton-Euler for- mulation yields amodel in a recursive form. It is composed of forward computation of velocities and accelerations of each link, followed by backward computation of forces and moments in each joint (Wu et al., 2010). This algorithm is computationally more efficient because it exploits the typically open structure of the manipulator kinematic chain (Djuric et al., 2012; Sicilia- no et al., 2009). On the other hand, the Newton-Euler procedure is very difficult to use in an advanced control application because of the closed structure, as the expense of calculation is considerably high (Gherman et al., 2012). Despite this, the Newton-Euler equations of motion 102 R. Lewkowicz, G. Kowaleczko are used tomodel the kinematics of centrifuge motion simulator due to the fact that these equ- ations incorporate all accelerations that act on the individual links of themotion system (Chen and Repperger, 1996; Dančuo et al., 2012a; Kvrgic et al., 2014; Vidaković et al., 2012). During kinematic modelling of the centrifuge, the small elastic deformation of the centrifuge links is neglected. 4.1. Coordinate frames and matrices determining relations for centrifuge links This Section defines coordinate frames for the centrifuge links (Fig. 2) and matrices that determine their relations. The centrifuge links and their coordinate frames are denoted by using theEuler angle convention. The centrifuge base coordinates are denoted by Ox0y0z0 (theEarth- -fixed system), the arm coordinates by AxAyAzA (link 1), the roll ring coordinates by RxRyRzR (link 2), the cabin coordinates by CxCyCzC (link 3) and the pilot head-fixed coordinates system by HxHyHzH. The pilot’s head is placed in the intersection of the cabin roll and pitch axes, therefore xC = xH, yC = yH, and zC = zH. To determine the mutual position of the defined coordinate systems, the following angles are used (Fig. 2): • ψA – the yaw angle between the x0-axis and the yA-axis. This angle, enlarged by 90◦, provides coverage of the y0-axis with the yA-axis, in this way defining the position of the Earth-fixed coordinate system Ox0y0z0 relative to the arm-fixed coordinate system AxAyAzA (link 1), • φR – the roll angle between the zA-axis and the zP-axis. This angle determines the position of the arm-fixed coordinate system AxAyAzA (link 1) relative to the ring-fixed coordinate system RxRyRzR (link 2), • θC – the pitch angle between the zP-axis and the zC-axis. This angle determines the position of the ring-fixed coordinate system RxRyRzR (link 2) relative to the cabin-fixed coordinate system CxCyCzC (link 3). To derive the kinematic equations for the motion simulator, the matrices for the relation between the centrifuge link coordinate frames are defined. These transformation matrices are obtained for the Z → Y → X rotation convention of the coordinate systems in the following way: — L0/A – matrix for transformation from the Earth-fixed coordinate system Ox0y0z0 to the arm-fixed coordinate system AxAyAzA L0/A =Lzr(90 ◦)Lz0(−ψA) (4.2) where Lz0(−ψA)=    cosψA −sinψA 0 sinψA cosψA 0 0 0 1    Lzr(90 ◦)=    0 −1 0 1 0 0 0 0 1    (4.3) —Lxr(φR) – matrix for transformation from the arm-fixed coordinate system AxAyAzA to the ring-fixed coordinate system RxRyRzR Lxr(φR)=    1 0 0 0 cosφR sinφR 0 −sinφR cosφR    (4.4) — Lyp(θC) – matrix for transformation from the ring-fixed coordinate system RxRyRzR to the cabin-fixed coordinate system CxCyCzC Lyp(θC)=    cosθC 0 −sinθC 0 1 0 sinθC 0 cosθC    (4.5) An inverse kinematic model of the human training centrifuge motion simulator 103 The matrix L0/C for transformation from the cabin-fixed coordinate system, CxCyCzC, to the Earth-fixed coordinate system, Ox0y0z0, is determined by multiplying transformation matrices (4.2)-(4.5), in the following way L0/C =Lyp(θC)Lxr(φR)L0/A (4.6) By using the convenient shorthand notation c =cos and s = sin, the components of the trans- formation matrix (4.6) L0/A =Lzr(90 ◦)Lz0(−ψA) (4.2) become the following L0/C =    sψAcθC + cψAsφRsθC cψAcφR −sψAsθC + cψAsφRcθC −cψAcθC +sψAsφRsθC sψAcφR cψAsθC +sψAsφRcθC cφRsθC −sφR cφRcθC    (4.7) Assuming that ψA =0, for further calculations matrix (4.7) is reduced to the form L0/C(ψA=0) =    sinφR sinθC cosφR sinφRcosθC −cosθC 0 sinθC cosφR sinθC −sinφR cosφRcosθC    (4.8) On the basis of these transformationalmatrices, the equations of forward kinematics that relate to the velocities and accelerations of the links and the end-effector-pilot’s head are developed in the next Section. 5. Linear acceleration acting on the pilot’s head The linear acceleration components at the intersection point of the roll (link 2) andpitch (link 3) axes are: the normal (radial) an, tangential at and gravitational g accelerations (Fig. 3). From these accelerations, the orthogonal components Gx0, Gy0, and Gz0 for the normal, tangential and vertical accelerations, respectively, are calculated as follows    Gx0 Gy0 Gz0    = 1 g    −an −at g    =    dAψ̇ 2 A/g −dAψ̈A/g 1    (5.1) where dA is the simulator arm length, ψ̇A and ψ̈A are angular velocity and acceleration of the arm (link 1), respectively. The link angles φR and θC, angular velocity ψ̇A and acceleration ψ̈A of the arm define the orthogonal components GxC, GyC and GzC of the resultant vector GC that are experienced by the pilot. Based on equations (4.8) and (5.1), the resultant vector GC experienced by the pilot can be found from [GxC ,GyC,GzC ] T =L−1 0/C(ψA=0) [Gx0,Gy0,Gz0] T (5.2) The transverse GxC, lateral GyC and longitudinal GzC components of the acceleration GC that act on the pilot’s head are GxC =sinθC(Gx0 sinφR +Gz0cosφR)−Gy0 cosθC GyC = Gx0cosφR −Gz0 sinφR GzC =cosθC(Gx0 sinφR +Gz0cosφR)+Gy0 sinθC (5.3) Equations (5.3) develop the inverse kinematics algorithmwhich determines the link angles that are required to generate a desired trajectory of the cabin-centrifuge and accelerations in the centrifuge axes. 104 R. Lewkowicz, G. Kowaleczko Fig. 3. The acceleration components at the intersection point of the roll and pitch axes 6. An inverse kinematics of the cabin simulator The inverse kinematics is defined as the problem of determining a set of appropriate joint configurations forwhich theendeffector (thepilot’s head)moves todesiredpositionsas smoothly, rapidly and as accurately as possible. The inverse kinematics of the centrifuge motion system is first based on calculating the angular displacement, velocity and acceleration of the links that are needed for the given linear acceleration of the simulator cabin. Then, taking into account the limitations of the motion system (Table 1), it is checked whether this system can achieve such accelerations. If it cannot, the maximum successive link angular accelerations that the motion system can achieve are calculated. The inverse kinematics for the centrifuge motion simulator can be described by the relationship [ψ̈A, ψ̇A,θC, θ̈C,φR, φ̈R] T = f−1[GxH,GyH ,GzH] T (6.1) where f−1 is a nonlinear, continuous and differentiable function that performs the inverse trans- formation to the function f, (4.1). There are twodistinctmethods for solvingEq. (6.1) of inverse kinematics, namely iterative and analytical. The iterative method gives the solution by solving an approximation of the system, and by updating the system with the output from the solver for each iteration until it converges. The analytical method solves the whole system at once; however, the complexity of it arises when large chains of joints attempt to be solved. Themost prominent among iterativemethods are based on the Jacobianmatrix, which describes the non- linear and configuration dependent transformation between velocities in the joint configuration coordinates and the task spaces. There are several versions of Jacobian-based methods, such as An inverse kinematic model of the human training centrifuge motion simulator 105 the Jacobian Transpose (Wolovich and Elliott, 1984), damped least squares (Wampler, 1986), damped least squares with singular value decomposition (Wampler, 1986), selectively damped least squares (Buss and Kim, 2005) and several extensions (Baillieul, 1985; Nakamura and Ha- nafusa, 1986). An iterative method can be also viewed as an optimization task solved with general-purposemethods (neural networks (Tejomurtula andKak, 1999) and genetic algorithms (Nearchou, 1998)), but those approaches are usually computationally ineffective. Due to small chains of joints (a centrifuge simulator is a three-jointmanipulator), direct and analytical computations of inverse kinematics are chosen. Based on the acceleration vector GC components, Eqs. (5.3), the link angles in each joint in the system are derived. To determine the angular accelerations and velocities of the arm (link 1), ring (link 2) and cabin (link 3), the following calculation algorithm for the inverse kinematics is used (Kvrgic et al., 2014). Step 1: Determination of the arm angular acceleration ψ̇A. The angular acceleration ψ̇A is derived fromEq. (5.1) that describes linear acceleration components at the intersection point of the roll andpitch axes (Fig. 3).The resultant acceleration of this point is a sumof thenormalan, tangential at, and gravitational g acceleration. This acceleration is as follows a2A = d 2 Aψ̇ 4 A +d 2 Aψ̈ 2 A +g 2 (6.2) For a positive angular acceleration ψ̈A, the angular velocity ψ̇A in the i-th moment of time is equal to ψ̇A(i)= ψ̇A(i−1)+ ψ̈A(i)dt (6.3) By substituting (6.3) to equation (6.2), the resultant acceleration takes the form a2A(i) = d 2 A[ψ̇A(i−1)+ ψ̈A(i)dt] 4+d2Aψ̈ 2 A(i)+g 2 (6.4) After calculations have been performed, this equation becomes a2A(i) = d 2 A[ψ̇ 4 A(i−1)+4ψ̇ 3 A(i−1)ψ̈A(i)dt+6ψ̇ 2 A(i−1)ψ̈ 2 A(i)dt 2]+d2Aψ̈ 2 A(i)+g 2 (6.5) and then a2A(i)−g 2 d2A = ψ̇4A(i−1)+4ψ̇ 3 A(i−1)ψ̈A(i)dt+6ψ̇ 2 A(i−1)ψ̈ 2 A(i)dt 2+ ψ̈2A(i) (6.6) By reducing equation (6.6) to the form of a quadratic equation [1+6ψ̇2A(i−1)dt 2]ψ̈2A(i)+4ψ̇ 3 A(i−1)ψ̈A(i)dt+ ψ̇ 4 A(i−1)− a2A(i)−g 2 d2A =0 (6.7) it is possible to obtain its solution in the form of two roots ψ̈A(i)= −2ψ̇3A(i−1) dt± √ −2ψ̇6A(i−1)dt2− ψ̇ 4 A(i−1)+ [1+6ψ̇ 2 A(i−1)dt2]k(i) 1+6ψ̇2A(i−1)dt2 (6.8) where k(i) = a2A(i)−g 2 d2A (6.9) Kvrgic et al. (2014) noted that equation (6.8) is valid for the movement that has a positive acceleration onset. For a negative acceleration onset, the discriminant−2ψ̇6A(i−1)dt 2− ψ̇4A(i− 1)+ [1+6ψ̇2A(i−1)dt 2]k(i) is mostly negative, which means that this equation cannot be used 106 R. Lewkowicz, G. Kowaleczko directly.Vidaković et al. (2012, 2013) proposeda solution in the formof aJacobi elliptic function, which describes the arm angular velocity as ψ̇A(t)= 4 √ ksn( 4 √ kt+ 4 √ kC1,−1) (6.10) where k is constant for every interpolation period of time and is given by Eq. (6.9), sn is a Jacobian elliptic function and C1 is the constant obtained from the value of angular velocity from the previous interpolation period. After equation (6.10) has been developed in Taylor series expansions of the Jacobi elliptic function (Wrigge, 1981), it becomes ψ̇A(i)= 4 √ k(i) ( t1(i)− t51(i) 10 + t91(i) 120 − 11t131 (i) 15600 + 211t171 (i) 3536000 ) (6.11) where t1(i)= 4 √ k(i)(dt+C1). Equation (6.11) describes the arm angular velocity, ψ̇A(i), for each i-th interpolation period of time. The angular acceleration, ψ̈A(i), of the arm for every interpolation period of time is calculated as ψ̈A(i)= ψ̇A(i+1)−dotψA(i) dt (6.12) Another approach to solve the problem of calculation of the negative acceleration was proposed by Liwen et al. (2015). The researchers generated both a trapezoidal G-load curve and three- -axis G-load commands using a real-time motion planning algorithm with two G-dimensional interpolation. Dančuo et al. (2018) and Vidaković et al. (2012) indicated that equation (6.8) could be also solved numerically for a small time increment dt → 0. Step 2:Determination of the angular velocity ψ̇A of the arm(6.3) and accelerations components (5.1) at the intersection point of the roll and pitch axes (Fig. 3). Step 3: Determination of the roll ring angle φR based on equation (5.3)2, which describes the lateral acceleration GyC. Expressing sinφR and cosφR by using the tangent function cosφR = 1 √ 1+tan2φR sinφR = tanφR √ 1+tan2φR (6.13) and by substituting these functions to equation (5.3)2, multiplying both sides by √ 1+tan2φR, and then raising to the power, the following expression is obtained G2x0−2Gx0Gz0 tanφR +G 2 z0 tan 2φR = G 2 yC (1+tan2φR) (6.14) After the next transformation, Eq. (6.14) is reduced to the form (G2z0−G 2 yC )tan2φR −2Gx0Gz0 tanφR +G2x0−G 2 yC =0 (6.15) for which the roots are the following tanφR = 2Gx0± √ 4G2x0−4(G 2 z0−G2yC)(G 2 x0−G2yC) 2(G2z0−G2yC) (6.16) By substituting Gz0 =1 to Eq. (6.16), and then performingmanipulations, the following result is obtained tanφR = Gx0±GyC √ 1+G2x0−G2yC 1−G2yC (6.17) An inverse kinematic model of the human training centrifuge motion simulator 107 For G2x0+1­ G 2 yC , the roll ring angle is equal to φR =arctan Gx0+GyC √ 1−G2yC +G 2 x0 1−G2yC (6.18) otherwise φR =arctan Gx0 1−G2yC (6.19) If GyC < 0, and G 2 yC > 1, the roll ring angle is equal to φR = φR +π. The angular velocity φ̇R and acceleration φ̈R of the ring (link 2) are determined as follows φ̇R(i) = φR(i)−φR(i−1) dt φ̈R(i)= φ̇R(i)− φ̇R(i−1) dt (6.20) Step 4: Completion of the pitch cabin angle θC calculation. This angle can be derived from equation (5.3)1 that describes the lateral acceleration GxC or based on Eq. (5.3)3 which defines the longitudinal acceleration GzC. Equations (5.3)1 and (5.3)3 indicate that it is not possible to obtain simultaneously the desired values of GxC and GzC acceleration, even if they do not exceed the limit ranges (Table 1). Therefore, to determine the pitch cabin angle, Eq. (5.3)1 is used. For the known lateral acceleration GxC, using similar substitution (6.13) for the pitch cabin angle θC, equation (5.3)1 takes the form GxC = tanθC √ 1+tan2θC (Gx0 sinφR +Gz0cosφR)−Gy0 1 √ 1+tan2θC (6.21) Performing analogous transformations (6.14) and (6.15) as for the angle of tilting the ring, the above equation takes the form G2xC(1+tan 2θC)= tan 2θC(Gx0 sinφR +Gz0cosφR) 2 −2Gy0 tanθC(Gx0 sinφR +Gz0cosφR)+G2y0 (6.22) By substituting Gx0 sinφR +Gz0cosφR = d, and then performing some manipulations, a qua- dratic equation is obtained (d2−G2xC)tan 2θC −2Gy0dtanθC +G2y0+G 2 xC =0 (6.23) for which the roots are the following tanθC = Gy0d±GxC √ d2−G2y0−G2xC d2−G2xC (6.24) For (d2+G2y0)­ G 2 xC , the pitch cabin angle θC is equal to θC =arctan Gy0d+GxC √ d2+G2y0−G2xC d2−G2xC (6.25) otherwise, if (d2+G2y0) < G 2 xC , it isnotpossible toobtain thedesired transverseacceleration GxC. Then equation (6.21) takes the form θC =arctan Gy0 d−G2xC/d (6.26) 108 R. Lewkowicz, G. Kowaleczko The angular velocity θ̇C and acceleration of the cabin θ̈C (link 3) are determined as follows θ̇C(i)= θC(i)−θC(i−1) dt θ̈C(i)= θ̇C(i)− θ̇C(i−1) dt (6.27) A centrifuge is capable of simulating all three load components, but simulating a pure GzC profile becomes a problem due to motor limitation. A pure GzC training profile is a profile without the transverse GxC and lateral GyC loads. The transverse load GxC in the centrifuge is a result of the tangential acceleration at (Fig. 3) and has a large value at the beginning and at the end of planetary arm motion. The greater the tangential acceleration is, the greater is the GxC acceleration. The emergence of the large tangential acceleration at decreases the angular velocity of the arm motion and has negative effects on the overall centrifuge perfor- mance (Dančuo et al., 2013). These effects are minimized by adjusting the pitch cabin angle by θC =arctan(at/ √ a2n +1). Equations (6.3), (6.8), (6.18), (6.20), (6.25) and (6.27) compose a system of 8 ordinary differential equations that describe an inverse kinematics model of the HTC motion simulator. Based on these equations, the centrifuge kinematic parameters: ψA, ψ̇A, ψ̈A, φR, φ̇R, φ̈R, θC, θ̇C, and θ̈C are calculated in three phases, according to the algorithm described byKvrgic et al. (2014). 7. Verification of the inverse kinematics model The presented inverse kinematics model of the HTC motion simulator has been tested using numerical calculations. The simulation was performed for the GxC, GyC and GzC acceleration forces profile (Fig. 5, solid line), whichwas generated by the software of theHTCcontrol system. This acceleration forces profile changes as follows: • starting from 1g with G-onset rate n =0.2g/s up to the baseline level (1.41g), • constant baseline level at 1.41g, • increase of the acceleration with G-onset rate n =3g/s up to 6g, • constant acceleration at 6g, • decrease of the acceleration with G-onset rate n =−3g/s up to the baseline level (1.41g). Additionally, for calculations the following data, namely dA = 8m, g = 9.81m/s 2, time step dt =0.005s andMatlab/SimulinkMathWorks software, have been used. Figures 4-7 present the results of numerical simulations (dotted line)plotted togetherwith the correspondingparameters which were recorded during operation of the HTC simulator (solid line). The figures show: • angular velocity ψ̇A and acceleration ψ̈A of the arm (link 1) (Fig. 4), • GxC, GyC, and GzC acceleration forces profile (Fig. 5), • angle φR, angular velocity φ̇R and acceleration φ̈R of the roll ring (link 2) (Fig. 6), • angle θC, angular velocity θ̇C and acceleration θ̈C of the pitch cabin (link 3) (Fig. 7). The calculated angular velocity of the arm ψ̇A (dotted line on the upper plot in the Fig. 4), which is responsible for generating the centripetal acceleration an, largely covers the envelope of the angular velocity obtained from the HTC control system (solid line). The difference between the two curves of ψ̇A (Fig. 4, model vs. HTC) is noticeable only during a decrease in the acceleration from 6g to 1.41g with G-onset rate n = −3g/s (Fig. 5). The maximum of this difference is approximately 0.2rad/s. A similar difference was found in the study by Vidaković et al. (2012). In Fig. 5, the GxH, GyH, and GzH components of the absolute acceleration force obtained by equations (5.3) are given. The curves are very close to each other, except for the phase of An inverse kinematic model of the human training centrifuge motion simulator 109 Fig. 4. Angular velocity and acceleration of the centrifuge arm Fig. 5. Components of the absolute acceleration 110 R. Lewkowicz, G. Kowaleczko deceleration (negative angular acceleration ψ̈A shown in Fig. 4), when a small difference has appeared. The maximum error of the absolute acceleration does not exceed the value of 0.2g (GzH). FromFig. 6, it is clear that thepresentedangleφR, angular velocity φ̇R andacceleration φ̈R of the roll ring (link 2) provide good results.Adifferencebetween thedesired (HTC)and calculated (model) angle φR is minimal. The angle φR is derived from equation (6.17) and depends on Gx0 acceleration (5.1). Thus, the observed difference in the calculated angle φR comes from angular velocity of the arm ψ̇A (Fig. 4) which affects Gx0 acceleration. Moreover, according to equation (6.20), the angular velocity φ̇R, and acceleration φ̈R of the roll ring are calculated on the basis of the angle φR. Therefore, these parameters curves (model vs. HTC) are also different. Fig. 6. Kinematic parameters of the roll ring: angle, angular velocity and acceleration Figure 7 shows the angle θC, angular velocity θ̇C and acceleration θ̈C of the pitch cabin (link 3) obtained by equations (6.25), and (6.27), respectively. A difference between the desired (solid line) and calculated (dotted line) angle θC appears only for the phase of deceleration (negative angular acceleration ψ̈A, shown in Fig. 4) of the arm. Similar to the kinematic para- meters calculated for the ring (link 2), the observed difference in the calculated angle θC (6.25) depends on the centrifuge armmovement (angular acceleration of the arm ψ̈A, which affects Gy0 acceleration (5.1)). The angular velocity θ̇C and acceleration θ̈C of the roll ring (Fig. 7) are calculated based on equation (6.27). Therefore, the differences between two curves (model vs. HTC) of these parameters are easily observed. The obtained kinematic model is not completely accurate, but the calculated link accelera- tions, velocities and angles do not differmuch from their actual values. It is a concern, especially for maximum andminimum values indicating whether the limit ranges (Table 1) have not been exceeded to achieve the desired values of the acceleration vector GC components. In order to eliminate the differences between the desired (HTC) and calculated (model) parameters for the An inverse kinematic model of the human training centrifuge motion simulator 111 Fig. 7. Kinematic parameters of the pitch cabin: angle, angular velocity and acceleration movement having a negative acceleration onset, Kvrgic et al. (2014) proposed a simple solution in which the values of the positive G-onset rate n of the samemagnitude are reversed. 8. Conclusions The purpose of the work is to present a way to solve the problem of correctly defining complex acceleration profiles that are recreated by a centrifugemotion simulator. The proposed solution, in the form of an inverse kinematic model of the centrifuge, indicates not only the exceeded limit values of the parameters, but also their changes over time. The simulation has revealed that the developed inverse kinematic model makes it possible to calculate the angular displace- ment, velocity and acceleration of the links, which is needed for the given linear acceleration of the simulator cabin. Simulation performed in Simulink proved the correctness of the presented expressions for angular displacement, velocity and acceleration of the centrifuge links. The pre- sented algorithm achieved the predefined profile of absolute acceleration in the centrifuge cabin where the onset rate of the absolute acceleration is constant. Thedevelopedmodel of the inverse kinematics can be used for computer simulation of motion of the centrifuge simulator system. Through an overview of the behaviour of the model under various operating conditions, it is possible to predict how a real systemwill behave. References 1. AMST-Systemtechnik GmbH, 2011, User manual Human Training Centrifuge HTC-07 2. Baillieul J., 1985,Kinematic programming alternatives for redundantmanipulators,Proceedings of 1985 IEEE International Conference onRobotics andAutomation, St. Louis,MO,USA, 722-728 112 R. Lewkowicz, G. Kowaleczko 3. Buss S.R., Kim J.S., 2005, Selectively damped least squares for inverse kinematics, Journal of Graphics Tools, 10, 37-49 4. ChenY.C.,ReppergerD.W., 1996,A study of the kinematics, dynamics and control algorithms for a centrifuge motion simulator,Mechatronics, 6, 829-852 5. Crosbie R.J., 1988, Dynamic Flight Simulator Control System, United States Patent, Number 4,751,662 6. Dančuo Z., Rasuo B., Bengin A., Zeljkovic V., 2018, Flight to Mars: Envelope simulation in a ground based high-performance human centrifuge, FME Transactions, 46, 1-9 7. Dančuo Z., RasuoB., Vidaković J., KvrgićV., BućanM., 2013,Onmechanics of a high-G human centrifuge,Proceedings in Applied Mathematics and Mechanics, 39-40 8. Dančuo Z., Vidaković J., Kvrgic V., Ferenc G., Lutovac M., 2012a, Modeling a human centrifuge as three-DoF robotmanipulator, 2012 Mediterranean Conference on Embedded Compu- ting, MECO, IEEE, 149-152 ISBN 9940943601 9. Dančuo Z., Zeljković V., Rašuo B., Dapić M., 2012b, High-G training profiles in a high performance human centrifuge, Scientific and Technical Review, 62, 64-69 10. Djuric A., Al Saidi R., Elmaraghy W., 2012, Dynamics solution of n-DOF globalmachinery model,Robotics and Computer-Integrated Manufacturing, 28, 621-630 11. GhermanB., PislaD., VaidaC., PliteaN., 2012,Development of inverse dynamicmodel for a surgical hybrid parallel robot with equivalent lumpedmasses,Robotics and Computer-Integrated Manufacturing, 28, 402-415 12. Grotjahn M., Heimann B., Kuehn J., Grendel H., 2004, Dynamics of robots with parallel kinematic structure,The 11th World Congress in Mechanism and Machine Science, 1689-1693 13. Kvrgic V.M., Vidaković J., Lutovac M.M., Ferenc G.Z., Cvijanovic V.B., 2014,A con- trol algorithm for a centrifugemotion simulator,Robotics andComputer-IntegratedManufacturing, 30, 399-412 14. Liwen G., Hui L.I.U., Meng F.U., 2015, Real-time motion planning algorithm for dynamic flight simulators (in Chinese),Tsinghua Science and Technology, 55, 709-715 15. Nakamura Y., Hanafusa H., 1986, Inverse kinematic solutions with singularity robustness for robot manipulator control, Journal of Dynamic Systems Measurement and Control, 108, 163-171 16. Nearchou A.C., 1998, Solving the inverse kinematics problem of redundant robots operating in complex environments via a modified genetic algorithm, Mechanism and Machine Theory, 33, 273-292 17. Newman D.G., 2015,High G Flight: Physiological Effects and Countermeasures, 1st ed., Ashgate Publishing Ltd., Monash University, Australia, ISBN 9781472414571 18. Siciliano B., Sciavicco L., Villani L., Oriolo G., 2009,Robotics: Modelling, Planning and Control, Soft Computing, Springer, ISBN 9781846286414 19. TejomurtulaS.,KakS., 1999, Inversekinematics in roboticsusingneuralnetworks, Information Sciences (Ny), 116, 147-164 20. Truszczyński O., Kowalczuk K., 2012, The Polish centrifuge as a dynamic flight simulator. New application and ideas,Polish Journal of Aviation Medicine and Psychology, 18, 71-80 21. Tsai L., 1999, Robot Analysis: the Mechanics of Serial and Parallel Manipulators, 1st ed., John Wiley & Sons, Inc., NewYork, NY, USA, ISBN 978-0-471-32593-2 22. Vidaković J., FerencG., LutovacM.,KvrgicV., 2012,Development and implementation of an algorithm for calculating angular velocity of main arm of human centrifuge, 15th International Power Electronics and Motion Control Conference (EPE/PEMC), Novi Sad, Serbia, DS2a.17-1- -DS2a.17-6 An inverse kinematic model of the human training centrifuge motion simulator 113 23. Vidaković J., LazarevicM., Kvrgic V.M., Dančuo Z., LutovacM.M., 2013, Comparison ofnumerical simulationmodels for open loopflight simulations in thehumancentrifuge,Proceedings in Applied Mathematics and Mechanics, 485-486 24. Wampler C.W., 1986,Manipulator inverse kinematic solutions based on vector formulations and damped least-squaresmethods, IEEE Transactions on Systems, Man, and Cybernetics, 16, 93-101 25. Wojtkowiak M., 1991, Human centrifuge training of men with lowered +Gz acceleration tole- rance,Physiologist, 34, 80-82 26. WolovichW.A., ElliottH., 1984,A computational technique for inverse kinematics,The 23rd IEEE Conference on Decision and Control, Las Vegas, Nevada, USA, 1359-1363 27. Wrigge S., 1981, Calculation of the Taylor series expansion coefficients of the Jacobian elliptic function sn(x,k), Mathematics of Computation, 36, 555-564 28. Wu J., Chen X., Li T., Wang L., 2013, Optimal design of a 2-DOF parallel manipulator with actuation redundancy considering kinematics and natural frequency,Robotics and Computer- -Integrated Manufacturing, 29, 80-85 29. Wu J., Wang J., Wang L., Li T., 2009, Dynamics and control of a planar 3-DOF parallel manipulator with actuation redundancy,Mechanism and Machine Theory, 44, 835-849 30. Wu J., Wang J., You Z., 2010, An overview of dynamic parameter identification of robots, Robotics and Computer-Integrated Manufacturing, 26, 414-419 31. ZhaoY., GaoF., 2009, Inverse dynamics of the 6-DOFout-parallelmanipulator bymeans of the principle of virtual work,Robotica, 27, 259-268 Manuscript received January 14, 2018; accepted for print July 26, 2018