An efficient and accurate linearization approach for hydraulically actuated multibody systems with holonomic and nonholonomic constraints

Hydraulics is often used to actuate mechanisms in the applications of heavy machinery. In this work, a linearization approach for hydraulically driven multibody systems is presented. The approach allows linearizing the equations of motion of general multibody systems with holonomic and nonholonomic constraints, augmented with the hydraulic equations of the hydraulic subsystem. The derivation of this linearization approach is of interest in many applications, such as the performance of linear stability analyses. The procedure is tested with a three-dimensional multibody model of a hydraulically actuated four-bar mechanism. The validation of the approach is performed by means of the forward dynamics simulation of the linear and nonlinear systems. The results show the power of the approach, obtaining the linearized equations of motion around the equilibrium position of the four-bar mechanism multibody model in terms of the mechanical and hydraulic parameters. A comparison of the proposed procedure with a conventional counterpart approach is included, demonstrating the great accuracy and computational efficiency of the approach developed in this work.


Introduction
Hydraulically actuated systems are widely used in heavy machinery, being present in a large variety of industries and research applications. Some examples of hydraulic systems used in daily life are aircrafts (activation and motion of landing gears or flaps), hydraulic lifts, mobile machines (cranes or excavators), hydraulic power steering, braking systems of vehicles, hydraulic jacks or shock absorbers. The modeling and simulation of hydraulic systems enables evaluating their behavior for a wide range of operating conditions. Positioning accuracy of the systems or identification of peak stresses are examples of important issues that can be addressed by means of simulation of hydraulic systems.
The use of the lumped fluid theory [1] allows 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 [2][3][4][5] and modeling complex hydraulic circuits for real-time simulation [6].
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 [7]. Nevertheless, including the dynamics of the hydraulic actuator is necessary in many applications [8,9]. From the standpoint of 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 [10,11]. Docquier et al. [12] presented a multibody model of a modern car equipped with a novel suspension system using a monolithic scheme. Ylinen et al. [13] 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. [14]. Rahikainen et al. [15] combined the index-3 semirecursive formulation [16] 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 [2]. Lastly, within the monolithic framework, the friction modeling in the hydraulic cylinder, which plays an important role in the accuracy of the simulations and can introduce numerical stiffness [17], was thoroughly studied by Jaiswal et al. [5]. Four friction modeling 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) [18,19] or resorting to a different software for each problem (co-simulation) [20][21][22]. Some relevant aspects have been addressed in the literature, such as co-simulation configuration [23], energy-based coupling error minimization [23,24] or the multirate cosimulation [25].
Different possibilities exist to carry out the linearization of the equations of motion of multibody systems, depending on the form of these nonlinear equations. Some approaches are based on the direct linearization of the Differential-Algebraic system of Equations (DAE) [26][27][28], while others resort to a coordinate partitioning to reduce the nonlinear DAE system to a nonlinear system of Ordinary Differential Equations (ODEs) [27,29,30]. Agúndez et al. [30] proposed a linearization approach, consisting in the linearization of the index-2 DAE system and then the reduction to a linear ODE system. The procedure, which showed an excellent accuracy and computational efficiency with complex multibody systems as the bicycle benchmark of Meijaard et al. [31], achieves the maximum reduction of the linearized equations of motion of constrained multibody systems. Nevertheless, to the best knowledge of the authors, there is no procedure for systematically obtaining the linearized equations of motion of hydraulically actuated multibody systems, which is required in several applications.
First, an important application is the building of state and input estimators, like Kalman filters, existing different works in the framework of hydraulically driven multibody systems. Khadim et al. [32] proposed a parameter estimation algorithm, consisting in the combination of the augmented discrete extended Kalman filter (ADEKF) with a curve-fitting method, and Jaiswal et al. [33] presented a state estimator based on an indirect Kalman filter. Secondly, another important application is the performance of linear stability analyses and the design of linear feedback controllers. Hydraulic steering systems play a key role in keeping the directional stability and tracking the steerhandling capability of articulated steering vehicles (ASVs). After some simplified models without considering the dynamics characteristics [34,35], where the hydraulic steering system was modeled as a torsion spring, the full-hydraulic steering control unit of the steering system model was modeled as a directional control valve by Pazooki et al. [36]. Several works are devoted to improving the yaw stability of ASVs, being some techniques increasing the damping at the articulation joint [37], introducing leakage across the cylinders [34], or other active strategies as differential braking [38], torque vectoring [35] and active steering [39]. Gao et al. [40] analyzed the stability of an ASV and designed a stability controller to avoid the oscillatory yaw motion, with the application of the optimal control theory. An optimal tuned cascade control strategy, based on feedback linearization, is proposed by Nedić et al. [41] to carry out the reference trajectory of a 6-DOF parallel robot platform. An optimal design of the parameters of the cascade load force controller was effectively performed. In these applications, a linearized version of the equations of motion is required.
The main objective of this work is to present an approach for systematically linearizing the equations of motion of hydraulically actuated multibody systems with holonomic and nonholonomic constraints. The procedure corresponds to an extension of the approach presented and validated by Agúndez et al. [30] for mechanical multibody systems. The main advantages of the proposed procedure are its computational efficiency and accuracy, which are compared with a linearization counterpart approach; the maximum reduction of the linearized equations of motion, which allows the elimination of spurious null eigenvalues in the performance of linear stability analyses; and the capacity of the proposed approach to generate the exact linearized equations of motion in terms of the mechanical and hydraulic parameters of the multibody system under study.
The paper is structured as follows. Following the Introduction, Sect. 2 summarizes the main aspects of the hydraulic modeling, using the lumped fluid method, and presents the nonlinear equations of motion of hydraulically actuated multibody systems. Next, Sect. 3 develops in detail the linearization approach and derives the resulting Jacobian matrix. In Sect. 4, the use of the procedure is illustrated with a hydraulically actuated three-dimensional four-bar mechanism model and validated by means of the forward dynamics simulation of the linear and nonlinear systems. Finally, Sect. 5 summarizes the main conclusions drawn from the present work.

