Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

Maximilian Gießler; Bernd Waltersberger

doi:10.1007/s10514-023-10092-x

Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

Gießler, Maximilian; Waltersberger, Bernd 2023-04-01 00:00:00 Maintaining stability while walking on arbitrary surfaces or dealing with external perturbations is of great interest in humanoid robotics research. Increasing the system’s autonomous robustness to a variety of postural threats during locomotion is the key despite the need to evaluate noisy sensor signals. The equations of motion are the foundation of all published approaches. In contrast, we propose a more adequate evaluation of the equations of motion with respect to an arbitrary moving reference point in a non-inertial reference frame. Conceptual advantages are, e.g., getting independent of global position and velocity vectors estimated by sensor fusions or calculating the imaginary zero-moment point walking on different inclined ground surfaces. Further, we improve the calculation results by reducing noise-amplifying methods in our algorithm and using speciﬁc characteristics of physical robots. We use simulation results to compare our algorithm with established approaches and test it with experimental robot data. Keywords Multi-body dynamics · Equation of motion · Moving reference point · Imaginary zero-moment point · Humanoid robot 1 Introduction examples are household chores, care of the elderly, and use as a service worker. Humanoid robots have a signiﬁcant advantage over all other The prerequisite of control algorithms for maintaining types of robots. Since they are designed to resemble the stability is the calculation of the robot’s dynamics. Euler’s human body, they are suitable for interacting in a human- theorem is often used in robotics as a starting point for made habitat. Consequently, the research ﬁeld of humanoid dynamics calculation, e.g. Henze et al. (2016), Buschmann robotics is currently gaining importance. A crucial aspect of (2011), Hirukawa et al. (2006) and Kajita et al. (2014)use increasing the autonomy of humanoid robots is to improve it. This theorem states that the absolute time derivative of dynamic stability during locomotion in challenging settings, the total angular momentum with respect to a ﬁxed reference e.g., the robot’s ability to effectively tolerate external pertur- point 0 equals the total resultant moment about the same (0) (0) bations during dynamic, omnidirectional locomotion. On the ˙ reference point (L = M ). other hand, locomotion in arbitrary environments is of great The approaches from Henze et al. (2016) and Buschmann relevance. In this context, the main focus is to enhance the (2011) require methods to determine at least the absolute robot’s ability to cope with uneven and unknown ground sur- position vector of one body segment. If this vector is not faces of the man-made habitat. Given that stability-enhancing explicitly measurable by sensors, it is estimated by sensor control algorithms are implemented, the robot’s autonomy fusion. E.g., Henze et al. (2016) use the kinematics, absolute will increase. This opens up new areas of application. Some position vector and the orientation provided by IMUs of all end effectors in contact. Buschmann (2011) uses the odome- Maximilian Gießler try. In general, sensor fusions or odometry can be sensitive to maximilian.giessler@hs-offenburg.de errors, because of signal noise or undeﬁned contact impacts Bernd Waltersberger the accuracy. Further, the proposed sensor fusion from Henze bernd.waltersberger@hs-offenburg.de et al. (2016) requires the assumption of rigid contact of at Department of Mechanical and Process Engineering, least one body segment. Likewise, a stand-alone calculation University of Applied Sciences Offenburg, 77652 Offenburg, of absolute position vectors from the kinematic chains of the Germany 123 Autonomous Robots robot is not feasible. Kajita et al. (2014) calculate the angular the calculation of the remaining kinematic quantities for the (0) momentum L and propose to use a numerical differentia- following kinematic chain using the relative derivative and (0) tion to determine its time derivative. To calculate L both the the kinematic state of all joints. absolute position and linear velocity vector are required. Due Moreover, we use simulation and experimental results of to the mentioned aspects, we regard the requirement of the the humanoid robot Sweaty to show our approach’s capa- absolute position or linear velocity vector as the main prob- bilities applied on the physical robots. Section 5 illustrates lem of these published approaches. This especially applies the outcomes. In the article’s ﬁnal parts, we discuss the dif- to the need to handle undeﬁned contacts and noisy sensor ferences and impacts on physical robots of our approach signals on real robots. compared to Henze et al. (2016), Buschmann (2011) and Also, a human being has no sense organs to directly per- Kajita et al. (2014) and give an outlook to future work. ceive its global position or absolute velocity with respect to a ﬁxed spatial point. Further, it does not know the position vec- tor of its center of mass (COM). In conclusion, human control 2 Previous work algorithms for stabilizing locomotion can operate without this information. Several approaches to describe the stability of humanoid Considering all mentioned aspects, we propose an robots during locomotion have been published. Mostly, their approach that does not require any information about the aim was to determine an indicator assessing stability in global position, velocity, or COM vector of the robot. We dynamic situations quantitatively. One of the most discussed achieve this by evaluating the Newton–Euler equations of and applied stability criteria uses the notion of the zero- motion (EOM) with respect to an arbitrary moving reference moment point (ZMP). If the ZMP lies within the support point in a non-inertial reference frame (NIF). The dynamics region, the current dynamic equilibrium of the robot can be calculation is then independent of absolute position vectors assumed stable in the sense that the robot cannot tip over of the robot. Similar approaches using a moving reference an edge of the support region. The ZMP is the spatial point point are common in multi-body dynamics, as seen e.g. in in the ground surface plane where the components of the Wittenburg (1977). moment vector resulting from the reaction force parallel to Furthermore, one goal is to provide an approach that is this plane are zero. Vukobratovic et al. (1970) and Vuko- usable on real robots without ﬁltering the sensor signal. bratovic´ and Stepanenko (1972) introduced this concept to Therefore, we propose to measure the angular acceleration the ﬁeld of robotics. Their approach can be regarded as a vector via angular velocity and linear acceleration measure- generalization of the classical static tilting stability concept ments to reduce the number of noise-amplifying numerical for a rigid body, e.g. Mehrtens (1887), where tipping over differentiation. So this work contributes to the practical use edges of the support region due to external forces is prevented of evaluating the EOMs on physical robots with noisy sensor when the resultant of the externally applied forces and the signals. body’s weight applied to the effective point of action inter- For this work, we choose the imaginary zero-moment sects the support region. Algorithms for calculating the ZMP point (IZMP) to evaluate the results of the EOM calcula- are shown by, e.g., Kajita et al. (2014). The work of Sardain tion. We ﬁnd the IZMP gives additional valuable insight and Bessonnet (2004) shows that the ZMP coincides with the into the stability evaluation of gait behaviors compared to COP under the assumption that the soles are in ﬂat contact the sole evaluation of the position of the center of pressure with one single surface and the robot is free of external per- (COP) measured with force/torque sensors. Especially for turbation forces and moments. If the theoretically calculated gait behaviors with possible external perturbations, calculat- ZMP differs from the COP then this is obviously caused by ing the distance of the IZMP from the support region of the perturbation forces or moments, since such perturbations are robot’s feet or as well from the COP could give additional not considered in the kinematical ZMP calculation. Vukobra- quantitative information as to the impact of the perturbation tovic et al. (2001) hence proposed the name IZMP for this or control strategies. conﬁguration of the ZMP, to which we adhere in our work. In the following Sect. 2, we review the state of the art of Several other concepts exist to quantify stability using referring control algorithms and stability criteria and point a point in space. For example, there is the foot-rotation- out where we see further contributions achieved through our indicator by Goswami (1999a, 1999b), the feasible-solution- approach. In Sect. 3, we explain the mathematical founda- of-wrench by Saida et al. (2003), the zero-rate-of-change tions. In Sect. 4, we discuss the methods used for providing angular momentum by Goswami and Kallem (2004)orthe the data and sensor signals required for the core algorithm capture point (CP) by Pratt et al. (2006). Furthermore, meth- (seen in Sect. 3). On the one hand, we describe a way to ods were developed to modulate contact state and control measure the total kinematic state, especially the angular the ground reaction forces. Hirukawa et al. (2006), Anitescu acceleration of the root segment. On the other hand, we show and Potra (1997), Pang and Trinkle (2000) propose numer- 123 Autonomous Robots ical methods to search for feasible solutions and, based on touchdown introduces a new reference frame. It is placed their characteristics, predictions are made about the stability below the stance foot. This avoids numerical inaccuracies proposed. through large distances from the initial reference point. All the mentioned proposals use the EOMs as a founda- Finally, Buschmann (2011) also uses a dynamics calcula- tion. The level of detail of the mechanical model is essential tion based on a ﬁxed reference point in the current reference for the accuracy. Different approaches have been established frame. over the last decades. On the one hand, there are concep- In contrast, we propose an algorithm to evaluate the EOMs tual models such as the linear inverted pendulum (LIP), e.g., using Newton–Euler equations for a rigid multi-body system shown in Kajita et al. (2010). The use case is to implement related to an arbitrary moving reference point in a NIF. To simpliﬁed control algorithms in combination with stability maintain generality, we use the ﬂoating base model. We use criteria. For example, Englsberger et al. (2011) present a the rigid body theory typically known from the ﬁeld of multi- bipedal walking control approach based on the LIP and the body dynamics to derive the EOMs. We further use the so- CP dynamics. Mesesan et al. (2021) introduced a trajectory called relative derivative to evaluate the EOMs referenced in adaptation for push and stumble recovery based on the CP. a NIF as well as to calculate all kinematic quantities for the Kamioka et al. (2018) are optimizing the ZMP and footsteps multi-body segments. based on the CP for a LIP. Also, Wang et al. (2018) and InJoon In combination, this will allow us to carry the reference et al. (2020) use the concepts of LIP and ZMP to stabilize the point and the NIF along with arbitrary locomotion of the gait of humanoid robots. On the other hand, modeling robots robot. E.g., the reference point can be an arbitrary spatial as rigid multi-body models has become established, see, e.g., point or ﬁxed to any material or immaterial point of the robot. in Vijaykumar Shah et al. (2013). Featherstone (2008)shows Thus, for different types of control algorithms, an adapted the ﬂoating-base model as a kind of rigid multi-body model reference point can be used. Also, we will show that our and it is deﬁned most generally. algorithm requires no absolute position and linear velocity A further distinctive criterion is the calculation method of vector of the root segment. Thus, our approach omits the mea- the time derivatives of the EOMs. Either numerical (Kajita surement or estimation by sensor fusion to determine those. et al., 2014) or analytical derivatives of the equations are To evaluate the effects when applying these approaches proposed. For the numerical derivatives of the linear and on physical robots, we will simulate the inﬂuence of realistic angular momentum, the knowledge of the linear velocity sensor noise behavior. We will also show the impact of feasi- of at least one segment is required. Wang et al. (2017)or ble sensor and actuator controller frequencies. In addition, we Kimura et al. (2018) show approaches to estimate the linear will demonstrate a way to measure the angular acceleration velocity vector. As mentioned before, mostly the reference of the root segment using gyroscope and linear acceleration point for the balance of angular momentum is referenced to data. We will compare this method against the numerical a ﬁxed point in the inertial frame (WF) or the immaterial calculation from the angular velocity signal. COM. Kajita et al. (2014) propose to calculate the rate of Furthermore, using a moving reference point as the origin change of angular momentum by applying a numerical time of a moving evaluation frame allows following the contour derivative. Hirukawa et al. (2006) use the rate of change of of the ground’s surface. IMUs located in the feet provide angular momentum without specifying an explicit equation. information about the inclination of the surface during the Only the equation of angular momentum is stated explicitly. stance phase. A transition for changing surface inclinations Further, Kajita et al. (2003) also propose a method to control is no longer necessary for calculating, e.g., the IZMP. the total linear and angular momentum. Besides combining the mentioned approaches with our The proposed analytical derivative methods require also work, we intend to fully outline the computation procedure the absolute position vector of the root segment from a spa- of the dynamics and the IZMP of a robot omitting noise- tially ﬁxed reference point as an input. Henze et al. (2016) amplifying numerical differentiation for critical terms in the use a dynamics calculation based on a ﬁxed reference point in procedure and without estimations of global vectors. For an a WF. They estimate the absolute position vector of the root easy application, the backgrounds and the algorithm we are segment by sensor fusion. This estimation requires infor- going to present are fully adapted to use the Kinematics and mation about the kinematics, absolute position vectors and Dynamics Library (KDL) and to integrate into the Robot orientations of all end effectors in contact. In their approach, Operating System 2 (ROS2). they calculate the absolute position vector of the root seg- ment from each end effector in rigid contact. The redundant results are combined in a sensor fusion to estimate global 3 Backgrounds position vectors. Buschmann (2011) describes a method to introduce a At ﬁrst, we introduce in Fig. 1 the frames and vectors. NIF is new reference frame for every gait cycle. Each swing leg the frame, in which the EOMs are evaluated. Its origin rep- 123 Autonomous Robots the moving reference frame. The term ω × p describes the i0 rate of change of the vector due to the rotation of frame i. For the implemented EOMs, Newton’s second law and angular momentum theorem of Euler, are taken as starting point. Hence P = F (5) (0) (0) L = M . (6) (0) apply. Here, P = mv is the linear and L the angular COM momentum, which we discuss in more detail in the follow- (0) ing section. F is the total external force and M the total moment with respect to the origin of WF applied on the body. Fig. 1 Deﬁnition of the vectors and frames used to describe the kine- Due to the linearity of the mathematical operators used, the matic description of a rigid multi-body system balance of momentum and angular momentum for the whole robot can be calculated as the sum of the single rigid-body resents the arbitrary moving reference point for the angular segments. momentum balance. Each rigid-body part is given a body- In the following, the velocity v is the ﬁrst and the accel- ﬁxed frame (BFF). The origin of a BFF is located at the eration a is the second total time derivative of r .Inthis X X center of the joint axis considering the mechanical joint struc- context, X represents all the following different subscripts ture. The equivalent transformations used in the following are used, e.g. O or COM . Substituting the second total time i i based on the rigid body theory used in multi-body dynamics. derivative of Eq. 3 into Eq. 5 and by taking into account the The vector from the origin of the NIF to the corresponding law of the COM, results in BFF of an arbitrary segment i is called r . The vector from an arbitrarily positioned inertial frame to the origin of a BFF F = m a + ω ˙ × s + ω × (ω × s ) . (7) i i O i0 i i0 i0 i is denoted as r . The vector from the origin of the BFF to any mass element dm is declared as r . The vector s points from k i (0) For L , we start from the following equation the origin of the BFF to the segment COM and x points from the COM to an arbitrary mass element dm. The deﬁnitions (0) L = m r × v + ω × s i i O O i0 i i i of these vectors are given by (8) (O ) +s × v + ω , i O i0 r = r + r (1) O k proposed by Wittenburg (1977), who therefrom derived a r = r + r (2) O NIF H i i formulation with respect to an arbitrary moving reference r = r + s . (3) COM O i i i point. However, our approach differs so far that we will derive a formulation that assumes an arbitrary moving reference Here, r is the position vector pointing from the origin NIF point for the moments while using another arbitrary reference of a arbitrary WF to the origin of NIF. The absolute time point, e.g., pivot points for the inertia tensors. The latter is derivative of an arbitrary vector p = p e is calculated i β i β β a common practice in the generic description of input data, according to the basic equation of kinematics as follows e.g., in C++ libraries. Form now, i indicates the reference of the quantities to the dp (O ) dp de i β i β i ith body of a multi-body system. In Eq. 8, is the inertia = e + p i β i β dt dt dt tensor of the ith body with respect to the reference point O . β=x ,y,z β=x ,y,z (4) We substitute Eq. 2 and the total time derivative of Eq. 3 into dp = + ω × p. i0 Eq. 8. Further, we take the total time derivative of Eq. 8.This dt gives Here, p is the coordinate in direction β of p and e is the i β i β (0) unit vector of direction β.The sufﬁx i indicates an arbitrary L = m r + r × a i i NIF H COM i i frame in which the vector is referenced and ω indicates the i0 +s × a i O dp (9) frame’s angular velocity. The term represents the so called dt (O ) (O ) i i + ω ˙ + ω × ω i0 i0 i0 relative derivative. This term can be interpreted as the vector’s (0) rate of change that one would observe if one were ﬁxed to = M . 123 Autonomous Robots We see all absolute velocity terms cancel out. forces and moments, e.g., caused by the gravity force or the Next, we subtract the term r × m a from Eq. 9. ground reaction force remain. This results in NIF i COM Further, for the moment on body i and by using Eq. 5 (P) (0) F = (m (a + ω ˙ × s ext i O i0 i M = M − r × F NIF i i i i =1 (10) (0) = M − r × m a NIF i COM i + ω × (ω × s ))) (12) i0 i0 i (P) (P) applies in general. Here, M describes the total moment M = m (r × a + r × (ω ˙ × s ) i H O H i0 i ext i i i i =1 applied on body i with respect to the arbitrary moving refer- ence point P. We get ﬁnally + r × (ω × (ω × s )) + s × a ) H i0 i0 i i O i i (O ) (O ) i i (P) + ω ˙ + ω × ω . (13) i0 i0 i0 M = m (r × a + r × (ω ˙ × s ) i i i H O H i0 i i i i + r × (ω × (ω × s )) + s × a ) (11) H i0 i0 i i O i i The algorithm is recursively designed and thus generic. The (O ) (O ) i i + ω ˙ + ω × ω . i0 i0 i0 mechanical structure and the physical properties of each rigid body segment of the robot can be provided via a Uniﬁed Robot Equation 11 describes the change in angular momentum Description Format ﬁle. for an arbitrary moving reference point. From Eq. 11,we For the implementation in the code, we switch from vec- see that all absolute position and velocity vectors cancel out. tor to coordinate matrix notation. E.g., the representation of The Eqs. 7 and 11 describe the Newton–Euler EOMs for an Eq. 13 is given by arbitrary moving reference point in a NIF. (P) NIF NIF i i i M = m R ( r × a i H O ext i i 4 Methods for a robust and eﬃcient i =1 evaluation of the inverse dynamics i i i + r × ω ˙ × s H i0 i (14) i i i i With our approach, we intend to provide an effective cal- + r × ( ω × ( ω × s )) H i0 i0 i culation of the inverse dynamics. We use this to calculate i i i (O)i + s × a ) + ω ˙ i O i i i0 the IZMP exemplarily. We believe the IZMP provides an i i (O) i additional beneﬁt in evaluating gait situations under exter- + ω × ω . i0 i i0 nal perturbations compared to a COP measured using 6-axis force/torque sensors. For example, the distance of the IZMP NIF i Here, R is a 3 × 3 rotation matrix that transforms the from the tipping edge or the measured COP provides a metric for the needed response to stabilize the robot in the presence coordinates with respect to the frame i to the coordinates with respect to the frame NIF. If a symbol is superscripted of external perturbations. We set the goal that our approach is efﬁcient and works stand-alone with noisy and unﬁltered on the left side by i, this means that the referring vector or tensor is expressed as its 3 × 1 (vector) or 3 × 3 (tensor) inertia sensor signals. Other inaccuracies or disturbances in the evaluation of the EOMs are, e.g., measurement errors of coordinate matrix with respect to frame i. This notation is adopted from Woernle (2016). Each summand of the right- encoder inadequacies in the mechanical design of the robot, the underlying rigid body hypothesis, or modeling errors of hand side of Eq. 14 is a 3 × 1 coordinate matrix. The number of n rigid bodies corresponds to the number of summands. the inertia data. However, these types of disturbance are not part of this work. By sensor fusion with 6-axis force/torque sensors and combined with an adapted Kalman ﬁlter, our 4.2 Analytically derived EOMs contribution is extendable to become less sensitive to the mentioned inaccuracies in the model assumption. We link the arbitrary moving reference point and the origin of the NIF to a material point of the robot. This point coincides 4.1 EOMs for a multi-body system with the origin of the root segment, called pelvis. As seen in Fig. 2,the x-axis of the NIF is deﬁned as pointing always To describe the dynamics of a multi-body system, we sum parallel to the sagittal plane of the robot and the z-axis points up the Eqs. 7 and 11 for all n rigid bodies. In this sum, the parallel to the normal vector of the ground surface. Also, it reaction forces and moments between two body segments i is feasible to place the NIF on the ground surface. E.g., the and i + 1 cancel out, due to Newton’s third law. The external origin could coincide with the moving immaterial point of 123 Autonomous Robots Due to the addition formula of angular velocities, for Eq. 15 follows dω dω dω i i0 i (i −1)0 i i (i −1) = + dt dt dt (16) dω i (i −1)0 := +¨ q n. dt Here, q ¨ n is the joint acceleration vector. The vector n describes the unit vector of the joint axis. In general, ω ˙ (i −1)0 can be calculated by dω i (i −1)0 ω ˙ = + ω × ω , (17) (i −1)0 i0 (i −1)0 dt using Eq. 4. Equation 17 can be rearranged accordingly dω i (i −1)0 = ω ˙ − ω × ω (i −1)0 i0 (i −1)0 dt = ω ˙ (i −1)0 (18) Fig. 2 Visualization of the vectors and frames WF and NIF in the − (ω + ω ) × ω (i −1)0 i (i −1) (i −1)0 mentioned application related to the robot Sweaty. The vector r pelvis describes the absolute position of the root segment with respect to the = ω ˙ −˙ q n × ω . (i −1)0 i (i −1)0 origin of the WF. r and r indicate the relative position vectors foot foot R L from the origin of the NIF to origin of each foot Here, q ˙ is the joint velocity of body i compared to i − 1. For ω ˙ ﬁnally follows i0 ω ˙ = ω ˙ + ω ×˙ q n +¨ q n. (19) the vertical projection of the origin of the BFF from the root i0 (i −1)0 (i −1)0 i i segment onto the ground surface. (0) i (0) ˙ The matrix ω ˙ is given by, i0 On the other hand, the approaches based on L = M evaluate the EOMs with respect to a ﬁxed spatial point. Fig- i i i i i ω ˙ = ω ˙ + ω ×˙ q n +¨ q n. (20) ure 2 displays the difference between the mentioned (Henze i0 (i −1)0 (i −1)0 i i et al., 2016; Buschmann, 2011; Kajita et al., 2014) and our 4.4 IZMP calculation approach. These approaches use an absolute position vec- tor pointing from the ﬁxed reference point to at least one To calculate the IZMP, the EOMs will be evaluated with mea- material point of the robot, e.g., r . With kinematics, it pelvis is then feasible to calculate the absolute position vectors of sured sensor and actuator information at discrete points in all segments with respect to the ﬁxed reference point. Since time. Due to the time derivatives being available as measured values, we can interpret the coupled nonlinear differential our chosen reference point coincides with one material point, no absolute position vector is needed. Thus we can use rela- equations as algebraic equations. We formulate all equations presented in Sect. 3 as continuous-time equations. They can tive vectors, e.g., r and r pointing from the reference foot foot R L point to the origin of each foot. be evaluated without adjustments at any discrete instant in time t . For the calculation of the IZMP, we use Eq. 13. We assume the ground reaction force F and the gravity force mg act gr 4.3 Relative derivative of angular velocity externally on the robot. Then the following applies to the sum (P) of all applied external moments M , ext We use the relative derivative to link the joint kinematics with the total angular velocity ω and angular acceleration i0 (P) (IZMP) M = M + r × F + r × mg (21) IZMP gr COM ext ω ˙ . Using the law of time derivative (Eq. 4), the following i0 relationship generally applies (IZMP) Here, M is the resultant vector of the moment at point r . Per the deﬁnition of the IZMP, the components of IZMP (IZMP) M are zero for the horizontal plane. When we substi- dω dω i i0 i i0 ω ˙ = + ω × ω = . (15) i0 i0 i0 tute Eq. 5 into 21, the ground reaction force can be replaced. dt dt 123 Autonomous Robots Further, substituting this into Eq. 13, the IZMP calculation equation is given by (IZMP) M + r × (P − mg) + r × mg IZMP COM = m (r × a + r × (ω ˙ × s ) i H O H i0 i i i i i =1 (22) + r × (ω × (ω × s )) + s × a ) H i0 i0 i i O i i (O ) (O ) i i Fig. 3 Location of the IMU and the accelerometers on a rigid body to + ω ˙ + ω × ω . i0 i0 i0 i i calculate the angular acceleration. r describes the position vector of the IMU. r are feasible relative position vectors of the accelerometers with respect to the IMU. r are linearly independent to each other This system of the coordinate equations from Eq. 22 is under- determined due to the skew-symmetric property of the cross product. To obtain a single unique solution of the equation behavior of the angular velocity signal and on the increment system, we calculate the IZMP coordinates upon the assump- t used. tion that the IZMP is located in the ground surface plane. For our approach, we propose an algebraic method to cal- Therefore, the z-coordinate of the IZMP is set equal to that culate the angular acceleration vector, similar to the proposed of the reference point in relation to the ground surface plane. methods, e.g., from Padgaonkar et al. (1975). Therefore, we place three additional accelerometers on the root segment and 4.4.1 Required inputs use the rigid body assumption. Figure 3 shows the schematic setup. With this setup, we measure the linear acceleration By analyzing Eq. 13 and regarding the underlying ﬂoating- vector of speciﬁc spatial points and the angular velocity vec- base model, the required input signals can be identiﬁed. From tor of the rigid body. Eq. 13 it follows the evaluation of the dynamics requires the In general, the locations of the sensors on the rigid body position, angular velocity, linear and angular acceleration must satisfy the condition that the connection vectors from vector. For all segments of the multi-body system except the the IMU to the accelerometers are linearly independent. This root segment, these vectors can be calculated through sensor prevents the sensors from lying in a spatial plane. If the sen- measurements of the actuators and the total angular velocities sors are placed with orthogonal connection vectors to each ω and ω ˙ of the previous segment known from the other, then uniaxial accelerometers can be used. (i −1)0 (i −1)0 recursive loop. As shown, these vectors are functions of the For the acceleration of an arbitrary ﬁxed point on a rigid joint angle q , joint velocity q ˙ , and joint acceleration q ¨ .The i i i body applies in general, commonly used actuator encoder in robotics delivers signals for q . We calculate q ˙ and q ¨ with numerical derivatives using i i i a = a + ω ˙ × r + ω × (ω × r ). (23) i 0 i i the actuator controller. Therefore, we can use the increment t corresponding to the frequency of the actuator controller. Eq. 23 can be rearranged as follows, Common frequencies of actuator controllers used in robotics are in the range of kilohertz. Therefore, we name this pro- ω ˙ × r = a − a − ω × (ω × r ). (24) i i 0 i cedure actuator coordinate differentiation (ACD). Due to the option of choosing an arbitrary moving reference point in To separate the term ω ˙ , the cyclic permutation of the vector a non-inertial frame and the analytical way of deriving the triple product is used. For this purpose, the scalar product of EOMs, we name the proposed method MR-AD. Eq. 24 and the connection vector r is formed, as follows Considering the root segment, these vectors must be pro- vided from sensor data or estimations. Since our chosen (ω × r ) · r = (a − a − ω × (ω × r )) · r . (25) i j i 0 i j reference point coincides with a material point of the root segment, we can use a relative body-ﬁxed position vector. Applying the cyclic permutation Eq. 25 gives The angular velocity is measured by a gyroscope and the linear acceleration vector by an accelerometer. (r × r ) · ω ˙ = (a − a − ω × (ω × r )) · r . (26) In perspective of the authors, direct angular accelera- i j i 0 i j tion measurements are not commonly used in robotics. One approach is to calculate the angular acceleration with a Eq. 26 leads to three scalar equations. Each equation contains numerical derivative of ω . In general, numerical derivatives one coefﬁcient from the angular acceleration vector ω ˙.For are noise-amplifying. So the accuracy depends on the noise i = 2, j = 3; i = 3, j = 1 and i = 1, j = 2 the following 123 Autonomous Robots applies, ω ˙ = ((a − a − ω × (ω × r )) · r ) (27) x 2 0 2 3 r r 2 3 a a 2z 0z = − − ω ω y z r r 2 2 ω ˙ = ((a − a − ω × (ω × r )) · r ) (28) y 3 0 3 1 r r 3 1 a a 3x 0x = − − ω ω z x r r 3 3 ω ˙ = ((a − a − ω × (ω × r )) · r ) (29) z 1 0 1 2 r r 1 2 a a 1y 0y = − − ω ω . x y r r 1 1 4.4.2 Momentum based dynamics calculation To put the generated results in perspective, the momentum- Fig. 4 The real humanoid robot Sweaty (left) and the simulation model based differentiation method (IR-MD) for calculating the (right), which was designed to interact autonomously with the human- linear and angular momentum balance proposed by Kajita et made habitat al. is also implemented. The implementation strictly follows the procedure given in Kajita et al. (2014). Due to the numeri- cal derivative calculation, only the equations for momentum CAN-BUS is clocked with 125 Hz. Sweaty was designed to and angular momentum are used. To determine the rate of interact autonomously with the human-made habitat. Each change of these quantities in time, a numerical time derivative sensor used has an equivalent to the human senses. Sweaty has based on the differential quotient is used. Since the calculated no sensors like GPS or linear velocity measurement installed. momentum and angular momentum vectors are subtracted The global position and linear velocity must be estimated. from each other in two different instants in time, they must Figure 4 (right) shows the simulation model of Sweaty. be referenced in the time-constant WF. The purpose of IR- This model has 31 DOF. The DOFs of both hands are not MD is to reduce the complexity of the required equations. implemented. Webots 2021b from the company Cyberbotics Also, any dependencies on the acceleration states of the seg- Ltd. is used as the simulation tool. In the simulation, however, ments and actuators are omitted. Kajita et al. use the ﬂoating all DOFs are simulated as rotational motors. This modiﬁ- base model for the mentioned algorithm. Consequently, the cation is made due to the closed kinematic chains of the MD method requires q and q ˙ and for the root segment the i i rod-and-lever systems leading to a decrease in simulation absolute position, the absolute linear and the angular veloc- performance. ity vector must also be speciﬁed. The increment t of the Sweaty’s software architecture is based on ROS2, where numerical differentiation is linked to the frequency of the information is exchanged in a network according to the pub- ﬁeld bus in physical robot systems. lisher and subscriber principle. Figure 5 shows the ROS2 nodes and their structural connections used to generate the 4.5 Physical robot Sweaty and simulation model results for this work. The interface between ROS2 and the simulator is the node Sweaty, shown in Fig. 4 (left), is 1.65 m tall and weighs 27.6 Webots-Controller-Client. Via this interface, the node Motion kg. The humanoid robot has 39 degrees of freedom (DOF). sends the desired motion pattern to the simulator. In the Spindle motors actuate the DOFs of the lower body and torso. simulator, these data are processed and the physics engine Rod-and-lever systems translate the linear motion of the simulates the dynamics for a ﬁxed time step. The simula- spindle nut into rotational motion. The mechanical structure tion results for joint angles, sensor data, and the current time of the rod-and-lever system leads to a variable gear ratio, are made available by the simulator to the ROS2 network which results in variable torque and angular velocity char- via the node Webots-Controller-Client. Webots 2021b does acteristics depending on the joint angle. This design also not provide information about the angular velocity q ˙ of the minimizes the joint backlash due to the high gear ratio of actuators. Therefore, the required input signals, e.g., angular the spindles. The Xsens MTi-20 VRU IMU is mounted on velocities q ˙ or the angular accelerations q ¨ of the actuators i i the hip of Sweaty. Two CAN-BUS lines provide communi- are calculated through a differential quotient with the simu- cation between actuator controllers and the control unit. The lation time step t.The TF-Framework of ROS2 makes it 123 Autonomous Robots Fig. 5 Nodegraph of the setup to generate the simulation results possible to request arbitrary coordinate transformations by Table 1 List of the main abbreviation using the joint angle information and other sensor data. Abbreviation Meaning Especially for this work, a GPS and an IMU are located ACD Actuator controller differentiation in the root segment. The IMU combines a gyroscope and an BFF Body-ﬁxed frame accelerometer. It delivers the signals for angular velocity, lin- IR-MD Inertial reference—momentum derivative ear acceleration, and the orientation of the root segment. The MR-AD Moving reference—analytically derived GPS signal is subscribed by the node Localizer and published NIF Non-inertial frame as the transformation from WF to the root segment in the TF- WF Inertial frame Framework. Furthermore, this node takes a time derivative of the GPS signal by using the differential quotient to provide a Method to calculate the time derivatives of the actuator position signal the linear velocity vector. Both mentioned calculation meth- q obtaining q ˙ and q ¨ i i i Dynamics evaluated by numerical differentiation of momentum based ods are integrated parallel as separate nodes. This ensures equations for a ﬁxed reference point in inertial frame compatibility, as the calculations access the same input sig- Dynamics evaluated by analytically derived EOMs for an arbitrary nals in each time step. Both calculation methods publish their moving reference point in non-inertial frame results, which are subscribed by the node IZMP Calculation. The calculated coordinates of the IZMP are transformed into choose the BFF of the root segment as the common evalua- a common evaluation frame via the TF-Framework and are tion frame, which corresponds to the introduced NIF (see recorded in ROS2 bag ﬁles. Sect. 4.2). A further quantitative evaluation is the relative | f (t )− f (t )| n proc k ref k 1 error . It is deﬁned as = ·100%· . k=0 | f (t )| n ref k 4.6 Evaluation methods Here, f (t ) represents the value at the time point t of the proc k k investigated procedure’s signal. Similarly, f (t ) represents ref k We verify our MR-AD approach, compare it to IR-MD and the reference signal for the comparison at the time point t . show its general scopes using simulation results. Further, To analyze the applicability to physical robots in real- we test our approach on the physical robot Sweaty. Table 1 world applications, we add white noise behavior to the ideal reviews all of the necessary main abbreviations for the fol- sensor signals to simulate typical noisy signals s ˆ(t ) from lowing sections. physical sensors. White noise behavior is added to linear Due to Sweaty having no GPS or sensors for measuring velocity, angular velocity, and linear acceleration signals. It the linear velocity installed, the global position and the linear generally applies s ˆ(t ) = s(t ) + n(t ). Here, s(t ) is the original velocity vector are not available during the execution of the ideal signal and n(t ) is the noise. n(t ) is generated normally experiment. An estimation by our localization and odome- distributed depending on the standard deviation σ and the try does not provide usable results. Therefore, we focus on mean μ of the physical sensor. We use σ = ν · f to calculate MR-AD for the evaluation of the EOMs on the real robot. the standard deviation for the gyroscope and accelerometer For the comparison in simulation, the mentioned approaches as mentioned in the data sheet of Xsens MTi-20 VRU. ν is the perform different test cases and the obtained results will be noise density and f is the sample rate. The mean μ is set to compared quantitatively through the IZMP coordinates. We zero. 123 Autonomous Robots Sweaty does not have a sensor to measure linear veloc- Table 2 Parameters and properties of the line diagrams ity. Therefore, we use the same standard deviation σ for Figure ACD/bus (Hz) Noise the noise term of the linear velocity signal that applies to Fig. 6 125 None the accelerometer. The following equation applies to cal- Fig. 7 62.5, 125, 200, 250 None culate the noise for each coordinate direction, s ˆ (t ) = β k a b Fig. 8 125 IMU ,GPS s (t ) + n (t ) = s (t ) + g (σ, μ). Here, t indicates the β k β k β k k discrete point in time. Fig. 12 1000/125 Real sensors For this work, we add noise for our approach to a for i Fig. 11 125 None in range i = 1, 2, 3, 4 and ω of the root segment. It applies If the column ACD/Bus contains one value, it belongs to both ACD and a =ˆs and ω =ˆs . Here, β indicates the coordi- i β aβ β ωβ ﬁeld bus frequency Noisy accelerometer and gyroscope signals nate axis in rage β = x , y, z. For IR-MD, we add noise Noisy velocity signal from Webots just to the velocity v of the root segment, in the same way, it applies to a . Thus for IR-MD, we expect the inﬂuence Table 3 Parameters and properties of the bar diagrams of noise as underestimated since velocity measurements are usually indirect measurements based on ﬁltered acceleration Setup Method ACD/bus (Hz) Ang. acc. calc and camera position data. A IR-MD 125 Diff. quotient In general, real robot systems offer the possibility to pro- B MR-AD 125 Algebraic vide higher update frequencies of sensor and actuator signals C IR-MD 500/125 Diff. quotient compared to the ﬁeld bus frequency. To analyze the impacts D MR-AD 500/125 Algebraic of this characteristic in calculating the dynamics, the sim- E MR-AD 500/125 Diff. quotient ulation runs with a time step of 0.002 s (500 Hz). However, the dynamics calculation runs with 0.008 s (125 Hz). Thus all If the column ACD/Bus contains only one value, it belongs to both ACD and ﬁeld bus frequency sensor and actuator signals are available in 500 Hz frequency One IMU and three accelerometers are mounted on the root segment for processing. Consequently, the frequency of the simulation to measure ω ˙ set to 500 Hz can be interpreted as the actuator controller fre- b Differential quotient is used quency, and the dynamics calculation frequency set to 125 Hz as the ﬁeld bus frequency of a physical robot. 5.1 Initial conditions for the simulation 5 Results At ﬁrst, we compare the results of MR-AD and IR-MD under ideal conditions. For this purpose, the robot walks in an arbi- In the following sections, we show the results of the EOM trary direction starting from the origin of the WF at a gait evaluation by the IZMP. Each line diagram listed in Table 2 speed of 0.3m/s. Optimal signal qualities from the simula- contains curve developments for the MR-AD and IR-MD tion are provided for each calculation method. Figure 6 shows method. The adjustable parameters are the ACD and ﬁeld exemplary curves of the x- and y-coordinate from the IZMP. bus frequencies and the signal noise that can be switched on The gait phase of the robot can be divided into two sub-step or off. If signal noise is present, it is applied, as explained in phases as shown in Fig. 6. The jump in both curves describes Sect. 4.6. On the one hand, the column ACD/Bus of Table 2 the switch of the support foot and thus the transition from contains the chosen ACD and ﬁeld bus frequencies for the sub-step one into two. From the data, it follows that a sub- simulated or real robot. On the other hand in the column step takes 0.448 s, and the entire step takes 0.896 s. In the Noise speciﬁes the present signal noise for each test case. development of the x-coordinate, the two sub-steps can be Likewise, the adjustable parameters for all bar diagrams identiﬁed via the consecutively similar curve sections. The are the ACD and ﬁeld bus frequencies and the angular accel- sub-step phase begins with the heel strike and ends with the eration vector determination approach. Table 3 contains the toe-off. The x-coordinate starts after the peaks of the heel chosen parameters and properties of each bar diagram. It strike at approx. − 0.08 m and increases in the development displays in the column ACD/Bus the ACD and ﬁeld bus fre- to approx. 0.1 m until the toe-off event. The x-coordinate, quencies used. In the column Ang. acc. calc. it speciﬁes the related to the foot, moves from the rear area to the front area underlying method of angular acceleration vector estimation during the duration of a sub-step. for each setup. Further, the signal noise is turned on for each The y-coordinate starts with negative values. The develop- setup (see Sect. 4.6). ment of the curve shows a relatively constant range between If the column ACD/Field bus of Tables 2 and 3 contain the time frame 0.1 s and 0.25 s. In the further development only one value, this value represents the frequency of both of the sub-step, the value of the y-coordinate moves towards ACD and ﬁeld bus. the zero line. Here, the robot tilts towards the swing leg. 123 Autonomous Robots Fig. 6 Comparison of the 0.2 0.1 calculated IZMP x-(left)and y-coordinate (right) of MR-AD 0.1 and IR-MD under ideal conditions in simulation. The calculations are executed in 0 0 parallel for the same step phase. The gait speed is0.3 m/s −0.1 MR-AD MR-AD −0.1 IR-MD IR-MD −0.2 00.10.20.30.40.50.60.70.8 00.10.20.30.40.50.60.70.8 Time t [s] Time t [s] 5.2 Simulated robot applications The output signals of physical sensors are generally noisy. To simulate the qualitative impact of this noise, noise behavior is added as explained in 4.6. Figure 8 shows the white noise effects of linear velocity, angular velocity, and linear accel- eration signals on the x- and y-coordinates of the IZMP. For Rel. error MR-AD comparison, the coordinate curves calculated by MR-AD and Rel. error IR-MD IR-MD with ideal input signals are also shown in Fig 8.All 45 8 16 displayed curves are generated with a simulation and dynam- Simulation time step Δt [ms] ics calculation frequency of 125 Hz. Figure 8 shows MR-AD Fig. 7 Convergence behavior for the IZMP y-coordinates of MR-AD and IR-MD do not provide similar curves compared to Fig. 6. and IR-MD for decreasing simulation time steps and ideal signals. The The noise behavior of the sensors has a marginal impact on reference curve for the corresponding method is generated with a time the MR-AD method compared to the curves with ideal sensor step of 0.002 s. The gait speed of 0.3 m/s are the same for each case signals. In contrast, noise behavior has a signiﬁcant impact on the IR-MD method. In contrast to Fig. 8,Fig. 9 shows the relative errors of the IZMP y-coordinate compared between MR-AD and IR- MD. Setup B and D use MR-AD and Setup A and C use The development of the y-coordinate curve indicates that the IR-MD. Further, Setup A and B use a frequency of 125 Hz analyzed step starts with the right leg as the support leg. for ACD and ﬁeld bus. In contrast, Setup C and D use sensor In Fig. 7 the convergence behavior of the IR-MD and MR- and actuator controller frequencies of 500 Hz and a modeled AD methods for the y-coordinate of IZMP are compared at ﬁeld bus frequency of 125 Hz. For IR-MD increased sensor different simulation time steps. The tested frequencies are and actuator controller frequencies impact the accuracy neg- 250 Hz, 200 Hz, 125 Hz and 62.5 Hz. All curves were gen- atively. The relative error increases from 65.2 Hz to 112.4%. erated with ideal signals. The reference curves have been For MR-AD the contrary is observed. Here, increased sen- determined with the corresponding method at a frequency sor and actuator controller frequencies increase the accuracy. of 500 Hz. Figure 7 shows a linear convergence behavior for The relative error decreases from 16.1 to 4.4%. both methods with respect to the simulation time step t. Fig. 8 Comparison between the 0.3 0.3 ideal MR-AD results of the IZMP x-and noised MR-AD 0.2 y-coordinate of MR-AD and 0.2 ideal IR-MD IR-MD with noised linear noised IR-MD 0.1 velocity, angular velocity and 0.1 linear acceleration signals at a simulation time step of 0.008 s −0.1 (corresponds to the ﬁeld bus ideal MR-AD −0.2 noised MR-AD frequency of the real robot −0.1 ideal IR-MD Sweaty). The gait speed is −0.3 noised IR-MD −0.2 0.3 m/s 00.20.40.60.800.20.40.60.8 Time t [s] Time t [s] Rel. error y -coord. [%] x -coordinate [m] x -coordinate [m] y -coordinate [m] y -coordinate [m] Autonomous Robots 120 112.42 0.2 0.1 65.22 Level ground Slope 10 16.05 4.44 −0.1 Setup A Setup B Setup C Setup D Fig. 9 Relative error comparison between MR-AD (Setup B and D) −0.2 and IR-MD (Setup A and C). Setup A and B use a frequency of 125 Hz for ACD and ﬁeld bus. Setup C and D use 500 Hz for ACD and 125 Hz −0.10 0.1 for ﬁeld bus. Setup B and D use the method for measuring of angular x -coordinate [m] acceleration vector of the root segment. The relative error is calculated compared to the corresponding method with ideal signals and a fre- Fig. 11 The IZMP x-and y-coordinate calculated by MR-AD for the quency of 500 Hz for ACD and ﬁeld bus xy-plane and referenced in the evaluation frame while walking on a level surface and slope. The gait speed is 0.15 m/s. The slope angle is set to 10 . The support polygon for each foot is represented with the shown 18.9 polygons 16.05 Next, we analyze the impact of walking onto a slope, espe- cially for the IZMP calculation. The robot starts on level ground (slope angle 0 ) and walks up an inclined plane with 4.44 a ﬁxed slope angle of 10 without stopping. For this purpose, the gait speed is 0.15 m/s. Figure 11 shows the IZMP coordinates referenced in the Setup B Setup E Setup D evaluation frame and displayed in the xy-plane of the ground Fig. 10 Comparison between the different calculation methods. Setup surface for one exemplary gait cycle. Here we use the time B, D and E use MR-AD. Setup B uses a frequency of 125 Hz for ACD and t as the parameter for the xy-plot. In Fig. 11 the solid curve ﬁeld bus. Setup D and E use 500 Hz for ACD and 125 Hz for ﬁeld bus. indicates a gait cycle on the level surface, the dotted curve Setup E uses a differential quotient to calculate the angular acceleration vector of the root segment. Setup B and D use the method for measuring indicates a gait cycle walking on a slope (slope angle 10 ). of angular acceleration vector of the root segment. The relative error The polygons in Fig. 11 represent the support polygon of is calculated compared to the corresponding method with ideal signals each foot referenced in the evaluation frame. Both support and a frequency of 500 Hz for ACD and ﬁeld bus polygons are calculated with the kinematics in the respective stance phase. Due to the support polygons being repre- sented in the NIF, both support polygons have the same Figure 10 compares the impact of the estimation of the angular acceleration on the IZMP y-coordinate. To analyze x-coordinates. For the gait cycle on the level surface as well as for the gait cycle on the slope, it can be seen that IZMP just the impact of the method providing the angular acceler- curves are located inside the support polygon. ation vector of the root segment, Setup B, D and E use the MR-AD. Setup E and B use 500 Hz for ACD and 125 Hz for the ﬁeld bus frequency. Setup E uses a differential quo- 5.3 Real robot application tient of the measured angular velocity and Setup B and D use the angular acceleration measurement method to deter- Figure 12 shows the development of the x- and y-coordinates mine this vector of the root segment. The relative error is of the IZMP calculate on the physical robot Sweaty. The actu- calculated compared to the corresponding method with the ator controllers return the actuator positions with a frequency same angular acceleration calculation method, ideal input of 1 kHz and the CAN-BUS frequency is set to 125 Hz. A signals and a frequency of 500 Hz for ACD and ﬁeld bus. differential quotient calculates the angular acceleration of The numerical differentiation leads to a higher relative error the root segment. The robot walks on ﬂat ground with a gait (18.9%) compared to the method using angular acceleration speed of 0.15 m/s. In comparison to the developments of the measurement (4.4%) under equal test case settings. Setup E x- and y-coordinates of MR-AD shown in Fig. 8, these curves has also a higher relative error compared to Setup B using have a smoother transition in the change of the stance foot. the ACD frequency of 125 Hz and the angular acceleration The x-coordinate curve shows a larger slope followed by measurement (16.1%), shown in Fig. 9. a relatively constant proﬁle at about 0.1 m, compared with Rel. error for y -coord. [%] Rel. error for y -coord. [%] y -coordinate [m] Autonomous Robots Fig. 12 IZMP x-and 0.2 y-coordinate of MR-AD 0.1 calculated on the physical robot 0.1 Sweaty. The gait speed is 0.15 m/s −0.1 −0.1 −0.2 −0.2 MR-AD real robot MR-AD real robot 00.20.40.60.800.20.40.60.8 Time t [s] Time t [s] Fig. 8. Overall, the y-coordinate curve shows a similar but differences in calculating time derivatives cause the devia- more inconstant development as seen in Fig. 8. tions in the development of IZMP coordinate curves shown in Fig. 6). As far as we consider realistic noise behavior of physical 6 Discussion sensor signals, e.g., for the linear velocity or the position vec- tor, the results of the MD method are no longer usable (see Before we discuss our main ﬁndings, we start with the inter- Figs. 8, 9). This is caused by the dependency on the quantities pretation of the simulation characteristics. The convergence from different instances in time. Due to this characteristic, analysis (Fig. 7) shows for AD and MD methods linear con- the white noise behavior of real sensors has a direct impact vergence with respect to the simulation time step t.The on the results. By using discrete sensor signal values, the simulator Webots 2021b uses a ﬁxed time step for the physics dynamics calculation inevitably includes at least two differ- engine. This time step is speciﬁed in milliseconds and can ent noise terms. Due to the randomness of the underlying be set to an integer. Therefore, simulations are limited at normal distribution of the noise, the numerical time deriva- t = 0.001 s (1000 Hz). With this time step, the simulation tives could determine physically unrealistic gradients of the delivers no usable results simulating the robot Sweaty due momentum. Figures 8 and 9 exactly show this behavior. to impaired communication between the Webots-Controller- Further, measuring linear velocity is impractical since a Client and the other ROS2 nodes. Regarding this, we set direct measurement method is not feasible. Thus, the analyt- simulation results generated with a time step of 0.002 s as ically obtained derivatives have an obvious advantage as far the reference. However, due to the linear convergence, we as noisy signals are concerned. Consequently, they are signif- consider it appropriate to deﬁne these curves as reference icantly better suited for calculating the dynamics of physical curves. Further, for the theoretical comparisons of all simu- robotic systems. lation results, we omit measurement inaccuracies of encoders Second, for physical robots, it is commonly feasible to and joint backlash as signal noise for q to exclusively analyze operate sensors and actuator controllers with signiﬁcantly the inﬂuence of white noise of inertial sensors and velocity higher sample rates than the ﬁeld bus. A numerical deriva- estimations. and joint acceleration tive can calculate the joint velocity q ˙ Our ﬁrst main ﬁndings are about the characteristics of the q ¨ using the joint position q provided by the encoder. Fig- i i MD method. Generally, this approach reduces the complexity ure 9 shows the impact of higher sample rates. Due to the of the EOMs since only total linear and angular momentum smaller time step, both methods beneﬁt theoretically by the have to be calculated. Consequently, IR-MD gets rid of joint increased accuracy of the calculated q ˙ and q ¨ (see. Figure 7). i i acceleration q ¨ , linear and angular acceleration vectors of Additionally, the measurement error of incremental encoders the root segment. However, it requires the linear velocity does not scale proportionally with the square root of the sam- vector of the root segment. Likewise, taking numerical time pling frequency compared to, e.g., inertial sensors used (see derivatives, the calculation requires quantities of at least two Sect. 4.6). discrete points in time. For MR-AD, the calculation requires However, the numerical time derivative of IR-MD cal- quantities of one instant in time and no linear velocity vectors. culates even larger gradients for the momentum in terms Instead, the complexity of the calculation increases. These of magnitude caused by the combination of the noise- amplifying numerical differentiation and the white noise in the sensor signals. Figure 9 shows ﬁnally that the inﬂuence of the larger inaccuracies due to the white noise in the sensor signals in combination with the numerical time derivatives x -coordinate [m] y -coordinate [m] Autonomous Robots dominates. The relative error increases signiﬁcantly com- between the accelerometers. For example, if the distance r pared to the ideal reference. In contrast, MR-AD beneﬁts tends to zero, the measured values a and a converge. The 2z 0z through higher sample rates for sensors and actuator con- noise behavior is then the dominant difference between these trollers due to increased accuracy of q ˙ and q ¨ . It is further measurements. Thus ω ˙ would be mostly inﬂuenced by the i i x signiﬁcantly less impacted through white noise than the MD noise. method, as seen in Fig. 8. Here, we want to discuss the use of three additional BasedontheresultsshowninFig. 9, we propose to execute accelerometers to determine all joint velocities q ˙ and accel- this calculation directly in the actuator controllers. In general, eration q ¨ of the robot. The above mentioned arguments the frequency of the ﬁeld bus from robots depends mostly against numerical differentiation apply analogously. One on the total number of motors and sensors connected to the consequence would be to omit numerical differentiation ﬁeld bus. Typically frequencies of ﬁeld buses (e.g., CAN or entirely when calculating joint kinematics. Further, it is thus EtherCAT) for robots with more than 30 degrees of freedom independent of the actuator controller and ﬁeld bus frequen- are in the range of 100–1000 Hz. cies. A disadvantage is the increasing number of installed For physical robots, the higher frequency of the actuator sensors and the transmitted data. controllers is an advantage regarding the accuracy of q ˙ and Fourth, MR-AD and all analytical derivative methods are q ¨ . The errors of the numerical derivative approximation of based on the same basic physical laws. Provided that consis- q ˙ and q ¨ decrease due to smaller t. However, numerical tent inputs are given for the calculation, the results must be i i differentiation with t tending to zero generally ampliﬁes identical. Compared to our MR-AD, the mentioned analyt- signal noise. Therefore, e.g., implemented ﬁltering methods ical derivative methods require the absolute position vector in the actuator controller have to compensate for the noise- of the root segment from a spatially ﬁxed reference point as amplifying effects. Through the signiﬁcantly higher actuator an input. controller frequency, we can design ﬁltering methods with Henze et al. (2016) estimate the absolute position vector adapted delays so that they have no effects on the dynamics of the root segment by sensor fusion. Therefore, this kind calculations executed in ﬁeld bus frequency. Consequently, of estimation requires the kinematics as well as the absolute in this case we do not use the t of the ﬁeld bus frequency position vector and orientation of each end effector in contact. in any calculation step. Therefore, it has no impact on the One assumption required for the proposed sensor fusion is a accuracy of the dynamics calculation. rigid contact of at least one body segment. Third, in contrast to the MD method, analytical deriva- Despite the redundant calculation and the introduced tive EOMs generally require the angular acceleration vector averaging method in the approach of Henze et al. (2016), of the root segment. In this work, we compare two different this sensor fusion always contains inaccuracies due to non- approaches to assess angular acceleration. On the one hand, a present rigid contact and white noise of the IMU sensor numerical derivative of the angular velocity signal, measured signals. Also, the continuous determination of the absolute by IMU, is used. On the other hand, we present a measure- position vectors of all end effectors contains inaccuracies, ment through three additional uniaxial accelerometers. especially on unknown, uneven terrain and through undeﬁned For calculating the angular acceleration of the root seg- contacts. This kind of estimation for the absolute position ment, a noise-amplifying numerical differentiation has the vector of the root segment is, in our view, error-prone. same disadvantages as mentioned before. The increment t The approach of Buschmann (2011) avoids numerical used depends on the sampling rate of the sensor. Noise in inaccuracies through large distances from the initial reference the estimated angular velocity signal combined with high point to the robot by shifting the reference frame depending sampling rates could result in unrealistic angular accelera- on the gait cycle. In general, this reference frame does not tion approximations. Figure 10 conﬁrms this. Low sampling coincide with the BFF of the stance foot and the origin is rates would lead to inaccuracy for angular acceleration due ﬁxed. This means the position vector of the COM cannot to the increasing t between two measured quantities. A be calculated only through kinematics on physical robots. decreasing t through higher sampling rates increases the Buschmann does not propose an estimation of the global noise-amplifying effects of the differential quotient used. position vector of the COM, since the work aims to simulate Since the calculation method from three additional robots. accelerometers does not involve numerical differentiation, In contrast, our approach leads to high ﬂexibility in choos- the accuracy gets independent of t or the physical sensor ing the reference point and the reference frame for dynamics update frequency. The disadvantages previously described do calculations. The reference point can be any material point, not apply. However, this includes four additional noisy sig- e.g., the origin of the root segment or any immaterial point. nals in the calculation, which affects the measurement of the The state of motion of this point compared to the other mate- angular acceleration. Equations 27, 28, and 29 show that the rial points of the robot is irrelevant. This leads to the omission components of angular acceleration depend on the distance of estimating absolute position vectors. We made this effort 123 Autonomous Robots to avoid all estimations since they can cause inaccuracies Sixth, if we compare the IZMP coordinates of the physical in the results and for saving computational capacity. E.g., robot shown in Fig. 12 with the simulation results shown in if we choose the origin of the root segments BFF as the Fig. 6, we recognize generally similar curve developments. moving reference point, we can calculate relative position The main differing characteristics of the coordinate curves vectors of arbitrary material points through kinematics with- from physical robot to simulation are, on the one hand, the out inaccuracies due to uneven ground or undeﬁned contact. more inconstant development and the noticeably smoother Consequently, conditions with no rigid body contact of at transition between the change of the stance foot. The joint least one segment or in ﬂight phases do not disturb our used backlash present and measurement errors of the encoders EOM evaluation approach. on physical robots, the elasticities, and manufacturing and Furthermore, the spatial distance between the evaluation assembly inaccuracies mainly cause these two characteris- frame and all BFFs stays always in a similar order of magni- tics. In conclusion, despite the mechanical inadequacies and tude. Similarly, the recommendation of Featherstone (2008) white noise of all sensor signals, MR-AD is suitable to calcu- for the distance of one to two radii of gyration is con- late the IZMP also for physical robots. If the proposed method tinuously satisﬁed. This reduces numerical inaccuracy. We for measuring angular acceleration vectors is applied to all further use an algorithm for calculating the relative derivative segments of the robot, we expect that the results of dynamics from numerically available values, which allows us to derive calculation on real robots will be further improved. physical quantities in moving frames. Speciﬁcally, we use it Further, Fig. 12 shows the justiﬁcation that omitting the to calculate the relative time derivative of the angular velocity measurement errors of the encoders and the joint backlash is for all body segments. With this, we evaluate the dynamics sensible in the theoretical comparisons of the applied meth- for each body segment in its BFF. This is an efﬁcient way to ods. The impact of the measurement errors of the encoders, evaluate the EOMs according to Featherstone (2008). the joint backlash, and other mechanical inaccuracies on the Fifth, a further advantage of using our approach is seen for IZMP curves is much lower than the impact of the numerical the IZMP calculation walking, e.g., on slopes. In general, one differentiation of the momentum equations with noisy linear axis of the evaluation frame must point parallel to the normal and angular velocity signals. vector of the ground surface for calculating the IZMP. This As we mentioned previously, there are disturbances in the is required to resolve the under-determination of the system dynamics calculation in general. These include inaccuracies of coordinate equations from Eq. 22 through determining the in the mechanical model compared to the real robot, e.g., due z-coordinate. This demand cannot be fulﬁlled by the com- to deviating mounting positions of the joint axes, inertia data pared approaches. Consequently, they do not provide usable or the underlying rigid body hypothesis, to name a few. To results without further adjustments in calculating the IZMP reduce the impact of inaccuracies, we propose combining a on slopes. For valid results, the compared approaches would sensor fusion from our approach with a 6-axis force/torque require a transition of the inertial frame. Likewise, the equa- sensor and adding a Kalman ﬁlter. A limitation of the sensor tion to calculate the IZMP, as, e.g., shown by Kajita et al. fusion is that it only provides valid data for a bipedal walking (2014), also has to be adapted to the rotated gravity vector. without external perturbations applied to the robot. For evaluating the EOMs in a MR, this condition for the IZMP calculation remains fulﬁlled. Here, the evaluation frame rotates in a suitable orientation to cover the slope to 7 Conclusion and future work ensure that the required axis parallelism remains fulﬁlled. IMUs placed in the feet measure the orientation of the sur- From the perspective of the authors, the combination of face during the stance phase in each gait cycle. evaluating the EOMs formulated through the Newton–Euler The comparison of the resulting IZMP coordinates for equations with respect to a moving reference point in a NIF the level ground and the slope shows that they are located and measuring the angular acceleration vector of the root within the support polygons as seen in Fig. 11. Our approach segment contributes to an efﬁcient and robust evaluation evaluates the locomotion on the slope as stable. of the inverse dynamics in robotics. The entire proposed Consequently, with our approach, the reference point is approach produces comparably similar results as the previ- adjustable to the implemented control algorithm. E.g., if we ous approaches evaluating the EOMs with respect to a ﬁxed use an IZMP-based control algorithm, the moving reference point in the inertial frame or the center of mass for a humanoid could coincide with the surface plane to be independent of robot, with the main beneﬁt that the proposed approach is less friction forces. Or, if we control the ground reaction forces, prone to noise than the previous approaches. it may be efﬁcient to place the reference point into the mov- We conﬁrm that in theory and by comparing it with the ing root segment. The calculated quantities can thus directly mentioned approaches using simulation results. In theory, trigger a reaction without additional spatial transformation. we show the conceptual advantages of our approach. E.g., 123 Autonomous Robots it gets rid on the dependency of any global position vector. tarity problems. Nonlinear Dynamics, 14(3), 231–247. https://doi. org/10.1023/A:1008292328909 Further, we show the ﬂexibility of choosing the reference Buschmann, T. (2011). Simulation and control of biped walking robots point and the evaluation frame concerning the preferences of (1 Auﬂ ed.). München: Verl. Dr. Hut. the control algorithm used. Englsberger, J. , Ott, C. , Roa, M.A. , Albu-Schaffer, A. , & Hirzinger, G. The comparison of simulation results clearly shows fur- (2011). Bipedal walking control based on Capture Point dynamics. In 2011 IEEE/RSJ international conference on intelligent robots ther advantages for the application to physical robots. On the and systems (pp. 4420–4427). CAIEEE. https://doi.org/10.1109/ one hand, we reduce the quantity of noise-amplifying numer- IROS.2011.6094435 ical differentiation. On the other hand, with our approach, Featherstone, R. (2008). Dynamics of rigid body systems. R. Feath- we can calculate the IZMP for the gait from the level to the erstone (Ed.), Rigid body dynamics algorithms (pp. 39–64). Springer. https://doi.org/10.1007/978-1-4899-7560-7_3 inclined surface. We also show that our algorithm calculates Goswami, A. (1999a). Foot rotation indicator (FRI) point: A new usable results applied to physical robots. gait planning tool to evaluate postural stability of biped robots. It is also worth mentioning that this approach is not lim- In Proceedings 1999 IEEE international conference on robotics ited to bipedal robots or humanoid robots in general. Due to and automation (Cat. No.99CH36288C) (vol. 1, pp. 47–52). USAIEEE. https://doi.org/10.1109/ROBOT.1999.769929 the generic structure, the EOMs can be calculated for multi- Goswami, A. (1999b). Postural stability of biped robots and the legged robots or human beings. The adaption to KDL and the foot-rotation indicator (FRI) point. The International Journal integration in ROS2 of this algorithm support straightforward of Robotics Research, 18(6), 523–533. https://doi.org/10.1177/ usability. 02783649922066376 Goswami, A. , & Kallem, V. (2004). Rate of change of angular Next, we will equip all segments of Sweaty with one sen- momentum and balance maintenance of biped robots. In IEEE sor unit combining an IMU and three additional uniaxial international conference on robotics and automation, 2004 pro- accelerometers to eliminate all numerical differentiation and ceedings (ICRA’04). 2004 (vol.4, pp. 3785–3790). USAIEEE. use our approach to evaluate the stability of the robot Sweaty. https://doi.org/10.1109/ROBOT.2004.1308858 Henze, B., Roa, M. A., & Ott, C. (2016). Passivity-based whole-body As a second step, we will setup up the sensor fusion with the balancing for torque-controlled humanoid robots in multi-contact self-developed 6-axis force/torque sensor and we will try to scenarios. The International Journal of Robotics Research, 35(12), implement a Kalman ﬁlter. Further, we will try to use the 1522–1543. https://doi.org/10.1177/0278364916653815 sensor unit also to evaluate the stability of human beings. Hirukawa, H. , Hattori, S. , Harada, K. , Kajita, S. , Kaneko, K. , Kane- hiro, F., & Morisawa, M. (2006). A universal stability criterion Finally, we will use this approach as the basis of a con- of the foot contact of legged robots—Adios ZMP. In Proceedings trol algorithm to stabilize both static situations and dynamic 2006 IEEE international conference on robotics and automation, locomotion. We want to align this controller speciﬁcally to 2006 (ICRA 2006) (pp. 1976–1983). USAIEEE. https://doi.org/ the control paradigms of human beings, taking beneﬁt of the 10.1109/ROBOT.2006.1641995 InJoon, M. , DongHa, Y. , MinSung, A. , & Jeakweon, H. (2020). Body identiﬁed characteristics of our approach. trajectory generation using quadratic programming in bipedal robots. In 2020 20th international conference on control, automa- Funding Open Access funding enabled and organized by Projekt tion and systems (ICCAS) (pp. 251–257). IEEE. https://doi.org/10. DEAL. The authors have no relevant ﬁnancial or non-ﬁnancial interests 23919/ICCAS50221.2020.9268204 to disclose. The project Sweaty of the University of Applied Sciences Kajita, S., Hirukawa, H., Harada, K., & Yokoi, K. (2014). Introduction Offenburg is a non-proﬁt project. The objective is to contribute to the to humanoid robotics (vol. 101). Springer. https://doi.org/10.1007/ research ﬁeld of humanoid robotics. Regional companies sponsor equip- 978-3-642-54536-8 ment. However, there is no inﬂuence on or exchange of research results. Kajita, S., Kanehiro, F., Kaneko, K., Fujiwara, K., Harada, K., Yokoi, K., & Hirukawa, H. (2003). Resolved momentum control: Humanoid Open Access This article is licensed under a Creative Commons motion planning based on the linear and angular momentum. In Attribution 4.0 International License, which permits use, sharing, adap- Proceedings 2003 IEEE/RSJ international conference on intelli- tation, distribution and reproduction in any medium or format, as gent robots and systems (IROS 2003) (cat. no.03ch37453) (vol. 2, long as you give appropriate credit to the original author(s) and the pp. 1644–1650). https://doi.org/10.1109/IROS.2003.1248880 source, provide a link to the Creative Commons licence, and indi- Kajita, S., Morisawa, M., Miura, K., Nakaoka, S., Harada, K., Kaneko, cate if changes were made. The images or other third party material K., & Yokoi, K. (2010). Biped walking stabilization based on linear in this article are included in the article’s Creative Commons licence, inverted pendulum tracking. In 2010 IEEE/RSJ international con- unless indicated otherwise in a credit line to the material. If material ference on intelligent robots and systems (pp. 4489–4496). IEEE. is not included in the article’s Creative Commons licence and your https://doi.org/10.1109/IROS.2010.5651082 intended use is not permitted by statutory regulation or exceeds the Kamioka, T., Kaneko, H., Takenaka, T., & Yoshiike, T. (2018). Simul- permitted use, you will need to obtain permission directly from the copy- taneous optimization of ZMP and footsteps based on the analytical right holder. To view a copy of this licence, visit http://creativecomm solution of divergent component of motion. In 2018 IEEE inter- ons.org/licenses/by/4.0/. national conference on robotics and automation (ICRA) (pp. 1763–1770). IEEE. https://doi.org/10.1109/ICRA.2018.8460572 Kimura, K. , Nozawa, S., Mizohana, H., Okada, K., & Inaba, M. (2018). Riding and speed governing for parallel two-wheeled scooter based References on sequential online learning control by humanoid robot. In 2018 IEEE/RSJ international conference on intelligent robots and sys- tems (IROS) (pp. 1–9). IEEE. https://doi.org/10.1109/IROS.2018. Anitescu, M., & Potra, F. A. (1997). Formulating dynamic multi-rigid- body contact problems with friction as solvable linear complemen- 8593685 123 Autonomous Robots Mehrtens, G. (1887). Technische Mechanik fester, ﬂüssiger und luftför- Wang, P., Liu, J., Zha, F., Guo, W., Wang, X., Li, M., & Sun, L. miger Körper(vol. 3). Kommissions-Verlag v. E. Toeche. (2017). A velocity estimation algorithm for legged robot. Advances Mesesan, G., Englsberger, J., & Ott, C. (2021). Online DCM trajectory in Mechanical Engineering, 9(12), 1687–8140. https://doi.org/10. adaptation for push and stumble recovery during humanoid loco- 1177/1687814017732736 motion. In 2021 IEEE international conference on robotics and Wittenburg, J. (1977). Dynamics of systems of rigid bod- automation (ICRA) (pp. 12780–12786). IEEE. https://doi.org/10. ies. Vieweg+Teubner Verlag. https://doi.org/10.1007/978-3-322- 1109/ICRA48506.2021.9560808 90942-8 Padgaonkar, A. J., Krieger, K. W., & King, A. I. (1975). Measurement Woernle, C. (2016). Mehrkörpersysteme: eine Einführung in die Kine- of angular acceleration of a rigid body using linear accelerometers. matik und Dynamik von Systemen starrer Körper (2., erweiterte Journal of Applied Mechanics, 42(3), 552–556. https://doi.org/10. Auﬂage ed.). Springer Vieweg. 1115/1.3423640 Pang, J. , & Trinkle, J. (2000). Stability characterizations of ﬁxtured Publisher’s Note Springer Nature remains neutral with regard to juris- rigid bodies with coulomb friction. In Proceedings 2000 icra. mil- dictional claims in published maps and institutional afﬁliations. lennium conference. IEEE international conference on robotics and automation (vol. 1, pp. 361–368). IEEE. https://doi.org/10. 1109/ROBOT.2000.844083 Pratt, J., Carff, J., Drakunov, S., & Goswami, A. (2006). Capture Maximilian Gießler received his Point: A step toward humanoid push recovery. In 2006 6th IEEE- B.Eng. in 2018 and M.Sc. in 2020 RAS international conference on humanoid robots (pp. 200–207). in Mechanical Engineering from University of Genova, IEEE. https://doi.org/10.1109/ICHR.2006. the Offenburg University of App- lied Sciences. He is currently a Saida, T., Yokokohji, Y., & Yoshikawa, T. (2003). FSW (fea- Ph.D. student and employed as a sible solution of wrench) for multi-legged robots. In 2003 scientiﬁc assistant at the Offen- IEEE international conference on robotics and automation (Cat. burg University of Applied Sci- No.03CH37422) (vol. 3, pp. 3815–3820). IEEE. https://doi.org/ ences. His research interests inclu- 10.1109/ROBOT.2003.1242182 de multibody dynamics and Sardain, P., & Bessonnet, G. (2004). Forces acting on a biped robot. humanoid robotics, with an Center of pressure–Zero moment point. IEEE Transactions on Sys- emphasis on bipedal locomotion tems, Man, and Cybernetics - Part A: Systems and Humans, 34(5), and stability theory. 630–637. https://doi.org/10.1109/TSMCA.2004.832811 Vijaykumar Shah, S., Kumar Saha, S., & Kumar Dutt, J. (2013). Dynam- ics of tree-type robotic systems (vol. 62). Springer. Vukobratovic, M., Borovac, B., & Šurdilovic, ´ D. (2001). Zero-moment Bernd Waltersberger has been point—Proper interpretation and new applications. In IEEE-RAS a professor at the Faculty of international conference on Humanodi robots (pp. 237–244). Mechanical and Process Engineer- IEEE. ing at the Offenburg University of Vukobratovic, M., Frank, A. A., & Juricic, D. (1970). On the stability Applied Sciences since 2010. He of biped locomotion. IEEE Transactions on Biomedical Engi- received his diploma in mechani- neeringBME, 17(1), 25–36. https://doi.org/10.1109/TBME.1970. cal engineering in 2002 from the University of Karlsruhe (now Vukobratovic, ´ M., & Stepanenko, J. (1972). On the stability of anthro- K.I.T) and was awarded his doc- pomorphic systems. Mathematical Biosciences, 15(1–2), 1–37. torate there in 2007 for a the- https://doi.org/10.1016/0025-5564(72)90061-2 sis on structural mechanics. He Wang, H. , Tian, Z., Hu, W., & Zhao, M. (2018). Human-like ZMP subsequently worked as a devel- generator and walking stabilizer based on divergent component opment engineer in the automo- of motion. In 2018 IEEE-RAS 18th international conference on tive ﬁeld and is co-inventor on 7 humanoid robots (Humanoids) (pp. 1–9). IEEE. https://doi.org/ patent applications. His research 10.1109/HUMANOIDS.2018.8624926 interests are in the dynamic interaction of machine elements and multibody dynamics, especially in the area of vibration engineering and stability theory. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Autonomous Robots Springer Journals http://www.deepdyve.com/lp/springer-journals/robust-inverse-dynamics-by-evaluating-newton-euler-equations-with-ArZZ06fYee

