Alexandr Klimchik

Anatol Pashkevich

Damien Chablat

Fundamentals of manipulator stiffness modeling using matrix structural analysis

Keywords: Manipulator stiffness modeling, Matrix Structural Analysis, Cartesian stiffness matrix, Hybrid architectures, Passive joints

published or not. The documents may come

Fundamentals of manipulator stiffness modeling using matrix structural analysis

Introduction

In many modern robotic applications, manipulators are subject to essential external loadings that affect the positioning accuracy and provoke non-negligible positioning errors caused by the compliance of mechanical components [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Slavkovic | A method for off-line compensation of cutting force-induced errors in robotic machining by tool path modification[END_REF]. For this reason, manipulator stiffness analysis becomes one of the most important issues in the design of robot mechanics and control algorithms. It allows the designer to achieve required balance between the dynamics and accuracy since the usual reduction of manipulator moving masses or cross-sections leads to increasing of achievable speed and acceleration but also deteriorating undesirable compliance errors. However, to make the stiffness analysis efficient, it should rely on a simple and computationally reasonable method that is able to deal with complex architectures including numerous closed-loops, flexible links, rigid connections, passive and elastic joints that are common for hybrid robots. At present, there exist three main techniques in this area, that are Finite Element Analysis (FEA), Matrix Structural Analysis (MSA) and Virtual Joint Modeling (VJM) [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF][START_REF] Yi | Geometric analysis of antagonistic stiffness in redundantly actuated parallel mechanisms[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF]. The most accurate of them is the FEA [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Wang | Finite Element Analysis and Comparison of Two Hybrid Robots-the Tricept and the TriVariant, Intelligent Robots and Systems[END_REF][START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF], which allows modeling links and joints with their true dimension and shape. However, this technique is usually applied at the final design stage because of the high computational expenses required for high order matrix inversion [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF][START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF]. In contrast, the VJM method is treated as the simplest one, it is based on the extension of the traditional rigid model by adding the virtual joints (localized springs), which describe the elastic deformations of the links, joints and actuators [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF][START_REF] Quennouelle | Instantaneous kinemato-static model of planar compliant parallel mechanisms[END_REF][START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF][START_REF] Klimchik | Stiffness matrix of manipulators with passive joints: computational aspects[END_REF]. This technique is widely used for serial and strictly parallel robots, but it can be hardly applied to the manipulators with more complex topology. The MSA is considered as a compromise technique, which incorporates the main ideas of the FEA, but operates with rather large elements such as flexible links connected by the actuated and passive joints in the overall manipulator structure [START_REF] Nagai | A systematic approach to stiffness analysis of parallel mechanisms[END_REF][START_REF] Marie | Elasto-geometrical modeling and calibration of robot manipulators: Application to machining and forming applications[END_REF][START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]. This obviously leads to the reduction of the computational expenses that are quite acceptable for robotics, but it requires some non-trivial actions for including of passive and elastic joints in the related mathematical model. From the other side, the MSA is very convenient for a description of complex structures with numerous closed-loops and cross-linkage. For this reason, this paper focuses on some enhancement and generalization of the MSA technique for robotic applications providing a compromise between complexity of the stiffness model generation and complexity of subsequent computations. Several reviews of existing works on manipulator stiffness analysis can be found in the literature [START_REF] Klimchik | Serial vs. quasi-serial manipulators: Comparison analysis of elasto-static behaviors[END_REF][START_REF] Gonçalves | A comparison of stiffness analysis methods for robotic systems[END_REF][START_REF] Yeo | Design and analysis of cable-driven manipulators with variable stiffness[END_REF][START_REF] Liu | Stiffness Modeling of Parallel Mechanisms at Limb and Joint/Link Levels[END_REF][START_REF] Cao | A method for stiffness modeling of 3R2T overconstrained parallel robotic mechanisms based on screw theory and strain energy[END_REF][START_REF] Subrin | New Redundant Architectures in Machining: Serial and Parallel Robots[END_REF][START_REF] Taghvaeipour | On the elastostatic analysis of mechanical systems[END_REF][START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF][START_REF] Pinto | A methodology for static stiffness mapping in lower mobility parallel manipulators with decoupled motions[END_REF][START_REF] Chen | Conservative congruence transformation for joint and Cartesian stiffness matrices of robotic hands and fingers[END_REF][START_REF] Klimchik | Efficiency evaluation of robots in machining applications using industrial performance measure[END_REF], where the authors addressed different aspects of the FEA, MSA and VJM and particularities of their application in robotics. These reviewers cover results starting from the early works of [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF] [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF] till nowadays, excluding last few years. Among recent contributions devoted to the MSA, it is worth mentioning the work of Cammarata [START_REF] Cammarata | Unified formulation for the stiffness analysis of spatial mechanisms[END_REF], who introduced the notion of the Condensed Stiffness Matrix that merges together flexibilities of the link and adjacent joints. Further, the obtained stiffness matrices of a two-node super-elements are used in a traditional for MSA way, which includes the manual merging of the lines and columns in the global stiffness matrix. Besides, despite the apparent simplicity, this technique does not allow direct computing of the Cartesian Stiffness Matrix, which is the main outcome of the manipulator stiffness analysis. In addition, here the influence of the external loading and buckling analysis cannot be executed in a simple way. Another contribution in this area [START_REF] Soares Júnior | Stiffness analysis of multibody systems using matrix structural analysis-MSA[END_REF] deals with MSA-based stiffness analysis of a parallel manipulator composed of several L-structures. The authors applied classical MSA, with the connectivity matrix for lines/columns merging. Although they took into account both the link and joint flexibility, there are still several open questions related to the definition of link/joint stiffness properties and invertibility of local stiffness matrices describing the L-structures. Another useful extension of the MSA for the case of links and joints with non-linear stiffness was proposed in [START_REF] Detert | Extended Procedure for Stiffness Modeling Based on the Matrix Structure Analysis[END_REF]. Here, a passive revolute joint with ball bearings was presented as an element with a rank-deficient force-dependent stiffness matrix, whose parameters were estimated experimentally. The relevant computational procedure included several iterations of conventional MSA linear model. This technique was applied to PARAGRIP handling system and validated by measurements of the end-effector deflection under vertical load. There are also several works the deal with the MSA application to the stiffness analysis of particular parallel and serial manipulators. In [START_REF] Shi | Static stiffness modelling of EAST articulated maintenance arm using matrix structural analysis method[END_REF] the MSA method was applied to EAST articulated maintenance arm with 11 degrees of freedoms (EAMA robot), which is used for remote inspection of inner components inside the vacuum vessel. In [START_REF] Mousavi | Model-Based Stability Prediction of a Machining Robot[END_REF] the MSA technique was employed to obtain a dynamic model of the industrial machining robot ABB IRB 6660 in order to predict vibration instability in machining (chatter). In [START_REF] Azulay | Comparative analysis of a new 3×PPRS parallel kinematic mechanism[END_REF] the MSA was applied to derive the static stiffness model of 9-dof redundant reconfigurable 3×PPPRS parallel manipulator for meso-Milling Machine Tool (RmMT). To our knowledge, the most essential contribution to the robot-oriented modification of the MSA was done by Deblaise and his co-authors [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]. They proposed a general technique that is able to take into account passive joints and rigid connections, which are common for parallel robots, without modification of 12 12  stiffness matrices describing structural elements (in contrast to some authors who integrated the passive joints and rigid connections in the link stiffness matrices). The properties of passive joints and rigid connections were described by means of matrix linear constraints, which complemented the classical MSA formulation containing the stiffness models of individual elements. This idea allowed to simplify assembling of the MSA global matrix and avoid tedious merging of the matrix rows and columns. To integrate these two sets of equations (conventional MSA relation and additional constraints) the authors used the energy approach with Lagrange multipliers to obtain the desired extended stiffness matrix describing all force-deflection relations for the entire robotic manipulator. Recent advances in computing power of commercially available facilities motivate re-thinking in the selection of the stiffness analysis technique for the robotic manipulators, making reasonable some increase of computational expenses related to the MSA method. This allows the user to apply the MSA to the stiffness analysis of complex manipulators with numerous closed-loops, preloadings and external forces/torques while operating with matrices of relatedly higher dimension that within the traditional MSA are reduced to the minimum size by means of rows/columns merging. This paper proposes an advancement in the MSA-based stiffness modeling of robotic manipulators allowing to analyze in a similar way both serial, parallel and hybrid architectures containing closed-loops, flexible links, and rigid connections, passive and elastic joints, flexible and rigid platforms, external wrenches and preloading, etc. The proposed approach leads to essential simplification of the stiffness model development while slightly increasing the time for the stiffness analysis stage. In addition, it also provides the user with additional data on the manipulator internal force/torque and corresponding deflections in links/joints.

MSA method background

MSA-based stiffness model of a separate link

In classical mechanics, the stiffness properties of the cantilever beam (fixed at one side, see Figure 1a) are described by the Hook's law that defines a linear relation between the applied external wrench (force/torque) W and corresponding deflection t at the free-end   W K t [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF] where K is 6 6  stiffness matrix. It should be mentioned that here t is 6-dimensional deflection vector that includes both In general case, the stiffness matrix K from the Hook's law is symmetric and positive definite but may include a number of off-diagonal elements [START_REF] Timoshenko | Theory of elasticity[END_REF]. For typical beams commonly used in practice (with regular cross-section) the stiffness matrix can be computed analytically as follows

uv ES L E Iz E Iz LL E Iy E Iy LL GJ L E Iy E Iy L L E Iz E Iz L L                                    K (pourquoi u, v ?) (2)
where L the beam length, S is the beam cross-section area, Iy , Iz are the second moments, J is the polar moment, E and G are Young's and Coulomb's modules of the beam material, respectively. It should be also stressed that for the cantilever beam, both the wrench W and deflection t are usually expressed with respect to the coordinate system attached to the beam's fixed-end. However, in general case, the vectors W and t may be presented in the global coordinate system, which requires some revision of ( 1) and ( 2). In fact, for small angular displacements, the deflection and wrench vectors in the global system can be presented as

global 33 global 3 3 66                   p R0 p φ 0 R φ (3)
global 33 global 3 3 66

              F R0 F M 0 R M (4)
where the orthogonal matrix R defines the orientation of the local coordinate system relative to the global one.

T T                             Fp R 0 R 0 M 0 R R K φ 0 (5)
that gives a simple rule for transforming the local stiffness matrices of all mechanical components to a single coordinate frame.

For the unsupported beams (with two non-fixed ends, see Figure 1b) that are used in Matrix Structural Analysis (MSA) as principal components, it is necessary to define the deflections and wrenches for both sides. The latter will be further referred to as "u" and "v" or "1" and "2". In this case, the stiffness model is presented in an extended form

                     W K K t W K K t (6)
that relates the deflections on both sides WW should satisfy the static equilibrium equation that defines linear dependence between the matrix rows in eq. (6).

To find the desired 66  blocks, let us consider two special cases that allow us to apply directly the cantilever beam equation ( 1) that includes 66  matrix K . In the first case, let us orient the global system axes in the standard way and locate its origin at the left-end of the beam. This corresponds to 1  t0 , 22   W K t and leads to simplification of eq. ( 6)

1 12 2 2 22 2     W K t W K t (7)
Relevant static equilibrium equation is written at the point "1" allows us to express the wrenches 12 , WW via the force 2 F and torque 2 M at the point "2" as

22 1 2 2 2 2 ;                FF W MM L W F ( 8
)
where L is the beam length vector directed from the point "1" to the point "2" (i.e. along the beam principal axis). Applying further the cantilever beam standard equation 22   W K t one can get explicit expressions for the matrices 12 22 , KK : ;

[] T           I0 K K K K LI (9)
where []  L denotes the 3 3  skew-symmetric matrix derived from the vector L . Similarly, in the second case, the global system is located at the right-end of the beam (with the same orientation as in the previous case), which yields 2  t0 and simplifies eq. ( 6) to

1 11 1 2 21 1     W K t W K t (10)
The static equilibrium equation leads to the following expressions

2 11 1 1 1 1 ;                  FF WW M M F L ( 11
)
where L is the same beam length vector (directed from the point "1" to point "2"). However, the cantilever beam equation should be slightly modified to take into account standard directions of the coordinate axes

3 3 3 3 3 1 3 1 33 zz zz                    R 0 R 0 W K t 0 R 0 R (12)
where ( 1, 1, 1)

z diag     R is 3 3
 rotation matrix around z axis by the angle  , which changes directions of the axes x, y. The latter allows us to obtain explicit expressions for 11 21 , KK   ;

zz z T z                           I0 R 0 R 0 K K K K LI 0 R 0 R (13)
The physical interpretation of the matrices 11 12 21 22 , , , K K K K can be done in the following way. The stiffness matrix 11 K describes the force/torque reaction at the beam left-end point caused by the left-end deflection 1 t . The stiffness matrix 12 K describes the force/torque reaction at the beam left-end point caused by the right-end deflection 2 t . The stiffness matrix 21 K describes the force/torque reaction at the beam right-end point caused by the left-end deflection 1 t . The stiffness matrix 22 K describes the force/torque reaction at the beam right-end point caused by the right-end deflection 2 t . The physical meaning of , , , K K K K is also illustrated for the simplest case of the 1D linear spring (Figure 2), where 66  matrices are reduced to the scalar stiffness coefficients

x x x 2 x  1 x  1 x  2 x  2 F 1 F 1 F 2 F 2 2 F k x    k k k 1 1 F k x        1 1 2 1 1 2 2 2 2 1 x x x k k x x x F k F k F k F k                                    1 1 1 2 2 2

Figure 2

Stiffness models of a simple spring with fixed and free ends

For the beam with a regular shape (with uniform cross-section) the above stiffness matrices can be computed analytically using the following expressions

ES L E Iz E I LL E Iy E Iy LL GJ L E Iy E Iy L L E Iz E Iz L L                                      K (14)
ES L E Iz E Iz LL E Iy E Iy LL GJ L E Iy E Iy L L E Iz E Iz L L                                        K ( 15
)
ES L E Iz E Iz LL E Iy E Iy LL GJ L E Iy E Iy L L E Iz E Iz L L                                        K ( 16
)
ES L E Iz E I LL E Iy E Iy LL GJ L E Iy E Iy L L E Iz E Iz L L                                      K ( 17
)
where all notations have the same meaning as in the cantilever beam stiffness matrix [START_REF] Slavkovic | A method for off-line compensation of cutting force-induced errors in robotic machining by tool path modification[END_REF]. Similar to the cantilever bean case, the stiffness matrix of the unsupported beam can be transformed to the global system using an extended version of eq. ( 5)

T T T T                          R 0 0 0 R 0 0 0 K K K K 0 R 0 0 0 R 0 0 0 0 R 0 K K K K 0 0 R 0 0 0 0 R 0 0 0 R (18)
where the orthogonal matrix R defines the orientation of the local coordinate system relative to the global one.

MSA-based stiffness models of simple linkages

Two-link serial system with rigid connection

Now let us apply the MSA method to a mechanical system composed of two beam elements of lengths L 1 , L 2 with the rigid connection between them (Figure 3). It is assumed that the beams are not aligned and the angle between them is equal to q . Here, each beam has two local coordinate systems: 0, 1 for the first beam and 2, 3 for the second one, as shown in Figure 3. Let us assume that the global coordinate system is located at the left-end of the first beam (point 0) and the global coordinate axes are aligned with the local ones. Under this assumption, the force deflection relations for the first link can be written as 11 00 11 12 11 11 [START_REF] Yeo | Design and analysis of cable-driven manipulators with variable stiffness[END_REF] 22

                  Wt KK Wt KK (19)
where the upper superscript denotes the link number and relevant matrices are computed using eqs. ( 14)- [START_REF] Marie | Elasto-geometrical modeling and calibration of robot manipulators: Application to machining and forming applications[END_REF]. For the second beam, the similar relation should be written taking into account that the local axes 22 ( , )

x z and 33 ( , )

x z are turned with respect to the global ones. This yields the following force-deflection relation 22 22 11 12 22 33 [START_REF] Yeo | Design and analysis of cable-driven manipulators with variable stiffness[END_REF] 22

TT TT                            Wt Q K Q Q K Q Wt Q K Q Q K Q
(avec K 1 la matrice de rigidité de la poutre 1 et K 2 la matrice de rigidité de la poutre 2) [START_REF] Gonçalves | A comparison of stiffness analysis methods for robotic systems[END_REF] where ( , )

qq yy diag  Q R R is the 6 6
 matrix that is expressed via the orthogonal rotation matrix q y R describing rotation around y-axis by the angle q . For further convenience let us denote 22

q ij i T j    K Q K Q
                                   K K 0 0 Wt K K 0 0 Wt Wt 0 0 K K Wt 0 0 K K
(doit on écrire la deuxième dimension pour un vecteur?)

Further, taking into account the boundary condition 0  t0 the above system can be reduced down to

                          K 0 0 Wt W 0 K K t Wt 0 K K (23)
Furthermore, the rigid connection between the beams yields relations 12    tt and 12  WW , where the last one directly follows from Newton's third law. The latter agrees on us to sum-up the first and second lines in the system [START_REF] Cao | A method for stiffness modeling of 3R2T overconstrained parallel robotic mechanisms based on screw theory and strain energy[END_REF] and obtain the final expression allowing us to compute deflections at the beam ends caused by external loading 3 W applied to the free end of the considered 2-beam system

qq qq                 t K K K 0 W t KK (24)
It is worth mentioning that here the 12 12  matrix is invertible, and equations can be easily solved with respect to the deflections 23 ,  tt . In particular, an analytical solution for the free-end deflection 3 t caused by 3 W may be presented in the following form

33 C   W K t ( 25
)
where C K is the Cartesian stiffness matrix of the two-link system under study  

qq C qq      K K K K K K ( 26
)
that depends on both beams stiffness parameters included in the matrices 12 ,

q ij ij
KKand relative orientations of the links described by the matrix Q . It is worth mentioning that here the matrix C K is non-singular and invertible. Hence, for this simplest case study, the MSA approach allowed rather easy compute the system stiffness matrix. It is clear that this approach can be also generalized for the multi-beam serial systems with rigid connections, however, the final analytical expression will include the matrix inversion of higher dimension. Besides, it is worth mentioning that the above-presented example implements the general MSA assembling technique that usually produces the matrices similar to [START_REF] Subrin | New Redundant Architectures in Machining: Serial and Parallel Robots[END_REF], which includes a number of matrix sums.

Two-link serial system with a passive joint

Let us apply the MSA method to a mechanical system composed of two beam elements of lengths L 1 , L 2 with a revolute passive joint between them (Figure 4). In structural mechanics, this type of connection is also called as a pin joint. It is assumed that the passive joint axis is directed along y 1 and y 2 , the beams are not aligned and the angle between them is equal to the passive joint coordinate q . It is worth mentioning that this problem cannot be solved straightforwardly using the standard MSA method, which was primary designed for rigid and elastic connections. However, it is possible to apply some special techniques to adjust MSA for the passive joint case. Figure 4 Two-link serial system with a passive joint It is clear that the basic MSA equations for this 2-beam system have the same structure as above

                          K 0 0 Wt W 0 K K t Wt 0 K K (27)
but here the matrices 2q ij K are not constant and depend on the passive joint coordinate q , which may vary under the influence of the external loading 3 W . However, in contrast to the previous case, where the connection is described by equations 12    tt and 12  WW , here these relations must be replaced by

    12 2 1 1 2 ;; r r p p             Λ t t 0 Λ W W 0 Λ W 0 Λ W 0 (28)
where the matrices (1,1,1,1, 0,1)

r diag  Λ
and (0, 0, 0, 0,1, 0)

p diag  Λ
describe the passive joint geometry, and the superscripts 'r' and 'p' are referred to the "rigid" and "passive" connections. The first of these relations takes into account that the connection 1 and 2 ensures equality of all components of 12 ,

 tt except of 12 , yy  (since the passive joint axis is directed along y axis). The second group of relations shows that the passive joint ensures the validity of the 3 rd Newton law for five wrench components only (except of y M ), but the torques 1 y M and 2 y M are equal to zero. It should be noted that in general case the matrices r Λ and p Λ may be non-diagonal, but their rank remains the same and is equal to 5 and 1, respectively (as in the case of diagonal matrices). Combining equations [START_REF] Pinto | A methodology for static stiffness mapping in lower mobility parallel manipulators with decoupled motions[END_REF] and [START_REF] Chen | Conservative congruence transformation for joint and Cartesian stiffness matrices of robotic hands and fingers[END_REF]

                                     Λ Λ 0 0 Λ K 0 0 t 0 0t 0 Λ K Λ K 0 t Λ K Λ K Λ K W 0 K K (29)
that should be solved with respect to 3 t . It should be mentioned that here the system matrix is non-square (of the size 30 18  ), but it contains 12 zero rows that are caused by zero diagonal components in r Λ and p Λ (each r Λ produces one zero-row and each p Λ produces five zero-rows). To eliminate the zero-rows that are useless here, let us introduce the modified non-singular matrices

  ** 1 0 0 0 0 0 0 1 0 0 0 0 ; 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 rp         ΛΛ ( 30
)
and rewrite the system (29) in the reduced form (with a square matrix of size [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]

q r q qq r                                           Λ Λ 0 0 Λ K 0 0 0 t 0 t 0 Λ K Λ K 0 t Λ K Λ K Λ K W 0 K K (31)
This system can be presented in more compact form after defining relevant block matrices

q q q qq                                Λ Λ 0 0 Λ K 0 AB ΛK 0 ΛK ΛK Λ K Λ K C 0 K D K (33)
Further, after elimination 12 ,  tt one can get the following expression for the wrench applied to the system end-point

33 C   W K t ( 34
)
where

1 C     K D C A B ( 35
)
is the desired Cartesian stiffness matrix of the two-link system under study. It has a similar structure as (26) corresponding to the case of the rigid connection. However, it is easy to verify numerically that here the stiffness matrix C K is rank-deficient, which is the result of the passive joint connection. Hence, for this case study, some modification of the MSA approach allowed us to compute the rank deficient stiffness matrix. However, relevant transformations include a number of non-trivial steps, which will be generalized in the following Sections.

Two-link serial system with an elastic joint

Let us apply the MSA method to a mechanical system composed of two beam elements of lengths L 1 , L 2 with an elastic joint between them (Figure 5). It is assumed that the elastic joint axis is directed along y 1 and y 2 , the beams are not aligned and the angle between them is equal to the elastic joint coordinate q . It is clear that the basic MSA equations for this 2-beam system have the same structure as in the cases (a) and (b)

qq qq                                K 0 0 Wt W 0 K K t Wt 0 K K (36)
Here the matrices 2q ij K are not constant and depend on the elastic joint coordinate q , which may vary under the influence of the external loading 3 W . However, in contrast to the cases of the rigid and passive connections, here the joint static equations must be replaced by

  12 2 2 1 2 1 ;
()

r ee q K              Λ t t 0 W W 0 Λ Λ t t W (37)
where q K is the stiffness coefficient of the elastic joint, the matrices (1,1,1,1, 0,1)

r diag  Λ and (0,0,0,0,1,0) e diag  Λ
describe the passive joint geometry, and the superscripts 'r' and 'e' are referred to "rigid" and "elastic" connections. The first of these relations takes into account that the connection 1 and 2 ensures equality of all components of 12 ,  tt except of 12 , yy  (since the elastic joint axis is directed along y). The second group of relations shows the elastic joint ensures validity of the 3 rd Newton law for all wrench components, and the third relation describes the Hooke's law for the elastic joint.

Base x0 y0 z 0 y1, y2 z 1 z3 x 3 y 3 z 2 x1 q L 1 x 2
                              Λ Λ 0 0 t Λ Λ K Λ Λ K 0 t 0 K K K t W 0 K K (38)
that should be solved with respect to 3 t . It should be mentioned that here the system matrix is non-square of the size [START_REF] Subrin | New Redundant Architectures in Machining: Serial and Parallel Robots[END_REF] 18  , but it contains 6 zero rows that are caused by zero diagonal components in r Λ and e Λ (each r Λ produces one zero-row and each e Λ produces five zero-rows). To eliminate the zero-rows that are useless here, let us introduce the modified non-singular matrices * r Λ of the size 5 6  (see eq. ( 30)) and

  * 0 0 0 0 1 0 e  Λ
. This allows us to rewrite the system [START_REF] Deblaise | Contribution à la modélisation et à l'étalonnage élasto-géométriques des manipulateurs à structure parallèle[END_REF] in the reduced form (with a square matrix of size [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]

                                  Λ Λ 0 0 t Λ Λ K Λ Λ K 0 t 0 K K K t W 0 K K (39)
This system can be simplified after defining relevant block matrices and presented as

KK                              0 ΛΛ A Λ Λ K Λ B Λ K K KK C 0 K D K (41)
Further, after elimination 12 ,  tt one can get the following expression for the wrench applied to the system end-point

33 C   W K t ( 42
)
where

1 C     K D C A B ( 43
)
is the desired Cartesian stiffness matrix of the two-link system under study. It has a similar structure as ( 26) and ( 35) corresponding to the cases of the rigid and passive connections. It is easy to verify numerically that here the stiffness matrix C K is full-rank, which is in a good agreement with the nature of the elastic joint.

Two-link parallel system with passive joints

Let us apply the MSA method to a closed-loop mechanical system composed of two beam elements with three revolute passive joints, which connect the beams to the rigid base at the left-hand side and to the end-effector on the right-hand side (Figure 6). It is assumed that the passive joint axes are directed along y and the angle between the beams is equal to q .

                                   K K 0 0 Wt K K 0 0 Wt Wt 0 0 K K Wt 0 0 K K (44)
where the matrices 2q ij K depend on the system geometric parameter q . However, because of the non-rigid connections between the beams and the base, here 0 0  t and 2 0  t that does not allow us to reduce the matrix dimension in a trivial way (as it was done in all previous examples, where the lines corresponding to 0 W were simply removed). For the considered system, the left-hand side connections between the beam and rigid base (via the passive joints) produce the following equations 02 02 ;

;

rr pp           Λ t 0 Λ t 0 Λ W 0 Λ W 0 (45)
where the first group describes the passive joint geometry and the second group follows from the 3 rd Newton's law (the matrices r Λ and p Λ are defined above in the example (b)). For the right-hand side connections between the beams and end-effector, it is necessary to define additional variables t and e W describing the end-effector deflection and wrench, respectively. Using these notations, the end-effector wrench may be expressed as follows 31 ()

r e    W Λ W W (46)
Besides, the passive joint geometry and the 3 rd Newton's law at the right-hand side yield equations

    13 13 ;
;

rr pp               Λ t t 0 Λ t t 0 Λ W 0 Λ W 0 (
                                               0 Λ 0 0 0 0 0 00 Λ 0 0 0 0 Λ 0 0 Λ 0 0 0 0 ΛΛ 0 Λ K Λ K 0 0 0 0 Λ K Λ K 0 0 0 0 00 ΛK 0 0 W 0 1 2 3 22
q p pp r r r qq r q qq                                   t t t t t Λ K 0 00 Λ K Λ K 0 Λ K Λ K Λ K Λ K 0 (48)
that should be solved with respect to t . It should be mentioned that here the system matrix is non-square of the size 54 30  , but it contains 24 zero rows that are caused by zero diagonal components in r Λ and p Λ (each r Λ produces one zero-row and each p Λ produces five zero-rows). To eliminate the zero-rows, let us use the modified non-singular matrices * r Λ and * p Λ (defined in example (b)) and rewrite the system (48) in the reduced form

                                             Λ 0 0 0 0 0 00 Λ 0 0 0 0 Λ 0 0 Λ 0 0 0 0 ΛΛ 0 0 Λ K Λ K 0 0 0 0 Λ K Λ K 0 0 0 0 00 Λ 0 W 0 1 2 3 22
                                  t t t t t K Λ K 0 00 Λ K Λ K 0 Λ K Λ K Λ K Λ K 0 (49)
where the main matrix of size 30 30  is obviously rank-deficient (since the 29 th column is equal to zero). The latter is in a good agreement with physical nature of the system, where the end-effector is connected to the beams via the passive joints. To obtain the final expression, the reduced system can be presented in a compact form where

                                   Λ 0 0 0 0 00 Λ0 0 0 Λ 0 0 Λ 0 0 0 Λ Λ AB Λ K Λ K 0 0 Λ K Λ K 0 0 00 Λ K Λ K 00 Λ K Λ K 6 *
1 C     K C A B ( 53
)
is the desired Cartesian stiffness matrix of the two-link closed-loop system under study. It has a similar structure as in the above examples. Also, it is easy to prove analytically that the stiffness matrix C K is rank-deficient (here ( ) 5 rank  B because of a zero column), which is the result of the passive joint connection to the end-effector. Summarising all above-presented case studies, one can conclude that some modifications of the MSA approach allowed us to obtain the desired Cartesian stiffness matrix that may be either a full-rank or rank-deficient one. However, relevant techniques include a number of non-trivial steps, which will be generalized in the following Section to be applied to more complex structures.

General methodology of classical MSA

As follows from the above case studies, basic ideas of the classical MSA can be successfully used in robotics but some enhancement is required in order to take into account numerous passive and actuated connections. Before presenting the proposed enhancement, let us remind briefly basic steps of the classical MSA that is perfectly suited for computer-aided analysis of complex structures such as trusses, bridges, high voltage towers, etc. In general, the MSA based stiffness analysis includes the following steps:

Step 1: Decoupling the original system into generic structural members such as beams, rods, etc. and presenting the system as a set of member elements connected at the nodes (MSA idealization).

Step 2: Creating "free-free" master stiffness models for each element using 12 12  stiffness matrices and presenting them in the global coordinate system (MSA members modeling).

Step 3: Creating the node-element connectivity matrix that defines links between model members via rigid, elastic or passive connections (MSA members connecting).

Step 4: Merging the individual stiffness matrices to the global stiffness matrix for the entire structure using the node-element connectivity matrix and overlapping technique (MSA assembling).

Step 5: Defining nodal loads and system supports, incorporating them in the global stiffness model, application of the boundary conditions/constraints and reducing the global stiffness matrix dimension (MSA constraints).

Step 6: Solving the resulting set of the reduced equations for unknown nodal displacements corresponding to given external loads (MSA solving).

Step 7: Computing reaction forces at the supports, the nodal internal forces and stresses for all system members (MSA post-processing).

It should be noted that practical computer implementation of the MSA method significantly differs from the manual technique presented in the previous Section. There are two main differences: (i) at the assembling stage the individual stiffness matrices are not expanded but are directly merged through the use of special "freedom pointer array"; (ii) both the individual and global stiffness matrices are stored using a special format that takes advantage of symmetry and sparseness. There are also some differences in the application of boundary conditions and the global stiffness matrix reduction. Nevertheless, relevant computer-oriented routines ensure satisfaction of two basic rules of structural mechanics: compatibility of displacement and force equilibrium, which are equivalent to the simple merging of 6-dimensional rows and columns in the global stiffness matrix. Hence, some modifications of these routines are required to apply MSA in robotics where the size of the global stiffness matrix is not very high but the connections between the structural members (passive and actuated joints) should be treated in a different manner.

MSA enhancement for robotic manipulators

In contrast to structural mechanics where the main interest is in the area of the nodal displacements, reaction/internal forces and stresses, the stiffness analysis in robotics concentrates on computing of the Cartesian stiffness matrix defining the force-deflection relation for the manipulator end-effector. Let us present the MSA based technique adapted for robotic applications, which is issues from previous works [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF][START_REF] Deblaise | Contribution à la modélisation et à l'étalonnage élasto-géométriques des manipulateurs à structure parallèle[END_REF] and also contains some new contributions and generalizations.

3.1

MSA models of manipulator links and platforms

Modeling of a flexible link

If the link flexibility is non-negligible, the 2-node "free-free" stiffness model should be used, which can be obtained either from the approximation of the link by a beam or using the CAD-based technique [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF] allowing to evaluate the stiffness parameters taking into account the link real shape and geometry. In both cases the link is described by the linear matrix equation, containing [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF] 12  "free-free" stiffness matrix ,,,

( ) ( ) 11 12 ( ) ( ) 21 22 1 12 2 ij ij ii ij ij jj                    Wt KK Wt KK (
ij ij ij ij
KKKK are 66  stiffness matrices. It should be noted that the system(54) ensures both the displacement and force equilibriums for the flexible link, it includes 12 scalar equations and 24 scalar variables contained in six-dimensional vectors , ij  tt and , ij WW. Besides, it can be proved the rank deficiency of this model is equal to 12. It is clear that in general case when the link is arbitrary oriented with respect to the global coordinate system, the above equations should be slightly modified by simply rotating the local stiffness matrices

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 12 11 12 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 21 22 21 22
g ij g ij ij ij ij ij ij ij g ij g ij ij ij ij ij T T i ij T T j                      K K Q K Q Q K Q K K Q K Q Q K Q (55) where ( ) ( ) ( ) ( ,
) ij ij ij diag  Q R R is composed of two similar orthogonal matrices () ij
R defining the orientation of the local coordinate system of the link (ij) relative to the global one and left superscript "g" indicates that matrix is presented in the global coordinate system.

Modeling of a rigid link

If the link flexibility is negligible, the above stiffness model should be replaced by two types of equations describing the displacement and force equilibriums. The first of them can be treated as a simple "rigidity constraint" that keeps constant the distance between the nodes i and j. Applying to this link the rigid body kinematic equations, one can get the following relations between the nodal displacements [ ; ]

i i i    p φ t , [ ; ] j j j    p φ t
expressed via the linear and angular components ()

ji ij j i i         φφ pp d φ (56)
where the vector ( )

( ) ( ) ( ) [ , , ] ij ij ij ij x y z T ddd  d
describes the link geometry and is directed from the i th to j th node. This constraint can be also rewritten in the matrix form

3 3 3 3 3 3 61 3 3 3 3 3 3 3 3 () ] [ i T j j i                    I t d I 0 I t I0 0 0 (57)
where () [] ij  d denotes the 3 3  skew-symmetric matrix derived from the vector () ij d . Further, after definition 6 6

 block matrix () 33
3 3 3 ) 3 6 ( 6 [ ] T j i i j         d I D 0I (58)
the displacement constraint for the rigid-link can be presented in the following form

() 6 6 6 1 ij i j            t 0 t DI (59)
that is convenient for further aggregation of the stiffness model components.

The second group of equations describing the force equilibrium can be derived using the rigid body static equations, which yields the following relations

() 0 0 ij j i j ji      F F d MM F (60)
where , ij FF and , ij MM denote the forces and torques applied at the nodes i and j respectively. In a matrix form, they can be presented as follows

() 0 0 ij ij        I WW I d (61)
Similarly, using the above-introduced definition for () ij D , this equation can be rewritten in a more compact form

() 61 ij T ij   D0 WW (62)
Hence, the rigid link produces 12 scalar equations describing the displacement and force equilibriums

() 6 6 6 1
() 6 6 6 1 ij i j ij T i j                   DI D t W W I0 t 0 (63)
which are written with respect to 24 scalar variables contained in six-dimensional vectors , ij  tt and , ij WW, similar to the flexible link case (the rank deficiency of the matrix is also equal to 12 here).

Modeling of a rigid mobile platform

If the platform can be treated as a rigid body, it can be included in the global stiffness model by means of several virtual rigid links connecting the nodes , , ,... i j k of the manipulator leg clamping points and the virtual rigid node e corresponding to the manipulator end-effector reference point (see Figure 7a). From the displacement equilibrium, one can derive several equations similar to (59) that may be aggregated in a common matrix equation

    D W D W 0 W D W (65)
which is an extended form of (62). In a matrix form, it can be presented as

( ) ( ) ( ) 6

Modeling of a non-rigid mobile platform

If the platform cannot be treated as a rigid body, it should be approximated by means of an equivalent structure consisting of several virtual flexible links connecting the nodes , , ,... i j k of the manipulator leg clamping and the virtual rigid node e corresponding to the manipulator end-effector reference point (see Figure 7b). These links are described by 12 12  stiffness matrix of the type (54) that should be aggregated in a common matrix equation

                                    K 0 0 K Wt Wt 0 K 0 K Wt 0 0 K K Wt K K K K K K (67)
where ,,

i j k
   t t t are the deflections at the leg clamping points, e t is the deflection at the end-effector reference point, ,

i j k
W W W are the wrenches at the leg clamping points, e W is the total wrench applied to the end-effector from the virtual links side. It is assumed here that all 66  stiffness matrices are presented in the global coordinate system. For this case of three attached legs, the model includes 24 scalar equations and 48 scalar variables, the rank deficiency of the model matrix is equal to 24, which is in a good agreement with four free-points of the considered elastic mechanical system.

MSA models of a manipulator's joints

In robotic manipulators, the robot base, legs/links and mobile platform are connected using either passive, actuated, elastic or rigid joints (Figure 8). Let us consider in detail their contributions to the global stiffness model of the manipulator.

Ac Ps

i j i j i j i j K e K e q q (a)

Modeling of a rigid joint

If adjacent links are connected by means of a rigid joint, the stiffness model must include relations describing two principal rules of structural mechanics: (a) displacement compatibility, (b) force equilibrium. Assuming that the nodes corresponding to the adjacent link denoted as , , ,... i j k , the first of these rules can be expressed as ...

i j k       t t t (68)
that allows the user to merge corresponding columns in the global stiffness matrix at the assembling stage. It is clear that the above equations can be also included in the stiffness model as an additional constraint without elimination redundant variables. In a matrix form, this constraint can be expressed as

  6 6 6 6 6 1 6 12 i j             t 0 t II (69)
in the case of two adjacent links with the rigid connection of the nodes , ij , and as

i j k                        I 0 I 0 t t 0I t I (70)
in the case of three adjacent links with the rigid connection of the nodes ,, i j k . In the case of more adjacent links, the constraint is modified by a simple extension of the relevant matrix. The second rule (the force equilibrium) produce the following relation ...

i j k     W W W 0 (71)
that directly follows from Newton's third law and allows the user to merge corresponding rows in the global stiffness matrix at the assembling stage. This relation can be also included in the stiffness model as an additional constraint without elimination redundant variables. In the case of two adjacent links with the rigid connection of the nodes , ij , the force equilibrium for the rigid connection of two links can be expressed as follows  

6 6 6 6 6 1 6 12 i j         W I I 0 W (72)
Similarly, the equation can be written for the case of three and more adjacent links with a rigid connection. For example, for three adjacent links, it will have the form   Hence, the rigid joint imposes some additional constraints on the model variables i t and i W . A number of these constraints is equal to 12 for the case of two links, 18 for three links, etc.

Modeling of a passive joint

If adjacent links are connected by means of a passive joint, weaker constraints should be applied compared to the rigid joint case. In particular, for two adjacent links (ai), (jb) with common nodes , ij , the displacement compatibility condition of , ij  tt must be replaced by the matrix equation

  r ij i j      Λ t t 0 (74)
where r ij Λ is a rank-deficient matrix that defines directions for the passive joint that do not admit free relative motions. For instance, for the revolute passive joint directed along the global y-axis, the matrix r Λ is expressed as

(1,1,1,1, 0,1) r Ry diag  Λ (75)
For the spherical joint, it is presented in the form (1,1,1, 0, 0, 0)

r Rxyz diag  Λ (76)
In general case, to derive the matrix r Λ , let us define the orthonormal basis 1 2 6 ,. , .., u u u associated with the passive joint in such way that the unit vectors 1 ,..., r uu describe the directions of the rigid connection and the unit vector corresponding to the passive connection allowing free relative motions of the links. This definition allows us to present the matrix r Λ as  

6 6 6 6 1 1 1 ,..., , | ,.., r r T    Λ u u 0 0 ( 77
)
that corresponds to the set of linear equations

( ) 0, 1, ij T k kr        u t t (78)
Using the matrix r ij Λ , the deflection compatibility constraint may be directly written in the form 61 6 12

rr ij i j j i             t 0 t ΛΛ (79)
which may be further reduced to a full-rank matrix equation

1 12 ** rr i r r j ij ij             Λ 0 t Λ t (80)
that contains a rectangular matrix * r ij Λ of the rank r . In the simple cases, the matrix * r ij Λ can be easily derived from r ij Λ by simple elimination of zero rows (see expression [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF]). In the general case, the reduced matrix *

r ij Λ is formed as follows   1 * 6 ,..., r r T r   Λ u u (81)
The force equilibrium condition for the passive joint should be also presented in the weaker form as

  r ij i j    Λ W W 0 (82)
that after elimination zero lines gives the second constraint

1 12 ** rr ij ij i r r j         W 0 W ΛΛ (83)
In addition, it is necessary to take into account that passive joints do not transmit the force/torque in the direction corresponding to zero columns in the matrix r ij Λ , which yields the following equations ;

pp ij i ij j     Λ W 0 Λ W 0 (84)
Here, the matrix p ij Λ defines the passive joint axes and contains the remaining unit vectors of the basis

Hence, for the passive joints, the force equilibrium produces two types of constraints ( 83), (84) that may be aggregated in the following matrix equation It should be noted that the above equations deal with the case when the passive joint connects two links only. In the case when more than two links are connected by a single passive joint, similar equations should be derived for each pair of the links, i.e. all passive connections should be treated separately (in contrast to the multiple rigid connections presented in the Section). It is also worth mentioning that the passive joint imposes the same number of constraints on the model variables i t and i W as in the case of the rigid joint, but their physical nature is different. In particular, the trivial (one dof) passive connection of two links produces 12 constraints, where 5 ones are imposed on i t and 7 others are imposed on i W (in contrast to the rigid joint case, where there are 6 constraints on i t and 6 constraints on i W ).

Modeling of an elastic joint

In some manipulators, the adjacent links are connected by means of an elastic joint that can be also treated as passive compliant joint with springs. In this case, the deflection compatibility condition (80) remains the same as above, it ensures equality for r components of deflection vectors , ij  tt . However, the force equilibrium condition must be slightly revised to take into account that some wrench components may be expressed via the Hooke's law:

                                 t II t W KK Λ W 0 0 0 Λ Λ Λ 0 W (92)
Hence, the elastic joint may be included in the global stiffness model similarly to the passive one by minor modifications of the corresponding matrix constraints.

Modeling of an actuated joint

The actuated joint ensures transmission of the force/torque between the manipulator links. In the frame of MSA technique, it can be presented either as a rigid, passive or elastic connection. In the first case, the actuating effort is transmitted in all direction 6 1 ,..., uu . In the second case, it is transmitted in the direction corresponding to the vectors 1 ,..., r uu while in the last case transmission is performed in the direction defined by the vectors 1 6 ,..., r  uu (see previous sections for details). Hence, the actuated joint may be included in the global stiffness model as a set of two linear matrix constraints describing the deflection compatibility and force equilibrium, which slightly differ depending on the accepted joint idealization. It is also worth mentioning that in the frame of the global stiffness model the actuated joint influences on the manipulator configuration, i.e. on the relative position and orientation of the links, which is obtained from the static equilibrium condition. Depending on the joint control type (position or torque based control), the actuator influence in the above-presented constraints is included either via the relation between , ij  tt or via the preloading 0 ij W .

Including of boundary conditions and external forces

To complete the global stiffness model, it is necessary to take into account specificity of the manipulator links connections to the robot base as well as the influence of the external forces (including reactions at the supports). Let us consider in detail contributions to the global stiffness model different types of the supports that may be implemented using either fixed, passive or elastic connections (Figure 10). In the frame of the MSA method, the equations produced by the supports are treated as the boundary conditions.

Base Ps Link

Base Link Base Link

0 1 0 1 0 1 K e q (a

Modeling of a rigid support

The rigid connection of the link (jb) to the robot base can be presented as a special case of the rigid joint with an eliminated link (ai) and zero deflection j  t0 . This simplifies the deflection compatibility constraint (69) down to

  6 6 6 1 j       I0 t (93)
that contains 6 linear equations to be included in the global stiffness model. Corresponding reaction wrench at the rigid support may be computed as

() 12
g jb jb   W K t ( 94
)
where the deflection b t is obtained from the solution of the global stiffness equations.

Modeling of a passive support

The passive connection of the link (jb) to the robot base can be presented as a special case of the passive joint with an eliminated link (ai) and zero deflection in the non-passive directions * r ij j    Λ t 0 . This simplifies the deflection compatibility constraint (80) down to

1 * 6 jr r ij r          t 0 Λ ( 95
)
that contributes r linear equations to the global stiffness model. In addition, it is necessary to take into account that the wrench components for the passive direction are equal to zero, which leads to the simplified form of (84)

1 * p ij j p  Λ W 0 (96)
that contributes another p equations to the global stiffness model. Totally, this gives 6 rp  independent linear constraints. Corresponding reaction wrench at the passive support may be computed as ( ) ( ) [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF] 12

g jb g jb j j b       W K t K t ( 97
)
where the deflections , jb  tt are obtained from the solution of the global stiffness equations.

Modeling of an elastic support

Similarly, the elastic connection of the link (jb) to the robot base can be presented as a special case of the elastic joint with an eliminated link (ai) and zero deflection in the non-elastic directions *

r ij j    Λ t 0
. This allows us to use the same deflection compatibility constraint (95) as above, which contributes r linear equations to the global stiffness model. In addition, it is necessary to take into account that the wrench components corresponding to the non-rigid directions are produced by the elastic forces satisfying the Hooke's law, i.e. that contributes another e equations to the global stiffness model. Totally, this gives 6 re  independent linear constraints. To compute the reaction wrench at the elastic support the same expression (97) can be used.

Including an external loading

In robotic manipulators, both serial and parallel, there is at least one node that is not connected directly to the robot base. It corresponds to the end-effector that interacts with robot environment by applying the force/torque to the external objects. For the global stiffness model, the end-effector produces the boundary conditions that should include the vector of the external wrench W , which is also necessary for computation of the Cartesian stiffness matrix.

To take into account the external loading ext W , the global stiffness model must be completed by the linear constraint derived from the force equilibrium at the node e, i.e. It is clear that this technique should be applied to all nodes with applied external wrenches.

Aggregation of MSA model components

In contrast to large mechanical structures consisting of a huge number of flexible components, a robotic manipulator is rather simple for the MSA analysis. Since the number of flexible links in manipulator is relatively small, the assembling stage can be simplified and the columns/rows merging operations can be avoided and replaced by adding relevant constraints [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]. Let us present this technique in detail.

As it was shown above, the MSA equations for robotic manipulator are derived from three main sources: (i) link models, (ii) joint models and (iii) boundary conditions. The first of them describes the force-displacement relations for all links (both flexible and rigid) yielding 12 scalar equations per link relating 24 variables { , , , }

i j i j
 t t W W . The second group of equations ensures the displacement compatibility and force/torque equilibrium for each internal connection (both rigid, passive and elastic), it includes two types of relations written for { , ,...} ij  tt and { , ,...} ij WW separately. Independent of the connection type, each joint provides 12 scalar equations if it connects two links, 18 scalar equations for the connection of three links, etc. The third group of equations is issued from the manipulator connections to the environment, they are presented as the force/displacement constraints for certain nodes that give up to 6 scalar equations for { , ,...} ij  tt or at least 6 scalar equations for { , ,...} ij WW depending on connection type (rigid, elastic or passive). In the general form, the aggregated stiffness model can be presented as follows [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF] (1)

(1) W0 (2) (2) ( t 2 W Δ 0 Δt ) {} {} ... ... ... i i                          t A A b W A A b (102)
where each line is issued from the equations/constraints presented in Sections 3.1-3.3, the matrices ()

W i A and ) Δt (i
A correspond to the terms with i W and i t respectively, the vectors () 0 i b represent the right-hand side of the corresponding equation/constraint. In this model, the contribution of the links can be expressed in the following way  flexible links:

12 12 12 12 {} { } { } {} ij i i          W IK t 0 (103)
 rigid links:

6 66 6 6 6 6 6 { , } { } {} { , } i iT i i j j                   0D 0 0 t IW I D 0 (104)
which gives 12 scalar equations per link (both flexible and rigid). The platform is described by one of the equations for flexible platform:

6 6 6 6 66 { } { } { } {} {} ij ij ex i i t                I K W 0 W 0K t (105)
for rigid platform:

  () 6 6 6 () 6 6 6 { , } {} {} ei i ei T i ext                 0I W 0 D W D 0 t (106)
whose dimension depends on the number of the legs. The joints contribute to the aggregated model the following constraints  rigid joints:

66 0 66 66 66 { , } { } {} { , } { } j i i i                      0 0 t I I W W I I 0 (107)
 passive joints:

** * * * * 0 { , } {} { , } { } {} {} rr ij ij r r r i ij ij ij i p ij ij                 0 0 W 0W 0 Λ Λ t Λ 0 Λ ΛΛ (108)
 elastic joints:

** 66 * * * * 0 66 0 { , } {} { , } { } {} } { , { } rr ij ij i i e e e e ij ij ee ij ij ij i ij j ij                         ΛΛ t Λ Λ Λ Λ 0 0 W I I 0 W KK W (109)
giving 12 scalar equations per internal joint (both rigid, passive and elastic). It should be noted that the external "joints" (connected to the base or end-effector) produce 6 scalar equations per connection only. So, their contribution to the aggregated model can be presented as

 rigid supports:   6 6 {} {} {} i i        W 0I t 0 (110)
 passive supports:

* * {} {} {} {} r ij i p i ij              0 W 0 0 0 Λ t Λ (111)
 elastic supports:

0 * ** {} {} {} {} { } { } r ij i ee i e ij ij ij ij              Λ t Λ 0 W K Λ 0 W (112)
And finally, the external loadings applied to some manipulator nodes are included in the model by means of the equations  external loadings:

  6 6 {} { } { } {} i i ext i        W I 0 W t (113)
It is worth mentioning that in practice the external loading applied to the manipulator end-effector is in the focus of stiffness analysis only.

For computational convenience, the equations ( 103)-( 113) can be arranged in the aggregated linear matrix equation of the following structure agr agr agr agr agr agr agr ext

                  S K W EF W b t (114)
where the matrices agr S , agr K , agr E , agr F , agr b , ext W are generated using relevant link/joint models or boundary conditions, while the vectors agr W and agr t contain all variables describing the wrenches and displacements, respectively. It can be proved that usually the system (114) cannot be solved straightforwardly because of the rank-deficiency of the matrix to be inverted. Nevertheless, this system allows us to obtain the Cartesian stiffness matrix defining the relation between the end-effector deflection and applied force/torque (which obviously always exists but can be singular). To find the desired Cartesian stiffness matrix, let us divide the node displacement variables agr t into two groups m t and e t corresponding to the manipulator internal nodes and the end-effector node, where the external wrench e W is applied. The latter allows us to rewrite the system (114) in the form

                              S K K W b E F C W W t ED t F (115)
and further, present it as

ee                 μb AB C t W D (116)
where all internal variables are included in the vector ( , )

                   S K b K A B C E F b C E F W (117)
Using the obtained system, the relation between the displacement e t and applied external wrench e W can be written as

0 e e e C    W K t W (118)
that includes the desired Cartesian stiffness matrix

1 C      KC A DB (119)
and the constant component

1 0 e     A W C b (120)
depending on the manipulator loadings and preloading included in the vector b .

It is worth mentioning that the matrix inversion in 1 

A usually exists (if the manipulator does not include redundant passive joints), otherwise the matrix inversion should be replaced by pseudo-inverse that distribute deflections between the redundant joints. Besides, the Cartesian stiffness matrix C K may be rank-deficient even for the non-singular matrix A , if the examined manipulator is under-constrained. Another important comment is related to the manipulator equilibrium configuration corresponding to given loadings and preloadings, which is implicitly assumed to be known here (i.e. computed before applying the MSA). However, for robotic manipulators, it is often not a trivial task [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF]. It should be also mentioned that the final expression (119) presents the manipulator stiffness model in the form of the classical Cartesian stiffness matrix C K of a size 6 6  corresponding to the fixed base and free end-effector. However, it can be easily transformed into the 'free-free' form that operates with 12 12  matrices using expressions similar to ( 9) and ( 13) taking into account that C K corresponds to 22 K .

Illustrative examples: MSA-based analysis of NAVARO manipulator

To demonstrate the utility of the MSA-based technique for robotics, let us apply it to the stiffness analysis of the NaVaRo robot (Figure 11a), a three-degree-of-freedom planar parallel manipulator with variable actuation schemes. It is composed of three identical legs and a moving platform formed of three segments rigidly linked at the central point. Each leg consists of four non-rigid links connected by five revolute joints to create a parallelogram linkage. Among them, there are four passive joints and one actuated joint connected to the motor via a double-clutching mechanism allowing to actuate one of two links adjacent to the motor axis (Figure 11b). The latter determines the main particularity of the NaVARo robot that has eight actuation modes that are switched from one to another in order to avoid kinematic singularities. To show the application of the MSA technique step-by-step, let us split the manipulator mechanism into four parts: three kinematically identical legs and a mobile platform. For convenience, let us divide the longest link of the leg into two rigidly connected parts. This allows us to present each leg (Figure 12) as a set of five flexible links (1,2), [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF], [START_REF] Yi | Geometric analysis of antagonistic stiffness in redundantly actuated parallel mechanisms[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF], [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Wang | Finite Element Analysis and Comparison of Two Hybrid Robots-the Tricept and the TriVariant, Intelligent Robots and Systems[END_REF], (9,e) with passive connections of the nodes <2,3>, <4,5>, <6,7> and rigid connection of the nodes <6,9>. Besides, depending on the actuation mode, the boundary conditions induced by the motor/clutch mechanics can be specified as the passive connection of the nodes <0,8> and elastic connection of the nodes <0,1>, or vice versa. Using the notation adopted in this paper, the "free-free" stiffness models of the links (54) can be presented in the form   ( , ) (1, 2), [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF], [START_REF] Yi | Geometric analysis of antagonistic stiffness in redundantly actuated parallel mechanisms[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF], [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Wang | Finite Element Analysis and Comparison of Two Hybrid Robots-the Tricept and the TriVariant, Intelligent Robots and Systems[END_REF], (9, )

ij ij ii ij ij jj i j e                     Wt KK Wt KK (121)
where all matrices are assumed to be already presented in the global coordinate system (here and below, the left superscript 'g' denoting the global system is omitted to improve readability). After assembling (121) into a single matrix one can get the following aggregated stiffness matrix Further, the stiffness model of the leg should be complemented by the second component that includes equation issued from the inter-link connections. For the considered manipulator, the links are connected by means of the rigid joint <6,9>, passive joints <2,3>, <4,5>, <6,7>, and elastic revolute joint <0,1>. The passive and elastic joints provide rotation around z-axis, which corresponds to the constraint matrices of the following form   * * * 1 0 0 0 0 0 0 1 0 0 0 0 ; 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

r p e           Λ Λ Λ (125)
Using this notation, the constraints imposed by the passive joints <2,3>, <4,5> can be presented as  

i i j j j ij                                             t W W Λ0 Λ Λ Λ Λ 00 t W W 0 Λ 0 (126)
The rigid joint constraint <6,9>, which ensures stiff connection of the links [START_REF] Yi | Geometric analysis of antagonistic stiffness in redundantly actuated parallel mechanisms[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF]

                                                  tt Λ 0 I I 0 W Λ tt W ΛΛ W 0 0 Λ Λ W 0 W ΛΛ W (127)
The third model component includes the boundary conditions. For the adopted actuation mode when the motor is connected to the link (1,2) via an elastic transmission and the link [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Wang | Finite Element Analysis and Comparison of Two Hybrid Robots-the Tricept and the TriVariant, Intelligent Robots and Systems[END_REF] is passively connected to the base, the boundary condition at the node <0,1,8> is presented in the form where 0 t is assumed to be equal to zero. The external loading can be introduced into the system using the following equation

e ext  WW (129)
After assembling all equations describing the constraints and boundary conditions the aggregated stiffness model of a single manipulator leg is presented as follows , which allow us to compute the desired stiffness matrices for all manipulator legs that will be further denoted as (1) (2) CC () C 3 ,, K K K . It should be also noted that here the Cartesian stiffness matrices of the legs are rank-deficient (their ranks are equal to 5), i.e. the presented enhanced MSA technique is able to deal with under-constrained structures. Nevertheless, assembling of three legs and the mobile platform in the NaVaRo robot provides the full-rank stiffness matrix (for non-singular configurations). To find the Cartesian stiffness matrix of the entire manipulator, let us consider the moving platform as a mechanical structure composed of three flexible links, which are rigidly connected to the end-effector on the right-hand side and to the manipulator's legs via a passive joint on the left-hand side. The stiffness model of this platform is described by the equation (67), that should be presented as [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF] (

                              KI A0 KS B 0 DC (
   KK I K K W 0 KK W 0 K K K K t (138)
where the nodes connected to the legs are denoted as <1>, <2>, <3>, the end-effector node is denoted as <e>, and

1 2 3
( , , )

agr col  W W W W , 1 2 3 ( , , ,
)

agr e col       t t t t t .
Relevant boundary conditions take into account connections to the flexible legs described by the stiffness matrices (1) (2) CC () C 3 ,, K K K via the passive joints with the matrices specified in (125). Using adopted notation, they are written as

  () ** C * ;0 1, 2,3 i i i r r i p i    Λ W Λ K t Λ W (139)
that can be presented in the extended matrix form as

 F K K K K (144)
For computational convenience, this system can be also converted into the form (114), where  b0 , which allow us to compute the desired stiffness matrix for the entire manipulator. Numerical results related to this example are presented in Tables 1-3, which contain the stiffness matrices of the manipulator and its legs corresponding to 3 different locations of the moving platform. The stiffness parameters of the manipulator links were estimated by means of FEA modeling in the CAD environment CATIA using the technique proposed in [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF]. Corresponding matrices are presented in Appendix. As follows from the obtained results, the proposed generalization of the MSA allows the user to obtain both full-rank and rank-deficient Cartesian stiffness matrices for complex architectures containing numerous passive, rigid and elastic joints. In particular, the stiffness matrices of the individual legs are rank-deficient for all considered manipulator configurations. In contrast, the stiffness matrices of the entire manipulator are non-singular in all examined cases. Moreover, these matrices are diagonal for the configurations 1 and 2, which are in good agreement with the manipulator postures symmetry in these locations. Hence, the presented example confirms the advantages of the developed techniques, which enlarges the MSA application area from conventional truss structures to complex robotic mechanisms.

                                   K (1) ( ) 5 C rank  K

Manipulator configuration

Il manqué la direction des axes x et y Couleur des jambes 1, 2, 3

Leg #2

                             K (2) ( ) 5 C rank  K Leg #3
                               K (3)
                                   K ( 1
C                                   K (2) ( ) 5 C rank  K Leg #3 2 4 3 4 5 4 5 4 4 4 3 3 4 (3) 33 3 4 3 6
                               K (3)
                               K (1) ( ) 5 C rank  K Manipulator configuration Leg #2
                                   K (2) ( ) 5 C rank  K Leg #3
                                  K (3) ( ) 5 C rank  K ( ) 6 C rank  K

Discussion

As follows from the above presented example, the proposed enhancement of the MSA method, which was originally designed for analysis of truss structures, allows user effectively analysing stiffness properties of serial and parallel robotic manipulator (i.e. with open and closed loop kinematics) that are composed of both flexible and rigid links connected by actuated, passive, elastic or rigid joints. This approach can be applied to under-constrained, fully-constrained and over-constrained structures both in non-singular and singular configurations, it is able to produce either full-rank or rank-deficient Cartesian stiffness matrices describing the manipulator stiffness properties. It presents the manipulator stiffness model as a set of two groups of matrix equations describing elasticity of separate links and connections between the links in the form of constraints. Its principal advantage is the simplicity of the model generation that includes straightforward aggregation of link/joint equations without conventional merging of rows and columns in the global stiffness matrix. Nevertheless, despite numerous advantages, the proposed technique has some limitations that should be discussed in detail. The most evident difficulty is related to the numerical precision of high-dimensional matrix inversion required for the computing of the desired Cartesian stiffness matrix. In fact, the matrix to be inverted contains both small and large non-homogeneous elements (with different dimensions [ ], [ ],[ ], ], [ ...[ / ] m rad N N m Nm ) whose magnitude essentially depends on the adopted units (m or mm, rad or deg, etc.). For example, the stiffness parameters of the NaVaRo manipulator links differs more than six orders of magnitude. The corresponding matrix of the aggregated stiffness model for a single leg (of size 114 114  ) has the condition number about 12 10 which requires special efforts in the inversion procedure (such as normalization, scaling, etc.) to ensure required precision. It is clear that for the entire manipulator, where the matrix size is 396 396  , this ill-conditioning problem is more severe. Nevertheless, modern linear algebra routines are able to handle even more difficult problems (arising in the FEA modeling, for instance). So, the above-mentioned difficulty is not critical here and does not affect essentially the final results. Another difficulty appears while preparing reliable input data for the model, which include the link geometric and stiffness parameters. As follows from our experience, there are two possible ways to tackle this problem. The simplest one is based on the approximation of the link by simple shapes (see Table 1) and using standard expressions for the stiffness matrix elements. Otherwise, the CAD-based technique [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF] is quite efficient provided that the link stiffness matrices are properly filtered from small elements that theoretically should be equal to zeros. Besides, the output result containing the desired Cartesian stiffness matrix should be also filtered by elimination non-significant elements. It should be also mentioned that in this paper we have not addressed the problem of the static equilibrium computing, which is implicitly assumed to be known while transforming the link stiffness matrices to the global coordinate system. To solve this problem, either a nominal extended geometric model of the manipulator can be used (providing an analytical solution) or an iterative technique proposed in our previous work [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF] can be adapted. On the other hand, further development of the static equilibrium computing methods and application of the proposed MSA-based technique allow the user to build a non-linear stiffness model, which will be able to detect such important phenomena as buckling effects.

Conclusions

The paper presents fundamentals of the MSA-based manipulator stiffness modeling. It proposes a novel unified and systematic approach that is suitable for serial, parallel and hybrid architectures containing closed-loops, flexible links, and rigid connections, passive and elastic joints, flexible and rigid platforms. This approach can be applied to both under-constrained, fully-constrained and over-constrained manipulators, which are subject to numerous external and internal loadings. It can produce either full-rank or rank-deficient Cartesian stiffness matrices corresponding to both generic and singular manipulator configurations. The developed approach is based on a unified mathematical formulation that presents the manipulator stiffness model as a set of two groups of matrix equations describing elasticity of separate links and connections between the links (in the form of constraints). Its principal advantage is the simplicity of the model generation that includes straightforward aggregation of link/joint equations without conventional merging of rows and columns in the global stiffness matrix. Compared to the previous works, which concentrate either on the classical MSA applied to truss systems or present some MSA extensions for the manipulator stiffness modelling via potential energy minimization, this work proposes a systematic approach allowing to describe in the same way both flexible and rigid mechanical elements and easily take into account passive connections, without using the Lagrange constraint optimisation technique. The developed technique has been successfully applied to the stiffness modeling of the planar parallel manipulator NaVaRo with three legs, which consists of 13 flexible links connected by 15 passive joints. The stiffness analysis was performed both for the manipulator legs separately (under-actuated quasi-serial architecture with a kinematic parallelogram) and entire parallel manipulator (parallel architecture with six closed-loops). The latter allowed to proof the method capability to produce both singular and full-rank Cartesian stiffness matrices for the complex architecture with numerous closed loops and passive joints. In future, this technique will be used for parametric optimization of the NaVaRo manipulator and applied to more complex architectures.

Acknowledgments

The work presented in this paper was supported by the grant of Russian Science Foundation №17-19-01740 and the ANR (Project ANR-2010-SEGI-003-02-COROUSSO), France.

Appendix: Compliance matrices of the NAVARO manipulator links

The stiffness parameters of the NAVARO manipulator links were estimated using CAD-based technique proposed in our previous work [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF]. To find them, the virtual forces/torques were applied to all links separately and corresponding deflections were obtained via the FEA modeling. These data were used for the identification of the link compliance coefficients describing the relation between the external force/torque and the linear/angular deflections. The obtained values of the compliance matrices are given in Table 5, where the ()

Figure 1

 1 Figure 1 Cantilever and unsupported beams: definition of frames, deflections and wrenches

(Figure 3

 3 Figure 3 Two-link serial system with rigid connection

Figure 6

 6 Figure 6 Two-link parallel system with passive joints Similar to all previous cases the basic elasto-static equations and are written as follows

  this example where the number of the leg clamping points is assumed to be equal to three, the model includes 24 scalar equations and 48 scalar variables. The rank deficiency of the model matrix is equal to 24, which agrees with the rigid body mechanics. It is clear that this model can be easily generalized to the case of 4, 5, … of attached legs by straightforward expanding the matrix sizes.

K 3 Figure 7

 37 Figure 7 Presentation of rigid (a) and non-rigid (b) platforms for the stiffness model of the parallel manipulator

Figure 8

 8 Figure 8 Connections between elements of a robotic manipulator

Figure 9

 9 Figure 9 Rigid linking of several joints

Figure 10

 10 Figure 10Connections between the manipulator links and base.

  is possible to derive the constraints for the case of three and more adjacent links, for instance

Figure 11

 11 Figure 11 Parallel robot NaVaRo and its geometry + axes et couleurs des jambs significations

Figure 12

 12 Figure 12 MSA-based representation of the NaVaRo leg (actuation of the link (1,2))

  137) Further, after separating the node variables {} i t in two groups corresponding to the internal nodes m t and to the end effector node e t one can get the linear systems (115) and (

22 ijC

 22 corresponds to the inverse of the stiffness matrix () 22 ij K in basic eq. (54). The stiffness coefficient of the actuator was estimated

  one can get the following linear system

		rr
				66
	p	1
		22	6 6	6 6	1
			22 qq pp
		6 6	11	12	2
	r	1 22	2 11 qq 2 12 rr	3 18 1
	3 30 1	6 6	22 21 22 q q	1 30 8

  It should be stressed that the above expressions are valid for the so-called "non-preloaded case", when equal deflections do not generate the elastic forces, i.e. ** . In the case when the springs are initially preloaded by the wrench

	6 6 * e ij ee 6 6 i j ij    6 6 * * 6 (6 66 e e ij i je       00 Λ Λ Λ II KK 0	) 24  e	24 1      t i j i j      W t W	0	(6 ) 1 e 	(90)
													ee ij i i   Λ W Λ W 0 0 j j   ij W , the
	force equilibrium equation must be presented in the form
	0 ** ; ) ( j i   W W 0 i e e e ij ij ij ij     () i j   Λ K Λ t t WW	(91)
	that leads to the following linear matrix constraint
															i
			6 6				6 6				6 66	6	j	6 1
			e	e				e	e			e		e	0
			ij	* ij			ij	* ij		* ij	e	6 (6 ) 2 4 ee * (6 ) 1 ij ij i
															j	24 1
	** ; i e g e e ij ij ij  W W 0 j i  () ij    Λ K Λ t t W	(88)
	where the matrix * e ij Λ of size	6 e ,		e	6  corresponds to the non-rigid directions of the joint (similar to * r p ij Λ in the passive joint
	case)
	Λ	* e			u	r	1 	,...,	u	6	 e T 	6	,		(89)
	and ge														as
	the following matrix constraint

ij

K is e e  stiffness matrix describing elastic properties of the joint. So, the force equilibrium condition can be presented

Table 1 .

 1 Stiffness matrices of NAVARO manipulator for platform location #1 (

	0, xy  0, 	0  o ) (unités)

Table 2 .

 2 Stiffness matrices of NAVARO manipulator for platform location #2 (

	K	C	3.75 ~0 0 0 0 10  ~0 0 0 0 ~0 3.75 0 0 0 ~0 0 2.92 ~0 ~0 0 0 ~0 4.46 ~0 0 0 ~~0 4.46 0 10 10 10 0 10            
			~0   ~0 0 0 0 8. 10 00 
			rank	( ) 6 C  K

Table 3 .

 3 Stiffness

						matrices of NAVARO manipulator for platform location #3 (	150, xy  150, 	0  o )
										Leg #1
		1.06	10	6	5.39	10	5	0				0	0	5.81	10	4
		5.39	10	5	2.74	10	5	0				0	0	2.96	10	4
	) (1	0				0				2.24	5 10	2.92	10	4	2.78	10	4	0
	C	0				0				2.92	10	4	5.43	10	3	4.69	10	3	0
		0				0				2.7 8	10	4	4.69	10	33 5.18 10	0
		5.81	10	43 4 2.96 0 1 0 0 3. 01 8 0 0

Corresponding author. Tel. +7 904 718 58 00; E-mail address: a.klimchik@innopolis.ru (A. Klimchik).

Link (i,e)