Formulation of the problem
In this section, the fundamental modeling aspects of a hydraulic linear actuator system are presented. More-over, the nonlinear equations of motion of hydraulically actuated multibody systems with holonomic and nonholonomic constraints are shown below.

Modeling of the hydraulic system
In this work, the lumped fluid theory is used to model the hydraulic system. The hydraulic circuit is divided into discrete volumes of uniformly distributed pressure. Considering a discrete control volume V , the evolution of the pressure p is given by the following first-order differential equation: where Q k represents the incoming (positive value) or outcoming (negative value) flows of the control volume, n f is the number of hydraulic flows going in or out of the volume, dV dt is the volume change term that usually represents the piston movement inside the cylinder and B e is the effective bulk modulus of the hydraulic volume, given by: In Eq. (2), B oil is the oil bulk modulus, n s is the number of subvolumes V k forming the volume V and B k is the bulk modulus of the subvolume V k . The semi-empirical modeling method [42] is used to describe the valves in the hydraulic circuit. First, the volume flow rate through a directional control valve Q d is: where C v d is the semi-empirical flow rate coefficient of the directional control valve; p is the pressure difference, which presents the same direction as the volume flow rate; and U is a normalized spool displacement, which determines the spool position in the control valve. The time evolution of U is given by the following first-order differential equation: where U ref is the reference normalized spool displacement and τ is the time constant of the valve.
In the case of a throttle valve, the flow rate Q t is expressed as: In Eq. (5), A t is the area of the throttle valve, C d is the coefficient of discharge and ρ oil is the oil density. The expressions of the volume flow rates in Eqs. (3) and (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. (3) and (5) present a vertical tangent for p → 0. Therefore, the laminar regime of the volume flow rates is considered for pressure drops lower than a predefined p lim , following a linear relation with the pressure difference. In this study, it is assumed that p lim = 2 bar. The volume flow rate through a directional control valve or a throttle valve can be written as: where γ 1 and γ 2 are a function of the valve parameters. To ensure the continuity of the laminar and turbulent regimes at p lim , corresponding to a volume flow rate Q lim (see Fig. 1), the relation 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 F cyl , is computed as follows:  where A 1 and A 2 are the piston and piston-rod side areas of the cylinder, respectively; p 1 is the pressure of the piston side chamber; p 2 is 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 2 shows a free-body diagram of the hydraulic cylinder. Following Ref. [5], the Brown and McPhee model [43] is numerically one of the most efficient approaches to describe the friction force in hydraulically driven multibody systems. This friction model incorporates the Coulomb, stiction and viscous friction, and is valid for both positive and negative relative tangential velocity. Therefore, the force F μ can be written as: 2.2 Nonlinear equations of motion of hydraulically actuated constrained multibody systems In the present work, a multibody system with n generalized coordinates, m holonomic constraints, l nonholonomic constraints, a hydraulic system with r discrete volumes and u normalized 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: In Eqs. (9)-(13), x is the n × 1 vector of generalized coordinates, which belongs to a domain of R n ; p is the r × 1 vector of pressures of the hydraulic system, defined in R r ; U is the u × 1 vector of normalized spool displacement variables existing in the hydraulic system; U ref is the u × 1 vector of reference normalized spool displacements; is the (m + l) × 1 vector of Lagrange multipliers; M (x) is the n × n mass matrix; Q x,ẋ, p is the n × 1 vector of generalized forces; H (x,ẋ, p, U) corresponds to the right-hand side of the pressure equations, which are built based on Eq. (1); W (U, U ref ) is the right-hand side of the valves dynamics equations, computed as defined in Eq. (4); C (x) is the m × 1 vector of holonomic constraints; and C nh x,ẋ is the l ×1 vector of nonholonomic constraints, linearly dependent on velocities. The matrices B (x) and D (x) are l × n and (m +l) × n, respectively, given by: where C x = ∂ 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:
First, the variationsx,ẋ,ẍ,˜ ,p,ṗ,Ũ andU with respect to this reference solution are defined as: Performing the Taylor expansion of the dynamic equilibrium equations: Equation (22) can be linearized by retaining up to firstorder terms. Using Eq. (16) yields: where the partial derivatives are evaluated for the reference solution. Similarly, the linearization of the velocity constraints in Eq. (15) leads to: where it has been used that D x 0 ẋ 0 = 0. Furthermore, the Taylor expansion of the hydraulic equations (12) with respect to the reference solution yields: Simplifying by using Eq. (19) and retaining up to firstorder terms in Eq. (25): Lastly, the Taylor expansion of the first-order differential equations (13) describing the valves dynamics results in: Using Eq. (20) and that the second-and higher-order derivatives in Eq. (27) are null yields:U To reduce the linearized equations of motion, the generalized coordinate partition of Ref. [30] is used. Given that the index-3 DAE system of equations (9)-(13) present m nonlinear holonomic constraints, the ncoordinates vector is split into m dependent coordinates x d and n − m admissible position coordinates x a : x = x a x d T . Moreover, the l nonholonomic constraints allow distinguishing between l dependent admissible velocitiesẋ ad and n − m − l independent admissible velocitiesẋ ai in the set of admissible velocitiesẋ a , and thereforeẋ a = ẋ aiẋad T . The same partition can be considered at position level: and for the vector of variationsx: The admissible dependent coordinatesx ad and the dependent coordinatesx d can be grouped in the set The following transformation matrix is defined: where is a (m + l)-square matrix, formed by the columns of matrix D (x) associated with the coordinatesx dd , and D ai (x) is built from the columns of D (x) associated with the independent coordinatesx ai . 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.
Step 1. Eliminate the Lagrange multipliers variations and reduce the linearized dynamic equations. By premultiplying Eq. (23) by T T x 0 , the Lagrange multipliers variations˜ are eliminated, given that T T x 0 D T x 0 = 0. By defining T 0 = T x 0 , the linearized dynamic equations (23) become: Step 2. Obtain a transformation at velocity level: express the velocitiesẋ in terms of the independent velocitiesẋ ai and positionsx. Next, the objective is to express Eq. (32) in terms ofx ai and its time derivatives. By using the linearized velocity constraints (24), the dependent velocities are written as follows: wherē From Eq. (33), the following transformation at velocity level is obtained: Step 3. Obtain a transformation at acceleration level: express the accelerationsẍ in terms of the independent accelerationsẍ ai , velocitiesẋ and positionsx. Secondly, the velocity constraints in Eq. (15) are differentiated with respect to time: The linearization of Eq. (36) yields: From Eq. (37), the dependent accelerationsẍ dd can be obtained as a function ofẍ ai , and the following relation at acceleration level is derived as: with Step 4. Obtain a transformation at position level: eliminate the dependent coordinatesx d by using the linearized holonomic constraints. To achieve the maximum reduction of the linearized equations of motion, the approach of Ref. [30] removes the dependent coordinatesx d by linearizing the holonomic constraints: Therefore, Eq. (41) allows expressing the dependent coordinatesx d as a function ofx ai andx ad : The use of Eq. (42) leads to the following transformation at position level: where Step 5. Express the velocitiesẋ ad in terms ofx ai , x ai andx ad : use of the time derivative of the holonomic constraints. Furthermore, an expression of the velocitiesẋ ad as a function ofx ai ,ẋ ai andx ad can be derived from the time derivative of the holonomic constraints and the nonholonomic constraints. Computing the time derivative of the holonomic constraints: the linearization of Eq. (45) with respect to the reference solution yields: Resorting to Eq. (46), the following transformation is obtained: where Step 6. Express the velocitiesẋ ad in terms ofx ai , x ai andx ad : use of the nonholonomic constraints. Finally, the nonholonomic constraints (11) are linearized: Substituting Eq. (47) in Eq. (49): and definingB x 0 = B x 0 T h 0 , the following expression for the velocitiesẋ ad is obtained: where 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. (35), (38) and (43), respectively, are used in the linearized dynamic equations (32) (result of Step 1) and the linearized hydraulic equations (26). Furthermore, including Eq. (51) (result of Step 5 and Step 6) and Eq. (28), yield the following linear ODE system of equations: whereV 0 =V x 0 ,ẋ 0 and the matrices m 0 , R 0 , S 0 , H 0 are given by: Note that the linear ODE system (53)-(56) is comprised of n − m − l linearized dynamic equations, l equations associated with the linearized nonholonomic constraints, r linearized hydraulic equations and u linear equations associated with the valves dynamics. 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.
By definingX = x aiẋaixadpŨ T , the linearized equations of motion (53)-(56) can be written as a first-order system of the formẊ = JX, where J is the Jacobian matrix: The size of the Jacobian matrix (61) is (2n−2m−l+r + u) × (2n − 2m −l +r + u). This Jacobian matrix represents the maximum possible reduction of the linearized equations of motion of a general hydraulically actuated multibody system with holonomic and nonholonomic constraints. In the particular case of a multibody system only with holonomic constraints (l = 0), the Steps 5 and 6 previously presented are not necessary and the linearized equations of motion (53)-(56) become:x and definingX = x aiẋaipŨ T , the first-order systemẊ = JX is obtained, with the Jacobian matrix:

Validation of the approach with a three-dimensional hydraulically actuated four-bar mechanism model
In this section, the procedure developed in Sect. 3 is validated with the forward dynamics simulation of a threedimensional hydraulically driven four-bar mechanism model. Moreover, a comparative analysis between the proposed approach and a conventional linearization counterpart procedure is performed. The comparison is made in terms of the preliminary steps required to compute the linearized equations of motion; the size of the resulting Jacobian matrix; the capability to analytically obtain the coefficients of the Jacobian matrix; the computational efficiency and the accuracy.

Description of the four-bar mechanism model with hydraulic actuator
The hydraulically actuated four-bar mechanism presents four rigid bodies: the ground link (inertial frame) is designated as body 1; the input link is denoted as body 2; the coupler is body 3, and the output link corresponds to body 4. The centres of mass G j , with j = {2 . . . 4}, correspond to the origins of the body frames. The origin of the inertial frame is located at O 1 , which also corresponds to the revolute joint connecting body 1 with body 2. The hydraulic actuation is performed by means of a hydraulic cylinder, consisting of the cylinder chamber and the piston-rod, whose masses are assumed to be negligible compared to those of the mechanism links. The endpoints of the hydraulic cylinder are placed at O c , which corresponds to a revolute joint between the cylinder and the ground link, and G 2 . The input link is connected to the coupler by means of the revolute joint C, and the revolute joint D allows the rotation of the output link with respect to the coupler. The output link is connected with the ground link by means of the revolute joint O 2 . A set of n = 10 generalized coordinates is used to describe the system, with the n × 1 vector of coordinates x given by: The position of G 2 is located by means of the coordinates x G 2 , y G 2 and z G 2 , and the orientation of the input link is determined with the angles ψ 2 , θ 2 and φ 2 , the latter being the rotation of the input link in the Y 1 Z 1plane; the angle φ 32 represents the relative rotation of the coupler with respect to the input link; and φ 43 corresponds to the relative rotation of the output link with respect to the coupler. Note that the constraints of the multibody system must ensure that the motion of the mechanism takes place in the Y 1 Z 1 -plane. Moreover, the length of the hydraulic cylinder (distance between the points O c and point G 2 ) is given by the coordinate s. Lastly, α corresponds to the angle between the hydraulic cylinder and the horizontal axis Y 1 . Figure 3 shows a three-dimensional view of the four-bar mechanism model, with the numbering of the bodies, all the coordinates and the body frames. The orientation matrices of the body frames, expressed as a function of the elemental rotation matrices, are given by:  Table 1 summarizes all the mechanical parameters (geometric and inertial) of the multibody system, including their numerical values, being shown in Fig. 4. The mechanical parameters shown in Table 1 are grouped in the set P m . Given that n = 10 coordinates are used, a total of m = 9 holonomic constraints are required, since the number of degrees of freedom of the multibody system is n g = n − m = 1. The set of holonomic constraints C (x) is given by: In Eq. (67), the first five constraints arise from the revolute joint connecting bodies 1 and 2 in O 1 , since this joint allows only the rotation of body 2 in the Y 1 Z 1 -plane. The absolute position vector of the rev- Fig. 4 View of the four-bar mechanism in the Y 1 Z 1 -plane: parameters L j and d 1 olute joint O 1 is r O 1 ; v is the unit vector along the X 2 -axis, expressed in the body frame 1, which is constrained to be parallel to the X 1 -axis by means of the fourth and fifth holonomic constraints; and u 1 , u 2 are unit vectors, perpendicular to the X 1 -axis. These vectors are computed as follows: In the next two holonomic constraints of Eq. (67), r Y O 2 and r Z O 2 are the Y and Z -components of the vector r O 2 , which is computed as follows: where 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

Description of the hydraulic circuit
The hydraulic circuit considered in this work consists of a double-acting hydraulic cylinder, a throttle valve, a directional control valve, a pump, a tank and the connecting hoses. For sake of simplicity, leakage in the hydraulic components is neglected. The hydraulic circuit is divided into three control volumes, denoted by V 1 , V 2 and V 3 . These control volumes are highlighted in Fig. (5) and are given by: In Eq. (72), V h 1 , V h 2 and V h 3 are the volumes of the hoses of the corresponding control volumes; A 2 and A 3 are the surfaces of the piston and piston-rod sides sections, respectively; and l 2 , l 3 are the lengths of the piston and piston-rod sides. These lengths are a function of the coordinate s and can be computed as: where c 1 is the length between the revolute joint O c and the base of the cylinder chamber, c 2 is the length of the piston and l is the length of the cylinder chamber. The dimensions c 1 , c 2 and l, and the lengths l 2 , l 3 are shown in detail in Fig. 5. The pressures of the control volumes p 1 , p 2 and p 3 can be computed, following Eq. (1), as follows: where B e 1 , B e 2 and B e 3 are the effective bulk moduli of the corresponding control volumes. These moduli are computed following Eq. (2): with B h and B c being the bulk moduli of the hoses and the hydraulic cylinder, respectively. The volume flow rates Q A 1 , Q A 2 and Q B 1 of Eq. (74) are defined as shown in Fig. 5. The parameters of the hydraulic circuit are listed in Table 2, with their corresponding numerical values being those of Ref. [5]. The hydraulic parameters shown in Table 2 are grouped in the set P h .
Equations (74) of the hydraulic circuit can be written in vector form as in Eq. (12):ṗ = H x,ẋ, p, U . Note that, in the hydraulic circuit of the present example, there is only one directional control valve, and thus only one normalized spool displacement U is required. Therefore, the equations of motion of the hydraulically actuated four-bar mechanism model are given by the dynamic equilibrium equations, which are computed following Ref. [44], the set of holonomic constraints (67), the hydraulic equations (74) and Eq. (13), which describes the dynamics of the directional control valve. The following index-3 DAE system is obtained: Note that, due to the absence of nonholonomic constraints, the matrix D (x) defined in Eq. (14) becomes C x (x).

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: where Given an equilibrium angle of the input link φ 0 2 , the values of the coordinates in the equilibrium configuration can be obtained from the holonomic constraints particularized for the equilibrium position: C x 0 = 0. The angles φ 0 32 and φ 0 43 are numerically obtained from nonlinear equations, while the remaining coordinates can be determined analytically with Eq. (78). Furthermore, by using the dynamic equilibrium equations and the hydraulic equations particularized for the equilibrium position: where U 0 = 0 andẋ 0 =ẍ 0 =ṗ 0 = 0 in the equilibrium, the pressures p 0 are obtained: with Note that, in Eq. (82), g( p 0 3 , x 0 , P m , P h ) allows obtaining the equilibrium pressure p 0 2 as a function of p 0 3 , and is given by: with the nonzero Lagrange multipliers being: In Eq. (85), μ 1 x 0 , P m and μ 2 x 0 , P m are functions of the coordinates in the equilibrium position x 0 and the mechanical parameters P m : Given that the system presents one degree of freedom, the input link angle φ 2 is chosen as independent coordinate. The sets of independent and depen-dent coordinatesx ai andx d , defined in Eq. (30), are therefore given by: The use of Eq. (64), particularized for the equilibrium of the four-bar mechanism model, leads to the following Jacobian matrix: where the coefficients ν k in Eq. (88), with k = {1 . . . 8}, are functions of the coordinates and pressures in the equilibrium, x 0 and p 0 , respectively, and the geometric and hydraulic parameters P m and P h : ν k = ν k x 0 , p 0 , P m , P h . 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. (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 (88) is (2n −2m −l +r +u)×(2n −2m −l +r +u) = 6×6.

Numerical results: validation of the approach
To validate the results obtained from the linearization, the forward dynamics simulations of the nonlinear fourbar mechanism multibody model and the linear system are performed. To carry out the numerical integration of the nonlinear system (performed with MATLAB ODE15s solver), the holonomic constraints in Eq. (76) are differentiated twice with respect to time, leading to the following index-1 DAE system: Note that the differentiation of the holonomic constraints may result in the violation of the constraints. Therefore, Baumgarte stabilization method [45] is used in Eq. (89) to avoid numerical drift, β and γ being the Baumgarte stabilization constants. To minimize the numerical drift of the constraints, while keeping the numerical cost of integration reasonably low, a numerical value of β = γ = 50 has been used. It has been verified that the results of the forward dynamics simulation are not dependent on the numerical values of these stabilization constants. Since the linearization has been performed around the equilibrium, corresponding to the neutral position of the directional control valve (U = 0), a reference normalized spool displacement U ref = 0 is considered in the simulation. The equilibrium position x 0 in Eq. (77) is considered as initial condition at position level: x (0) = x 0 , with the following numerical values: The exact numerical values in Eq. (90), corresponding to φ 2 (0) = 60 • , can be computed by solving the holonomic constraints for this equilibrium angle, as detailed after Eq. (78). Moreover,ẋ (0) = 0 10×1 at velocity level and the equilibrium pressures p 0 given by Eq. (82) are considered as initial pressures: p (0) = p 0 , with p 0 3 = 3.5 MPa: Note that the exact numerical value of p 2 (0) is computed by using Eq. (82). The term Q M (x) in Eq. (89) is the generalized force vector associated with the external torque M 2 , which is applied in the input link of the mechanism: The external torque M 2 , whose time evolution is described in Eq. (92) and shown in Fig. 6, linearly increases with time until t = t 1 , maintaining its maximum value M 0 until t = t 2 , when it is suddenly removed. Note that this torque deviates the multibody system from its equilibrium position. In the present case, t 1 = 1 s, t 2 = 1.5 s and M 0 = 250 N · m. Furthermore, the forward dynamics simulation of the linear system is performed. From Eq. (88), the linearized equations of motion are given by: The numerical integration of the linear system of equations (93) is performed from t = t 2 onwards, which corresponds to the instant of time when the external torque disappears. Therefore,φ 2 (t 2 ),ω 2 (t 2 ),p (t 2 ) andŨ (t 2 ) are considered as initial conditions for the numerical integration, with: where the numerical values of φ 2 (t 2 ),φ 2 (t 2 ) and p (t 2 ) are obtained from the forward dynamics simulation of the nonlinear system. The time evolution of the linear system can be predicted from the linear stability analysis. The eigenvalues are computed for three different scenarios, corresponding to different values of the friction force parameters F c , F s and σ 2 . Denoting as F * c , F * s and σ * 2 the numerical values of F c , F s and σ 2 used in Ref. [5], with F * c = 210 N, F * s = 830 N and σ * 2 = 330 Ns/m, Table 3 shows the numerical values of the eigenvalues λ k , with k = {1 . . . 6}, obtained from the Jacobian matrix (88). The results are shown with fourteen decimal digits for those readers who may wish to compare their results.
In the linear system of equations (93), it can be seen that the linear hydraulic equation ofṗ 2 can be obtained from a linear combination of the linearized equations ofṗ 1 andṗ 3 , which leads to one of the null eigenvalues shown in all the scenarios of Table 3. The second null eigenvalue arises from the linear dependence between the linear equation ofṗ 3 andφ 2 . Furthermore, the eigenvalue λ 6 = − 1 τ , which is associated with the linearized equation of the valve dynamics in Eq. (93), is obtained in all the cases, since this equation is decoupled. Lastly, the eigenvalues λ 3 , λ 4 and λ 5 correspond to the remaining three linearized equations. The numerical values of λ 3 , λ 4 and λ 5 are highly sensitive to the variation of the Coulomb friction F c , the static friction F s and the coefficient of viscous friction σ 2 .
Overdamped scenario. In the first case, with F c = F * c , F s = F * s and σ 2 = σ * 2 , the system is overdamped and λ 3 , λ 4 and λ 5 are negative real eigenvalues. Figure 7a shows the time evolution of the variation of the coordinate φ 2 (t) with respect to the equilibrium position φ 0 2 = 60 • , in the nonlinear and linear dynamic simulations. According to the time evolution of the external torque M 2 in Eq. (92), the angle φ 2 (t) increases and deviates from φ 0 2 for t < t 1 , with t 1 = 1 s. Next, for t 1 ≤ t ≤ t 2 , with t 2 = 1.5 s, M 2 is constant and the angle maintains its value over time. Lastly, for t > t 2 , the external torque vanishes and the input link approximately returns to its equilibrium position. Figure 7b shows in detail the time evolutions, in the nonlinear and linear cases, of the variations of φ 2 (t) with Overdamped scenario   respect to φ 0 2 = 60 • for t > t 2 . Moreover, Fig. 8a represents the evolution of p 1 and p 2 over time, and Fig. 8b shows the time evolution of p 3 . In all the figures, the nonlinear responses are denoted by NL and the linear results are labelled with L. It can be seen that the linear responses accurately reproduce the nonlinear results.
Underdamped scenario (high friction). Secondly, the numerical values of the friction parameters are reduced In this scenario, the system is underdamped: as shown in Table 3, the eigenvalue λ 5 remains real and negative, while λ 3 and λ 4 are, in this case, a complex conjugate pair of eigenvalues with negative real parts. Figure 9a shows the time evolution of the variation of the coordinate φ 2 (t) with respect to the equilibrium position φ 0 2 = 60 • , in the nonlinear and linear dynamic simulations. Note that, until t = t 2 = 1.5 s, the time evolution of φ 2 (t) is similar to that shown in the overdamped scenario in Fig. 7. Nevertheless, for t > t 2 , the system exhibits its underdamped behavior by oscillating around the equilibrium position φ 0 2 . These oscil- F * s and σ 2 = 1 100 σ * 2 . As in the previous case, the system is underdamped, with the corresponding eigenvalues shown in Table 3. Note that, despite there is very small variation in λ 5 , the real parts of λ 3 and λ 4 substantially decrease. The underdamped behavior is shown in Fig. 11a, b, with the time evolution of the variation of the coordinate φ 2 (t) with respect to the equilibrium position φ 0 2 = 60 • . Furthermore, the pressures p 1 (t) and p 2 (t) are shown in Fig. 12a, and p 3 (t) is plotted in Fig. 12b. Due to the reduction of the friction, the amplitude and frequency of the oscillations in Figs. 11a, b and 12a, b significantly increase with respect to those of the second scenario. In the particular case of no friction force (F c = F s = σ 2 = 0), the time evolutions of φ 2 (t) and the pressures are qualitatively similar to those of this scenario have 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. Moreover, a decrease in the numerical value of the maximum external torque M 0 in Eq. (92) 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 13 shows the friction force F μ presented in Eq. (8), with F c = 210 N, F s = 830 N and σ 2 = 330 Ns/m, together with its first-order Taylor approximation, denoted as T 1 μ . It can be seen that the linear approximation is excellent for low values of the velocityṡ, as is the case of the overdamped scenario. Nevertheless, the approximation becomes worse for larger values ofṡ, which is the case of the underdamped scenarios.

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 (15). The velocity constraints D (x)ẋ = 0 can be differentiated once with respect to time, leading to the following index-1 DAE system: where Eqs. (95) can be rewritten as the following DAE system: . (96) Eq. (95) can be expressed as follows: where the matrix A and the vector f are defined in Eq. (96). Denoting the first 2n + r + u components of F (X, U ref ) as F (X, U ref ), the following ODE system is obtained: Note that the last m + l equations of Eq. (97) provide the expressions of the Lagrange multipliers . The linearization of the ODE system (98) with respect to an arbitrary reference solution X 0 = x 0ẋ0 p 0 U 0 T leads to the following Jacobian matrix: In Eq. (99), N = 2n + r + u and J j represents the jth column of the Jacobian matrix, with j = {1 · · · N }. The term J j is computed as follows: where the partial derivative in Eq. (100) can be numerically computed by using central finite differences. In Eq. (100), h is the step used in the computation of the finite differences, and e j corresponds to the jth column of the identity matrix I N . Note that the Jacobian matrix (99) is (2n + r + u) × (2n + r + u).
A detailed comparison between the linearization approach proposed in Sect. 3 and this counterpart procedure is performed by means of the hydraulically actuated four-bar mechanism model described in Subsects. 4.1 and 4.2. First, a qualitative comparative analysis is carried out. The preliminary steps required to compute the Jacobian matrices (61) and (99) differ significantly. The counterpart approach first computes the index-1 form of the nonlinear DAE system (96); then obtains the first-order ODE system (98) by removing the equations associated with the Lagrange multipliers and finally performs the linearization, obtaining the Jacobian matrix (99). In contrast, the proposed linearization approach follows the steps described in Sect. 3. The DAE system (15) is first linearized along the reference solution; then, a coordinate partition in terms of independent and dependent coordinates is used to reduce the linearized dynamic equations and eliminate the variations of the Lagrange multipliers (see Eq. (32)); next, the linearized constraints at velocity and acceleration levels are used to express the complete set of velocity and acceleration variationsẋ and x in terms of the independent ones,ẋ ai ,ẍ ai ; and finally, the linearized holonomic constraints are used to express the variations of the dependent coordinatesx d in terms of the independent onesx ai . 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 coordinatesx ai and their time derivatives.
An important advantage of the proposed approach is that it leads to the maximum possible reduction of the linearized equations of motion. For the general case of a hydraulically driven multibody system with holonomic and nonholonomic constraints, the Jacobian matrix (61) is (2n − 2m − l + r + u) × (2n − 2m − l + r + u), while the conventional counterpart is a bulky linearization procedure that provides an augmented Jacobian matrix (see Eq. (99)), which is (2n + r + u) × (2n + r + u). Therefore, the proposed approach provides 2n − 2m − l + r + u eigenvalues, among which 2 (n − m − l) + r + u correspond to the real spectrum of the problem, and l spurious null eigenvalues are obtained, associated with the linearized nonholonomic constraints (54). In the particular case of a holonomic multibody system, the proposed approach leads to 2 (n − m) + r + u eigenvalues, which correspond to the real spectrum, and no spurious eigenvalues are obtained. In contrast, the use of the counterpart linearization approach results in 2n + r + u eigenvalues, with 2 (n − m − l) + r + u eigenvalues corresponding to the real spectrum, and 2 (m + l) spurious null eigenvalues associated with the m+l dependent coordinates. In the particular case of the four-bar mechanism model, the proposed approach leads to the Jacobian matrix (88), which is 6 × 6, obtaining the eigenvalues λ k , with k = {1 . . . 6}, of Table 3. Conversely, the Jacobian matrix (99) is 24×24, obtaining the six eigenvalues λ k of Table 3 and eighteen additional spurious null eigenvalues. Therefore, the use of the proposed approach allows for the reduction of the linearized equations of motion and the elimination of the spurious null eigenvalues associated with the dependent coordinates, while retaining all the stability information.
Another important advantage to highlight is the power of the proposed approach. This procedure allows generating 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 model, the coefficients ν k of the Jacobian matrix (88), with k = {1 . . . 8}, are computed analytically as a function of the coordinates and pressures in the equilibrium, x 0 and p 0 , respectively, and the geometric and hydraulic parameters P m and P h : ν k = ν k x 0 , p 0 , P m , P h . This significantly eases 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 (99) 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 (x) in Eq. (97) is required. For this reason, the Jacobian matrix (99) is numerically computed by using finite differences.
A quantitative comparative analysis is also included. The computational efficiency and accuracy of the proposed and counterpart approaches are compared. The computational efficiency is assessed by means of the required time of computation of the Jacobian matri-ces (88) and (99). For the particular case of the four-bar mechanism model, the average time of computation of the Jacobian matrix (88) is 0.005 s. In contrast, the required time for computing the Jacobian matrix (99) 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. Therefore, the proposed approach is more efficient, with a time of computation approximately three times lower for the four-bar mechanism.
Regarding the accuracy, the conventional counterpart provides an approximation of the Jacobian matrix, whose accuracy is highly sensitive to the step h used in the computation of the finite differences in Eq. (99). In contrast, the proposed approach leads to the exact linearized equations of motion. To compare the accuracy of both approaches, the relative errors 3 and 4 in the computation of the eigenvalues λ 3 and λ 4 of the overdamped scenario (see Table 3) are computed. Figure 14a, b show the variations of the relative errors 3 and 4 , respectively, with the step h. To ease the visualization, a logarithmic scale is used for the X-axis. In Fig. 14a, b, it can be seen that the relative errors 3 and 4 significantly increase for very low and high values of the step h. The counterpart approach provides the most accurate results of the eigenvalues λ 3 and λ 4 for an intermediate value of h 10 −5 , with 3 0 and 4 0. Therefore, the choice of an appropriate step is crucial for the accuracy of the alternative procedure. On the other hand, this problem does not arise in the proposed approach, since the linearization procedure presents a high capability to compute the exact Jacobian matrix. The use of a first-order finite difference formula (as for example the forward finite difference) to numerically calculate the Jacobian matrix could certainly decrease the computational burden of the counterpart approach due to a smaller number of function evaluations required, at the cost of decreasing the accuracy of the approximation.

Conclusions and future work
In this work, a linearization approach for systematically linearizing the equations of motion of hydraulically actuated multibody systems with holonomic and nonholonomic constraints has been presented. The proposed procedure provides the exact linearized equations of motion of the multibody system, and enables the computation of the Jacobian matrix analytically or numerically. The procedure has been validated by means of the forward dynamics simulation of a threedimensional four-bar mechanism model. By comparing the time evolution of the input link angle and the pressures of the hydraulic circuit in the nonlinear and linear systems, it has been shown that the linearized model accurately reproduces the results of the nonlinear multibody system. The procedure presents a great power, which is demonstrated by generating the Jacobian matrix as a function of the mechanical and hydraulic parameters of the three-dimensional four-bar mechanism model. This allows performing a linear stability analysis of the system, computing the eigenvalues for different scenarios. In this particular case, the eigenvalues have been computed for three-different set of values of the Brown and McPhee friction model parameters. The results show that these eigenvalues are highly sensitive to the Coulomb and static friction forces and the coefficient of viscous friction. Despite the linearized model reproduces well the results of the nonlinear system in the vicinity of the equilibrium configuration for all the scenarios, the results of the overdamped case improve compared to those of the underdamped cases, given that the linear approximation of the friction force becomes better in the overdamped scenario. The proposed approach is also compared with a counterpart procedure, presenting several advantages. First, the maximum reduction of the linearized equations of motion is obtained, which simplifies working with the resulting Jacobian matrix. In the performance of linear stability analyses, this leads to the elimination of the spurious null eigenvalues associated with the dependent coordinates of the multibody system model. Secondly, the proposed approach is more efficient than its counterpart. The computational efficiency is assessed by means of the required time of computation of both Jacobian matrices for the hydraulically driven four-bar mechanism. A time of computation of approximately three times lower was obtained with the proposed approach. Lastly, the accuracy of the procedure was also demonstrated. While the procedure developed in this work provides the exact Jacobian matrix of the system, the counterpart approach is highly sensitive to the step used in the numerical computation of the partial derivatives required in the Jacobian matrix. To illustrate this, the variation of the relative errors with the step, in the computation of the eigenvalues of the system, was shown. It is important to note that, in this study, 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 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. In future work, the proposed procedure will be used for developing state observers and linear feedback controllers of hydraulically actuated multibody systems.