Modelling and linear stability analysis of highly mobile nonholonomic multibody systems

Universidad de Sevilla Doctoral Thesis Modelling and Linear Stability Analysis of Highly Mobile Nonholonomic Multibody Systems Author: Alfonso Garc´ıa-Ag´undez Blanco Supervisors: Daniel Garc´ıa Vallejo Emilio Freire Mac´ıas A thesis submitted in fulfilment of the requirements for the degree of Doctor of Philosophy in the PhD Programme of Mechanical and Industrial Engineering Department of Mechanical Engineering and Manufacturing September 19, 2023 This PhD thesis has received funding from the Ministry of Science, Innovation and Universities of Spain under the Training Program for Academic Staff (FPU, by its Spanish acronym) [grant FPU18/05598], and from the mobility grants for research stays for FPU recipients (years 2021 and 2023) [grants EST21/00635 and EST23/00475]. Esta tesis ha sido financiada por el Ministerio de Ciencia, Innovaci´on y Universidades en el marco del programa para la Formaci´on de Profesorado Universitario (FPU) [ayuda FPU18/05598], y por las ayudas de movilidad para realizaci´on de estancias breves para beneficiarios FPU (a˜nos 2021 y 2023) [ayudas EST21/00635 y EST23/00475]. Acknowledgements Quiero dar las gracias a mis padres, por regalarme una vida llena de felicidad, y a mis hermanos, con los que he crecido y he compartido todo. Este trabajo ha sido posible gracias a vosotros. A Paula, por tu ayuda continua y la oportunidad de haber compartido contigo la estancia en Eindhoven. Es una experiencia que siempre recordar´e. A mis amigos, gracias por el d´ıa a d´ıa, buenos ratos y tantos momentos vividos. Quiero estar con vosotros siempre. Gracias a Dani y a Emilio, por haberme regalado los a˜nos m´as felices de mi vida acad´emica. No voy a olvidar nunca nuestras reuniones y buenos momentos. Haber podido aprender de vosotros, tanto desde el punto de vista investigador como personal, ha sido un privilegio. Sin duda, lo mejor de esta tesis doctoral es vuestra amistad. I would like to thank Professor Aki Mikkola for his warm welcome and his attention during my research stay at Lappeenranta University of Technology. The gratitude also extends to the LUT Department of Mechanical Engineering, where I was able to enjoy an excellent working atmosphere from the very beginning. I would also like to thank Professor Alessandro Saccon for giving me the opportunity to visit the Dynamics and Control group of the Eindhoven University of Technology, and for his assistance during my research stay. Abstract This PhD thesis aims to develop a methodology that contributes to the modelling of nonholonomic multibody systems, the study of their linear stability and the design of linear feedback controllers. The ultimate objective is on single-person nonholonomic vehicles, particularly in terms of safety in urban transportation. By understanding how critical design parameters affect stability, it will be possible to carry out design modifications that result in safer vehicles and reduce accidents. These modifications will also cater to a wider range of potential users, including elderly and physically impaired individuals. From a regulatory perspective, understanding stability conditions will play a crucial role in formulating rules for the use of these vehicles within cities. Moreover, advancements focused on enhancing safety will promote the adoption of electric single-person vehicles, leading to a reduction in environmental degradation and aiding in the accomplishment of pollution reduction objectives. The methodology consists of the following phases. First, multibody models of a class of nonholonomic systems are developed. Subsequently, in order to carry out the linear stability analyses of these multibody systems, efficient and accurate linearization approaches are required. To this end, novel linearization procedures, for multibody systems with holonomic and nonholonomic constraints, are proposed. The correctness of the linearization approaches is validated with the linear stability results of a well-acknowledged bicycle benchmark. The results show that the procedures are completely accurate, efficient, valid for any multibody system (regardless of its complexity) and powerful, obtaining the linearized equations of motion along an arbitrary reference solution as a function of the design parameters of the multibody system. Next, by varying these parameters over a wide range of values, a detailed sensitivity analysis of the eigenvalues can be performed to assess the linear stability of the multibody system. Lastly, linear feedback controllers can be designed by using the linearized equations of motion. These linearization approaches are applied to study the linear stability of several nonholonomic systems. First, linear time-invariant (LTI) systems are addressed. In particular, the linear stability of the steady forward motions of classical nonholonomic systems, as the skateboard or the rolling hoop, and highly mobile nonholonomic multibody systems, such as the bicycle, the waveboard and the electric kickscooter, is assessed. Next, the case of periodic linear time-varying (LTV) systems is studied, including a detailed linear stability analysis of the circular steady motion of the rolling toroidal wheel. Another application is the use of the linearization procedures in multiphysics scenarios, and, in particular, with hydraulically actuated multibody systems. Finally, the approaches are applied in the design of linear feedback controllers for nonholonomic multibody systems. As an example of application, the optimal control of the well-known ball-plate system, using a Linear-Quadratic Regulator (LQR), is shown. Contents 1 Introduction 17 1.1 Multibody System Dynamics . . . . . . . . . . . . . . . . . . . . . 17 1.2 Nonholonomic multibody systems . . . . . . . . . . . . . . . . . . . 21 1.3 Linearization in Multibody System Dynamics . . . . . . . . . . . . 31 1.4 General motivation and research objectives . . . . . . . . . . . . . 36 1.5 By-chapter summary . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.6 Publications of this PhD thesis . . . . . . . . . . . . . . . . . . . . 42 2 Modelling of constrained multibody systems 47 2.1 Basic concepts in multibody systems . . . . . . . . . . . . . . . . . 48 2.1.1 Degrees of freedom of a multibody system . . . . . . . . . . 48 2.1.2 Types of coordinates . . . . . . . . . . . . . . . . . . . . . . 48 2.1.3 Global and local reference frames: the orientation matrix ofabody ............................ 52 2.2 Kinematics of multibody systems . . . . . . . . . . . . . . . . . . . 56 2.2.1 Velocity vectors: linear and angular velocity . . . . . . . . . 57 2.2.2 Acceleration vectors: linear and angular acceleration . . . . 59 2.2.3 Velocity and acceleration vectors in terms of the generalized coordinates ........................... 60 2.2.4 Kinematic constraint equations . . . . . . . . . . . . . . . . 62 2.2.5 Kinematic simulation . . . . . . . . . . . . . . . . . . . . . 63 2.3 Dynamics of multibody systems . . . . . . . . . . . . . . . . . . . . 66 2.3.1 Dynamics of a rigid body: Newton-Euler equations . . . . . 66 2.3.2 Dynamic equations of multibody systems . . . . . . . . . . 67 2.3.3 Equations of motion of constrained multibody systems . . . 69 2.3.4 Inverse dynamics analysis . . . . . . . . . . . . . . . . . . . 69 2.3.5 Forward dynamics analysis . . . . . . . . . . . . . . . . . . 70 2.4 Forms of the equations of motion . . . . . . . . . . . . . . . . . . . 71 2.4.1 Index-3 and index-1 DAE form . . . . . . . . . . . . . . . . 71 2.4.2 ODE system with coordinate partition in terms of independent and dependent coordinates . . . . . . . . . . . . . . . 72 2.4.3 Equations of motion in DAE form without Lagrange multipliers .............................. 73 5 Contents 3 An overview of linearization approaches in Multibody System Dynamics 77 3.1 Group 1: procedures based on the linearization of the index-3 DAE system, using a redundant set of coordinates . . . . . . . . . . . . 79 3.2 Group 2: procedures based on the linearization of the index-1 DAE system, using a redundant set of coordinates . . . . . . . . . . . . 86 3.3 Group 3: procedures based on coordinate partitioning in terms of independent and dependent coordinates . . . . . . . . . . . . . . . 88 3.4 Group 4: procedures based on coordinate partitioning in terms of independent and dependent coordinates, with holonomic constraintselimination........................... 91 4 Novel linearization approaches for constrained multibody systems 97 4.1 Procedures based on symbolic computation . . . . . . . . . . . . . 99 4.1.1 First symbolic procedure . . . . . . . . . . . . . . . . . . . 99 4.1.2 Second symbolic procedure . . . . . . . . . . . . . . . . . . 101 4.2 Procedures based on numerical computation . . . . . . . . . . . . . 105 4.2.1 First numerical procedure . . . . . . . . . . . . . . . . . . . 105 4.2.2 Second numerical procedure . . . . . . . . . . . . . . . . . . 109 4.2.3 Third numerical procedure . . . . . . . . . . . . . . . . . . 112 4.3 Validation of the linearization approaches with the bicycle benchmark117 4.3.1 Multibody model of the bicycle benchmark . . . . . . . . . 117 4.3.2 Description of the steady forward motion of the bicycle . . 122 4.3.3 Use of the first numerical procedure . . . . . . . . . . . . . 124 4.3.4 Use of the first symbolic procedure and second numerical procedure............................127 4.3.5 Use of the second symbolic procedure and third numerical procedure............................131 4.3.6 Comparison with the benchmark results . . . . . . . . . . . 135 4.4 Comparison: advantages and disadvantages of the approaches . . . 138 4.5 The role of the dependent coordinates in the linearization . . . . . 142 4.5.1 Test example: bicycle multibody model with toroidal wheels 143 4.5.2 The role of the dependent coordinates in the linear stability analysis of the bicycle . . . . . . . . . . . . . . . . . . . . . 146 5 Linear stability analysis of nonholonomic multibody systems: study of linear time-invariant systems 155 5.1 Stability analysis of linear time-invariant systems: Lyapunov theory 156 5.2 Skateboard: Chaplygin sleigh . . . . . . . . . . . . . . . . . . . . . 158 5.2.1 Description of the multibody model and equations of motion 158 5.2.2 Linearization with the third numerical procedure . . . . . . 160 5.2.3 Linear stability analysis of the skateboard: discussion of the results..............................161 5.3 Hoop rolling without slipping . . . . . . . . . . . . . . . . . . . . . 163 6 Contents 5.3.1 Description of the multibody model and equations of motion 163 5.3.2 Linearization with the third numerical procedure . . . . . . 166 5.3.3 Linear stability analysis of the rolling hoop: discussion of theresults............................167 5.4 Bicycle with elliptic toroidal wheels . . . . . . . . . . . . . . . . . . 170 5.4.1 Linearization with the third numerical procedure . . . . . . 170 5.4.2 Linear stability analysis of the bicycle with toroidal wheels: discussion of the results . . . . . . . . . . . . . . . . . . . . 173 5.5 Thewaveboard.............................176 5.5.1 Description of the multibody model and equations of motion 177 5.5.2 Linearization with the third numerical procedure . . . . . . 182 5.5.3 Linear stability analysis of the waveboard: discussion of the results..............................184 5.5.4 Inverse dynamics analysis: influence of the design parameters on the external actuations . . . . . . . . . . . . . . . . 189 5.6 Electric kickscooter with toroidal wheels and suspensions . . . . . 192 5.6.1 Description of the e-scooter multibody model . . . . . . . . 193 5.6.2 Description and characterisation of the suspensions . . . . . 197 5.6.3 List of the e-scooter multibody model parameters . . . . . . 200 5.6.4 Equations of motion of the e-scooter multibody model . . . 202 5.6.5 Linearization with the third numerical procedure . . . . . . 203 5.6.6 Linear stability analysis of the e-scooter: discussion of the results..............................207 6 Linear stability analysis of nonholonomic multibody systems: study of linear periodic time-varying systems 219 6.1 Linear Time-Varying systems . . . . . . . . . . . . . . . . . . . . . 220 6.2 Stability analysis of periodic trajectories: Floquet-Lyapunov theory . . . . . . . . . . . . . . . . . . . . . . . 221 6.3 The rolling toroidal wheel: stability analysis of the circular steady motion..................................224 6.3.1 Description of the model and equations of motion . . . . . 225 6.3.2 Description of the circular steady motion . . . . . . . . . . 228 6.3.3 Linearization of the equations of motion . . . . . . . . . . . 230 6.3.4 Computation of the time-independent Jacobian matrix . . . 234 6.3.5 Results and discussion: linear stability of the steady straightlinerollingmotion .......................237 6.3.6 Results and discussion: linear stability of the steady spinning motion around a vertical diameter . . . . . . . . . . . 242 6.3.7 Results and discussion: linear stability of the circular steady motion..............................247 6.3.8 Validation of the results using Floquet-Lyapunov theory . . 255 7 Contents 7 Linearization of the equations of motion of multibody systems in multiphysics scenarios 257 7.1 Formulation of the problem . . . . . . . . . . . . . . . . . . . . . . 259 7.1.1 Modelling of the hydraulic system . . . . . . . . . . . . . . 259 7.1.2 Equations of motion of hydraulically actuated constrained multibodysystems.......................262 7.2 Linearization approach for hydraulically actuated multibody systems263 7.3 Validation of the approach with a hydraulically actuated four-bar mechanismmodel............................271 7.3.1 Description of the four-bar mechanism model with hydraulic actuator.............................271 7.3.2 Description of the hydraulic circuit . . . . . . . . . . . . . . 275 7.3.3 Computation of the Jacobian matrix . . . . . . . . . . . . . 277 7.3.4 Numerical results: validation of the approach . . . . . . . . 281 7.3.5 Comparison of the proposed linearization procedure with a conventional counterpart approach . . . . . . . . . . . . . . 289 8 Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems 295 8.1 Extension of the third numerical linearization approach to controlled multibody systems . . . . . . . . . . . . . . . . . . . . . . . 297 8.2 Time-varying discrete linear systems . . . . . . . . . . . . . . . . . 299 8.3 Optimal control of nonholonomic multibody systems: discrete time finite-horizon Linear-Quadratic Regulator . . . . . . . . . . . . . . 302 8.4 Example of study: ball-plate system . . . . . . . . . . . . . . . . . 303 8.4.1 Description of the ball-plate multibody model and equationsofmotion.........................303 8.4.2 Computation of the reference solution: circular steady motionoftheball .........................306 8.4.3 Linearization along the circular steady motion . . . . . . . 307 8.4.4 Controllability of the linear time-varying system . . . . . . 309 8.4.5 Controllable system of linear equations . . . . . . . . . . . . 314 9 Conclusions and future lines of research 321 9.1 Original contributions of this PhD thesis . . . . . . . . . . . . . . . 321 9.2 Possible future lines of research . . . . . . . . . . . . . . . . . . . . 327 8 List of Figures 1.1 Schematic representation of a generic multibody system. . . . . . . 19 1.2 Classical examples of nonholonomic systems: the Chaplygin sleigh andtherattleback............................ 22 1.3 Classical nonholonomic system: the rolling disk. . . . . . . . . . . . 23 1.4 Snakeboard and its propelling mechanism. . . . . . . . . . . . . . . 25 1.5 Sequence of configurations in the snakeboard motion. . . . . . . . . 25 1.6 Components of the waveboard and propelling mechanism. . . . . . 26 1.7 Meandering trajectory of the waveboard. . . . . . . . . . . . . . . . 27 1.8 Single-person vehicles in urban mobility: the electric bicycle and motorcycle. ............................... 29 1.9 SEAT eXS kickscooter [1]. . . . . . . . . . . . . . . . . . . . . . . . 30 1.10 Main objectives of the PhD thesis. . . . . . . . . . . . . . . . . . . 39 2.1 Four-bar mechanism described by using relative coordinates. . . . . 49 2.2 Four-bar mechanism described by using absolute coordinates. . . . 50 2.3 Four-bar mechanism described by using natural coordinates. . . . . 52 2.4 Rigid body iundergoing planar motion: global and local reference frames................................... 53 2.5 Examples of rotation sequences: proper Euler angles and Tait- Bryanangles............................... 55 2.6 Rigid body iundergoing spatial motion: global and local reference frames................................... 56 3.1 Existing linearization approaches in Multibody System Dynamics. 78 4.1 Novel linearization approaches presented in Chapter 4. . . . . . . . 98 4.2 Three-dimensional multibody model of the bicycle benchmark. . . 118 4.3 Lateral view of the bicycle benchmark with the main geometric parameters. ...............................120 4.4 Root locus of the bicycle benchmark model. . . . . . . . . . . . . . 137 4.5 Bicycle benchmark model: evolution of the eigenvalues with the forwardspeed. .............................138 9 List of Tables 16 Chapter 1 Introduction In this chapter, a description and brief state-of-the-art review of Multibody System Dynamics are performed. Special emphasis is placed on mechanical systems with nonholonomic constraints, exploring classical nonholonomic systems, such as the Chaplygin sleigh or the rolling disk, as well as complex nonholonomic multibody systems like the bicycle or the electric kickscooter. Since the linearization of the equations of motion of multibody systems is required to assess their linear stability, a summary with the main linearization procedures developed in the context of Multibody System Dynamics is also included. Next, the general motivation and the research objectives, a by-chapter summary and the main publications of this PhD thesis are presented in detail. 1.1 Multibody System Dynamics Multibody dynamics is the branch of mechanics devoted to the study of the dynamic behavior of multibody systems. A multibody system is an assembly of several bodies, connected to each other by joints and subject to the action of forces and motion constraints [3, 4]. Multibody system formalisms typically resort to computer-aided techniques to efficiently represent and simulate the dynamics of these systems. A wide range of engineering fields make use of Multibody System Dynamics. Some of the key applications include robotics [5, 6], to model and simulate the motion of robot manipulators, humanoid robots, wheeled mobile and aerial robots; design and analysis of mechanisms [7, 8]; biomechanics [9, 10, 11, 12], to study the human motion, design prosthetic devices or analyse sports performance; the aerospace industry, to simulate the motion of complex systems as aircraft [13, 14], helicopters [15], satellites [16] or spacecraft [17]; vehicles [18, 19, 20] and railway dynamics [21, 22, 23], for ride comfort analysis and design optimization for improved safety and efficiency; or virtual prototyping and simulation [24, 25], which leads to cost and time savings in the product development by virtually testing and optimizing the performance of systems, without 17 Chapter 1. Introduction building physical prototypes. In a multibody system, the bodies describe large translational and rotational motions, and can be either rigid or flexible. A body is said to be rigid when its deformations are assumed to be negligible, and therefore the distance between two given points of the body remains constant over time, independently of the external forces or torques acting on it. Although solids are not completely rigid, in many practical applications their flexibility can be disregarded. In contrast, a flexible body presents an elastic structure and can undergo significant deformations. A free rigid body, moving in the three-dimensional space, presents six degrees of freedom, and therefore can be fully described by using six generalized coordinates. Conversely, in the case of a flexible body, additional generalized coordinates are required to describe its deformations [3, 26, 27]. This thesis focuses on multibody systems consisting of rigid bodies. These bodies are interconnected by means of essential elements, called joints, which allow for certain relative motions and constrain others. A joint between two bodies introduces kinematical constraints, which limit the number of relative motions of one body relative to the other. Some of the most usual kinematic joints are the revolute, spherical or prismatic joints. Moreover, the forces acting on a multibody system can be the result of passive elements, such as springs and dampers, or active elements as actuators. The external applied forces present different nature and level of complexity [28, 29]. Figure 1.1 shows a schematic representation of a generic multibody system with its main components: rigid and flexible bodies, different kinematic joints (revolute and spherical), force elements (spring, damper and actuator) and external actuations (applied forces and torques). A review of the roots and the state-of-the-art of Multibody System Dynamics can be found in the work of Schiehlen [30]. The origins of Multibody System Dynamics can be traced back to the 17th and 18th centuries. Newton’s laws of motion, presented in his work “Philosophiae Naturalis Principia Mathematica” [31] in 1687, laid the foundation for classical mechanics. Newton introduced the most simple element of a multibody system: a free particle. In 1776, Euler presented, in “Nova methodus motum corporum rigidarum determinandi” [32], the concept of rigid body, which consists of a collection of material points and constitutes the paramount element of a multibody system. D’Alembert published, in 1743, the most famous of his scientific works: “Trait´e de Dynamique” [33]. This contribution contains the first statement of what is known as d’Alembert’s principle and deals with constrained rigid bodies. In [33], d’Alembert made the difference between applied and reaction forces, and labeled the reaction forces as “lost forces” [30]. Next, in 1788, Lagrange published his work “M´ecanique analytique” [34], where a consistent mathematical formulation of d’Alembert’s principle and the systematic analysis of constrained mechanical systems were established. By considering the total kinetic and potential energy of the system, along with its kinematical constraints and the generalized coordinates, Lagrange applied a variational principle [34] to obtain the Lagrangian equations of the first and second kind. The former constitutes a system of Differential-Algebraic Equations (DAE), while the latter leads to a minimal set of Ordinary Differential Equations 18 1.1. Multibody System Dynamics Figure 1.1: Schematic representation of a generic multibody system with its main components. (ODE). Lagrange, with all these contributions, is one of the founders of analytical mechanics. Another remarkable figure in the context of analytical mechanics is Hamilton, who made significant contributions with the development of the Hamiltonian formulation. In 1833, his work “On a General Method in Dynamics” [35] was published, where the concepts of canonical transformations and the Hamiltonian equations of motion were developed. These outstanding early contributions established the groundwork of classical and analytical mechanics, which form the basis of modern Multibody System Dynamics. The advances in the dynamics of constrained mechanical systems continued with notable contributions, such as the work of Jourdain [36], which presented an extension of d’Alembert’s principle valid for holonomic systems, and Kane and Levinson [37], where a method for forming the equations of motion of constrained multibody systems, including nonholonomic constraints, was developed. Some of the first applications of rigid body dynamics are found in gyrodynamics [38, 39], mechanism theory [40, 41, 42] and biomechanics [43, 44, 45, 46]. Nevertheless, the use of the dynamics of rigid bodies in practical applications faced significant limitations in terms of model complexity until the 1960s. These limitations arised from less efficient numerical methods for solving differential equations, the nonlinearity of large rotations and the highly nonlinear gyroscopic coupling in the equations of motion [30]. Computational advances at both software and hardware levels, driven by a growing demand for complex models of aircraft and satellites, gave rise to what is known today as Multibody System Dynamics. Some 19 Chapter 1. Introduction multibody formalisms, as those of Hooker and Margulies [47], and Roberson and Wittenburg [48], appeared for developing multibody systems of satellites. Moreover, the progress in computational capacity lead to the development of symbolic computation, which constitutes a powerfull tool for computing the equations of motion of complex multibody systems. A contribution in this line was that of Levinson [49]. Schwertassek and Roberson [50] described software systems for the modeling, simulation and animation offered in the 1980s. The increasing number of research groups devoted to the study of mechanical systems lead to the need for standardization. Therefore, in 1977, supported by the IUTAM (International Union of Theoretical and Applied Mechanics), the first international symposium on what began to be called Multibody System Dynamics was organized in Munich by Magnus [39]. In the same year, the first book on multibody dynamics was published by Wittenburg [51]. This book addressed the basics of rigid body kinematics and dynamics, and explored classical problems of rigid body motion as well as general multibody systems. Subsequent international events, devoted to Multibody System Dynamics, took place in the next years, being organized by NATO Advanced Study Institute [52, 53], the International Association of Vehicle System Dynamics (IAVSD) [54] and IUTAM [55, 56, 57]. Schiehlen’s book [58], published in 1986, presented multibody systems, finite element and continuous systems as equivalent models in an unified manner. In 1988, Nikravesh published the first book on computer-aided analysis of multibody systems [28]. Haug’s book [59] illustrated basic methods of computer-aided kinematics and dynamics of mechanical systems, for planar and spatial systems. Relevant textbooks, devoted to flexible multibody systems, were also published, as the well-known work of Shabana [26]. The development of computer science and numerical techniques played a key role for the consolidation of Multibody System Dynamics. Amirouche’s book [60] addressed computational methods, with special emphasis on matrix methods. Efficient methods for the kinematic and dynamic simulation to meet the real time challenge were presented by Garcia de Jal´on and Bayo [61], and Shabana’s book [62] focused on computational dynamics of rigid multibody systems. Moreover, several general-purpose programs for the analysis and simulation of multibody systems, such as ADAMS [63], RECURDYN [64] or DADS [65], were released to the market. The first number of the international journal Multibody System Dynamics was released in 1997, with two review articles authored by Shabana [66] and Schiehlen [30]. In particular, [30] presented a historial overview of the development of Multibody System Dynamics and highlighted some of the early applications that generated interest in the field. In a retrospective study published in 2007, Schiehlen [67] examined the first decade of the journal Multibody System Dynamics, outlining the research trends of that period. Presently, there are two international journals fully dedicated to this field: Multibody System Dynamics and Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. Furthermore, recent journals, as International Journal of Mechanical System Dynamics (IJMSD), have multibody dynamics as an essen- 20 1.2. Nonholonomic multibody systems tial part of their scope; many others accept contributions related to this topic; and most mechanical engineering conferences include a multibody dynamics session. Nowadays, Multibody System Dynamics is a well-established branch of theoretical, computational and applied mechanics. This field follows a continuous development, with many institutions focused on this research area [68], and offers tools in an excellent state of maturity for the development of modern engineering, science and technology. This PhD thesis is focused on the modelling, stability analysis and control of nonholonomic multibody systems. These systems are of interest in many engineering fields, as vehicle dynamics and mobile robots, and they also attract mathematicians, looking for analytical solutions in a reduced minimal state space. The next section is devoted to nonholonomic multibody systems, presenting a historial overview of the development of nonholonomic mechanics and highlighting some of the most notable nonholonomic systems. 1.2 Nonholonomic multibody systems Nonholonomic systems are physical systems subject to constraints that involve both positions and velocities and are non-integrable. A classical example of nonholonomic system is a body with rolling contact in spatial motion, as a wheel rolling without slipping. The interest in nonholonomic systems dates back to the late 19th century. The equations of motion of a nonholonomic system were derived by Ferrers [69] in 1872. These equations were given in the form of Euler–Lagrange equations, including some additional terms arising from the constraints (but without Lagrange multipliers). Next, Voronetz [70] performed, in 1901, the formal derivation of these equations. Chaplygin obtained, for the case of cyclic configuration variables, the so-called Chaplygin equations [71, 72], which engendered the technique of nonholonomic reduction. In his fundamental work Principles of Mechanics [73], Hertz first coined the term ’nonholonomic’ and showed the inapplicability of variational principles and Lagrange equations to nonholonomic systems. Hertz’s observations were developed by Poincar´e in his well-known paper Hertz’s Ideas in Mechanics [74]. Moreover, the work of Korteweg [75] clarified the confusion existing in the literature between nonholonomic mechanical systems and vakonomic systems. The derivation of the Gibbs-Appell equations of motion [76], to describe the evolution of nonholonomic systems without using Lagrange multipliers, was also a significant step in the development of nonholonomic mechanics. Hamel’s work [77] also included nonholonomic systems, and Kane and Levinson [37] presented a method for forming the equations of motion of constrained multibody systems, including nonholonomic constraints. Kalaba and Udwadia [78] obtained the equations of nonholonomic dynamical systems by using Gauss’s principle. Routh’s equations were applied to nonholonomic systems with variable mass by Yao-huang Luo and Yong-da Zhao [79]. Rismantab-Sany and Shabana [80, 81] dealt with nonholonomic deformable systems and their nu- 21 Chapter 1. Introduction merical solution by means of Differential-Algebraical Equations. The book of Neimark and Fufaev [82] was a major landmark for the development of nonholonomic mechanics, addressing the kinematics, dynamics and theory of stability of nonholonomic systems. Furthermore, since the beginning of the eighties, nonholonomic mechanics was studied from the standpoint of Geometric Mechanics. Some notable works in this framework are those of Vershik and Faddeev [83] and Koiller [84]. The work of Koiller addressed the reduction of nonholonomic systems with symmetries and inspired Bloch et al. in [85], where it is discussed that, in contrast to mechanical systems with standard configuration space constraints (holonomic), the existence of symmetries in the nonholonomic case does not necessarily lead to conserved quantities. In the context of Geometric Mechanics, the works of Bates et al. [86, 87, 88, 89], Saunders et al. [90, 91], Marle [92, 93], Ehlers et al. [94] and Ibort et al. [95] must also be highlighted. The connections of control theory with the theory of nonholonomic systems can be found in the book of Bloch [96]. Furthermore, control aspects of nonholonomic multibody systems were addressed by Schiehlen [97]. The works of Borisov and Mamaev [98], Bloch et al. [99] and Manuel de Le´on [100] can be followed to find more details about the development of nonholonomic mechanics. (a) The Chaplygin sleigh on a horizontal plane (b) Carved wooden rattleback [101] Figure 1.2: Classical examples of nonholonomic systems: the Chaplygin sleigh and the rattleback. Some of the most studied nonholonomic systems in the literature are the wellknown Chaplygin sleigh [102, 103, 104] and the rattleback [105, 106, 107, 108], also known as wobblestone or Celtic stone. The Chaplygin sleigh is a classical mechanical system, consisting of a rigid body that resembles a sleigh, moving on a horizontal plane supported at three points. The sleigh presents a knife edge, which constrains its motion so that the knife slides only along its longitudinal axis. Figure 1.2 (a) shows an example of the Chaplygin sleigh, with point A 22 1.2. Nonholonomic multibody systems representing the knife edge. Despite its simplicity, the Chaplygin sleigh constitutes an interesting example, widely used in the field of Multibody System Dynamics to study nonholonomic systems. Figure 1.2 (b) shows a rattleback, which is a small asymmetric toy with an intriguing behavior when spun. This behavior arises from the shape and asymmetry of the toy, which causes the system to exhibit a preferred direction of rotation. The system spins steadily in the preferred direction of rotation, but when it is spun in the opposite direction, it starts to wobble and eventually stops spinning, reversing its rotation to the preferred direction of spin. Figure 1.3: Classical nonholonomic system: disk rolling without slipping on a horizontal plane. Nevertheless, the rolling contact of bodies constitutes the most widely investigated nonholonomic example. The research on rolling bodies went hand in hand with the development of nonholonomic mechanics. The first known work devoted to the dynamics of a rolling rigid body is attributed to Euler [109], in which small oscillations of a body rolling without slipping were studied. The work of Poisson [110], which addressed the integrability of a heavy top moving on a smooth surface, gave rise to multiple papers investigating the integrability of the equations of motion of a disk rolling without slipping. Slesser [111] studied for the first time, in 1861, the motion of a rolling symmetrical disk on a horizontal rough plane. Vierkandt [112], in 1892, revealed that, on an appropriate symmetryreduced space, all motions of the rolling disk are periodic. Figure 1.3 shows a disk in the upright position, rolling without slipping on a horizontal plane. Chaplygin [113] showed the integrability of the rolling motion of an axisymmetric body on a horizontal plane. In the particular case of the rolling disk, the reduction of the problem to the analysis of hypergeometric quadratures was presented. Independently from Chaplygin, in 1900, Appel [114] and Korteweg [75] carried out the integration of the equations of motion of a disk with hyperelliptic functions, and Gallop [115] obtained the same results by resorting to Legendre functions. Furthermore, Carvallo [116] and Routh [117] studied the stability of the rolling disk, although these works did not include any information on the bifurcations of the system. 23 Chapter 1. Introduction Other works devoted to analyse the stability of the steady motions of the rolling disk, taking advantage of the classical integrability results, are those of Mindlin [118, 119], Duvakin [120] and Karapetyan [121]. Kolesnikov [122] and Fedorov [123] discussed qualitative properties of the disk motion, and Kozlov et al. [124] studied the nonintegrable problem of a disk moving on an inclined plane. O’Reilly [125] analysed the bifurcations and stability of the steady motions for both rolling and sliding disks, obtaining two-parameter families of integrable second order differential equations. The results showed that the bifurcations of the steady motions of the disks are either of the pitchfork or saddle-node type. The study of bifurcations and stability of the disk’s stationary motions was also performed by Cushman et al. [126] and Kuleshov [127]. Borisov et al. [128] showed different trajectories drawn by the contact point in the body-fixed and relative frames of references, and presented a complete three-dimensional bifurcation diagram in the space of the first integrals. Zenkov et al. [129] analysed, by using an energy-momentum analysis for nonholonomic systems, the stability of the relative equilibria of the rolling disk. Paris et al. [130] studied the circular rolling steady motion of the disk and derived the expression for the frequency of small oscillations about this reference motion. The book of Cushman et al. [131] provided a modern differential geometric treatment of linearly nonholonomically constrained systems and, in particular, analyzed the stability and the corresponding bifurcation diagram of the rolling disk. Moreover, a detailed linear stability analysis of the straight-line rolling, spinning about a vertical diameter and tumbling solutions of the rolling disk was presented by Przybylska [132]. The dynamics of an unbalanced disk with a single nonholonomic constraint was studied in a recent work by Kilin et al. [133]. In addition, there exist several works devoted to other nonholonomic rolling bodies. Borisov et al. [134] considered cases of existence of invariant measure, additional first integrals and Poisson structure in the problem of rigid body’s rolling without sliding on a plane and a sphere. The particular examples of nonholonomic systems addressed in [134], in addition to the rolling disk, are: a dynamically symmetric ball with the displaced center of mass rolling on a plane; the rolling of balanced, dynamically nonsymmetric ball (Chaplygin ball); the rolling of an unbalanced, dynamically nonsymmetric ball on a plane and the rolling of an ellipsoid on a plane. Moreover, different examples of bodies rolling on a sphere were included in [134]: the problem of rolling of round disk; a ball with displaced center; a balanced, dynamically nonsymmetric ball on sphere; the unbalanced, dynamically nonsymmetric ball on sphere; and the rolling of a body with partially flat surface on sphere. Other important works devoted to the rolling ball are those of Hennans [135] and Borisov et al. [136, 137, 138, 139]. A recent work of Antali et al. [140] addressed the kinematics and dynamics of a rigid body in contact with two rigid surfaces in the presence of dry friction. In the context of nonholonomic systems, mechanics and control of robotics and biological locomotion have also been widely studied by several authors [141, 142]. In addition to the skateboard [143, 144], the snakeboard is a very well-known nonholonomic example. This mechanism consists of a longitudinal bar that connects two sets of wheels, being allowed to rotate independently. In the snakeboard, the 24 1.2. Nonholonomic multibody systems (a) Snakeboard: rotations of the rear and front wheelsets (b) Maneuvering of the snakeboard Figure 1.4: Snakeboard and its propelling mechanism: the rotations of the rear and front wheelsets are ϕrand ϕf, respectively; the twist angle of the rider’s torso is ψ; and the torques exerted by the rider with their ankles and torso are Γr, Γf and Γψ, respectively. forward propelling is achieved without kicking the ground. Figure 1.4 (a) shows a view of the snakeboard, where the angles ϕrand ϕfrepresent the rotations of the rear and front wheelsets, respectively. Figure 1.4 (b) illustrates the snakeboard propelling mechanism, where the twist angle of the rider’s torso is denoted as ψ. To achieve a forward motion, the rider carries out alternate rotations with their ankles (exerting the torques Γrand Γfon the rear and front wheelsets, respectively, shown in Fig. 1.4 (b)) and torso (torque Γψ). The sequence of motion of the snakeboard is shown in Fig. 1.5. First, both platforms are turned in, and the twist of the rider’s torso through an angle ψ allows the motion depicted by configuration (1) in Fig. 1.5. Subsequently, the rider turns both feet out and moves the torso in the opposite direction, obtaining the configuration depicted as (2) in Fig. 1.5. The repetition of this sequence over time leads to the forward propelling of the snakeboard. Figure 1.5: Sequence of configurations in the snakeboard motion. Ostrowski et al. [145] addressed the motion planning of the snakeboard, using steering techiques introduced by Murray and Sastry [141], and extended their work to perform optimal gait analysis [146]. Then, Kuleshov [147] carried out a further development of the mathematical model, deriving the equations of motion 25 Chapter 1. Introduction the stability of motorcycles models [211] by manually analysing the eigenvalues of a set of linearized equations at the equilibrium configuration. When the stability of periodic motions has to be analysed, Floquet theory is the preferred methodology [208]. This kind of analysis requires the evaluation of the so-called Floquet multipliers, based on a monodromy matrix which follows by numerical integration of a linearized set of equations. This procedure has been successfully used in the context of multibody systems, see for example [212]. Obtaining linearized equations of motion of unconstrained systems is a very well-known task that is usually learnt by science and engineering students. Such equations are often used to study linear mechanical vibrations [213] or to analyse the linear stability of mechanical systems [214]. However, when the equations of motion are constrained either by holonomic or nonholonomic constraints, the linearization procedure is not simple at all. The equations of motion of constrained multibody systems are usually given by nonlinear Differential-Algebraic Equations (DAE) systems, whose correct linearization poses a challenge when the multibody system model presents a large number of coordinates and constraints. The DAE system can be reduced to a system of Ordinary Differential Equations (ODE) by resorting to a coordinate partition in terms of independent and dependent coordinates. In this case, the linearization must be conducted carefully and, in general, the dependent coordinates cannot be ignored in the computation of the correct linearized equations. In the field of Multibody System Dynamics, the number of works devoted to the linearization of the equations of motion has increased in recent years. These linearization approaches differ from each other depending on the initial form of the nonlinear equations of motion, the selection of the generalized coordinates that they use (a redundant or a minimal set) or the structure of the linearized equations of motion, which is highly dependent on the selection of coordinates. A first group of works intends to linearize the index-3 DAE equations of motion. The work of Negrut et al. [215] presented a linearization approach suitable for heterogeneous systems, which include nonholonomic constraints, friction, flexible bodies or user-defined differential equations. In this procedure, the governing equations are augmented and a set of sensitivities that provide the linearization of interest are computed. The approach was numerically tested with a nonuniform helicopter blade, a washing machine, an all terrain vehicle and a rotor-stator multibody models. Escalona et al. [216, 217] also linearized the index-3 DAE equations of motion. In [216], the stability analysis of the steady curving of vehicles was addressed. An alternative stability analysis to the use of Floquet’s theory, based on two coordinate projections, was presented. In particular, the linear stability of the circular steady motion of a simple wheeled mechanism was analysed. The same approach was used in [217] to linearize the equations of motion of railroad vehicles moving on curved tracks. Another procedure devoted to the linearization of the index-3 DAE system was described by Gonz´alez et al. [218]. This approach, known as Redundant Coordinate Set (RCS) Formulation in [218], leads to the same structure of the linearized equations motion obtained in [216, 217]. The properties of the approach were shown with the linearization of two numerical 32 1.3. Linearization in Multibody System Dynamics examples in static equilibrium: a multiple loop four-bar linkage and a flexible double pendulum. The work of Gonz´alez et al. [218] presented two other linearization procedures, whose properties were also assessed with these numerical examples. Recently, Pappalardo et al. [219] presented an approach to linearize the index-3 DAE form of the equations of motion of multibody systems with holonomic and nonholonomic constraints. The use of the procedure was illustrated with different numerical examples of multibody systems: a two-dimensional physical pendulum and wheeled inverted pendulum, a spinning top, a disk rolling without slipping on a horizontal plane and a Watt centrifugal governor model. The approach was also revisited in [220] and applied to study the stability of a bicycle multibody model [221]. Moreover, Xiong et al. [222] linearized the index-3 DAE equations of motion of the Whipple bicycle, moving on a revolution surface, along its circular steady motion. The resulting linearized equations of motion present the same structure obtained with the RCS approach of Gonz´alez et al. [218] and the approach of Escalona et al. [216, 217]. Nevertheless, since these approaches result in a structure of the linearized system of equations that hinders the performance of a standard eigenvalue analysis, a reduced system of first-order ODEs was derived in [222]. All the linearization procedures of this group make use of a redundant set of coordinates. A second group of approaches is based on the linearization of the index-1 form of the DAE system. The index-1 form is obtained by computing the constraints at acceleration level and assembling them with the dynamic equilibrium equations. A straightforward way to carry out the linearization of these nonlinear equations of motion is based on the isolation of the acceleration vector. Once the accelerations are obtained, a linear ODE system can be built. Nevertheless, this approach is impractical for complex multibody systems with large number of bodies, coordinates and constraints, due to the high computational cost required. Van Khang et al. [223] proposed a symbolic linearization procedure using redundant coordinates. The correctness of the approach was verified with a four-bar mechanism with an elastic connecting link. The reduction of the DAE system to an ODE system, by using a coordinate partition in terms of independent and dependent coordinates, gives rise to a third group of linearization approaches. Some works describing the transformation of the nonlinear equations of motion of a constrained dynamical system in DAE form into an ODE system, expressed in terms of the independent coordinates, are the works of Wehage and Haug [224] and Garc´ıa de Jal´on and Bayo [61]. The linearization of this ODE system of equations was performed by Cuadrado et al. [225, 226] to compute the state-space Jacobian matrix, required for building continuous extended Kalman filters (CEKF). In [225, 226, 227], the derivatives of the independent accelerations with respect to the independent coordinates and velocities are required, and the size of the state-space Jacobian matrix is twice the number of degrees of freedom of the multibody system. These derivatives are also necessary in the construction of the state-transition matrix of discrete extended Kalman filters (DEKF) [227] and error-state extended Kalman filters [227, 228, 229]. A similar result to that of [225, 226] was obtained by Gonz´alez et al. [218] with 33 Chapter 1. Introduction the procedure known as Minimal Coordinate Set (MCS) Formulation: Velocity Projection. In this work, a further development of the derivatives of the independent accelerations with respect to the independent positions and velocities is presented. In all these works, the state-space Jacobian matrix includes only the partial derivatives with respect to the independent positions and velocities. Nevertheless, some explanatory comments regarding the role of the dependent coordinates in these linearization approaches are required to compute the correct linearized equations of motion. As stated by Peterson et al. [230], where a linearization procedure was developed for equations of motion generated by Kane’s method [37], the dependent coordinates and dependent velocities are functions of the independent coordinates and independent velocities. Therefore, the dependent states cannot be eliminated in general without considering this dependency with the independent ones. To the best of the authors’ knowledge, no previous author makes this explicitly clear when presenting procedures for symbolically linearizing the equations of motion of constrained dynamical systems [230]. In some particular cases, one can intuitively find a set of independent coordinates that allows for eliminating the derivatives with respect to the dependent ones in the linearization. For example, Escalona et al. [171] analysed the linear stability of the steady forward motion of a bicycle multibody model by resorting to a similar procedure to those of [225, 226], considering only the derivatives with respect to the independent coordinates in the linearization. In [171], it is specified that the dependent coordinates are considered as null in the linearization, since, in the particular case of the bicycle multibody model, they do not have influence on the linear stability. Lastly, a fourth group of linearization approaches can be considered. As in the third group, these procedures make use of a coordinate partition in terms of independent and dependent coordinates. Nevertheless, these approaches achieve a linear ODE system of equations that is only expressed in terms of the independent coordinates. In this case, the influence of the dependent coordinates on the linearization is considered by resorting to the linearized constraints, which allows obtaining an expression of the dependent coordinates in terms of the independent ones. In the case of a multibody system constrained only by holonomic constraints, the Jacobian matrix only includes the derivatives with respect to the independent variables, without ignoring the dependent ones, and presents a size of twice the number of degrees of freedom of the multibody system (as in [225, 226, 227]). Bae [231] presented a linearization procedure in which the perturbed contraint equations are solved for variations of all the variables in terms of the independent ones. Subsequently, these relations are substituted into the variational dynamic equations to obtain the linearized equations of motion in terms of the independent variations. This strategy was also followed by Cossalter et al. [190], where a procedure for evaluating the stability of steady state solutions in straight running and on a curve of sports motorcycles was presented. The linearization approaches of this group achieve the maximum possible reduction of the Jacobian matrix. In addition to the approaches already mentioned, there exist alternative pro- 34 1.3. Linearization in Multibody System Dynamics cedures to linearize the equations of motion of multibody systems. An example is the third procedure presented by Gonz´alez et al. in [218], so-called Unconstrained Coordinate Set (UCS) Formulation. This approach is based on the Penalty method, in which a penalty-based relaxation of the constraints allows transforming a DAE system into a set of ODEs [232]. Moreover, Gonz´alez et al. [233] identified two types of linearization problems: heavily constrained multibody systems, with a number of degrees of freedom much lower than the number of kinematic variables, and systems in which both numbers are similar. Other works are based on the linearization of equations of motion without Lagrange multipliers. Xiong et al. [172] derived the Gibbs–Appell equations [76, 234] for the dynamic analysis of the Whipple bicycle, obtaining a minimal set of equations of motion without Lagrange multipliers. These equations are linearized and the stability analyses of the steady straight-line motion and the circular steady solution are performed. Xiong et al. [235] also applied Voronets equations [70] to the Whipple bicycle moving on a revolution surface, to obtain a seven-dimensional reduced dynamic system on the reduced constraint space. Two types of relative equilibria (static and dynamic) were studied and the stability of the bicycle was discussed. As previously mentioned, a linearization approach based on Kane’s equations of motion [37] was presented by Peterson et al. [230]. Moreover, Saccon et al. [236] proposed, from the standpoint of Geometric Mechanics, a recursive and singularity-free linearization procedure for moving-based robot systems. Additional works address the linearization and include stability analyses of multibody systems as bicycles [169, 170, 176, 178, 177, 179, 180], motorcycles [182, 187, 188, 189, 190, 191, 237, 238, 239] or electric kickscooters [2, 202]. While the theory of modelling constrained multibody systems is well developed and many references in the literature include formalisms to generate the nonlinear equations of motion, the linearization of those equations has not been detailed in a general procedure that could be applied to any kind of holonomic or nonholonomic multibody system, irrespective of size and complexity. As it was shown in this chapter, there are many works devoted to the linearization of multibody systems. Nevertheless, the literature includes examples of linearized equations of motion of multibody systems based on symbolic computation, which is not an option when the multibody system model is sufficiently large in terms of equations and coordinates. Efficiency of the linearization procedures has not been an issue so far since symbolic computation of the linearized equations results in very cumbersome mathematical expressions, which in some cases cannot be handled by today’s personal computers. In this PhD thesis, different numerical linearization approaches for general multibody systems are developed. To validate the proposed procedures, a well-acknowledged benchmark bicycle model has been used. According to the results, the numerical linearization procedures are completely general, more efficient and very accurate as compared to their symbolic counterpart approaches. 35 Chapter 1. Introduction 1.4 General motivation and research objectives The general motivation of this PhD thesis is the development of a methodology that contributes to the modelling of nonholonomic multibody systems, the study of their linear stability and the design of linear feedback controllers. Particular emphasis is placed on single-person nonholonomic vehicles. The computational modelling and stability analysis of these systems are important from the point of view of safety in urban transportation. Understanding how the critical design parameters influence stability will enable the implementation of design modifications that will result in safer vehicles, leading to a decrease in accidents in urban transportation. At the same time, these safer designs will be oriented to a wider range of potential users, as elderly and physically impaired individuals. The development of an analysis methodology for these vehicles will open the door to new driver assistance systems that increase safety while driving. Just as the vast majority of cars today have active safety systems, such as ABS or traction control, it is possible to develop driving assistance systems for single-person vehicles that act on power management or braking systems. From a regulatory point of view, this work will also have a significant impact. Knowledge of the stability conditions of these vehicles will be of great use to experts who are part of regulatory committees at global levels and have the task of drafting regulations for the circulation of these vehicles within cities. It is important to emphasize, as previously mentioned, that there is currently no uniform regulation on the use of these vehicles. From an environmental perspective, the developments aimed at improving safety will increase the number of electric single-person vehicle users. The increasing use of this eco-friendly means of transportation will help reduce the environmental degradation in cities due to the reduction in the number of combustion engine vehicles. In addition, it will promote the achievement of pollution reduction goals set by the world’s major cities. All these design modifications could be developed by means of computationally efficient models and analyses carried out in this study. To this end, the following research objectives are defined: 1. Objective 1. Development of nonholonomic multibody models. The computational modelling and analysis of multibody systems with nonholonomic constraints are considerably more complex than in holonomic multibody systems. Nonholonomic systems are addressed less frequently by researchers, existing a greater lack of knowledge about their characteristics. This PhD thesis aims to increase knowledge on highly mobile nonholonomic multibody systems, with nonlinear constraints in position but linear in velocity, which are obtained when modeling rolling contact. In this study, vehicles’ wheels are modelled as rigid tori rolling without slipping on the ground surface. The classical Newton-Euler formalism, including the holonomic and nonholonomic constraints by means of the method of Lagrange multipliers, is used to systematically obtain the nonlinear equations of motion of nonholonomic multibody systems. The models obtained 36 1.4. General motivation and research objectives through this approach are of particular importance for two reasons. First, it is generally accepted that this approach leads to the correct modelling of nonholonomic systems. On the other hand, the resulting equations of motion constitute a Differential-Algebraic Equations (DAE) system, containing differential equations and algebraic constraint equations. With the DAE system, it can be analysed how the presence of constraints affects the study of stability. 2. Objective 2. Overview of existing linearization approaches in Multibody System Dynamics. The linearization of the equations of motion of nonholonomic multibody systems is of paramount importance in this study. Therefore, a comprehensive overview of the linearization approaches existing in the field of Multibody System Dynamics is required. As already mentioned in the previous section, there exist numerous works devoted to linearization in the context of multibody dynamics. A compilation of the main linearization procedures, illustrated by using a common notation and classified in different groups, will be very useful for the multibody dynamics community. The classification of the approaches could be performed by considering the initial form of the nonlinear equations of motion, the use of a redundant or a minimal set of generalized coordinates, or the structure of the linearized equations of motion, which is highly dependent on the selection of coordinates. 3. Objective 3. Development of novel linearization approaches for constrained multibody systems. Numerous works in the literature include examples of linearized equations of motion of multibody systems based on symbolic computation, which is not feasible for complex multibody systems, with high number of coordinates and long kinematic chains. Obtaining cumbersome nonlinear equations of motion constitutes a great limitation to efficiently derive the symbolic linearized equations, and in many cases, this task cannot be handled by today’s personal computers. General linearization approaches, which could be implemented symbolically or numerically and applied to any kind of holonomic or nonholonomic multibody system, irrespective of size and complexity, need to be detailed. Therefore, this PhD thesis aims to develop novel linearization procedures for general constrained multibody systems, which will be a significant contribution for the multibody dynamics community. By using these approaches, the linearization can be performed along predefined trajectories, which sometimes may be expressed analytically or, in general, numerically. Some examples are the steady forward motion of two-wheeled vehicles or their circular motions. 4. Objective 4. Linear stability analysis of highly mobile nonholonomic multibody systems. As previously mentioned, the main objective of this PhD thesis is the performance of detailed stability analyses of nonholonomic multibody systems. The study of the linear stability of these systems has not yet been approached in a methodological way. The use 37 Chapter 1. Introduction of powerful and efficient linearization approaches allows for generating the Jacobian matrix of the system under study in terms of the main design parameters. This enables obtaining the eigenvalues of the multibody system parameterized in terms of the geometric and dynamic parameters that influence the stability, as the aspect ratios of the toroidal wheels, steering angle, trail, positions of centres of mass, etc. Therefore, by varying these parameters over a wide range of values, a detailed sensitivity analysis can be carried out to understand changes in the vehicle behavior. This will result in improved vehicle safety conditions and a wider range of potential users, as elderly people, whose motor skills and sensory faculties may decline, or users with lower skills for driving. As examples of application, detailed multibody systems of the bicycle, the waveboard or the electric kickscooter are considered in this work. In particular, the linear stability of the steady forward motion of these vehicles is studied. In these cases, linear time-invariant (LTI) systems are obtained, and the study of the linear stability can be directly performed by means of the computation of the eigenvalues of the system. Moreover, the stability analysis of more complex trajectories, as periodic ones, are also addressed in this study. In this case, the linearization leads to periodic linear time-varying (LTV) systems. The linear stability of the rolling toroidal wheel along the circular steady motion is analysed as a case study. Finally, another application included in this study is the use of the novel linearization procedures in multiphysics scenarios, and, in particular, with hydraulically actuated multibody systems. 5. Objective 5. Design of linear feedback controllers for nonholonomic multibody systems. The linearization procedures developed in this document constitute an excellent tool for combining multibody dynamics and linear control theory. Taking advantage of the linearized equations of motion, linear feedback controllers can be designed by using widely extended optimal control methods, such as the Linear-Quadratic Regulator (LQR). As a case study, the optimal control of the ball-plate system, which constitutes a classical nonholonomic benchmark problem in the context of nonprehensile manipulation, is considered. Figure 1.10 shows a schematic representation of the objectives of this PhD thesis and the relationships between them, with the corresponding chapters where they are developed. Objectives 1 and 2 provide the basis for developing Objective 3, which constitutes the main contribution of this document. Once this objective is fulfilled, the novel linearization approaches can be used in several applications. First, the linear stability analysis of nonholonomic multibody systems (Objective 4) is performed. Furthermore, the linearized equations of motions are used to design linear feedback controllers (Objective 5). Objective 1 is developed in Chapter 2; Chapter 3 covers Objective 2; Chapter 4 is devoted to Objective 3; Chapters 5, 6 and 7 develop Objective 4 and Chapter 8 corresponds to Objective 5. 38 1.5. By-chapter summary Figure 1.10: Schematic representation of the main objectives of the PhD thesis, with the corresponding chapter of the document. 1.5 By-chapter summary A by-chapter summary with the work done throughout this PhD thesis is presented. The first chapter corresponds to this Introduction, and the following chapters are: -Chapter 2. Modelling of constrained multibody systems. The main objective of Chapter 2 is the development of nonholonomic multibody models from the broadest possible perspective, with the sole particularity of employing nonholonomic constraints linear in velocities, obtained by modelling rolling without slipping. To this end, the basics of rigid body kinematics and dynamics, using the Multibody System Dynamics formalism and with the aim of computational implementation, are summarized. The main tools, used in the following chapters of this study, are reviewed. First, some basic concepts in multibody dynamics are presented. Next, the most important aspects of kinematics, including some comments of holonomic and nonholonomic constraints, are addressed. The dynamic equilibrium equations of multibody systems, using the Newton-Euler formalism, are derived. Different forms of the nonlinear equations of motion, depending on the generalized coordinates used (a minimal or a redundant set) or the structure of the equations, are outlined. Finally, the inverse and forward dynamic simulations of constrained multibody systems are summarized. 39 Chapter 1. Introduction -Chapter 3. An overview of linearization approaches in Multibody System Dynamics. A comprehensive overview of the linearization approaches for constrained multibody systems existing in the literature is shown. The procedures are classified into four different groups and presented by using a common notation. This provides context to the previous work done in linearization of constrained multibody systems, and justifies the novel work carried out in this PhD thesis regarding linearization. Furthermore, the handling of the dependent coordinates in each of the reviewed approaches is clarified. -Chapter 4. Novel linearization approaches for constrained multibody systems. This chapter presents five novel linearization approaches for constrained multibody systems. All the steps required for the derivation of the linearized equations of motion are shown in detail. Among the five novel approaches, two of them are labeled as symbolic procedures, since the linearized equations can only be computed by means of symbolic computation. The other three approaches allow for symbolic and numerical implementation, and are labeled as numerical procedures. The correctness of the approaches is validated by using the well-acknowledged benchmark bicycle model of Meijaard et al. [168]. According to the results, the numerical linearization procedures are completely general, suitable for constrained multibody systems, irrespective of size and complexity, accurate and more efficient. Finally, the handling of the dependent coordinates in these novel linearization approaches is illustrated with the linear stability results of a bicycle multibody model with toroidal wheels, which constitutes an extension of the bicycle benchmark. -Chapter 5. Linear stability analysis of nonholonomic multibody systems: study of linear time-invariant systems. This chapter is devoted to the main goal of the PhD thesis. Using the linearization approaches developed in Chapter 4, the linearized equations of motion of different nonholonomic multibody systems are obtained. The linear stability analyses of this chapter correspond to relative equilibria that lead to linear-time invariant (LTI) systems. A brief summary of Lyapunov theory for the linear stability analysis of LTI systems is included. In particular, the linear stability of the steady forward motion of the nonholonomic multibody systems is assessed. First, some classical nonholonomic systems, as the skateboard or the rolling disk, are analysed. Next, the linear stability of highly mobile nonholonomic vehicles, as the bicycle, the waveboard or the electric kickscooter, is studied. In the case of a linear time-invariant (LTI) system, the study of the linear stability can be directly performed by means of the computation of the eigenvalues of the system. By sweeping the main design parameters across a wide range of values, their influence on the stability of the nonholonomic multibody systems is studied. -Chapter 6. Linear stability analysis of nonholonomic multibody systems: study of linear periodic time-varying systems. This chapter aims to study the linear stability of periodic linear time-varying (LTV) systems. First, basic aspects of LTV systems are presented. Next, a summary of Floquet- 40 1.5. By-chapter summary Lyapunov theory for the stability analysis of periodic trajectories is included. As example of application, a detailed and comprehensive linear stability analysis of a rolling toroidal wheel is performed. The circular steady motion and the linearized equations of motion along this relative equilibrium are obtained, for both the solid and hollow tori. By virtue of Floquet’s theorem, the periodic time-dependent Jacobian matrix is transformed into a Jacobian matrix with constant coefficients, which allows for the computation of the eigenvalues. These eigenvalues are obtained as a function of the geometric parameters of the rolling toroidal wheel. The change of variable to carry out this transformation can be intuitively found due to the symmetry of the reference solution. Finally, the stability results are validated using the Floquet theory for periodic linear differential equations, by computing the Floquet multipliers of the monodromy matrix of the system. -Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios. The linearization approaches developed in Chapter 4 are applied to multiphysics problems. A particularly interesting case is that of hydraulically actuated multibody systems, since hydraulics is often used to actuate mechanisms in several applications, such as heavy machinery. In this case, the equations of motion of multibody systems with holonomic and nonholonomic constraints are augmented with the differential equations of the hydraulic subsystem. A linearization approach for hydraulically driven multibody systems is presented and tested with a hydraulically actuated four-bar mechanism. A comparison with a conventional counterpart procedure is included, demonstrating the accuracy and efficiency of the proposed approach. -Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems. The linearization approaches proposed in Chapter 4 can be used to design linear feedback controllers. In particular, the third numerical linearization approach is extended to the case of controlled multibody systems. The expression of the input matrix is provided. From the continuos-time description of the linear system, a discrete-time representation is derived. Next, a linear feedback controller, designed by using a widely extended optimal control method as the Linear-Quadratic Regulator (LQR), is developed. As a case study, the optimal control of a ball rolling without slipping on a plate, which constitutes a well-known nonholonomic benchmark problem in the context of nonprehensile manipulation, is considered. -Chapter 9. Conclusions and future lines of research. This chapter highlights the original contributions of this PhD thesis and presents some possible lines of research for future work. 41 Chapter 2. Modelling of constrained multibody systems 2.1 Basic concepts in multibody systems In this section, some basic concepts regarding rigid multibody systems are summarized: the number of degrees of freedom of a multibody system, the types of coordinates that can be used in a multibody model and basic ideas of global and local reference frames, including the orientation matrix of a body. 2.1.1 Degrees of freedom of a multibody system The motion of a body with respect to a global reference frame is mathematically defined by the time evolution of a set of variables. These description variables, called generalized coordinates, must be able to uniquely define the position of all the material points composing the body at any instant of time. In the case of rigid bodies with spatial motion, which are the scope of this work, a minimum of six coordinates is required. Nevertheless, a larger set of coordinates can be used to describe the motion of a rigid body. In this case, these coordinates must satisfy certain relationships, which from a mathematical point of view are constraint equations that are solved together with the dynamic equations of the multibody system. Knowing the number of degrees of freedom is crucial for the modelling and analysis of multibody systems. The number of degrees of freedom of a multibody system can be defined as the minimum number of independent generalized coordinates to uniquely describe the configuration of the system [3, 251]. The kinematic joints interconnecting the bodies of a multibody system introduce different constraints, which must be independent from one another. The number of constraints depends on the number and type of joints existing in the multibody system. Therefore, the number of degrees of freedom, henceforth denoted as ng, can be computed as the difference between the number of coordinates chosen to model the system and the number of independent constraints [252]: ng=n−nc,(2.1) where nis the number of generalized coordinates and ncis the number of position constraints. 2.1.2 Types of coordinates The dynamics of the rigid bodies that form a certain multibody system can be described by different sets of coordinates. A common way to classify generalized coordinates is into dependent and independent coordinates [224]. In some cases, it is possible to use the Lagrange equations of the second kind, using a minimum set of coordinates equal to the number of degrees of freedom [253]. In this scenario, one talks about “independent coordinates”. Nevertheless, the use of this type of coordinates is not systematic and, therefore, the programming of the equations of motion is reduced to some particular cases [27]. Using a larger set of variables to obtain the equations of motion usually leads to less cumbersome equations. 48 2.1. Basic concepts in multibody systems In this case, the variables must satisfy certain relationships that constitute the constraint equations. The set of coordinates used in this approach is referred to as “dependent”, and the dynamic problem is formulated according to Lagrange’s equations of the first kind, which constitutes a system of Differential-Algebraic Equations (DAE). Within the dependent coordinates, a classification in relative coordinates [254], absolute or reference point coordinates [255] and natural coordinates [61] can be made. A discussion with the main types of coordinates used in multibody dynamics, highlighting their advantages and drawbacks, can be found in [3, 28, 61, 256]. The choice of the appropriate type of coordinates to be used depends on the type of problem to be analyzed. To comprehensively illustrate the use of relative, absolute and natural coordinates, a planar multibody model of a four-bar mechanism is used below. Relative coordinates The relative coordinates, also known as joint coordinates, are used to describe the position and orientation of a body relative to the preceding one in a multibody system. These coordinates are typically associated with the relative degrees of freedom permitted by the joints that link the bodies. Relative coordinates can be associated with linear or angular displacements. Figure 2.1: Four-bar mechanism described by using relative coordinates. Consider the planar four-bar mechanism shown in Fig. 2.1. The multibody system presents four rigid bodies: the ground link is body 1; the input link is body 2; the coupler is denoted as body 3, and the output link is body 4. These links are connected to each other by means of revolute joints, located at points O2,C,Dand O4. The following set of relative coordinates can be used to model the four-bar mechanism: x=θ2θ32 θ43 T,(2.2) where θ2represents the rotation of the input link with respect to the ground link; 49 Chapter 2. Modelling of constrained multibody systems θ32 describes the rotation of the coupler with respect to the input link; and θ43 is the rotation of the output link with respect to the coupler. It is well-known that the number of degrees of freedom of a planar four-bar mechanism is one. Therefore, given that n= 3 relative coordinates are used in Eq. (2.2) to describe the system, two constraints relating these coordinates are required. These constraints correspond to the loop equations of the mechanism, and are given by: L2cos (θ2) + L3cos (θ2+θ32) + L4cos (θ2+θ32 +θ43)−L1= 0,(2.3) L2sin (θ2) + L3sin (θ2+θ32) + L4sin (θ2+θ32 +θ43)=0,(2.4) where Li, with i= 1 . . . 4, is the length of link i. Absolute coordinates On the other hand, the absolute or reference point coordinates refer the position and orientation of each body to a global reference frame. Considering an arbitrary body of the multibody system, the location of its center of mass is typically defined. The components of its position vector, with respect to the global reference frame, are the translational absolute coordinates. Furthermore, the orientation of the body is defined by a set of angles or angular parameters, describing the orientation of the body reference frame in the global frame. In the three-dimensional space, the location of a rigid body using absolute coordinates can be defined by using three translational and three rotational coordinates. In the case of planar motion, three absolute coordinates are required to completely describe a given body i: two Cartesian coordinates, xiand yi, and one angle θi, measured in the counterclockwise direction with respect to the X-axis of the global frame. Figure 2.2: Four-bar mechanism described by using absolute coordinates. Consider now the previous four-bar mechanism example. The following set of absolute coordinates can be used to describe the system: x=x2y2θ2x3y3θ3x4y4θ4T,(2.5) 50 2.1. Basic concepts in multibody systems where, in this case, a set of n= 9 coordinates is chosen (see Fig. 2.2). Note that, since the multibody system is planar, only three absolute coordinates are used for each of the moving links. Given that the number of degrees of freedom of the system is one, eight constraints (nc= 8) are required: x2−L2 2cos (θ2) = 0,(2.6) y2−L2 2sin (θ2)=0,(2.7) x2+L2 2cos (θ2)−x3+L3 2cos (θ3)=0,(2.8) y2+L2 2sin (θ2)−y3+L3 2sin (θ3)=0,(2.9) x3+L3 2cos (θ3)−x4−L4 2cos (θ4)=0,(2.10) y3+L3 2sin (θ3)−y4−L4 2sin (θ4)=0,(2.11) x4−L4 2cos (θ4)−L1= 0,(2.12) y4−L4 2sin (θ4)=0.(2.13) The use of absolute coordinates is straightforward, being an appropriate choice to systematically model complex multibody systems. Natural coordinates An alternative to absolute and relative coordinates, so-called natural coordinates, was developed by Garc´ıa de Jal´on et al. [257, 258] during the 1980s for modelling rigid multibody systems. Unlike absolute coordinates, which include angles or angular parameters, natural coordinates present the great advantage that there is no need for any angular variables, being composed by position vectors of material points or vectors that define trajectories in solids, that is, exclusively Cartesian coordinates. The points usually correspond to relevant positions of the multibody system components, such as joints. Given that the use of natural coordinates allows for sharing coordinates between different bodies, a significant reduction in the total number of coordinates required to model a mechanism is achieved. A comparative analysis between the use of absolute and natural coordinates was carried out by Unda et al. [259], concluding that the number of natural coordinates is always lower than the number of absolute coordinates. 51 Chapter 2. Modelling of constrained multibody systems Figure 2.3: Four-bar mechanism described by using natural coordinates. Returning to the example of the four-bar mechanism, the Cartesian components of the position vectors of the revolute joints Cand Dcan be chosen as natural coordinates: x=xCyCxDyDT.(2.14) Three constraint equations are required. The three conditions are obtained by using the concept of rigid link, since the distances between points O2and C;C and D;Dand O4are constant. (xC−xO2)2+ (yC−yO2)2−L2 2= 0,(2.15) (xD−xC)2+ (yD−yC)2−L2 3= 0,(2.16) (xD−xO4)2+ (yD−yO4)2−L2 4= 0.(2.17) In Eqs. (2.15)-(2.17), (xO2, yO2) and (xO4, yO4) are known and define the position of the fixed revolute joints O2and O4, respectively. These three types of coordinates are used throughout this document. Given that the absolute coordinates have the great virtue of being quite straightforward, even for systems with high level of complexity, it is the preferred option for modelling multibody systems in this work. 2.1.3 Global and local reference frames: the orientation matrix of a body Prior to expose the basics of rigid body kinematics, some aspects of global and local reference frames are presented. Given a mechanical system with nbbodies, consider body i, with i= 1 . . . nb, which is part of the multibody system. This body has a local frame ⟨Xi, Yi, Zi⟩, whose origin can be located at any point (it may be even a point that does not belong to the body). The origin of this body frame is denoted as Oi. The global reference frame is denoted as ⟨X1, Y1, Z1⟩, with origin O1. 52 2.1. Basic concepts in multibody systems As previously mentioned, when absolute coordinates are used to describe the body, its position and orientation are defined by means of a set of translational and rotational coordinates. Therefore, the location of the origin of the body frame Oiis defined in the global frame with the vector rOi. Moreover, if a multibody system undergoing planar motion is considered, the orientation of the body-fixed coordinate system is specified by means of the angle of rotation θiof this system of coordinates with respect to the X1-axis of the global frame. The global and body reference frames with their respective origins, the vector rOiand the angle θiare shown in Fig. 2.4 for a general rigid body undergoing planar motion. Figure 2.4: Rigid body iundergoing planar motion. Global reference frame ⟨X1, Y1, Z1⟩with origin O1; body reference frame ⟨Xi, Yi, Zi⟩with origin Oi; position vector rOiand angle of rotation θi. Consider now an arbitrary point Pthat belongs to body i. Let ¯rOiPbe the vector from Oito P, expressed in the body reference frame ⟨Xi, Yi, Zi⟩, and rOiP the same vector expressed in the global frame. The components of these vectors read as follows: ¯rOiP=¯rX OiP¯rY OiPT,(2.18) rOiP=rX OiPrY OiPT.(2.19) In Eq. (2.18), the X- and Y-components of ¯rOiP, denoted as ¯rX OiPand ¯rY OiP, are constant (time-independent), since the vector ¯rOiPis expressed in the body reference frame. From now on, any vector expressed in a local body frame will present a bar, as is the case of ¯rOiP. In contrast, the absolute components of the vector rOiP, denoted as rX OiPand rY OiPin Eq. (2.19), are expressed in the global reference frame, and thus change as the body goes through motion, not being constant. To express the vector ¯rOiP 53 Chapter 2. Modelling of constrained multibody systems in the global frame, the orientation matrix of body i, denoted as Ri, is used as follows: rOiP=Ri¯rOiP.(2.20) If the motion is planar, the orientation matrix Riis a function of the angle of rotation θi, given by: Ri(θi) = cos (θi)−sin (θi) sin (θi) cos (θi).(2.21) To carry out the transformation in reverse, from the global to the body reference frame: ¯rOiP=RT irOiP.(2.22) So far, only the case of a body with planar motion has been considered. In a multibody system undergoing spatial motion, the location of a rigid body ican be defined, by using absolute coordinates, with three translational and three rotational coordinates that describe the origin and orientation of the body reference frame with respect to the global frame. These rotational coordinates are known as Euler angles. The use of three composed elemental rotations, about the axes of a coordinate system, always allows for reaching any orientation. This is also equivalent to say that any rotation can be decomposed as a product of three elemental rotation matrices. There exist twelve possible sequences of rotation axes [260], divided in two groups: -Proper Euler angles. These angles use the same type of axis (X,Yor Z) for both the first and third elemental rotations. The sequences within this group are: X1-Y2-X3,X1-Z2-X3,Y1-X2-Y3,Y1-Z2-Y3,Z1-X2-Z3and Z1-Y2-Z3. -Tait-Bryan angles. Also known as Cardan or nautical angles, they represent rotations about three distinct axes. The following sequences fall within this group: X1-Y2-Z3,Y1-Z2-X3,Z1-X2-Y3,X1-Z2-Y3,Z1-Y2-X3and Y1-X2-Z3. Two examples are shown below, starting from the global reference frame and moving to the body reference frame. From now on, when it is indicated that a rotation is counterclockwise or clockwise, it is considered to be viewed from above. First, among the proper Euler angles, the most common sequence is the Z1-X2-Z3 convention, which is shown in Fig. 2.5 (a). In this sequence, a counterclockwise rotation ψaround the Z1-axis is first performed. Next, a counterclockwise rotation ϕaround the intermediate Xψ-axis is made, and finally the counterclockwise rotation γaround the Zϕ-axis is carried out. The elemental rotation matrices are given by: Rψ=  cψ−sψ0 sψcψ0 0 0 1  ,Rϕ=  1 0 0 0 cϕ−sϕ 0 sϕcϕ ,Rγ=  cγ−sγ0 sγcγ0 0 0 1  , (2.23) 54 2.1. Basic concepts in multibody systems where the abbreviations c∗= cos (∗) and s∗= sin (∗) are introduced for the sake of simplicity. The orientation matrix of body iis given by: Ri(ψ, ϕ, γ) = RψRϕRγ=  cψcγ−cϕsψsγ−cψsγ−cϕcγsψsψsϕ cγsψ+ cψcϕsγcψcϕcγ−sψsγ−cψsϕ sϕsγcγsϕcϕ . (2.24) (a) Proper Euler angles: rotations in the Z1-X2-Z3conven- tion (b) Tait-Bryan angles: rotations in the Z1-X2-Y3conven- tion Figure 2.5: Examples of rotation sequences: proper Euler angles in the Z1-X2-Z3 convention and Tait-Bryan angles in the Z1-X2-Y3convention. Secondly, the Tait-Bryan angles in the Z1-X2-Y3convention are shown in Fig. 2.5 (b). In this case, a counterclockwise rotation ψaround the Z1-axis is first made. Next, a counterclockwise rotation ϕaround the intermediate Xψ-axis is performed, and finally, the counterclockwise rotation θaround the Yϕ-axis is carried out. The elemental rotation matrices are: Rψ=  cψ−sψ0 sψcψ0 0 0 1  ,Rϕ=  1 0 0 0 cϕ−sϕ 0 sϕcϕ ,Rθ=  cθ0 sθ 0 1 0 −sθ0 cθ , (2.25) and the orientation matrix of body iis, in this case: Ri(ψ, ϕ, θ) = RψRϕRθ=  cψcθ−sψsϕsθ−cϕsψcψsθ+ cθsψsϕ cθsψ+ cψsϕsθcψcϕsψsθ−cψcθsϕ −cϕsθsϕcϕcθ . (2.26) To concentrate on the rotational coordinates of the body, the global and body frames have been drawn with the same origin in Figs. 2.5 (a) and (b). The rows and columns of the orientation matrices Ri(ψ, ϕ, γ) in Eq. (2.24) and Ri(ψ, ϕ, θ) in Eq. (2.26) have a geometric interpretation. The columns of 55 Chapter 2. Modelling of constrained multibody systems Rirepresent the X-, Y- and Z-components of the body frame unit vectors in the global reference frame. In contrast, the rows correspond to the X-, Y- and Z-components of the global frame unit vectors, expressed in the body reference frame. Figure 2.6: Rigid body iundergoing spatial motion. Global frame ⟨X1, Y1, Z1⟩ with origin O1; body frame ⟨Xi, Yi, Zi⟩with origin Oiand position vector rOi. Figure 2.6 shows the global frame ⟨X1, Y1, Z1⟩with origin O1; the body frame ⟨Xi, Yi, Zi⟩with origin Oiand the position vector of the origin of the body frame rOi, expressed in the global frame. In contrast to the planar case, the orientation of ⟨Xi, Yi, Zi⟩is not determined with a single rotation (θiin Fig. 2.4), and three rotations are required. The absolute position vector of an arbitrary point Pof body i, expresed in the global frame, is given by: rP=rOi+Ri¯ rOiP.(2.27) As explained in the planar case, the orientation matrix Rican be used to transform the components of a vector, expressed in the body frame, into the global frame by using Eq. (2.20). Equation (2.22) also holds for the reverse transformation, from the global to the body reference frame. 2.2 Kinematics of multibody systems The fundamental aspects of rigid body kinematics are reviewed. Without loss of generality, it is assumed that, in the next chapters of this PhD thesis, the origin of the body reference frame Oicoincides, for convenience, with the center of mass of 56 2.2. Kinematics of multibody systems the body, which is denoted as Gi. Given the set of generalized coordinates x∈Rn, the position of the center of mass of body iis a function of these coordinates: rGi=rGi(x).(2.28) 2.2.1 Velocity vectors: linear and angular velocity The velocity vector of an arbitrary point Pof the body, expressed in the global frame, can be obtained by differentiating with respect to time Eq. (2.27): vP=vGi+˙ Ri¯rGiP,(2.29) where vGiis the velocity of the center of mass Gi, given by vGi=˙rGi, and ˙ Riis the time derivative of the orientation matrix Ri, which is computed as follows: ˙ Ri= n X j=1 ∂Ri ∂xj ˙xj.(2.30) Alternatively, by using the traditional expression relating the velocity of two points in a rigid body, vPis given by: vP=vGi+ωi×rGiP,(2.31) where ωiis the angular velocity of body i, and rGiPis the position vector from Gito P. Note that both vectors are expressed in the global frame. By comparing Eqs. (2.29) and (2.31), the following terms are identified: - Velocity of the origin of the body reference frame (in this case, the center of mass): vGi=˙rGi. - Relative velocity of point Pwith respect to Gi:˙ Ri¯ rGiP=ωi×rGiP. The concept of angular velocity is developed below [250]. The angular velocity is a vector that describes the angular speed at which the orientation of the rigid body is changing and the instantaneous axis about it is rotating. All the points of a rigid body experience the same angular velocity. Consider a vector r, which is fixed in the body frame ⟨Xi, Yi, Zi⟩that comoves with the body. Therefore, the components of r, expressed in the body frame ⟨Xi, Yi, Zi⟩, are constants, with ¯r=¯rX¯rY¯rZT, and therefore: ˙ ¯r=0.(2.32) Nevertheless, the components of r, expressed in the global frame ⟨X1, Y1, Z1⟩, vary with time. Therefore, a vector is necessary to take into account the change in orientation of the local frame ⟨Xi, Yi, Zi⟩with respect to the global frame. As a result, the angular velocity of body i,ωi, verifies the following relation: ˙r=ωi×r.(2.33) 57 Chapter 2. Modelling of constrained multibody systems Cnh (x,˙x) = B(x)˙x=0.(2.69) Computing the time derivative of the holonomic constraints (2.68), one obtains the constraints at velocity level: Cx(x)˙x=0,(2.70) Cnh (x,˙x) = B(x)˙x=0,(2.71) where Cx(x) = ∂C(x) ∂x. To perform a kinematic simulation, it is necessary to prescribe the time evolution of n−m−lgeneralized coordinates, hereinafter referred to as the set of independent coordinates xai. On the other hand, the remaining m+lcoordinates constitute the set of dependent coodinates xdd. The time evolution of the dependent coordinates is subject to the evolution of xai and the fulfilment of the constraints. Therefore, the following coordinate partition at velocity level is considered: ˙x=˙xai ˙xdd T.(2.72) Similarly, the corresponding partition at position level is defined: x=xai xdd T.(2.73) The time derivative of the holonomic constraints (2.70), as well as the nonholonomic (2.71), can be written in terms of the dependent and independent coordinates to solve the velocity problem: Cxdd (x)˙xdd +Cxai (x)˙xai =0,(2.74) Bdd (x)˙xdd +Bai (x)˙xai =0.(2.75) Therefore, Cxdd (x) Bdd (x) | {z } Ddd(x) ˙xdd =−Cxai (x) Bai (x) | {z } Dai(x) ˙xai →Ddd (x)˙xdd =−Dai (x)˙xai. (2.76) Since the independent coordinates xai are known time functions in the kinematic simulation, i.e., xai =xai (t), the following first-order system of differential equations can be obtained by isolating the dependent velocities in Eq. (2.76): ˙xdd =−(Ddd (xdd,xai (t)))−1Dai (xdd,xai (t)) ˙xai (t).(2.77) Note that the choice of the independent and dependent coordinates is subject to the existence of the inverse matrix D−1 dd . Given an initial condition xdd (t0) that satisfies the constraints, the system of differential equations (2.77) can be integrated along time by using any numerical algorithm for first-order ODE systems. Nevertheless, the numerical integration might bring along some errors in 64 2.2. Kinematics of multibody systems the fulfilment of position constraints. To avoid this, one can resort to Baumgarte’s stabilization method of constraints [261, 262], by imposing the relation: ˙ C(x) + αeC(x) = 0,(2.78) where αeis a positive constant. Therefore, the velocity problem becomes: ˙xdd =−Cxdd (xdd,xai) Bdd (xdd,xai)−1Cxai (xdd,xai)˙xai +αeC(xdd,xai) Bai (xdd,xai)˙xai , (2.79) with the initial condition xdd0=xdd (t0).(2.80) Note that, in Eq. (2.79), xai =xai (t) is a known function of time. Once the integration of the first ODE system (2.79) is performed, the time evolutions of the dependent coordinates xdd (t) are known. This result can be used in Eq. (2.77) to directly compute the time evolutions of the dependent velocities ˙xdd (t). Finally, it would be necessary to compute the time evolution of the dependent accelerations ¨xdd (t). To this end, the following steps are followed. The time derivative of the holonomic constraints (2.70) and the nonholonomic constraints (2.71) can be assembled to express the velocity constraints as follows: D(x)˙x=0,(2.81) where D(x) is a (m+l)×nmatrix, given by: D(x) = Cx(x) B(x).(2.82) Note that this matrix was already used in Eq. (2.76) in the velocity problem. Differentiating with respect to time Eq. (2.81), the constraints at acceleration level are obtained: D(x)¨x=Qd(x,˙x),(2.83) where Qd(x,˙x) = −∂(D(x)˙x) ∂x˙x. Equation (2.83) can be written in terms of the dependent and independent coordinates to solve the acceleration problem: Ddd (x)¨xdd =Qd(x,˙x)−Dai (x)¨xai.(2.84) The dependent accelerations can be obtained from Eq. (2.84), since the time evolution of all the generalized coordinates at position and velocity level, x(t) and ˙x(t), are already known, and the independent accelerations ¨xai (t) are known functions of time: ¨xdd = (Ddd (x(t)))−1(Qd(x(t),˙x(t)) −Dai (x(t)) ¨xai (t)) .(2.85) 65 Chapter 2. Modelling of constrained multibody systems 2.3 Dynamics of multibody systems As in Sect. 2.2 with the kinematics, the main dynamics aspects of a rigid body are presented first, and then extended to the case of a multibody system with nb bodies from a computational perspective. 2.3.1 Dynamics of a rigid body: Newton-Euler equations Consider a generic body iof the multibody system, with mass mi, center of mass Gi, and subject to the action of NFforces and NMtorques. The traslational motion of the body is governed by Newton’s second law: NF X j=1 Fj=miaGi,(2.86) where Fjis each of the external forces acting on the body. On the other hand, the rotational motion is given by Euler’s second law: NM X k=1 Mk+ NF X j=1 rj×Fj=dLGi dt,(2.87) where Mkis each of the external torques acting on the body; the terms rj×Fj correspond to the torques generated by the forces Fj, exerted on points different to the center of mass of the body; and LGiis the angular momentum of the body with respect to the center of mass, expressed in the global reference frame. To obtain the angular momentum LGi, the inertia tensor of the body IGi, also expressed in the global frame, is used: LGi=IGiωi.(2.88) Nevertheless, it is convenient to apply Euler’s second law using the local reference frame, since the inertia tensor, observed in the body reference frame (and denoted as ¯ IGi), does not change as the body rotates. The angular momentum ¯ LGi, expressed in the local reference frame, is given by: ¯ LGi=¯ IGi¯ωi.(2.89) The absolute time derivative of the angular momentum ¯ LGican be computed as: d¯ LGi dt=d¯ LGi dtrot +¯ωiׯ LGi,(2.90) where d¯ LGi dtrot represents the time derivative of ¯ LGiin the rotating basis, corresponding to the “regular” time derivative acting on the components of ¯ LGi. 66 2.3. Dynamics of multibody systems Note that the additional term ¯ωiׯ LGiin Eq. (2.90) arises from the rotation of the local reference frame ⟨Xi, Yi, Zi⟩. The time derivative d¯ LGi dtrot is given by: d¯ LGi dtrot =d¯ IGi dt¯ωi+¯ IGi d¯ωi dt.(2.91) Since d¯ IGi dt=0and ¯αi=d¯ωi dt, Eq. (2.91) becomes: d¯ LGi dtrot =¯ IGi¯αi.(2.92) Substituting Eq. (2.89) and (2.92) in Eq. (2.90), yields: d¯ LGi dt=¯ IGi¯αi+¯ωi×¯ IGi¯ωi.(2.93) The assembly of Eqs. (2.86) and (2.87), and the use of Eq. (2.93), provide the Newton-Euler equations of a rigid body: mi03 03¯ IGi aGi ¯αi=Fi ¯ Mi+03×1 −¯ωi×¯ IGi¯ωi,(2.94) where mi=miI3,I3is the 3 ×3 identity matrix, 03is a 3 ×3 null matrix, and 03×1is a 3 ×1 null matrix. In Eq. (2.94), Fiis the resultant force acting on body i, expressed in the global reference frame, and ¯ Miis the resultant torque, expressed in the body reference frame. Fiand ¯ Mican be computed as: Fi= NF X j=1 Fj,(2.95) ¯ Mi=RT iMi,(2.96) where Miis the resultant torque, expressed in the global reference frame, given by: Mi= NM X k=1 Mk+ NF X j=1 rj×Fj.(2.97) 2.3.2 Dynamic equations of multibody systems Consider now a multibody system with nbbodies. The Newton-Euler equations (2.94) of all the bodies can be assembled, leading to: ˆ Mˆa=ˆ Qa+ˆ Qv,(2.98) 67 Chapter 2. Modelling of constrained multibody systems where ˆawas defined in Eq. (2.63), and ˆ M,ˆ Qaand ˆ Qvare given by: ˆ M=                m203··· 030303··· 03 03m3··· 030303··· 03 . . .. . ..... . .. . .. . .. . .. . . 0303··· mnb0303··· 03 0303··· 03¯ IG203··· 03 0303··· 0303¯ IG3··· 03 . . .. . .. . .. . .. . .. . ..... . . 0303··· 030303··· ¯ IGnb                ,(2.99) ˆ Qa=           F2 . . . Fnb ¯ M2 . . . ¯ Mnb           ,ˆ Qv=            03×1 . . . 03×1 −¯ω2×¯ IG2¯ω2 . . . −¯ωnb×¯ IGnb¯ωnb            .(2.100) Substituting Eq. (2.65) in Eq. (2.98), one obtains: ˆ M(L¨x+l˙x) = ˆ Qa+ˆ Qv.(2.101) Equation (2.101) can be premultiplied by LT(matrix Lwas defined in Eq. (2.64)). Rearranging terms, the dynamic equations of the multibody system can be written as: M(x)¨x=Qa(x,˙x) + Qv(x,˙x),(2.102) where M=LTˆ ML,(2.103) Qa=LTˆ Qa,(2.104) Qv=LTˆ Qv−ˆ Ml ˙x.(2.105) In Eq. (2.102), the following terms can be identified: -M(x)¨xis the generalized inertia term, being Mthe n×nmass matrix. Note that, although the matrix ˆ Mis not a function of the generalized coordinates, the generalized mass matrix Mis a function of x. -Qvis the generalized quadratic-velocity force vector, associated with the inertia forces that are quadratic with respect to the system velocities (centrifugal and Coriolis forces). -Qais the generalized force vector associated with the external forces and torques acting on the multibody system. 68 2.3. Dynamics of multibody systems In Eq. (2.102), the generalized force vector arising from the reaction forces, associated with the constraints that the multibody system satisfies during its motion, must be included. This vector, which is denoted as Qreac, can be obtained by using the method of Lagrange multipliers. Considering a multibody system with mholonomic and lnonoholonomic constraints, Qreac is given by: Qreac =−DT(x)Λ,(2.106) where Λ∈Rm+lis the (m+l)×1 vector of Lagrange multipliers, and D(x) was defined in Eq. (2.82). The Lagrange multipliers represent the reaction forces associated with the constraints, with opposite sign. Including the generalized force vector given by Eq. (2.106) in Eq. (2.102), the dynamic equilibrium equations of a constrained multibody system are obtained: M(x)¨x+DT(x)Λ=Qa(x,˙x) + Qv(x,˙x).(2.107) The dynamic equations (2.107) must be solved together with the constraints of the multibody system. 2.3.3 Equations of motion of constrained multibody systems The equations of motion of a multibody system, modeled with ncoordinates, m holonomic and lnonholonomic constraints, are given by the dynamic equilibrium equations (2.107), augmented with the holonomic constraints (2.66) and nonholonomic constraints (2.67). The equations of motion constitute a nonlinear index-3 Differential-Algebraic Equations (DAE) system: M(x)¨x+DT(x)Λ=Q(x,˙x),(2.108) C(x) = 0,(2.109) Cnh (x,˙x) = B(x)˙x=0.(2.110) In Eq. (2.108), the vector of generalized forces Q(x,˙x) = Qa(x,˙x) + Qv(x,˙x) has been defined for the sake of simplicity. All the terms in Eq. (2.108), (2.109) and (2.110) have been presented in the previous sections. Once the equations of motion of a constrained multibody system are obtained, some comments about the inverse and forward dynamic analysis are included. 2.3.4 Inverse dynamics analysis Inverse dynamic analysis constitutes an extension of the kinematic simulation presented in Sect. 2.2, but it also allows for obtaining the forces involved in the desired motion of the multibody system. Consider a multibody system with holonomic and nonholonomic constraints, n−m−ldegrees of freedom and subject to n−m−lexternal actuations. The 69 Chapter 2. Modelling of constrained multibody systems (n−m−l)×1 vector of actuations is denoted as τ. The dynamic equilibrium equations (2.108), including these external actuations, can be written as: M(x)¨x+DT(x)Λ=Q(x,˙x) + Aτ(x)τ,(2.111) where Aτ(x) is a n×(n−m−l) projection matrix. The kinematic simulation provides the time evolution of the generalized coordinates x(t), velocities ˙x(t) and accelerations ¨x(t). In Eq. (2.111), the m+lLagrange multipliers Λ(each of them associated with a constraint) and the n−m−l actuations τcan be obtained as follows: −DT(x(t)),Aτ(x(t)) | {z } Aan×n Λ τ | {z } Xan×1 =M(x(t))¨x(t)−Q(x(t),˙x(t)) | {z } ban×1 .(2.112) Note that, to obtain Λand τ, a system of equations of the form AaXa=baneeds to be solved. It is important to remember that the Lagrange multipliers are equal to the required reaction forces for satisfying the constraints, with opposite sign. 2.3.5 Forward dynamics analysis Forward dynamic analysis corresponds to the true dynamic simulation, and allows for knowing the motion of the multibody system from the applied forces by integrating the equations of motion. To perform a forward dynamic simulation, the index-3 DAE system (2.108)-(2.110) is transformed into the index-1 form. First, the time derivative of the holonomic constraints and the nonholonomic constraints are assembled, leading to the velocity constraints (2.81). The following index-2 DAE system is obtained: M(x)¨x+DT(x)Λ=Q(x,˙x),(2.113) D(x)˙x=0.(2.114) Next, the velocity constraints (2.114) can be differentiated once with respect to time, leading to the index-1 DAE system: M(x)¨x+DT(x)Λ=Q(x,˙x),(2.115) D(x)¨x=Qd(x,˙x),(2.116) where the definition of Qd(x,˙x), which was already presented in Eq. (2.83), is repeated here for the reader’s convenience. Qd(x,˙x) = −∂(D(x)˙x) ∂x˙x.(2.117) Note that working with the constraints at acceleration level of Eq. (2.116) may result in the violation of the holonomic and nonholonomic constraints. Therefore, 70 2.4. Forms of the equations of motion Baumgarte’s stabilization method [261] can be used to avoid numerical drift. Equation (2.116) becomes: D(x)¨x=Qd(x,˙x)−2βe˙ C(x) Cnh (x,˙x) | {z } Q1(x,˙x) −γ2 eC(x) 0 | {z } Q2(x) ,(2.118) where βeand γeare two Baumgarte stabilization positive constants, and the vectors Q1(x,˙x) and Q2(x) have been defined for convenience. Assembling the dynamic equations (2.115) and Eq. (2.118), one obtains: M(x)DT(x) D(x)0(m+l) ¨x Λ=Q(x,˙x) Qd(x,˙x)−2βeQ1(x,˙x)−γ2 eQ2(x), (2.119) where 0(m+l)is a null square matrix of size m+l. Given the initial conditions at position and velocity level x0=x(t0) and ˙x0=˙x(t0), the index-1 DAE system (2.119) can be numerically integrated to obtain the time evolution of the generalized coordinates and the Lagrange multipliers. In the particular case of a multibody system with only holonomic constraints, Eq. (2.119) is simplified to: M(x)CT x(x) Cx(x)0m ¨x Λ=Q(x,˙x) Qd(x,˙x)−2βe˙ C(x)−γ2 eC(x). (2.120) Finally, in the case of a nonholonomic multibody system without holonomic constraints, Eq. (2.119) becomes: M(x)DT(x) D(x)0l ¨x Λ=Q(x,˙x) Qd(x,˙x)−2βeCnh (x,˙x).(2.121) In addition to Baumgarte’s stabilization method, there exist other stabilization techniques developed to control constraint drift. A review of these techniques can be found in the work of Bauchau et al. [263]. Among them, Penalty based stabilization techniques can be highlighted, being an example the Augmented Lagrangian formulation developed by Bayo et al. [232]. 2.4 Forms of the equations of motion This section presents a summary with different forms of the nonlinear equations of motion of multibody systems. The equations of motion can be classified according to their structure and the coordinates they use (a redundant or a minimal set). 2.4.1 Index-3 and index-1 DAE form So far, the index-3 and index-1 DAE forms of the equations of motion of a multibody system with ncoordinates, mholonomic and lnonholonomic constraints 71 Chapter 2. Modelling of constrained multibody systems have been presented. Both DAE systems of equations are shown again for convenience, since they are used immediately below to derive additional forms of the equations of motion. The index-3 DAE system is given by: M(x)¨x+DT(x)Λ=Q(x,˙x),(2.122) C(x) = 0,(2.123) Cnh (x,˙x) = B(x)˙x=0.(2.124) The index-1 DAE form is: M(x)¨x+DT(x)Λ=Q(x,˙x),(2.125) D(x)¨x=Qd(x,˙x),(2.126) or, alternatively: M(x)DT(x) D(x)0(m+l) ¨x Λ=Q(x,˙x) Qd(x,˙x),(2.127) with Qd(x,˙x) = −∂(D(x)˙x) ∂x˙x. 2.4.2 ODE system with coordinate partition in terms of independent and dependent coordinates The equations of motion (2.122)-(2.124) can be reduced by means of a coordinate partition at velocity level, based on n−m−lindependent velocities ˙xai and m+l dependent velocities ˙xdd: ˙x=˙xai ˙xdd T.(2.128) This coordinate partition was already presented to perform the kinematic simulation in Sect. 2.2. Similarly, the corresponding partition at position level is defined: x=xai xdd T.(2.129) By using this coordinate partition in the velocity constraints (2.81) (D(x)˙x=0), the following velocity transformation is obtained: ˙x=T(x)˙xai,(2.130) where Tis a transformation matrix, given by: T(x) = I(n−m−l) −(Ddd (x))−1Dai (x).(2.131) In Eq. (2.131), I(n−m−l)is the identity matrix of size n−m−l;Ddd (x) is a (m+l)-square matrix, formed by the columns of matrix D(x) corresponding to the dependent coordinates; and Dai (x)isa(m+l)×(n−m−l) matrix, defined 72 2.4. Forms of the equations of motion by the columns of D(x) associated with the independent coordinates. Note that T(x) is the null-space of D(x), and therefore: TT(x)DT(x) = 0.(2.132) Differentiating with respect to time the velocity transformation (2.130), the following acceleration transformation is obtained: ¨x=T(x)¨xai +˙ T(x)˙xai,(2.133) with ˙ T(x)˙xai =∂(T˙xai) ∂x˙x. Premultiplying the dynamic equations (2.122) by TT(x), and using the velocity and aceleration transformations of Eqs. (2.130) and (2.133), respectively, in the dynamic equations (2.122), the following reduced equations of motion in ODE form are obtained: m(xai,xdd)¨xai =f(xai,xdd,˙xai),(2.134) ˙xdd =H(xai,xdd)˙xai,(2.135) with m(xai,xdd) = TTMT , f(xai,xdd,˙xai) = TTQ−M˙ T˙xai, H(xai,xdd) = −(Ddd)−1Dai. (2.136) Note that the Lagrange multipliers Λof Eq. (2.122) are eliminated in Eq. (2.134) by virtue of Eq. (2.132). An important clarification regarding the ODE equations of motion (2.134)- (2.135) must be made. Note that, despite expressing the complete set of velocities and accelerations as a function of the independent velocities and accelerations by means of Eqs. (2.130) and (2.133), the equations of motion (2.134)-(2.135) are in general a function of the dependent coordinates xdd. In general, this dependency at position level with the dependent coordinates cannot be eliminated, since the holonomic constraints are nonlinear algebraic equations that hinder obtaining an explicit expression of the dependent coordinates as a function of the independent ones. Moreover, the nonholonomic constraints are non-integrable first-order differential equations, which also impedes obtaining the dependent coordinates as a function of the independent ones. Therefore, the first-order differential equations (2.135) are required. 2.4.3 Equations of motion in DAE form without Lagrange multipliers Due to the mnonlinear holonomic constraints (2.123) in the DAE system (2.122)- (2.124), the n-coordinates vector can be partitioned in n−madmissible coordinates xaand mdependent coordinates xd, and therefore x=xaxdT. 73 Chapter 3. An overview of linearization approaches in Multibody System Dynamics Defining ˜ X=˜x˙ ˜x˜ ΛT, the linear DAE system of equations (3.5)-(3.6) can be expressed as the first-order ODE system: A0˙ ˜ X=B0˜ X,(3.7) where A0and B0are the following (2n+m)×(2n+m) matrices: A0=   In0n0n×m 0nMx00n×m 0m×n0m×n0m   ,(3.8) B0=       0nIn0n×m B21 0 ∂Q(x,˙x) ∂˙x0−CT xx0 −Cxx00m×n0m        ,(3.9) with B21 0=∂Q(x,˙x) ∂x0−∂M(x)¨x0 ∂x0−∂CT x(x)Λ0 ∂x0 . Matrices A0and B0are (2n+m+l)×(2n+m+l) if the multibody system presents nonholonomic constraints. As can be seen from the structure of matrix A0in Eq. (3.7), this matrix is singular, since the last mrows are null. Therefore, as stated by Escalona et al. [216, 217] and Gonz´alez et al. [218], the system of equations (3.7) cannot be solved with a standard eigenvalue analysis. Considering a solution of the form ˜ X(t) = X0eλt, where X0is a constant vector, the generalized eigenvalue problem associated with the matrices A0and B0is: (B0−λA0)X0=0.(3.10) The generalized eigenvalue analysis of the 2n+mlinearized equations (3.7) leads to 3mspurious eigenvalues [216, 217, 218], among which 2mare related to the mconstrained dependent variables, and mare associated with the mLagrange multipliers. The resulting 2 (n−m) eigenvalues correspond to the spectrum of the constrained problem. In the particular case that the multibody system includes l nonholonomic constraints, a total of 2n+m+leigenvalues are obtained, among which 2 (n−m−l) correspond to the spectrum of the constrained problem and 3 (m+l) are spurious eigenvalues. This linearization approach makes use of the redundant set of coordinates. As can be seen in Eq. (3.7), the variations of the complete set of coordinates ˜x and their time derivatives ˙ ˜x,¨ ˜xappear. In this case, no coordinate partition in terms of independent and dependent coordinates is made. Therefore, there is no ambiguity in the handling of the dependent coordinates in this approach. Xiong et al. [222] follow a similar approach to linearize the index-3 DAE equations of motion of a Whipple bicycle multibody model, moving on a revolution 80 3.1. Group 1: procedures based on the linearization of the index-3 DAE system, using a redundant set of coordinates surface, along its circular steady motion. The resulting linearized equations of motion present the same structure obtained with the RCS approach of Gonz´alez et al. [218] and Escalona et al. [216, 217]. Nevertheless, since these linearized systems of equations do not allow the performance of a standard eigenvalue analysis, a reduced system of first-order ODEs is derived in [222]. Procedure by Pappalardo et al. [219, 220] The procedure presented by Pappalardo et al. [219, 220] performs the linearization of multibody systems with holonomic and nonholonomic constraints. Computing the Taylor expansion of the nonlinear index-3 DAE system of equations (2.122)-(2.124), retaining up to first-order terms and using Eqs. (3.1)- (3.3), yields: Mx0¨ ˜x+∂M(x)¨x0 ∂x0 ˜x+DTx0˜ Λ+∂DT(x)Λ0 ∂x0 ˜x= ∂Q(x,˙x) ∂x0 ˜x+∂Q(x,˙x) ∂˙x0 ˙ ˜x, ∂C(x) ∂x0 ˜x=0,(3.11) Bx0˙ ˜x+∂B(x)˙x0 ∂x0 ˜x=0. The linear DAE system of equations (3.11) can be expressed as: Mx0¨ ˜x+∂M(x)¨x0 ∂x0 ˜x+DTx0˜ Λ+∂DT(x)Λ0 ∂x0 ˜x= ∂Q(x,˙x) ∂x0 ˜x+∂Q(x,˙x) ∂˙x0 ˙ ˜x,(3.12) Ex|0˜x+E˙x|0˙ ˜x=0,(3.13) where Ex|0and E˙x|0are (m+l)×nmatrices, given by: Ex|0=     ∂C(x) ∂x0 ∂B(x)˙x0 ∂x0      ,E˙x|0=0m×n Bx0.(3.14) In the absence of nonholonomic constraints, Ex|0=∂C(x) ∂x0 ,E˙x|0=0m×n and Dx0=Cxx0. Defining ˜p=˜x˜ ΛT, the linear system given by Eqs. (3.12) and (3.13) can be rewritten as: ¯ M0¨ ˜p+¯ C0˙ ˜p+¯ K0˜p=0,(3.15) 81 Chapter 3. An overview of linearization approaches in Multibody System Dynamics where ¯ M0,¯ C0and ¯ K0are the n+m+lcomposite mass, damping and stiffness matrices, respectively, given by: ¯ M0=Mx00n×(m+l) 0(m+l)×n0(m+l),(3.16) ¯ C0= −∂Q(x,˙x) ∂˙x0 0n×(m+l) E˙x|00(m+l) ,(3.17) ¯ K0=  −∂Q(x,˙x) ∂x0 +∂M(x)¨x0 ∂x0 +∂DT(x)Λ0 ∂x0 DTx0 Ex|00(m+l)   . (3.18) As stated in [219, 220], to carry out the stability analysis of the reference solution, a state-space reformulation of the linearized equations of motion (3.15) is necessary. To this end, by defining ˜ X=˜p˙ ˜pT, Eq. (3.15) can be expressed as the following first-order ODE system: C0˙ ˜ X=D0˜ X,(3.19) where the matrices C0and D0are given by: C0=¯ C0¯ M0 ¯ M00(n+m+l),D0=−¯ K00(n+m+l) 0(n+m+l)¯ M0.(3.20) As the matrix A0in Eq. (3.7), the matrix C0of Eq. (3.19) is singular, since the last m+lrows of ¯ M0are zero. Therefore, the system of equations (3.19) cannot be solved with a standard eigenvalue analysis. Considering a solution of the form ˜ X(t) = X0eλt, where X0is a constant vector, the generalized eigenvalue problem associated with the matrices C0and D0is: (D0−λC0)X0=0.(3.21) The solution of Eq. (3.21) provides 2 (n+m+l) eigenvalues. Given that the spectrum of the constrained problem consists of 2 (n−m−l) eigenvalues, this approach results in 4 (m+l) spurious eigenvalues. Among these spurious eigenvalues, 2 (m+l) are associated with the m+ldependent coordinates of the multibody model, and the remaining 2 (m+l) are related to the m+lLagrange multipliers. As in the linearization approach of Escalona et al. [216, 217] and Gonz´alez et al. [218], the linearized equations of motion are expressed in terms of the variations of the redundant set of coordinates ˜xand their time derivatives ˙ ˜x,¨ ˜x. There is no ambiguity in the handling of the dependent coordinates in this approach, since no coordinate partition in terms of independent and dependent coordinates is made. 82 3.1. Group 1: procedures based on the linearization of the index-3 DAE system, using a redundant set of coordinates Procedure by Negrut et al. [215] Negrut et al. [215] also present a linearization approach for the index-3 DAE equations of motion. In this procedure, the governing equations are augmented and a set of sensitivities of interest that provide the linearization are computed. The approach presented in [215] is developed for heterogeneous multibody systems, including nonholonomic constraints, friction, flexible bodies or user-defined differential equations. Nevertheless, as in [215], it is assumed that the set of constraints are holonomic to keep the presentation of the approach simpler. Moreover, in [215], the vector of generalized velocities, denoted as vin this document, is not assumed to be the direct time derivative of the vector of position coordinates x, and therefore: v=T(x)˙x,(3.22) where T(x) is a transformation matrix. In this work, to present the linearization procedure in a clearer way (without losing any relevant information), it is assumed that T=In, and therefore v=˙x. Negrut et al. [215] present the augmented equations of motion of the multibody system as follows: M(x)¨x+CT x(x)Λ=Q(x,˙x), ¨ C(x,˙x,¨x) = 0, ˙ C(x,˙x) = 0,(3.23) C(x) = 0, B1˙x=˙xai, B0x=xai, where the constraints at velocity and acceleration levels are included. In Eqs. (3.23), xai and ˙xai are n−mindependent positions and velocities, respectively, arising from a coordinate partition of the complete set of coordinates xin dependent and independent coordinates. Moreover, B1and B0are (n−m)×nBoolean matrices. The nonlinear equations of motion (3.23) can be written as: F(x,˙x,¨x,Λ) = b(xai,˙xai),(3.24) where F(x,˙x,¨x,Λ) corresponds to the left-hand side of Eqs. (3.23), and b(˙xai,xai) to the right-hand side. The main objective of this approach is the computation of several first-order sensitivities. In particular, the partial derivatives of the complete set of accelerations and velocities with respect to the independent velocities and positions are of interest in [215]: ∂¨x ∂˙xai ,∂¨x ∂xai ,∂˙x ∂˙xai and ∂˙x ∂xai . To this end, the partial derivatives 83 Chapter 3. An overview of linearization approaches in Multibody System Dynamics ∂F ∂˙xai and ∂F ∂xai are computed: ∂F ∂˙xai =                             M(x)∂¨x ∂˙xai +CT x(x)∂Λ ∂˙xai +Υ˙x ∂˙x ∂˙xai +Υx ∂x ∂˙xai ¨ C¨x ∂¨x ∂˙xai +¨ C˙x ∂˙x ∂˙xai +¨ Cx ∂x ∂˙xai ˙ C˙x ∂˙x ∂˙xai +˙ Cx ∂x ∂˙xai Cx ∂x ∂˙xai B1 ∂˙x ∂˙xai B0 ∂x ∂˙xai                             ,(3.25) ∂F ∂xai =                             M(x)∂¨x ∂xai +CT x(x)∂Λ ∂xai +Υ˙x ∂˙x ∂xai +Υx ∂x ∂xai ¨ C¨x ∂¨x ∂xai +¨ C˙x ∂˙x ∂xai +¨ Cx ∂x ∂xai ˙ C˙x ∂˙x ∂xai +˙ Cx ∂x ∂xai Cx ∂x ∂xai B1 ∂˙x ∂xai B0 ∂x ∂xai                             ,(3.26) where Υ˙xand Υxare given by: Υ˙x=−∂Q(x,˙x) ∂˙x,Υx=∂M(x)¨x+CT x(x)Λ−Q(x,˙x) ∂x.(3.27) 84 3.1. Group 1: procedures based on the linearization of the index-3 DAE system, using a redundant set of coordinates These partial derivatives verify: ∂F ∂˙xai =         0n×(n−m) 0m×(n−m) 0m×(n−m) 0m×(n−m) I(n−m) 0(n−m)         ,(3.28) ∂F ∂xai =         0n×(n−m) 0m×(n−m) 0m×(n−m) 0m×(n−m) 0(n−m) I(n−m)         .(3.29) Equations (3.28) and (3.29) can be assembled to build the following linear system: A1(x)X1=b1,(3.30) where A1(x) =            M(x)CT x(x)Υ˙xΥx ¨ C¨x0m¨ C˙x¨ Cx 0m×n0m˙ C˙x˙ Cx 0m×n0m0m×nCx 0(n−m)×n0(n−m)×mB10(n−m)×n 0(n−m)×n0(n−m)×m0(n−m)×nB0            ,(3.31) X1=                 ∂¨x ∂˙xai ∂¨x ∂xai ∂Λ ∂˙xai ∂Λ ∂xai ∂˙x ∂˙xai ∂˙x ∂xai ∂x ∂˙xai ∂x ∂xai                 ,b1=         0n×(n−m)0n×(n−m) 0m×(n−m)0m×(n−m) 0m×(n−m)0m×(n−m) 0m×(n−m)0m×(n−m) I(n−m)0(n−m) 0(n−m)I(n−m)         .(3.32) In Eq. (3.30), the assembly of Eqs. (3.28) and (3.29) is possible since both equations have the matrix A1(x) in common. Note that the linear system of equations (3.30) provide the sensitivities of the velocities and accelerations with respect to the independent coordinates and velocities. 85 Chapter 3. An overview of linearization approaches in Multibody System Dynamics 3.2 Group 2: procedures based on the linearization of the index-1 DAE system, using a redundant set of coordinates The procedures summarized in this section are based on the linearization of the index-1 DAE system (2.127). Symbolic linearization of the index-1 DAE system A straightforward way to carry out the linearization of the DAE system (2.127) is based on obtaining, first, the acceleration vector ¨xfrom Eq. (2.127): ¨x Λ=M(x)DT(x) D(x)0(m+l)−1Q(x,˙x) Qd(x,˙x)=F(x,˙x),(3.33) where F(x,˙x)isa(n+m+l)×1 vector. Given that the symbolic computation of the inverse matrix in Eq. (3.33) requires a high computational cost, this procedure is not valid for moderately complex multibody systems. Denoting the first ncomponents of F(x,˙x) as F1(x,˙x), and the last m+l components by F2(x,˙x), Eq. (3.33) is rewritten as: ¨x=F1(x,˙x),(3.34) Λ=F2(x,˙x).(3.35) Note that Eq. (3.35) corresponds to the expressions of the Lagrange multipliers Λ, which do not provide information about the stability of the multibody system. In the particular case that the mass matrix M(x) is not singular, an explicit expression can be found for F1(x,˙x) and F2(x,˙x): F1(x,˙x) = M−1Q+DTDM−1DTQd−DM−1Q,(3.36) F2(x,˙x) = −DM−1DT−1Qd−DM−1Q.(3.37) Performing the Taylor expansion of Eq. (3.34), retaining up to first-order terms and using that the reference motion verifies ¨x0=F1x0,˙x0, yield: ¨ ˜x=∂F1 ∂x0 ˜x+∂F1 ∂˙x0 ˙ ˜x.(3.38) Defining ˜ X=˜x˙ ˜xT, the following first-order linear ODE system is obtained: ˙ ˜ X=J˜ X,with J=  0nIn ∂F1 ∂x0 ∂F1 ∂˙x0  .(3.39) The size of Jis, in this case, 2n×2n, obtaining 2neigenvalues. Given that the real spectrum consists of 2 (n−m−l) eigenvalues, associated with the ng=n−m−l 86 3.2. Group 2: procedures based on the linearization of the index-1 DAE system, using a redundant set of coordinates degrees of freedom of the multibody system, this approach leads to 2 (m+l) spurious eigenvalues, arising from the m+ldependent coordinates. Note that, as in the approaches based on the linearization of the index-3 DAE system (Sect. 3.1), the Jacobian matrix (3.39) includes the partial derivatives with respect to the redundant set of coordinates x. Therefore, no distinction between independent and dependent coordinates is made in this linearization approach. Procedure by Van Khang et al. [223] The approach of Van Khang et al. [223] is also based on the linearization of the index-1 DAE system (2.127). In this approach, the coordinate partition in n−m−lindependent and m+ldependent coordinates, described in Sect. 2.4.2, is used. Using the transformation matrix T(x) presented in Eq. (2.130), and premultiplying the dynamic equations (2.113) by TT(x), the index-1 DAE system (2.127) is transformed into the following ODE system: TT(x)M(x)¨x=TT(x)Q(x,˙x),(3.40) D(x)¨x=Qd(x,˙x).(3.41) The Lagrange multipliers Λin Eq. (2.127) are eliminated since TT(x)DT(x) = 0. Presenting the following definitions: f1(x,¨x) = TT(x)M(x)¨x,k1(x,˙x) = TT(x)Q(x,˙x), f2(x,¨x) = D(x)¨x,k2(x,˙x) = Qd(x,˙x),(3.42) the nonlinear ODE system (3.40)-(3.41) can be written as: f1(x,¨x)−k1(x,˙x) = 0,(3.43) f2(x,¨x)−k2(x,˙x) = 0.(3.44) The linearization of Eqs. (3.43)-(3.44) leads to: ∂f1 ∂¨x0 ¨ ˜x−∂k1 ∂˙x0 ˙ ˜x+∂f1 ∂x0−∂k1 ∂x0˜x=0,(3.45) ∂f2 ∂¨x0 ¨ ˜x−∂k2 ∂˙x0 ˙ ˜x+∂f2 ∂x0−∂k2 ∂x0˜x=0.(3.46) The partial derivatives ∂f1 ∂x0 and ∂f2 ∂x0 include the second derivative with respect to time of the reference solution, ¨x0, and ∂k1 ∂x0 ,∂k2 ∂x0 include the first derivative with respect to time of the reference solution, ˙x0. The system of linear equations (3.45)-(3.46) can be expressed as: ¯ M¨ ˜x+¯ C˙ ˜x+¯ K˜x=0,(3.47) 87 Chapter 3. An overview of linearization approaches in Multibody System Dynamics where ¯ M,¯ Cand ¯ Kare n×nmatrices, representing the composite mass, damping and stiffness matrices, respectively, given by: ¯ M=    ∂f1 ∂¨x0 ∂f2 ∂¨x0    ,¯ C=    −∂k1 ∂˙x0 −∂k2 ∂˙x0    ,¯ K=    ∂f1 ∂x0−∂k1 ∂x0 ∂f2 ∂x0−∂k2 ∂x0    .(3.48) The symbolic computation of matrices (3.48) is further developed in the work of Van Khang et al. [223]. Defining ˜ X=˜x˙ ˜xT, Eq. (3.47) can be written as the following first-order ODE system: ˙ ˜ X=J˜ X,with J=0nIn −¯ M−1¯ K−¯ M−1¯ C.(3.49) As in the previous approach, the resulting Jacobian matrix (3.49) is 2n×2n. This leads to 2neigenvalues, among which 2 (n−m−l) corresponds to the real spectrum of the problem, and 2 (m+l) are spurious eigenvalues due to the m+l dependent coordinates. Note that, despite this approach makes use of the coordinate partition in terms of independent and dependent coordinates described in Sect. 2.4.2, the resulting linearized equations of motion (3.47) are expressed in terms of the variations of the redundant set of coordinates ˜xand their time derivatives ˙ ˜xand ¨ ˜x. Therefore, the dependent coordinates are not eliminated in this approach. For this reason, and given that the procedure is based on the linearization of the index-1 DAE system, this approach has been included in Group 2 (Sect. 3.2) instead of Group 3 (Sect. 3.3). 3.3 Group 3: procedures based on coordinate partitioning in terms of independent and dependent coordinates The approaches included in this section are based on the linearization of the ODE system of equations (2.134)-(2.135), expressed in terms of the independent and dependent coordinates. By virtue of the coordinate partition of Eq. (2.129), the reference solution x0(t) and its time derivatives ˙x0(t), ¨x0(t), presented in Sect. 3.1, can be expressed as: x0 ai(t),˙x0 ai(t),¨x0 ai(t),x0 dd(t),˙x0 dd(t).(3.50) The following variations with respect to the reference solution (3.50) are introduced: ˜xai =xai −x0 ai,˙ ˜xai =˙xai −˙x0 ai, ¨ ˜xai =¨xai −¨x0 ai,˜xdd =xdd −x0 dd,(3.51) 88 3.3. Group 3: procedures based on coordinate partitioning in terms of independent and dependent coordinates ˙ ˜xdd =˙xdd −˙x0 dd. Procedure by Cuadrado et al. [225, 226] Cuadrado et al. [225, 226] presented a linearization approach for multibody systems with holonomic constraints. In [225, 226], the dynamic equations (2.134) are written as the first-order ODE system: ˙xai =vai,(3.52) ˙vai = (m(xai,xdd))−1f(xai,xdd,vai) = g(xai,xdd,vai),(3.53) where the matrix mand the vector fwere defined in Eq. (2.136). Defining X=xai vai T, Eqs. (3.52)-(3.53) can be expressed in compact form as: ˙ X=G1(xai,xdd,vai).(3.54) Introducing the vector of variations ˜ X=˜xai ˜vai T, the linearization of the first-order system (3.54) along the reference solution (3.50) leads to a first-order linear ODE system of the form: ˙ ˜ X=J˜ X,with J=  0(n−m)I(n−m) ∂g ∂xai 0 ∂g ∂˙xai 0  .(3.55) The state-space Jacobian matrix Jof Eq. (3.55) is 2 (n−m)×2 (n−m), leading to 2 (n−m) eigenvalues. These eigenvalues correspond to the real spectrum of the multibody system, which presents ng=n−mdegrees of freedom. This procedure can be found in several works. Cuadrado et al. [225, 226] and Sanjurjo et al. [227] followed this approach to construct the state-space Jacobian matrix, in order to build continuous extended Kalman filters. Procedure by Gonz´alez et al. [218] (MCS Formulation) A similar procedure is presented by Gonz´alez et al. [218], under the name of Minimal Coordinate Set (MCS) Formulation: Velocity Projection. As in the case of Cuadrado et al. [225, 226], the approach is presented for multibody systems with holonomic constraints. The dynamic equations (2.134) are written as: H1=m(xai,xdd)¨xai −f(xai,xdd,˙xai) = 0,(3.56) where, according to Eq. (2.136), m=TTMT and f=TTQ−M˙ T˙xai. The linearization of Eq. (3.56) along the reference solution yields: ∂H1 ∂xai 0 ˜xai+∂H1 ∂˙xai 0 ˙ ˜xai +∂H1 ∂¨xai 0 ¨ ˜xai =0,(3.57) where the partial derivatives in Eq. (3.57) are given by: ∂H1 ∂xai =∂TT ∂xai MT +TT∂M ∂xai T+TTM∂T ∂xai ¨xai −∂TT ∂xai Q−M˙ T˙xai 89 Chapter 3. An overview of linearization approaches in Multibody System Dynamics 96 Chapter 4 Novel linearization approaches for constrained multibody systems The objective of this chapter is to present general linearization approaches, which could be implemented symbolically or numerically, and applied to any kind of holonomic or nonholonomic multibody system. First, two symbolic linearization approaches are proposed. The first procedure transforms the original nonlinear DAE system into a nonlinear ODE system of equations by using a coordinate partition based on independent and dependent coordinates. Next, the linearization of this ODE system is performed, including the blocks associated with the derivatives with respect to the dependent coordinates in the resulting Jacobian matrix. The second symbolic approach carries out a further reduction of the Jacobian matrix, and the resulting linearized equations of motion are obtained only in terms of the variations of the independent coordinates and their time derivatives. The dependency of the dependent coordinates with the independent ones is considered by using the linearized holonomic constraints. These two procedures are labeled as symbolic because they involve certain steps that can only be carried out symbolically. The first of them falls within Group 3 of Sect. 3.3, and the second belongs to Group 4 of Sect. 3.4. On the other hand, three approaches, which can be implemented numerically and symbolically, are presented. The first of them is based on the direct linearization of the index-1 DAE system, being an excellent alternative to the approaches presented in Group 2 of Sect. 3.2. The second numerical approach is a counterpart procedure to the first symbolic linearization approach presented in this document, and significantly reduces the computational cost for the computation of the linearized equations of motion. While the approaches found in [171, 218, 225, 226, 227] transform the original nonlinear DAE system into a nonlinear ODE system of equations by using the coordinate partition, and sub- 97 Chapter 4. Novel linearization approaches for constrained multibody systems sequently perform the linearization, this approach first linearizes the DAE system and then carries out the reduction to a linear ODE system. The resulting Jacobian matrix also includes the blocks associated with the derivatives with respect to the dependent coordinates. Finally, the third numerical procedure follows the steps of the second numerical approach, and achieves the maximum possible reduction of the linearized equations of motion by considering the relation between the dependent and independent coordinates with the linearized holonomic constraints. The correctness of all the linearization approaches is validated with the linear stability results of the bicycle benchmark of Meijaard et al. [168]. Novel linearization approaches First numerical approach (Sect. 4.2.1) Second numerical approach (Sect. 4.2.2) Third numerical approach (Sect. 4.2.3) Second symbolic approach (Sect. 4.1.2) First symbolic approach (Sect. 4.1.1) Figure 4.1: Novel linearization approaches presented in Chapter 4. The novel linearization procedures are presented by using the same notation as in Chapter 3. The linearization is performed along a known reference solution of the DAE system (2.122)-(2.124), and therefore Eqs. (3.1)-(3.3) are fulfilled. The time evolution of the generalized coordinates in this reference solution, at position, velocity and acceleration levels, are defined by x0(t), ˙x0(t), ¨x0(t). Furthermore, the Lagrange multipliers in the reference solution are denoted as Λ0(t). Thus, the 98 4.1. Procedures based on symbolic computation reference solution can be expressed as: x0(t),˙x0(t),¨x0(t),Λ0(t).(4.1) 4.1 Procedures based on symbolic computation The two linearization procedures based on symbolic computation are shown below. 4.1.1 First symbolic procedure The first symbolic approach is based on the coordinate partition (2.129), i.e., x=xai xdd T, and the linearization of the ODE system of equations (2.134)- (2.135), which is shown again for convenience: m(xai,xdd)¨xai =f(xai,xdd,˙xai),(4.2) ˙xdd =H(xai,xdd)˙xai,(4.3) with m(xai,xdd) = TTMT , f(xai,xdd,˙xai) = TTQ−M˙ T˙xai, H(xai,xdd) = −(Ddd)−1Dai. (4.4) This linearization approach belongs to Group 3 of Sect. 3.3. The reference solution (4.1) can be expressed, by virtue of the coordinate partition (2.129), as: x0 ai(t),˙x0 ai(t),¨x0 ai(t),x0 dd(t),˙x0 dd(t).(4.5) The following variations with respect to the reference solution (4.5) are introduced: ˜xai =xai −x0 ai,˙ ˜xai =˙xai −˙x0 ai, ¨ ˜xai =¨xai −¨x0 ai,˜xdd =xdd −x0 dd,(4.6) ˙ ˜xdd =˙xdd −˙x0 dd. The Taylor expansion of the equations of motion (4.2)-(4.3) yields:  m0+ n−m−l X j=1 ∂m ∂xaij0 ˜xaij+ m+l X j=1 ∂m ∂xddj0 ˜xddj+... ¨ ˜xai +¨x0 ai= f0+∂f ∂xai 0 ˜xai +∂f ∂xdd 0 ˜xdd +∂f ∂˙xai 0 ˙ ˜xai +..., (4.7) 99 Chapter 4. Novel linearization approaches for constrained multibody systems ˙ ˜xdd +˙x0 dd = H0+ n−m−l X j=1 ∂H ∂xaij0 ˜xaij+ m+l X j=1 ∂H ∂xddj0 ˜xddj+... ˙ ˜xai +˙x0 ai, (4.8) where the dependence of m,fand Hwith (xai,xdd), detailed in Eqs. (4.4), is omitted for simplicity and the following definitions are used: m0=mx0 ai,x0 dd,(4.9) f0=fx0 ai,x0 dd,˙x0 ai,(4.10) H0=Hx0 ai,x0 dd.(4.11) Since the reference motion fulfils the equations of motion (4.2)-(4.3), the following relations are verified: m0¨x0 ai =f0,(4.12) ˙x0 dd =H0˙x0 ai.(4.13) Using Eqs. (4.12) and (4.13), and retaining up to first-order tems in the system of equations (4.7)-(4.8), the following linear ODE system is obtained: m0¨ ˜xai +∂m¨x0 ai ∂xai 0 ˜xai +∂m¨x0 ai ∂xdd 0 ˜xdd = ∂f ∂xai 0 ˜xai +∂f ∂xdd 0 ˜xdd +∂f ∂˙xai 0 ˙ ˜xai,(4.14) ˙ ˜xdd =H0˙ ˜xai +∂H˙x0 ai ∂xai 0 ˜xai +∂H˙x0 ai ∂xdd 0 ˜xdd,(4.15) where the reference velocity and acceleration vectors ˙x0 ai and ¨x0 ai have been grouped with Hand m, respectively, for a simpler notation. Rearranging terms in Eqs. (4.14)-(4.15), one obtains: m0¨ ˜xai =∂f−m¨x0 ai ∂xai 0 ˜xai +∂f ∂˙xai 0 ˙ ˜xai +∂f−m¨x0 ai ∂xdd 0 ˜xdd, ˙ ˜xdd =∂H˙x0 ai ∂xai 0 ˜xai +H0˙ ˜xai +∂H˙x0 ai ∂xdd 0 ˜xdd. (4.16) Defining ˜ X=˜xai ˙ ˜xai ˜xdd T, Eqs. (4.16) can be rewritten as a first- 100 4.1. Procedures based on symbolic computation order linear ODE system of the form ˙ ˜ X=J˜ X, where Jis the Jacobian matrix: J=         0(n−m−l)I(n−m−l)0(n−m−l)×(m+l) m−1 0 ∂f−m¨x0 ai ∂xai 0 m−1 0 ∂f ∂˙xai 0 m−1 0 ∂f−m¨x0 ai ∂xdd 0 ∂H˙x0 ai ∂xai 0 H0 ∂H˙x0 ai ∂xdd 0         .(4.17) Observe that the size of Jis (2n−m−l)×(2n−m−l). Note that the blocks associated with the dependent coordinates have been considered in the Jacobian matrix (4.17). The linearized equations of motion (4.16) are not only expressed in terms of the independent variations ˜xai and their time derivatives ˙ ˜xai and ¨ ˜xai, but also include the variations of the dependent coordinates ˜xdd and the time derivative ˙ ˜xdd. The approach is labeled as symbolic since the Jacobian matrix (4.17) requires the symbolic computation of derivatives of inverse matrices, as matrix ˙ Tin Eq. (4.4). 4.1.2 Second symbolic procedure The second symbolic approach is based on the coordinate partition (2.137), i.e., x=xai xad xdT, and the linearization of the equations of motion presented in Sect. 2.4.3. These equations of motion, which are shown again for convenience, constitute a nonlinear DAE system without Lagrange multipliers: m(xai,xad,xd)¨xai =f(xai,xad,xd,˙xai),(4.18) C(xai,xad,xd) = 0,(4.19) ˙xad =G(xai,xad,xd)˙xai,(4.20) where mand fcan be found in Eqs. (4.4), and Gwas presented after Eq. (2.144). This linearization approach belongs to Group 4 of Sect. 3.4. The known reference motion (4.1) can be expressed, by virtue of the coordinate partition (2.137), as: x0 ai(t),˙x0 ai(t),¨x0 ai(t),x0 ad(t),˙x0 ad(t),x0 d(t).(4.21) The following variations with respect to the reference solution (4.21) are defined: ˜xai =xai −x0 ai,˙ ˜xai =˙xai −˙x0 ai, ¨ ˜xai =¨xai −¨x0 ai,˜xad =xad −x0 ad,(4.22) ˙ ˜xad =˙xad −˙x0 ad,˜xd=xd−x0 d. The Taylor expansion of the equations of motion (4.18)-(4.20) yields:  m0+ n−m−l X j=1 ∂m ∂xaij0 ˜xaij+ l X j=1 ∂m ∂xadj0 ˜xadj 101 Chapter 4. Novel linearization approaches for constrained multibody systems + m X j=1 ∂m ∂xdj0 ˜xdj+... ¨ ˜xai +¨x0 ai=f0+∂f ∂xai 0 ˜xai +∂f ∂xad 0 ˜xad +∂f ∂xd0 ˜xd+∂f ∂˙xai 0 ˙ ˜xai +..., (4.23) C0+∂C ∂xai 0 ˜xai +∂C ∂xad 0 ˜xad +∂C ∂xd0 ˜xd=0,(4.24) ˙ ˜xad +˙x0 ad = G0+ n−m−l X j=1 ∂G ∂xaij0 ˜xaij+ l X j=1 ∂G ∂xadj0 ˜xadj + m X j=1 ∂G ∂xdj0 ˜xdj+... ˙ ˜xai +˙x0 ai,(4.25) where the dependence of m,fand Gwith (xai,xad,xd), detailed in Eqs. (4.18)- (4.20), is omitted for simplicity and the following definitions have been used: m0=mx0 ai,x0 ad,x0 d,(4.26) f0=fx0 ai,x0 ad,x0 d,˙x0 ai,(4.27) C0=Cx0 ai,x0 ad,x0 d,(4.28) G0=Gx0 ai,x0 ad,x0 d.(4.29) Since the reference motion fulfils the equations of motion (4.18)-(4.20), the following relations are verified: m0¨x0 ai =f0,(4.30) C0=0,(4.31) ˙x0 ad =G0˙x0 ai.(4.32) Using Eqs. (4.30)-(4.32), and neglecting second-order and higher tems in the system of equations (4.23)-(4.25), the following linear DAE system is obtained: m0¨ ˜xai +∂m¨x0 ai ∂xai 0 ˜xai +∂m¨x0 ai ∂xad 0 ˜xad +∂m¨x0 ai ∂xd0 ˜xd= ∂f ∂xai 0 ˜xai +∂f ∂xad 0 ˜xad +∂f ∂xd0 ˜xd+∂f ∂˙xai 0 ˙ ˜xai,(4.33) ∂C ∂xai 0 ˜xai +∂C ∂xad 0 ˜xad +∂C ∂xd0 ˜xd=0,(4.34) ˙ ˜xad =G0˙ ˜xai +∂G˙x0 ai ∂xai 0 ˜xai +∂G˙x0 ai ∂xad 0 ˜xad +∂G˙x0 ai ∂xd0 ˜xd,(4.35) 102 4.1. Procedures based on symbolic computation where the reference velocity and acceleration vectors ˙x0 ai and ¨x0 ai have been grouped with Gand m, respectively. Therefore, Eqs. (4.33)-(4.35) are further expressed as m0¨ ˜xai =∂f−m¨x0 ai ∂xai 0 ˜xai +∂f−m¨x0 ai ∂xad 0 ˜xad +∂f−m¨x0 ai ∂xd0 ˜xd+∂f ∂˙xai 0 ˙ ˜xai,(4.36) ∂C ∂xai 0 ˜xai +∂C ∂xad 0 ˜xad +∂C ∂xd0 ˜xd=0,(4.37) ˙ ˜xad =G0˙ ˜xai +∂G˙x0 ai ∂xai 0 ˜xai +∂G˙x0 ai ∂xad 0 ˜xad +∂G˙x0 ai ∂xd0 ˜xd.(4.38) The linearized holonomic constraints in Eq. (4.37) can be used to eliminate the dependent coordinates ˜xdin Eqs. (4.36) and (4.38), obtaining the following relation: ˜xd=−∂C ∂xd0−1∂C ∂xai 0 ˜xai −∂C ∂xd0−1∂C ∂xad 0 ˜xad.(4.39) Substituting Eq. (4.39) in Eqs. (4.36) and (4.38), yields m0¨ ˜xai = ∂f−m¨x0 ai ∂xai −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 ˜xai + ∂f−m¨x0 ai ∂xad −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 ˜xad +∂f ∂˙xai 0 ˙ ˜xai, (4.40) ˙ ˜xad =G0˙ ˜xai + ∂G˙x0 ai ∂xai −∂G˙x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 ˜xai + ∂G˙x0 ai ∂xad −∂G˙x0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 ˜xad.(4.41) Note that, in contrast to the dependent variations ˜xd, the set ˜xad cannot be removed, since the nonholonomic constraints are first-order differential equations that cannot be eliminated. Finally, defining ˜ X=˜xai ˙ ˜xai ˜xad T, the system of equations (4.40)- (4.41) can be rewritten as a first-order system of the form ˙ ˜ X=J˜ X, with the Jacobian matrix J: J=   0(n−m−l)I(n−m−l)0(n−m−l)×l J21 J22 J23 J31 J32 J33   ,(4.42) 103 Chapter 4. Novel linearization approaches for constrained multibody systems where J21 =m−1 0 ∂f−m¨x0 ai ∂xai −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 , J22 =m−1 0 ∂f ∂˙xai 0 , J23 =m−1 0 ∂f−m¨x0 ai ∂xad −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 , J31 = ∂G˙x0 ai ∂xai −∂G˙x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 , J32 =G0, J33 = ∂G˙x0 ai ∂xad −∂G˙x0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 . (4.43) Observe that, in this case, the Jacobian matrix Jis (2n−2m−l)×(2n−2m−l). This size is the minimum that can be achieved and cannot be further reduced. In the particular case of a holonomic multibody system, the linear equations of motion in ODE form (4.40)-(4.41) can be expressed in terms of only the independent variations ˜xai and their time derivatives. The linearization of Eqs. (4.18)- (4.19) with respect to the reference solution (4.21) yields the following linear DAE system: m0¨ ˜xai =∂f−m¨x0 ai ∂xai 0 ˜xai +∂f−m¨x0 ai ∂xd0 ˜xd+∂f ∂˙xai 0 ˙ ˜xai,(4.44) ∂C ∂xai 0 ˜xai +∂C ∂xd0 ˜xd=0.(4.45) The linearized holonomic constraints (4.45) can be used to express the variations of the dependent coordinates ˜xdas a function of the independent ones ˜xai: ˜xd=−∂C ∂xd0−1∂C ∂xai 0 ˜xai.(4.46) The substitution of Eq. (4.46) in the linearized dynamic equations (4.44) leads to: m0¨ ˜xai = ∂f−m¨x0 ai ∂xai −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 ˜xai +∂f ∂˙xai 0 ˙ ˜xai.(4.47) 104 4.2. Procedures based on numerical computation Defining ˜ X=˜xai ˙ ˜xai T, the linear second-order ODE of Eq. (4.47) can be written as a first-order ODE of the form ˙ ˜ X=J˜ X, with the Jacobian matrix J: J= 0(n−m)I(n−m) J21 J22 !,(4.48) where the blocks J21 and J22 in Eq. (4.48) are given by: J21 =m−1 0 ∂f−m¨x0 ai ∂xai −∂f−m¨x0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 , J22 =m−1 0 ∂f ∂˙xai 0 . (4.49) In this case, the Jacobian matrix (4.48) is 2 (n−m)×2 (n−m). As in the first symbolic approach, this procedure is also labeled as symbolic, since the linearization of the equations of motion requires the symbolic computation of the matrix ˙ Tin Eq. (4.4), which entails the computation of derivatives of inverse matrices. Nevertheless, it is unlikely that these symbolic linearization procedures can be applied in general, given that the computation of ˙ Tdemands a high computational cost. While for simple multibody systems (low number of generalized coordinates and short kinematic chains) they might be convenient, the symbolic procedures might be unaffordable for moderately large systems. 4.2 Procedures based on numerical computation The following three novel procedures are presented below. It is important to clarify that, although these approaches are called numerical, they also allow for computing the linearized equations of motion analytically. 4.2.1 First numerical procedure The first numerical linearization approach is based on the direct linearization of the index-1 DAE system (2.125)-(2.126) (or, alternately, Eq. (2.127)), and therefore it is classified within Group 1 of Sect. 3.1. The main results of the intermediate steps are shown in a box. First, the variations with respect to the reference solution (8.7) are presented again for convenience: ˜x=x−x0,˙ ˜x=˙x−˙x0, ¨ ˜x=¨x−¨x0,˜ Λ=Λ−Λ0.(4.50) The Taylor expansion of the dynamic equilibrium equations (2.125), with respect to the reference solution, leads to: 105 Chapter 4. Novel linearization approaches for constrained multibody systems The Jacobian matrices of Eqs. (4.17) and (4.95) are (2n−m−l)×(2n−m−l). Among the 2n−m−leigenvalues obtained with these approaches, 2 (n−m−l) correspond to the real spectrum of the problem, and m+lare associated with the mholonomic constraints differentiated with respect to time and the lnonholonomic constraints. Note that the resulting linearized equations of motion (4.89)-(4.90) present the same structure as Eq. (4.16). In this approach, the dependent coordinates are handled in the same way as in the symbolic counterpart. As in Eq. (4.16), the linearized equations of motion (4.89)-(4.90) are not only expressed in terms of the independent variations ˜xai and their time derivatives ˙ ˜xai and ¨ ˜xai, but also include the variations of the dependent coordinates ˜xdd and the time derivatives ˙ ˜xdd. 4.2.3 Third numerical procedure A counterpart procedure to the second symbolic linearization approach, which allows the symbolic and numerical computation of the Jacobian matrix (4.42), is presented. The coordinate partition of Eq. (2.137) is used. In the linearization approach, this partition is considered for the vector of variations ˜x: ˜x=˜xai ˜xad ˜xdT.(4.96) The admissible dependent coordinates ˜xad and the dependent coordinates ˜xdcan be grouped in the set ˜xdd =˜xad ˜xdT. Since obtaining the linearized equations with this approach involves several steps, they are listed below, with the main result of each step shown in a box. Step 1. Eliminate the variations of the Lagrange multipliers and reduce the number of linearized dynamic equations. This step was also done in the second numerical procedure. The variations of the Lagrange multipliers ˜ Λcan be eliminated by premultiplying Eq. (4.52) by TTx0, given that TTx0DTx0=0. Defining T0=Tx0, the reduced linearized dynamic equations are obtained: TT 0Mx0¨ ˜x+TT 0 ∂M(x)¨x0 ∂x0 ˜x+TT 0 ∂DT(x)Λ0 ∂x0 ˜x= TT 0∂Q ∂x0 ˜x+∂Q ∂˙x0 ˙ ˜x. (4.97) Step 2. Obtain a transformation at velocity level: express the velocities ˙ ˜xin terms of the independent velocities ˙ ˜xai and positions ˜x. The objective is to express Eq. (4.97) in terms of the variations of the independent coordinates ˜xai and their time derivatives ˙ ˜xai,¨ ˜xai. By using the linearized 112 4.2. Procedures based on numerical computation velocity constraints (4.78), the variations of the dependent velocities are written as follows: ˙ ˜xdd =Tdd x0˙ ˜xai +¯ Vx0,˙x0˜x,(4.98) where Tdd x0=−Ddd x0−1Dai x0, ¯ Vx0,˙x0=−Ddd x0−1∂D(x)˙x0 ∂x0 .(4.99) From Eq. (4.98), the following transformation at velocity level is obtained: ˙ ˜x=Tx0˙ ˜xai +¯ ¯ Vx0,˙x0˜x,(4.100) where ¯ ¯ Vx0,˙x0=0(n−m−l)×n ¯ Vx0,˙x0. Step 3. Obtain a transformation at acceleration level: express the accelerations ¨ ˜xin terms of the independent accelerations ¨ ˜xai, velocities ˙ ˜xand positions ˜x. Next, the linearized acceleration constraints (4.79) are used to express the dependent accelerations ¨ ˜xdd as a function of ¨ ˜xai, and the following relation at acceleration level is derived: ¨ ˜x=Tx0¨ ˜xai +Ux0,˙x0˙ ˜x+Vx0,˙x0,¨x0˜x,(4.101) with Ux0,˙x0=  0(n−m−l)×n Ddd x0−1∂Qd ∂˙x0  ,(4.102) Vx0,˙x0,¨x0=   0(n−m−l)×n −Ddd x0−1 ∂D(x)¨x0 ∂x0−∂Qd ∂x0!  .(4.103) Steps 2 and 3 were already followed in the second numerical approach. Step 4. Obtain a transformation at position level: eliminate the dependent coordinates ˜xdby using the linearized holonomic constraints. To achieve the maximum reduction of the linearized equations of motion, the dependent coordinates ˜xdare eliminated by linearizing the holonomic constraints (4.19): ∂C ∂xai 0 ˜xai +∂C ∂xad 0 ˜xad +∂C ∂xd0 ˜xd=0.(4.104) Therefore, Eq. (4.104) allows expressing the dependent coordinates ˜xdas a function of ˜xai and ˜xad: ˜xd=−∂C ∂xd0−1∂C ∂xai 0 ˜xai +∂C ∂xad 0 ˜xad.(4.105) 113 Chapter 4. Novel linearization approaches for constrained multibody systems Note that Eq. (4.105) was already used in the second symbolic approach. The use of Eq. (4.105) leads to the following transformation at position level: ˜x=Vh,ai 0˜xai +Vh,ad 0˜xad,(4.106) where Vh,ai 0=   I(n−m−l) 0l×(n−m−l) −∂C ∂xd0−1∂C ∂xai 0   ,(4.107) Vh,ad 0=   0(n−m−l)×l Il −∂C ∂xd0−1∂C ∂xad 0   .(4.108) Step 5. Express the velocities ˙ ˜xad in terms of ˜xai,˙ ˜xai and ˜xad: use of the time derivative of the holonomic constraints. Furthermore, an expression of the velocities ˙ ˜xad as a function of ˜xai,˙ ˜xai and ˜xad can be derived from the time derivative of the holonomic constraints and the nonholonomic constraints. Computing the time derivative of the holonomic constraints: Cx(x)˙x=0,(4.109) the linearization of Eq. (4.109) with respect to the reference solution yields: Cxx0˙ ˜x+∂Cx(x)˙x0 ∂x0 ˜x=0.(4.110) Resorting to Eq. (4.110), the following transformation is obtained: ˙ ˜x=Th 0˙ ˜xa+Vh 0˜x,(4.111) where Th 0=I(n−m) −Cxdx0−1Cxax0,(4.112) Vh 0=   0(n−m)×n −Cxdx0−1∂Cx(x)˙x0 ∂x0   .(4.113) Step 6. Express the velocities ˙ ˜xad in terms of ˜xai,˙ ˜xai and ˜xad: use of the nonholonomic constraints. Finally, the nonholonomic constraints (2.124) are linearized: Bx0˙ ˜x+∂B(x)˙x0 ∂x0 ˜x=0.(4.114) 114 4.2. Procedures based on numerical computation Substituting Eq. (4.111) in Eq. (4.114): Bx0Th 0˙ ˜xa+ Bx0Vh 0+∂B(x)˙x0 ∂x0!˜x=0,(4.115) and defining ¯ Bx0=Bx0Th 0, the following expression for the velocities ˙ ˜xad is obtained: ˙ ˜xad =Unh,ai 0˙ ˜xai +Vnh,ai 0˜xai +Vnh,ad 0˜xad,(4.116) where Unh,ai 0=−¯ Bad x0−1¯ Bai x0, Vnh,ai 0=−¯ Bad x0−1 Bx0Vh 0+∂B(x)˙x0 ∂x0!Vh,ai 0, Vnh,ad 0=−¯ Bad x0−1 Bx0Vh 0+∂B(x)˙x0 ∂x0!Vh,ad 0. (4.117) Step 7. Obtain the linearized equations of motion and the Jacobian matrix. The transformations at velocity level (result of Step 2), acceleration level (result of Step 3) and position level (result of Step 4), given by Eqs. (4.100), (4.101) and (4.106), respectively, are used in the linearized dynamic equations (4.97) (result of Step 1). Furthermore, including Eq. (4.116) (result of Step 5 and Step 6), yields the following linear ODE system of equations: m0¨ ˜xai =R0¯ ¯ V0+S0Vh,ai 0˜xai +R0T0˙ ˜xai +R0¯ ¯ V0+S0Vh,ad 0˜xad,(4.118) ˙ ˜xad =Vnh,ai 0˜xai +Unh,ai 0˙ ˜xai +Vnh,ad 0˜xad,(4.119) where ¯ ¯ V0=¯ ¯ Vx0,˙x0. The matrices m0,R0and S0in Eqs. (4.118)-(4.119) are given by: m0=TT 0Mx0T0,(4.120) R0=TT 0∂Q ∂˙x0−Mx0Ux0,˙x0,(4.121) S0=TT 0∂Q ∂x0−Mx0Vx0,˙x0,¨x0 −∂M(x)¨x0 ∂x0−∂DT(x)Λ0 ∂x0!.(4.122) Note that the linear ODE system (4.118)-(4.119) presents the same structure as Eqs. (4.40)-(4.41). Equation (4.118) corresponds to the n−m−llinearized 115 Chapter 4. Novel linearization approaches for constrained multibody systems dynamic equations, and Eq. (4.119) to the llinearized nonholonomic constraints. It must be pointed out that, in contrast to the holonomic constraints, which are nonlinear algebraic equations that can be eliminated after the linearization, the nonholonomic constraints cannot be eliminated due to their non-integrable nature. Finally, defining ˜ X=˜xai ˙ ˜xai ˜xad T, the system of equations (4.118)- (4.119) can be rewritten as a first-order system of the form ˙ ˜ X=J˜ X,(4.123) with the Jacobian matrix J: J=   0(n−m−l)I(n−m−l)0(n−m−l)×l J21 J22 J23 J31 J32 J33   ,(4.124) where J21 =m−1 0R0¯ ¯ V0+S0Vh,ai 0, J22 =m−1 0R0T0, J23 =m−1 0R0¯ ¯ V0+S0Vh,ad 0, J31 =Vnh,ai 0, J32 =Unh,ai 0, J33 =Vnh,ad 0. (4.125) As in Eq. (4.42), the Jacobian matrix (4.124) is (2 (n−m)−l)×(2 (n−m)−l). This size is the minimum that can be achieved, not being possible a further reduction, since the nonholonomic constraints are first-order differential equations that cannot be eliminated. Among the 2 (n−m)−leigenvalues obtained with this approach, 2 (n−m−l) correspond to the real spectrum of the problem, and lare spurious eigenvalues associated with the llinearized nonholonomic constraints. In the particular case of a holonomic multibody system, without nonholonomic constraints, the linearized dynamic equations (4.118) are expressed in terms of the independent variations ˜xai and their time derivatives ˙ ˜xai,¨ ˜xai: m0¨ ˜xai =R0¯ ¯ V0+S0Vh,ai 0˜xai +R0T0˙ ˜xai.(4.126) Defining ˜ X=˜xai ˙ ˜xai T, a first-order ODE system of the form ˙ ˜ X=J˜ Xis obtained, with the Jacobian matrix J: J= 0(n−m)I(n−m) J21 J22 !,(4.127) 116 4.3. Validation of the linearization approaches with the bicycle benchmark where the blocks J21 and J22 in Eq. (4.127) are those of Eq. (4.125). The Jacobian matrix (4.127) is 2 (n−m)×2 (n−m). Therefore, 2 (n−m) eigenvalues are obtained, corresponding to the real spectrum of the problem. Note that, despite the linearized equations of motion (4.126) are only expressed in terms of the independent variations ˜xai and their time derivatives, the dependent coordinates are not ignored in the linearization. The dependency with the independent ones has been taken into accout by means of the relations at velocity, acceleration and position levels (4.100), (4.101) and (4.106), respectively. 4.3 Validation of the linearization approaches with the bicycle benchmark The accuracy of the linearization approaches developed in Sects. 4.1 and 4.2 is validated with the linear stability results of the well-acknowledged bicycle benchmark of Meijaard et al. [168]. This benchmark has been widely used in several experimental [169] and theoretical works [170, 171, 172] devoted to multibody dynamics. The linear stability of the steady forward and circular motions of the bicycle benchmark was meticulously analyzed by Basu-Mandal et al. [170] and Xiong et al. [172, 173]. Furthermore, more advanced models of frames, tyres and riders were introduced in different extensions of the bicycle benchmark [176, 177, 178, 179, 180]. Therefore, a bicycle multibody model, based on this bicycle benchmark, is built by using the Newton-Euler formalism, including holonomic and nonholonomic constraints. The resulting nonlinear DAE equations of motion are linearized by using all the novel symbolic and numerical approaches along the steady forward motion of the bicycle. 4.3.1 Multibody model of the bicycle benchmark The multibody model consists of five rigid bodies. The rear body and frame assembly, which includes the rider, is denoted as B; the front handlebar and fork assembly is represented by H; and the rear and front wheels, which are modelled as hoops rolling without slipping, are denoted as Rand F, respectively. Lastly, the global frame corresponds to body 1. The origins of the body reference frames are located at the respective centres of mass, denoted by Gi, with i={B, H, R, F}. The following n×1 vector of coordinates, with n= 11, is used to describe the system: x=xbybzbψbϕbθbδ θRθFξRξFT.(4.128) In Eq. (4.128), xb,yband zbare Cartesian coordinates that locate the centre of mass of the rear body and frame assembly B, while the triplet {ψb, ϕb, θb}are the orientation angles of body Bin space. In particular, ψbis the yaw angle, ϕbis the lean angle and θbrepresents the pitch. The steering angle δcorresponds to the rotation of the handlebar with respect to the rear body, and θR,θFare used 117 Chapter 4. Novel linearization approaches for constrained multibody systems to describe the rotations of the rear and front wheels with respect to bodies B and H, respectively. Lastly, ξRand ξFare two angular coordinates that allow for locating the radii which contain the contact points Pand Qin the body reference frames of the rear and front wheels, respectively. Figure 4.2 shows the numbering of the bodies, the generalized coordinates of the system and its body reference frames. The orientation matrices of the body frames, expressed as a function of the elemental rotation matrices, are given by: RB(x) = RψbRϕbRθb, RR(x) = RBRθR, RH(x) = RBRνRδ, RF(x) = RHRθF, (4.129) where Rνconsiders the steer axis tilt angle, denoted by ν. Figure 4.2: Three-dimensional multibody model of the bicycle benchmark: numbering of the bodies, generalized coordinates and body reference frames. To build the multibody model, the geometric and dynamic parameters of Meijaard et al. [168] are used, being listed in Table 4.1 with their corresponding numerical values. Considering the orientation of the global and body reference 118 4.3. Validation of the linearization approaches with the bicycle benchmark frames shown in Fig. 4.2, the inertia tensors of the bodies Band F, denoted as ¯ IBand ¯ IF, respectively, can be expressed in their corresponding body frames as: ¯ IB=RT πIBRπ, ¯ IH=RT νRT πIHRπRν.(4.130) In Eq. (4.130), IBand IHare the inertia tensors about a local reference frame fixed to GBand GH, parallel to the reference frame of Meijaard et al. [168]. Given that this reference frame uses an axes convention with the positive Z-axis pointing downwards and the positive Y-axis aiming at the right of the rider, the transformations in Eq. (4.130) are required, where Rπrepresents a counterclockwise rotation of 180◦about the X-axis. Therefore, the components of ¯ IBand ¯ IF can be expressed as a function of the moments of inertia provided by Meijaard et al. [168] as follows: ¯ IBxx =IBxx , ¯ IByy =IByy , ¯ IBzz =IBzz , ¯ IBxz =−IBxz , ¯ IHxx =IHxx cos2(ν) + IHzz sin2(ν)−IHxz sin (2ν), ¯ IHyy =IHyy , ¯ IHzz =IHxx sin2(ν) + IHzz cos2(ν) + IHxz sin (2ν), ¯ IHxz = (IHzz −IHxx ) sin (ν) cos (ν)−IHxz cos (2ν). (4.131) Similarly, the vertical distances zBand zHhave been taken positive in Table 4.1. Note that the axes of the body reference frames of bodies Band H,⟨XBYBZB⟩ and ⟨XHYHZH⟩, respectively, are written with capital letters in Fig. 4.2, whereas lower case letters have been used for the geometric parameters xB,zB,xHand zHof Table 4.1. Next, the absolute position vectors of the centres of mass rGiand the contact points, rPand rQ, can be obtained as: rGB=xbybzbT,rI=rGB+RB¯rGBI, rGR=rGB+RB¯rGBGR,rGH=rI+RH¯rIGH, rGF=rI+RH¯rIGF,rP=rGR+RRRξR¯ρR,(4.132) rQ=rGF+RFRξF¯ρF, where ¯rGBI,¯rGBGR,¯rIGH,¯rIGF,¯ρRand ¯ρFare the following vectors, expressed in the corresponding body reference frames: ¯rGBI=xI0 0 T,¯rGBGR=−xB0ρR−zBT, ¯rIGH=a0−bT,¯rIGF=d0−eT,(4.133) 119 Chapter 4. Novel linearization approaches for constrained multibody systems Figure 4.3: Lateral view of the bicycle benchmark with the main geometric parameters. ¯ρR=ρR0 0 T,¯ρF=ρF0 0 T. The auxiliar point Iand the dimensions xI,a,b,d,ein Eqs. (4.133) are shown in Fig. 4.3. These magnitudes, in terms of the geometric parameters of Table 4.1, are given by: xI=w+c−xB−zBtan (ν),(4.134) a= (xH−w−c) cos (ν) + zHsin (ν),(4.135) b=1 2 cos (ν)((xH−w−c) sin (2ν)−zHcos (2ν)+2zB−zH),(4.136) d=RFsin (ν)−ccos (ν),(4.137) e=1 2 cos (ν)(2zB−RF−csin (2ν)−RFcos (2ν)) .(4.138) On the other hand, four holonomic and four nonholonomic constraints are considered in the multibody model. The holonomic constraints arise from the contact between the bicycle and the ground. Thus, the Z-component of the position vectors of the contact points, rPZand rQZ, must be zero, and the tangent vectors to these contact points (denoted as tPand tQ) must be orthogonal to the normal to the ground n. Therefore, the set of holonomic constraints is: C(x) = rPZrQZn·tPn·tQT=04×1,(4.139) with the tangent vectors tPand tQ: tR=RR ∂¯rGRP ∂ξR ,tF=RF ∂¯rGFQ ∂ξF .(4.140) 120 4.3. Validation of the linearization approaches with the bicycle benchmark Table 4.1: List of the main geometric and inertial parameters of the bicycle multibody model and their numerical values, taken from the benchmark by Meijaard et al. [168]. Parameter [unit] Symbol Numerical value Wheel base [m] w1.02 Trail [m] c0.08 Steer axis tilt [rad / ◦]νπ 10 /18 Horizontal distance P,GB [m] xB0.3 Vertical distance P,GB[m] zB0.9 Horizontal distance P,GH [m] xH0.9 Vertical distance P,GH[m] zH0.7 Rear body and frame assembly (B) Mass [kg] mB85 with rider’s mass Moment of inertia tensor [kg m2] ¯ IB=  ¯ IBxx 0¯ IBxz 0¯ IByy 0 ¯ IBxz 0¯ IBzz    9.2 0 −2.4 0 11 0 −2.4 0 2.8  Front handlebar and fork assembly (H) Mass [kg] mH4 Moment of inertia tensor [kg m2] ¯ IH=  ¯ IHxx 0¯ IHxz 0¯ IHyy 0 ¯ IHxz 0¯ IHzz    5.841 0 −0.912 060 −0.912 0 0.759  10−2 from Eqs. (4.131) Rear wheel (R) Radius [m] ρR0.3 Mass [kg] mR2 Moment of inertia tensor [kg m2] ¯ IR=  ¯ IRxx 0 0 0¯ IRyy 0 0 0 ¯ IRxx    0.0603 0 0 0 0.12 0 0 0 0.0603  Front wheel (F) Radius [m] ρF0.35 Mass [kg] mF3 Moment of inertia tensor [kg m2] ¯ IF=  ¯ IFxx 0 0 0¯ IFyy 0 0 0 ¯ IFxx    0.1405 0 0 0 0.28 0 0 0 0.1405  121 Chapter 4. Novel linearization approaches for constrained multibody systems The linear system of equations (4.16) obtained with the first symbolic procedure can be particularized for a relative equilibrium as the one defined in Eq. (4.164), which verifies ˙x0 ai(t) = V0 ai and ¨x0 ai(t) = 0. Therefore, the linearized equations of motion (4.16) become: m0¨ ˜xai =∂f ∂xai 0 ˜xai +∂f ∂xdd 0 ˜xdd +∂f ∂˙xai 0 ˙ ˜xai,(4.166) ˙ ˜xdd =H0˙ ˜xai +∂HV 0 ai ∂xai 0 ˜xai +∂HV 0 ai ∂xdd 0 ˜xdd.(4.167) Defining ˜ X=˜xai ˙ ˜xai ˜xdd T, Eqs. (4.166)-(4.167) can be written as a firstorder ODE system ˙ ˜ X=J˜ X, with the Jacobian matrix: J=       0(n−m−l)I(n−m−l)0(n−m−l)×(m+l) m−1 0 ∂f ∂xai 0 m−1 0 ∂f ∂˙xai 0 m−1 0 ∂f ∂xdd 0 ∂HV 0 ai ∂xai 0 H0 ∂HV 0 ai ∂xdd 0        .(4.168) Observe that, in this case, the size of Jis (2n−m−l)×(2n−m−l) = 14 ×14. On the other hand, the use of the second numerical approach leads to the Jacobian matrix: J=    0(n−m−l)I(n−m−l)0(n−m−l)×(m+l) m−1 0R0¯ ¯ V0+S0Eai m−1 0R0T0m−1 0R0¯ ¯ V0+S0Edd ¯ V0Eai Tdd x0¯ V0Edd    . (4.169) The substitution of the numerical values of the bicycle geometric and dynamic parameters in the Jacobian matrices (4.168) and (4.169) leads to: 128 4.3. Validation of the linearization approaches with the bicycle benchmark J=                           0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9.4898 0.8912 v2+ 0.5715 0 −0.1055 v0.3305 v0 0 0 0 0 0 0 0 0−11.7195 30.9088 −1.9717 v20−3.6768 v−3.0849 v0 0 0 0 0 0 0 0 0 0 0.2797 v0−0.9 0.0224 0 0 v00000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9324 v0 0 0.0746 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3.3333 0 0 0 0 0 0 0 0 0 0 0 0 0 2.8571 0 0 0 0 0 0 0 0 0 0 0 0 0 −3.3333 0 0 0 0 0 0 0 0 0 0 0 0 0 −2.8571 0 0 0 0 0 0 0 0 0 0                           .(4.170) 129 Chapter 4. Novel linearization approaches for constrained multibody systems Note that, as in Eq. (4.159), the Jacobian matrix (4.170) is obtained as a function of the forward velocity v. The Jacobian matrices derived from the first symbolic approach (Eq. (4.168)) and the second numerical approach (Eq. (4.169)) coincide in the first 32 numerical digits. This demonstrates the accuracy of the numerical procedure. Equation (4.170) shows the result with four decimal places because of space limitations. Note that the Jacobian matrix (4.159) contains linear terms in the forward speed v, while Eq. (4.170) exhibits quadratic terms. The Jacobian matrix (4.170) does not depend explicitly on time, and presents constant coefficients. The characteristic equation is: q(λ) = λ10 λ4+ζ1vλ3+ζ2v2−ζ3λ2+ζ4v2−ζ5vλ +ζ6−ζ7v2= 0, (4.171) where the coefficients of the quartic polynomial ζkin Eq. (4.171), with k= 1 . . . 7, coincide with those given in Eq. (4.161). Given that Jis 14 ×14, Eq. (4.171) provides fourteen eigenvalues, among which ten are zero. This number of null eigenvalues arises from the two eigenvalues corresponding to the forward motion, and m+l= 8 spurious null eigenvalues, associated with the time derivative of the m= 4 holonomic constraints and the l= 4 nonholonomic constraints. As noted in Sect. 4.3.3, the four remaining nonzero eigenvalues arise from the roots of the quartic polynomial of Eq. (4.171), associated with the lean and steering angles. Concerning the first symbolic approach developed in Sect. 4.1.1, the main advantages are: •The Lagrange multipliers Λare eliminated and one works with the reduced equations of motion. •The Jacobian matrix (4.17) is (2n−m−l)×(2n−m−l), which leads to m+lfewer spurious eigenvalues than the first numerical procedure. This results from the fact that, in the first numerical approach, the presence of m+ldependent coordinates leads to 2 (m+l) spurious eigenvalues, while in the first symbolic approach, the dependency between the coordinates is considered by means of the m+lconstraints, which implies only m+l spurious eigenvalues. On the other hand, the first symbolic approach presents the following disadvantages: •When working with moderately complex multibody systems, the symbolic computation of the inverse matrix (Ddd (x))−1in Eq. (2.131) cannot be performed by most of personal computers. Denoting a generic element of the matrix Ddd (x) as Dij (x), in the particular case of the bicycle, this obstacle has been overcome in this study by doing dij =Dij (x), with dij being constant. Thus, the matrix inversion is carried out using the elements dij, which significantly reduces the computational cost, and then the elements Dij (x) are back-substituted. 130 4.3. Validation of the linearization approaches with the bicycle benchmark •In addition, the construction of the time derivative of the transformation matrix, ˙ T, in the term fof Eq. (4.4) involves a high computational cost, since it entails the computation of the time derivative of an inverse matrix. •Lastly, computing the partial derivatives in the Jacobian matrix (4.17) is computationally expensive and not feasible for complex multibody systems. All these drawbacks make the second numerical linearization approach (developed in Sect. 4.2.2) more appropriate, being its principal advantages: •The approach provides the exact linearized equations of motion, obtaining by means of Eq. (4.95) the same Jacobian matrix of Eq. (4.17). The results coincide in the first 32 numerical digits. •It is much more efficient than the first symbolic approach, with a substantial reduction in the computational cost. This results from the fact that the computation of the inverse matrix in Tx0, defined in Eq. (4.66), is performed after the evaluation for the reference solution x0. Furthermore, the numerical procedure avoids the computation of derivatives of inverse matrices, as occurred in the symbolic case with the matrix ˙ T. As disadvantages, the Lagrange multipliers in the reference motion Λ0are required and the construction of the Jacobian matrix (4.95) entails the computation of a larger number of derivatives with respect to the first symbolic approach. 4.3.5 Use of the second symbolic procedure and third numerical procedure Lastly, the second symbolic approach (Sect. 4.1.2) and the third numerical approach (Sect. 4.2.3) are used. According to the coordinate partition presented in Eq. (2.137), the following sets of coordinates are considered: xai =xbϕbδT, xad =ybψbθRθFT, xd=zbθbξRξFT. (4.172) Note that the independent coordinates in Eq. (4.172) are the same as those of Eq. (4.162), corresponding to the forward motion, the lean and steering angles. Similarly, the set of dependent coordinates xdd in Eq. (4.162) is splitted in the admissible dependent xad and dependent xdsets. In the third numerical approach, the same coordinate partition is considered for the variations of the coordinates, as shown in Eq. (4.96), and therefore: ˜xai =˜xb˜ ϕb˜ δT, ˜xad =˜yb˜ ψb˜ θR˜ θFT,(4.173) 131 Chapter 4. Novel linearization approaches for constrained multibody systems ˜xd=˜z2˜ θ2˜ ξR˜ ξFT. The steady forward motion (4.148) can be expressed, according to this coordinate partition, as: x0=x0 ai x0 ad x0 dT, x0 ai(t) = P0 ai +V0 ait, x0 ad(t) = P0 ad +V0 adt, x0 d(t) = P0 d+V0 dt, (4.174) where P0 ai,V0 ai,P0 ad,V0 ad,P0 dand V0 dare constant vectors, given by: P0 ai =03×1,V0 ai =v0 0 T, P0 ad =04×1,V0 ad =0 0 v RR ,v RFT ,(4.175) P0 d=zB0π 2ν+π 2T,V0 d=0 0 −v ρR−v ρFT . First, the second symbolic procedure is used. The linear system of equations (4.40)-(4.41) can be particularized for the relative equilibrium defined in Eq. (4.174), which verifies ˙x0 ai(t) = V0 ai and ¨x0 ai(t) = 0. Therefore, the linearized equations of motion (4.16) become: m0¨ ˜xai = ∂f ∂xai −∂f ∂xd∂C ∂xd−1∂C ∂xai !0 ˜xai + ∂f ∂xad −∂f ∂xd∂C ∂xd−1∂C ∂xad !0 ˜xad +∂f ∂˙xai 0 ˙ ˜xai,(4.176) ˙ ˜xad =G0˙ ˜xai + ∂GV 0 ai ∂xai −∂GV 0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 ˜xai + ∂GV 0 ai ∂xad −∂GV 0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 ˜xad.(4.177) Defining ˜ X=˜xai ˙ ˜xai ˜xad T, Eqs. (4.176)-(4.177) can be written as a first-order ODE system ˙ ˜ X=J˜ X, with the Jacobian matrix: J=   0(n−m−l)I(n−m−l)0(n−m−l)×l J21 J22 J23 J31 J32 J33   ,(4.178) 132 4.3. Validation of the linearization approaches with the bicycle benchmark where J21 =m−1 0 ∂f ∂xai −∂f ∂xd∂C ∂xd−1∂C ∂xai !0 , J22 =m−1 0 ∂f ∂˙xai 0 , J23 =m−1 0 ∂f ∂xad −∂f ∂xd∂C ∂xd−1∂C ∂xad !0 , J31 = ∂GV 0 ai ∂xai −∂GV 0 ai ∂xd∂C ∂xd−1∂C ∂xai !0 , J32 =G0, J33 = ∂GV 0 ai ∂xad −∂GV 0 ai ∂xd∂C ∂xd−1∂C ∂xad !0 . (4.179) On the other hand, the third numerical linearization approach leads to the Jacobian matrix of Eq. (4.124). The substitution of the numerical values of the bicycle geometric and dynamic parameters in the Jacobian matrices of Eqs. (4.178) and (4.124) results in: 133 Chapter 4. Novel linearization approaches for constrained multibody systems J=                  0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9.4898 0.5715 + 0.8912 v20−0.1055 v0.3305 v0 0 0 0 0−11.7195 30.9088 −1.9717 v20−3.6768 v−3.0849 v0 0 0 0 0 0 0.2797 v0−0.9 0.0224 0 v0 0 0 0 0.9324 v0 0 0.0746 0 0 0 0 0 0 0 3.3333 0 0 0 0 0 0 0 0 0 2.8571 0 0 0 0 0 0                  .(4.180) 134 4.3. Validation of the linearization approaches with the bicycle benchmark As in Sect. 4.3.4, the Jacobian matrices obtained with the second symbolic approach and the third numerical procedure coincide in the first 32 numerical digits. The Jacobian matrix of Eq. (4.180) is (2n−2m−l)×(2n−2m−l) = 10 ×10. As Eq. (4.170), it does not explicitly depend on time and contains quadratic terms with the forward speed v. Its characteristic equation is: h(λ) = λ6λ4+ζ1vλ3+ζ2v2−ζ3λ2+ζ4v2−ζ5vλ +ζ6−ζ7v2= 0, (4.181) where the coefficients of the quartic polynomial ζkin Eq. (4.181), with k= 1 . . . 7, coincide with those given in Eqs. (4.161) and (4.171). Equation (4.181) provides 2n−2m−l= 10 eigenvalues, among which six are zero. In this way, four out of the six null eigenvalues are spurious, associated with the l= 4 linearized nonholonomic constraints, and the other two correspond to the forward motion. Note that the elimination of the m= 4 holonomic constraints results in the reduction of the number of spurious null eigenvalues from ten (see Eq. (4.171)) to six. The second symbolic approach presents the same drawbacks of the first symbolic procedure, described in Sect. 4.1.1. The construction of fand Gin Eqs. (4.18)- (4.20), respectively, involves the inversion of matrices that, for moderately large multibody systems, cannot be addressed by the majority of personal computers. Furthermore, the computation of frequires the performance of time derivatives of inverse matrices, which is computationally expensive. Finally, obtaining the derivatives of these terms with respect to the generalized coordinates requires large amounts of time. Therefore, despite one obtains the maximum possible reduction of the linearized equations of motion, this approach might be unaffordable. In contrast, the third numerical procedure avoids these problems, since it performs the inversion of matrices evaluated for the reference solution, greatly simplifying the proccess. Moreover, no derivatives of inverse matrices are required, resulting in much greater efficiency. The procedure is completely accurate, obtaining the same Jacobian matrix (first 32 numerical digits) as in the symbolic counterpart approach. On the negative side, the construction of the Jacobian matrix in Eq. (4.124) involves more steps and derivatives than the Jacobian matrix derived with the second symbolic approach, shown in Eq. (4.42). 4.3.6 Comparison with the benchmark results Meijaard et al. [168] wrote the linearized lean and steer equations of the bicycle as: M¨q+vC1˙q+gK0+v2K2q=f,(4.182) where q=ϕ δ Tis the set of generalized coordinates, including the lean and steering angles; f=TϕTδTrepresents the lean disturbance and steering torques exerted by the rider; and M,C1,K0and K2correspond to the inertia, velocity-dependent damping, gravitational stiffness and gyroscopic and centrifugal effects matrices, respectively. The detailed expressions of these matrices can be 135 Chapter 4. Novel linearization approaches for constrained multibody systems found in Ref. [168]. The mass matrix Mis a symmetric positive definite matrix; the velocity-dependent matrix C1presents a gyroscopic part (skew-symmetric) and a symetric part; the gravity-dependent stiffness matrix K0is symmetric; and the velocity-dependent stiffness matrix K2presents a symmetric part, corresponding to the centrifugal forces, and a non-conservative positional part (skewsymmetric). In this case, the symmetric part of the velocity-dependent matrix C1 and the skew-symmetric part of the velocity-dependent stiffness matrix K2are associated with the nonholonomic nature of the system. The equation corresponding to the forward motion is decoupled from the system of equations (4.182). The linear stability of the uncontrolled bicycle will be considered, which corresponds to f=0. As stated in [168], given that the eigenvalues are independent of coordinate choice and equation form, they allow the comparison between different approaches. Thus, the characteristic equation associated with Eq. (4.182) is: det Mλ2+vC1λ+gK0+v2K2= 0,(4.183) which can be written as the following quartic polynomial: p0(λ) = λ4+ζ∗ 1vλ3+ζ∗ 2v2−ζ∗ 3λ2+ζ∗ 4v2−ζ∗ 5vλ +ζ∗ 6−ζ∗ 7v2.(4.184) Note that the structure of the polynomial in Eq. (4.184) is the same as that of the quartic polynomials obtained in Eqs. (4.160), (4.171) and (4.181). The coefficients of Eq. (4.184) incorporate the superscript ∗to indicate that they correspond to the bicycle benchmark. Substituting the numerical values of the geometric and dynamic parameters of Ref. [168] in Eq. (4.183), one obtains the numerical values of the coefficients in Eq. (4.184). Table 4.2: Comparison of the characteristic equations coefficients. Characteristic equations coefficients Benchmark Present work ζ∗ 1= 3.190 387 977 238 80 ζ1= 3.190 387 977 238 81 ζ∗ 2= 3.512 478 401 483 87 ζ2= 3.512 478 401 483 87 ζ∗ 3= 40.398 527 840 014 26 ζ3= 40.398 527 840 014 28 ζ∗ 4= 3.484 820 713 307 97 ζ4= 3.484 820 713 307 98 ζ∗ 5= 26.561 398 462 011 76 ζ5= 26.561 398 462 011 81 ζ∗ 6= 300.015 050 749 111 96 ζ6= 300.015 050 749 112 22 ζ∗ 7= 8.266 760 150 418 84 ζ7= 8.266 760 150 418 82 Table 4.2 shows the numerical values of the coefficients of the characteristic polynomial equations, obtained with the linearization procedures of the present work, and those of Meijaard et al. [168]. The results are shown with fourteen decimal places. One can see that the first fourteen significant digits coincide, proving the accuracy of the linearization approaches presented in this study. 136 4.3. Validation of the linearization approaches with the bicycle benchmark The computation of the eigenvalues of the uncontrolled bicycle shows that it is statically unstable at low speeds, asymptotic stable for medium speeds and oscillatory unstable at high speeds. Figure 4.4 shows the root locus of the bicycle benchmark, for the velocity range 0 < v < 10 m/s. The eigenvalues corresponding to v= 0 m/s are highlighted with crosses, and the eigenvalues associated with v= 10 m/s are represented with squares. Figure 4.5 shows the evolution of the eigenvalues in the forward speed range 0 < v < 10 m/s, obtaining the same continuation diagram as Meijaard et al. [168]. The real parts of the eigenvalues have been plotted with solid lines, while the imaginary parts are shown with dashed line. It can be seen that, for the forward speed vd, the oscillatory unstable weave motion appears after the splitting of the two real eigenvalues λdinto a complex conjugate pair. Once the weave speed vwis reached, the complex conjugate pair traverses the imaginary axis at λw, so the weave mode is stabilized. The asymptotic stability of the bicycle benchmark continues until the capsize speed vc, when the capsize real eigenvalue becomes positive, resulting in the instability of the bicycle. In this way, the asymptotic stability is achieved for the range vw< v < vc. The points corresponding to vwand vcare also highlighted in the root locus of Fig. 4.4. The numerical values of the speeds vd,vwand vcare shown in Table 4.3, being the corresponding eigenvalues λdand λwequal to those obtained in Ref. [168] up to the fourteenth significant digit. Therefore, it is concluded that the proposed linearization procedures are completely accurate. Figure 4.4: Root locus of the bicycle benchmark model, for the velocity range 0< v < 10 m s. 137 Chapter 4. Novel linearization approaches for constrained multibody systems (a) Three-dimensional view of the toroidal wheel (b) Cross-section of the toroidal wheel Figure 4.9: Model of the toroidal wheel: angular coordinates and geometric parameters. where µiis the torus aspect ratio, given by the ratio of the minor to the major radius of the wheel, and σiis the aspect ratio of the elliptical cross-section, defined as the ratio of the minor to the major radius of the ellipse. Note that the hoopshaped wheel case of Meijaard et al. [168] is also captured with this model, and is obtained by degenerating the tori for µi=σi= 0, leading to ρi=Ri. The inertia tensors of the rear and front wheels used in Ref. [168] are substituted by the inertia tensors of an elliptic torus, parameterized in terms of the cross-section parameters ρi,aiand bi, whose expressions can be found in Diaz et al. [266]: Iixx =Iizz =1 8mi4ρ2 i+ 3a2 i+ 2b2 i,(4.189) Iiyy =1 4mi4ρ2 i+ 3a2 i,(4.190) with i={R, F}. The set of holonomic constraints, which arise from the contact of the bicycle wheels with the ground, is given by: C(x) = rPZrQZn·tLRn·tTRn·tLFn·tTFT=0,(4.191) where rPZand rQZare the Z-components of the position vectors of the contact points, and were already defined in Eq. (4.139); tLiand tTiare the longitudinal and transversal tangent vectors to the contact points; and nis the normal vector to the ground surface. These vectors are given by: rP=rGR+RR¯rGRP,rQ=rGF+RF¯rGFQ, 144 4.5. The role of the dependent coordinates in the linearization tLR=RR ∂¯rGRP ∂ξR ,tTR=RR ∂¯rGRP ∂ηR ,(4.192) tLF=RF ∂¯rGFQ ∂ξF ,tTF=RF ∂¯rGFQ ∂ηF . In Eqs. (4.192), RRand RFare the orientation matrices of bodies Rand F, respectively. Moreover, ¯rGRPand ¯rGFQare position vectors, expressed in the body reference frames of bodies Rand F, respectively, computed as follows: ¯ rGRP=RξR(ρR+RηRrR),¯ rGFQ=RξF(ρF+RηFrF), ρR=ρR0 0 T,ρF=ρF0 0 T,(4.193) rR=aR0 0 T,rF=aF0 0 T, where Rξiand Rηiare rotation matrices, associated with the coordinates ξiand ηi. The nonholonomic constraints arise from the assumption of rolling without slipping, and are computed as in Eq. (4.145): Cnh (x,˙x) = vPxvPyvQxvQyT=0.(4.194) In Eq. (4.194), vPand vQare the velocity of the contact points Pand Q, respectively, computed as follows: vP=vGR+RR(¯ ωRׯ rGRP),vQ=vGF+RF(¯ ωFׯ rGFQ),(4.195) where vGRand vGFare the absolute velocities of GRand GF, and ¯ωRand ¯ωF are the angular velocities of the rear and front wheels, respectively, expressed in the body reference frames of bodies Rand F. Note that, in this case, n= 13, m= 6 and l= 4, and as in the hoop-shaped wheel case, the number of degrees of freedom of the system is ng= 3. The equations of motion of the bicycle with toroidal wheels are given by the index-3 DAE of Eq. (4.147). The steady forward motion of the system is given by: x0=x0 by0 bz0 bψ0 bϕ0 bθ0 bδ0θ0 Rθ0 Fξ0 Rξ0 Fη0 Rη0 FT,(4.196) with x0 b=vt, ψ0 b= 0, y0 b= 0, ϕ0 b= 0, z0 b=zB, θ0 b= 0, δ0= 0, θ0 R=v RR t, (4.197) θ0 F=v RF t, ξ0 R=π 2−θ0 R, ξ0 F=ν+π 2−θ0 F, η0 R= 0, η0 F= 0. 145 Chapter 4. Novel linearization approaches for constrained multibody systems 4.5.2 The role of the dependent coordinates in the linear stability analysis of the bicycle The role of the dependent coordinates in the linearization is illustrated with the bicycle model described in Sect. 4.5.1. Groups 1 and 2: as already shown, the linearization approaches described in Sects. 3.1 and 3.2 (devoted to the linearization of the index-3 and index-1 DAE equations of motion, respectively) make use of the redundant set of coordinates. In these approaches, no coordinate partition in terms of independent and dependent coordinates is considered, and no distinction in the handling of the coordinates is made. The resulting Jacobian matrices include the partial derivatives with respect to the complete set of coordinates. Among the linearization approaches of Group 1, the procedure of Pappalardo et al. [219, 220], summarized in Sect. 3.1, was used to analyse the stability of the steady forward motion of the bicycle with hoop-shaped wheels. In the same way, the first numerical procedure proposed in this PhD thesis, exposed in Sect. 4.2.1, was used to linearize the index-1 DAE equations of motion of the bicycle multibody model along the reference solution (4.148). Group 3: the approaches of the third group make use of a coordinate partition in terms of independent and dependent coordinates. In this case, it is important to clarify the handling of the dependent coordinates in the linearization. Consider the vector of coordinates (4.185). The following coordinate partition is chosen: ˜xai =˜xb˜ ϕb˜ δT, ˜xdd =˜yb˜zb˜ ψb˜ θb˜ θR˜ θF˜ ξR˜ ξF˜ηR˜ηFT, (4.198) where n−m−l= 3 independent coordinates (the same number as degrees of freedom of the system) are considered. The forward motion, the lean and steering angles have been chosen as independent coordinates. This partition is possible, since the submatrix Ddd (x) in Eq. (2.131) is non-singular. The Jacobian matrices (4.17) and (4.95) are computed for the steady forward motion of the bicycle, leading to: 146 4.5. The role of the dependent coordinates in the linearization J=                              0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0000000 0α1α2+α3v20α4v α5v00 0 0000000 0α6α7+α8v20α9v α10 v00 0 0000000 0 0 α11 v0α12 α13 0 0 α14v0000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 α15 v0 0 α16 00 0 0000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 α17 0 0 00 0 0000000 0 0 0 α18 0 0 00 0 0000000 0 0 0 −α17 0 0 00 0 0000000 0 0 0 −α18 0 0 00 0 0000000 0 0 0 0 α19 0 00 0 0000000 0 0 0 0 α19 α20 00 0 0000000                              ,(4.199) 147 Chapter 4. Novel linearization approaches for constrained multibody systems where the coefficients αk, with k= 1 . . . 20, are functions of the geometric and dynamic parameters of the bicycle. The computation of the eigenvalues associated with the Jacobian matrix (4.199), for the hoop-shaped wheels case (µi= 0), yields the velocity continuation digram of Fig. 4.5. Moreover, ten spurious null eigenvalues are obtained, associated with the m+l= 10 dependent coordinates in Eq. (4.198). The linear stability results of the toroidal case will be shown in Sect. 5.4. With the Jacobian matrix (4.199), the following question arises: can the Jacobian matrix (4.199) be reduced by eliminating the blocks associated with the dependent coordinates? Note that, in the Jacobian matrices (4.17) and (4.95), the blocks including the derivatives with respect to the dependent coordinates are J23 =m−1 0 ∂f−m¨x0 ai ∂xdd 0 =m−1 0R0¯ ¯ V0+S0Edd. In the particular case of the bicycle, with the selection of independent coordinates of Eq. (4.198), it can be seen that these blocks correspond to null submatrices (highlighted with blue in Eq. (4.199)), and thus the dependent coordinates can be ignored. In this scenario, the Jacobian matrices (4.17) and (4.95) can be simplified to: J=   0(n−m−l)I(n−m−l) ∂f−m¨x0 ai ∂xai 0 m−1 0 ∂f ∂˙xai 0   = 0(n−m−l)I(n−m−l) m−1 0R0¯ ¯ V0+S0Eai m−1 0R0T0!. (4.200) In the bicycle example, Eq. (4.200) leads to: J=          0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0α1α2+α3v20α4v α5v 0α6α7+α8v20α9v α10 v          .(4.201) In the linear stability analysis of the bicycle performed by Escalona et al. [171], it is expressly mentioned that the dependent coordinates of Eq. (4.198) can be considered as null in the linearization, and therefore they can be ignored. This is demonstrated by obtaining the null submatrix in the block associated with the dependent coordinates, as shown in Eq. (4.199). The linearization approaches by Cuadrado et al. [225, 226] and Gonz´alez et al. [218] described in Sect. 3.3, which do not include the blocks associated with the dependent coordinates, also lead to the Jacobian matrix (4.201). Despite the reduction of the Jacobian matrix (4.199) to (4.201), all the relevant stability information is retained. The same velocity continuation diagram of Fig. 4.5 is obtained. The Jacobian matrix (4.199) is (2n−m−l)×(2n−m−l) = 16 ×16, while the Jacobian matrix (4.201) is 148 4.5. The role of the dependent coordinates in the linearization (2 (n−m−l)) ×(2 (n−m−l)) = 6 ×6, leading to the elimination of the ten spurious null eigenvalues obtained with Eq. (4.199). Note that the third row of the Jacobian matrix (4.201) is null. This row corresponds to the equation of the bicycle’s forward motion, ¨ ˜xb= 0, which is decoupled from the others. Therefore, the Jacobian matrix (4.201) can be further reduced. The fifth and sixth rows of Eq. (4.201) correspond to the linearized lean and steer equations, which are equivalent to those derived by Meijaard et al. [168]. Defining ˜q=ϕ δ T, the linearized lean and steer equations are given by: ¨ ˜q=vα4α5 α9α10 ˙ ˜q+α1α2 α6α7˜q+v20α3 0α8˜q.(4.202) Equations (4.202) contain all the information about the linear stability of the lateral motion of the bicycle. Nevertheless, the dependent coordinates cannot generally be ignored. Consider now the following coordinate partition: ˜xai =˜xb˜ ψb˜ ϕbT, ˜xdd =˜yb˜zb˜ θb˜ δ˜ θR˜ θF˜ ξR˜ ξF˜ηR˜ηFT, (4.203) where the forward motion, the yaw and lean angles are selected as independent coordinates. The coordinate partition of Eq. (4.203) is also possible, since the submatrix Ddd (x) in Eq. (2.131) is non-singular. The Jacobian matrices (4.17) and (4.95), using the coordinate partition (4.203), become: 149 Chapter 4. Novel linearization approaches for constrained multibody systems J=                              0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 000 0 000000 0 0 β10β2v β3v000 β4+β5v2000000 0 0 β60β7v β8v000β9+β10v2000000 0β11v0 0 β12 β13 000 0 000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 β14 0 0 0 0 β15v000000 0 0 0 β16 0 0 000 0 000000 0 0 0 β17 0 0 000 0 000000 0 0 0 −β16 0 0 000 0 000000 0 0 0 −β17 0 0 000 0 000000 0 0 0 0 0 β18 000 0 000000 0 0 0 0 β19 β20 000 β21 v000000                              .(4.204) 150 4.5. The role of the dependent coordinates in the linearization In Eq. (4.204), the coefficients βk, with k= 1 . . . 21, are functions of the geometric and dynamic parameters of the bicycle. Despite a different coordinate partition to that of Eq. (4.198) is used, the Jacobian matrix (4.204) yields the same velocity continuation diagram of Fig. 4.5. Moreover, as in Eq. (4.199), ten spurious null eigenvalues are obtained. The following question is now posed: can the Jacobian matrix (4.204) be reduced by eliminating the blocks associated with the dependent coordinates? Note that, in this case, this block is a non-null submatrix, highlighted with blue in Eq. (4.204). If the derivatives with respect to the dependent coordinates are ignored in Eq. (4.204), the following reduced Jacobian matrix is obtained: J=          0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 β10β2v β3v 0 0 β60β7v β8v          .(4.205) The resulting velocity continuation diagram associated with the Jacobian matrix (4.205) does not correspond to that of Fig. 4.5, obtaining completely incorrect results. Due to the symmetry of the bicycle system, the selection of independent coordinates in Eq. (4.198) allowed obtaining the correct linearized equations of motion despite ignoring the dependent coordinates. Nevertheless, with the coordinate partition of Eq. (4.203), the dependent coordinates cannot be ignored and they are required in the computation of the exact linearized equations of motion. Therefore, the Jacobian matrix (4.204) cannot be reduced to Eq. (4.205). Group 4: the use of the approaches described in Sect. 3.4 allows the elimination of mdependent coordinates. In the particular case of the bicycle multibody model, m= 6. This is accomplished by using Eq. (4.105), which makes use of the linearized holonomic constraints to express the variations of the dependent coordinates ˜xdin terms of ˜xai and ˜xad. Choosing as independent coordinates the forward motion, the lean and steering angles, a possible coordinate partition, based on Eq. (2.137), is: ˜xai =˜xb˜ ϕb˜ δT, ˜xad =˜yb˜ ψb˜ θR˜ θFT,(4.206) ˜xd=˜z2˜ θ2˜ ξR˜ ξF˜ηR˜ηFT. The use of Eq. (4.124), considering the coordinate partition (4.206), leads to the 151 Chapter 4. Novel linearization approaches for constrained multibody systems following Jacobian matrix: J=                  0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0γ1γ2+γ3v20γ4v γ5v0 0 0 0 0γ6γ7+γ8v20γ9v γ10 v0 0 0 0 0 0 γ11 v0γ12 γ13 0γ14v0 0 0 0 γ15v0 0 γ16 0 0 0 0 0 0 0 γ17 0 0 0 0 0 0 0 0 0 γ18 0 0 0 0 0 0                  ,(4.207) where γk, with k= 1 . . . 18, are functions of the geometric and dynamic parameters of the bicycle. The eigenvalues associated with Eq. (4.207) correspond to those of the velocity continuation diagram of Fig. 4.5. In this case, the Jacobian matrix (4.207) is (2 (n−m)−l)×(2 (n−m)−l) = 10 ×10, and the number of spurious null eigenvalues is reduced from ten to four, associated with the l= 4 linearized nonholonomic constraints. In Eq. (4.207), the dependent coordinates ˜xdof Eq. (4.206) are eliminated with the use of the linearized holonomic constraints. Moreover, the block associated with the coordinates ˜xad is a null submatrix, highlighted with color blue. This block corresponds to J23 in the Jacobian matrices (4.42) and (4.124), obtained with the second symbolic and third numerical approaches, respectively. As in Eq. (4.199), the Jacobian matrix (4.207) can be reduced to Eq. (4.201). This reduced Jacobian matrix leads to the correct eigenvalues of the system. Nevertheless, as in the approaches of Group 3, the reduction ignoring the coordinates ˜xad cannot always be done. Consider the following coordinate partition: ˜ xai =˜xb˜ ψb˜ ϕbT, ˜ xad =˜yb˜ δ˜ θR˜ θFT,(4.208) ˜ xd=˜zb˜ θb˜ ξR˜ ξF˜ηR˜ηFT, where the forward motion, the yaw and lean angles are selected as independent 152 4.5. The role of the dependent coordinates in the linearization coordinates. In this case, the Jacobian matrices (4.42) and (4.124) are given by: J=                  0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 δ10δ2v δ3v0δ4+δ5v20 0 0 0 δ60δ7v δ8v0δ9+δ10v20 0 0δ11v0 0 δ12 δ13 0 0 0 0 0 0 0 0 δ14 0 0 δ15v0 0 0 0 0 δ16 0 0 0 0 0 0 0 0 0 δ17 0 0 0 0 0 0                  .(4.209) In the Jacobian matrix (4.209), the coefficients δk, with k= 1 . . . 17, are functions of the geometric and dynamic parameters of the bicycle. In contrast to Eq. (4.207), the block associated with the coordinates ˜xad is a non-null submatrix, highlighted with color blue in Eq. (4.209). The Jacobian matrix (4.209) leads to the same velocity continuation diagram of Fig. 4.5. However, eliminating this block results in completely incorrect results, as in Eq. (4.205). This analysis clearly shows that the dependent coordinates need to be carefully handled in the linearization of the equations of motion of constrained multibody systems. If one of the approaches of Group 1 (Sect. (3.1)) and Group 2 (Sect. (3.2)) is used, no distinction between coordinates is considered and the handling of the dependent coordinates does not require further attention. As a drawback, a bulky Jacobian matrix is obtained. In the approaches of Group 3 (Sect. (3.3)), it is shown that, in general, the dependent coordinates cannot be ignored. In the particular case of the bicycle, the selection of the forward motion, the lean and steering angles as independent coordinates allows ignoring the derivatives with respect to the dependent coordinates, due to the symmetry of the bicycle system along the steady forward motion. However, a different admissible choice of independent coordinates, as the one shown in Eq. (4.203), inhibits ignoring the derivatives with respect to the dependent coordinates. Lastly, by using the linearization approaches of Group 4 (Sect. (3.4)), the linearized holonomic constraints can be used to eliminate mdependent coordinates. In the particular case of a multibody system with only holonomic constraints, this allows obtaining the linearized equations of motion only in terms of the n−mindependent coordinates and their time derivatives, eliminating the dependent ones. Nevertheless, in a multibody system with nonholonomic constraints as the bicycle, the linearized nonholonomic constraints cannot be used to express the remaining dependent coordinates as a function of the independent ones, since these constraints are non-integrable. The approaches of Group 4 lead to the maximum possible reduction of the linearized equations of motion. As shown with the bicycle, in some multibody systems with a low number of coordinates and in presence of symmetries, a selection of independent coordinates can be intuitively found to ignore the dependent coordinates. However, this is not the usual case, and the sys- 153 Chapter 6. Linear stability analysis of nonholonomic multibody systems: study of linear periodic time-varying systems λ2 cof Eq. (6.130). Since |ρ2 c|>1, the periodic solution is unstable. Next, the monodromy matrix ¯ X(T) is obtained for ϕ0=−5◦, ¯ω= 0.75 and ¯ Ω≃0.0794, which leads to the Floquet multipliers: ρ1 c≃ −0.9797 + 0.2006i, ρ2 c≃ −0.9797 −0.2006i. (6.133) As with the case of ¯ω= 0.5, the Floquet exponents associated with the Floquet multipliers ρ1 cand ρ2 cof Eq. (6.132) correspond to the eigenvalues λ1 cand λ2 cof Eq. (6.131). For ¯ω= 0.75, the magnitude of the Floquet multipliers ρ1 cand ρ2 c verify |ρ1 c|= 1 and |ρ2 c|= 1, and the periodic solution is stable. Figure 6.20 shows a qualitative representation of the evolution of the Floquet exponents λj c and Floquet multipliers ρj c, with j= 1,2, varying the angular velocity ¯ω. The cross markers correspond to an initial speed ¯ωfor which the system is unstable, and the square markers correspond to a value of ¯ωfor which the system is stable. Note that, once the gyroscopic stabilization is achieved, the Floquet exponents are located on the imaginary axis. This is equivalent to the Floquet multipliers being situated on the unit circle. Figure 6.20: Qualitative representation of the evolution of the Floquet exponents and Floquet multipliers with ¯ω. Therefore, the use of Floquet theory to assess the linear stability of the periodic solution provides the same results as those obtained with the time-independent Jacobian matrix of Eq. (6.76), validating the linear stability analysis. The third numerical linearization approach, developed in Sect. 4.2.3, showed an excellent power to obtain the linearized equations of motion of the rolling toroidal wheel along the circular steady motion and to perform the corresponding linear stability analysis. The results have been analytically obtained in terms of the geometric parameters of the toroidal wheel. Furthermore, these results confirm that the procedure is suitable for performing the linearization along reference solutions that lead to linear time-varying (LTV) systems, as the circular steady motion of the rolling wheel shown in this chapter. 256 Chapter 7 Linearization of the equations of motion of multibody systems in multiphysics scenarios The linearization approaches proposed in Chapter 4 can be applied in multiphysics scenarios and, in particular, with hydraulically actuated multibody systems, which have extensive applications in heavy machinery and are utilized in a diverse range of industries and research areas. Hydraulically driven multibody systems play a vital role in several applications such as aircraft (activation and motion of landing gears or flaps), hydraulic lifts, mobile machines (cranes or excavators), hydraulic power steering, vehicle braking systems, hydraulic jacks, and shock absorbers. Modelling and simulating hydraulic systems allow for the evaluation of their behavior across a wide range of operating conditions. Simulation enables addressing crucial concerns like positioning accuracy of the systems or identification of peak stresses. The use of the lumped fluid theory [272] allows for obtaining computationally efficient models of hydraulic circuits. This theory can be applied to hydraulic systems in which the effect of acoustic waves is negligible, which is generally the case of mobile working machine applications, where the motions are relatively slow and pipelines are usually short. The lumped fluid method is widely used in several works with hydraulically actuated multibody systems [273, 274, 275, 276] and modelling complex hydraulic circuits for real-time simulation [277]. Different approaches exist to include the behavior of hydraulic actuators in the simulation of multibody systems. One simplified technique is the kinematic guidance of the variable corresponding to the endpoints of the hydraulic actuator [278]. In a hydraulically actuated multibody system, there is a coupling between the 257 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios dynamics of the mechanical system and the hydraulic actuator. Including the dynamics of the hydraulic actuator is necessary in many applications [279, 280]. From the standpoint of numerical integration of the multiphysics problem, two main approaches exist in the literature: the monolithic approach and multirate integration. The first one, known as monolithic approach or unified scheme, consists of the combination of the multibody system and the hydraulic equations, obtaining a single system of differential equations that can be integrated in time [281, 282]. Docquier et al. [283] presented a multibody model of a modern car equipped with a novel suspension system using a monolithic scheme. Ylinen et al. [284] proposed a hydraulic cylinder model for dynamic simulations, coupling the hydraulic and mechanical variables in a monolithic way. A monolithic formulation, for combined simulation of multibody and hydraulic systems, based on the index-3 augmented Lagrangian, was presented by Naya et al. [285]. Rahikainen et al. [286] combined the index-3 semi-recursive formulation [287] and the lumped fluid theory for hydraulically driven multibody systems, and an improvement of the proposed monolithic formulation was proposed by introducing the singular perturbation method [273]. Jaiswal et al. [275] compared the performance of the double-step semi-recursive approach and the penalty-based semi-recursive approach with hydraulically driven multibody systems, using a hydraulically actuated four-bar and quick-return mechanisms as case studies. Lastly, within the monolithic framework, the friction modelling in the hydraulic cylinder, which plays an important role in the accuracy of the simulations and can introduce numerical stiffness [288], was thoroughly studied [276]. Four friction modelling approaches were compared in terms of the work cycle, friction force, energy balance and numerical efficiency. As an alternative to the monolithic schemes, in a second approach, known as multirate integration, there exist different subsystems that are integrated separately, exchanging information between them in predetermined time intervals. This can be carried out by using a single environment where the different problems are integrated separately (co-integration) [289, 290] or resorting to a different software for each problem (co-simulation) [291, 280, 292, 293]. Some relevant aspects have been addressed in the literature, such as co-simulation configuration [294, 295], energy-based coupling error minimization [294, 296, 297] or the multirate cosimulation [298]. Chapter 3 showed an overview of the main linearization procedures for constrained multibody systems that exist in the literature. Next, Chapter 4 presented novel linearization approaches for multibody systems with holonomic and nonholonomic constraints. Nevertheless, to the best knowledge of the author, there is no procedure for systematically obtaining the linearized equations of motion of hydraulically actuated multibody systems, which is required in several applications. Firstly, an important application is the development of state and input estimators, such as Kalman filters, which have been extensively studied in the context of hydraulically driven multibody systems. Khadim et al. [299] proposed a parameter estimation algorithm that combines the augmented discrete extended Kalman 258 7.1. Formulation of the problem filter (ADEKF) with a curve-fitting method, and Jaiswal et al. [228] presented a state estimator based on an indirect Kalman filter. Secondly, linear stability analyses and the design of linear feedback controllers are other significant applications. Hydraulic steering systems play a crucial role in maintaining directional stability and tracking the steer-handling capability of articulated steering vehicles (ASVs). Previous works [300, 301, 302] have used simplified models that do not consider the dynamic characteristics, where the hydraulic steering system is represented as a torsion spring. Pazooki et al. [303] modelled the full-hydraulic steering control unit of the steering system as a directional control valve. Numerous studies have focused on improving the yaw stability of ASVs using techniques such as increasing damping at the articulation joint [304, 305], introducing leakage across the cylinders [300, 306], and employing active strategies such as differential braking [307, 308], torque vectoring [301] and active steering [309]. Gao et al. [310] analyzed the stability of an ASV and designed a stability controller to mitigate oscillatory yaw motion using optimal control theory. Nedi´c et al. [311] addressed the modelling and simulation of a check valve in a hydraulic transmission. Moreover, an optimal tuned cascade control strategy, based on feedback linearization, was proposed by Nedi´c et al. [312] to carry out the reference trajectory of a 6-DOF parallel robot platform. Djordjevic et al. [313] tackled the control problem of hydraulic servo actuators (HSAs) using a linear discrete model and an online learning data-driven controller based on measured input and output data. In all these applications, a linearized version of the equations of motion is required. Therefore, the main objective of this chapter is to apply the third numerical linearization approach (Sect. 4.2.3) for linearizing the equations of motion of hydraulically actuated multibody systems with holonomic and nonholonomic constraints. First, the main aspects of hydraulic modelling, using the lumped fluid method, are summarized and the equations of motion of hydraulically actuated multibody systems are presented. Next, the linearization approach is detailed, and its use is illustrated with a hydraulically actuated three-dimensional fourbar mechanism model. The linearization procedure is validated by means of the forward dynamics simulation of the linear and nonlinear systems. 7.1 Formulation of the problem The fundamental modelling aspects of a hydraulic linear actuator system are presented. Moreover, the equations of motion of hydraulically actuated multibody systems with holonomic and nonholonomic constraints are shown below. 7.1.1 Modelling of the hydraulic system In this chapter, the lumped fluid theory is used to model the hydraulic system. The hydraulic circuit is divided into discrete volumes of uniformly distributed pressure. The hydraulic system analyzed in this study comprises a double-acting 259 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios hydraulic cylinder, a throttle valve, a directional control valve, a pump, a tank and the connecting hoses. For the sake of simplicity, any leakage in the hydraulic components is disregarded. Figure 7.1 shows a scheme of the hydraulic circuit with its main components. Figure 7.1: Scheme of the hydraulic circuit considered in this study: main components. Considering a discrete control volume V, the evolution of the pressure pis given by the following first-order differential equation: ˙p=Be V nf X k=1 Qk−dV dt!,(7.1) where Qkrepresents the incoming (positive value) or outcoming (negative value) flows of the control volume, nfis the number of hydraulic flows going in or out of the volume, dV dtis the volume change term that usually represents the piston movement inside the cylinder and Beis the effective bulk modulus of the hydraulic volume, given by: 1 Be =1 Boil + ns X k=1 1 Bk Vk V.(7.2) where V=Pns k=1 Vk. In Eq. (7.2), Boil is the oil bulk modulus, nsis the number of subvolumes Vkforming the volume Vand Bkis the bulk modulus of the subvolume Vk. 260 7.1. Formulation of the problem The semi-empirical modelling method [314] is used to describe the valves in the hydraulic circuit. First, the volume flow rate through a directional control valve Qdis: Qd=CvdUp|∆p|∆p |∆p|,(7.3) where Cvdis the semi-empirical flow rate coefficient of the directional control valve; ∆pis the pressure difference, which presents the same direction as the volume flow rate; and Uis a normalized spool displacement, which determines the spool position in the control valve. Equation (7.3) corresponds to a turbulent flow. The time evolution of Uis given by the following first-order differential equation: ˙ U=Uref −U τv ,(7.4) where Uref is the reference normalized spool displacement and τvis the time constant of the valve. In the case of a throttle valve, the turbulent flow rate Qtis expressed as: Qt=Cvtp|∆p|∆p |∆p|,with Cvt=CdAtr2 ρoil .(7.5) In Eq. (7.5), Atis the area of the throttle valve, Cdis the coefficient of discharge and ρoil is the oil density. The expressions of the volume flow rates in Eqs. (7.3) and (7.5) correspond to a turbulent flow, with a Reynolds number of Re > 2400. The derivatives of these expressions with respect to the pressure drop may lead to numerical problems for small pressure differences, since the square root functions of Eqs. (7.3) and (7.5) present a vertical tangent for ∆p→0. For low Reynolds number, the laminar regime is considered, and the laminar-turbulent transition depends on the value of the Reynolds number. In this study, the laminar regime of the volume flow rates is considered for pressure drops lower than a predefined ∆plim, following a linear relation with the pressure difference. It is assumed that ∆plim = 2 bar. Therefore, the volume flow rate through a directional control valve or a throttle valve can be written as: Q=   γ1∆p|∆p| ≤ ∆plim γ2p|∆p|∆p |∆p||∆p|>∆plim (7.6) where γ1and γ2are functions of the valve parameters. To ensure the continuity of the laminar and turbulent regimes at ∆plim, corresponding to a volume flow rate Qlim (see Fig. 7.2 (a)), the relation γ1=√∆plim ∆plim γ2is verified. From Eq. (7.3), γ2=CvdUand γ1=CvdU ∆plim √∆plim in the case of a directional control valve, while γ2=Cvtand γ1=Cvt ∆plim √∆plim for throttle valves. 261 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios (a) Volume flow rate Qas a function of drop pressure ∆p: laminar and turbulent regimes (b) Free-body diagram of the hydraulic cylinder Figure 7.2: Laminar and turbulent regimes of the volume flow rate and free-body diagram of the hydraulic cylinder. The hydraulic actuation is performed by means of a hydraulic cylinder, which transforms hydraulic pressure into mechanical force. The force in the extension direction of the cylinder, denoted by Fcyl, is computed as follows: Fcyl =A1p1−A2p2−Fµ,(7.7) where A1and A2are the piston and piston-rod side areas of the cylinder, respectively; p1is the pressure of the piston side chamber; p2is the pressure of the piston-rod side chamber and Fµis the friction force, arising from the contact between the seal material with the cylinder wall and cylinder rod. Figure 7.2 (b) shows a free-body diagram of the hydraulic cylinder. Following Ref. [276], the Brown and McPhee model [315] is numerically one of the most efficient approaches to describe the friction force in hydraulically driven multibody systems. This friction model incorporates the Coulomb, Stribeck and viscous friction, and is valid for both positive and negative relative tangential velocity. Therefore, the force Fµcan be written as: Fµ=Fctanh 4˙s vs+ (Fs−Fc) ˙s vs 1 4˙s vs2 +3 4!2+σ2˙s, (7.8) where Fcand Fsare the Coulomb and static friction forces, respectively, vsis the transition velocity and σ2is the coefficient of viscous friction. 7.1.2 Equations of motion of hydraulically actuated constrained multibody systems A multibody system with ngeneralized coordinates, mholonomic constraints, l nonholonomic constraints, a hydraulic system with rdiscrete volumes and unor- 262 7.2. Linearization approach for hydraulically actuated multibody systems malized spool displacement variables, which determine the spool positions of the directional control valves, is considered. The equations of motion of a hydraulically actuated constrained multibody system constitute an index-3 Differential- Algebraic system of Equations (DAE), given by the dynamic equilibrium equations, the holonomic and nonholonomic constraints, the hydraulic equations and the valves dynamics equations: M(x)¨x+DT(x)Λ=Q(x,˙x,p),(7.9) C(x) = 0,(7.10) Cnh (x,˙x) = B(x)˙x=0,(7.11) ˙p=H(x,˙x,p,U),(7.12) ˙ U=W(U,Uref ).(7.13) The notation employed in the previous chapters is used. In Eqs. (7.9)-(7.13), x is the n×1 vector of generalized coordinates, which belongs to a domain of Rn; pis the r×1 vector of pressures of the hydraulic system, defined in Rr;Uis the u×1 vector of normalized spool displacement variables existing in the hydraulic system; Uref is the u×1 vector of reference normalized spool displacements; Λis the (m+l)×1 vector of Lagrange multipliers; M(x) is the n×nmass matrix; Q(x,˙x,p) is the n×1 vector of generalized forces; H(x,˙x,p,U) corresponds to the right-hand side of the pressure equations, which are built based on Eq. (7.1); W(U,Uref ) is the right-hand side of the valves dynamics equations, computed as defined in Eq. (7.4); C(x) is the m×1 vector of holonomic constraints; and Cnh (x,˙x) is the l×1 vector of nonholonomic constraints, linearly dependent on velocities. The matrices B(x) and D(x) are l×nand (m+l)×n, respectively, given by: B(x) = ∂Cnh (x,˙x) ∂˙x,D(x) = Cx(x) B(x),(7.14) where Cx=∂C ∂x. The time derivative of the holonomic constraints can be assembled with the nonholonomic constraints, resulting in the following nonlinear index-2 DAE system: M(x)¨x+DT(x)Λ=Q(x,˙x,p), D(x)˙x=0, ˙p=H(x,˙x,p,U), ˙ U=W(U,Uref ). (7.15) 7.2 Linearization approach for hydraulically actuated multibody systems Consider a reference solution of the system of equations (7.9)-(7.13), given by x0(t), ˙x0(t), ¨x0(t), Λ0(t), p0(t), ˙p0(t), U0(t) and ˙ U0(t). Therefore, the following 263 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios relations are verified: Mx0¨x0+DTx0Λ0=Qx0,˙x0,p0,(7.16) Cx0=0,(7.17) Cnh x0,˙x0=Bx0˙x0=0,(7.18) ˙p0=Hx0,˙x0,p0,U0,(7.19) ˙ U0=WU0,Uref ,(7.20) where the time dependence has been omitted for simplicity. First, the variations ˜x,˙ ˜x,¨ ˜x,˜ Λ,˜p,˙ ˜p,˜ Uand ˙ ˜ Uwith respect to this reference solution are defined: ˜x=x−x0,˜ Λ=Λ−Λ0, ˙ ˜x=˙x−˙x0,˜p=p−p0,(7.21) ¨ ˜x=¨x−¨x0,˙ ˜p=˙p−˙p0, ˜ U=U−U0,˙ ˜ U=˙ U−˙ U0. Performing the Taylor expansion of the dynamic equilibrium equations:  Mx0+ n X j=1 ∂(M(x)) ∂xj0 ˜xj+... ¨ ˜x+¨x0 + DTx0+ n X j=1 ∂DT(x) ∂xj0 ˜xj+... ˜ Λ+Λ0= (7.22) Qx0,˙x0,p0+∂Q ∂x0 ˜x+∂Q ∂˙x0 ˙ ˜x+∂Q ∂p0 ˜p... Equation (7.22) can be linearized by retaining up to first-order terms. Using Eq. (7.16), yields: Mx0¨ ˜x+∂M(x)¨x0 ∂x0 ˜x+DTx0˜ Λ+∂DT(x)Λ0 ∂x0 ˜x= ∂Q ∂x0 ˜x+∂Q ∂˙x0 ˙ ˜x+∂Q ∂p0 ˜p, (7.23) where the partial derivatives are evaluated for the reference solution. Similarly, the linearization of the velocity constraints in Eq. (7.15) leads to: Dx0˙ ˜x+∂D(x)˙x0 ∂x0 ˜x=0,(7.24) where it has been used that Dx0˙x0=0. 264 7.2. Linearization approach for hydraulically actuated multibody systems Furthermore, the Taylor expansion of the hydraulic equations (7.12) with respect to the reference solution yields: ˙ ˜p+˙p0=Hx0,˙x0,p0,U0+∂H ∂x0 ˜x+∂H ∂˙x0 ˙ ˜x+∂H ∂p0 ˜p+∂H ∂U0 ˜ U... (7.25) Simplifying by using Eq. (7.19) and retaining up to first-order terms in Eq. (7.25): ˙ ˜p=∂H ∂x0 ˜x+∂H ∂˙x0 ˙ ˜x+∂H ∂p0 ˜p+∂H ∂U0 ˜ U.(7.26) Lastly, the Taylor expansion of the first-order differential equations (7.13) describing the valves dynamics results in: ˙ ˜ U+˙ U0=WU0,Uref +∂W ∂U0 ˜ U... (7.27) Using Eq. (7.20) and that the second- and higher-order derivatives in Eq. (7.27) are null, yields: ˙ ˜ U=∂W ∂U0 ˜ U.(7.28) To reduce the number of linearized equations of motion, the generalized coordinate partition of Eq. (4.96) is used. Given that the index-3 DAE system of equations (7.9)-(7.13) presents mnonlinear holonomic constraints, the n-coordinates vector is split into mdependent coordinates xdand n−madmissible position coordinates xa:x=xaxdT. Moreover, the lnonholonomic constraints allow distinguishing between ldependent admissible velocities ˙xad and n−m−l independent admissible velocities ˙xai in the set of admissible velocities ˙xa, and therefore ˙xa=˙xai ˙xad T. The same partition can be considered at position level: x=xai xad xdT,(7.29) and for the vector of variations ˜x: ˜x=˜xai ˜xad ˜xdT.(7.30) The admissible dependent coordinates ˜xad and the dependent coordinates ˜xdcan be grouped in the set ˜xdd =˜xad ˜xdT. The steps of the approach to obtain the linearized equations of motion are listed below, with the main result of each step shown in a box. Although these steps are the same as those listed in Sect. (4.2.3), they are repeated here for convenience since the equations include additional terms associated with the hydraulic system. Step 1. Eliminate the variations of the Lagrange multipliers and reduce the number of linearized dynamic equations. By premultiplying Eq. (7.23) by TTx0(the transformation matrix Tx0was defined in 265 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios The coordinates xG2,yG2and zG2determine the position of G2, while the orientation of the input link is described by the angles ψ2,θ2and ϕ2. Among these angles, ϕ2represents the rotation of the input link in the Y1Z1-plane. The angle ϕ32 denotes the relative orientation between the coupler and the input link, while ϕ43 corresponds to the relative rotation between the output link and the coupler. It is important to note that the system’s constraints ensure that the motion of the mechanism takes place in the Y1Z1-plane. Additionally, the length of the hydraulic cylinder, which is the distance between points Ocand G2, is denoted by the coordinate s. Lastly, αis the angle between the hydraulic cylinder and the horizontal axis Y1. Figure 7.3 provides a three-dimensional representation of the four-bar mechanism model, displaying the body numbering, all the generalized coordinates and the body reference frames. The orientation matrices of the body frames, expressed as functions of elemental rotation matrices, are defined as follows: R2(x) = Rψ2Rθ2Rϕ2,R3(x) = R2Rϕ32 , R4(x) = R3Rϕ43 ,Rc(x) = Rψ2Rθ2Rα,(7.73) where Rc(x) determines the orientation of the hydraulic cylinder. Figure 7.3: Three-dimensional four-bar mechanism model with hydraulic actuation: numbering of the bodies, coordinates and body frames. 272 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model The lengths and masses of the links are represented as Ljand mj, respectively, with j= 1 . . . 4. The distance between the revolute joints Ocand O1is denoted as d1. The moments of inertia of the bodies, assuming they are slender rods, are defined relative to their centers of mass and are given by ¯ Ijxx =¯ Ijzz =1 12mjL2 j and ¯ Ijyy = 0. The products of inertia are zero due to the symmetries of the bodies. Table 7.1 summarizes all the mechanical parameters (geometric and inertial) of the multibody system, including their numerical values, being shown in Fig. 7.4. The mechanical parameters shown in Table 7.1 are grouped in the set Pm. Table 7.1: List of the mechanical parameters of the four-bar mechanism model, with the corresponding numerical values. Parameter Symbol Numerical value Length of link 1 L110 m Length of link 2 L22 m Length of link 3 L38 m Length of link 4 L45 m Distance OcO1d11 m Mass of link 2 m22 kg Mass of link 3 m38 kg Mass of link 4 m45 kg Figure 7.4: View of the four-bar mechanism in the Y1Z1-plane: parameters Lj and d1. With n= 10 coordinates being used, the multibody system requires m= 9 holonomic constraints, since the system’s degrees of freedom is ng=n−m= 1. 273 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios The holonomic constraints C(x) can be defined as follows: C(x) =           rO1 u1·v u2·v rY O2−L1 rZ O2 d1+rY OcG2−rY O1G2 rZ OcG2−rZ O1G2           =09×1.(7.74) In Eq. (7.74), the first five constraints correspond to the revolute joint between bodies 1 and 2, located at point O1. This joint constrains the rotation of body 2 to the Y1Z1-plane. The absolute position vector of this revolute joint is represented by rO1. The unit vector valong the X2-axis is defined in the body reference frame 1 and is constrained to be parallel to the X1-axis due to the fourth and fifth holonomic constraints. Additionally, u1and u2are unit vectors perpendicular to the X1-axis. These vectors are computed as follows: rO1=rG2+R2¯ rG2O1, rG2=xG2yG2zG2T, ¯rG2O1=0−L2 20T ,(7.75) u1=010T, u2=001T, v=cos (ψ2) cos (θ2) sin (ψ2) cos (θ2)−sin (θ2)T. Figure 7.3 depicts vectors v,u1and u2. In the next two holonomic constraints of Eq. (7.74), rY O2and rZ O2are the Y- and Z-components of the vector rO2, which is computed as follows: rO2=r2+r3+r4,(7.76) where r2=rG2+R2¯rG2C,¯rG2C=0L2 20T , r3=R3¯rCD,¯rCD =0L30T,(7.77) r4=R4¯rDO4,¯rDO4=0L40T. Finally, the last two holonomic constraints arise from the Y- and Z-components of the loop formed by the hydraulic cylinder, the ground link and the input link, which involves the vectors rOcG2,rO1Ocand rO1G2: rOcG2=Rc¯rOcG2,¯rOcG2=0s0T, 274 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model rO1G2=R2¯rO1G2,¯rO1G2=0L2 20T ,(7.78) rO1Oc=0−d10T. 7.3.2 Description of the hydraulic circuit As stated in Sect. 7.1.1, the hydraulic system analyzed in this study comprises a double-acting hydraulic cylinder, a throttle valve, a directional control valve, a pump, a tank and the connecting hoses. The hydraulic system is divided into three distinct control volumes, labeled as V1,V2and V3. These control volumes are clearly illustrated in Fig. 7.5 and are given by: V1=Vh1, V2=Vh2+A2l2,(7.79) V3=Vh3+A3l3. In Eq. (7.79), Vh1,Vh2and Vh3are the volumes of the hoses of the corresponding control volumes; A2and A3are the surfaces of the piston and piston-rod sides sections, respectively; and l2,l3are the lengths of the piston and piston-rod sides. These lengths are a function of the coordinate s, and can be computed as: l2(s) = s−c1−c2, l3(s) = l−l2(s),(7.80) where c1is the length between the revolute joint Ocand the base of the cylinder chamber, c2is the length of the piston and lis the length of the cylinder chamber. The dimensions c1,c2and l, and the lengths l2,l3are shown in detail in Fig. 7.5. The pressures of the control volumes p1,p2and p3can be computed, following Eq. (7.1), as follows: ˙p1=Be1 V1 (QA1−QA2), ˙p2=Be2 V2 (QA2−A2˙s),(7.81) ˙p3=Be3 V3 (A3˙s−QB1), where Be1,Be2and Be3are the effective bulk moduli of the corresponding control volumes. These moduli are computed following Eq. (7.2): 1 Be1 =1 Boil +Vh1 V1Bh , 1 Be2 =1 Boil +Vh2 V2Bh +A2l2 V2Bc ,(7.82) 1 Be3 =1 Boil +Vh3 V3Bh +A3l3 V3Bc , 275 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios Figure 7.5: Scheme of the hydraulic circuit and dimensions of the hydraulic cylinder. 276 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model with Bhand Bcbeing the bulk moduli of the hoses and the hydraulic cylinder, respectively. The volume flow rates QA1,QA2and QB1of Eqs. (7.81) are defined as shown in Fig. 7.5. The parameters of the hydraulic circuit are listed in Table 7.2, with their corresponding numerical values being those of Ref. [276]. The hydraulic parameters shown in Table 7.2 are grouped in the set Ph. Equations (7.81) of the hydraulic circuit can be written in vector form as in Eq. (7.12): ˙p=H(x,˙x,p, U). It should be noted that, in the hydraulic system of the current example, there is only one directional control valve, requiring only a single normalized spool displacement U. As a result, the equations of motion of the hydraulically actuated four-bar mechanism can be obtained by considering the dynamic equilibrium equations, computed in accordance with Chapter 2, the set of holonomic constraints (7.74), the hydraulic equations (7.81) and Eq. (7.13), which describes the dynamics of the directional control valve. The following index-3 DAE system is obtained: M(x)¨x+CT x(x)Λ=Q(x,˙x,p), C(x) = 0, ˙p=H(x,˙x,p, U), ˙ U=Uref −U τ. (7.83) Note that, due to the absence of nonholonomic constraints, the matrix D(x) defined in Eq. (7.14) becomes Cx(x). 7.3.3 Computation of the Jacobian matrix Prior to the linearization, the equilibrium of the four-bar mechanism multibody model is defined. In the equilibrium configuration, which corresponds to the neutral position of the directional control valve (U= 0), the coordinates are given by: x0=x0 G2y0 G2z0 G2ψ0 2θ0 2ϕ0 2ϕ0 32 ϕ0 43 s0α0T,(7.84) where x0 G2= 0, y0 G2=L2 2cos ϕ0 2, z0 G2=L2 2sin ϕ0 2, ψ0 2= 0,(7.85) θ0 2= 0, s0=sd2 1+L2 22 +d1L2cos (ϕ0 2), tan α0=sin ϕ0 2 cos (ϕ0 2) + 2d1 L2 . 277 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios Table 7.2: List of the hydraulic parameters and numerical values. Hydraulic parameter [unit] Symbol Numerical value Pressure of the pump [MPa] pP7.6 Pressure of the tank (atmospheric pressure) [MPa] pT0.1 Density of the oil [kg/m3]ρoil 850 Bulk modulus of the oil [MPa] Boil 1500 Bulk modulus of the hoses [MPa] Bh550 Bulk modulus of the hydraulic cylinder [MPa] Bc31500 Volume of the hose (control volume 1) [m3]Vh14.71 ·10−5 Volume of the hose (control volume 2) [m3]Vh23.14 ·10−5 Volume of the hose (control volume 3) [m3]Vh37.85 ·10−5 Area of the throttle valve [m2]At2.83 ·10−5 Semi-empirical flow rate constant of the directional control valve hm3/s√Pai Cvd2.138 ·10−8 Flow discharge coefficient of the throttle valve Cd0.8 Diameter of the piston [mm] D280 Diameter of the piston-rod [mm] D335 Length of the cylinder chamber [m] l0.9 Constant 1 [m] c10.43 Constant 2 (length of the piston) [m] c21 Coulomb friction [N] Fc210 Static friction [N] Fs830 Stribeck velocity [m/s] vs1.25 ·10−2 Coefficient of viscous friction [Ns/m] σ2330 Time constant of directional control valve [s] τ0.0045 278 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model Given an equilibrium angle ϕ0 2for the input link, the coordinates in the equilibrium configuration can be determined by solving the holonomic constraints, which are particularized for the equilibrium position: Cx0=0. The angles ϕ0 32 and ϕ0 43 are obtained numerically by solving nonlinear equations, while the remaining coordinates can be analytically computed using Eqs. (7.85). Furthermore, by using the dynamic equilibrium equations and the hydraulic equations particularized for the equilibrium position: Mx0¨x0+DTx0Λ0=Qx0,˙x0,p0,(7.86) ˙p0=Hx0,˙x0,p0, U0,(7.87) where U0= 0 and ˙x0=¨x0=˙p0=0in the equilibrium, the pressures p0are obtained: p0=p0 1p0 2p0 3T,(7.88) with p0 1=p0 2, p0 2=g(p0 3,x0,Pm,Ph).(7.89) Note that, in Eq. (7.89), g(p0 3,x0,Pm,Ph) allows obtaining the equilibrium pressure p0 2as a function of p0 3, and is given by: g(p0 3,x0,Pm,Ph) = A3 A2 p0 3−g(m2+ 2m3+m4) sin ϕ0 2+ϕ0 43 2A2sin (α0−ϕ0 2) sin (ϕ0 43) −g(m2+m3) sin ϕ0 2−ϕ0 43+ (m3+m4) sin ϕ0 2+ 2ϕ0 32 +ϕ0 43 2A2sin (α0−ϕ0 2) sin (ϕ0 43).(7.90) Lastly, the Lagrange multipliers in the equilibrium Λ0are also derived from Eqs. (7.86), leading to: Λ0=0 Λ0 2Λ0 300Λ0 6Λ0 70 0 T,(7.91) with the nonzero Lagrange multipliers being: Λ0 2=g 8 sin (α0−ϕ0 2) sin (ϕ0 43)µ1x0,Pm, Λ0 3=g 8 sin (α0−ϕ0 2) sin (ϕ0 43)µ2x0,Pm,(7.92) Λ0 6=−g(m3+m4) 4 sin (ϕ0 43)cos 2ϕ0 2+ 2ϕ0 32 +ϕ0 43+ cos ϕ0 43, Λ0 7=−g 4 sin (ϕ0 43)(m3+m4) sin 2ϕ0 2+ 2ϕ0 32 +ϕ0 43+ (m3+ 3m4) sin ϕ0 43. In Eq. (7.92), µ1x0,Pmand µ2x0,Pmare functions of the coordinates in the equilibrium position x0and the mechanical parameters Pm: µ1x0,Pm=−(m3+m4)sin α0+ϕ0 2+ 2ϕ0 32 +ϕ0 43279 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios + sin 3ϕ0 2−α0+ 2ϕ0 32 +ϕ0 43 + (2m2+ 4m3+ 2m4) sin α0+ϕ0 2+ϕ0 43 + (2m2+ 3m3+m4) sin ϕ0 2−α0+ϕ0 43 + (2m2+ 3m3+m4) sin α0−ϕ0 2+ϕ0 43 −2 (m2+m3) sin ϕ0 2+α0−ϕ0 43 −2 (m3+m4) sin ϕ0 2−α0+ 2ϕ0 32 +ϕ0 43,(7.93) µ2x0,Pm= (m3+m4)cos α0+ϕ0 2+ 2ϕ0 32 +ϕ0 43 + cos 3ϕ0 2−α0+ 2ϕ0 32 +ϕ0 43 −(2m2+ 4m3+ 2m4) cos α0+ϕ0 2+ϕ0 43 + (m3+m4−2m2) cos ϕ0 2−α0+ϕ0 43 + (2m2+m3+m4) cos α0−ϕ0 2+ϕ0 43 + 2 (m2+m3) cos ϕ0 2+α0−ϕ0 43 −2 (m3+m4) cos ϕ0 2−α0+ 2ϕ0 32 +ϕ0 43.(7.94) Given that the system presents one degree of freedom, the input link angle ϕ2is chosen as independent coordinate. The sets of independent and dependent coordinates ˜xai and ˜xd, defined in Eq. (4.96), are therefore given by: ˜xai =˜ ϕ2,˜xd=˜xG2˜yG2˜zG2˜ ψ2˜ θ2˜ ϕ32 ˜ ϕ43 ˜s˜αT.(7.95) The use of Eq. (7.71), particularized for the equilibrium of the four-bar mechanism model, leads to the following Jacobian matrix: J=           0 1 0 0 0 0 ν1ν20ν3ν40 0 0 ν5−ν50 0 0ν6ν7−ν70 0 0ν80 0 0 0 0 0 0 0 0 −1 τ           ,(7.96) where the coefficients νkin Eq. (7.96), with k= 1 . . . 8, are functions of the coordinates and pressures in the equilibrium, x0and p0, respectively, and the geometric and hydraulic parameters, Pmand Ph:νk=νkx0,p0,Pm,Ph. In the computation of the linearized hydraulic equations, the laminar regime of the volume flow rates through the directional control valve and the throttle valve described in Eq. (7.6) has been considered. Note that, in this case, n= 10, m= 9, l= 0, r= 3 and u= 1, and therefore the Jacobian matrix (7.96) is (2n−2m−l+r+u)×(2n−2m−l+r+u)=6×6. 280 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model 7.3.4 Numerical results: validation of the approach To verify the accuracy of the linearization results, forward dynamic simulations are conducted for both the nonlinear four-bar mechanism multibody model and the linearized system. For the numerical integration of the nonlinear system, the holonomic constraints in Eq. (7.83) are differentiated twice with respect to time. This yields the following index-1 DAE, which is solved using the MATLAB ODE15s solver: M(x)¨x+CT x(x)Λ=Q(x,˙x,p) + QM(x), Cx(x)¨x=−˙ Cx(x)˙x−2βe˙ C(x)−γ2 eC(x), ˙p=H(x,˙x,p, U), ˙ U=Uref −U τv . (7.97) Note that the differentiation of the holonomic constraints may result in the violation of the constraints. Therefore, as described in Sect. 2.3.5, Baumgarte stabilization method [261] is used in Eq. (7.97) to avoid numerical drift, βeand γebeing the Baumgarte stabilization constants. To minimize the numerical drift of the constraints while maintaining a reasonable integration cost, a numerical value of βe=γe= 50 has been chosen. It has been verified that the results of the forward dynamics simulation are not affected by the specific numerical values of these stabilization constants. Since the linearization has been performed around the equilibrium state, which corresponds to the neutral position of the directional control valve (U= 0), a reference normalized spool displacement Uref = 0 is used in the simulation. The equilibrium position x0in Eq. (7.84) is considered as the initial condition at position level: x(0) = x0, with the following numerical values: xG2(0) = 0 m, yG2(0) = 0.5 m, zG2(0) ≃0.866 m, ψ2(0) = 0, θ2(0) = 0, ϕ2(0) = 60◦,(7.98) ϕ32 (0) ≃322.1◦, ϕ43 (0) ≃266.4◦, s(0) = √3 m, α (0) = 30◦. The exact numerical values in Eq. (7.98), corresponding to ϕ2(0) = 60◦, can be computed by solving the holonomic constraints for this equilibrium angle, as detailed after Eq. (7.85). Moreover, ˙x(0) = 010×1and the equilibrium pressures, p0, given by Eq. (7.89), are considered as initial pressures: p(0) = p0, with p0 3= 3.5 MPa: p1(0) = p2(0) ≃2.82 MPa, p3(0) = 3.5 MPa.(7.99) Note that the exact numerical value of p2(0) is computed by using Eq. (7.89). The term QM(x) in Eq. (7.97) is the generalized force vector associated with the 281 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios with respect to those of the second scenario. In the particular case of no friction force (Fc=Fs=σ2= 0), the time evolutions of ϕ2(t) and the pressures are qualitatively similar to those of this scenario (Fc=1 100F∗ c, Fs=1 100F∗ sand σ2= 1 100σ∗ 2), given that there exists a second source of dissipation due to the throttle valve, and the eigenvalues λ3and λ4present negative real parts. The numerical values of Fc=1 100F∗ c, Fs=1 100F∗ sand σ2=1 100σ∗ 2have been chosen to improve the visualization of the oscillations of the system. As in the overdamped case, the linear responses of the underdamped scenarios accurately reproduce the nonlinear results for small amplitude oscillations. Morover, a decrease in the numerical value of the maximum external torque, M0, in Eq. (7.100) further improves the results of the linear model, since the multibody system is closer to its equilibrium position. In contrast, when the external torque is too high, the multibody system deviates significantly from the equilibrium position, and therefore the linear model ceases to be valid. Despite the linearized system accurately reproduces the nonlinear system in all the cases, the results in the overdamped scenario are better than in the underdamped cases. This can be explained with the linearization of the friction force. Figure 7.13 shows the friction force Fµpresented in Eq. (7.8), with Fc= 210 N, Fs= 830 N and σ2= 330 Ns/m, together with its first-order Taylor approximation, denoted as T1 µ. It can be seen that the linear approximation is excellent for low values of the velocity ˙s, as is the case of the overdamped scenario. Nevertheless, the approximation becomes worse for larger values of ˙s, which is the case of the underdamped scenarios. 0246810 _s(mm/s) 0 500 1000 1500 2000 F7(N) F7 T1 7 Figure 7.13: First-order Taylor approximation of the friction force Fµ. 288 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model 7.3.5 Comparison of the proposed linearization procedure with a conventional counterpart approach The proposed linearization approach is compared with a conventional counterpart. This counterpart procedure is based on the linearization of an ODE system and the numerical computation of the required partial derivatives by using finite differences. The main steps of this counterpart procedure are shown below. Consider the index-2 DAE system (7.15). The velocity contraints D(x)˙x=0 can be differentiated once with respect to time, leading to the following index-1 DAE system: M(x)¨x+DT(x)Λ=Q(x,˙x,p), D(x)¨x=Qd(x,˙x), ˙p=H(x,˙x,p,U), ˙ U=W(U,Uref ), (7.103) where Qd(x,˙x) was defined in Eq. (2.127). Defining ˙x=vand X=x v p U T, Eqs. (7.103) can be rewritten as the following DAE system: Ah(x)˙ X Λ=fh(X,Uref ),(7.104) with Ah(x) =          In0n0n×r0n×u0n×(m+l) 0nM(x)0n×r0n×uDT(x) 0(m+l)×nD(x)0(m+l)×r0(m+l)×u0(m+l) 0r×n0r×nIr0r×u0r×(m+l) 0u×n0u×n0u×rIu0u×(m+l)          ,(7.105) fh(X,Uref ) =       v Q(x,v,p) Qd(x,v) H(x,v,p,U) W(U,Uref )       .(7.106) Alternatively, Eq. (7.104) can be expressed as follows: ˙ X Λ=A−1 h(x)fh(X,Uref ) = F(X,Uref ).(7.107) Denoting the first 2n+r+ucomponents of F(X,Uref ) as F(X,Uref ), the following ODE system is obtained: ˙ X=F(X,Uref ).(7.108) 289 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios Note that the last m+lequations of Eq. (7.107) provide the expressions of the Lagrange multipliers Λ. The linearization of the ODE system of equations (7.108) with respect to an arbitrary reference solution X0=x0˙x0p0U0Tleads to the following Jacobian matrix: J=∂F ∂X0 =J1··· Jj··· JN.(7.109) In Eq. (7.109), N= 2n+r+uand Jjrepresents the jth column of the Jacobian matrix, with j= 1 . . . N. The term Jjis computed as follows: Jj=∂F ∂Xj0≃FX0+hej−FX0−hej 2h,(7.110) where the partial derivative in Eq. (7.110) can be numerically computed by using central finite differences. In Eq. (7.110), his the step used in the computation of the finite differences, and ejcorresponds to the jth column of the N×Nidentity matrix IN. The Jacobian matrix of Eq. (7.109) is (2n+r+u)×(2n+r+u). A thorough comparison is conducted between the linearization approach proposed in Sect.7.2 and its counterpart procedure using the hydraulically actuated four-bar mechanism model described in Sects. 7.3.1 and 7.3.2. First, a qualitative comparative analysis is performed. The preliminary steps required to obtain the Jacobian matrices of Eqs. (7.62) and (7.109) differ significantly. In the counterpart approach, the index-1 DAE of Eq. (7.103) is first transformed into Eq. (7.107). Then, the equations associated with the Lagrange multipliers are removed to obtain the first-order ODE system of equations (7.108). Finally, the linearization is carried out, resulting in the Jacobian matrix of Eq. (7.109). On the other hand, the proposed linearization approach follows the steps outlined in Sect.7.2. The DAE system of Eq. (7.15) is initially linearized around the reference solution. Subsequently, a coordinate partitioning in terms of independent and dependent coordinates is used, reducing the number of linearized dynamic equations and eliminating the variations of the Lagrange multipliers (as shown in Eq. (7.31)). Then, the linearized constraints at velocity and acceleration levels are utilized to express the complete set of velocity and acceleration variations ˙ ˜xand ¨ ˜xin terms of the independent variations ˙ ˜xai and ¨ ˜xai. Finally, the linearized holonomic constraints are employed to express the dependent coordinate variations ˜xdin terms of the independent ones ˜xai. All these steps lead to a reduced system of linear equations, which in the case of a multibody system with only holonomic constraints as the four-bar mechanism, is only expressed in terms of the variations of the independent coordinates ˜xai and their time derivatives. A significant advantage of the proposed approach is that it allows obtaining the maximum possible reduction in the number of linearized equations of motion. The Jacobian matrix of Eq. (7.62) is (2n−2m−l+r+u)×(2n−2m−l+r+u) in the general case of a hydraulically driven multibody system with holonomic and nonholonomic constraints, while the conventional counterpart results in a bulky 290 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model Jacobian matrix (see Eq.(7.109)) with dimensions (2n+r+u)×(2n+r+u). As a result, the proposed approach yields 2n−2m−l+r+ueigenvalues, among which 2 (n−m−l)+r+ucorrespond to the real spectrum of the problem, and lspurious eigenvalues are obtained, which are associated with the linearized nonholonomic constraints (7.54). In the case of a holonomic multibody system, the proposed approach leads to 2 (n−m) + r+ueigenvalues, all of which correspond to the real spectrum, without any spurious eigenvalues. On the other hand, the use of the counterpart linearization approach results in 2n+r+ueigenvalues, with 2 (n−m−l)+r+ueigenvalues corresponding to the real spectrum, and 2 (m+l) spurious eigenvalues associated with the m+ldependent coordinates. In the specific case of the four-bar mechanism, the proposed approach yields the 6 ×6 Jacobian matrix of Eq. (7.96), from which the eigenvalues λk, with k= 1 . . . 6, are obtained as shown in Table 7.3. In contrast, the counterpart approach leads to the 24 ×24 Jacobian matrix of Eq. (7.109), resulting in the six eigenvalues λkof Table 7.3 along with eighteen additional spurious null eigenvalues. Hence, the utilization of the proposed approach enables the reduction of the number of linearized equations of motion and eliminates the spurious null eigenvalues associated with the dependent coordinates, while preserving all the stability information. Another key advantage to highlight is the power of the proposed approach to generate the exact linearized equations of motion as a function of the geometric, dynamic and hydraulic parameters of the multibody system under study. In the particular case of the four-bar mechanism, the coefficients νkof the Jacobian matrix of Eq. (7.96), with k= 1 . . . 8, are computed analytically as a function of the coordinates and pressures in the equilibrium, x0and p0, respectively, and the geometric and hydraulic parameters Pmand Ph:νk=νkx0,p0,Pm,Ph. This greatly simplifies the performance of detailed eigenvalues sensitivity analyses, since the eigenvalues of the system are parameterized in terms of the design parameters. In contrast, the analytical computation of the Jacobian matrix of Eq. (7.109) is not possible for complex multibody systems with large number of coordinates, constraints and long kinematic chains, since the symbolic computation of the inverse matrix A−1 h(x) in Eq. (7.107) is required. For this reason, the Jacobian matrix of Eq. (7.109) is numerically computed by using finite differences. A quantitative comparative analysis between the proposed approach and its counterpart procedure is also included. The comparison focuses on two aspects: efficiency and accuracy. Efficiency is assessed by measuring the computation time required for calculating the Jacobian matrices of Eqs. (7.96) and (7.109). In particular, for the four-bar mechanism model, the average computation time for the Jacobian matrix of Eq. (7.96) is found to be 0.005 s. In contrast, the computation time for the Jacobian matrix of Eq. (7.109) is 0.016 s. The assessment was carried out by using a computer HP with Intel(R) Core(TM) i7-6700HQ 2.6 GHz and 12 GB of RAM. Consequently, the proposed approach exhibits higher efficiency, with a computation time approximately three times lower for the four-bar mechanism. In terms of accuracy, the conventional counterpart approach provides an ap- 291 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios proximation of the Jacobian matrix. Its accuracy is strongly influenced by the step size hused in the computation of the finite differences in Eq. (7.109). In contrast, the proposed approach leads to the exact linearized equations of motion. To compare the accuracy of both approaches, the relative errors ϵ3and ϵ4 in the computation of the eigenvalues λ3and λ4of the overdamped scenario (see Table 7.3) are computed. Figures 7.14 (a) and (b) show the variations of the relative errors ϵ3and ϵ4, respectively, with the step h. To ease the visualization, a logarithmic scale is used for the X-axis. The relative errors ϵ3and ϵ4exhibit significant increases when the step his very small or large, as shown in Figs. 7.14 (a) and (b). The counterpart approach achieves the highest accuracy in determining the eigenvalues λ3and λ4when using an intermediate step value of approximately h≃10−5, resulting in ϵ3≃0 and ϵ4≃0. Hence, selecting an appropriate step size is crucial for maintaining accuracy in the alternative procedure. However, this issue does not arise in the proposed approach since the linearization procedure has the capability to accurately compute the exact Jacobian matrix. -10 -8 -6 -4 -2 Step, log10(h) 0 1 2 3 4 5 Rel. error 63,03(%) (a) Relative error ϵ3in the computation of λ3 -10 -8 -6 -4 -2 Step, log10(h) 0 1 2 3 4 5 Rel. error 64,04(%) (b) Relative error ϵ4in the computation of λ4 Figure 7.14: Relative errors ϵ3and ϵ4in the computation of the eigenvalues λ3 and λ4: influence of the step hused in the central finite differences scheme of Eq. (7.110). Therefore, it is concluded that the main advantages of the proposed linearization procedure are its efficiency and accuracy; the maximum reduction of the number of linearized equations of motion, which allows the elimination of spurious eigenvalues in the performance of linear stability analyses; and the power to generate the exact linearized equations of motion in terms of the mechanical and hydraulic parameters of the multibody system under study. It is important to note that, in this chapter, the linearization has been carried out around the equilibrium position of the four-bar mechanism multibody model. Therefore, the resulting Jacobian matrix presents constant (time-independent) coefficients, which are a function of the multibody system parameters. In contrast, in other applications, as state observers based on the linearized equations of 292 7.3. Validation of the approach with a hydraulically actuated four-bar mechanism model motions like Kalman filters, a linear approximation of the dynamics of the multibody system is required at each instant of time. In this case, the Jacobian matrix must be updated and its coefficients vary throughout the reference solution. 293 Chapter 7. Linearization of the equations of motion of multibody systems in multiphysics scenarios 294 Chapter 8 Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems The linearization approaches developed in Chapter 4 can be used to design linear feedback controllers for nonholonomic multibody systems. The problem of nonholonomic rolling is usual in nonprehensile manipulation tasks, which are relevant in areas such as mobile robotics, robotic manipulation and autonomous systems. Nonprehensile manipulation involves manipulating objects in tasks where grasping them may not be feasible or desirable [316, 317]. Instead, control inputs such as forces or torques are applied to the system to achieve the desired motion. The planning and stabilization of the trajectory of a smooth body rolling on the surface of another [318] is usual within nonprehensile manipulation. Specific examples are the butterfly robot [319, 320, 321] or controlling a sphere rolling on a motion-controlled flat plate [322, 323]. The works of Woodruff and Lynch [318, 324] must be highlighted. In [318], a motion planner and feedback controller for pure-rolling motions between general smooth rigid bodies were presented. The motion planning and control carried out in [318] was limited to first-order kinematics, without dealing with the dynamic equations. The methods were demonstrated for a sphere rolling on a sphere and an ellipsoid rolling on an ellipsoid. Next, their work was extended in [324]. The rolling dynamic equations of smooth rigid bodies were derived, and their results were validated by means of simulations of the rolling sphere on a plane. Two scenarios were simulated: the motion of the sphere on a horizontal rotating plane 295 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems and on a tilted rotating plane. In both cases, the simulation results were in correspondence with the analytical solutions of the problems. Next, the motion planner developed in [318], based on a direct collocation method, was used with the dynamic equations to plan robotic juggling motions. This motion planner was employed to reorientate a sphere on a plate. Lastly, Woodruff and Lynch [318] demonstrated the use of a Linear-Quadratic Regulator (LQR) to stabilize planned trajectories in both simulations (in the previous problem of reorienting the sphere) and experiments (3-DOF robotic arm manipulating an object). In particular, the ball-plate system constitutes a well-studied example of a nonholonomic system in the context of nonprehensile manipulation. In this system, the objective is to control the ball’s motion by manipulating the plate without directly grasping the ball. The nonholonomic constraints arise due to the rolling contact between the ball and the plate. The number of references devoted to the ball-plate system is very large. The work of Knuplez et al. [325] focused on the synthesis of a controller for a two-dimensional electromechanical system consisting of a ball and plate. The study aimed to investigate the system dynamics and perform laboratory experiments with different control approaches based on classical and modern control theory. The control objective was to maintain the freely rolling ball in a specific position on the plate. Ker et al. [326] proposed a controller, based on the backstepping control procedure, whose performance was assessed through tests involving oscillatory stabilization and circular trajectory tracking of the ball-plate system. In the paper by Date et al. [327], the control objective was to move a pair of horizontal plates and a ball held between the plates to a desired position with a desired orientation, maintaining rolling contact without slip. The control was achieved by applying horizontal motion to the plates, allowing both translation and rotation. The problem was divided into two subproblems: position control, which focuses on controlling the positions of the ball and plates while ignoring the ball’s orientation, and orientation control, which aims to control the ball’s orientation while keeping the positional relationship between the ball and upper plate constant. The combined control of the entire configuration is referred to as simultaneous control. Fan et al. [323] proposed a hierarchical fuzzy control scheme for trajectory planning and tracking in the ball-plate system, using genetic algorithms to optimize trajectory planning. Simulation results demonstrated effective control of the ball without hitting obstacles and in a short time. Oriolo et al. [328] presented an iterative control scheme that guarantees asymptotic stability with exponential convergence towards any desired set point, given appropriate conditions. The proposed strategy was illustrated by designing a stabilizing controller for the ball-plate manipulation system. A sufficient stabilizability condition for general nonlinear systems was provided by Casagrande et al. [329] and illustrated by means of the ball-plate system. The Linear-Quadratic Regulator (LQR) algorithm is used in several works devoted to control of the ball-plate system. In the paper by Cheng et al. [330], a skillful robotic wrist system, focused on the application of table tennis, was presented. A ball-plate system was used in the initial stage of their work. For changing the plate’s orientation, an LQR control algorithm was designed, demonstrating that 296 8.1. Extension of the third numerical linearization approach to controlled multibody systems the ball can be balanced at specific positions with different experiments. Kassem et al. [331] addressed the control of a system consisting of a plate mounted on a Stewart platform with 6 degrees of freedom. To tackle the static and dynamic position tracking problem, the paper proposed four strategies: PID, LQR, Sliding Mode and Fuzzy controllers. The authors presented simulation results as well as validation on a real system. A comparison of the proposed strategies was included in their work. Furthermore, other control methods, as Model Predictive Control (MPC), have been applied to the ball-plate system. In the paper by Fan et al. [332], a novel nonlinear Model Predictive Control (NMPC) strategy based on Gaussian particle swarm optimization (GPSO) was presented. The authors utilized the Taylor expansion to transform the NMPC problem into a quadratic programming problem with unknown parameters. Finally, the proposed control strategy was applied to the ball-plate system to verify the effectiveness. Oravec et al. [333] also implemented a control algorithm, based on MPC, verified with nonlinear simulations and reference trajectory tracking in a real ball-plate laboratory model. In the paper by Bang et al. [334], a method for implementing MPC on a microcontroller was proposed, focusing on enhancing tracking performance in the ball-plate system. To improve computation efficiency and facilitate implementation, the MPC quadratic programming problem was converted into an equivalent nonnegative least-squares problem. The MPC proposed by Bang et al. was applied to a tracking control problem using a laboratory-built ball-plate system. In this chapter, the third numerical linearization approach (Sect. 4.2.3) is extended to the case of controlled multibody systems. The expression of the input matrix is derived. Starting from the continuos-time description of the linear system, a discrete-time representation is obtained. Subsequently, a linear feedback controller is developed using a widely extended optimal control method as the Linear-Quadratic Regulator (LQR). As an example of application, the optimal control of the well-known ball-plate system is considered. 8.1 Extension of the third numerical linearization approach to controlled multibody systems The third numerical approach (Sect. 4.2.3) is used to obtain the linearized equations of motion of controlled multibody systems. Consider a multibody system with ngeneralized coordinates, mholonomic constraints, lnonholonomic constraints and nuactuations. Denoting by uthe nu×1 vector of actuations, the equations of motion given by the index-3 DAE (2.122)-(2.124) can be written as: M(x,u)¨x+DT(x)Λ=Q(x,˙x,u),(8.1) C(x) = 0,(8.2) Cnh (x,˙x) = B(x)˙x=0.(8.3) 297 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems Figure 8.1: Multibody model of the ball-plate system. The Cartesian coordinates xCand yCare defined in the local reference frame of body 2 and locate the contact point Cof the sphere on the plate; ψis the yaw angle of the sphere; ϕrepresents the lean angle; and θis the roll angle. The angles ψ, ϕ and θdescribe the orientation of the ball in space. The set of generalized coordinates of the multibody system and the body reference frames are shown in Fig. 8.1. The orientation matrices of the body frames, expressed as a function of the elemental rotation matrices, are given by: R2(x) = RαRβ,R3(x) = R2RψRϕRθ.(8.46) The position vectors of the contact point Cand the centre of mass G3, expressed in the global reference frame, are given by: rC=R2¯rC,¯rC=xCyC0T,(8.47) rG3=rC+R2¯rCG3,¯rCG3=0 0 RT.(8.48) The set of nonholonomic constraints is denoted as Cnh (x,˙x). These nonholonomic constraints result from the assumption that the sphere (body 3) rolls without slipping on the plate surface (body 2): ¯v2 C32 =0,(8.49) where ¯v2 C32 is expressed in the body reference frame 2. The l×1 vector of nonholonomic constraints, with l= 2, is given by: Cnh (x,˙x) = ¯vX C32 ¯vY C32 T=02×1,(8.50) 304 8.4. Example of study: ball-plate system where ¯vX C32 and ¯vY C32 are the X- and Y-components of ¯v2 C32 , which can be obtained as: ¯v2 C32 =RT 2vC32 .(8.51) The relative velocity vC32 , expressed in the global reference frame, can be computed as follows: vC32 =vC31 −vC21 ,(8.52) with vC31 =vG31 +R3(¯ ω3ׯrG3C),(8.53) vC21 =vG21 +R2(¯ ω2ׯrC).(8.54) In Eqs. (8.53) and (8.54), vG21 and vG31 are the velocities of the centres of mass of bodies 2 and 3, respectively; ¯ ω2and ¯ ω3are the angular velocities of bodies 2 and 3, respectively, expressed in the corresponding body frames; and ¯rG3C=−¯rCG3, where ¯rCG3was defined in Eq. (8.48). The expression of velocity ¯v2 C32 , obtained with symbolic computation, is given by: ¯v2 C32 =  ¯vX C32 ¯vY C32 ¯vZ C32  =  ˙xC−R˙ θcos (ψ) cos (ϕ)−R˙ ϕsin (ψ) ˙yC−R˙ θsin (ψ) cos (ϕ) + R˙ ϕcos (ψ) 0 .(8.55) Note that, in Eq. (8.55), the Z-component of ¯v2 C32 is already null and therefore not included in the nonholonomic constraints (8.50). Since the rolling sphere is described by using five coordinates (n= 5), the number of degrees of freedom of the system is three: ng=n−l= 3. As shown in Eq. (8.44), the number of control inputs of the system is nu= 2. Since nu< ng, the system is underactuated. Denoting the nonholonomic constraints as C1 nh and C2 nh: C1 nh = ˙xC−R˙ θcos (ψ) cos (ϕ)−R˙ ϕsin (ψ)=0,(8.56) C2 nh = ˙yC−R˙ θsin (ψ) cos (ϕ) + R˙ ϕcos (ψ)=0,(8.57) the time derivatives of the angles ˙ ϕand ˙ θcan be obtained by combining Eqs. (8.56) and (8.57) as follows: C1 nh sin (ψ)−C2 nh cos (ψ)=0 →˙ ϕ=˙xCsin (ψ)−˙yCcos (ψ) R,(8.58) C1 nh cos (ψ) + C2 nh sin (ψ)=0 →˙ θ=˙xCcos (ψ) + ˙yCsin (ψ) Rcos (ϕ).(8.59) Equations (8.58) and (8.59) will be useful in the discussion of the results obtained with the design of the linear feedback controller. The equations of motion of the ball-plate multibody system are given by the dynamic equations, which are obtained with symbolic computation as explained in Chapter 2, and the nonholonomic constraints (8.50). The following nonlinear index-2 DAE system is obtained: M(x,u)¨x+DT(x)Λ=Q(x,˙x,u),(8.60) 305 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems Cnh (x,˙x) = D(x)˙x=0.(8.61) Note that, in the dynamic equations (8.60), the n×nmass matrix M(x,u) and the n×1 vector of generalized forces Q(x,˙x,u) depend on the vector of control inputs u. The expressions of the equations of motion are overly cumbersome to be shown in this document. 8.4.2 Computation of the reference solution: circular steady motion of the ball The objective of Sect. 8.4 is to develop a linear feedback controller, based on the LQR algorithm, to perform the circular tracking of the ball on the plate. To this end, the first step is to obtain the reference solution, which corresponds to the circular steady motion of the ball. This reference solution can be expressed as: x0(t) = x0 C(t)y0 C(t)ψ0(t)ϕ0(t)θ0(t)T,(8.62) with x0 C(t) = xA+ρcos (ωt), y0 C(t) = yA+ρsin (ωt), ψ0(t) = π 2+ωt, ϕ0(t) = ϕ0,(8.63) θ0(t)=Ωt. Note that Eqs. (8.63) represent a circular trajectory of radius ρ, described by the contact point Cand centered at point A. The Cartesian coordinates of Aare given by xAand yA. To satisfy the nonholonomic constraints, ωand Ω must fulfil: Ω = ρω Rcos (ϕ0).(8.64) Equation (8.64) can be derived by substituting the reference solution given by Eq. (8.62) in Eq. (8.59). From the dynamic equations particularized for Eq. (8.62): Mx0,u0¨x0+DTx0Λ0=Qx0,˙x0,u0,(8.65) the Lagrange multipliers in the reference solution Λ0are obtained: Λ0(t) = Λ0 1(t) Λ0 2(t)T,(8.66) with Λ0 1=ρmeω2cos (ωt) + megcos α0(t)sin β0(t),(8.67) Λ0 2=ρmeω2sin (ωt)−megsin α0(t).(8.68) In the same way, the control inputs in the reference solution, denoted as u0, are obtained from Eq. (8.65): u0(t) = α0(t)β0(t)T,(8.69) 306 8.4. Example of study: ball-plate system where sin α0(t)=(ρmeωR +IΩ cos (ϕ0)) ω megR sin (ωt),(8.70) sin β0(t)=−(ρmeωR +IΩ cos (ϕ0)) ω megRs1−ω2ρmeωR +IΩ cos (ϕ0) megR 2 sin2(ωt) cos (ωt). (8.71) 8.4.3 Linearization along the circular steady motion The linearized equations of motion are obtained using the third numerical linearization approach (Sect. 4.2.3). Choosing as independent coordinates ˜xai in the coordinate partition of Eq. (4.96) the three orientation angles of the sphere: ˜xai =˜ ψ˜ ϕ˜ θT,(8.72) ˜xad =˜xC˜yCT,(8.73) and defining ˜ X=˜xai ˙ ˜xai ˜xad T, the following linear time-varying (LTV) system is obtained: ˙ ˜ X=A(t)˜ X+B(t)˜u.(8.74) Note that the coordinate partition of Eqs. (8.72) and (8.73) is the same as that used with the rolling hoop (Eq. (5.38)) and the rolling toroidal wheel (Eq. (6.54)). The power of the linearization approach allows for obtaining the expressions of matrices A(t) and B(t) analytically: 307 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems A(t) =              0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 a10 0 0 a20 0 0 0a30a40a50 0 a60 0 0 a70 0 0 a8cos (ωt)a9sin (ωt) 0 0 a10 cos (ωt)a11 sin (ωt) 0 0 a12 sin (ωt)−a13 cos (ωt) 0 0 a14 sin (ωt)−a15 cos (ωt) 0 0              ,(8.75) 308 8.4. Example of study: ball-plate system B(t) =              0 0 0 0 0 0 b1(t)b2(t) b3(t)b4(t) b5(t)b6(t) 0 0 0 0              .(8.76) The expressions of the coefficients aj, with j= 1 . . . 15, and bk, with k= 1 . . . 6, are overly long to be included in the document. The system of linear equations (8.74) is given by: ˙ ˜ ψ= ˜ωψ,(8.77) ˙ ˜ ϕ= ˜ωϕ,(8.78) ˙ ˜ θ= ˜ωθ,(8.79) ˙ ˜ωψ=a1˜ ψ+a2˜ωϕ+b1(t) ˜α+b2(t)˜ β, (8.80) ˙ ˜ωϕ=a3˜ ϕ+a4˜ωψ+a5˜ωθ+b3(t) ˜α+b4(t)˜ β, (8.81) ˙ ˜ωθ=a6˜ ψ+a7˜ωϕ+b5(t) ˜α+b6(t)˜ β, (8.82) ˙ ˜xC=a8cos (ωt)˜ ψ+a9sin (ωt)˜ ϕ+a10 cos (ωt) ˜ωϕ+a11 sin (ωt) ˜ωθ,(8.83) ˙ ˜yC=a12 sin (ωt)˜ ψ−a13 cos (ωt)˜ ϕ+a14 sin (ωt) ˜ωϕ−a15 cos (ωt) ˜ωθ.(8.84) 8.4.4 Controllability of the linear time-varying system The controllability of a linear time-varying (LTV) system along a reference trajectory can be determined with the controllability Gramian [318, 339]: Wc(t0, t1) = Zt1 t0 Φ(t1, τ)B(τ)BT(τ)ΦT(t1, τ) dτ, (8.85) where Φ(t1, τ) is the state transition matrix, which can be computed for LTV systems as explained in Sect. 6.1, by means of the fundamental matrix solution or other methods as the Peano-Baker series [268]. The linear system is controllable if the N×Nmatrix Wc(t0, t1) is non-singular on the interval t∈[t0, t1]. For a multibody system with ngeneralized coordinates, mholonomic constraints and l nonholonomic constraints, N= 2 (n−m)−l. In case the linear system verifies rank (Wc(t0, t1)) < N, it is rank deficient and therefore uncontrollable. For the ball-plate system, n= 5, m= 0 and l= 2, and thus N= 2n−l= 8. The numerical computation of the controllability Gramian in Eq. (8.85), for t0= 0 and t1=2π ω, with ω= 1 rad/s, results in rank (Wc(t0, t1)) = 7 (rank deficient), which means that the system is uncontrollable. The state transition matrix in 309 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems Eq. (8.85) has been computed numerically by means of the fundamental matrix solution of the system. Alternatively, for linear time-varying (LTV) systems, the controllability criterion of Silverman and Meadows [340, 341] can be used. By means of this criterion, the total controllability and total observability can be characterized in terms of the matrices A(t), B(t) and their time derivatives. The time-varying matrices A(t) and B(t) must be piecewise differentiable at least N−2 and N−1 times, respectively, which is verified in the case of the linear system of Eq. (8.74). Since the power of the linearization approach allows the analytical computation of these matrices, this result provides an easier way to test the controllability of the ballplate system along the circular steady trajectory. Define the matrix QN−1(t) as follows: QN−1(t) = P0(t)P1(t). . . PN−1(t),(8.86) where Pk(t), with k= 0 . . . N −1, is computed recursively as: Pk(t) = −A(t)Pk−1(t) + ˙ Pk−1(t),(8.87) with P0(t) = B(t). According to [340], the system is completely state controllable at t0if there exists a finite interval [t1, t2], with t1> t0, such that rank (QN−1) = Nfor almost all t∈[t1, t2]. The analytical and numerical computation of QN−1(t) for the linear system of Eq. (8.74) leads to rank (QN−1) = 7, being very ill-conditioned and therefore rank-deficient. This results confirms that the linear system (8.74) is not controllable. Although it has been shown that the system is uncontrollable, the forward dynamics simulation of the ball-plate system is performed. The reference solution and the results of the open-loop and closed-loop systems, in the linear and nonlinear scenarios, are shown. The following numerical values for the ball-plate system and the circular trajectory are considered: me= 0.1 kg, R = 0.05 m, I =2 5meR2, ρ= 0.5 m, ω = 1 rad/s, ϕ0= 5◦.(8.88) From Eq. (8.64), with the numerical values of Eq. (8.88), Ω ≃10.0382 rad/s. In the dynamic simulation, the following initial conditions, at position and velocity levels, are considered: x(0) = xC(0) yC(0) ψ(0) ϕ(0) θ(0) T,(8.89) ˙x(0) = ˙xC(0) ˙yC(0) ˙ ψ(0) ˙ ϕ(0) ˙ θ(0) T,(8.90) with xC(0) = ρ+εx, yC(0) = εy, ψ (0) = π 2+εψ, ϕ(0) = ϕ0+εϕ, θ (0) = εθ, ˙xC(0) = vx+εvx,˙yC(0) = ρω +εvy,˙ ψ(0) = ω+εωψ, ˙ ϕ(0) = εωϕ,˙ θ(0) = Ω + εωθ. (8.91) 310 8.4. Example of study: ball-plate system In Eqs. (8.91), εx,εy,εψ,εϕ,εθ,εvx,εvy,εωψ,εωϕand εωθare initial deviations of the generalized coordinates and velocities with respect to their values in the reference solution. When these initial deviations are set equal to zero, the ball describes the reference solution presented in Eq. (8.62) and (8.63). The following numerical values are considered: εx=−0.1 m, εy= 0.1 m, εψ= 5◦≃0.0873 rad, εϕ= 0◦, εθ= 10◦≃0.1745 rad, εvx= 0.1 m/s, εvy= 0.1 m/s, εωψ= 0.1 rad/s, εωϕ≃3.0383 rad/s, εωθ≃1.7868 rad/s. (8.92) The initial conditions (8.89) and (8.90) satisfy the nonholonomic constraints (8.61). -0.5 0 0.5 1 1.5 xC(m) -0.5 0 0.5 1 yC(m) Reference LQR;L case Open loop;NL case LQR;NL case (a) Trajectory of the contact point C, with ϕ0= 5◦ 0 2 4 6 8 10 12 Time t, (s) -10 -5 0 5 10 15 20 ?(deg) Reference LQR;L case LQR;NL case Open loop;NL case (b) Time evolution of ϕ, for ϕ0= 5◦ Figure 8.2: Closed-loop system with first coordinate partition: trajectory of the ball on the plate, for ϕ0= 5◦. The contact point apparently follows the reference path in both the linear (L) and nonlinear (NL) scenarios, while the reference value of the lean angle ϕ0is not achieved. The results of the open-loop simulation are also shown. Figure 8.2 (a) shows the trajectory of the contact point Cof the ball on the plate, in both the closed-loop and open-loop cases. In the closed-loop case, it can be seen that Capparently follows the reference trajectory in the linear (L) and nonlinear (NL) scenarios. On the other hand, Fig. 8.2 (b) shows that the lean angle ϕis not controllable. The reference value of ϕis not achieved, obtaining a value close to the initial condition ϕ(0). Furthermore, Figs. 8.3 (a) and (b) show the time evolution of the angles ψand θ, respectively, in the closed-loop and open-loop cases. With the LQR controller, the reference values are apparently reached in the linear (L) and nonlinear (NL) scenarios. Nevertheless, the fact that the lean angle ϕis not controllable leads to the following reasoning. From the expression of the nonholonomic constraint (8.59), it can be seen that if this constraint is fulfilled, and ˙xC, ˙yCand ψachieve their references, and ϕdoes not, ˙ θcannot either. However, as shown in Fig. 8.3 (b), 311 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems it seems that θfollow its reference. When the dynamic simulation is performed for a small reference value of the lean angle (in this case ϕ0= 5◦), the lean angle ϕachieves a value close to ϕ0, with cos (ϕ)≃cos (ϕ0). The error made by the controller, which is small in this case, is extended to ˙xC, ˙yC,ψand ˙ θ, leading to negligible differences with respect to their references. This is the reason why xC, yC,ψand θseem to follow their references in Figs. 8.2 and 8.3, when in reality they do not. Consider now a larger reference value of the lean angle, ϕ0= 35◦, which corresponds to Ω ≃12.2077 rad/s. In this case, with the same initial conditions given by Eqs. (8.89) and (8.90), none of the generalized coordinates and their time derivatives follow their references and the errors are perceptible in all the generalized coordinates, at position and velocity level. These results are consistent with the fact that the system is not completely linearly controllable, according to the controllability Gramian and Silverman-Meadows criterion previously applied. Figure 8.4 (a) shows the trajectory of the contact point Cof the ball on the plate, in both the controlled and open-loop cases, for ϕ0= 35◦. In this case, the contact point Cdoes not follow the reference trajectory in both the linear (L) and nonlinear (NL) scenarios. Moreover, Fig. 8.4 (b) shows the time evolution of the lean angle ϕ, which is further from the reference than in Fig. 8.2 (b). Lastly, Fig. 8.5 shows the time evolution of the plate angles αand βin the controlled case. It can be seen that the control inputs do not achieve their reference values α0and β0. 0 2 4 6 8 10 12 Time t, (s) 0 5 10 15 A(rad) Reference LQR;L case LQR;NL case Open loop;NL case (a) Time evolution of the yaw angle ψ, with ϕ0= 5◦ 0 2 4 6 8 10 12 Time t, (s) 0 20 40 60 80 100 120 3(rad) Reference LQR;L case LQR;NL case Open loop;NL case (b) Time evolution of the roll angle θ, for ϕ0= 5◦ Figure 8.3: Time evolution of the angles ψand θ, in the closed-loop and open-loop cases, for ϕ0= 5◦. With the LQR controller, the reference values are apparently reached in the linear (L) and nonlinear (NL) scenarios. Nevertheless, as previously shown, the system is not linearly controllable. 312 8.4. Example of study: ball-plate system -0.5 0 0.5 1 1.5 xC(m) 0 0.5 1 yC(m) Reference LQR;L case Open loop;NL case LQR;NL case (a) Trajectory of the contact point C, with ϕ0= 35◦ 0 2 4 6 8 10 12 Time t, (s) 20 30 40 50 60 ?(deg) Reference LQR;L case LQR;NL case Open loop;NL case (b) Time evolution of ϕ, for ϕ0= 35◦ Figure 8.4: Closed-loop system with first coordinate partition: trajectory of the ball and time evolution of the lean angle, for ϕ0= 35◦. The results show that the system is not controllable. 0 2 4 6 8 10 12 Time t, (s) -30 -20 -10 0 10 20 30 ,,-(deg) ,0 -0 ,LQR -LQR Figure 8.5: Time evolution of the control inputs αand β, for ϕ0= 35◦: comparison with their values in the reference solution α0and β0. 313 Chapter 8. Application of the linearization approaches in the design of linear feedback controllers for nonholonomic multibody systems 320 Chapter 9 Conclusions and future lines of research This chapter summarizes the main contributions of this PhD thesis and presents potential lines of future research. 9.1 Original contributions of this PhD thesis The main objective of this PhD thesis is the development of a methodology that contributes to the modelling of nonholonomic multibody systems, the study of their linear stability and the design of linear feedback controllers. The ultimate focus centers on single-person nonholonomic vehicles, particularly with regard to enhancing safety in urban transportation. By understanding how crucial design parameters affect stability, it is possible to carry out design modifications that result in safer vehicles and reduce accidents. These modifications leads to a broader range of potential users, encompassing the elderly and individuals with physical impairments. From a regulatory standpoint, understanding stability conditions plays an indispensable role in formulating guidelines for these vehicles in urban settings. Furthermore, efforts concentrated on increasing safety also propels the adoption of electric single-person vehicles, thereby contributing to decreased environmental degradation and pollution reduction objectives. The methodology presented in this study consists of the following phases. First, multibody models of the nonholonomic systems are developed. Subsequently, in order to carry out the linear stability analyses of these multibody systems, efficient and accurate linearization approaches are proposed. Next, the linear stability of the multibody systems is assessed by performing detailed eigenvalues sensitivity analysis. Lastly, linear feedback controllers can be designed by using the linearized equations of motion. These phases have been addressed in detail in the chapters of this PhD thesis. The main contributions of each chapter are highlighted below: 321 Chapter 9. Conclusions and future lines of research - Chapters 1 and 2 do not contain original contributions, but they established the foundation for the development of subsequent chapters. Chapter 1 is the introductory chapter and carried out a description and brief state-of-the-art review of Multibody System Dynamics. Emphasis was placed on nonholonomic multibody systems, exploring classical nonholonomic systems, such as the Chaplygin sleigh or the rolling disk, as well as highly mobile nonholonomic multibody systems as the bicycle, the waveboard or the electric kickscooter. A summary with the main linearization procedures developed in the context of Multibody System Dynamics was included, since the linearization of constrained multibody systems is a key matter in this study, required for conducting the linear stability analyses. Chapter 2 is devoted to the modelling of constrained multibody systems. This chapter summarized the basics of rigid body kinematics and dynamics, using the Multibody System Dynamics formalism and with the aim of computational implementation, to develop nonholonomic multibody models from the broadest possible perspective, with the sole particularity of employing nonholonomic constraints linear in velocities. The nonholonomic constraints arised from modelling rolling without slipping. The main tools, used in the following chapters of this study, were reviewed, and the equations of motion of constrained multibody systems were derived using the Newton-Euler formalism. - In Chapter 3, a comprehensive overview of the most notable linearization approaches in Multibody System Dynamics was presented. The procedures were illustrated by using a common notation and classified into four groups. To the best knowledge of the author, this has not been done in any previous work and can be considered an useful contribution for the multibody dynamics community. The classification of the approaches was performed by considering the initial form of the equations of motion and the structure of the linearized equations, which is highly dependent on the use of a redundant or a minimal set of generalized coordinates. A clarification on how the procedures handle the dependent coordinates in the linearization, which is a fundamental issue for correctly obtaining the linearized equations of motion, was also included for each of the approaches. This chapter provided context to the previous work done in linearization of constrained multibody systems before this PhD thesis. - Chapter 4 contains the main contribution of this PhD thesis. In this chapter, five novel linearization approaches for general constrained multibody systems were developed, which constitutes a significant contribution for the multibody dynamics community. By using these approaches, the linearization can be performed along predefined trajectories, which sometimes may be expressed analytically or, in general, numerically. Among these five approaches, two are classified as symbolic procedures, as the linearized equations can only be computed using symbolic toolboxes. The remaining three approaches, labeled as numerical procedures, allow for both symbolic and numerical implementation. A detailed step-by-step explanation of the process for deriving the linearized equations of 322 9.1. Original contributions of this PhD thesis motion was provided. To validate the correctness of these approaches, the wellacknowledged benchmark bicycle model proposed by Meijaard et al. [168] was utilized. The results demonstrated that the numerical linearization procedures are accurate. The eigenvalues obtained from the linear stability analysis of the steady forward motion of the bicycle multibody model coincided with those of the bicycle benchmark of Meijaard et al. [168] in the first fourteen significant digits. Furthermore, the Jacobian matrices derived with the numerical procedures were identical to those of the symbolic procedures, with the same first 32 significant digits. In terms of efficiency, the numerical linearization approaches are more efficient than the symbolic ones. Given that the symbolic approaches involve a very high computational cost and cannot be handled by most of the personal computers, the numerical procedures are much more attractive. In the particular case of the bicycle multibody model, the average time of computation of the Jacobian matrices using the numerical approaches was 0.0005 s. The assessment was carried out by using a computer HP with Intel(R) Core(TM) i7-6700HQ 2.6 GHz and 12 GB of RAM. Moreover, the power of the numerical approaches allowed for generating these Jacobian matrices in terms of the geometric and dynamic parameters of the bicycle, being valid for any constrained multibody system, regardless of size and complexity. Despite presenting a greater number of preparatory steps, the use of the third numerical linearization approach (Sect. 4.2.3) is recommended, as it leads to the minimum number of linearized equations of motion and spurious eigenvalues in the performance of linear stability analyses. For a multibody system with ngeneralized coordinates, mholonomic constraints and lnonholonomic constraints, a (2n−2m−l)×(2n−2m−l) Jacobian matrix is obtained. In this case, l spurious eigenvalues are obtained, associated with the linearized nonholonomic constraints. In the particular case of a holonomic multibody system, no spurious eigenvalues are obtained using this approach. Lastly, another contribution of this chapter is the demonstration of the importance of considering the dependent coordinates in the linearization, which is illustrated with the proposed linearization approaches. As shown with the bicycle, in some multibody systems with a low number of coordinates and in presence of symmetries, a selection of independent coordinates can be intuitively found to ignore the dependent coordinates. However, this is not the usual case, and the systematic linearization of constrained multibody systems requires considering these dependent coordinates. - Chapter 5 addresses the main objective of this PhD thesis, which is the performance of detailed stability analyses of nonholonomic multibody systems. The paramount contribution of this chapter is conducting the linear stability analyses of these systems in a methodological way. This is possible due to the use of the efficient and powerful linearization approaches presented in Chapter 4, which allows for generating the Jacobian matrix of the system under study in 323 Chapter 9. Conclusions and future lines of research terms of the main design parameters. The eigenvalues of the multibody system are parameterized in terms of the geometric and dynamic parameters that influence the stability, as the aspect ratios of the toroidal wheels, steering angle, trail, positions of centres of mass, etc. Therefore, by varying these parameters over a wide range of values, a detailed sensitivity analysis can be carried out to understand changes in the vehicle behavior. The linear stability analyses performed in this chapter correspond to relative equilibria that lead to linear-time invariant (LTI) systems. In particular, the linear stability of the steady forward motion of the nonholonomic systems was studied. A LTI system was obtained since the angular coordinate of the nonholonomic system that evolves periodically is cyclic and is absent from the dynamic equations. Therefore, the study of the linear stability was directly performed by means of the computation of the eigenvalues of the system. The linear stability of highly mobile nonholonomic multibody systems, as the bicycle, the waveboard or the electric kickscooter, was studied. This led to original contributions for each of the systems. First, in the case of the bicycle, the linearized equations of motion of a bicycle multibody model with elliptic toroidal wheels were obtained, along the steady forward motion, as a function of the geometric, dynamic and wheels’ parameters. The results of the linear stability analysis showed that increasing the rear torus aspect ratio broadens the self-stability velocity range, and increasing the front torus aspect ratio has the opposite effect, shrinking the width of the stability regions. When the toroidal geometry is considered in both wheels, the destabilizing effect of the front torus aspect ratio rules over the stabilizing effect of the rear aspect ratio. Furthermore, the consideration of an elliptic cross-section in the wheels results in a stabilizing effect, and the self-stability velocity range increases as the steering axis tilt angle increases. Secondly, the linear stability analysis of the waveboard was addressed, to the best knowledge of the authors, for the first time. To this end, the minimal set of linear equations of the system was obtained as a function of the main design parameters. The influence of the forward speed, the casters’ inclination angles and the tori aspect ratios on the stability was shown, concluding that the waveboard is completely unstable along the steady forward motion. The excellent performance of the third numerical linearization approach with the waveboard confirmed its adequacy for the linearization of complex multibody systems with holonomic and nonholonomic constraints. The influence of the design parameters on the external actuations exerted by the rider in the waveboard maneuvering was also studied by means of an inverse dynamic analysis. Lastly, the linear stability of a detailed e-scooter multibody model, with toroidal wheels and suspensions, was addressed for the first time. The reduced linearized equations along the steady forward motion were obtained with the third numerical linearization approach. The linearized equations were derived in terms of the design parameters of the e-scooter. The multibody model presented was validated with the bicycle benchmark of Meijaard et al. [168], which can be 324 9.1. Original contributions of this PhD thesis used with other two-wheeled vehicles as the e-scooter. In contrast to bicycles, considering the numerical values of the SEAT eXS Kickscooter ES2 [1], the escooter multibody model was shown to be completely unstable. The variation of the moments of inertia of the front handlebar led to significant changes in the stability of these vehicles. In particular, a reduction of the inertia tensor of the front handlebar (associated, for instance, with a removal of the battery from the handlebar) gave rise to the appearance of a self-stability velocity range. Next, the influence of the rider model, considered by modifying the inertia tensor of body B, was studied, and no self-stability velocity range was obtained. The combined effect of the inertia tensor of body Bfound in Paudel et al. [2], together with the reduction of the inertia tensor of the handlebar, led to a evolution of the eigenvalues with the forward speed qualitatively similar to that of Ref. [2]. Furthermore, the inclusion of toroidal wheels resulted in no self-stability velocity range. Regarding the geometry of the cross-section of the toroidal wheel, a lower value of the ellipse aspect ratio led to a stabilizing effect. Similarly, the unstable behavior remained in spite of the modification of the steer axis tilt angle. Lastly, the introduction of the rear and front suspensions led to the emergence of two complex conjugate pairs of eigenvalues. The real parts of these eigenvalues, which were shown to be independent of the forward speed, are highly dependent on the values of the damping coefficients of the suspensions, while the imaginary parts are highly dependent on the stiffness constants of the suspensions’ springs and the mass of the vehicle. The results presented with the e-scooter multibody model demonstrates that the methodology proposed in this PhD thesis can be used for the systematic analysis of the stability of highly mobile nonholonomic multibody systems. - The main contribution of Chapter 6 is the application of the linearization approaches developed in Chapter 4 to assess the linear stability of periodic linear time-varying (LTV) systems. First, basic aspects of LTV systems were presented. Next, a summary of Floquet-Lyapunov theory for the stability analysis of periodic trajectories was included. As example of application, a detailed and comprehensive linear stability analysis of the rolling toroidal wheel was performed. This study was not found in any previous work in the literature of nonholonomic systems. The circular steady motion and the linearized equations of motion along this relative equilibrium were obtained, for both the solid and hollow tori. By virtue of Floquet’s theorem, the periodic time-dependent Jacobian matrix was transformed into a Jacobian matrix with constant coefficients, which allows for the computation of the eigenvalues. These eigenvalues were obtained as a function of the geometric parameters of the rolling toroidal wheel. The change of variable to carry out this transformation was intuitively found due to the symmetry of the reference solution. Finally, the stability results were validated using the Floquet theory for periodic linear differential equations, by computing the Floquet multipliers of the monodromy matrix of the system. The results confirmed that the third numerical linearization procedure is suitable for performing the linearization along reference solutions that lead to 325 Chapter 9. Conclusions and future lines of research linear time-varying (LTV) systems, as the circular steady motion of the rolling wheel. - Chapter 7 addresses the application of the proposed linearization approaches in multiphysics scenarios. The main contribution of this chapter is the development of a linearization approach for hydraulically driven multibody systems. The approach was tested with a hydraulically actuated four-bar mechanism. A comparison with a conventional counterpart linearization procedure was included. In terms of accuracy, the conventional counterpart approach provides an approximation of the Jacobian matrix, its accuracy being strongly influenced by the step size used in the computation of the finite differences. In contrast, the proposed approach leads to the exact linearized equations of motion and exhibits an excellent power, generating the linearized equations of motion as a function of the geometric, dynamic and hydraulic parameters of the multibody system under study. The proposed approach is also more efficient, with a computation time approximately three times lower for the four-bar mechanism. The average computation time for the Jacobian matrix of the four-bar mechanism was found to be 0.005 s with the proposed approach and 0.016 s with the counterpart procedure. The assessment was carried out by using a computer HP with Intel(R) Core(TM) i7-6700HQ 2.6 GHz and 12 GB of RAM. Lastly, a significant advantage of the procedure developed in this chapter is that it allows for obtaining the maximum possible reduction in the number of linearized equations of motion for hydraulically actuated multibody systems with holononomic and nonholonomic constraints. - The main contribution of Chapter 8 is the application of the linearization approaches developed in Chapter 4 in the design of linear feedback controllers for nonholonomic multibody systems. In particular, the third numerical linearization approach was extended to the case of controlled multibody systems, obtaining the expression of the input matrix and deriving a continuos and discrete-time representation of the linear system. A linear feedback controller, designed by using an extended optimal control method as the Linear-Quadratic Regulator (LQR), was developed. The use of the third numerical linearization approach of Chapter 4 also led to original contributions in the optimal control of the ball-plate system, which is a well-known nonholonomic benchmark problem in the context of nonprehensile manipulation. First, the circular steady motion of the rolling ball on the plate was derived analytically, obtaining the time evolution of all the states and control inputs. Next, the exact linearized equations of motion along this reference solution were obtained. The coefficient and input matrices were derived analytically. This significantly facilitates the application of linear controllability criteria for LTV systems, such as the Silverman and Meadows criterion [340, 341]. In contrast to Refs. [325, 331, 333, 334], where a simplified model of the ball-plate system can be found, the results presented in this study have been obtained for a three-dimensional model. Moreover, in Refs. [325, 331, 333, 334], the equa- 326 9.2. Possible future lines of research tions of motion were linearized with respect to the equilibrium position of the sphere on the horizontal plate. Therefore, the control input angles αand βwere very small, so that the linear model remained valid. In contrast, the linearized equations of motion obtained in this study are valid for a circular steady motion with an arbitrary radius, without any restriction in the values of the orientation angles αand βof the plate. 9.2 Possible future lines of research The methodology developed in this PhD thesis provides interesting research opportunities for future works. Some of the possible future research directions are listed below. 1. Application of the novel linearization approaches in flexible multibody systems. The scope of this study has been limited to rigid multibody systems. Nevertheless, the linearization approaches developed in Chapter 4 could be applied to linearize the equations of motion of flexible multibody systems. This could be of interest, for example, to perform linear stability analyses or develop linear feedback controllers. 2. Linear stability analysis of more advanced multibody models. Despite the methodology proposed in this PhD thesis has been successfully applied with highly mobile nonholonomic multibody systems, such as the electric kickscooter with toroidal wheels and suspensions, more advanced multibody models could be considered. The vehicle models could incorporate, for instance, advancements in tires or rider modelling. 3. Linear stability analysis of more complex trajectories. The linear stability analysis of several linear time-invariant (LTI) systems was carried out in this study. Furthermore, the linear stability of periodic linear time-varying (LTV) systems was also addressed with the rolling toroidal wheel as a case study. Nevertheless, the stability of complex nonholonomic vehicles, such as the bicycle or the electric kickscooter with toroidal wheels and suspensions, along more complex trajectories as the circular steady motion, could be analyzed in future work. 4. Development of state observers based on the linearized equations of motion. The linearization approaches of Chapter 4 could be used to develop state observers based on the linearized equations of motions, like Kalman filters. Some examples of works in this line are those of Cuadrado et al. [225, 226], Jaiswal et al. [228] and Qhadim et al. [299]. In this application, a linear approximation of the dynamics of the multibody system is required at each instant of time. The Jacobian matrices obtained with the proposed linearization approaches must be updated throughout the reference solution. The proposed procedures are suitable for this application, as they can be 327 Chapter 9. Conclusions and future lines of research used to provide the linearization along an arbitrary reference solution. Furthermore, they have already demonstrated an excellent performance with linear time-varying (LTV) systems, as shown in Chapters 6 with the rolling toroidal wheel and in Chapter 8 with the ball-plate system. 5. Development of linear feedback controllers for nonholonomic vehicles. In Chapter 8, a Linear-Quadratic Regulator (LQR) controller was developed for the optimal control of the ball-plate system, which is a well-acknowledged nonholonomic problem. In future works, linear feedback controllers for complex nonholonomic multibody systems, such as the bicycle or the e-scooter, could be designed following the methodology proposed in this PhD thesis. The development of controllers that ensure the vehicle’s operation within a stable domain is crucial from the point of view of safety in urban transportation. 6. Development of controllers using other control techniques. Despite the LQR controller has been chosen as feedback controller in this study, the linearized equations obtained with the linearization approaches developed in Chapter 4, whether in continuous or discrete form, could be used with other widely extended control techniques, such as model predictive control [337, 338]. 328 Bibliography [1] “Electric kickscooter’s webpage.” https://www. seat.co.uk/about-seat/news-events/corporate/ seat-launches-electric-urban-mobility-vehicles. Accessed: July 3, 2023. [2] M. Paudel and F. Fah Yap, “Front steering design guidelines formulation for e-scooters considering the influence of sitting and standing riders on self-stability and safety performance,” Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, vol. 235, no. 9, pp. 2551–2567, 2021. [3] P. Flores, Concepts and formulations for spatial multibody dynamics. Springer, 2015. [4] P. E. Nikravesh, Planar Multibody Dynamics: Formulation, Programming with MATLAB®, and Applications. CRC press, 2018. [5] W.-H. Zhu, J.-C. Piedboeuf, and Y. Gonthier, “A dynamics formulation of general constrained robots,” Multibody System Dynamics, vol. 16, no. 1, pp. 37–54, 2006. [6] A. Jain, Robot and multibody dynamics: analysis and algorithms. Springer Science & Business Media, 2010. [7] P. Flores, “A methodology for quantifying the kinematic position errors due to manufacturing and assembly tolerances,” Strojniˇski vestnik-Journal of Mechanical Engineering, vol. 57, no. 6, pp. 457–467, 2011. [8] S. Sanchez-Salinas, C. Nunez-Torres, J. Lopez-Martinez, D. Garcia-Vallejo, and J. M. Muyor, “Design and analysis of a constant-force bench press,” Mechanism and Machine Theory, vol. 142, p. 103612, 2019. [9] M. Silva and J. Ambr´osio, “Kinematic data consistency in the inverse dynamic analysis of biomechanical systems,” Multibody System Dynamics, vol. 8, pp. 219–239, 2002. 329 Bibliography [97] W. Schiehlen, “Dynamics and control of nonholonomic mobile robot systems,” IFAC Proceedings Volumes, vol. 27, no. 14, pp. 329–334, 1994. [98] A. V. Borisov and I. S. Mamaev, “On the history of the development of the nonholonomic dynamics,” arXiv preprint nlin/0502040, 2005. [99] A. M. Bloch, J. E. Marsden, and D. V. Zenkov, “Nonholonomic dynamics,” Notices of the AMS, vol. 52, no. 3, pp. 324–333, 2005. [100] M. De Leon, “A historical review on nonholomic mechanics,” Revista de la Real Academia de Ciencias Exactas, Fisicas y Naturales. Serie A. Matematicas, vol. 106, no. 1, pp. 191–224, 2012. [101] “Rattleback’s page in wikipedia.” https://en.wikipedia.org/wiki/ Rattleback. Accessed: July 3, 2023. [102] M. J. Coleman and P. J. Holmes, “Motions and stability of a piecewise holonomic system: the discrete Chaplygin sleigh,” Regular and Chaotic Dynamics, vol. 4, no. 2, pp. 55–77, 1999. [103] A. M. Bloch, “Asymptotic Hamiltonian dynamics: the Toda lattice, the three-wave interaction and the non-holonomic Chaplygin sleigh,” Physica D: Nonlinear Phenomena, vol. 141, no. 3-4, pp. 297–315, 2000. [104] I. A. Bizyaev, A. V. Borisov, and I. S. Mamaev, “The Chaplygin sleigh with parametric excitation: chaotic dynamics and nonholonomic acceleration,” Regular and Chaotic Dynamics, vol. 22, pp. 955–975, 2017. [105] H. Bondi, “The rigid body dynamics of unidirectional spin,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, vol. 405, no. 1829, pp. 265–274, 1986. [106] H. Moffatt and T. Tokieda, “Celt reversals: a prototype of chiral dynamics,” Proceedings of the Royal Society of Edinburgh Section A: Mathematics, vol. 138, no. 2, pp. 361–368, 2008. [107] A. V. Borisov and I. S. Mamaev, “Strange attractors in rattleback dynamics,” Physics-Uspekhi, vol. 46, no. 4, p. 393, 2003. [108] A. V. Borisov, A. O. Kazakov, and S. P. Kuznetsov, “Nonlinear dynamics of the rattleback: a nonholonomic model,” Physics-Uspekhi, vol. 57, no. 5, p. 453, 2014. [109] L. Euler, “De minimis oscillationibus corporum tam rigidorum quam exililium, methodus nova et facilis,” Commentarii Academiae Scientiarum Imperialis Petropolitanae, vol. 7, pp. 99–122, 1734. [110] S.-D. Poisson, A treatise of mechanics, vol. 2. Longman and Company, 1842. 336 Bibliography [111] G. Slesser, “Notes on rigid dynamics,” Quart. J. of Math, vol. 4, p. 6567, 1861. [112] A. Vierkandt, “On sliding and rolling motion,” Monatshefte der Mathematik and Physik, vol. 3, pp. 31–54 and 97–134 (in German), 1892. [113] S. A. Chaplygin, “On a motion of a heavy body of revolution on a horizontal plane,” Regular and Chaotic Dynamics, vol. 7, no. 2, pp. 119–130, 2002. [114] P. Appell, “Sur l’int´egration des ´equations du mouvement d’un corps pesant de r´evolution roulant par une arˆete circulaire sur un plan horizontal; cas particulier du cerceau,” Rendiconti del Circolo Matematico di Palermo, vol. 14, pp. 1–6, 1900. [115] E. Gallop, “On the rise of a spinning top,” Trans. Camb. Phil. Soc, vol. 19, pp. 356–373, 1904. [116] M. Carvallo, Theory of motion of a monocycle and a bicycle: Part 1. Hoop and monocycle, vol. 5. 1900. [117] E. J. Routh, The Advanced Part of A Treatise on the Dynamics of a System of Rigid Bodies: Being Part II. of a Treatise on the Whole Subject, vol. 229. Dover Publications, 1905. [118] I. M. Mindlin, On the stability of motion of a heavy body of revolution on a horizontal plane, vol. 4. 1964. [119] I. Mindlin and G. Pozharitskii, “On the stability of steady motions of a heavy body of revolution on an absolutely rough horizontal plane,” Journal of Applied Mathematics and Mechanics, vol. 29, no. 4, pp. 879–883, 1965. [120] A. P. Duvakin, On the stability of motions of a disk, vol. 5. 1965. [121] A. Karapetyan and V. Rumyantsev, “Stability of conservative and dissipative systems,” Itogi Nauki Tekh. Ser. Obshch. Mekh, vol. 6, 1983. [122] S. N. Kolesnikov, “Rolling of a disk along a horizontal plane,” Vestnik Moskovskogo Universiteta. Seriya 1. Matematika. Mekhanika, no. 2, pp. 55– 60, 1985. [123] Y. N. Fedorov, “On disk rolling on absolutely rough surface,” Izv. Akad. Nauk SSSR. Mekh. Tverd. Tela, vol. 4, pp. 67–75, 1987. [124] V. Kozlov, “On motion of disk on an inclined plain.,” Proc. of Russian Acad. of Sciences, Rigid body mech., no. 5, pp. 29–35, 1996. [125] O. O’reilly, “The dynamies of rolling disks and sliding disks,” Nonlinear Dynamics, vol. 10, no. 3, pp. 287–305, 1996. [126] R. Cushman, J. Hermans, and D. Kemppainen, “The rolling disc,” in Nonlinear dynamical systems and chaos, pp. 21–60, Springer, 1996. 337 Bibliography [127] A. S. Kuleshov, “On stationary rollings of disk on rough surface,” Appl. Math. & Mech., vol. 65, no. 1, pp. 173–175, 2001. [128] A. V. Borisov, I. S. Mamaev, and A. A. Kilin, “Dynamics of rolling disk,” Regular and Chaotic Dynamics, vol. 8, no. 2, pp. 201–212, 2003. [129] D. V. Zenkov, A. M. Bloch, and J. E. Marsden, “The energy-momentum method for the stability of non-holonomic systems,” Dynamics and stability of systems, vol. 13, no. 2, pp. 123–165, 1998. [130] P. Paris and L. Zhang, “A disk rolling on a horizontal surface without slip,” Mathematical and computer modelling, vol. 36, no. 7-8, pp. 855–860, 2002. [131] R. H. Cushman, H. Duistermaat, and J. Sniatycki, Geometry of nonholonomically constrained systems, vol. 26. World Scientific, 2009. [132] M. Przybylska and S. Rauch-Wojciechowski, “Dynamics of a rolling and sliding disk in a plane. Asymptotic solutions, stability and numerical simulations,” Regular and Chaotic Dynamics, vol. 21, no. 2, pp. 204–231, 2016. [133] A. A. Kilin and E. N. Pivovarova, “Dynamics of an unbalanced disk with a single nonholonomic constraint,” Regular and Chaotic Dynamics, vol. 28, no. 1, pp. 78–106, 2023. [134] A. Borisov and I. Mamaev, “The rolling motion of a rigid body on a plane and a sphere,” Hierarchy of dynamics Regul. Chaotic Dyn, vol. 7, pp. 177– 200, 2002. [135] J. Hermans, “A symmetric sphere rolling on a surface,” Nonlinearity, vol. 8, no. 4, p. 493, 1995. [136] A. V. Borisov, I. S. Mamaev, and I. A. Bizyaev, “The hierarchy of dynamics of a rigid body rolling without slipping and spinning on a plane and a sphere,” Regular and Chaotic Dynamics, vol. 18, no. 3, pp. 277–328, 2013. [137] A. V. Borisov, I. S. Mamaev, and A. A. Kilin, “The rolling motion of a ball on a surface. New integrals and hierarchy of dynamics,” arXiv preprint nlin/0303024, 2003. [138] A. V. Borisov, I. Mamaev, I. Bizyaev, et al., “The Jacobi integral in nonholonomic mechanics,” Regular and Chaotic Dynamics, vol. 20, no. 3, pp. 383– 400, 2015. [139] A. V. Borisov, T. B. Ivanova, Y. L. Karavaev, and I. S. Mamaev, “Theoretical and experimental investigations of the rolling of a ball on a rotating plane (turntable),” European Journal of Physics, vol. 39, no. 6, p. 065001, 2018. [140] M. Antali and G. Stepan, “Slipping–rolling transitions of a body with two contact points,” Nonlinear Dynamics, vol. 107, no. 2, pp. 1511–1528, 2022. 338 Bibliography [141] R. M. Murray and S. S. Sastry, “Nonholonomic motion planning: Steering using sinusoids,” IEEE transactions on Automatic Control, vol. 38, no. 5, pp. 700–716, 1993. [142] A. M. Bloch and P. E. Crouch, “Kinematics and dynamics of nonholonomic control systems on Riemannian manifolds,” in Proceedings of the 32nd IEEE Conference on Decision and Control, pp. 1–5, 1992. [143] M. Hubbard, “Mechanics of skate boards,” Journal of Applied Mechanics, vol. 46, p. 931, 1979. [144] A. V. Kremnev and A. S. Kuleshov, “Nonlinear dynamics and stability of the skateboard,” Discrete & Continuous Dynamical Systems-S, vol. 3, no. 1, p. 85, 2010. [145] J. Ostrowski, A. Lewis, R. Murray, and J. Burdick, “Nonholonomic mechanics and locomotion: the snakeboard example,” in Proceedings of the 1994 IEEE International Conference on Robotics and Automation, pp. 2391– 2397, IEEE, 1994. [146] J. P. Ostrowski, J. P. Desai, and V. Kumar, “Optimal gait selection for nonholonomic locomotion systems,” The International journal of robotics research, vol. 19, no. 3, pp. 225–237, 2000. [147] A. S. Kuleshov, “Further development of the mathematical model of a snakeboard,” Regular and Chaotic Dynamics, vol. 12, no. 3, pp. 321–334, 2007. [148] E. Shammas and M. de Oliveira, “An analytic motion planning solution for the snakeboard,” in Robotics: Science and Systems (RSS), p. 297, MIT Press Cambridge, MA, 2012. [149] T. Dear, S. D. Kelly, M. Travers, and H. Choset, “Snakeboard motion planning with viscous friction and skidding,” in 2015 IEEE International Conference on Robotics and Automation (ICRA), pp. 670–675, IEEE, 2015. [150] I. I. Kosenko, M. S. Loginova, Y. P. Obraztsov, and M. S. Stavrovskaya, “Multibody systems dynamics: Modelica implementation and bond graph representation,” in Proceedings of the 5th International Modelica Conference, arsenal research, Vienna, Austria, pp. 213–223, 2006. [151] F. Boyer and S. Ali, “Recursive inverse dynamics of mobile multibody systems with joints and wheels,” IEEE Transactions on Robotics, vol. 27, no. 2, pp. 215–228, 2011. [152] A. M¨uller and Z. Terze, “Differential-geometric modelling and dynamic simulation of multibody systems,” Strojarstvo: Journal for Theory and Application in Mechanical Engineering, vol. 51, no. 6, pp. 597–612, 2009. 339 Bibliography [153] Z. Terze and J. Naudet, “Structure of optimized generalized coordinates partitioned vectors for holonomic and non-holonomic systems,” Multibody System Dynamics, vol. 24, pp. 203–218, 2010. [154] K. Kinugasa, M. Ishikawa, Y. Sugimoto, and K. Osuka, “Modeling and control of casterboard robot,” IFAC Proceedings Volumes, vol. 46, no. 23, pp. 785–790, 2013. [155] A. Agrawal, H. M. Zaini, T. Dear, and H. Choset, “Experimental gait analysis of waveboard locomotion,” in ASME 2016 Dynamic Systems and Control Conference, pp. V002T22A011–V002T22A011, American Society of Mechanical Engineers, 2016. [156] B. Su, T. Wang, J. Wang, and S. Kuang, “Kinematic mechanism and path planning of the essboard,” Science China Technological Sciences, vol. 56, no. 6, pp. 1499–1516, 2013. [157] T. Wang, B. Su, S. Kuang, and J. Wang, “On kinematic mechanism of a two-wheel skateboard: The essboard,” Journal of Mechanisms and Robotics, vol. 5, no. 3, p. 034503, 2013. [158] Z. Zeng, D. Chen, T. Zhang, and Y. Guan, “Kinematic modelling and analysis of an ess-board-like robot,” in 2018 IEEE International Conference on Robotics and Biomimetics (ROBIO), pp. 1371–1376, IEEE, 2018. [159] S. Ito, K. Niwa, S. Sugiura, and R. Morita, “An autonomous mobile robot with passive wheels propelled by a single motor,” Robotics and Autonomous Systems, vol. 122, p. 103310, 2019. [160] S. Ito, S. Sugiura, Y. Masuda, S. Nohara, and R. Morita, “Mechanism and control of a one-actuator mobile robot incorporating a torque limiter,” Journal of Intelligent & Robotic Systems, vol. 97, no. 2, pp. 431–448, 2020. [161] R. Fukai, K. Yagi, and Y. Mori, “Dynamic model for using casterboard by a humanoid robot,” Advanced Robotics, vol. 34, no. 10, pp. 648–660, 2020. [162] M. M. Gadzhiev, A. S. Kuleshov, and A. I. Bukanov, “Geometric constraints in the problem of motion of a caster board,” Journal of Mathematical Sciences, pp. 1–5, 2020. [163] A. DasGupta, “Dynamics of a waveboard simplified,” Proceedings of the Royal Society A, vol. 476, no. 2244, p. 20200486, 2020. [164] M. Weiss, P. Dekker, A. Moro, H. Scholz, and M. K. Patel, “On the electrification of road transportation–a review of the environmental, economic, and social performance of electric two-wheelers,” Transportation Research Part D: Transport and Environment, vol. 41, pp. 348–366, 2015. [165] F. J. Whipple, “The stability of the motion of a bicycle,” Quarterly Journal of Pure and Applied Mathematics, vol. 30, no. 120, pp. 312–321, 1899. 340 Bibliography [166] K. J. Astrom, R. E. Klein, and A. Lennartsson, “Bicycle dynamics and control: adapted bicycles for education and research,” IEEE Control Systems Magazine, vol. 25, no. 4, pp. 26–47, 2005. [167] D. J. Limebeer and R. S. Sharp, “Bicycles, motorcycles, and models,” IEEE Control Systems Magazine, vol. 26, no. 5, pp. 34–61, 2006. [168] J. P. Meijaard, J. M. Papadopoulos, A. Ruina, and A. L. Schwab, “Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and review,” Proceedings of the Royal society A: mathematical, physical and engineering sciences, vol. 463, no. 2084, pp. 1955–1982, 2007. [169] J. Kooijman, A. Schwab, and J. P. Meijaard, “Experimental validation of a model of an uncontrolled bicycle,” Multibody System Dynamics, vol. 19, no. 1-2, pp. 115–132, 2008. [170] P. Basu-Mandal, A. Chatterjee, and J. M. Papadopoulos, “Hands-free circular motions of a benchmark bicycle,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 463, no. 2084, pp. 1983–2003, 2007. [171] J. L. Escalona and A. M. Recuero, “A bicycle model for education in multibody dynamics and real-time interactive simulation,” Multibody System Dynamics, vol. 27, no. 3, pp. 383–402, 2012. [172] J. Xiong, N. Wang, and C. Liu, “Stability analysis for the whipple bicycle dynamics,” Multibody System Dynamics, vol. 48, no. 3, pp. 311–335, 2020. [173] J. Xiong, N. Wang, and C. Liu, “Bicycle dynamics and its circular solution on a revolution surface,” Acta Mechanica Sinica, vol. 36, no. 1, pp. 220–233, 2020. [174] E. X. Wang, J. Zou, G. Xue, L. Yuan, Z. Zeng, G. Xue, Y. Liu, G. Zhang, and Q. Fan, “Symbolic derivation of bicycle kinematics with toroidal wheels,” in MATEC Web of Conferences, vol. 34, p. 04001, EDP Sciences, 2015. [175] G. Frosali and F. Ricci, “Kinematics of a bicycle with toroidal wheels,” Communications in Applied and Industrial Mathematics, vol. 3, no. 1, 2012. [176] J. Meijaard and A. Schwab, “Linearized equations for an extended bicycle model,” in III European Conference on Computational Mechanics, pp. 772– 772, Springer, 2006. [177] A. Schwab, J. Meijaard, and J. Kooijman, “Some recent developments in bicycle dynamics,” in Proceedings of the 12th World Congress in Mechanism and Machine Science, pp. 1–6, Russian Academy of Sciences Moscow, Russia, 2007. 341 [Document text truncated for crawler view.]