Loading next page...

References (32)

A. Padgaonkar, K. Krieger, A. King (1975)
Measurement of Angular Acceleration of a Rigid Body Using Linear Accelerometers
Journal of Applied Mechanics, 42
Riding and speedgoverning for parallel two - wheeled scooter based on sequential online learning control by humanoid robot
H. Hirukawa, Shizuko Hattori, K. Harada, S. Kajita, K. Kaneko, F. Kanehiro, K. Fujiwara, M. Morisawa (2006)
A universal stability criterion of the foot contact of legged robots - adios ZMP
Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006.
A. Goswami, V. Kallem (2004)
Rate of change of angular momentum and balance maintenance of biped robots
IEEE International Conference on Robotics and Automation, 2004. Proceedings. ICRA '04. 2004, 4
J. Pang, J. Trinkle (2000)
Stability characterizations of fixtured rigid bodies with Coulomb friction
Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065), 1
T. Saida, Y. Yokokohji, T. Yoshikawa (2003)
FSW (feasible solution of wrench) for multi-legged robots
2003 IEEE International Conference on Robotics and Automation (Cat. No.03CH37422), 3
M. Anitescu, F. Potra (1997)
Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementarity Problems
Nonlinear Dynamics, 14
M. Vukobratovic, A. Frank, D. Juricic (1970)
On the stability of biped locomotion.
IEEE transactions on bio-medical engineering, 17 1
S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, H. Hirukawa (2003)
Resolved momentum control: humanoid motion planning based on the linear and angular momentum
Proceedings 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003) (Cat. No.03CH37453), 2
Pengfei Wang, Jikai Liu, F. Zha, Wei Guo, Xin Wang, Mantian Li, Lining Sun (2017)
A velocity estimation algorithm for legged robot
Advances in Mechanical Engineering, 9
M Vukobratović, J Stepanenko (1972)
On the stability of anthropomorphic systems
Mathematical Biosciences, 15
(2016)
Mehrkörpersysteme: eine Einführung in die Kinematik und Dynamik von Systemen starrer Körper
D. Chevallier, J. Lerbet (2018)
Dynamics of Rigid Body Systems
Johannes Englsberger, C. Ott, M. Roa, A. Albu-Schäffer, G. Hirzinger (2011)
Bipedal walking control based on Capture Point dynamics
2011 IEEE/RSJ International Conference on Intelligent Robots and Systems
Haitao Wang, Zhongyuan Tian, Wenbin Hu, Mingguo Zhao (2018)
Human-Like ZMP Generator and Walking Stabilizer Based on Divergent Component of Motion
2018 IEEE-RAS 18th International Conference on Humanoid Robots (Humanoids)
M. Vukobratovic, J. Stepanenko (1972)
On the stability of anthropomorphic systems
Bellman Prize in Mathematical Biosciences, 15
J. Pratt, John Carff, S. Drakunov, A. Goswami (2006)
Capture Point: A Step toward Humanoid Push Recovery
2006 6th IEEE-RAS International Conference on Humanoid Robots
Takumi Kamioka, Hiroyuki Kaneko, T. Takenaka, T. Yoshiike (2018)
Simultaneous Optimization of ZMP and Footsteps Based on the Analytical Solution of Divergent Component of Motion
2018 IEEE International Conference on Robotics and Automation (ICRA)
(2001)
Zero - moment point — Proper interpretation and new applications
J. Wittenburg, P. Likins (1977)
Dynamics of systems of rigid bodies
S. Kajita, H. Hirukawa, K. Harada, K. Yokoi (2014)
Introduction to Humanoid Robotics
, 101
A. Goswami (1999)
Foot rotation indicator (FRI) point: a new gait planning tool to evaluate postural stability of biped robots
Proceedings 1999 IEEE International Conference on Robotics and Automation (Cat. No.99CH36288C), 1
Min. InJoon, Yoo. DongHa, Ahn. MinSung, Han. Jeakweon (2020)
Body Trajectory Generation Using Quadratic Programming in Bipedal Robots
2020 20th International Conference on Control, Automation and Systems (ICCAS)
T. Buschmann (2010)
Simulation and Control of Biped Walking Robots
A. Goswami (1999)
Postural Stability of Biped Robots and the Foot-Rotation Indicator (FRI) Point
The International Journal of Robotics Research, 18
P. Sardain, G. Bessonnet (2004)
Forces acting on a biped robot. Center of pressure-zero moment point
IEEE Transactions on Systems, Man, and Cybernetics - Part A: Systems and Humans, 34
George Mesesan, Johannes Englsberger, C. Ott (2021)
Online DCM Trajectory Adaptation for Push and Stumble Recovery during Humanoid Locomotion
2021 IEEE International Conference on Robotics and Automation (ICRA)
M Anitescu (1997)
231
Nonlinear Dynamics, 14
Bernd Henze, M. Roa, C. Ott (2016)
Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios
The International Journal of Robotics Research, 35
(1994)
On the Stability of
S. Kajita, M. Morisawa, K. Miura, Shin'ichiro Nakaoka, K. Harada, K. Kaneko, F. Kanehiro, K. Yokoi (2010)
Biped walking stabilization based on linear inverted pendulum tracking
2010 IEEE/RSJ International Conference on Intelligent Robots and Systems
Technische Mechanik fester, flüssiger und luftförmiger Körper

