Dynamic inversion and optimal tracking control on the ball-plate system based on a linearized nonholonomic multibody model

Contents lists available at ScienceDirect Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmt Research paper Dynamic inversion and optimal tracking control on the ball-plate system based on a linearized nonholonomic multibody model A.G. Agúndez a,b,∗, A. Saccon c, D. García-Vallejo a, E. Freire d aDepartment of Mechanical Engineering and Manufacturing, Universidad de Sevilla, Seville, Spain bDepartment of Mechanical Engineering, Lappeenranta University of Technology, Lappeenranta, Finland cDepartment of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands dDepartment of Applied Mathematics II, Universidad de Sevilla, Seville, Spain ARTICLE INFO Keywords: Ball-plate system Trajectory tracking Linearization Multibody dynamics Nonholonomic system ABSTRACT This paper addresses the optimal control of the ball-plate system, a well-known nonholonomic system in the context of nonprehensile manipulation, using a multibody dynamics approach. The trajectory tracking control of a steady-state circular motion of the ball on the plate, for any radius and potentially off-centric with respect to the plate’s pivoting point, is achieved by designing a Linear-Quadratic Regulator. A spatial multibody model of the ball-plate system is considered. A key contribution is the analytical computation of the circular steady motion of the ball by dynamic inversion, including the control actions to achieve this reference solution. This enables the analytical computation of the linearized equations along this reference motion, resulting in a periodic linear time-varying (LTV) system, and the application of linear controllability criteria for LTV systems. A controllable linear system, involving the Cartesian coordinates of the contact point and the yaw angle of the sphere, is obtained using a convenient coordinate partition in the linearization. Compared to existing results on the same problem, closed-loop stability about the desired trajectory is achieved for any radius of the circular trajectory. 1. Introduction 1.1. Background The problem of nonholonomic rolling is usual in nonprehensile manipulation [1] and can be found 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 [1–3] and can be found in several works of the literature [1,4]. The planning and stabilization of the trajectory of a smooth body rolling on the surface of another [5] is usual within nonprehensile manipulation. Specific examples are the butterfly robot [2,6–9] or controlling a sphere rolling on a motion-controlled flat plate (known as ball-plate system [10–13]). The ball-plate system constitutes a well-studied example of 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 without slipping between the ball and the plate. There are numerous ∗Corresponding author at: Department of Mechanical Engineering, Lappeenranta University of Technology, Lappeenranta, Finland. E-mail addresses: [email protected] (A.G. Agúndez), [email protected] (A. Saccon), [email protected] (D. García-Vallejo), [email protected] (E. Freire). https://doi.org/10.1016/j.mechmachtheory.2024.105795 Received 7 June 2024; Received in revised form 11 September 2024; Accepted 12 September 2024 Mechanism and Machine Theory 203 (2024) 105795 Available online 21 September 2024 0094-114X/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ). A.G. Agúndez et al. Nomenclature Symbol Meaning 𝑛Number of generalized coordinates 𝑚Number of holonomic constraints 𝑙Number of nonholonomic constraints 𝑛𝑢Number of actuations 𝑡Time 𝒙Vector of generalized coordinates  𝒙, 𝒙Velocity and acceleration vectors, respectively 𝜦Vector of Lagrange multipliers 𝒖Vector of actuations 𝑰𝑝Identity matrix of size 𝑝 𝟎𝑝×𝑞𝑝×𝑞zero matrix 𝑴(𝒙,𝒖)Mass matrix 𝑸(𝒙, 𝒙,𝒖)Vector of generalized forces 𝑪(𝒙)Vector of holonomic constraints 𝑪𝑛ℎ (𝒙, 𝒙)Vector of nonholonomic constraints 𝑪𝒙(𝒙)Jacobian matrix of the holonomic constraints with respect to the generalized coordinates 𝑩(𝒙)Jacobian matrix of the nonholonomic constraints with respect to the time derivative of the generalized coordinates 𝑫(𝒙)Jacobian matrix of the time derivative of the holonomic constraints, and the nonholonomic constraints, with respect to the time derivative of the generalized coordinates 𝒙𝑎𝑖 Vector of independent coordinates 𝒙𝑎𝑑 Vector of dependent admissible coordinates 𝒙𝑑Vector of dependent coordinates 𝑫𝑎𝑖 (𝒙)Submatrix formed by the columns of matrix 𝑫(𝒙)corresponding to the independent coordinates 𝑫𝑑𝑑 (𝒙)Submatrix formed by the columns of matrix 𝑫(𝒙)corresponding to the dependent coordinates 𝑻(𝒙)Transformation matrix used in the coordinate partitioning method 𝒙0(𝑡)Time evolution of the vector of generalized coordinates in the reference solution  𝒙0(𝑡), 𝒙0(𝑡)Velocity and acceleration vectors, respectively, in the reference solution 𝜦0(𝑡)Time evolution of the vector of Lagrange multipliers in the reference solution 𝒖0(𝑡)Time evolution of the actuations in the reference solution  𝒙Variations of the generalized coordinates with respect to the reference solution  𝒙, 𝒙Variations of the velocity and acceleration vectors, respectively, with respect to the reference solution  𝜦Variations of the Lagrange multipliers with respect to the reference solution  𝒖Variations of the actuations with respect to the reference solution  𝒙𝑎𝑖 Variations of the independent coordinates with respect to the reference solution  𝒙𝑎𝑖, 𝒙𝑎𝑖 Variations of the independent velocities and accelerations, respectively, with respect to the reference solution  𝒙𝑎𝑑 Variations of the dependent admissible coordinates with respect to the reference solution  𝒙𝑑Variations of the dependent coordinates with respect to the reference solution 𝑻0Transformation matrix 𝑻(𝒙)particularized for the reference solution Jacobian matrix in the linearized equations of motion Input matrix in the linearized equations of motion 𝑐Jacobian matrix of the controllable linear system 𝑐Input matrix of the controllable linear system 𝑾𝑐(𝑡0, 𝑡1)Controllability Gramian of the linear time-varying system 𝜱(𝑡1, 𝜏)State-transition matrix of the linear system 𝑸𝑁−1 (𝑡)Controllability matrix in the controllability criterion of Silverman and Meadows 𝑁𝑡Time horizon of the discrete-time finite-horizon LQR  𝒖∗Optimal control sequence in the discrete-time finite-horizon LQR 𝐽∗Quadratic cost function to minimize in the discrete-time finite-horizon LQR  𝑿𝑘State at 𝑡𝑘  𝑿𝑘+1 State at 𝑡𝑘+1 𝑘Jacobian matrix particularized at 𝑡𝑘 Mechanism and Machine Theory 203 (2024) 105795 2 A.G. Agúndez et al. 𝑘Input matrix particularized at 𝑡𝑘  𝒖𝑘Actuations at 𝑡𝑘 𝑲𝑘Time-varying gain matrix of the LQR at 𝑡𝑘 𝑷𝑘Matrix solution of the discrete-time Riccati equation at 𝑡𝑘 𝑸𝑤State-cost weighted matrix 𝑹𝑤Input-cost weighted matrix  𝑿𝑁𝑡State at the final instant of time 𝑸𝑁𝑡State-cost weighted matrix corresponding to the final instant of time 𝐸𝑠 max Maximum acceptable values of the errors in the states 𝐸𝑢 max Maximum acceptable values of the errors in the inputs 𝐶Contact point of the sphere with the plate 𝑚𝑒Mass of the sphere 𝑅Radius of the sphere 𝐼Moment of inertia of the sphere, taken about any diameter 𝜌Radius of the circular trajectory of the ball on the plate 𝜔Constant angular velocity at which the contact point 𝐶describes the circular steady motion 𝛺Rotational speed of the sphere in the circular steady motion works devoted to the ball-plate in the literature, with many different control objectives. A summary with the main ones is presented in Section 1.3, for then highlighting the contributions of this study within this field. 1.2. Formulation of the problem of interest for this investigation This paper focuses on the problem of trajectory tracking of the ball-plate system. In particular, the trajectory tracking control of a steady-state circular motion of the ball on the plate, for any radius and potentially off-centric with respect to the plate’s pivoting point, is addressed. Existing models have limitations due to their reliance on linear approximations around the trivial equilibrium position of the ball on the plate in horizontal position, which restrict the radius of the circular motion and the control input angles. This paper seeks to overcome these limitations by using a nonlinear spatial multibody model of the ball-plate system, to derive a more accurate linearized model that can handle arbitrary radii and unrestricted plate orientation angles. The study also aims to design a linear feedback controller using a novel linearization approach, extending methodologies previously applied to uncontrolled systems to controlled multibody systems. 1.3. Literature survey Numerous works in the literature are devoted to nonprehensile dynamic manipulation and can be found in the survey presented by Ruggiero et al. [1], which focuses on planning and control in the nonprehensile dynamic manipulation domain. The problem of nonholonomic rolling in nonprehensile manipulation tasks was also reviewed by Gutierrez-Giles [4] through two illustrative examples: the robotic hula-hoop and the ballbot system. Many of these works fall within the areas of mobile robotics and robotic manipulation, aimed at motion planning and controller design [14–21]. In particular, Woodruff et al. [5] presented a motion planner and feedback controller for pure-rolling motions between general smooth rigid bodies. The methods were demonstrated for a sphere rolling on a sphere and an ellipsoid rolling on an ellipsoid. Next, Woodruff et al. [22] extended their work, deriving the rolling dynamic equations of smooth rigid bodies, and validated their results by means of simulations of the rolling sphere on a plane. A Linear-Quadratic Regulator (LQR) to stabilize planned trajectories in both simulations (reorientation of a sphere on a plate) and experiments (3-DOF robotic arm manipulating an object) was used in [22]. It should also be noted that there are works in the literature devoted to the dynamic simulation of nonholonomic mechanical systems subject to a general set of holonomic and/or nonholonomic constraints, as Pappalardo et al. [23]. The ball-plate system constitutes a well-known nonholonomic system in the domain of nonprehensile manipulation. Some of the works devoted to the control of the ball-plate system aim to move and maintain the ball in a desired position on the plate. This is the case of Knuplez et al. [24], where the control objective is to maintain the freely rolling ball in a specific position. To this end, a controller for a two-dimensional electromechanical system was synthesized. The study aimed to investigate the system dynamics and perform laboratory experiments with different control approaches based on classical and modern control theory. In Oriolo et al. [25], an iterative control scheme that guarantees asymptotic stability with exponential convergence towards any desired set point was presented. The control objective of achieving a desired position with a desired orientation can also be found in Date et al. [26], where a ball is held between a pair of horizontal plates maintaining rolling contact without slip. The control was achieved by applying horizontal motion to the plates, allowing both translation and rotation. Azizi et al. [27] also addressed the stabilization of a sphere from any initial configuration to any desired position and orientation. To this end, a nonlinear model predictive control (NMPC) was implemented for the point stabilization of a nonholonomic spherical mobile robot actuated by two internal rotors. In the same vein, Mukherjee et al. [10] proposed two different algorithms for the motion planning and reconfiguration of a rolling Mechanism and Machine Theory 203 (2024) 105795 3 A.G. Agúndez et al. sphere, and Serra et al. [28] presented a method to reconfigure the pose (position and orientation) of a sphere rolling on a plate in a nonprehensile manner. The proposed method in Serra et al. [28] was experimentally validated using a humanoid-like robot. In Cheng et al. [29], an LQR control algorithm was designed to adjust the plate’s orientation in order to balance the ball at specific positions with different experiments. A sufficient stabilizability condition for general nonlinear systems was provided by Casagrande et al. [30] and illustrated by means of the ball-plate system. The trajectory tracking control of the ball-plate system can also be found in numerous works in the literature. A controller, based on the backstepping control procedure, for the oscillatory stabilization and circular trajectory tracking of the ball-plate system, was presented by Ker et al. [31]. Fan et al. [11] proposed a hierarchical fuzzy control scheme for trajectory planning and tracking in the ball-plate system, using genetic algorithms. Simulation results demonstrated effective control of the ball without hitting obstacles and in a short time. Kassem et al. [32] addressed the control of a system consisting of a plate mounted on a Stewart platform with six 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, and a comparison of the proposed strategies was included. Furthermore, other control methods, as Model Predictive Control (MPC) [33,34], have been applied to address the trajectory tracking of the ball-plate system. In Fan et al. [35], a novel nonlinear Model Predictive Control (NMPC) strategy 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. [36] also implemented a control algorithm, based on MPC, verified with nonlinear simulations and reference trajectory tracking in a real ball-plate laboratory model. In Bang et al. [37], a method for implementing MPC on a microcontroller was proposed, focusing on enhancing tracking performance in the ball-plate system. A laboratory-built ball-plate system was used in Bang et al. [37]. As highlighted in Section 1.4, in contrast to the previously summarized works [11,31,32,35–37] devoted to the trajectory tracking of the ball-plate system, this paper makes use of a spatial nonholonomic multibody model of the ball-plate and achieves the trajectory tracking of the circular steady motion of the ball on the plate for an arbitrarily large radius. 1.4. Scope and contributions of this study This paper presents several contributions compared to existing results on the same problem, both from the perspective of the ball-plate system control and from the standpoint of the methodology used: 1. First, the trajectory tracking of the circular steady motion of the ball on the plate, for an arbitrarily large radius and potentially off centric with respect to the plate pivoting point, is achieved. Previous works [11,24,31,32,36,37] devoted to control of the ball-plate system made use of a linear model about the trivial equilibrium position of the sphere on the horizontal plate. Therefore, the radius of the circular motion described by the sphere in the circular trajectory tracking and the control input angles of the plate had to be very small for the linear model to remain valid. In contrast, the present work obtains, from a nonlinear multibody model, a linearized model of the ball-plate system along the circular steady motion, which is valid for a circular motion with an arbitrary radius, without any restriction in the values of the plate orientation angles. 2. A key contribution is the analytical computation of the reference solution by dynamic inversion, as well as the necessary control actions on the plate to achieve the circular steady motion. This enables the analytical derivation of the linearized equations of motion as a function of the geometric and dynamic parameters of the ball-plate. 3. In contrast to Refs. [11,24,31,32,36,37], where a simplified model of the ball-plate system can be found, the results presented in this study have been obtained for a spatial nonholonomic multibody model of the ball-plate, which considers the location of the contact point and the orientation angles of the sphere. 4. Regarding the methodology, the study aims to design a linear feedback controller using a novel linearization approach, previously applied to uncontrolled systems, extending its use to controlled multibody systems. The design of the linear feedback controller is carried out using the reduced linearized equations of the system, which are computed using the linearization approach proposed in [38]. This linearization approach has previously been applied to uncontrolled multibody systems, such as the bicycle benchmark [38] or the waveboard [39]. The Linear-Quadratic Regulator (LQR) is chosen as control method, due to its simplicity of implementation, computational efficiency and its presence in several works devoted to the control of the ball-plate system. 5. Lastly, the importance of properly choosing a coordinate partition in the linearization, to achieve a controllable linear system, is demonstrated with the ball-plate system. 1.5. Organization of the paper The following structure is considered. Following the Introduction, Section 2describes the methodology used in the paper. The nonlinear equations of motion and the reduced linearized equations of nonholonomic controlled multibody systems are presented. Next, a discrete-time representation is derived from the continuous-time description of the linear system, and the optimal control of nonholonomic multibody systems is addressed, summarizing the steps to design a discrete-time finite-horizon Linear-Quadratic Regulator. Section 3applies this methodology to the trajectory tracking of the circular steady motion of the ball-plate. A spatial multibody model is described, the circular steady motion of the ball on the plate is computed and the linearized model along this reference solution is derived. The linear controllability is tested and a controllable system of linear equations is obtained. Next, Section 4shows the numerical results of the circular trajectory tracking of the ball-plate system. Finally, Section 5summarizes the main conclusions derived from the present work. Mechanism and Machine Theory 203 (2024) 105795 4 A.G. Agúndez et al. 2. Methodology In this section, the methodology applied subsequently to the ball-plate multibody model is summarized. First, the equations of motion of a controlled multibody system are presented. Next, the linearization approach [38] is applied to controlled multibody systems. By using this approach, the linearization can be performed along predefined trajectories, which sometimes may be expressed analytically or, in general, numerically. From the continuous-time description of the linear system, a discrete-time representation is obtained. Next, a linear feedback controller, designed by using a widely extended optimal control method as the Linear-Quadratic Regulator (LQR), is developed. 2.1. Equations of motion of controlled multibody systems Consider a multibody system with 𝑛generalized coordinates, 𝑚holonomic constraints, 𝑙nonholonomic constraints and 𝑛𝑢 actuations. The equations of motion of the multibody system are given by the index-3 DAE: 𝑴(𝒙,𝒖) 𝒙+𝑫T(𝒙)𝜦=𝑸(𝒙, 𝒙,𝒖),(1) 𝑪(𝒙)=𝟎,(2) 𝑪𝑛ℎ (𝒙, 𝒙)=𝑩(𝒙) 𝒙=𝟎.(3) In Eqs. (1)–(3),𝒙∈R𝑛is the 𝑛× 1 vector of generalized coordinates, 𝜦∈R𝑚+𝑙is the (𝑚+𝑙) × 1 vector of Lagrange multipliers, 𝒖∈R𝑛𝑢is the 𝑛𝑢× 1 vector of actuations, 𝑴(𝒙,𝒖)is the 𝑛×𝑛mass matrix, 𝑸(𝒙, 𝒙,𝒖)is the 𝑛× 1 vector of generalized forces, 𝑪(𝒙) is the 𝑚× 1 vector of holonomic constraints, and 𝑪𝑛ℎ (𝒙, 𝒙)is the 𝑙× 1 vector of nonholonomic constraints, which in this work is assumed to have a linear dependence on the system velocities. The matrices 𝑩(𝒙)and 𝑫(𝒙)are 𝑙×𝑛and (𝑚+𝑙) × 𝑛, respectively, given by: 𝑩(𝒙)=𝜕𝑪𝑛ℎ (𝒙, 𝒙) 𝜕 𝒙,𝑫(𝒙)=(𝑪𝒙(𝒙) 𝑩(𝒙)),(4) where 𝑪𝒙=𝜕𝑪 𝜕𝒙. It should be noted that, in this work, it is assumed that the constraints are scleronomic, and therefore the explicit time dependence does not appear in the constraint equations. In Eqs. (1)–(3), it is assumed that the control actions 𝒖can be, in addition to forces or torques, an angle or position, and therefore, they appear in both the mass matrix 𝑴(𝒙,𝒖)and the generalized force vector 𝑸(𝒙, 𝒙,𝒖). In the particular case where control actions were only forces or torques, they would only appear in the generalized force vector. Finally, in this work, it is assumed that the constraints do not depend on the control actions. 2.2. Reduced linearized equations of motion of controlled nonholonomic multibody systems The third numerical linearization approach proposed in [38] has previously been applied to uncontrolled multibody systems, such as the bicycle benchmark [38] or the waveboard [39]. In this section, the approach is applied to controlled multibody systems, deriving the expression of the input matrix. This linearization approach considers a coordinate partitioning in dependent and independent coordinates [40]. Other methods that allow for coordinate reduction, reducing the DAE system to an ordinary differential equations (ODE) system, are the SVD, QR and Schur decompositions or the Udwadia-Kalaba formulation [41,42]. Consider a known reference solution of the index-3 DAE (1)–(3), defined by the time evolution of the states and their time derivatives: 𝒙0(𝑡), 𝒙0(𝑡), 𝒙0(𝑡), and the Lagrange multipliers 𝜦0(𝑡). Moreover, the actuations in the reference solution are denoted as 𝒖0(𝑡). Therefore, the following holds: 𝑴(𝒙0,𝒖0) 𝒙0+𝑫T(𝒙0)𝜦0=𝑸(𝒙0, 𝒙0,𝒖0),(5) 𝑪(𝒙0)=𝟎,(6) 𝑪𝑛ℎ (𝒙0, 𝒙0)=𝑩(𝒙0) 𝒙0=𝟎,(7) where the time dependence has been omitted for sake of conciseness. The following variations with respect to the reference solution are introduced:  𝒙=𝒙−𝒙0, 𝒙= 𝒙− 𝒙0,  𝒙= 𝒙− 𝒙0, 𝜦=𝜦−𝜦0,(8)  𝒖=𝒖−𝒖0. Following Ref. [38], given a multibody system with 𝑚holonomic constraints, a coordinate partition of the 𝑛-coordinates vector 𝒙into 𝑛−𝑚admissible position coordinates, grouped in 𝒙𝑎, and 𝑚dependent coordinates 𝒙𝑑, is made. The set of admissible position coordinates 𝒙𝑎spans the domain of admissible positions of the multibody system. This partition is represented as 𝒙=(𝒙T 𝑎𝒙T 𝑑)T. Furthermore, the 𝑙nonholonomic constraints allow the partition of the time derivative of the admissible position coordinates  𝒙𝑎 in 𝑙dependent and 𝑛−𝑚−𝑙independent admissible velocities, given by the vectors  𝒙𝑎𝑑 and  𝒙𝑎𝑖, respectively, and therefore Mechanism and Machine Theory 203 (2024) 105795 5 A.G. Agúndez et al.  𝒙𝑎=( 𝒙T 𝑎𝑖  𝒙T 𝑎𝑑 )T. The same partition can be used at position level, with the coordinate vector finally partitioned as in [38]: 𝒙=(𝒙T 𝑎𝑖 𝒙T 𝑎𝑑 𝒙T 𝑑)T.(9) Therefore, the vector of variations  𝒙can be written as:  𝒙=( 𝒙T 𝑎𝑖  𝒙T 𝑎𝑑  𝒙T 𝑑)T.(10) By using this coordinate partition, the following velocity transformation is obtained:  𝒙=𝑻(𝒙) 𝒙𝑎𝑖,(11) where 𝑻is a transformation matrix, given by: 𝑻(𝒙)=(𝑰(𝑛−𝑚−𝑙) −(𝑫𝑑𝑑 (𝒙))−1 𝑫𝑎𝑖 (𝒙)).(12) In Eq. (12),𝑰(𝑛−𝑚−𝑙)is the identity matrix of size 𝑛−𝑚−𝑙;𝑫𝑑𝑑 (𝒙)is a (𝑚+𝑙)-square matrix, formed by the columns of matrix 𝑫(𝒙) corresponding to the dependent coordinates; and 𝑫𝑎𝑖 (𝒙)is a (𝑚+𝑙)×(𝑛−𝑚−𝑙)matrix, defined by the columns of 𝑫(𝒙)associated with the independent coordinates. The selection of independent coordinates can be obtained, in practice, using Gaussian elimination with column pivoting. The transformation matrix of Eq. (12), particularized for the reference solution, is denoted as 𝑻0. The computation of the linearized equations of motion can be performed by following the steps detailed in [38]. Performing the Taylor expansion of the dynamic Eqs. (1), retaining up to first-order terms and using Eq. (5), the linearized dynamic equations are obtained: 𝑴(𝒙0,𝒖0) 𝒙+𝜕(𝑴(𝒙,𝒖) 𝒙0) 𝜕𝒙|||||0  𝒙+𝑫T(𝒙0) 𝜦+𝜕(𝑫T(𝒙)𝜦0) 𝜕𝒙|||||0  𝒙= 𝜕𝑸 𝜕𝒙||||0  𝒙+𝜕𝑸 𝜕 𝒙||||0  𝒙+(𝜕𝑸 𝜕𝒖||||0 −𝜕(𝑴(𝒙,𝒖) 𝒙0) 𝜕𝒖|||||0) 𝒖, (13) where (⋅)|0is henceforth used to indicate that the partial derivatives are evaluated for the reference solution. Note that Eq. (13) includes, with respect to Ref. [38], the term associated with the variations of the actuations  𝒖. Premultiplying Eq. (13) by 𝑻T 0, the number of linearized dynamic equations is reduced, leading to: 𝑻T 0𝑴(𝒙0,𝒖0) 𝒙+𝑻T 0 𝜕(𝑴(𝒙,𝒖) 𝒙0) 𝜕𝒙|||||0  𝒙+𝑻T 0 𝜕(𝑫T(𝒙)𝜦0) 𝜕𝒙|||||0  𝒙= 𝑻T 0(𝜕𝑸 𝜕𝒙||||0  𝒙+𝜕𝑸 𝜕 𝒙||||0  𝒙)+𝑻T 0(𝜕𝑸 𝜕𝒖||||0 −𝜕(𝑴(𝒙,𝒖) 𝒙0) 𝜕𝒖|||||0) 𝒖. (14) Considering Eq. (12), it can be seen that the variations of the Lagrange multipliers,  𝜦, have been eliminated from Eq. (13), since 𝑻T 0𝑫T(𝒙0)=𝟎. Following the steps summarized in Ref. [38], Eq. (14) can be expressed in terms of the variations of the independent coordinates,  𝒙𝑎𝑖, and their time derivatives. The following linear ODE system is obtained:  𝒙𝑎𝑖 =𝑨21  𝒙𝑎𝑖 +𝑨22  𝒙𝑎𝑖 +𝑨23  𝒙𝑎𝑑 +𝑩1 𝒖,(15)  𝒙𝑎𝑑 =𝑨31  𝒙𝑎𝑖 +𝑨32  𝒙𝑎𝑖 +𝑨33  𝒙𝑎𝑑 ,(16) where the expressions of the matrices 𝑨21,𝑨22,𝑨23,𝑨31,𝑨32 and 𝑨33 can be found in detail in Ref. [38]. Matrix 𝑩1is a (𝑛−𝑚−𝑙)×𝑛𝑢 matrix, given by: 𝑩1=𝒎−1 0𝑻T 0(𝜕𝑸 𝜕𝒖||||0 −𝜕(𝑴(𝒙,𝒖) 𝒙0) 𝜕𝒖|||||0).(17) In Eq. (17),𝒎0=𝑻T 0𝑴(𝒙0,𝒖0)𝑻0. Eqs. (15) correspond to the linearized dynamic equations, and Eqs. (16) to the linearized nonholonomic constraints. Defining  𝑿=( 𝒙T 𝑎𝑖  𝒙T 𝑎𝑖  𝒙T 𝑎𝑑 )T, the linearized equations of motion can be written as:   𝑿=(𝑡) 𝑿+(𝑡) 𝒖,(18) where is the 𝑁×𝑁Jacobian matrix presented in Ref. [38], with 𝑁= 2 (𝑛−𝑚)−𝑙, given by: =⎛⎜⎜⎝ 𝟎(𝑛−𝑚−𝑙)𝑰(𝑛−𝑚−𝑙)𝟎(𝑛−𝑚−𝑙)×𝑙 𝑨21 𝑨22 𝑨23 𝑨31 𝑨32 𝑨33 ⎞⎟⎟⎠ .(19) Mechanism and Machine Theory 203 (2024) 105795 6 A.G. Agúndez et al. Regarding the input matrix , its expression is given by: =⎛⎜⎜⎝ 𝟎(𝑛−𝑚−𝑙)×𝑛𝑢 𝑩1 𝟎𝑙×𝑛𝑢⎞⎟⎟⎠ .(20) Note that, in Eq. (18), time dependence in matrices and is considered, since the linearization can be carried out along an arbitrary reference solution. 2.3. Time-varying discrete linear systems In this section, a discrete-time representation is derived from the continuous-time description of the linear system given by Eq. (18). The following definitions are introduced for the sake of simplicity:  𝑿𝑘= 𝑿(𝑡𝑘), 𝑿𝑘+1 = 𝑿(𝑡𝑘+1), 𝑘=(𝑡𝑘),𝑘=(𝑡𝑘),(21)  𝒖𝑘= 𝒖(𝑡𝑘). Let 𝛥𝑡 > 0be the time step between two time instants, 𝑡𝑘and 𝑡𝑘+1. Assuming 𝛥𝑡 ≪ 1, Eq. (18) can be approximated as follows:   𝑿=𝑘 𝑿+𝑘 𝒖𝑘,(22) with  𝑿(𝑡𝑘)= 𝑿𝑘being known. Note that the time-dependent matrices (𝑡),(𝑡)and the input vector  𝒖in Eq. (18) have been approximated by their values at 𝑡𝑘in Eq. (22). Eq. (22) can be solved exactly at (𝑡𝑘, 𝑡𝑘+1), leading to:  𝑿𝑘+1 = e𝑘𝛥𝑡  𝑿𝑘+(e𝑘𝛥𝑡 −𝑰)−1 𝑘𝑘 𝒖𝑘.(23) Eq. (23) is a discrete-time representation of the linear system (18). 2.4. Optimal control of nonholonomic multibody systems: discrete time finite-horizon linear-quadratic regulator The Linear-Quadratic Regulator (LQR) algorithm is a widely used control technique that allows for the optimal control of linear dynamical systems [43]. It combines the concepts of state feedback and optimal control to design controllers that minimize a specified performance criterion. By formulating the control problem as an optimization task, the LQR finds the optimal control law that minimizes a cost function while satisfying system dynamics. In this study, the discrete time finite-horizon LQR is used. Given the linear system (23), the objective is to find the optimal control sequence  𝒖∗= { 𝒖∗ 0, 𝒖∗ 1… 𝒖∗ 𝑁𝑡−1}as a function of the state, which minimizes the following quadratic cost function over the time horizon 𝑁𝑡[43]: 𝐽∗= 𝑿T 𝑁𝑡𝑸𝑁𝑡  𝑿𝑁𝑡+ 𝑁𝑡−1 ∑ 𝑘=0 ( 𝑿T 𝑘𝑸𝑤 𝑿𝑘+ 𝒖T 𝑘𝑹𝑤 𝒖𝑘),(24) starting from the initial condition  𝑿(0)= 𝑿0. Matrix 𝑸𝑤is the state-cost weighted matrix (symmetric, positive semidefinite), 𝑹𝑤is the input-cost weighted matrix (symmetric, positive definite),  𝑿𝑁𝑡is the state in the final instant of time, and 𝑸𝑁𝑡is the state-cost weighted matrix corresponding to the final instant of time (symmetric, positive semidefinite) [43]. Considering a multibody system with 𝑛generalized coordinates, 𝑚holonomic constraints, 𝑙nonholonomic constraints and 𝑛𝑢actuations, and defining 𝑁= 2 (𝑛−𝑚)−𝑙, 𝑸𝑤and 𝑸𝑁𝑡are 𝑁×𝑁matrices, and 𝑹𝑤is a 𝑛𝑢×𝑛𝑢matrix. Matrices 𝑸𝑤,𝑸𝑁𝑡and 𝑹𝑤can be regarded as tuning parameters in order to meet design requirements. These matrices are often diagonal: 𝑸𝑤=⎛⎜⎜⎜⎜⎝ 𝑄1,10⋯0 0𝑄2,2⋯0 ⋮0⋱ ⋮ 0 0 ⋯𝑄𝑁,𝑁 ⎞⎟⎟⎟⎟⎠ ,𝑹𝑤=⎛⎜⎜⎜⎜⎝ 𝑅1,10⋯0 0𝑅2,2⋯0 ⋮0⋱ ⋮ 0 0 ⋯𝑅𝑛𝑢,𝑛𝑢 ⎞⎟⎟⎟⎟⎠ ,(25) where the values of the diagonal elements 𝑄𝑖,𝑖 and 𝑅𝑗,𝑗 , with 𝑖= 1 … 𝑁and 𝑗= 1 … 𝑛𝑢, can be set by means of Bryson’s rule [44]: 𝑄𝑖,𝑖 =1 (𝐸𝑠 max)2, 𝑅𝑗,𝑗 =1 (𝐸𝑢 max)2.(26) In Eq. (26),𝐸𝑠 max and 𝐸𝑢 max are the maximum acceptable values of the errors in the states and inputs, respectively. The optimal control sequence is given by the feedback law:  𝒖∗ 𝑘= −𝑲𝑘 𝑿𝑘,(27) where the time-varying gain matrix 𝑲𝑘is computed following [45]: 𝑲𝑘=(𝑹𝑤+T 𝑘𝑷𝑘+1𝑘)−1 T 𝑘𝑷𝑘+1𝑘.(28) Mechanism and Machine Theory 203 (2024) 105795 7 A.G. Agúndez et al. In Eq. (28),𝑷𝑘is a symmetric and positive-definite matrix, solution of the discrete time Riccati equation, evolving backwards in time according to [45]: 𝑷𝑘=𝑸𝑤+T 𝑘𝑷𝑘+1𝑘−T 𝑘𝑷𝑘+1𝑘(𝑹𝑤+T 𝑘𝑷𝑘+1𝑘)−1 T 𝑘𝑷𝑘+1𝑘,(29) for 𝑘= 0,1,…, 𝑁𝑡− 1, with the terminal condition 𝑷𝑁𝑡=𝑸𝑁𝑡. Once the optimal sequence  𝒖∗is computed, the control inputs are, by virtue of Eq. (8): 𝒖=𝒖0+ 𝒖∗.(30) Note that, in this study, the Linear-Quadratic Regulator (LQR) has been chosen as feedback controller. Nevertheless, the linear system derived in this section, whether in continuous (Eq. (18)) or discrete form (Eq. (23)), can be used with other widely extended control techniques such as Model Predictive Control [33,34]. 3. Trajectory tracking of the ball-plate multibody system The methodology described in Section 2is applied to perform the trajectory tracking of the circular steady motion of the ballplate. A spatial multibody model of the system is described, the circular steady motion is computed and the linearized model along this reference solution is derived. The linear controllability is tested and forward dynamics simulations of the closed-loop system are shown. 3.1. Description of the ball-plate multibody model and nonholonomic constraints The ball-plate multibody model consists of two rigid bodies: the plate is denoted as body 2 and the sphere, which rolls without slipping on the plate, is body 3. The global reference frame is designated as body 1. The centres of mass 𝐺𝑗, with 𝑗= 2,3, correspond to the origins of the body reference frames. The origin of the global reference frame is located at 𝑂1. The contact point of the sphere with the plate is 𝐶. The sphere presents a mass 𝑚𝑒, radius 𝑅and moment of inertia, taken about any diameter, 𝐼. It is assumed that the ball and the plate cannot separate from each other. The angles that orientate the plate in space are denoted by 𝛼and 𝛽, and constitute the inputs of the multibody system. The 𝑛𝑢×1 vector of inputs 𝒖∈R𝑛𝑢, with 𝑛𝑢= 2, is: 𝒖=(𝛼 𝛽 )T.(31) The 𝑛× 1 vector of generalized coordinates 𝒙∈R𝑛, with 𝑛= 5, is given by: 𝒙=(𝑥𝐶𝑦𝐶𝜓 𝜙 𝜃 )T.(32) The Cartesian coordinates 𝑥𝐶and 𝑦𝐶are defined in the local reference frame of body 2 and locate the contact point 𝐶of the sphere on the plate. The angles 𝜓, 𝜙 and 𝜃describe the orientation of the ball relative to the plate: 𝜓is the yaw angle; 𝜙represents the lean angle; and 𝜃is the roll angle. The numbering of the bodies, the local reference frames, the set of generalized coordinates and the plate control inputs are shown in detail in Fig. 1. The orientation matrices of the body reference frames, expressed as a function of the elemental rotation matrices, are given by: 𝑹2=𝑹𝛼𝑹𝛽,𝑹3=𝑹2𝑹𝜓𝑹𝜙𝑹𝜃.(33) The position vectors of the contact point 𝐶and the centre of mass 𝐺3, expressed in the global reference frame, are given by: 𝒓𝐶=𝑹2 𝒓𝐶, 𝒓𝐶=(𝑥𝐶𝑦𝐶0)T,(34) 𝒓𝐺3=𝒓𝐶+𝑹2 𝒓𝐶𝐺3, 𝒓𝐶𝐺3=(0 0 𝑅)T.(35) The set of nonholonomic constraints is denoted as 𝑪𝑛ℎ (𝒙, 𝒙). These nonholonomic constraints result from the assumption that the sphere (body 3) rolls without slipping on the plate surface (body 2). Denoting the relative velocity at the contact point 𝐶by 𝒗𝐶32 , the condition for rolling without slipping can be written as:  𝒗𝐶32 =𝟎,(36) where  𝒗𝐶32 is expressed in the local reference frame of body 2. The 𝑙× 1 vector of nonholonomic constraints, with 𝑙= 2, is given by: 𝑪𝑛ℎ (𝒙, 𝒙)=(𝑣𝑋 𝐶32 𝑣𝑌 𝐶32 )T =𝟎2×1,(37) where 𝑣𝑋 𝐶32 and 𝑣𝑌 𝐶32 are the 𝑋- and 𝑌-components of  𝒗𝐶32 , which can be obtained as:  𝒗𝐶32 =𝑹T 2𝒗𝐶32 .(38) The relative velocity 𝒗𝐶32 , expressed in the global reference frame, can be computed as follows: 𝒗𝐶32 =𝒗𝐶31 −𝒗𝐶21 ,(39) Mechanism and Machine Theory 203 (2024) 105795 8 A.G. Agúndez et al. Fig. 1. Multibody model of the ball-plate system: numbering of the bodies, local reference frames, generalized coordinates of the multibody system and control input angles. Only a part of body 2 is represented. with 𝒗𝐶31 =𝒗𝐺31 +𝑹3( 𝝎3× 𝒓𝐺3𝐶),(40) 𝒗𝐶21 =𝒗𝐺21 +𝑹2( 𝝎2× 𝒓𝐶).(41) In Eqs. (40) and (41),𝒗𝐺21 and 𝒗𝐺31 are the velocities of the centres of mass of bodies 2and 3, respectively;  𝝎2and  𝝎3are the angular velocities of bodies 2and 3, respectively, expressed in the corresponding body frames; and  𝒓𝐺3𝐶= − 𝒓𝐶𝐺3, where  𝒓𝐶𝐺3was defined in Eq. (35). The expression of velocity  𝒗𝐶32 is given by:  𝒗𝐶32 =⎛⎜⎜⎜⎜⎝ 𝑣𝑋 𝐶32 𝑣𝑌 𝐶32 𝑣𝑍 𝐶32 ⎞⎟⎟⎟⎟⎠ =⎛⎜⎜⎜⎝ 𝑥𝐶−𝑅 𝜃cos (𝜓)cos (𝜙)−𝑅 𝜙sin (𝜓) 𝑦𝐶−𝑅 𝜃sin (𝜓)cos (𝜙)+𝑅 𝜙cos (𝜓) 0⎞⎟⎟⎟⎠ .(42) Note that, in Eq. (42), the 𝑍-component of  𝒗𝐶32 is already null and therefore not included in the nonholonomic constraints (37). Since the rolling sphere is described by using five coordinates (𝑛= 5), the number of degrees of freedom of the system is three: 𝑛𝑔=𝑛−𝑙= 3. As shown in Eq. (31), the number of control inputs of the system is 𝑛𝑢= 2. Since 𝑛𝑢< 𝑛𝑔, the system is underactuated. Denoting the nonholonomic constraints as 𝐶1 𝑛ℎ and 𝐶2 𝑛ℎ: 𝐶1 𝑛ℎ =𝑥𝐶−𝑅 𝜃cos (𝜓)cos (𝜙)−𝑅 𝜙sin (𝜓)= 0,(43) 𝐶2 𝑛ℎ =𝑦𝐶−𝑅 𝜃sin (𝜓)cos (𝜙)+𝑅 𝜙cos (𝜓)= 0,(44) the time derivatives of the angles  𝜙and  𝜃can be obtained by combining Eqs. (43) and (44) as follows: 𝐶1 𝑛ℎ sin (𝜓)−𝐶2 𝑛ℎ cos (𝜓)= 0 → 𝜙=𝑥𝐶sin (𝜓)−𝑦𝐶cos (𝜓) 𝑅,(45) 𝐶1 𝑛ℎ cos (𝜓)+𝐶2 𝑛ℎ sin (𝜓)= 0 → 𝜃=𝑥𝐶cos (𝜓)+𝑦𝐶sin (𝜓) 𝑅cos (𝜙).(46) The equations of motion of the ball-plate multibody system are given by the dynamic equations, which are obtained as, and the nonholonomic constraints (37). The following nonlinear index-2 DAE system is obtained: 𝑴(𝒙,𝒖) 𝒙+𝑫T(𝒙)𝜦=𝑸(𝒙, 𝒙,𝒖),(47) 𝑪𝑛ℎ (𝒙, 𝒙)=𝑫(𝒙) 𝒙=𝟎.(48) Note that, in the dynamic Eqs. (47), the 𝑛×𝑛mass matrix 𝑴(𝒙,𝒖)and the 𝑛× 1 vector of generalized forces 𝑸(𝒙, 𝒙,𝒖)depend on the vector of control inputs 𝒖. The symbolic expressions of the equations of motion are overly cumbersome to be shown in this manuscript, being available to the reader upon reasonable request. 3.2. Computation of the reference solution: circular steady motion of the ball The objective of this section is to develop a linear feedback controller to perform the circular trajectory tracking of the ball-plate system. To this end, the first step is to obtain the reference solution, which corresponds to the circular steady motion of the ball. Mechanism and Machine Theory 203 (2024) 105795 9 A.G. Agúndez et al. Fig. 5. Controllable system: forward dynamics simulation of the closed-loop system, considering the linearly controllable system (71), in both the linear (L) and nonlinear (NL) scenarios. Fig. 6. Controllable system: forward dynamics simulation of the closed-loop and open-loop systems, considering the linearly controllable system (71). Fig. 7. Controllable system: forward dynamics simulation considering the linearly controllable system (71): time evolution of the Cartesian coordinate 𝑦𝐶and the yaw angle 𝜓, in both the linear (L) and nonlinear (NL) scenarios. The results of the open-loop simulation are also shown. Mechanism and Machine Theory 203 (2024) 105795 16 A.G. Agúndez et al. Fig. 8. Time evolution of the angles 𝜙and 𝜃, in the open-loop and closed-loop scenarios. Fig. 9. Controllable system: time evolution of the plate orientation angles 𝛼and 𝛽, and fulfilment of the nonholonomic constraints. adjusting so that the ball, which is initially outside the circular trajectory, can follow the reference solution. Once on the reference solution, the angles 𝛼and 𝛽follow the evolution of 𝛼0and 𝛽0. The designed LQR controller satisfies the control objectives, while the nonholonomic constraints are also fulfilled, as shown in Fig. 9(b). The maximum values of the control input angles of the plate, 𝛼(𝑡)and 𝛽(𝑡), can be reduced by increasing the numerical values of the input-cost weighted matrix 𝑹𝑤of Eq. (25) and tuning the numerical values of the state-cost weighted matrix 𝑸𝑤and the state-cost weighted matrix corresponding to the final instant of time, 𝑸𝑁𝑡. Given the numerical values of matrices 𝑸𝑤,𝑸𝑁𝑡and 𝑹𝑤 detailed in Section 4.1, the following values are now considered: 𝑸′ 𝑤= 100 ⋅𝑸𝑤, 𝑸′ 𝑁𝑡= 100 ⋅𝑸′ 𝑤,(84) 𝑹′ 𝑤= 104⋅𝑹𝑤. Fig. 10(a) represents the time evolution of the plate orientation angles, 𝛼and 𝛽, considering the values detailed in Eq. (84) in the design of the LQR controller. It can be seen that the maximum values of 𝛼and 𝛽are significantly smaller than those of Fig. 9(a), and these angles reach their references 𝛼0and 𝛽0, respectively. The reduction in the initial kicks of the orientation angles results in a lower peak in velocity and acceleration at the initial instants of time, which would be advantageous for working with a real experimental setup of the ball-plate system and for ensuring that the ball does not separate from the plate during these initial instants of time before reaching the reference circular trajectory. The time derivatives of the control input angles, 𝛼 and  𝛽, are shown in Mechanism and Machine Theory 203 (2024) 105795 17 A.G. Agúndez et al. Fig. 10. Controllable system: time evolution of the plate orientation angles, 𝛼and 𝛽, and their time derivatives 𝛼 and  𝛽. Fig. 10(b). Despite this reduction in the maximum values of 𝛼and 𝛽, the LQR controller designed considering Eq. (84) satisfies the control objectives, with the contact point 𝐶of the sphere following the circular reference path. The linearization approach developed in Ref. [38], and applied to controlled multibody systems in Section 2.2, demonstrated excellent performance in designing linear feedback controllers for a nonholonomic multibody system as the ball-plate. The power of the approach to analytically compute the Jacobian matrix and the input matrix in Eq. (18) significantly facilitates the application of linear controllability criteria for linear-time varying (LTV) systems, as illustrated with the linearized model of the ball-plate system along the circular steady motion. The approach also enables the numerical computation of these matrices along any arbitrary reference solution. An acceptable computational efficiency was demonstrated in the performance of the circular trajectory tracking of the closed-loop nonlinear system, taking 11.5 s to perform a 12.56-second simulation, which would allow for real-time simulation with a ball-plate experimental setup. 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. The practical implementation of circular trajectory tracking can be addressed in future work by applying the methodology proposed in this study. Due to the method’s computational efficiency, real-time trajectory tracking would be feasible. To achieve this, a linear approximation of the system dynamics would need to be computed at each time step, enabling the calculation of the required control actions for the subsequent step. Updating the Jacobian matrix, , and the input matrix, , is straightforward since these matrices have been analytically derived as time-dependent functions. Moreover, in the experimental setup, measuring the lean and roll angles of the sphere is not required to track the time evolution of the generalized coordinates. This is because, in the coordinate partition of Eq. (64), these angles are treated as dependent coordinates during the linearization process. Instead, only the position of the contact point and the yaw angle, which are easier to measure, need to be recorded. The lean and roll angles can then be computed from these independent coordinates using the nonholonomic constraints. Although an experimental setup will not be used in this work, the behaviour of the controller is shown below when noise is introduced into the simulation of the closed-loop system, emulating the noise that would be encountered in measurements in a practical implementation. The same values of the LQR controller presented in Eq. (84) and the same initial conditions given by Eqs. (79) and (80) are considered. To this end, a random noise with normal distribution has been included in the states at each time step. Fig. 11(a) shows that despite the injection of noise into the simulation, the controller is able to stabilize the system, following the circular reference trajectory. Fig. 11(b) represents the time evolution of the Cartesian coordinates, 𝑥𝐶and 𝑦𝐶, the time evolution of the Cartesian coordinates including noise (which would be the measurements provided by the sensors in the experimental setup), denoted by 𝑥𝑛 𝐶and 𝑦𝑛 𝐶in Fig. 11(b), and their references 𝑥0 𝐶and 𝑦0 𝐶. The time derivatives of the Cartesian coordinates, 𝑥𝐶and 𝑦𝐶, are shown in Fig. 12(a); and Fig. 12(b) and Fig. 13(a) represent the time evolution of the yaw angle and its time derivative. Lastly, Fig. 13(b) depicts the time evolution of the control input angles, 𝛼and 𝛽, when the noise is injected in the simulation. It can be seen that, although not with the same accuracy as in the dynamic simulation without noise, in all cases the generalized coordinates, their time derivatives, and the control input angles follow their reference values. Finally, it is important to emphasize that the coordinate partition defined in Eq. (64) not only facilitates the derivation of a fully controllable reduced system of linear equations, but also offers a significant advantage over the coordinate partition described in Eq. (59). This advantage stems from the simplicity of measuring the contact point’s position compared to the orientation angles of the sphere. 5. Conclusions This paper addressed the optimal control of the ball-plate system. The trajectory tracking of the circular steady motion of the ball on the plate, for an arbitrarily large radius and potentially off centric with respect to the plate pivoting point, was achieved. Mechanism and Machine Theory 203 (2024) 105795 18 A.G. Agúndez et al. Fig. 11. Forward dynamics simulation of the closed-loop system, including a random noise with normal distribution in the simulation: trajectory of the contact point 𝐶and time evolution of the Cartesian coordinates, 𝑥𝐶and 𝑦𝐶. Fig. 12. Forward dynamics simulation of the closed-loop system, including a random noise with normal distribution in the simulation: time evolution of the time derivatives of the Cartesian coordinates, 𝑥𝐶and 𝑦𝐶, and the yaw angle, 𝜓. Previous works [11,24,31,32,36,37] devoted to the control of the ball-plate system employed a linear model around the trivial equilibrium position of the sphere on the horizontal plate. Consequently, the radius of the circular motion described by the sphere in circular trajectory tracking, as well as the control input angles of the plate, had to be kept very small to ensure the validity of the linear model. In contrast, the present work obtained, from a nonlinear multibody model, a linearized model of the ball-plate system along the circular steady motion, which is valid for a circular motion with an arbitrary radius, without any restriction in the values of the plate orientation angles. In contrast to Refs. [11,24,31,32,36,37], where a simplified model of the ball-plate system was presented, a spatial nonholonomic multibody model of the ball-plate, which includes the location of the contact point and the orientation angles of the sphere, was considered and the equations of motion were derived. A key contribution of the present work is the analytical computation of the circular steady motion by dynamic inversion, as well as the necessary control actions on the plate to achieve this reference solution. This enabled the analytical derivation of the linearized equations of motion as a function of the geometric and dynamic parameters of the ball-plate. To perform the linearization, the approach proposed in Ref. [38], which enables obtaining the reduced linearized equations of motion of constrained multibody systems [38,39], was extended to controlled multibody systems. The expression of the input matrix was derived and continuous and discrete-time representations of the linear system were shown. A linear time-varying (LTV) system was obtained, and a linear feedback controller was designed by using a discrete-time finite-horizon Linear-Quadratic Regulator (LQR). Mechanism and Machine Theory 203 (2024) 105795 19 A.G. Agúndez et al. Fig. 13. Forward dynamics simulation of the closed-loop system, including a random noise with normal distribution in the simulation: time evolution of the time derivative of the yaw angle, 𝜓, and the control input angles, 𝛼and 𝛽. The application of linear controllability criteria for LTV systems showed that the ball-plate system was not completely linearly controllable. To solve this, an alternative coordinate partition, which allowed for decoupling the linearized dynamic equations from the linearized nonholonomic constraints, was used. The linearized dynamic equations, expressed in terms of the Cartesian coordinates of the contact point and the yaw angle of the sphere, together with their time derivatives, constituted a controllable linear system. The importance of properly choosing a coordinate partition to achieve a controllable linear system was demonstrated with the ball-plate system. The forward dynamics simulation of the controllable system showed that the Cartesian coordinates of the contact point and the yaw angle of the sphere accurately followed the references, in both the linear and nonlinear scenarios. The plate orientation angles also achieved their reference values. The remaining orientation angles of the sphere (lean and roll angles) were computed from the nonholonomic constraints. With respect to previous works, where the radius of the circular motion described by the sphere and the control input angles of the plate had to be very small for the linear model to remain valid, the circular trajectory tracking presented in this work could be performed for large radii (𝜌= 0.5 m), without any restriction in the values of the plate orientation angles. An acceptable computational efficiency was demonstrated in the performance of the circular trajectory tracking of the closed-loop nonlinear system, which would allow for real-time simulation with a ball-plate experimental setup. In future work, the trajectory tracking of more complex trajectories in the ball-plate system could be carried out. Although this work has focused on circular trajectory tracking, the methodology presented in the present work is applicable to any trajectory, as it is based on a linearization along the reference solution, obtaining a linear time-varying system that can be updated at each time instant. Furthermore, validation of the trajectory tracking in a real ball-plate laboratory model could be addressed. Lastly, development of controllers using other control techniques could be considered. Despite the LQR has been chosen as feedback controller in this study, the linearized equations derived in Section 2.2, whether in continuous or discrete form, could be used with other widely extended control techniques, such as Model Predictive Control. Nonlinear control could also be explored, avoiding the need to calculate a linear approximation of the system dynamics at each instant of time, as is done in the present work. CRediT authorship contribution statement A.G. Agúndez: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. A. Saccon: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. D. García- Vallejo: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. E. Freire: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Mechanism and Machine Theory 203 (2024) 105795 20 A.G. Agúndez et al. Data availability Data will be made available on request. Acknowledgements This work was supported by Grant FPU18/05598 of the Spanish Ministry of Science, Innovation and Universities and project US-1380740 of Consejería de Universidad, Investigación e Innovación, Junta de Andalucía. Appendix A. Approximation of the time-varying discrete linear system Consider the discrete-time representation of the linear system (18), given by Eq. (23):  𝑿𝑘+1 = e𝑘𝛥𝑡  𝑿𝑘+(e𝑘𝛥𝑡 −𝑰)−1 𝑘𝑘 𝒖𝑘. The exponential of a matrix can be introduced by an infinite series of the same form that is used to define the scalar exponential: e𝑘𝛥𝑡 =𝑰+𝑘𝛥𝑡 +1 2! 2 𝑘𝛥𝑡2+⋯= ∞ ∑ 𝑗=0 1 𝑗!𝑗 𝑘𝛥𝑡𝑗.(85) The series of Eq. (85) is as nicely behaved as in the scalar case: it converges absolutely ∀𝑘∈R𝑁×𝑁and ∀𝛥𝑡 ∈R, and it can be differentiated or integrated term by term. Through the definition of 𝜳𝑘as follows [50]: 𝜳𝑘=(e𝑘𝛥𝑡 −𝑰)−1 𝑘= ∞ ∑ 𝑗=0 𝑗 𝑘𝛥𝑡𝑗+1 (𝑗+ 1)! ,(86) Eq. (23) can be written as in [50]:  𝑿𝑘+1 = e𝑘𝛥𝑡  𝑿𝑘+𝜳𝑘𝑘 𝒖𝑘.(87) In practice, e𝑘𝛥𝑡 and 𝜳𝑘are approximated by the corresponding truncated series. For example, truncating the series of Eq. (85) to 𝑗= 1, yields: e𝑘𝛥𝑡 ≃𝑰+𝑘𝛥𝑡. (88) Using the approximation of Eq. (88) in Eq. (23) leads to:  𝑿𝑘+1 =(𝑰+𝑘𝛥𝑡) 𝑿𝑘+𝛥𝑡𝑘 𝒖𝑘.(89) Appendix B. Derivation of the equations of motion of the ball-plate multibody system This appendix details the derivation of the equations of motion of the ball-plate multibody system as explained in Schiehlen [51, 52], where the reader can find further details. Consider the ball-plate multibody model described in Section 3.1, which presents two bodies and 𝑛= 5 generalized coordinates. First, the velocities and accelerations are expressed as a function of the generalized coordinates of the multibody system and their time derivatives. Consider body 𝑖, with 𝑖= 2,3, in the ball-plate system. The velocity of the centre of mass 𝐺𝑖,𝒗𝐺𝑖, can be written as: 𝒗𝐺𝑖=d𝒓𝐺𝑖 d𝑡= 𝜕𝒓𝐺𝑖 𝜕𝒙 𝒙,(90) where 𝒓𝐺𝑖is the position vector of the centre of mass 𝐺𝑖. In Eq. (90), it is assumed that 𝒓𝐺𝑖does not depend explicitly on time. Eq. (90) can be written as a function of the generalized coordinates as: 𝒗𝐺𝑖=𝑯𝑖 𝒙,(91) where 𝑯𝑖is the following 3 × 𝑛Jacobian matrix: 𝑯𝑖= 𝜕𝒓𝐺𝑖 𝜕𝒙.(92) Likewise, differentiating Eq. (91) with respect to time, the expression for the acceleration of the centre of mass is obtained: 𝒂𝐺𝑖=d𝒗𝐺𝑖 d𝑡=𝑯𝑖 𝒙+ 𝑯𝑖 𝒙.(93) The time derivative of the Jacobian matrix  𝑯𝑖is computed as follows:  𝑯𝑖=𝜕𝑯𝑖 𝜕𝒙 𝒙=𝜕(𝑯𝑖 𝒙) 𝜕𝒙= 𝜕𝒗𝐺𝑖 𝜕𝒙=𝒉𝑖.(94) Mechanism and Machine Theory 203 (2024) 105795 21 A.G. Agúndez et al. Therefore, the acceleration 𝒂𝐺𝑖can be expressed in terms of the generalized coordinates and their time derivatives as: 𝒂𝐺𝑖=d𝒗𝐺𝑖 d𝑡=𝑯𝑖 𝒙+𝒉𝑖 𝒙.(95) Note that 𝒉𝑖is a 3 × 𝑛matrix. Regarding the angular acceleration  𝜶𝑖, the same steps in the computation of the linear acceleration 𝒂𝐺𝑖are followed, leading to:  𝜶𝑖=d 𝝎𝑖 d𝑡= 𝑮𝑖 𝒙+  𝑮𝑖 𝒙,(96) where  𝝎𝑖is the angular velocity of body 𝑖, expressed in the body reference frame. The angular velocity  𝝎𝑖can be written as:  𝝎𝑖= 𝑮𝑖 𝒙,(97) with matrix  𝑮𝑖:  𝑮𝑖=𝜕 𝝎𝑖 𝜕𝒙.(98) It is assumed that  𝝎𝑖does not dependent explicitly on time. In Eq. (96), the time derivative of   𝑮𝑖can be computed as:   𝑮𝑖=𝜕 𝑮𝑖 𝜕𝒙 𝒙=𝜕( 𝑮𝑖 𝒙) 𝜕𝒙=𝜕 𝝎𝑖 𝜕𝒙= 𝒈𝑖.(99) Therefore, the expression of the angular acceleration  𝜶𝑖, as a function of the generalized coordinates and their time derivatives, is given by:  𝜶𝑖=d 𝝎𝑖 d𝑡= 𝑮𝑖 𝒙+ 𝒈𝑖 𝒙.(100) As 𝑯𝑖and 𝒉𝑖, 𝑮𝑖and  𝒈𝑖are 3 × 𝑛matrices. The matrices  𝑮𝑖and  𝒈𝑖have been computed in the body reference frame. The linear and angular accelerations 𝒂𝐺𝑖and  𝜶𝑖, with 𝑖= 2,3, and expressed in the global and local reference frames, respectively, can be assembled by using Eqs. (97) and (100) as follows:  𝒂=⎛⎜⎜⎜⎜⎝ 𝒂𝐺2 𝒂𝐺3  𝜶2  𝜶3 ⎞⎟⎟⎟⎟⎠ =⎛⎜⎜⎜⎜⎝ 𝑯2 𝑯3  𝑮2  𝑮3 ⎞⎟⎟⎟⎟⎠  𝒙+⎛⎜⎜⎜⎜⎝ 𝒉2 𝒉3  𝒈2  𝒈3 ⎞⎟⎟⎟⎟⎠  𝒙.(101) Defining the matrices 𝑳and 𝒍as: 𝑳=⎛⎜⎜⎜⎜⎝ 𝑯2 𝑯3  𝑮2  𝑮3 ⎞⎟⎟⎟⎟⎠ ,𝒍=⎛⎜⎜⎜⎜⎝ 𝒉2 𝒉3  𝒈2  𝒈3 ⎞⎟⎟⎟⎟⎠ ,(102) Eq. (101) is rewritten as:  𝒂=𝑳 𝒙+𝒍 𝒙.(103) Consider now the Newton-Euler equations of a rigid body: (𝒎𝑖𝟎3 𝟎3 𝑰𝐺𝑖)( 𝒂𝐺𝑖  𝜶𝑖)=(𝑭𝑖  𝑴𝑖)+⎛⎜⎜⎝ 𝟎3×1 − 𝝎𝑖×( 𝑰𝐺𝑖 𝝎𝑖)⎞⎟⎟⎠ ,(104) where 𝒎𝑖=𝑚𝑖𝑰3(with 𝑖= 2,3in the ball-plate), 𝑚𝑖is the mass of body 𝑖,𝑰3is the identity matrix of size 3,𝟎3is the null square matrix of size 3,𝟎3×1 is a 3 ×1 null matrix, and  𝑰𝐺𝑖is the inertia tensor of body 𝑖, observed in the body reference frame. In Eq. (104), 𝑭𝑖is the resultant force acting on body 𝑖, expressed in the global reference frame, and  𝑴𝑖is the resultant torque, expressed in the body reference frame. The Newton-Euler Eqs. (104) of bodies 2(plate) and 3(ball) can be assembled, leading to:  𝑴 𝒂= 𝑸𝑎+ 𝑸𝑣,(105) where  𝒂was defined in Eq. (101), and  𝑴, 𝑸𝑎and  𝑸𝑣are given by:  𝑴=⎛⎜⎜⎜⎜⎝ 𝒎2𝟎3𝟎3𝟎3 𝟎3𝒎3𝟎3𝟎3 𝟎3𝟎3 𝑰𝐺2𝟎3 𝟎3𝟎3𝟎3 𝑰𝐺3 ⎞⎟⎟⎟⎟⎠ ,(106) Mechanism and Machine Theory 203 (2024) 105795 22 A.G. Agúndez et al.  𝑸𝑎=⎛⎜⎜⎜⎜⎝ 𝑭2 𝑭3  𝑴2  𝑴3 ⎞⎟⎟⎟⎟⎠ , 𝑸𝑣=⎛⎜⎜⎜⎜⎜⎜⎝ 𝟎3×1 𝟎3×1 − 𝝎2×( 𝑰𝐺2 𝝎2) − 𝝎3×( 𝑰𝐺3 𝝎3) ⎞⎟⎟⎟⎟⎟⎟⎠ .(107) Substituting Eq. (103) in Eq. (105), one obtains:  𝑴(𝑳 𝒙+𝒍 𝒙)= 𝑸𝑎+ 𝑸𝑣.(108) Eq. (108) can be premultiplied by 𝑳T(matrix 𝑳was defined in Eq. (102)). Rearranging the terms, the dynamic equations of the ball-plate multibody system can be written as: 𝑴 𝒙=𝑸𝑎+𝑸𝑣,(109) where 𝑴=𝑳T 𝑴𝑳,(110) 𝑸𝑎=𝑳T 𝑸𝑎,(111) 𝑸𝑣=𝑳T( 𝑸𝑣− 𝑴𝒍  𝒙).(112) In Eq. (109), the following terms can be identified: -𝑴 𝒙is the generalized inertia term, being 𝑴the 𝑛×𝑛mass matrix. -𝑸𝑣is the generalized quadratic-velocity forces vector, associated with the inertia forces that are quadratic with respect to the system velocities (centrifugal and Coriolis forces). -𝑸𝑎is the generalized forces vector, which can be separated into two terms: 𝑸𝑎=𝑸apl +𝑸reac.(113) In Eq. (113),𝑸apl is associated with the external forces and torques acting on the multibody system, while 𝑸reac arises from the reaction forces associated with the constraints that the multibody system satisfies during its motion. The vector 𝑸reac can be obtained by using the method of Lagrange multipliers. Considering a multibody system with 𝑚holonomic and 𝑙nonoholonomic constraints, 𝑸reac can be written as: 𝑸reac = −𝑫T(𝒙)𝜦,(114) where 𝜦∈R𝑚+𝑙is the (𝑚+𝑙)×1 vector of Lagrange multipliers, and 𝑫(𝒙)was defined in Eq. (4). The Lagrange multipliers represent the reaction forces associated with the constraints, with opposite sign. In the particular case of the ball-plate multibody model, 𝑚= 0 and 𝑙= 2. Using Eq. (114) in Eq. (109), the dynamic equilibrium equations of the ball-plate are given by: 𝑴(𝒙,𝒖) 𝒙+𝑫T(𝒙)𝜦=𝑸(𝒙, 𝒙,𝒖).(115) Note that, in the dynamic Eqs. (115) of the ball-plate, 𝑸=𝑸𝑣+𝑸𝑔, with 𝑸𝑔being the generalized gravity force vector. Furthermore, the 𝑛×𝑛mass matrix 𝑴(𝒙,𝒖)and the 𝑛× 1 vector of generalized forces 𝑸(𝒙, 𝒙,𝒖)depend on the vector of control inputs 𝒖. The dynamic Eqs. (115) must be solved together with the constraints of the multibody system, leading to the equations of motion (47)–(48) presented in Section 3.1. 𝑴(𝒙,𝒖) 𝒙+𝑫T(𝒙)𝜦=𝑸(𝒙, 𝒙,𝒖), 𝑪𝑛ℎ (𝒙, 𝒙)=𝑫(𝒙) 𝒙=𝟎. Appendix C. Coefficients of the linearized equations of motion The expressions of the coefficients 𝐴𝑗in Eq. (61), with 𝑗= 1 … 10, are given by: 𝐴1= −𝛺𝜔 sin (𝜙0), 𝐴2=𝜔sin (𝜙0)−𝛺 cos (𝜙0), 𝐴3=𝛺cos (𝜙0), 𝐴4=𝜔cos (𝜙0), 𝐴5=𝛺𝜔, 𝐴6=𝛺sin (𝜙0)−𝜔 cos (𝜙0),(116) 𝐴7= −𝛺𝑅 cos (𝜙0), 𝐴8=𝛺𝑅 sin (𝜙0), 𝐴9=𝑅, 𝐴10 = −𝑅cos (𝜙0). Mechanism and Machine Theory 203 (2024) 105795 23 A.G. Agúndez et al. The coefficients 𝐵𝑟(𝑡)in Eq. (62), with 𝑟= 1 … 6, are the following periodic time-dependent functions: 𝐵1(𝑡)=𝑚𝑒𝑔𝑅 tan (𝜙0) 𝐼+𝑚𝑒𝑅2(cos (𝜔𝑡)cos (𝛼0)− sin (𝜔𝑡)sin (𝛼0)sin (𝛽0)), 𝐵2(𝑡)=𝑚𝑒𝑔𝑅 tan (𝜙0) 𝐼+𝑚𝑒𝑅2sin (𝜔𝑡)cos (𝛼0)cos (𝛽0), 𝐵3(𝑡)= − 𝑚𝑒𝑔𝑅 𝐼+𝑚𝑒𝑅2(sin (𝜔𝑡)cos (𝛼0)+ cos (𝜔𝑡)sin (𝛼0)sin (𝛽0)),(117) 𝐵4(𝑡)=𝑚𝑒𝑔𝑅 𝐼+𝑚𝑒𝑅2cos (𝜔𝑡)cos (𝛼0)cos (𝛽0), 𝐵5(𝑡)= − 𝑚𝑒𝑔𝑅 (𝐼+𝑚𝑒𝑅2)cos (𝜙0)(cos (𝜔𝑡)cos (𝛼0)− sin (𝜔𝑡)sin (𝛼0)sin (𝛽0)), 𝐵6(𝑡)= − 𝑚𝑒𝑔𝑅 (𝐼+𝑚𝑒𝑅2)cos (𝜙0)sin (𝜔𝑡)cos (𝛼0)cos (𝛽0). The coefficients 𝑎𝑗in Eq. (65), with 𝑗= 1 … 5, are given by: 𝑎1= −𝛺2, 𝑎2= −𝑎5(𝛺−𝜔√1 − (𝜌𝜔 𝛺𝑅 )2), 𝑎3=𝜌𝜔 𝑅, 𝑎4=1 𝑅,(118) 𝑎5=𝛺 𝜌𝜔 , 𝑎6=𝑎5𝛺𝑅√1 − (𝜌𝜔 𝛺𝑅 )2 . The coefficients 𝑏𝑟(𝑡)in Eq. (66), with 𝑟= 1 … 5, are periodic time-dependent functions. Defining the functions 𝑓1(𝑡)and 𝑓2(𝑡) as follows: 𝑓1(𝑡)=√ √ √ √1 − 𝜌2𝜔4(𝐼+𝑚𝑒𝑅2)2 𝑚2 𝑒𝑔2𝑅4sin2(𝜔𝑡),(119) 𝑓2(𝑡)=√ √ √ √ √1 + 𝜌2𝜔4(𝐼+𝑚𝑒𝑅2)2 𝜌2𝜔4(𝐼+𝑚𝑒𝑅2)2sin2(𝜔𝑡)−𝑚2 𝑒𝑔2𝑅4 cos2(𝜔𝑡),(120) the coefficients 𝑏𝑟(𝑡)are given by: 𝑏1(𝑡)=𝜌2𝜔4(𝐼+𝑚𝑒𝑅2) 𝑚𝑒𝑔𝑅2𝑓1(𝑡)sin (𝜔𝑡)cos (𝜔𝑡), 𝑏2(𝑡)=𝑚𝑒𝑔𝑅2 𝐼+𝑚𝑒𝑅2𝑓1(𝑡)𝑓2(𝑡), 𝑏3(𝑡)= − 𝑚𝑒𝑔𝑅2 𝐼+𝑚𝑒𝑅2𝑓1(𝑡),(121) 𝑏4(𝑡)=𝑚𝑒𝑔𝑅 sin (𝜙0) (𝐼+𝑚𝑒𝑅2)cos (𝜙0)𝑓1(𝑡)cos (𝜔𝑡), 𝑏5(𝑡)=𝑚𝑒𝑔𝑅 sin (𝜙0)sin (𝜔𝑡) (𝐼+𝑚𝑒𝑅2)cos (𝜙0)𝑓1(𝑡)𝑓2(𝑡). References [1] F. Ruggiero, V. Lippiello, B. Siciliano, Nonprehensile dynamic manipulation: A survey, IEEE Robot. Autom. Lett. 3 (3) (2018) 1711–1718. [2] D. Serra, Robot control for nonprehensile dynamic manipulation tasks, in: International Conference on Informatics in Control, Automation and Robotics Doctoral Consortium, 2016, pp. 3–12. [3] B. Siciliano, F. Ruggiero, Robot Dynamic Manipulation: Perception of Deformable Objects and Nonprehensile Manipulation Control, vol. 144, Springer Nature, 2022. [4] A. Gutierrez-Giles, A.C. Satici, A. Donaire, F. Ruggiero, V. Lippiello, B. Siciliano, Nonholonomic rolling nonprehensile manipulation primitive, Robot Dyn. Manip. Percept. Deform. Objects Nonprehensile Manip. Control. (2022) 159–205. [5] J.Z. Woodruff, S. Ren, K.M. Lynch, Motion planning and feedback control of rolling bodies, IEEE Access 8 (2020) 31780–31791. [6] K.M. Lynch, N. Shiroma, H. Arai, K. Tanie, The roles of shape and motion in dynamic manipulation: The butterfly example, in: Proceedings. 1998 IEEE International Conference on Robotics and Automation (Cat. No. 98CH36146), vol. 3, IEEE, 1998, pp. 1958–1963. [7] M. Cefalo, L. Lanari, G. Oriolo, Energy-based control of the butterfly robot, IFAC Proc. Vol. 39 (15) (2006) 1–6. [8] D.O. Morales, P. La Hera, S.U. Réhman, Generating periodic motions for the butterfly robot, in: 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, IEEE, 2013, pp. 2527–2532. [9] M. Surov, A. Shiriaev, L. Freidovich, S. Gusev, L. Paramonov, Case study in non-prehensile manipulation: planning and orbital stabilization of one-directional rollings for the ‘‘butterfly’’ robot, in: 2015 IEEE International Conference on Robotics and Automation, ICRA, IEEE, 2015, pp. 1484–1489. Mechanism and Machine Theory 203 (2024) 105795 24 A.G. Agúndez et al. [10] R. Mukherjee, M.A. Minor, J.T. Pukrushpan, Motion planning for a spherical mobile robot: Revisiting the classical ball-plate problem, J. Dyn. Syst. Meas. Control 124 (4) (2002) 502–511. [11] X. Fan, N. Zhang, S. Teng, Trajectory planning and tracking of ball and plate system using hierarchical fuzzy control scheme, Fuzzy Sets and Systems 144 (2) (2004) 297–312. [12] L. Morales, M. Gordón, O. Camacho, A. Rosales, D. Pozo, A comparative analysis among different controllers applied to the experimental ball and plate system, in: 2017 International Conference on Information Systems and Computer Science, INCISCOS, IEEE, 2017, pp. 108–114. [13] H. Bang, Y.S. Lee, Implementation of a ball and plate control system using sliding mode control, IEEE access 6 (2018) 32401–32408. [14] J.Z. Woodruff, K.M. Lynch, Planning and control for dynamic, nonprehensile, and hybrid manipulation tasks, in: 2017 IEEE International Conference on Robotics and Automation, ICRA, IEEE, 2017, pp. 4066–4073. [15] W.H. Huang, G.F. Holden, Nonprehensile palmar manipulation with a mobile robot, in: Proceedings 2001 IEEE/RSJ International Conference on Intelligent Robots and Systems. Expanding the Societal Role of Robotics in the the Next Millennium (Cat. No. 01CH37180), vol. 1, IEEE, 2001, pp. 114–119. [16] C. Fan, S. Shirafuji, J. Ota, Modal planning for cooperative non-prehensile manipulation by mobile robots, Appl. Sci. 9 (3) (2019) 462. [17] M. Costanzo, M. De Simone, S. Federico, C. Natale, Non-prehensile manipulation actions and visual 6D pose estimation for fruit grasping based on tactile sensing, Robotics 12 (4) (2023) 92. [18] J. Shi, J.Z. Woodruff, P.B. Umbanhowar, K.M. Lynch, Dynamic in-hand sliding manipulation, IEEE Trans. Robot. 33 (4) (2017) 778–795. [19] F. Bertoncelli, F. Ruggiero, L. Sabattini, Linear time-varying MPC for nonprehensile object manipulation with a nonholonomic mobile robot, in: 2020 IEEE International Conference on Robotics and Automation, ICRA, IEEE, 2020, pp. 11032–11038. [20] H. Ebel, P. Eberhard, Non-prehensile cooperative object transportation with omnidirectional mobile robots: Organization, control, simulation, and experimentation, in: 2021 International Symposium on Multi-Robot and Multi-Agent Systems, MRS, IEEE, 2021, pp. 1–10. [21] Y. Farid, B. Siciliano, F. Ruggiero, Review and descriptive investigation of the connection between bipedal locomotion and non-prehensile manipulation, Annu. Rev. Control 53 (2022) 51–69. [22] J.Z. Woodruff, K.M. Lynch, Robotic contact juggling, 2021, arXiv preprint arXiv:2102.10421. [23] C.M. Pappalardo, D. Guida, On the dynamics and control of underactuated nonholonomic mechanical systems and applications to mobile robots, Arch. Appl. Mech. 89 (2019) 669–698. [24] A. Knuplez, A. Chowdhury, R. Svecko, Modeling and control design for the ball and plate system, in: IEEE International Conference on Industrial Technology, 2003, vol. 2, IEEE, 2003, pp. 1064–1067. [25] G. Oriolo, M. Vendittelli, A framework for the stabilization of general nonholonomic systems with an application to the plate-ball mechanism, IEEE Trans. Robot. 21 (2) (2005) 162–175. [26] H. Date, M. Sampei, M. Ishikawa, M. Koga, Simultaneous control of position and orientation for ball-plate manipulation problem based on time-state control form, IEEE Trans. Robot. Autom. 20 (3) (2004) 465–480. [27] M.R. Azizi, J. Keighobadi, Point stabilization of nonholonomic spherical mobile robot using nonlinear model predictive control, Robot. Auton. Syst. 98 (2017) 347–359. [28] D. Serra, J. Ferguson, F. Ruggiero, A. Siniscalco, A. Petit, V. Lippiello, B. Siciliano, On the experiments about the nonprehensile reconfiguration of a rolling sphere on a plate, in: 2018 26th Mediterranean Conference on Control and Automation, MED, IEEE, 2018, pp. 13–20. [29] C.-C. Cheng, C.-H. Tsai, Visual servo control for balancing a ball-plate system, Int. J. Mech. Eng. Robot. Res. 5 (1) (2016) 28. [30] D. Casagrande, A. Astolfi, T. Parisini, Switching-driving Lyapunov function and the stabilization of the ball-and-plate system, IEEE Trans. Autom. Control 54 (8) (2009) 1881–1886. [31] C.C. Ker, C.E. Lin, R.T. Wang, Tracking and balance control of ball and plate system, J. Chin. Inst. Eng. 30 (3) (2007) 459–470. [32] A. Kassem, H. Haddad, C. Albitar, Commparison between different methods of control of ball and plate system with 6dof stewart platform, IFAC-PapersOnLine 48 (11) (2015) 47–52. [33] P. Falcone, M. Tufo, F. Borrelli, J. Asgari, H.E. Tseng, A linear time varying model predictive control approach to the integrated vehicle dynamics control problem in autonomous systems, in: 2007 46th IEEE Conference on Decision and Control, IEEE, 2007, pp. 2980–2985. [34] E.F. Camacho, C.B. Alba, Model predictive control, Springer science & business media, 2013. [35] J. Fan, M. Han, Nonliear model predictive control of ball-plate system based on gaussian particle swarm optimization, in: 2012 IEEE Congress on Evolutionary Computation, IEEE, 2012, pp. 1–6. [36] M. Oravec, A. Jadlovská, Model predictive control of a ball and plate laboratory model, in: 2015 IEEE 13th International Symposium on Applied Machine Intelligence and Informatics, SAMI, IEEE, 2015, pp. 165–170. [37] H. Bang, Y.S. Lee, Embedded model predictive control for enhancing tracking performance of a ball-and-plate system, IEEE Access 7 (2019) 39652–39659. [38] A. García-Agúndez, D. García-Vallejo, E. Freire, Linearization approaches for general multibody systems validated through stability analysis of a benchmark bicycle model, Nonlinear Dynam. 103 (2021) 557–580. [39] A. Agúndez, D. García-Vallejo, E. Freire, A. Mikkola, A reduced and linearized high fidelity waveboard multibody model for stability analysis, J. Comput. Nonlinear Dyn. 17 (5) (2022) 051010. [40] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems, ASME J. Mech. Des. 104 (1) (1982) 247–255. [41] E. Pennestrı, P. Valentini, Coordinate reduction strategies in multibody dynamics: A review, in: Proceedings of the Conference on Multibody System Dynamics, 2007. [42] A.A. Shabana, Dynamics of multibody systems, Cambridge University Press, 2020. [43] S.L. Brunton, J.N. Kutz, Data-driven science and engineering: Machine learning, dynamical systems, and control, Cambridge University Press, 2022. [44] A.E. Bryson, Applied optimal control: optimization, estimation and control, CRC Press, 1975. [45] A. Shaiju, I.R. Petersen, Formulas for discrete time LQR, LQG, LEQG and minimax LQG optimal control problems, IFAC Proc. Vol. 41 (2) (2008) 8773–8778. [46] W.L. Brogan, Modern Control Theory, Prentice-Hall, Inc., 1991. [47] M. Baake, U. Schlaegel, The peano-baker series, Proc. Steklov Inst. Math. 275 (2011) 155–159. [48] L.M. Silverman, H. Meadows, Controllability and observability in time-variable linear systems, SIAM J. Control 5 (1) (1967) 64–73. [49] H. Leiva, S. Siegmund, A necessary algebraic condition for controllability and observability of linear time-varying systems, IEEE Trans. Autom. Control 48 (12) (2003) 2229–2232. [50] L. Pyrhönen, S. Jaiswal, A. Garcia-Agundez, D. García Vallejo, A. Mikkola, Linearization-based state-transition model for the discrete extended Kalman filter applied to multibody simulations, Multibody Syst. Dyn. 57 (1) (2023) 55–72. [51] W. Schiehlen, Multibody system dynamics: roots and perspectives, Multibody Syst. Dyn. 1 (1997) 149–188. [52] J.L. Escalona, A.M. Recuero, A bicycle model for education in multibody dynamics and real-time interactive simulation, Multibody Syst. Dyn. 27 (2012) 383–402. Mechanism and Machine Theory 203 (2024) 105795 25