Open access peer-reviewed chapter

Forward and Inverse Dynamics and Quasi-Static Analysis of Mechanizes with MATLAB®

By E. Corral, J. Meneses and J.C. García-Prada

Submitted: September 29th 2015Reviewed: March 30th 2016Published: July 7th 2016

DOI: 10.5772/63372

Downloaded: 1110

Abstract

There are many potential advantages of direct and inverse dynamic and quasi-static analysis of mechanisms, namely control the risk of slippage, improve stability, better adaptation to the environment, obtaining smooth movements and optimizing energy consumption. This chapter proposes new analysis methods and algorithms to bring new solutions to the mechanics of the machines under consideration. The methodology has been developed in modular programs thanks to the flexibility of MATLAB®.

Keywords

  • mechanism
  • dynamics
  • quasi-static
  • robot
  • vehicle

1. Introduction

Nowadays, walking robots, service robotics, and unmanned ground vehicles are considered one of the main areas of research. The employment of robots for dangerous tasks makes its design a crucial point. The interest in the development of humanoid robotics (personal assistance, social task, etc.) is rising, and it is being studied by a great number of research groups. The main issue is current mobile robots are not adapted to be used in domestic environments due to their lack of maneuverability but also to their large volume and/or weight.

These days humanoid biped have a high number of actuators that are used to control the high degrees of freedom (DOF) they possess. Nonetheless, one of the biggest drawbacks in humanoids is that both the weight and the power consumption still make this technology not suitable for many tasks and/or environments. In the majority of cases, around 30% of the total weight is related to the actuators and wires, and more than 25% to the reduction systems. That is why our work is focused on finding a new dynamics analysis to simplify the design of new mechanisms and kinematic chains which, maintaining the robot functionality, does not require such a high number of actuators. This would reduce the robot mass and hence its power consumption and total cost [14].

In recent years, different research groups have developed robots based on passive walking techniques. Biped locomotion has been studied from several perspectives. Much effort is devoted to the design of optimal trajectories and stabilized walking cycles via control programs. However, researchers have not focused enough their efforts in solving the slip problem, which has been largely ignored [5, 6].

In this chapter, the kinematics and dynamics analysis of PASIBOT (a biped walking robot designed and built by the Maqlab Groups of the Universidad Carlos III de Madrid, shown in Figure 2) are presented. The methodology to do the completed study from a theoretical point of view is explained. The study objective is to calculate all the forces and torques between links, as well as the linear and angular position coordinates, velocities, and accelerations for all links, for any time. The equations have been implemented in MATLAB® code, and the corresponding results have been contrasted [712].

The program accepts the initial kinematic state and the motor torque (as a function of time) as inputs and returns the bipedal movement, including the sliding of the supporting foot. The sliding is taken into account by adding one degree of freedom. Thus, we focus on the kinematics and dynamics of the sliding supporting foot [13].

In addition, there is no doubt about all advantages unmanned vehicles have. For this reason, the kinematics and dynamics analysis is one of the principal research lines in robotics. The Tallinn University of Technology is researching in the design and development of the unmanned ground vehicle (UGV) shown in Figure 20. This UGV is an all-terrain vehicle equipped with an engine for each of its wheels. The novel aspect of this vehicle is that each wheel is attached to the body by a leg so that the angle between the latter and the body may vary thanks to the use of an attached actuator [1416].

In this chapter, a parametric quasi-static half vehicle model is implemented on MATLAB® following the explained methodology. The program calculates the variation in configuration angles that optimized desired criteria as the vehicle passes on a particular track profile. This algorithm makes it possible to find the variation of configuration angles along the track profile that keeps the applied torque constant and/or minimized and/or satisfying certain criterion.

2. Methodology

The sketch of the methodology is shown in Figure 1.

The first step is to define the mechanism. It is basically naming variables and parameters (with the correct nomenclature the implementation in MATLAB® code is made easier) and classifying the type of movement (degrees of freedom, class of joints, etc.)

Figure 1.

Sketch of the methodology.

It is possible to solve the kinematics and dynamics through the inverse dynamics or the forward dynamics algorithm depending on the type of mechanisms, the known input, and the desired output. Also, if the dynamics are complicated, the quasi-static approach can be useful. The biggest advantage of this approach is that it is easily optimized. This is a valuable tool to improve design and control with mathematics software.

In this chapter, the methodology of inverse and forward dynamics is applied to the biped walking robot PASIBOT and the quasi-static approach to an unmanned ground vehicle.

Note that the input and output of each approach are different and the equations are also different. However, the methodology can be simplified in the Figure 1 for both, and the correctly chosen nomenclature helps to generate a MATLAB® code. The developed MATLAB® code is totally parametric, giving the possibility of changing the values, adding new analysis and new degrees of freedom (like slippage), or looking for analysis of sensibility.

Another advantage of doing the whole analysis with MATLAB® code is that the result can be used as an input for redesigning or as an optimizing function.

3. Kinematics of the biped robot PASIBOT

Nowadays, certain commercial software of mechanical simulation provides the dynamics of a mechanism with a small error. But, in some cases, like a biped robot actuated by small number of actuators, it is also possible to obtain the kinematics and apply this methodology and obtain the dynamics with a low degree of error. In this chapter, the kinematics of the biped PASIBOT is developed. In next chapters, the inverse and forward dynamics are addressed using the relations obtained in the kinematics.

The biped PASIBOT is a one-degree-of-freedom mechanical system based on a combination of classical mechanisms that emulates human walking. PASIBOT is shown in Figure 2.

Figure 2.

PASIBOT.

The biped PASIBOT (see Figures 24) is a mechanism composed of three sub-mechanisms of conventional type, each of which is designed to perform a different function:

  1. The mechanism “Chebyshev” is responsible for generating a quasi-straight line.

  2. A pantograph handles extend the “Chebyshev” coupling curve.

  3. The stability of the biped is achieved by parallel extensions.

In Figure 3, the sub-mechanisms Chebyshev and pantograph as well as the trajectories for the most relevant points are shown.

The only engine of the biped conveys a full rotation to the crank of the Chebyshev mechanism, so that the end of its connecting rod (point C in Figure 3) makes a cyclical movement, one of its tracts being a quasi-straight line. This point is connected to one end of the pantograph, so that its other end (point E in Figure 3) carries an inverted and amplified from the previous path. The corresponding amplification ratio depends on the lengths of the links of the pantograph. The amplification ratio is two for the design of PASIBOT presented here.

Figure 3.

The sub-mechanisms Chebyshev and pantograph with the trajectories of their notable points.

Points A, B, and D in Figure 3 are attached to the hip, as can be seen in Figure 4. The stabilization system consists of a series of articulated parallelograms which are based on the two longest links of the pantograph. These parallelograms guarantee the parallelism between the foot in contact with the floor and the stabilizing link, the end of which slides on a slot at the hip. Since that slot is aligned with the linear segment of the Chebyshev trajectory, the supporting foot remains parallel to the slot during the whole period of support. In order to provide the opposite leg with the appropriate movement, the corresponding crank is phased out π rad (see Figure 4) in the same motor axis. In fact, both cranks are part of the same rigid element.

Figure 4.

Sub-mechanisms of PASIBOT; nomenclature and numeration for the supporting leg and angular positions for the links.

As shown in Figure 4, any link belonging to the flying leg has the same name and number as the corresponding link from the supporting leg, but with a prime to distinguish between them. Each leg comprises 12 links, apart from the engine crank (link number 8), which is shared with both legs (no link number 8’ exists), so the biped PASIBOT has a total of 24 links, including the single hip (link number 13).

Listed below are the parameters and variables that describe the kinematics and dynamics of the biped:

  1. li : length of the link i (mm)

  2. mi: mass of the link i (kg)

  3. ϴi: angle between the link i and the hip (rad)

  4. ωi: rotational velocity of link i (rad/s)

  5. αi: rotational acceleration of link i (rad/s2)

  6. Ii: inertia moment for the link i (kg mm2)

  7. rij: position vector of the ij joint from the link i center of mass (mm)

  8. rijx: x projection of the position vector (mm)

  9. rijy: y projection of the position vector (mm)

  10. fij: force exerted by the link i on the link j (N)

  11. fijx: x projection of the fij (N)

  12. fijy: y projection of the fij (N)

In Figure 5, a sequence for one step of PASIBOT is presented, as simulated with a commercial program. Note that one step corresponds to a half rotation (π rad) of the motor crank.

Figure 5.

PASIBOT gait along one step (from ϴ8= π to 3π/2 rad).

After defining the whole type of movement and the nomenclature, the kinematical study of one PASIBOT step is presented here. The kinematics is developed for the phase of “simple support,” in which the supporting foot is in contact with the horizontal ground, whereas the other leg is flying. First, no sliding between the supporting foot and the ground is considered, so it can be considered part of the ground. Hence, the biped PASIBOT is a one DOF planar mechanism, and we can refer the angular positions of any link to the angular position of the motor crank (ω8):

θi=θi(θ8),i =1,2,1,2,E1

Then, the x,y coordinates for its center of mass can be easily expressed with respect to that angle:

xi= xi(θ8);yi=yi(θ8),i =1,2,1,2,E2

The angular velocities, accelerations and the center of mass linear velocities and accelerations are obtained by taking the first and second derivatives in Eqs. (1) and (2).

Actually, the biped kinematics is divided into three closed-loop kinematic chains:

  1. 1. Chebyshev chain (links number 7, 8, 9, and 13)

The distance between motor crank and rocker fixed points (A and B in Figure 3, respectively) in a Chebyshev mechanism is lAB = 2l8, the rocker arm length is l9 = 2.5l8, the connecting rod length is l7 = 5l8, and the rocker arm and connecting rod are joined at the middle point of the latter. The link lengths have been particularized for the designed biped, and normalized to the crank length, l8 = 1. Taking into account these lengths, the Chebyshev closed-loop kinematic chain provides (see Figure 6):

2.5ejϑ72.5ejϑ9ejϑ8+2E3

Figure 6.

Chebyshev chain (lengths in units of l8).

In Eqs. (3)–(5), both projections (vertical and horizontal) for each closed-loop equation are written in a compact form following the Euler’s formula, where j is the imaginary unit.

  1. 2. Pantograph chain (links number 9, 7, 3, 6, and 13)

The tendons length is l4 = l6 = 6l8, whereas the distance between the connecting rod-femur and upper tendon-femur joints (points C and F, respectively) is lCF = 3l8, and the distance between rocker arm-hip and upper tendon-hip joints (points B and D, respectively) is lBD = 12l8, so the pantograph closed-loop kinematic chain provides (see Figure 7):

6ejϑ6+3ejϑ3+2.5(ejϑ7+ejϑ9)12j=0E4

Figure 7.

Pantograph chain (lengths in units of l8).

  1. 3. Stability chain (links number 8, 7, 10, and 13)

In our model, the stabilizing link length is l10 = 4.2l8. The vertical distance between the motor crank joint and the slot at the hip is 4l8. The horizontal projection distance between the motor crank joint and the end of the stabilizing link is called x. The stability closed-loop kinematic chain is as follows (see Figure 8):

24.ejϑ105ejϑ7+ejϑ8x+4j=0E5

Figure 8.

Stabilization chain (lengths in units of l8).

As stated below, these equations determine the angles for all the links as functions of that for the motor crank, ϑ8, which is also a function of time. Solving Eq. (3), the following expressions for the connecting rod and rocker arm angles are found:

{ϑ7=acos[4cos2ϑ8+13cosϑ810+sinϑ816cos2ϑ860cosϑ8+1002520cosϑ8]ϑ9=acos[4cos2ϑ8+13cosϑ810sinϑ816cos2ϑ860cosϑ8+10025+20cosϑ8]E6

From Eq. (4), the femur and tibia angles are found as functions of the previous ones:

{ϑ6=acos[(2.5(cosϑ7+cosϑ9)(27A)2.5(sinϑ7+sinϑ912)144A(27A)212A]ϑ3=acos[(2.5cosϑ7+cosϑ9)(27A)+2.5(sinϑ7+sinϑ9)12)36A(27A)26A]E7

where A=(2.5(cosϑ7+cosϑ9))2+(2.5(sinϑ7+sinϑ9)12)2

Finally, Eq. (5) gives the solution for the stabilizing angle:

ϑ10=asin(5sinϑ7sinϑ844.2)E8

As can be seen in Figure 3, the rest of the angles involved are identical to one of the given ones in Eqs. (6)–(8), in particular:

ϑ1=ϑ5=ϑ10ϑ12=ϑ4=ϑ3ϑ2=ϑ13=ϑ6E9

For the links belonging to the opposite leg, we apply a phase out of π radians on ϑ8:

ϑi'(ϑ8)=ϑi(ϑ8+π)E10

In order to reference all values to the ground, this corresponding base change must be applied:

ϑiground=ϑiϑ1,E11

where ϑigroundϑis the angle related to the ground system, and ϑi is the corresponding one related to the reference system fixed at link 14 (hip).

The positions of the center of mass for every link are obtained using trigonometric relations (e.g. x2=L2cosϑ2/2, y2=L2sinϑ2/2; x3=L2cosϑ2+L3cosϑ3/2, y3=L2sinϑ2+L3sinϑ3/2; etc.). Then, by time differentiating once and then twice, the angular velocity and acceleration as well as linear velocity and acceleration, respectively, for any link are calculated.

Thanks to apply this equations on MATLAB program, the kinematics of the biped robot PASIBOT can be solved for one step, considering a motor crank constant angular velocity,

ω8:ϑ8(t) =ω8·tE12

The PASIBOT possessed a single DOF, so the positions of the center of mass of link i can be referred to the angular position of the motor crank (ϑ8):

ϑi=ϑi(ϑ8),i=2,31,2,(i8)E13
Xi=Xi(ϑ8),i1E14
yi=yi(ϑ8),i1E15

4. Slipping kinematics of the biped robot PASIBOT

If the supporting foot is allowed to slip, the PASIBOT becomes a 2-DOF mechanism (the biped moves across a plane, and the supporting foot supposedly remains horizontal). Eqs. (13) and (15) remain valid, while Eq. (14) increases by the value of the supporting foot slippage x1 as follows:

xi=x1+Xi(ϑ8);i1E16

The first and second time derivatives of Eq. (16) are as follows:

x˙i=dx1dt+dXidϑ8dϑ8dt=x˙1+Xi'(ϑ8)ϑ˙8x¨i=x¨1+Xi''(ϑ8)ϑ˙82+Xi'(ϑ8)ϑ¨8E17

where a prime denotes explicit derivative with respect to ϑ8, and dots denote time derivatives.

5. Non-slipping inverse dynamics of the biped robot PASIBOT

The inputs for the dynamical problem are the kinematic magnitudes (angular acceleration, αi, and its center of mass acceleration, (aix, aiy)). The dynamical equations for the motion of every link, using Newton action-reaction law, are exposed in Eq. (18):

iFi=miaifij=fjiiTon i=Iiα}{j<ifjixk>ifikx=mix¨ij<ifjiyk>ifiky=mig+miy¨iTi+j<i(rijxfjiyrijyfjix)k>i(rikxfikyrikyfikx)=Iiθ¨ii= 2, 3, , 13, 1,2,,7,9,, 12E18

There are three equations for each link (there are 23 links, excluding the supporting foot), and the system describing the dynamics of the whole mechanism consists of 69 linear equations. The system (Eq. 18) is expressed in a matrix form (Eq. 19), and then solved with a MATLAB® via matrix inversion:

[a11a12a21a22]·[f12xf12yT8]=[m2x¨2m2g+m2y¨2I2θ¨2][A(coefficient)]·[F(force)]=[I(inertia)][F]=[A]1·[I]E19

The code of the matrix that must be written in MATLAB® is shown in Figure 9.

Figure 9.

The matrix A (coefficient).

Using MATLAB® code the kinematical and dynamical equations have been implemented in order to obtain solutions depending on a set of parameters (link dimensions, masses and densities, motor angular velocity) entered by the user. In Figure 10, the MATLAB flow chart of the kinematics and dynamics algorithm is shown.

Figure 10.

MATLAB® flow chart for the PASIBOT kinematics and dynamics calculus code.

The MATLAB® program first finds the corresponding value of ϑ8(t), and using Eqs. (6) to (11), it obtains the corresponding values of the rest of the angles and the positions of the centers of mass. Then, it calculates the kinematics. These data, which define the state of the biped at the time t, form the inertia matrix, (I), in Eq (19). Finally the MATLAB® program inverts the coefficient matrix, (A), by means of a matrix inversion subroutine and multiplies both matrices to provide the forces and torques between links at this time step. These values are stored to be plotted and the calculations are restarted for the next time step.

The results have been validated by comparison with others programs (working model and ADAMS code). The main advantage of the program developed via MATLAB® is that it lets us perform fast modifications, making the final robot design easier by changing parameters.

As the first result, the program implemented in MATLAB® has calculated the motor torque required to perform the movement. Figure 11 shows the actuator torque in the crank (link number 8) related to time, for different values of the motor angular velocity, ω8.

Figure 11.

Torque for different crank velocities. T is the period for one step.

In Figure 11, the same shape for each case can be appreciated apart from the torque obtained from the highest velocity value (ω8 = 5 rad/s). With this velocity, the dynamical effect of the inertia forces becomes important. However, for low speeds (below 3 rad/s) torque graphs hardly differ from one to another.

In Figure 12, the torque is represented again but for different loads (5, 10, and 15 kg) added to the hip. It is an interesting result. It shows that the required motor torque depends slightly on the added load. This is because the hip remains at almost the same level in a course of a step.

Figure 12.

Actuator torque for different hip extra loads.

6. Forward dynamics for biped robot PASIBOT

The dynamics of mechanical systems can be modeled in two ways: inverse dynamics, which calculates the forces and torques that produce kinematics (movement), and forward dynamics, which computes the movement from known applied forces and torques.

When addressing forward dynamics, the kinematics is unknown. However, the angular position of the motor crank, ϑ8, defines the position of the remaining links by Eqs. (20)–(22). These functions were defined in Eqs. (6)–(12). The corresponding angular velocities and accelerations as well as the center of mass linear velocities and accelerations are obtained from the corresponding first and second time derivatives:

ϑ˙i=dϑidϑ8dϑ8dt=ϑi'(ϑ8)ϑ˙8ϑ¨i=ϑi''(ϑ8)ϑ˙82+ϑi'(ϑ8)ϑ¨8E20
X˙i=Xi'(ϑ8)ϑ˙8X¨i=Xi''(ϑ8)ϑ82+Xi'(ϑ8)ϑ¨8E21
y˙i=yi'(ϑ8)ϑ˙8y¨i=yi''(ϑ8)ϑ82+yi'(ϑ8)ϑ¨8E22

When the kinematics is unknown, Eq. (19) becomes a system of second-order differential equations. To solve it numerically, in addition to time discretization, a motor crank angle ϑ8 discretization is proposed. In this way, the derivatives of the known functions, ϑi = ϑi(ϑ8), Xi = Xi(ϑ8) and yi = yi(ϑ8), are computed with respect to ϑ8 as follows:

{ϑi'(ϑ8)1Δϑ8[ϑi(ϑ8+Δϑ8)ϑi(ϑ8)]ϑi''(ϑ8)1(Δϑ8)2[ϑi(ϑ8+2Δϑ8)2ϑi(ϑ8)+ϑi(ϑ8)]E23
{Xi'(ϑ8)1Δϑ8[Xi(ϑ8+Δϑ8)Xi(ϑ8)]Xi''(ϑ8)1(Δϑ8)2[Xi(ϑ8+2Δϑ8)2Xi(ϑ8)+Xi(ϑ8)]E24
{yi'(ϑ8)1Δϑ8[yi(ϑ8+Δϑ8)yi(ϑ8)]yi''(ϑ8)1(Δϑ8)2[yi(ϑ8+2Δϑ8)2yi(ϑ8)+yi(ϑ8)]E25

In fact, to implement the forward dynamic problem, Eqs. (23)–(25) are inserted into Eqs. (20)–(22) and then into Eq. (18). Thus, we obtain a system of equations in which the first and second time derivatives of ϑ8 are unknowns, while the torque is now a known function of time, T8 = T8(t). The resulting system of equations is as follows:

{j<ifjixk>ifikxmiXi'(ϑ8)ϑ¨8=miXi''(ϑ8)ϑ˙82j<ifjiyk>ifikymiyi'(ϑ8)ϑ¨8=mig+miyi''(ϑ8)ϑ˙82j<iTjik>iTik+j<i(rijxfjiyrijyfjix)k>i(rikxfikyrikyfikx)Iiϑi'(ϑ8)ϑ¨8=Iiϑi''(ϑ8)ϑ˙82E26

Figure 13.

Evolution of the motor crank angle for different motor torques, when the biped starts walking from rest, according to the program described in this chapter (a) and according to Working Model 2D simulations (b).

To solve this system of differential equations, Eq. (26), a time discretization is used. For each time step, a linear inhomogeneous system is calculated, where the forces between links, and the angular acceleration of the motor crank, ϑ¨8, are unknowns, while the angular velocity of the motor crank, ϑ˙8, is known. In fact, ϑ˙8is assigned an initial input, ϑ˙8(t=0), which is updated after solving Eq. (26) in the previous time step, regarding the determined angular acceleration as constant during ∆t. Then the updated angular velocity is as follows:

ϑ˙8[nΔt]=ϑ˙8[(n1)Δt]+ϑ¨8ΔtE27

Thus, the coefficient matrix is obtained from Eq. (19) by eliminating the column corresponding to the torque coefficients T8 (previously, it was unknown), and adding a column representing the coefficients of the motor crank angular acceleration, ϑ¨8. The torque T8 must now be included in the right-hand side (RHS) column vector. The forward dynamics system is obtained as

A=[a11a12a1,710a21a22a2,710a31a32a3,71m2x2'(ϑ8)a41a42a4,71m2y2'(ϑ8)a51a52a5,71I2ϑ2'(ϑ8)a71,1a71,2a71,71I12ϑ12'(ϑ8)];U=[f01xf01yf12xf12y.f10',12'yϑ¨8];C=[0m1gm2X2''(ϑ8)ϑ˙82m2g+m2y2''(ϑ8)ϑ˙82I2ϑ2''(ϑ8)ϑ˙82T8][A(coeff)]·[U(forces,ϑ¨8)]=[C(cts)][U]=[A]1·[C]E28

Note that the torque T8 now appears in the RHS column vector (constants column). Now some results from the developed forward dynamics model are presented.

As shown in Figure 13(a), the torque above which the biped can start walking is 0.84 Nm. The initial value of ϑ8is π/2 (both feet are on the floor).

These results were compared with those obtained with the software Working Model 2D. Comparing Figure 13(a) and (b), we can say that the results are similar enough to each other to validate the program described in this chapter.

The MATLAB® program allows changing the density. In Figure 14, the motor crank angle is plotted for a constant torque, T8 = 1 Nm and varying total weight (obtained by varying the density of all links).

Figure 14.

Evolution of the motor crank angle, for a constant motor torque, T8 = 1 Nm, and for different total weights.

7. Dynamics of the biped robot PASIBOT including slippage between the ground and the supporting foot

7.1. Inverse dynamics with slippage

In the previous sections, we have applied inverse dynamics to parametrically calculate the required torque at the sole motor for PASIBOT to walk at a steady state (constant speed) with no sliding between the supporting foot and the floor. However, when sliding between the supported foot and the floor is allowed the kinematics is unknown and other approaches must be applied. In fact, three more equations regarding the supporting foot dynamics must be considered, as well as the kinematics-statics friction condition.

The forces acting on the supporting foot are

f01xf12xf112x=m1x¨1;f12yf12yf112y=m1gE29

Figure 15.

Forces acting on the supporting foot (link 1). Link 1 is connected to links 2 and 12. The floor is considered as link 0.

Since r10xandf01yrare both unknown, Eq. (30) is the non-linear torque equation for the supporting foot (link 1; see Figure 15):

r10xf01yr10yf01x(r12xf12yr12yf12x)(r1,12xf1,12yr1,12yf1,12x)=0E30

Eq. (30) is non-linear. If Eq. (18) is solved without Eq. (30), the latter can be used to obtain the instantaneous zero moment point (ZMP) relative to the center of mass, ZMPx = r10xr10x, which determines whether the biped topples.

In summary, the “inverse dynamic static friction equation” is obtained as follows in a matrix form:

[a11a12a1,71a21a22a2,71a71,1a71,2a71,71]·[f01xf01yf12xf12yT8]=[0m1gm2X¨2m2g+m2y¨2I2ϑ¨2I12'ϑ¨12'][A(coefficients)]·[F(force)]=[I(inertia)][F]=[A]1·[I]E31

The forces and the motor torque in each time step are computed by solving Eq. (31) via matrix inversion encoded in MATLAB®.

When the supporting foot is allowed to slide, x1 becomes an independent variable (in general, x10,x˙10,x¨10). To consider sliding between two bodies involves three possible scenarios: (1) no sliding (static friction), (2) imminent sliding, and (3) actual sliding. To determine the sliding status at each time step, conditional branching was incorporated into the code.

Initially, no sliding is assumed and the state is static friction (x1=0;x˙1=0;x¨1=0). Thus, solving Eq. (31), the values of friction force (f01x) and normal force (f01y) are calculated. Then, the static friction condition

|f01x|μs|f01y|E32

is evaluated for a specified static friction coefficient, μs. If Eq. (31) is satisfied, the time is incremented by Δt and Eq. (31) is recalculated using the updated values of rij,x¨i,y¨i,ϑ¨i. This step closes the “static friction” conditional loop.

If Eq. (32) is false, the PASIBOT enters the state of imminent slipping. This system has the following conditions: First, since the acceleration of the supporting foot,x¨1, is no longer zero but unknown, it must appear in the column vector of unknowns. Second, the kinetic frictional relationship between normal and tangential components of the floor-foot force must be considered:

|f01x|=μk|f01y|,E33

where μkis the kinetic friction coefficient.

The new matrix of coefficients is then obtained from its predecessor by adding the following:

  • A column of mi elements in positions corresponding to the x components of Newton’s equations, with zeros elsewhere.

  • An additional row incorporating Eq. (34).

Therefore, the final matrix form of the “inverse dynamic sliding friction equation” is as follows:

[a11a12a1,71m1a21a22a2,710a31a32a3,71m2a71,1a71,2a71,7101±μk00]·[f01xf01yf12xf12yT8f10',12'yx¨1]=[0m1gm2X¨2m2g+m2y¨2I2ϑ¨2I12'ϑ¨12'0][A(Coefficients)]·[U(Unknowns)]=[C(Constants)][U]=[A]1·[C]E34

The sign of the friction force is depending on the sliding status in the previous calculation. The friction force opposes the horizontal component of the other forces acting on the foot when the previous state was imminent sliding. The friction force opposes the velocity of the supporting foot when the previous state was actual sliding.

With Eq. (34), the MATLAB® program calculates the acceleration with which the supporting foot has begun sliding, x¨1. In the current time interval ((n1)ΔTnΔt), the supporting foot is assumed to move with the calculated uniform acceleration. The velocity and position of the supporting foot are then updated by Eq. (35):

x¨1=constant in the interval Δtx˙1[nΔt]=x˙1[(n1)Δt]+x¨1Δtx1[nΔt]=x1[(n1)Δt]+x˙1[(n1)Δt]·Δt+12x¨1(Δt)2E35

The MATLAB® program obtain the kinematics and dynamics data (torque and force data) for all links, and then it increments the time by ∆t and re-solves Eq. (34). Note that the square brackets show the dependence.

When the PASIBOT is having a slippage, the friction force is against the supporting foot movement. This force is considered constant during this time interval, and it can stop the sliding of the PASIBOT but not change the direction of the movement. This is obtained with the stopping time, ts, and compares it with the time increment, ∆t, as shown in Eq. (36):

ts=x˙1x¨1E36

If ts is positive and less than ∆t, then friction stop the PASIBOT sliding before the end of the time interval. Else, if it becomes negative or exceeds ∆t, the PASIBOT continues sliding in the time interval. After that, the MATLAB® program updates the results using ts instead of ∆t and returns to the beginning to solve the case of static friction, as provided in Eq. (31).

Figure 16.

MATLAB® program flowchart with slippage.

In Figure 16, it is shown the MATLAB® program inverse dynamics flowchart with slippage. The forces between links and torque (dynamical variables) are calculated during the “STORING DATA” task. The position, velocity, and acceleration (kinematic unknowns) for the supporting foot are updated according to Eq. (34).

7.2. Forward dynamics with slippage

In the previous sections, we have applied this methodology to design the MATLAB® program of inverse dynamics with slippage. To obtain the forward dynamics program, the slippage is treated as in 30, while the dynamics is formulated as in 34. The MATLAB® flowchart showed in Figure 16 is mostly maintained, but the systems of equations are different:

  1. – The equation system to be calculated in the state of “STATIC FRICTION” is the forward dynamics system of Eq. (28): static system (ST).

  2. – The equation system that describes the slippage of PASIBOT (Eq. 35) in the state of “SLIDING FRICTION,” is added to Eq. (28). Also, the motor torque appears in the constants’ column and the motor acceleration and the slippage acceleration become the penultimate and final elements of the unknowns’ column. Using Eqs. (20)–(25), the first and second derivatives with respect to ϑ8 of the position, velocity, and acceleration of every link are calculated, with a sufficiently fine discretization of ϑ8. The resulting system of equations (36) is referred to as the “sliding system (SL).”

    [a11a12a1,710m1a21a22a2,7100a31a32a3,71m2x2pf'(ϑ8)m2a41a42a4,71m2y2'(ϑ8)0a51a52a5,71I2ϑ2'(ϑ8)0I80a71,1a71,2a71,71I12'ϑ12'(ϑ8)01±μk000].[f01xf01yf12xf12y.f10',12'yϑ¨8x¨1]=[0m1gm2X2''(ϑ8)ϑ˙82m2g+m2y2''(ϑ8)ϑ˙82I2ϑ2''(ϑ8)ϑ˙82T8.][A(coeff)]·[U(forces,ϑ¨8,x¨1)]=[C(cts)][U]=[A]1·[C]E37

  3. – In “STORING DATA” mode, the kinematics of the supporting foot is updated by Eq. (37). The motor crank position and velocity is updated according to Eq. (35).

Some results are presented in the following figures, in which different friction coefficients and motor crank velocities have been considered. Figure 17 plots the horizontal supporting foot position as a function of time for a constant motor crank angular velocity, ϑ˙8=3rad/s,and varying friction coefficient (here μs = μk μ). From this plot, we can deduce the time course of the supporting foot sliding.

Figure 17.

Supporting foot versus time (in units of one period, T) for constant ϑ8 = 3 rad/s, for a set of different friction coefficient.

In Figure 17, it is shown that the minimum friction coefficient that prevents the slippage is 0.08. Also, it can be observed that the slippage occurs during preferred phases. The sliding starts at mid-step until when the swinging leg has reached its highest point. If two slippages occur, one of them is again invoked at mid-step, while the other occurs at the first quarter step. For μ = 0.03, slippage occurs repeatedly at various phases.

Figure 18 shows the sliding characteristics for constant friction coefficient μ = 0.1 but different motor crank angular velocities.

Figure 18.

Supporting foot versus time (in units of one period, T) for μ = 0.1, varying ϑ8.

Because the program is parametric it is easy to set different values for static and kinetic friction coefficients (kinetic friction coefficient is smaller than static friction coefficient). Figure 19 shows the supporting foot of the PASIBOT with slippage for a static friction coefficient, μs = 0.2, and three different kinetic friction coefficient, μk = 0.2, 0.1, and 0.05.

Note that the bigger the kinetic friction coefficient, the smaller the sliding distance slippage occurs. Also, if there is slippage, it occurs at the same point for all three cases.

Figure 19.

Supporting foot versus time for constant ϑ8 = 5 rad/s, for the same static friction coefficient and for three different kinetic friction coefficients.

8. Applied quasi-static approach methodology to UGV

The quasi-static methodology is applied in this chapter. The main advantage of this approach is that it is easily optimized. It is applied to a vehicle in order to optimize the navigation capabilities. This methodology is applied to an UGV shown in Figure 20, which was designed and developed by the Tallinn University of Technology. This unmanned ground vehicle can change the angle between the body and the legs to improve the capabilities of passing obstacles or navigation. It changes the position of the center of mass (CoM) relative to the ground-wheel contacts, as well as the distance between the ground and the body [14].

In order to explain how to apply the methodology to implement this analysis in MATLAB® code, the nomenclature and geometry of the vehicle are presented. The position and/or trajectories of centers of mass, joints, and ground-wheel contact points are defined. Then the quasi-static model is developed and the equations to calculate the forces and torques involved are implemented in MATLAB®. The algorithm with the quasi-static equations obtains the position along the track for any configuration angles, and then calculates the optimal values of those angles that satisfy a given condition. If the vehicle slips or overturns at any point of the track, it is also calculated by the program [15, 16].

Figure 20.

The unmanned ground vehicle (UGV) of Tallinn University of Technology.

The nomenclature is shown in the following list:

  1. Flw: force on the rear wheel exerted by the rear leg

  2. Fwl: force on the front leg exerted by the front wheel

  3. Fr: friction force exerted by the ground on the rear wheel

  4. Ff: friction force exerted by the ground on the front wheel

  5. Nr: normal force exerted by the ground on the rear wheel

  6. Nf: normal force exerted by the ground on the front wheel

  7. M: torque on the wheels (supposedly the same on front and rear wheels)

  8. D: ground-front wheel contact point; T: ground-back wheel contact point

  9. L = 0.83 m: body length; l = 0.35m: leg length; C: Center of mass

  10. φr: angle between rear leg and body; φf: angle between front leg and body

  11. βr,βf: rear wheel contact angle, front wheel contact angle

  12. f(x): function defining the track profile

  13. g(x): function defining the trajectory for the centers of the wheels

  14. mw = 50 kg: wheel mass; mb = 300 kg: body mass; ml = 20 kg: leg mass

First, the problem of positions is resolved. For a given ground function, with the front wheel position, xf and the configuration angles, φr and φf, the back wheel can be located. The slope of the ground at x is obtained as: βx = tan-1[f'(x)]. After the locus of the centers of the wheels, g(s) are calculated as follows (see Figure 21):

s=xrsinβxf(x)g(s)=f(x)+rcosβxE38

Figure 21.

Position problem.

The front wheel center position (sf,g(sf)) is obtained from Eq. (38) and the wheelbase is solved as following:

R=(Llcoscosφrlcoscosφf)2+l2(sinsinφrsinsinφf)2E39

The intersection between the function g(s) and the circumference of radius R centered on the front wheel (xfrsinsinβf,f(xf)+rcoscosβf)locates the back wheel center (see Figure 22).

Figure 22.

Scheme for locating the back wheel center.

Once the position of the vehicle is established, the locations of the CoMs are solved. The quasi-static approach can be calculated with three subsystems: (see Figure 23): (1) back wheel, (2) back wheel, both legs and body, and (3) the whole vehicle.

Figure 23.

Vehicle subsystems.

For each subsystem, equilibrium requires two force equations and one moment equation:

a){FrcosβrNrsinβr+Flwx=0Frsinβr+Nrcosβr+Flwy=mwgFrr+M=0E40
b){FrFrcosβrBrNrNrsinβrBr+Fwlx=0FrFrfsinβrBr+NrNrcosβrBr+Fwly=(mw+2ml+mb)g(C2T)x(Frsinβr+Nrcosβr)(C2T)y(FrcosβrcosβrNrsinβrsinβr)++(C2E)yFwlx(C2E)xFwlyM=0E41
c){FrcosβrNrsinβr+FfcosβfNfsinβf=0Frsinβr+Nrcosβr+Ffsinβf+Nfcosβf=(2mw+2ml+mb)g(C3T)x(Frsinβr+Nrcosβr)(C3T)y(FrcosβrNrsinβr)++(C3D)x(Ffsinβf+Nfcosβf)(C3D)y(FfcosβfNfsinβf)=0E42

where the unknowns are: Fr, Nr, Ff, Nf, FlwxFlwx, FlwyFlwy, FlwxFwlx, FlwyFwly and M. The torque is the same on front and back wheels. This system of nine equations can be simplified into a new system of five equations, where the unknowns are the torque, normal forces, and friction forces. And in a matrix form: A⋅F = C ⇒ F = A-1 C, where

A= [coscosβrsinsinβrcoscosβfsinsinβf0sinsinβrcoscosβrsinsinβfcoscosβf0r0001(C2T)xsinsinβr(C2T)ycoscosβr)(C2T)xcoscosβr+(C2T)ysinsinβr(C2E)xsinsinβf(C2E)ycoscosβf(C2E)xcoscosβf+(C2E)ysinsinβf1(C3T)xsinsinβr(C3T)ycoscosβr(C3T)xcoscosβr+(C3T)ysinsinβr(C3D)xsinsinβf(C3D)ycoscosβf(C3D)xcoscosβf+(C3D)ysinsinβf0]
 F=[FrNrFfNfM]C=[0(2mw+2ml+mb)g0(C2E)x(mbg)0]E43

MATLAB® code is used to solve the system of five equations. The program calculate for a set of discretized values of the front wheel position, the needed torque for any combination of a set of discretized values of the angles ϕr and ϕd. Thus, different criterions can be applied: minimize the energy to be supplied to the wheels, minimize the instantaneous torque or maximization the grip, the ground-wheel normal force, etc.

Some results are presented for a soft bump: square exponential profile, f(x)=0.1e(x4)2(in meters). Figure 24 shows the torque function that must be applied for any static configurations angles passing the soft bump.

Figure 24.

Torque needed passing a soft bump for different static configuration angles.

Figure 25(a) shows configuration angles variation needed to pass over the soft bump, with the minimum variation of torque and the corresponding torque function. Figure 25(b) shows the sequence of the UGV passing through the soft bump.

Figure 25.

(a) Configuration angles and torque and (b) sequence of the UGV.

9. Conclusions

The methodology provided in this chapter can be applied to mechanisms, vehicles, or robots for the complete mechanical study. The kinematics and dynamics are solved using Newton-Euler equations, from the movement of the actuator to iterating during the time of initial condition as well as external forces, and with the quasi-static approach.

The programs have been designed so that all parameters can be modified. It was possible to automate these calculations creating an algorithm implemented in a programming language to simply find the solutions and the results of the analysis.

The methodology has been applied to design the biped robot PASIBOT. The kinematics and dynamics (both forward and inverse) of the biped robot “PASIBOT,” taking into account for support foot slippage are encoded in MATLAB® code.

The great advantage of creating a parametric MATLAB® code following this methodology is that the algorithm can be modify to obtain the results in a parametric way or even changing the conditions easily. For example, it can calculate the motion of the biped from the torque function given by the biped’s sole motor or the torques required for starting and braking as well as defining the conditions that prevent or control slippage.

Because the program remains parametric the lengths, densities and masses, motor velocities and torque, friction coefficients and other parameters can be modified by the user.

The methodology was also applied to another machine, a UGV vehicle, obtaining navigation optimization results. A numerical program based on a quasi-static half vehicle model is presented. For a given profile that could be read by sensors the program calculates how the angles between the body and the legs must vary, in order to fulfill the criteria like maintain as constant as possible the torque for example. The program created with this methodology in MATLAB® code also can calculate the values of normal and friction forces, checking if the UGV rolls over or slip at any point.

In conclusion, this methodology can help to generate MATLAB® programs that will be valuable tools to optimize some navigation capabilities, dynamics analysis, quasi-static analysis, and slippage control among other.

Following this link the reader can find some examples of MATLAB codes done with the methodology of the chapter: http://www.mathworks.com/matlabcentral/profile/authors/7854464-eduardo-corral

The code that calculates the inverse dynamic of the biped PASIBOT with slippage (using the methodology that we explain in the chapter) and the code that calculates the torque of the UGV and that optimized the best route are in the previous links.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

E. Corral, J. Meneses and J.C. García-Prada (July 7th 2016). Forward and Inverse Dynamics and Quasi-Static Analysis of Mechanizes with MATLAB®, Applications from Engineering with MATLAB Concepts, Jan Valdman, IntechOpen, DOI: 10.5772/63372. Available from:

chapter statistics

1110total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Design, Simulation, and Control of a Hexapod Robot in Simscape Multibody

By Claudio Urrea, Luis Valenzuela and John Kern

Related Book

First chapter

Distributionally Robust Optimization

By Jian Gao, Yida Xu, Julian Barreiro-Gomez, Massa Ndong, Michalis Smyrnakis and Hamidou Tembine

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us