Publisher: Springer Journals
Copyright: Copyright © The Author(s) 2023
ISSN: 0929-5593
eISSN: 1573-7527
DOI: 10.1007/s10514-023-10092-x
Publisher site: See Article on Publisher Site

Abstract

Maintaining stability while walking on arbitrary surfaces or dealing with external perturbations is of great interest in humanoid robotics research. Increasing the system’s autonomous robustness to a variety of postural threats during locomotion is the key despite the need to evaluate noisy sensor signals. The equations of motion are the foundation of all published approaches. In contrast, we propose a more adequate evaluation of the equations of motion with respect to an arbitrary moving reference point in a non-inertial reference frame. Conceptual advantages are, e.g., getting independent of global position and velocity vectors estimated by sensor fusions or calculating the imaginary zero-moment point walking on different inclined ground surfaces. Further, we improve the calculation results by reducing noise-amplifying methods in our algorithm and using speciﬁc characteristics of physical robots. We use simulation results to compare our algorithm with established approaches and test it with experimental robot data. Keywords Multi-body dynamics · Equation of motion · Moving reference point · Imaginary zero-moment point · Humanoid robot 1 Introduction examples are household chores, care of the elderly, and use as a service worker. Humanoid robots have a signiﬁcant advantage over all other The prerequisite of control algorithms for maintaining types of robots. Since they are designed to resemble the stability is the calculation of the robot’s dynamics. Euler’s human body, they are suitable for interacting in a human- theorem is often used in robotics as a starting point for made habitat. Consequently, the research ﬁeld of humanoid dynamics calculation, e.g. Henze et al. (2016), Buschmann robotics is currently gaining importance. A crucial aspect of (2011), Hirukawa et al. (2006) and Kajita et al. (2014)use increasing the autonomy of humanoid robots is to improve it. This theorem states that the absolute time derivative of dynamic stability during locomotion in challenging settings, the total angular momentum with respect to a ﬁxed reference e.g., the robot’s ability to effectively tolerate external pertur- point 0 equals the total resultant moment about the same (0) (0) bations during dynamic, omnidirectional locomotion. On the ˙ reference point (L = M ). other hand, locomotion in arbitrary environments is of great The approaches from Henze et al. (2016) and Buschmann relevance. In this context, the main focus is to enhance the (2011) require methods to determine at least the absolute robot’s ability to cope with uneven and unknown ground sur- position vector of one body segment. If this vector is not faces of the man-made habitat. Given that stability-enhancing explicitly measurable by sensors, it is estimated by sensor control algorithms are implemented, the robot’s autonomy fusion. E.g., Henze et al. (2016) use the kinematics, absolute will increase. This opens up new areas of application. Some position vector and the orientation provided by IMUs of all end effectors in contact. Buschmann (2011) uses the odome- Maximilian Gießler try. In general, sensor fusions or odometry can be sensitive to maximilian.giessler@hs-offenburg.de errors, because of signal noise or undeﬁned contact impacts Bernd Waltersberger the accuracy. Further, the proposed sensor fusion from Henze bernd.waltersberger@hs-offenburg.de et al. (2016) requires the assumption of rigid contact of at Department of Mechanical and Process Engineering, least one body segment. Likewise, a stand-alone calculation University of Applied Sciences Offenburg, 77652 Offenburg, of absolute position vectors from the kinematic chains of the Germany 123 Autonomous Robots robot is not feasible. Kajita et al. (2014) calculate the angular the calculation of the remaining kinematic quantities for the (0) momentum L and propose to use a numerical differentia- following kinematic chain using the relative derivative and (0) tion to determine its time derivative. To calculate L both the the kinematic state of all joints. absolute position and linear velocity vector are required. Due Moreover, we use simulation and experimental results of to the mentioned aspects, we regard the requirement of the the humanoid robot Sweaty to show our approach’s capa- absolute position or linear velocity vector as the main prob- bilities applied on the physical robots. Section 5 illustrates lem of these published approaches. This especially applies the outcomes. In the article’s ﬁnal parts, we discuss the dif- to the need to handle undeﬁned contacts and noisy sensor ferences and impacts on physical robots of our approach signals on real robots. compared to Henze et al. (2016), Buschmann (2011) and Also, a human being has no sense organs to directly per- Kajita et al. (2014) and give an outlook to future work. ceive its global position or absolute velocity with respect to a ﬁxed spatial point. Further, it does not know the position vec- tor of its center of mass (COM). In conclusion, human control 2 Previous work algorithms for stabilizing locomotion can operate without this information. Several approaches to describe the stability of humanoid Considering all mentioned aspects, we propose an robots during locomotion have been published. Mostly, their approach that does not require any information about the aim was to determine an indicator assessing stability in global position, velocity, or COM vector of the robot. We dynamic situations quantitatively. One of the most discussed achieve this by evaluating the Newton–Euler equations of and applied stability criteria uses the notion of the zero- motion (EOM) with respect to an arbitrary moving reference moment point (ZMP). If the ZMP lies within the support point in a non-inertial reference frame (NIF). The dynamics region, the current dynamic equilibrium of the robot can be calculation is then independent of absolute position vectors assumed stable in the sense that the robot cannot tip over of the robot. Similar approaches using a moving reference an edge of the support region. The ZMP is the spatial point point are common in multi-body dynamics, as seen e.g. in in the ground surface plane where the components of the Wittenburg (1977). moment vector resulting from the reaction force parallel to Furthermore, one goal is to provide an approach that is this plane are zero. Vukobratovic et al. (1970) and Vuko- usable on real robots without ﬁltering the sensor signal. bratovic´ and Stepanenko (1972) introduced this concept to Therefore, we propose to measure the angular acceleration the ﬁeld of robotics. Their approach can be regarded as a vector via angular velocity and linear acceleration measure- generalization of the classical static tilting stability concept ments to reduce the number of noise-amplifying numerical for a rigid body, e.g. Mehrtens (1887), where tipping over differentiation. So this work contributes to the practical use edges of the support region due to external forces is prevented of evaluating the EOMs on physical robots with noisy sensor when the resultant of the externally applied forces and the signals. body’s weight applied to the effective point of action inter- For this work, we choose the imaginary zero-moment sects the support region. Algorithms for calculating the ZMP point (IZMP) to evaluate the results of the EOM calcula- are shown by, e.g., Kajita et al. (2014). The work of Sardain tion. We ﬁnd the IZMP gives additional valuable insight and Bessonnet (2004) shows that the ZMP coincides with the into the stability evaluation of gait behaviors compared to COP under the assumption that the soles are in ﬂat contact the sole evaluation of the position of the center of pressure with one single surface and the robot is free of external per- (COP) measured with force/torque sensors. Especially for turbation forces and moments. If the theoretically calculated gait behaviors with possible external perturbations, calculat- ZMP differs from the COP then this is obviously caused by ing the distance of the IZMP from the support region of the perturbation forces or moments, since such perturbations are robot’s feet or as well from the COP could give additional not considered in the kinematical ZMP calculation. Vukobra- quantitative information as to the impact of the perturbation tovic et al. (2001) hence proposed the name IZMP for this or control strategies. conﬁguration of the ZMP, to which we adhere in our work. In the following Sect. 2, we review the state of the art of Several other concepts exist to quantify stability using referring control algorithms and stability criteria and point a point in space. For example, there is the foot-rotation- out where we see further contributions achieved through our indicator by Goswami (1999a, 1999b), the feasible-solution- approach. In Sect. 3, we explain the mathematical founda- of-wrench by Saida et al. (2003), the zero-rate-of-change tions. In Sect. 4, we discuss the methods used for providing angular momentum by Goswami and Kallem (2004)orthe the data and sensor signals required for the core algorithm capture point (CP) by Pratt et al. (2006). Furthermore, meth- (seen in Sect. 3). On the one hand, we describe a way to ods were developed to modulate contact state and control measure the total kinematic state, especially the angular the ground reaction forces. Hirukawa et al. (2006), Anitescu acceleration of the root segment. On the other hand, we show and Potra (1997), Pang and Trinkle (2000) propose numer- 123 Autonomous Robots ical methods to search for feasible solutions and, based on touchdown introduces a new reference frame. It is placed their characteristics, predictions are made about the stability below the stance foot. This avoids numerical inaccuracies proposed. through large distances from the initial reference point. All the mentioned proposals use the EOMs as a founda- Finally, Buschmann (2011) also uses a dynamics calcula- tion. The level of detail of the mechanical model is essential tion based on a ﬁxed reference point in the current reference for the accuracy. Different approaches have been established frame. over the last decades. On the one hand, there are concep- In contrast, we propose an algorithm to evaluate the EOMs tual models such as the linear inverted pendulum (LIP), e.g., using Newton–Euler equations for a rigid multi-body system shown in Kajita et al. (2010). The use case is to implement related to an arbitrary moving reference point in a NIF. To simpliﬁed control algorithms in combination with stability maintain generality, we use the ﬂoating base model. We use criteria. For example, Englsberger et al. (2011) present a the rigid body theory typically known from the ﬁeld of multi- bipedal walking control approach based on the LIP and the body dynamics to derive the EOMs. We further use the so- CP dynamics. Mesesan et al. (2021) introduced a trajectory called relative derivative to evaluate the EOMs referenced in adaptation for push and stumble recovery based on the CP. a NIF as well as to calculate all kinematic quantities for the Kamioka et al. (2018) are optimizing the ZMP and footsteps multi-body segments. based on the CP for a LIP. Also, Wang et al. (2018) and InJoon In combination, this will allow us to carry the reference et al. (2020) use the concepts of LIP and ZMP to stabilize the point and the NIF along with arbitrary locomotion of the gait of humanoid robots. On the other hand, modeling robots robot. E.g., the reference point can be an arbitrary spatial as rigid multi-body models has become established, see, e.g., point or ﬁxed to any material or immaterial point of the robot. in Vijaykumar Shah et al. (2013). Featherstone (2008)shows Thus, for different types of control algorithms, an adapted the ﬂoating-base model as a kind of rigid multi-body model reference point can be used. Also, we will show that our and it is deﬁned most generally. algorithm requires no absolute position and linear velocity A further distinctive criterion is the calculation method of vector of the root segment. Thus, our approach omits the mea- the time derivatives of the EOMs. Either numerical (Kajita surement or estimation by sensor fusion to determine those. et al., 2014) or analytical derivatives of the equations are To evaluate the effects when applying these approaches proposed. For the numerical derivatives of the linear and on physical robots, we will simulate the inﬂuence of realistic angular momentum, the knowledge of the linear velocity sensor noise behavior. We will also show the impact of feasi- of at least one segment is required. Wang et al. (2017)or ble sensor and actuator controller frequencies. In addition, we Kimura et al. (2018) show approaches to estimate the linear will demonstrate a way to measure the angular acceleration velocity vector. As mentioned before, mostly the reference of the root segment using gyroscope and linear acceleration point for the balance of angular momentum is referenced to data. We will compare this method against the numerical a ﬁxed point in the inertial frame (WF) or the immaterial calculation from the angular velocity signal. COM. Kajita et al. (2014) propose to calculate the rate of Furthermore, using a moving reference point as the origin change of angular momentum by applying a numerical time of a moving evaluation frame allows following the contour derivative. Hirukawa et al. (2006) use the rate of change of of the ground’s surface. IMUs located in the feet provide angular momentum without specifying an explicit equation. information about the inclination of the surface during the Only the equation of angular momentum is stated explicitly. stance phase. A transition for changing surface inclinations Further, Kajita et al. (2003) also propose a method to control is no longer necessary for calculating, e.g., the IZMP. the total linear and angular momentum. Besides combining the mentioned approaches with our The proposed analytical derivative methods require also work, we intend to fully outline the computation procedure the absolute position vector of the root segment from a spa- of the dynamics and the IZMP of a robot omitting noise- tially ﬁxed reference point as an input. Henze et al. (2016) amplifying numerical differentiation for critical terms in the use a dynamics calculation based on a ﬁxed reference point in procedure and without estimations of global vectors. For an a WF. They estimate the absolute position vector of the root easy application, the backgrounds and the algorithm we are segment by sensor fusion. This estimation requires infor- going to present are fully adapted to use the Kinematics and mation about the kinematics, absolute position vectors and Dynamics Library (KDL) and to integrate into the Robot orientations of all end effectors in contact. In their approach, Operating System 2 (ROS2). they calculate the absolute position vector of the root seg- ment from each end effector in rigid contact. The redundant results are combined in a sensor fusion to estimate global 3 Backgrounds position vectors. Buschmann (2011) describes a method to introduce a At ﬁrst, we introduce in Fig. 1 the frames and vectors. NIF is new reference frame for every gait cycle. Each swing leg the frame, in which the EOMs are evaluated. Its origin rep- 123 Autonomous Robots the moving reference frame. The term ω × p describes the i0 rate of change of the vector due to the rotation of frame i. For the implemented EOMs, Newton’s second law and angular momentum theorem of Euler, are taken as starting point. Hence P = F (5) (0) (0) L = M . (6) (0) apply. Here, P = mv is the linear and L the angular COM momentum, which we discuss in more detail in the follow- (0) ing section. F is the total external force and M the total moment with respect to the origin of WF applied on the body. Fig. 1 Deﬁnition of the vectors and frames used to describe the kine- Due to the linearity of the mathematical operators used, the matic description of a rigid multi-body system balance of momentum and angular momentum for the whole robot can be calculated as the sum of the single rigid-body resents the arbitrary moving reference point for the angular segments. momentum balance. Each rigid-body part is given a body- In the following, the velocity v is the ﬁrst and the accel- ﬁxed frame (BFF). The origin of a BFF is located at the eration a is the second total time derivative of r .Inthis X X center of the joint axis considering the mechanical joint struc- context, X represents all the following different subscripts ture. The equivalent transformations used in the following are used, e.g. O or COM . Substituting the second total time i i based on the rigid body theory used in multi-body dynamics. derivative of Eq. 3 into Eq. 5 and by taking into account the The vector from the origin of the NIF to the corresponding law of the COM, results in BFF of an arbitrary segment i is called r . The vector from an arbitrarily positioned inertial frame to the origin of a BFF F = m a + ω ˙ × s + ω × (ω × s ) . (7) i i O i0 i i0 i0 i is denoted as r . The vector from the origin of the BFF to any mass element dm is declared as r . The vector s points from k i (0) For L , we start from the following equation the origin of the BFF to the segment COM and x points from the COM to an arbitrary mass element dm. The deﬁnitions (0) L = m r × v + ω × s i i O O i0 i i i of these vectors are given by (8) (O ) +s × v + ω , i O i0 r = r + r (1) O k proposed by Wittenburg (1977), who therefrom derived a r = r + r (2) O NIF H i i formulation with respect to an arbitrary moving reference r = r + s . (3) COM O i i i point. However, our approach differs so far that we will derive a formulation that assumes an arbitrary moving reference Here, r is the position vector pointing from the origin NIF point for the moments while using another arbitrary reference of a arbitrary WF to the origin of NIF. The absolute time point, e.g., pivot points for the inertia tensors. The latter is derivative of an arbitrary vector p = p e is calculated i β i β β a common practice in the generic description of input data, according to the basic equation of kinematics as follows e.g., in C++ libraries. Form now, i indicates the reference of the quantities to the dp (O ) dp de i β i β i ith body of a multi-body system. In Eq. 8, is the inertia = e + p i β i β dt dt dt tensor of the ith body with respect to the reference point O . β=x ,y,z β=x ,y,z (4) We substitute Eq. 2 and the total time derivative of Eq. 3 into dp = + ω × p. i0 Eq. 8. Further, we take the total time derivative of Eq. 8.This dt gives Here, p is the coordinate in direction β of p and e is the i β i β (0) unit vector of direction β.The sufﬁx i indicates an arbitrary L = m r + r × a i i NIF H COM i i frame in which the vector is referenced and ω indicates the i0 +s × a i O dp (9) frame’s angular velocity. The term represents the so called dt (O ) (O ) i i + ω ˙ + ω × ω i0 i0 i0 relative derivative. This term can be interpreted as the vector’s (0) rate of change that one would observe if one were ﬁxed to = M . 123 Autonomous Robots We see all absolute velocity terms cancel out. forces and moments, e.g., caused by the gravity force or the Next, we subtract the term r × m a from Eq. 9. ground reaction force remain. This results in NIF i COM Further, for the moment on body i and by using Eq. 5 (P) (0) F = (m (a + ω ˙ × s ext i O i0 i M = M − r × F NIF i i i i =1 (10) (0) = M − r × m a NIF i COM i + ω × (ω × s ))) (12) i0 i0 i (P) (P) applies in general. Here, M describes the total moment M = m (r × a + r × (ω ˙ × s ) i H O H i0 i ext i i i i =1 applied on body i with respect to the arbitrary moving refer- ence point P. We get ﬁnally + r × (ω × (ω × s )) + s × a ) H i0 i0 i i O i i (O ) (O ) i i (P) + ω ˙ + ω × ω . (13) i0 i0 i0 M = m (r × a + r × (ω ˙ × s ) i i i H O H i0 i i i i + r × (ω × (ω × s )) + s × a ) (11) H i0 i0 i i O i i The algorithm is recursively designed and thus generic. The (O ) (O ) i i + ω ˙ + ω × ω . i0 i0 i0 mechanical structure and the physical properties of each rigid body segment of the robot can be provided via a Uniﬁed Robot Equation 11 describes the change in angular momentum Description Format ﬁle. for an arbitrary moving reference point. From Eq. 11,we For the implementation in the code, we switch from vec- see that all absolute position and velocity vectors cancel out. tor to coordinate matrix notation. E.g., the representation of The Eqs. 7 and 11 describe the Newton–Euler EOMs for an Eq. 13 is given by arbitrary moving reference point in a NIF. (P) NIF NIF i i i M = m R ( r × a i H O ext i i 4 Methods for a robust and eﬃcient i =1 evaluation of the inverse dynamics i i i + r × ω ˙ × s H i0 i (14) i i i i With our approach, we intend to provide an effective cal- + r × ( ω × ( ω × s )) H i0 i0 i culation of the inverse dynamics. We use this to calculate i i i (O)i + s × a ) + ω ˙ i O i i i0 the IZMP exemplarily. We believe the IZMP provides an i i (O) i additional beneﬁt in evaluating gait situations under exter- + ω × ω . i0 i i0 nal perturbations compared to a COP measured using 6-axis force/torque sensors. For example, the distance of the IZMP NIF i Here, R is a 3 × 3 rotation matrix that transforms the from the tipping edge or the measured COP provides a metric for the needed response to stabilize the robot in the presence coordinates with respect to the frame i to the coordinates with respect to the frame NIF. If a symbol is superscripted of external perturbations. We set the goal that our approach is efﬁcient and works stand-alone with noisy and unﬁltered on the left side by i, this means that the referring vector or tensor is expressed as its 3 × 1 (vector) or 3 × 3 (tensor) inertia sensor signals. Other inaccuracies or disturbances in the evaluation of the EOMs are, e.g., measurement errors of coordinate matrix with respect to frame i. This notation is adopted from Woernle (2016). Each summand of the right- encoder inadequacies in the mechanical design of the robot, the underlying rigid body hypothesis, or modeling errors of hand side of Eq. 14 is a 3 × 1 coordinate matrix. The number of n rigid bodies corresponds to the number of summands. the inertia data. However, these types of disturbance are not part of this work. By sensor fusion with 6-axis force/torque sensors and combined with an adapted Kalman ﬁlter, our 4.2 Analytically derived EOMs contribution is extendable to become less sensitive to the mentioned inaccuracies in the model assumption. We link the arbitrary moving reference point and the origin of the NIF to a material point of the robot. This point coincides 4.1 EOMs for a multi-body system with the origin of the root segment, called pelvis. As seen in Fig. 2,the x-axis of the NIF is deﬁned as pointing always To describe the dynamics of a multi-body system, we sum parallel to the sagittal plane of the robot and the z-axis points up the Eqs. 7 and 11 for all n rigid bodies. In this sum, the parallel to the normal vector of the ground surface. Also, it reaction forces and moments between two body segments i is feasible to place the NIF on the ground surface. E.g., the and i + 1 cancel out, due to Newton’s third law. The external origin could coincide with the moving immaterial point of 123 Autonomous Robots Due to the addition formula of angular velocities, for Eq. 15 follows dω dω dω i i0 i (i −1)0 i i (i −1) = + dt dt dt (16) dω i (i −1)0 := +¨ q n. dt Here, q ¨ n is the joint acceleration vector. The vector n describes the unit vector of the joint axis. In general, ω ˙ (i −1)0 can be calculated by dω i (i −1)0 ω ˙ = + ω × ω , (17) (i −1)0 i0 (i −1)0 dt using Eq. 4. Equation 17 can be rearranged accordingly dω i (i −1)0 = ω ˙ − ω × ω (i −1)0 i0 (i −1)0 dt = ω ˙ (i −1)0 (18) Fig. 2 Visualization of the vectors and frames WF and NIF in the − (ω + ω ) × ω (i −1)0 i (i −1) (i −1)0 mentioned application related to the robot Sweaty. The vector r pelvis describes the absolute position of the root segment with respect to the = ω ˙ −˙ q n × ω . (i −1)0 i (i −1)0 origin of the WF. r and r indicate the relative position vectors foot foot R L from the origin of the NIF to origin of each foot Here, q ˙ is the joint velocity of body i compared to i − 1. For ω ˙ ﬁnally follows i0 ω ˙ = ω ˙ + ω ×˙ q n +¨ q n. (19) the vertical projection of the origin of the BFF from the root i0 (i −1)0 (i −1)0 i i segment onto the ground surface. (0) i (0) ˙ The matrix ω ˙ is given by, i0 On the other hand, the approaches based on L = M evaluate the EOMs with respect to a ﬁxed spatial point. Fig- i i i i i ω ˙ = ω ˙ + ω ×˙ q n +¨ q n. (20) ure 2 displays the difference between the mentioned (Henze i0 (i −1)0 (i −1)0 i i et al., 2016; Buschmann, 2011; Kajita et al., 2014) and our 4.4 IZMP calculation approach. These approaches use an absolute position vec- tor pointing from the ﬁxed reference point to at least one To calculate the IZMP, the EOMs will be evaluated with mea- material point of the robot, e.g., r . With kinematics, it pelvis is then feasible to calculate the absolute position vectors of sured sensor and actuator information at discrete points in all segments with respect to the ﬁxed reference point. Since time. Due to the time derivatives being available as measured values, we can interpret the coupled nonlinear differential our chosen reference point coincides with one material point, no absolute position vector is needed. Thus we can use rela- equations as algebraic equations. We formulate all equations presented in Sect. 3 as continuous-time equations. They can tive vectors, e.g., r and r pointing from the reference foot foot R L point to the origin of each foot. be evaluated without adjustments at any discrete instant in time t . For the calculation of the IZMP, we use Eq. 13. We assume the ground reaction force F and the gravity force mg act gr 4.3 Relative derivative of angular velocity externally on the robot. Then the following applies to the sum (P) of all applied external moments M , ext We use the relative derivative to link the joint kinematics with the total angular velocity ω and angular acceleration i0 (P) (IZMP) M = M + r × F + r × mg (21) IZMP gr COM ext ω ˙ . Using the law of time derivative (Eq. 4), the following i0 relationship generally applies (IZMP) Here, M is the resultant vector of the moment at point r . Per the deﬁnition of the IZMP, the components of IZMP (IZMP) M are zero for the horizontal plane. When we substi- dω dω i i0 i i0 ω ˙ = + ω × ω = . (15) i0 i0 i0 tute Eq. 5 into 21, the ground reaction force can be replaced. dt dt 123 Autonomous Robots Further, substituting this into Eq. 13, the IZMP calculation equation is given by (IZMP) M + r × (P − mg) + r × mg IZMP COM = m (r × a + r × (ω ˙ × s ) i H O H i0 i i i i i =1 (22) + r × (ω × (ω × s )) + s × a ) H i0 i0 i i O i i (O ) (O ) i i Fig. 3 Location of the IMU and the accelerometers on a rigid body to + ω ˙ + ω × ω . i0 i0 i0 i i calculate the angular acceleration. r describes the position vector of the IMU. r are feasible relative position vectors of the accelerometers with respect to the IMU. r are linearly independent to each other This system of the coordinate equations from Eq. 22 is under- determined due to the skew-symmetric property of the cross product. To obtain a single unique solution of the equation behavior of the angular velocity signal and on the increment system, we calculate the IZMP coordinates upon the assump- t used. tion that the IZMP is located in the ground surface plane. For our approach, we propose an algebraic method to cal- Therefore, the z-coordinate of the IZMP is set equal to that culate the angular acceleration vector, similar to the proposed of the reference point in relation to the ground surface plane. methods, e.g., from Padgaonkar et al. (1975). Therefore, we place three additional accelerometers on the root segment and 4.4.1 Required inputs use the rigid body assumption. Figure 3 shows the schematic setup. With this setup, we measure the linear acceleration By analyzing Eq. 13 and regarding the underlying ﬂoating- vector of speciﬁc spatial points and the angular velocity vec- base model, the required input signals can be identiﬁed. From tor of the rigid body. Eq. 13 it follows the evaluation of the dynamics requires the In general, the locations of the sensors on the rigid body position, angular velocity, linear and angular acceleration must satisfy the condition that the connection vectors from vector. For all segments of the multi-body system except the the IMU to the accelerometers are linearly independent. This root segment, these vectors can be calculated through sensor prevents the sensors from lying in a spatial plane. If the sen- measurements of the actuators and the total angular velocities sors are placed with orthogonal connection vectors to each ω and ω ˙ of the previous segment known from the other, then uniaxial accelerometers can be used. (i −1)0 (i −1)0 recursive loop. As shown, these vectors are functions of the For the acceleration of an arbitrary ﬁxed point on a rigid joint angle q , joint velocity q ˙ , and joint acceleration q ¨ .The i i i body applies in general, commonly used actuator encoder in robotics delivers signals for q . We calculate q ˙ and q ¨ with numerical derivatives using i i i a = a + ω ˙ × r + ω × (ω × r ). (23) i 0 i i the actuator controller. Therefore, we can use the increment t corresponding to the frequency of the actuator controller. Eq. 23 can be rearranged as follows, Common frequencies of actuator controllers used in robotics are in the range of kilohertz. Therefore, we name this pro- ω ˙ × r = a − a − ω × (ω × r ). (24) i i 0 i cedure actuator coordinate differentiation (ACD). Due to the option of choosing an arbitrary moving reference point in To separate the term ω ˙ , the cyclic permutation of the vector a non-inertial frame and the analytical way of deriving the triple product is used. For this purpose, the scalar product of EOMs, we name the proposed method MR-AD. Eq. 24 and the connection vector r is formed, as follows Considering the root segment, these vectors must be pro- vided from sensor data or estimations. Since our chosen (ω × r ) · r = (a − a − ω × (ω × r )) · r . (25) i j i 0 i j reference point coincides with a material point of the root segment, we can use a relative body-ﬁxed position vector. Applying the cyclic permutation Eq. 25 gives The angular velocity is measured by a gyroscope and the linear acceleration vector by an accelerometer. (r × r ) · ω ˙ = (a − a − ω × (ω × r )) · r . (26) In perspective of the authors, direct angular accelera- i j i 0 i j tion measurements are not commonly used in robotics. One approach is to calculate the angular acceleration with a Eq. 26 leads to three scalar equations. Each equation contains numerical derivative of ω . In general, numerical derivatives one coefﬁcient from the angular acceleration vector ω ˙.For are noise-amplifying. So the accuracy depends on the noise i = 2, j = 3; i = 3, j = 1 and i = 1, j = 2 the following 123 Autonomous Robots applies, ω ˙ = ((a − a − ω × (ω × r )) · r ) (27) x 2 0 2 3 r r 2 3 a a 2z 0z = − − ω ω y z r r 2 2 ω ˙ = ((a − a − ω × (ω × r )) · r ) (28) y 3 0 3 1 r r 3 1 a a 3x 0x = − − ω ω z x r r 3 3 ω ˙ = ((a − a − ω × (ω × r )) · r ) (29) z 1 0 1 2 r r 1 2 a a 1y 0y = − − ω ω . x y r r 1 1 4.4.2 Momentum based dynamics calculation To put the generated results in perspective, the momentum- Fig. 4 The real humanoid robot Sweaty (left) and the simulation model based differentiation method (IR-MD) for calculating the (right), which was designed to interact autonomously with the human- linear and angular momentum balance proposed by Kajita et made habitat al. is also implemented. The implementation strictly follows the procedure given in Kajita et al. (2014). Due to the numeri- cal derivative calculation, only the equations for momentum CAN-BUS is clocked with 125 Hz. Sweaty was designed to and angular momentum are used. To determine the rate of interact autonomously with the human-made habitat. Each change of these quantities in time, a numerical time derivative sensor used has an equivalent to the human senses. Sweaty has based on the differential quotient is used. Since the calculated no sensors like GPS or linear velocity measurement installed. momentum and angular momentum vectors are subtracted The global position and linear velocity must be estimated. from each other in two different instants in time, they must Figure 4 (right) shows the simulation model of Sweaty. be referenced in the time-constant WF. The purpose of IR- This model has 31 DOF. The DOFs of both hands are not MD is to reduce the complexity of the required equations. implemented. Webots 2021b from the company Cyberbotics Also, any dependencies on the acceleration states of the seg- Ltd. is used as the simulation tool. In the simulation, however, ments and actuators are omitted. Kajita et al. use the ﬂoating all DOFs are simulated as rotational motors. This modiﬁ- base model for the mentioned algorithm. Consequently, the cation is made due to the closed kinematic chains of the MD method requires q and q ˙ and for the root segment the i i rod-and-lever systems leading to a decrease in simulation absolute position, the absolute linear and the angular veloc- performance. ity vector must also be speciﬁed. The increment t of the Sweaty’s software architecture is based on ROS2, where numerical differentiation is linked to the frequency of the information is exchanged in a network according to the pub- ﬁeld bus in physical robot systems. lisher and subscriber principle. Figure 5 shows the ROS2 nodes and their structural connections used to generate the 4.5 Physical robot Sweaty and simulation model results for this work. The interface between ROS2 and the simulator is the node Sweaty, shown in Fig. 4 (left), is 1.65 m tall and weighs 27.6 Webots-Controller-Client. Via this interface, the node Motion kg. The humanoid robot has 39 degrees of freedom (DOF). sends the desired motion pattern to the simulator. In the Spindle motors actuate the DOFs of the lower body and torso. simulator, these data are processed and the physics engine Rod-and-lever systems translate the linear motion of the simulates the dynamics for a ﬁxed time step. The simula- spindle nut into rotational motion. The mechanical structure tion results for joint angles, sensor data, and the current time of the rod-and-lever system leads to a variable gear ratio, are made available by the simulator to the ROS2 network which results in variable torque and angular velocity char- via the node Webots-Controller-Client. Webots 2021b does acteristics depending on the joint angle. This design also not provide information about the angular velocity q ˙ of the minimizes the joint backlash due to the high gear ratio of actuators. Therefore, the required input signals, e.g., angular the spindles. The Xsens MTi-20 VRU IMU is mounted on velocities q ˙ or the angular accelerations q ¨ of the actuators i i the hip of Sweaty. Two CAN-BUS lines provide communi- are calculated through a differential quotient with the simu- cation between actuator controllers and the control unit. The lation time step t.The TF-Framework of ROS2 makes it 123 Autonomous Robots Fig. 5 Nodegraph of the setup to generate the simulation results possible to request arbitrary coordinate transformations by Table 1 List of the main abbreviation using the joint angle information and other sensor data. Abbreviation Meaning Especially for this work, a GPS and an IMU are located ACD Actuator controller differentiation in the root segment. The IMU combines a gyroscope and an BFF Body-ﬁxed frame accelerometer. It delivers the signals for angular velocity, lin- IR-MD Inertial reference—momentum derivative ear acceleration, and the orientation of the root segment. The MR-AD Moving reference—analytically derived GPS signal is subscribed by the node Localizer and published NIF Non-inertial frame as the transformation from WF to the root segment in the TF- WF Inertial frame Framework. Furthermore, this node takes a time derivative of the GPS signal by using the differential quotient to provide a Method to calculate the time derivatives of the actuator position signal the linear velocity vector. Both mentioned calculation meth- q obtaining q ˙ and q ¨ i i i Dynamics evaluated by numerical differentiation of momentum based ods are integrated parallel as separate nodes. This ensures equations for a ﬁxed reference point in inertial frame compatibility, as the calculations access the same input sig- Dynamics evaluated by analytically derived EOMs for an arbitrary nals in each time step. Both calculation methods publish their moving reference point in non-inertial frame results, which are subscribed by the node IZMP Calculation. The calculated coordinates of the IZMP are transformed into choose the BFF of the root segment as the common evalua- a common evaluation frame via the TF-Framework and are tion frame, which corresponds to the introduced NIF (see recorded in ROS2 bag ﬁles. Sect. 4.2). A further quantitative evaluation is the relative | f (t )− f (t )| n proc k ref k 1 error . It is deﬁned as = ·100%· . k=0 | f (t )| n ref k 4.6 Evaluation methods Here, f (t ) represents the value at the time point t of the proc k k investigated procedure’s signal. Similarly, f (t ) represents ref k We verify our MR-AD approach, compare it to IR-MD and the reference signal for the comparison at the time point t . show its general scopes using simulation results. Further, To analyze the applicability to physical robots in real- we test our approach on the physical robot Sweaty. Table 1 world applications, we add white noise behavior to the ideal reviews all of the necessary main abbreviations for the fol- sensor signals to simulate typical noisy signals s ˆ(t ) from lowing sections. physical sensors. White noise behavior is added to linear Due to Sweaty having no GPS or sensors for measuring velocity, angular velocity, and linear acceleration signals. It the linear velocity installed, the global position and the linear generally applies s ˆ(t ) = s(t ) + n(t ). Here, s(t ) is the original velocity vector are not available during the execution of the ideal signal and n(t ) is the noise. n(t ) is generated normally experiment. An estimation by our localization and odome- distributed depending on the standard deviation σ and the try does not provide usable results. Therefore, we focus on mean μ of the physical sensor. We use σ = ν · f to calculate MR-AD for the evaluation of the EOMs on the real robot. the standard deviation for the gyroscope and accelerometer For the comparison in simulation, the mentioned approaches as mentioned in the data sheet of Xsens MTi-20 VRU. ν is the perform different test cases and the obtained results will be noise density and f is the sample rate. The mean μ is set to compared quantitatively through the IZMP coordinates. We zero. 123 Autonomous Robots Sweaty does not have a sensor to measure linear veloc- Table 2 Parameters and properties of the line diagrams ity. Therefore, we use the same standard deviation σ for Figure ACD/bus (Hz) Noise the noise term of the linear velocity signal that applies to Fig. 6 125 None the accelerometer. The following equation applies to cal- Fig. 7 62.5, 125, 200, 250 None culate the noise for each coordinate direction, s ˆ (t ) = β k a b Fig. 8 125 IMU ,GPS s (t ) + n (t ) = s (t ) + g (σ, μ). Here, t indicates the β k β k β k k discrete point in time. Fig. 12 1000/125 Real sensors For this work, we add noise for our approach to a for i Fig. 11 125 None in range i = 1, 2, 3, 4 and ω of the root segment. It applies If the column ACD/Bus contains one value, it belongs to both ACD and a =ˆs and ω =ˆs . Here, β indicates the coordi- i β aβ β ωβ ﬁeld bus frequency Noisy accelerometer and gyroscope signals nate axis in rage β = x , y, z. For IR-MD, we add noise Noisy velocity signal from Webots just to the velocity v of the root segment, in the same way, it applies to a . Thus for IR-MD, we expect the inﬂuence Table 3 Parameters and properties of the bar diagrams of noise as underestimated since velocity measurements are usually indirect measurements based on ﬁltered acceleration Setup Method ACD/bus (Hz) Ang. acc. calc and camera position data. A IR-MD 125 Diff. quotient In general, real robot systems offer the possibility to pro- B MR-AD 125 Algebraic vide higher update frequencies of sensor and actuator signals C IR-MD 500/125 Diff. quotient compared to the ﬁeld bus frequency. To analyze the impacts D MR-AD 500/125 Algebraic of this characteristic in calculating the dynamics, the sim- E MR-AD 500/125 Diff. quotient ulation runs with a time step of 0.002 s (500 Hz). However, the dynamics calculation runs with 0.008 s (125 Hz). Thus all If the column ACD/Bus contains only one value, it belongs to both ACD and ﬁeld bus frequency sensor and actuator signals are available in 500 Hz frequency One IMU and three accelerometers are mounted on the root segment for processing. Consequently, the frequency of the simulation to measure ω ˙ set to 500 Hz can be interpreted as the actuator controller fre- b Differential quotient is used quency, and the dynamics calculation frequency set to 125 Hz as the ﬁeld bus frequency of a physical robot. 5.1 Initial conditions for the simulation 5 Results At ﬁrst, we compare the results of MR-AD and IR-MD under ideal conditions. For this purpose, the robot walks in an arbi- In the following sections, we show the results of the EOM trary direction starting from the origin of the WF at a gait evaluation by the IZMP. Each line diagram listed in Table 2 speed of 0.3m/s. Optimal signal qualities from the simula- contains curve developments for the MR-AD and IR-MD tion are provided for each calculation method. Figure 6 shows method. The adjustable parameters are the ACD and ﬁeld exemplary curves of the x- and y-coordinate from the IZMP. bus frequencies and the signal noise that can be switched on The gait phase of the robot can be divided into two sub-step or off. If signal noise is present, it is applied, as explained in phases as shown in Fig. 6. The jump in both curves describes Sect. 4.6. On the one hand, the column ACD/Bus of Table 2 the switch of the support foot and thus the transition from contains the chosen ACD and ﬁeld bus frequencies for the sub-step one into two. From the data, it follows that a sub- simulated or real robot. On the other hand in the column step takes 0.448 s, and the entire step takes 0.896 s. In the Noise speciﬁes the present signal noise for each test case. development of the x-coordinate, the two sub-steps can be Likewise, the adjustable parameters for all bar diagrams identiﬁed via the consecutively similar curve sections. The are the ACD and ﬁeld bus frequencies and the angular accel- sub-step phase begins with the heel strike and ends with the eration vector determination approach. Table 3 contains the toe-off. The x-coordinate starts after the peaks of the heel chosen parameters and properties of each bar diagram. It strike at approx. − 0.08 m and increases in the development displays in the column ACD/Bus the ACD and ﬁeld bus fre- to approx. 0.1 m until the toe-off event. The x-coordinate, quencies used. In the column Ang. acc. calc. it speciﬁes the related to the foot, moves from the rear area to the front area underlying method of angular acceleration vector estimation during the duration of a sub-step. for each setup. Further, the signal noise is turned on for each The y-coordinate starts with negative values. The develop- setup (see Sect. 4.6). ment of the curve shows a relatively constant range between If the column ACD/Field bus of Tables 2 and 3 contain the time frame 0.1 s and 0.25 s. In the further development only one value, this value represents the frequency of both of the sub-step, the value of the y-coordinate moves towards ACD and ﬁeld bus. the zero line. Here, the robot tilts towards the swing leg. 123 Autonomous Robots Fig. 6 Comparison of the 0.2 0.1 calculated IZMP x-(left)and y-coordinate (right) of MR-AD 0.1 and IR-MD under ideal conditions in simulation. The calculations are executed in 0 0 parallel for the same step phase. The gait speed is0.3 m/s −0.1 MR-AD MR-AD −0.1 IR-MD IR-MD −0.2 00.10.20.30.40.50.60.70.8 00.10.20.30.40.50.60.70.8 Time t [s] Time t [s] 5.2 Simulated robot applications The output signals of physical sensors are generally noisy. To simulate the qualitative impact of this noise, noise behavior is added as explained in 4.6. Figure 8 shows the white noise effects of linear velocity, angular velocity, and linear accel- eration signals on the x- and y-coordinates of the IZMP. For Rel. error MR-AD comparison, the coordinate curves calculated by MR-AD and Rel. error IR-MD IR-MD with ideal input signals are also shown in Fig 8.All 45 8 16 displayed curves are generated with a simulation and dynam- Simulation time step Δt [ms] ics calculation frequency of 125 Hz. Figure 8 shows MR-AD Fig. 7 Convergence behavior for the IZMP y-coordinates of MR-AD and IR-MD do not provide similar curves compared to Fig. 6. and IR-MD for decreasing simulation time steps and ideal signals. The The noise behavior of the sensors has a marginal impact on reference curve for the corresponding method is generated with a time the MR-AD method compared to the curves with ideal sensor step of 0.002 s. The gait speed of 0.3 m/s are the same for each case signals. In contrast, noise behavior has a signiﬁcant impact on the IR-MD method. In contrast to Fig. 8,Fig. 9 shows the relative errors of the IZMP y-coordinate compared between MR-AD and IR- MD. Setup B and D use MR-AD and Setup A and C use The development of the y-coordinate curve indicates that the IR-MD. Further, Setup A and B use a frequency of 125 Hz analyzed step starts with the right leg as the support leg. for ACD and ﬁeld bus. In contrast, Setup C and D use sensor In Fig. 7 the convergence behavior of the IR-MD and MR- and actuator controller frequencies of 500 Hz and a modeled AD methods for the y-coordinate of IZMP are compared at ﬁeld bus frequency of 125 Hz. For IR-MD increased sensor different simulation time steps. The tested frequencies are and actuator controller frequencies impact the accuracy neg- 250 Hz, 200 Hz, 125 Hz and 62.5 Hz. All curves were gen- atively. The relative error increases from 65.2 Hz to 112.4%. erated with ideal signals. The reference curves have been For MR-AD the contrary is observed. Here, increased sen- determined with the corresponding method at a frequency sor and actuator controller frequencies increase the accuracy. of 500 Hz. Figure 7 shows a linear convergence behavior for The relative error decreases from 16.1 to 4.4%. both methods with respect to the simulation time step t. Fig. 8 Comparison between the 0.3 0.3 ideal MR-AD results of the IZMP x-and noised MR-AD 0.2 y-coordinate of MR-AD and 0.2 ideal IR-MD IR-MD with noised linear noised IR-MD 0.1 velocity, angular velocity and 0.1 linear acceleration signals at a simulation time step of 0.008 s −0.1 (corresponds to the ﬁeld bus ideal MR-AD −0.2 noised MR-AD frequency of the real robot −0.1 ideal IR-MD Sweaty). The gait speed is −0.3 noised IR-MD −0.2 0.3 m/s 00.20.40.60.800.20.40.60.8 Time t [s] Time t [s] Rel. error y -coord. [%] x -coordinate [m] x -coordinate [m] y -coordinate [m] y -coordinate [m] Autonomous Robots 120 112.42 0.2 0.1 65.22 Level ground Slope 10 16.05 4.44 −0.1 Setup A Setup B Setup C Setup D Fig. 9 Relative error comparison between MR-AD (Setup B and D) −0.2 and IR-MD (Setup A and C). Setup A and B use a frequency of 125 Hz for ACD and ﬁeld bus. Setup C and D use 500 Hz for ACD and 125 Hz −0.10 0.1 for ﬁeld bus. Setup B and D use the method for measuring of angular x -coordinate [m] acceleration vector of the root segment. The relative error is calculated compared to the corresponding method with ideal signals and a fre- Fig. 11 The IZMP x-and y-coordinate calculated by MR-AD for the quency of 500 Hz for ACD and ﬁeld bus xy-plane and referenced in the evaluation frame while walking on a level surface and slope. The gait speed is 0.15 m/s. The slope angle is set to 10 . The support polygon for each foot is represented with the shown 18.9 polygons 16.05 Next, we analyze the impact of walking onto a slope, espe- cially for the IZMP calculation. The robot starts on level ground (slope angle 0 ) and walks up an inclined plane with 4.44 a ﬁxed slope angle of 10 without stopping. For this purpose, the gait speed is 0.15 m/s. Figure 11 shows the IZMP coordinates referenced in the Setup B Setup E Setup D evaluation frame and displayed in the xy-plane of the ground Fig. 10 Comparison between the different calculation methods. Setup surface for one exemplary gait cycle. Here we use the time B, D and E use MR-AD. Setup B uses a frequency of 125 Hz for ACD and t as the parameter for the xy-plot. In Fig. 11 the solid curve ﬁeld bus. Setup D and E use 500 Hz for ACD and 125 Hz for ﬁeld bus. indicates a gait cycle on the level surface, the dotted curve Setup E uses a differential quotient to calculate the angular acceleration vector of the root segment. Setup B and D use the method for measuring indicates a gait cycle walking on a slope (slope angle 10 ). of angular acceleration vector of the root segment. The relative error The polygons in Fig. 11 represent the support polygon of is calculated compared to the corresponding method with ideal signals each foot referenced in the evaluation frame. Both support and a frequency of 500 Hz for ACD and ﬁeld bus polygons are calculated with the kinematics in the respective stance phase. Due to the support polygons being repre- sented in the NIF, both support polygons have the same Figure 10 compares the impact of the estimation of the angular acceleration on the IZMP y-coordinate. To analyze x-coordinates. For the gait cycle on the level surface as well as for the gait cycle on the slope, it can be seen that IZMP just the impact of the method providing the angular acceler- curves are located inside the support polygon. ation vector of the root segment, Setup B, D and E use the MR-AD. Setup E and B use 500 Hz for ACD and 125 Hz for the ﬁeld bus frequency. Setup E uses a differential quo- 5.3 Real robot application tient of the measured angular velocity and Setup B and D use the angular acceleration measurement method to deter- Figure 12 shows the development of the x- and y-coordinates mine this vector of the root segment. The relative error is of the IZMP calculate on the physical robot Sweaty. The actu- calculated compared to the corresponding method with the ator controllers return the actuator positions with a frequency same angular acceleration calculation method, ideal input of 1 kHz and the CAN-BUS frequency is set to 125 Hz. A signals and a frequency of 500 Hz for ACD and ﬁeld bus. differential quotient calculates the angular acceleration of The numerical differentiation leads to a higher relative error the root segment. The robot walks on ﬂat ground with a gait (18.9%) compared to the method using angular acceleration speed of 0.15 m/s. In comparison to the developments of the measurement (4.4%) under equal test case settings. Setup E x- and y-coordinates of MR-AD shown in Fig. 8, these curves has also a higher relative error compared to Setup B using have a smoother transition in the change of the stance foot. the ACD frequency of 125 Hz and the angular acceleration The x-coordinate curve shows a larger slope followed by measurement (16.1%), shown in Fig. 9. a relatively constant proﬁle at about 0.1 m, compared with Rel. error for y -coord. [%] Rel. error for y -coord. [%] y -coordinate [m] Autonomous Robots Fig. 12 IZMP x-and 0.2 y-coordinate of MR-AD 0.1 calculated on the physical robot 0.1 Sweaty. The gait speed is 0.15 m/s −0.1 −0.1 −0.2 −0.2 MR-AD real robot MR-AD real robot 00.20.40.60.800.20.40.60.8 Time t [s] Time t [s] Fig. 8. Overall, the y-coordinate curve shows a similar but differences in calculating time derivatives cause the devia- more inconstant development as seen in Fig. 8. tions in the development of IZMP coordinate curves shown in Fig. 6). As far as we consider realistic noise behavior of physical 6 Discussion sensor signals, e.g., for the linear velocity or the position vec- tor, the results of the MD method are no longer usable (see Before we discuss our main ﬁndings, we start with the inter- Figs. 8, 9). This is caused by the dependency on the quantities pretation of the simulation characteristics. The convergence from different instances in time. Due to this characteristic, analysis (Fig. 7) shows for AD and MD methods linear con- the white noise behavior of real sensors has a direct impact vergence with respect to the simulation time step t.The on the results. By using discrete sensor signal values, the simulator Webots 2021b uses a ﬁxed time step for the physics dynamics calculation inevitably includes at least two differ- engine. This time step is speciﬁed in milliseconds and can ent noise terms. Due to the randomness of the underlying be set to an integer. Therefore, simulations are limited at normal distribution of the noise, the numerical time deriva- t = 0.001 s (1000 Hz). With this time step, the simulation tives could determine physically unrealistic gradients of the delivers no usable results simulating the robot Sweaty due momentum. Figures 8 and 9 exactly show this behavior. to impaired communication between the Webots-Controller- Further, measuring linear velocity is impractical since a Client and the other ROS2 nodes. Regarding this, we set direct measurement method is not feasible. Thus, the analyt- simulation results generated with a time step of 0.002 s as ically obtained derivatives have an obvious advantage as far the reference. However, due to the linear convergence, we as noisy signals are concerned. Consequently, they are signif- consider it appropriate to deﬁne these curves as reference icantly better suited for calculating the dynamics of physical curves. Further, for the theoretical comparisons of all simu- robotic systems. lation results, we omit measurement inaccuracies of encoders Second, for physical robots, it is commonly feasible to and joint backlash as signal noise for q to exclusively analyze operate sensors and actuator controllers with signiﬁcantly the inﬂuence of white noise of inertial sensors and velocity higher sample rates than the ﬁeld bus. A numerical deriva- estimations. and joint acceleration tive can calculate the joint velocity q ˙ Our ﬁrst main ﬁndings are about the characteristics of the q ¨ using the joint position q provided by the encoder. Fig- i i MD method. Generally, this approach reduces the complexity ure 9 shows the impact of higher sample rates. Due to the of the EOMs since only total linear and angular momentum smaller time step, both methods beneﬁt theoretically by the have to be calculated. Consequently, IR-MD gets rid of joint increased accuracy of the calculated q ˙ and q ¨ (see. Figure 7). i i acceleration q ¨ , linear and angular acceleration vectors of Additionally, the measurement error of incremental encoders the root segment. However, it requires the linear velocity does not scale proportionally with the square root of the sam- vector of the root segment. Likewise, taking numerical time pling frequency compared to, e.g., inertial sensors used (see derivatives, the calculation requires quantities of at least two Sect. 4.6). discrete points in time. For MR-AD, the calculation requires However, the numerical time derivative of IR-MD cal- quantities of one instant in time and no linear velocity vectors. culates even larger gradients for the momentum in terms Instead, the complexity of the calculation increases. These of magnitude caused by the combination of the noise- amplifying numerical differentiation and the white noise in the sensor signals. Figure 9 shows ﬁnally that the inﬂuence of the larger inaccuracies due to the white noise in the sensor signals in combination with the numerical time derivatives x -coordinate [m] y -coordinate [m] Autonomous Robots dominates. The relative error increases signiﬁcantly com- between the accelerometers. For example, if the distance r pared to the ideal reference. In contrast, MR-AD beneﬁts tends to zero, the measured values a and a converge. The 2z 0z through higher sample rates for sensors and actuator con- noise behavior is then the dominant difference between these trollers due to increased accuracy of q ˙ and q ¨ . It is further measurements. Thus ω ˙ would be mostly inﬂuenced by the i i x signiﬁcantly less impacted through white noise than the MD noise. method, as seen in Fig. 8. Here, we want to discuss the use of three additional BasedontheresultsshowninFig. 9, we propose to execute accelerometers to determine all joint velocities q ˙ and accel- this calculation directly in the actuator controllers. In general, eration q ¨ of the robot. The above mentioned arguments the frequency of the ﬁeld bus from robots depends mostly against numerical differentiation apply analogously. One on the total number of motors and sensors connected to the consequence would be to omit numerical differentiation ﬁeld bus. Typically frequencies of ﬁeld buses (e.g., CAN or entirely when calculating joint kinematics. Further, it is thus EtherCAT) for robots with more than 30 degrees of freedom independent of the actuator controller and ﬁeld bus frequen- are in the range of 100–1000 Hz. cies. A disadvantage is the increasing number of installed For physical robots, the higher frequency of the actuator sensors and the transmitted data. controllers is an advantage regarding the accuracy of q ˙ and Fourth, MR-AD and all analytical derivative methods are q ¨ . The errors of the numerical derivative approximation of based on the same basic physical laws. Provided that consis- q ˙ and q ¨ decrease due to smaller t. However, numerical tent inputs are given for the calculation, the results must be i i differentiation with t tending to zero generally ampliﬁes identical. Compared to our MR-AD, the mentioned analyt- signal noise. Therefore, e.g., implemented ﬁltering methods ical derivative methods require the absolute position vector in the actuator controller have to compensate for the noise- of the root segment from a spatially ﬁxed reference point as amplifying effects. Through the signiﬁcantly higher actuator an input. controller frequency, we can design ﬁltering methods with Henze et al. (2016) estimate the absolute position vector adapted delays so that they have no effects on the dynamics of the root segment by sensor fusion. Therefore, this kind calculations executed in ﬁeld bus frequency. Consequently, of estimation requires the kinematics as well as the absolute in this case we do not use the t of the ﬁeld bus frequency position vector and orientation of each end effector in contact. in any calculation step. Therefore, it has no impact on the One assumption required for the proposed sensor fusion is a accuracy of the dynamics calculation. rigid contact of at least one body segment. Third, in contrast to the MD method, analytical deriva- Despite the redundant calculation and the introduced tive EOMs generally require the angular acceleration vector averaging method in the approach of Henze et al. (2016), of the root segment. In this work, we compare two different this sensor fusion always contains inaccuracies due to non- approaches to assess angular acceleration. On the one hand, a present rigid contact and white noise of the IMU sensor numerical derivative of the angular velocity signal, measured signals. Also, the continuous determination of the absolute by IMU, is used. On the other hand, we present a measure- position vectors of all end effectors contains inaccuracies, ment through three additional uniaxial accelerometers. especially on unknown, uneven terrain and through undeﬁned For calculating the angular acceleration of the root seg- contacts. This kind of estimation for the absolute position ment, a noise-amplifying numerical differentiation has the vector of the root segment is, in our view, error-prone. same disadvantages as mentioned before. The increment t The approach of Buschmann (2011) avoids numerical used depends on the sampling rate of the sensor. Noise in inaccuracies through large distances from the initial reference the estimated angular velocity signal combined with high point to the robot by shifting the reference frame depending sampling rates could result in unrealistic angular accelera- on the gait cycle. In general, this reference frame does not tion approximations. Figure 10 conﬁrms this. Low sampling coincide with the BFF of the stance foot and the origin is rates would lead to inaccuracy for angular acceleration due ﬁxed. This means the position vector of the COM cannot to the increasing t between two measured quantities. A be calculated only through kinematics on physical robots. decreasing t through higher sampling rates increases the Buschmann does not propose an estimation of the global noise-amplifying effects of the differential quotient used. position vector of the COM, since the work aims to simulate Since the calculation method from three additional robots. accelerometers does not involve numerical differentiation, In contrast, our approach leads to high ﬂexibility in choos- the accuracy gets independent of t or the physical sensor ing the reference point and the reference frame for dynamics update frequency. The disadvantages previously described do calculations. The reference point can be any material point, not apply. However, this includes four additional noisy sig- e.g., the origin of the root segment or any immaterial point. nals in the calculation, which affects the measurement of the The state of motion of this point compared to the other mate- angular acceleration. Equations 27, 28, and 29 show that the rial points of the robot is irrelevant. This leads to the omission components of angular acceleration depend on the distance of estimating absolute position vectors. We made this effort 123 Autonomous Robots to avoid all estimations since they can cause inaccuracies Sixth, if we compare the IZMP coordinates of the physical in the results and for saving computational capacity. E.g., robot shown in Fig. 12 with the simulation results shown in if we choose the origin of the root segments BFF as the Fig. 6, we recognize generally similar curve developments. moving reference point, we can calculate relative position The main differing characteristics of the coordinate curves vectors of arbitrary material points through kinematics with- from physical robot to simulation are, on the one hand, the out inaccuracies due to uneven ground or undeﬁned contact. more inconstant development and the noticeably smoother Consequently, conditions with no rigid body contact of at transition between the change of the stance foot. The joint least one segment or in ﬂight phases do not disturb our used backlash present and measurement errors of the encoders EOM evaluation approach. on physical robots, the elasticities, and manufacturing and Furthermore, the spatial distance between the evaluation assembly inaccuracies mainly cause these two characteris- frame and all BFFs stays always in a similar order of magni- tics. In conclusion, despite the mechanical inadequacies and tude. Similarly, the recommendation of Featherstone (2008) white noise of all sensor signals, MR-AD is suitable to calcu- for the distance of one to two radii of gyration is con- late the IZMP also for physical robots. If the proposed method tinuously satisﬁed. This reduces numerical inaccuracy. We for measuring angular acceleration vectors is applied to all further use an algorithm for calculating the relative derivative segments of the robot, we expect that the results of dynamics from numerically available values, which allows us to derive calculation on real robots will be further improved. physical quantities in moving frames. Speciﬁcally, we use it Further, Fig. 12 shows the justiﬁcation that omitting the to calculate the relative time derivative of the angular velocity measurement errors of the encoders and the joint backlash is for all body segments. With this, we evaluate the dynamics sensible in the theoretical comparisons of the applied meth- for each body segment in its BFF. This is an efﬁcient way to ods. The impact of the measurement errors of the encoders, evaluate the EOMs according to Featherstone (2008). the joint backlash, and other mechanical inaccuracies on the Fifth, a further advantage of using our approach is seen for IZMP curves is much lower than the impact of the numerical the IZMP calculation walking, e.g., on slopes. In general, one differentiation of the momentum equations with noisy linear axis of the evaluation frame must point parallel to the normal and angular velocity signals. vector of the ground surface for calculating the IZMP. This As we mentioned previously, there are disturbances in the is required to resolve the under-determination of the system dynamics calculation in general. These include inaccuracies of coordinate equations from Eq. 22 through determining the in the mechanical model compared to the real robot, e.g., due z-coordinate. This demand cannot be fulﬁlled by the com- to deviating mounting positions of the joint axes, inertia data pared approaches. Consequently, they do not provide usable or the underlying rigid body hypothesis, to name a few. To results without further adjustments in calculating the IZMP reduce the impact of inaccuracies, we propose combining a on slopes. For valid results, the compared approaches would sensor fusion from our approach with a 6-axis force/torque require a transition of the inertial frame. Likewise, the equa- sensor and adding a Kalman ﬁlter. A limitation of the sensor tion to calculate the IZMP, as, e.g., shown by Kajita et al. fusion is that it only provides valid data for a bipedal walking (2014), also has to be adapted to the rotated gravity vector. without external perturbations applied to the robot. For evaluating the EOMs in a MR, this condition for the IZMP calculation remains fulﬁlled. Here, the evaluation frame rotates in a suitable orientation to cover the slope to 7 Conclusion and future work ensure that the required axis parallelism remains fulﬁlled. IMUs placed in the feet measure the orientation of the sur- From the perspective of the authors, the combination of face during the stance phase in each gait cycle. evaluating the EOMs formulated through the Newton–Euler The comparison of the resulting IZMP coordinates for equations with respect to a moving reference point in a NIF the level ground and the slope shows that they are located and measuring the angular acceleration vector of the root within the support polygons as seen in Fig. 11. Our approach segment contributes to an efﬁcient and robust evaluation evaluates the locomotion on the slope as stable. of the inverse dynamics in robotics. The entire proposed Consequently, with our approach, the reference point is approach produces comparably similar results as the previ- adjustable to the implemented control algorithm. E.g., if we ous approaches evaluating the EOMs with respect to a ﬁxed use an IZMP-based control algorithm, the moving reference point in the inertial frame or the center of mass for a humanoid could coincide with the surface plane to be independent of robot, with the main beneﬁt that the proposed approach is less friction forces. Or, if we control the ground reaction forces, prone to noise than the previous approaches. it may be efﬁcient to place the reference point into the mov- We conﬁrm that in theory and by comparing it with the ing root segment. The calculated quantities can thus directly mentioned approaches using simulation results. In theory, trigger a reaction without additional spatial transformation. we show the conceptual advantages of our approach. E.g., 123 Autonomous Robots it gets rid on the dependency of any global position vector. tarity problems. Nonlinear Dynamics, 14(3), 231–247. https://doi. org/10.1023/A:1008292328909 Further, we show the ﬂexibility of choosing the reference Buschmann, T. (2011). Simulation and control of biped walking robots point and the evaluation frame concerning the preferences of (1 Auﬂ ed.). München: Verl. Dr. Hut. the control algorithm used. Englsberger, J. , Ott, C. , Roa, M.A. , Albu-Schaffer, A. , & Hirzinger, G. The comparison of simulation results clearly shows fur- (2011). Bipedal walking control based on Capture Point dynamics. In 2011 IEEE/RSJ international conference on intelligent robots ther advantages for the application to physical robots. On the and systems (pp. 4420–4427). CAIEEE. https://doi.org/10.1109/ one hand, we reduce the quantity of noise-amplifying numer- IROS.2011.6094435 ical differentiation. On the other hand, with our approach, Featherstone, R. (2008). Dynamics of rigid body systems. R. Feath- we can calculate the IZMP for the gait from the level to the erstone (Ed.), Rigid body dynamics algorithms (pp. 39–64). Springer. https://doi.org/10.1007/978-1-4899-7560-7_3 inclined surface. We also show that our algorithm calculates Goswami, A. (1999a). Foot rotation indicator (FRI) point: A new usable results applied to physical robots. gait planning tool to evaluate postural stability of biped robots. It is also worth mentioning that this approach is not lim- In Proceedings 1999 IEEE international conference on robotics ited to bipedal robots or humanoid robots in general. Due to and automation (Cat. No.99CH36288C) (vol. 1, pp. 47–52). USAIEEE. https://doi.org/10.1109/ROBOT.1999.769929 the generic structure, the EOMs can be calculated for multi- Goswami, A. (1999b). Postural stability of biped robots and the legged robots or human beings. The adaption to KDL and the foot-rotation indicator (FRI) point. The International Journal integration in ROS2 of this algorithm support straightforward of Robotics Research, 18(6), 523–533. https://doi.org/10.1177/ usability. 02783649922066376 Goswami, A. , & Kallem, V. (2004). Rate of change of angular Next, we will equip all segments of Sweaty with one sen- momentum and balance maintenance of biped robots. In IEEE sor unit combining an IMU and three additional uniaxial international conference on robotics and automation, 2004 pro- accelerometers to eliminate all numerical differentiation and ceedings (ICRA’04). 2004 (vol.4, pp. 3785–3790). USAIEEE. use our approach to evaluate the stability of the robot Sweaty. https://doi.org/10.1109/ROBOT.2004.1308858 Henze, B., Roa, M. A., & Ott, C. (2016). Passivity-based whole-body As a second step, we will setup up the sensor fusion with the balancing for torque-controlled humanoid robots in multi-contact self-developed 6-axis force/torque sensor and we will try to scenarios. The International Journal of Robotics Research, 35(12), implement a Kalman ﬁlter. Further, we will try to use the 1522–1543. https://doi.org/10.1177/0278364916653815 sensor unit also to evaluate the stability of human beings. Hirukawa, H. , Hattori, S. , Harada, K. , Kajita, S. , Kaneko, K. , Kane- hiro, F., & Morisawa, M. (2006). A universal stability criterion Finally, we will use this approach as the basis of a con- of the foot contact of legged robots—Adios ZMP. In Proceedings trol algorithm to stabilize both static situations and dynamic 2006 IEEE international conference on robotics and automation, locomotion. We want to align this controller speciﬁcally to 2006 (ICRA 2006) (pp. 1976–1983). USAIEEE. https://doi.org/ the control paradigms of human beings, taking beneﬁt of the 10.1109/ROBOT.2006.1641995 InJoon, M. , DongHa, Y. , MinSung, A. , & Jeakweon, H. (2020). Body identiﬁed characteristics of our approach. trajectory generation using quadratic programming in bipedal robots. In 2020 20th international conference on control, automa- Funding Open Access funding enabled and organized by Projekt tion and systems (ICCAS) (pp. 251–257). IEEE. https://doi.org/10. DEAL. The authors have no relevant ﬁnancial or non-ﬁnancial interests 23919/ICCAS50221.2020.9268204 to disclose. The project Sweaty of the University of Applied Sciences Kajita, S., Hirukawa, H., Harada, K., & Yokoi, K. (2014). Introduction Offenburg is a non-proﬁt project. The objective is to contribute to the to humanoid robotics (vol. 101). Springer. https://doi.org/10.1007/ research ﬁeld of humanoid robotics. Regional companies sponsor equip- 978-3-642-54536-8 ment. However, there is no inﬂuence on or exchange of research results. Kajita, S., Kanehiro, F., Kaneko, K., Fujiwara, K., Harada, K., Yokoi, K., & Hirukawa, H. (2003). Resolved momentum control: Humanoid Open Access This article is licensed under a Creative Commons motion planning based on the linear and angular momentum. In Attribution 4.0 International License, which permits use, sharing, adap- Proceedings 2003 IEEE/RSJ international conference on intelli- tation, distribution and reproduction in any medium or format, as gent robots and systems (IROS 2003) (cat. no.03ch37453) (vol. 2, long as you give appropriate credit to the original author(s) and the pp. 1644–1650). https://doi.org/10.1109/IROS.2003.1248880 source, provide a link to the Creative Commons licence, and indi- Kajita, S., Morisawa, M., Miura, K., Nakaoka, S., Harada, K., Kaneko, cate if changes were made. The images or other third party material K., & Yokoi, K. (2010). Biped walking stabilization based on linear in this article are included in the article’s Creative Commons licence, inverted pendulum tracking. In 2010 IEEE/RSJ international con- unless indicated otherwise in a credit line to the material. If material ference on intelligent robots and systems (pp. 4489–4496). IEEE. is not included in the article’s Creative Commons licence and your https://doi.org/10.1109/IROS.2010.5651082 intended use is not permitted by statutory regulation or exceeds the Kamioka, T., Kaneko, H., Takenaka, T., & Yoshiike, T. (2018). Simul- permitted use, you will need to obtain permission directly from the copy- taneous optimization of ZMP and footsteps based on the analytical right holder. To view a copy of this licence, visit http://creativecomm solution of divergent component of motion. In 2018 IEEE inter- ons.org/licenses/by/4.0/. national conference on robotics and automation (ICRA) (pp. 1763–1770). IEEE. https://doi.org/10.1109/ICRA.2018.8460572 Kimura, K. , Nozawa, S., Mizohana, H., Okada, K., & Inaba, M. (2018). Riding and speed governing for parallel two-wheeled scooter based References on sequential online learning control by humanoid robot. In 2018 IEEE/RSJ international conference on intelligent robots and sys- tems (IROS) (pp. 1–9). IEEE. https://doi.org/10.1109/IROS.2018. Anitescu, M., & Potra, F. A. (1997). Formulating dynamic multi-rigid- body contact problems with friction as solvable linear complemen- 8593685 123 Autonomous Robots Mehrtens, G. (1887). Technische Mechanik fester, ﬂüssiger und luftför- Wang, P., Liu, J., Zha, F., Guo, W., Wang, X., Li, M., & Sun, L. miger Körper(vol. 3). Kommissions-Verlag v. E. Toeche. (2017). A velocity estimation algorithm for legged robot. Advances Mesesan, G., Englsberger, J., & Ott, C. (2021). Online DCM trajectory in Mechanical Engineering, 9(12), 1687–8140. https://doi.org/10. adaptation for push and stumble recovery during humanoid loco- 1177/1687814017732736 motion. In 2021 IEEE international conference on robotics and Wittenburg, J. (1977). Dynamics of systems of rigid bod- automation (ICRA) (pp. 12780–12786). IEEE. https://doi.org/10. ies. Vieweg+Teubner Verlag. https://doi.org/10.1007/978-3-322- 1109/ICRA48506.2021.9560808 90942-8 Padgaonkar, A. J., Krieger, K. W., & King, A. I. (1975). Measurement Woernle, C. (2016). Mehrkörpersysteme: eine Einführung in die Kine- of angular acceleration of a rigid body using linear accelerometers. matik und Dynamik von Systemen starrer Körper (2., erweiterte Journal of Applied Mechanics, 42(3), 552–556. https://doi.org/10. Auﬂage ed.). Springer Vieweg. 1115/1.3423640 Pang, J. , & Trinkle, J. (2000). Stability characterizations of ﬁxtured Publisher’s Note Springer Nature remains neutral with regard to juris- rigid bodies with coulomb friction. In Proceedings 2000 icra. mil- dictional claims in published maps and institutional afﬁliations. lennium conference. IEEE international conference on robotics and automation (vol. 1, pp. 361–368). IEEE. https://doi.org/10. 1109/ROBOT.2000.844083 Pratt, J., Carff, J., Drakunov, S., & Goswami, A. (2006). Capture Maximilian Gießler received his Point: A step toward humanoid push recovery. In 2006 6th IEEE- B.Eng. in 2018 and M.Sc. in 2020 RAS international conference on humanoid robots (pp. 200–207). in Mechanical Engineering from University of Genova, IEEE. https://doi.org/10.1109/ICHR.2006. the Offenburg University of App- lied Sciences. He is currently a Saida, T., Yokokohji, Y., & Yoshikawa, T. (2003). FSW (fea- Ph.D. student and employed as a sible solution of wrench) for multi-legged robots. In 2003 scientiﬁc assistant at the Offen- IEEE international conference on robotics and automation (Cat. burg University of Applied Sci- No.03CH37422) (vol. 3, pp. 3815–3820). IEEE. https://doi.org/ ences. His research interests inclu- 10.1109/ROBOT.2003.1242182 de multibody dynamics and Sardain, P., & Bessonnet, G. (2004). Forces acting on a biped robot. humanoid robotics, with an Center of pressure–Zero moment point. IEEE Transactions on Sys- emphasis on bipedal locomotion tems, Man, and Cybernetics - Part A: Systems and Humans, 34(5), and stability theory. 630–637. https://doi.org/10.1109/TSMCA.2004.832811 Vijaykumar Shah, S., Kumar Saha, S., & Kumar Dutt, J. (2013). Dynam- ics of tree-type robotic systems (vol. 62). Springer. Vukobratovic, M., Borovac, B., & Šurdilovic, ´ D. (2001). Zero-moment Bernd Waltersberger has been point—Proper interpretation and new applications. In IEEE-RAS a professor at the Faculty of international conference on Humanodi robots (pp. 237–244). Mechanical and Process Engineer- IEEE. ing at the Offenburg University of Vukobratovic, M., Frank, A. A., & Juricic, D. (1970). On the stability Applied Sciences since 2010. He of biped locomotion. IEEE Transactions on Biomedical Engi- received his diploma in mechani- neeringBME, 17(1), 25–36. https://doi.org/10.1109/TBME.1970. cal engineering in 2002 from the University of Karlsruhe (now Vukobratovic, ´ M., & Stepanenko, J. (1972). On the stability of anthro- K.I.T) and was awarded his doc- pomorphic systems. Mathematical Biosciences, 15(1–2), 1–37. torate there in 2007 for a the- https://doi.org/10.1016/0025-5564(72)90061-2 sis on structural mechanics. He Wang, H. , Tian, Z., Hu, W., & Zhao, M. (2018). Human-like ZMP subsequently worked as a devel- generator and walking stabilizer based on divergent component opment engineer in the automo- of motion. In 2018 IEEE-RAS 18th international conference on tive ﬁeld and is co-inventor on 7 humanoid robots (Humanoids) (pp. 1–9). IEEE. https://doi.org/ patent applications. His research 10.1109/HUMANOIDS.2018.8624926 interests are in the dynamic interaction of machine elements and multibody dynamics, especially in the area of vibration engineering and stability theory.

Journal

Autonomous Robots – Springer Journals

Published: Apr 1, 2023

Keywords: Multi-body dynamics; Equation of motion; Moving reference point; Imaginary zero-moment point; Humanoid robot

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

Robust inverse dynamics by evaluating Newton–Euler equations with respect to a moving reference and measuring angular acceleration

References (32)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies