Anatol Pashkevich
email: anatol.pashkevich@emn.fr

Alexandr Klimchik

Damien Chablat

Enhanced stiffness modeling of manipulators with passive joints

Keywords: Stiffness, parallel mechanisms, external and internal loading, equilibrium configuration, stability

The paper presents a methodology to enhance the stiffness analysis of serial and parallel manipulators with passive joints. It directly takes into account the loading influence on the manipulator configuration and, consequently, on its Jacobians and Hessians. The main contributions of this paper are the introduction of a non-linear stiffness model for the manipulators with passive joints, a relevant numerical technique for its linearization and computing of the Cartesian stiffness matrix which allows rank-deficiency. Within the developed technique, the manipulator elements are presented as pseudo-rigid bodies separated by multidimensional virtual springs and perfect passive joints. Simulation examples are presented that deal with parallel manipulators of the Ortholide family and demonstrate the ability of the developed methodology to describe nonlinear behavior of the manipulator structure such as a sudden change of the elastic instability properties (buckling).

Introduction

Mechanical stiffness 1 of a manipulator is one of the most important indicators in performance evaluation of robotic systems [START_REF] Park | Safe link mechanism based on nonlinear stiffness for collision safety[END_REF][START_REF] Angeles | Performance Evaluation and Design Criteria[END_REF][START_REF] De Luca | Robots with Flexible Elements[END_REF]. In particular, for industrial robots where the primary target is the precise manipulation of a technological tool, the manipulator stiffness defines the positioning errors due to the external loading arising during the workpiece processing. Similarly, in industrial pick-and-place applications which are intended for simple but fast manipulations, the stiffness defines admissible velocity/acceleration while approaching the target point, in order to avoid undesirable displacements due to inertia forces [START_REF]Handbook of industrial robotics[END_REF]. Other examples include large robotic manipulators for a patient positioning in medical treatment, where elastic deformations of mechanical components under the task load (and under own link weight) is the primary source of positioning errors [START_REF] Meggiolaro | Geometric and elastic error calibration of a high accuracy patient positioning system[END_REF]. It is obvious that in all of these cases, the desired stiffness should be high enough to meet the requirements of the relevant application.

In contrast, for service robots, which interact directly with humans, a rather low stiffness is required to eliminate collisions causing an operator injury. This leads to special "soft" (bionic or biologically inspired) manipulator architectures that are based on utilization of numerous elastic links, joints and actuators [START_REF] Filippini | Toward soft robots you can depend on[END_REF][START_REF] Scarfogliero | The use of compliant joints and elastic energy storage in bio-inspired legged robots[END_REF]. Some more recent examples incorporating elastic elements, are wire-driven (cable-based) robots that are potentially interesting for medical rehabilitation and rescue operations [START_REF] Merlet | Analysis of Wire Elasticity for Wire-driven Parallel Robots[END_REF]. Another active research area that needs sophisticated stiffness analysis is micromanipulation where the motions are generated due to piezoelectric actuators but conventional passive joints are replaced with elastic hinges [START_REF] Kang | Analysis and design of parallel mechanisms with flexure joints[END_REF][START_REF] Yong | Kinetostatic modeling of 3-RRR compliant micro-motion stages with flexure hinges[END_REF][START_REF] Lobontiu | Analytical model of displacement amplification and stiffness optimization for a class of flexure-based compliant mechanisms[END_REF].

Related works

At the preliminary design stage, the stiffness analysis can be performed using the virtual joint method (VJM) that describes all types of flexibility (both distributed and lumped) by localized virtual springs located in the manipulator joints. Then, for this compliant mechanism, it is computed the Cartesian stiffness matrix that depends on current configuration (posture) of the manipulator. Mathematical background for this computation originates from the work of Salisbury [START_REF] Salisbury | Active Stiffness Control of a Manipulator in Cartesian Coordinates[END_REF] who derived a closed-form expression for the Cartesian stiffness matrix of a serial manipulator assuming that the mechanical elasticity is concentrated in the actuated joints. Retaining this assumption, Gosselin [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF] extended this result for the case of parallel manipulators (where the links were assumed to be rigid, and the passive joints to be perfect). Further development of this approach allowed taking into account the elasticity of the links, which were presented as rigid beams supplemented by linear and torsional springs [START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF]. At present, there are a number of variations and simplifications of VJM, which differ in modeling assumptions and numerical techniques. In particular, it was applied to the CaPAMan, Orthoglide and H4 robots, specific variants of Stewart-Gough platform, manipulators with US/UPS legs and other kinematic architectures [START_REF] Pashkevich | Nonlinear effect in stiffness modeling of robotic manipulators[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF][START_REF] Ceccarelli | A stiffness analysis for CaPaMan (Cassino Parallel Manipulator)[END_REF][START_REF] Company | Modelling and preliminary design issues of a 4-axis parallel machine for heavy parts handling[END_REF][START_REF] Arumugam | Stiffness analysis of a class of parallel mechanisms for micro-positioning applications[END_REF][START_REF] Vertechy | Static and stiffness analyses of a class of over-constrained parallel manipulators with legs of type US and UPS[END_REF].

In the first works, the stiffness parameters of the virtual springs describing the link elasticity were evaluated rather approximately, using rough presentation of the link shape by rectangular beams. Besides, it was assumed that all linear and angular deflections (compression/tension, bending, torsion) are decoupled and are presented by independent onedimensional springs that produce a diagonal stiffness matrix of size 66 for each link. Afterwards, this elasticity model was enhanced by using complete 66 non-diagonal stiffness matrix of the cantilever beam that is known from structural mechanics. This allowed taking into account all types of the translational/rotational compliance and relevant coupling between different deflections. Further advance in this direction (applicable to the links of complicated shape) led to the FEA-based identification technique that involves virtual loading experiments in CAD environment and stiffness matrix estimation using dedicated numerical routines [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]. The latter essentially increased accuracy of the VJM-modeling while preserving its high computational efficiency. It is worth mentioning that usual high computational expenses of FEA is not a critical issue here, because it is applied only once for each link (in contrast to the straightforward FEA-modeling for the entire manipulator, which requires complete re-computing for each manipulator posture).

Another important research issue is related to passive joints which are widely used in parallel manipulators. In the simplest case, when number/type of geometrical constraints perfectly corresponds to number/type of passive joints, the redundant variables (passive joint coordinates) may be just eliminated from the kinematic model. This allows computing the conventional Jacobian and enables a direct application of the Salisbury formula. However, in the case of overconstrained or under-constrained manipulators, the elimination technique can not be used directly.

For serial kinematic chains with passive joints, the problem was solved for the general case [START_REF] Pashkevich | Stiffness analysis of 3-d.o.f. overconstrained translational parallel manipulators[END_REF]. In particular, it was proposed an algorithmic solution that extends the Salisbury formula and is able to produce the rank-deficient stiffness matrices describing "zero-resistance" of the end-effector to certain type of displacements, which do not require deflections in the virtual springs (due to presence of passive joints and/or kinematic singularity of the examined posture). Relevant technique involves inversion of dedicated square matrix of size ( 6) ( 6)

nn    
which is composed of the links stiffness matrices and kinematic Jacobians of both virtual springs and passive joints (here n  is the number of passive joints). Then, the desired Cartesian stiffness matrix is obtained by simple extraction of an appropriate 66  sub-matrix from the computed inverse. The main advantage of this method is its computational simplicity, since the number of the virtual springs do not influence on the size of the matrix to be inverted 3 . Besides, the method does not require manual elimination of the redundant spring corresponding to the passive joints, since this operation is inherently included in the numerical algorithm.

However, for parallel manipulators with passive joints, solutions were obtained for less general cases. They include "pure" parallel architectures where the base and the end platform are connected by strictly serial kinematic chains. Here, the total stiffness matrix can be presented as the sum of partial matrices corresponding to separate chains (computed using the above described technique), so the passive joints are taken into account easily. Besides, in this case the overconstraining of the mechanism does not create additional difficulties. For instance, for the over-constrained manipulator of Orthoglide family [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF], each of the parallel chains yields the stiffness matrix of rank 4 while their aggregation gives the matrix of full rank 6. But if there exists a cross-linking between the parallel chains (like in kinematic parallelograms, for instance), this method can not be applied directly. For this case, some interesting results are presented in [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF][START_REF] Quennouelle | Instantaneous Kinemato-Static Model of Planar Compliant Parallel Mechanisms[END_REF][START_REF] Quennouelle | Modelisation geometrico-statique des mecanismes paralleles compliants[END_REF] where the geometrical constraints were treated in a general way but detailed computational techniques were not developed.

The majority of related works implicitly assume that the stiffness is evaluated in a quasi-static configuration with no external or internal loading. There are very limited number of publications that directly address the loaded working mode (or so-called "large deflections" case), where in addition to the conventional "elastic stiffness" in the joints it is necessary to take into account the "geometrical stiffness" due to the change in the manipulator configuration under the load. Although the existence of this additional stiffness component for elastic structures has been known for a long time [START_REF] Przemieniecki | Theory of Matrix Structural Analysis[END_REF], the importance of this problem for robotic manipulators has been highlighted rather recently. The most essential results in this area were obtained in [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Sh | Geometrical Approach to The Conservative Congruence Transformation (CCT) for Robotic Stiffness Control[END_REF][START_REF] Li | Stiffness Control and Transformation for Robotic Systems With Coordinate and Non-Coordinate Bases[END_REF] where there are presented both some theoretical issues and several case studies for serial and parallel manipulators. Several authors [START_REF] Chakarov | Study of passive compliance of parallel manipulators[END_REF][START_REF] Chakarov | Study of antagonistic stiffness of parallel manipulators with actuation redundancy[END_REF][START_REF] Wei | Design of planar parallel robots with preloaded flexures for guaranteed backlash prevention[END_REF] addressed the problem of stiffness analysis for the manipulators with internal preloading or antagonistic actuating, but in relevant equations some of the second order kinematic derivates were neglected.

To out knowledge, the first detailed study of the second-order coefficients in stiffness analysis was done in [START_REF] Yi | Synthesis of actively adjustable springs by antagonistic redundant actuation[END_REF][START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF], who took into account external loads, gravity, active and passive stiffness in actuated and unactuated joints and also considered antagonistic redundant actuation. However, the derived model requires solution of non-linear matrix equation (where the joint stiffnesses and external/internal loadings are parameters), which includes a number of first-and secondorder matrix derivatives that obviously must be computed in the neighborhood of loaded equilibrium that also depends on the loadings. But this issue was not considered in details and, as result, any nonlinear effects were no detected in stiffness behavior of the examined manipulators.

Table 1

Summary of the related works and expressions for the Cartesian stiffness matrix

Publications

Model & assumptions Stiffness matrix

Salisbury [START_REF] Salisbury | Active Stiffness Control of a Manipulator in Cartesian Coordinates[END_REF] Serial manipulator, elasticity in actuators [START_REF] Ciblak | Synthesis of stiffness by springs[END_REF][START_REF] Ciblak | Synthesis of Cartesian Stiffness for Robotic Applications[END_REF],

1 T c      K J K J Gosselin [41], Ciblak & Lipkin
Pigoski et al., [START_REF] Pigoski | Stiffness mappings employing different frames of reference[END_REF] Parallel manipulator, elasticity in actuators, non over-constrained

1 T c      K J K J
Zhang et al. [START_REF] Zhang | Analysis of parallel kinematic machine with kinetostatic modeling method[END_REF][START_REF] Xi | Global kinetostatic modelling of tripod-based parallel kinematic machine[END_REF] Serial kinematic chain without passive joints, elasticity in virtual joints

1 1 T c i i i i           K J K J
Pashkevich et al. [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] Serial kinematic chain with passive joints, elasticity in virtual joints [START_REF] Zhang | Kinetostatic modeling of parallel mechanisms with a passive constraining leg and revolute actuators[END_REF] Parallel manipulator without cross-linking between kinematic chains

1 1 T q c T q                J K J J K J0 Zhang & Gosselin,
1 1 ; T c ci ci j j j ij           K K K J K J
Pashkevich et al., [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] Parallel manipulator without cross-linking between kinematic chains

1 1 ; T i i i qi ci c ci T i qi                 J K J J K KK J0
Alici & Shirinzadeh, [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF] Chen & Kao [START_REF] Chen | Conservative Congruence Transformation for Joint and Cartesian Stiffness Matrices of Robotic Hands and Fingers[END_REF][START_REF] Li | Stiffness Control and Transformation for Robotic Systems With Coordinate and Non-Coordinate Bases[END_REF],

Marlet & Gosselin [START_REF] Merlet | Parallel mechanisms and robots[END_REF] Serial or parallel manipulator with external loading (non over-constrained)

  1 T cF      K J K K J
Quennouelle & Gosselin [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF][START_REF] Quennouelle | Instantaneous Kinemato-Static Model of Planar Compliant Parallel Mechanisms[END_REF] Parallel manipulator with external loading and supplementary geometric constraints (cross-linkings)

  1 T c F I        K J K K K J
Yi & Freeman [START_REF] Yi | Synthesis of actively adjustable springs by antagonistic redundant actuation[END_REF][START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF] Parallel manipulator with external loading, inertia and gravity loads, joint stiffness, actuation redundancy Therefore, in most of the related works the problem of finding the loaded equilibrium was omitted, so the Jacobian and Hessian were computed in a traditional way, i.e. for the unloaded configuration. The latter yielded essential computational simplification but also imposed crucial limitations, not allowing detecting the buckling and other non-linear phenomena known from general theory of elastic structures. On the other side, these issues become more and more critical in robotics applications where the geometrical buckling of the manipulator structure is more likely than the buckling in separate links.

1 T cu    K J K J u K -
The above mentioned publications are briefly summarised in Table 1, where all notations are adopted to those used in this paper. As it follows from their analysis, the stiffness modeling in loaded working mode for both serial and parallel manipulators with passive joints needs further improvement. In particular, it is required to develop dedicated numerical techniques that are well adopted for robotic applications and provide a designer with capability of computing the loaded equilibriums (which may be non-unique), their stability analysis, local linearization of the force-deflection relation, and also evaluation of possible non-linear phenomena (such as buckling) that can be potentially dangerous in practical applications. This motivates enhancement of our previous results [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF][START_REF] Pashkevich | Stiffness Analysis of Parallel Manipulators with Preloaded Passive Joints[END_REF] and extending them for the case of the external loading.

3

Problem of stiffness modeling

Basic Assumptions

The examined class of manipulators that is in the focus of this paper includes serial and parallel architectures with flexible actuated joints and flexible links that may be separated by a number of passive joints. In the case of a parallel manipulator, its architecture can be decomposed into several kinematic chains that should be studied separately. Thus, the stiffness analysis of the serial chain with passive joints (taking into account external and internal loading) is a key issue of this study. It is obvious that such kinematic chains are statically under-constrained and their stiffness analysis can not be performed by direct application of the standard methods.

To describe most of existing parallel architectures, it is assumed that the general serial kinematic chain under study consists of a fixed "Base", a number of flexible actuated joints "Ac", a serial chain of flexible "Links", a number of passive joints "Ps" and a moving "Platform" at the end (Fig. 1). It is also assumed that all links are separated by the joints (actuated or passive, rotational or translational) and their order is arbitrary. Besides, it is admitted that some links may be separated by actuated and passive joints simultaneously. Typical examples of the examined kinematic chains can be found in 3-PUU translational parallel kinematic machine [START_REF] Li | Stiffness analysis for a 3-PUU parallel kinematic machine[END_REF], in Delta parallel robot [START_REF] Clavel | DELTA, a fast robot with parallel geometry[END_REF] or in parallel manipulators of the Orthoglide family [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF] and other manipulators [START_REF] Merlet | Parallel Robots[END_REF]. It is worth mentioning that here a specific spatial arrangement of under-constrained chains yields the over-constrained mechanism that posses a high structural rigidity with respect to the external force. To evaluate stiffness of the considered serial manipulator, let us apply a modification of the virtual joint method (VJM), which is based on the lump modeling approach [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF]. According to this approach, the original rigid model should be extended by adding virtual joints (localized springs), which describe elastic deformations of the links. Besides, virtual springs are included in the actuating joints, to take into account stiffness of the control loop. Detailed description of such presentation as well as related assumptions are given in our previous publications [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF], where the kinematic model is presented as a product of homogeneous transformation matrices that after extraction of position and orientation of endeffector is transformed into the vector function [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF] ( , )  t g q θ [START_REF] Park | Safe link mechanism based on nonlinear stiffness for collision safety[END_REF] Here the vector

., )

T n q q q  q contains "perfect passive joint" coordinates (i.e., 6 without internal preloading), the vector 12 ( , , ..., )

T m     θ
collects coordinates of all virtual joints and "preloaded passive joints" (with auxiliary internal springs); n , m are the sizes of q and θ respectively.

Physical nature of passive joints with internal preloading is illustrated in Fig. 2. Such joints include internal springs, so their statics is described by the following expression  

0 i i i i K        (2)
where i   is the torque/force caused by deviation of the joint coordinate i  from its non-loaded ("zero") value 0 i  , and coefficient i K  defines the spring stiffness. For the purpose of generality, let us introduce similar "zero" values 0 i  for the virtual springs that described flexibility of the links (obviously they are equal to zero for this subset of θ ). This allows to define vector 0 θ of the same size as θ and to present the static equations corresponding to all variables (corresponding to perfect and preloaded passive joints, virtual springs of links and actuators) in general form

  0 , q      τ K θ θ τ 0 .
Here , q  ττ are the generalized torque/force in joints corresponding to the variables θ and q , the matrix  K collects stiffness coefficients of all springs of the kinematic chain.

It worth mentioning that in the case without internal preloading (see for example [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF]) the vector θ describes only flexibility of manipulator links/actuators that are presented by virtual springs, while vector q collects entire set of passive joint coordinates. In contrast, in this paper, the passive joint coordinates are divided into two subsets: (i) "perfect passive joints" included in q , and (ii) "preloaded passive joints" included in θ together with traditional virtual springs. Besides, if a passive joint includes a nonlinear spring (see Fig. 2), the corresponding joint variable may be included either in θ or q , depending on current configuration of manipulator. However, for each configuration, this assignment is strictly unique.

Problem statement

The designed stiffness model describes the resistance of manipulator to deformations caused by an external force or torque. It can be defined by the relation

() f  F Δt , where (...)
f is a so-called "force-deflections" function that associates a deflection Δt with an external force F that causes the transition. It is worth mentioning that the function (...) f can de determined even for the singular configurations (or redundant kinematics) while the inverse statement is not generally true. Hence, enhanced stiffness analysis must include computing of this function and detailed analysis of its singularities that may provoke various nonlinear phenomena (such as buckling). In the unloaded case, this function is usually defined through the "stiffness matrix" K , which describes the linear relation 00 ( ,

)  F K q θ Δt (3)
between small six-dimensional translational/rotational displacements Δt , and the external forces/torques F causing this transition. Here, it is assumed that Δt includes thee positional components   ,, x y z    describing displacement in Cartesian space and three angular components   ,,

x y z      
that describe the end-platform rotation around the Cartesian axes, while the vectors 00 , q θ correspond to the manipulator equilibrium for which the loadings (both internal and external) are equal to zero. However, for the loaded mode, similar linear relation is defined in the neighborhood of another static equilibrium, which corresponds to a different manipulator configuration ( , ) q θ , that is modified by external forces/torques F . Respectively, in this case, the stiffness model describes the relation between the increments of the force δF and the position δt ( , ) [START_REF]Handbook of industrial robotics[END_REF] where  0 qq Δq and  0 θ θ Δθ denote the new configuration of the manipulator, and Δq , Δθ are the deviations in the coordinates , q θ respectively. Hence, the considered problem of enhanced stiffness analysis may be divided into several sequential subtasks: (i) non-linear stiffness modeling that includes computing full-scale "force-deflections relation" for the externally loaded manipulator (taking into account internal preloading) and checking stability of corresponding equilibrium configurations;

F  δF K q θ δt
(ii) linearization of the relevant force-deflection relations in the neighborhood of this equilibrium and computing corresponding stiffness matrix that in general case can be singular due to presence of passive joints;

(iii) determining the critical forces that may cause undesired buckling phenomena or sudden change of current configuration of the loaded manipulator.

Static equilibrium for the loaded mode

Computing of the static equilibrium is a key issue for the non-linear stiffness analysis, since it defines the manipulator configuration ( , ) q θ required for the linearization of the "load-deflection" relation. In previous works, this issue was usually ignored and the linearization was performed in the neighborhood of the unloaded configuration assuming that the external load is small enough. It is obvious that the latter essentially limits relevant results and does not allow detecting non-linear effects such as buckling. From mathematical point of view, the problem is reduced to solving of a system of non-linear static equilibrium equations that may produce unique or non-unique, stable or unstable solutions.

Configuration of loaded manipulator

For computational reasons, let us consider the dual problem that deals with determining the external force F and the manipulator equilibrium configuration ( , ) q θ that correspond to the end-effector location t taking into account internal preloading in the joints. Let us assume that the joints are given small, arbitrary virtual displacements ,  q θ in the neighborhood of ( , ) q θ . According to the principle of virtual displacements, the virtual work of the external force F applied to the end-effector along the corresponding displacement

q         tJ θ J q is equal to the sum     TT q         FJ θ F J
q (here  J and q J are the kinematic Jacobians with respect to the coordinates , q). Since the passive joints do not produce the force/torque reactions, the virtual work corresponding to the generalized forces/torques , q  ττ includes only one component

T    
τθ (the minus sign takes into account the force-displacement directions for the virtual spring). In the static equilibrium, the total virtual work of all forces is equal to zero for any virtual displacement, therefore the equilibrium conditions may be written as T   JF τ ; T q  J F 0 , and taking into account assumptions and notations from the previous section, the static equilibrium conditions can be presented as

  0 ; ; ( , ) TT q         J F K θ θ J F 0 t g q θ (5)
where the vector 0 θ defines internal preloading in the joints, the matrix  K describes stiffness of all springs in the adopted VJM model, while the external loading F and the configuration ( , ) q θ are treated as an unknown for given end-effector location ( , )  t g q θ . Hence, the designed equilibrium configuration must satisfy the system of nonlinear equations [START_REF] Meggiolaro | Geometric and elastic error calibration of a high accuracy patient positioning system[END_REF]. It is evident that there is no general method for analytical solution of this system and it is required to apply numerical techniques. To derive the numerical algorithm, let us linearize the kinematic equation in the neighborhood of the current position ( ,

) ii q θ 11 ( , ) ( , ) ( ) ( , ) ( ) i i q i i i i i i i i           t g q θ
J q θ q q J q θ θ θ [START_REF] Filippini | Toward soft robots you can depend on[END_REF] and rewrite the static equilibrium equations as

  1 1 0 1 ( , ) ; ( , ) TT i i i i q i i i           Jq θ F K θ θ J q θ F 0 , (7)
which leads to a system of linear algebraic equations with respect to 111 ( , , )

iii  q θF that includes the Jacobians ( , )

ii  Jqθ , ( , )

q i i

Jqθ and geometrical location function ( , )

ii gq θ computed in the previous point ( ,

)
ii

q θ : 11 1 1 1 0 ( , ) ( , ) ( , ) ( , ) ( , ) ; ( , ) ( , )
q i i i i i i i i q i i i i i i T q i i i T i i i i                           Jqθ q J q θ θ t g q θ J q θ q J q θ θ Jq θ F 0 Jq θ F K θ K θ (8)
This gives the following iterative scheme

  1 1 1 10 ( , ) ( , ) , ( , ) ( , ) ( , ) 0 ( , )
i q i i i i i i q i i i i i i T i q i i T i i i                                          
F J q θ J q θ t g q θ J q θ q J q θ θ q J q θ θ J q θ K K θ [START_REF] Kang | Analysis and design of parallel mechanisms with flexure joints[END_REF] that may be reduced down to

1 1 1 1 1 1 0 1 ( , ) ( , ) ( , ) ; ( , ) ( , ) 0 0 T i i i i i q i i i T i i i i T i q i i                             F Jqθ K J q θ J q θ ε θ K J q θ F θ q
Jq θ [START_REF] Yong | Kinetostatic modeling of 3-RRR compliant micro-motion stages with flexure hinges[END_REF] where ( , ) ( , ) ( , )

i i i q i i i i i i        ε
t g q θ J q θ q J q θ θ . The latter is more convenient computationally, since it requires inversion of a lower dimension matrix    

66 nn    instead of     66 n m n m     
, where n , m are the sizes of the vectors q and θ respectively. For instance, for the kinematic chains of the Orthoglide manipulator (see example in section 6), expression (9) requires inversion of 34 34  matrix, while iterative scheme [START_REF] Yong | Kinetostatic modeling of 3-RRR compliant micro-motion stages with flexure hinges[END_REF]

K

is computed only ones, outside of the iterative loop. Similar to other iterative schemes, convergence of this algorithm highly depends on the starting point. However, due to physical nature of the considered problem, it is possible to start iterations from the non-loaded configuration ( 00 , q θ ). Besides, it is useful to modify the target point for each iteration in accordant with expression

  0 1 i i i       t t
t using scalar variable i  that is monotonically increasing from 0 up to 1. Another approach can be used for computing the forcedeflection curve. Here, the starting point can be taken from previously computed loaded configuration corresponding to another value of deflection that is very close to the target one. As it follows from computational experiments, for typical values of deformations the proposed iterative algorithm possesses rather good convergence (3-5 iterations are usually enough if a configuration is stable). In our simulation studies, the convergence was evaluated as the weighted sum of residual norms corresponding to equations ( 5) and the iterations were terminated when this criterion achieved prescribe value.

However, some computational difficulties may arise in the case of buckling or in the area of multiple equilibriums, where convergence problem becomes rather critical and highly depends on the initial guess. Here, number of iterations increase significantly and computational time becomes non-negligible. To overcome these difficulties, it is proposed to modify the developed iterative scheme and repeat computations several times, with slightly modified initial points that are obtained by adding small random noise to 00 , q θ . Another option is adding small disturbances to  

, ii θq at each iteration.

These techniques were used in case studies presented section 6.

The proposed iterative scheme can be also slightly modified to solve the dual problem: computing an equilibrium configuration corresponding to given external loading F (instead of given t ). In this case, expressions ( 9), [START_REF] Yong | Kinetostatic modeling of 3-RRR compliant micro-motion stages with flexure hinges[END_REF] are used in internal loop, while the designed algorithm is supplemented by external loop, which provides iterative searching for t corresponding to given

F   1 1 i i i i       t t K F F (11)
where , ii tF and i K are respectively the location, loading and stiffness matrix at the i-th iteration. It worth mentioning that the dual problem is meaningful only if the stiffness matrix i K is non-singular. It is obvious that for a separate serial chain with passive joints the matrix i K is always singular, while for the parallel manipulator it is usually non-singular due to specific assembling of kinematic chains. On the other hand, the original problem considered in this subsection (i.e. computing F corresponding to t ) is always physically meaningful and has at least one solution. More details concerning the developed algorithms are contained in web-appendix to this paper, which includes pseudo-code versions of relevant procedures [START_REF]Web-appendix[END_REF].

Stability of the loaded configuration

To investigate possible non-linear effects that may be caused by external and internal loading, it is necessary to extend the notion of stability associated with the stiffness analysis. Traditionally, the stability of compliant mechanical systems (including manipulators) is defined as resistance of end-point location t with respect to the "disturbing" effects of external force F applied at this point. In such formulation, the stability is completely defined by the stiffness matrix F K that describes linear relations (4) between the force and deflection deviations , δ F δ t with respect to the values F , t . It is obvious that the matrix F K is positive definite for the stable location t . However, in compliant manipulators with passive joints, the configuration ( , ) q θ corresponding to the same location t can be not unique. Moreover, these configurations may be both "stable" and "unstable" and correspond to different values of potential energy stored in the springs. From this point of view, it is worth to distinguish stability of end-point location t and stability of the corresponding configuration ( , ) q θ . This issue becomes extremely important for the loaded mode, when due to kinematic redundancy caused by passive joints and excessive number of virtual springs, small disturbances in ( , ) q θ may provoke essential change of current equilibrium configuration that leads to reduction of the potential energy and transition to another equilibrium state, while keeping the same end-location. Hence, it is necessary to evaluate internal properties of the kinematic chain in the state of the loaded equilibrium that may correspond either to minimum or maximum of potential energy for fixed value of t . Section 6 contains a number of illustrative examples that demonstrate equilibrium bifurcations while incrementing the external loading, buckling phenomena and difference between two types of stability notions mentioned above. Also, additional animated illustrations describing this issue are available in webappendix [START_REF]Web-appendix[END_REF].

To evaluate stability of the static equilibrium configuration ( , ) q θ , let us assume that the manipulator end-effector is fixed at the point ( , ) T  tp φ corresponding to the external load F , but the joint coordinates are given small virtual displacements  q , θ satisfying the geometrical constraint (1), i.e.

( , );

( , )       t g q θ t g q q θ θ [START_REF] Timoshenko | Theory of elasticity[END_REF] For these assumptions, let us compute the total virtual work in the joints that must be positive for a stable equilibrium and negative for an unstable one. To achieve the virtual configuration ( , )     qq θθ and restore the equilibrium conditions, each of the joints must include virtual motors that generate the generalized forces/torques q τ ,  τ which satisfies the equations:

  00 ; ( ) ( ) 0; ( ) TT TT q q q q                            J F K θ θ J J F K θ θ θ τ J F J J F τ (13)
After relevant transformations, the virtual torques may be expressed as

( ) ; ( ) TT qq               τ J F K θ τ J F (14)
where (...)



denotes the differential with respect to  q , θ that may be expanded via Hessians of the scalar function

( , ) T    gq θF : ( ) ; ( ) T F F T F F q q qq q                     J F H q H θ J F H q H θ ( 15
)
provided that

2 2 2 2 2 / ; / ; / F F F F qq q q                  H q H θ H H q θ (16)
Further, taking into account that the virtual displacement from ( , )

q θ to ( , )    
qq θθ leads to a gradual change of the virtual torques from (0, 0) to ( , ) q   ττ , the virtual work may be computed as a half of the corresponding scalar products  

1 2 TT q W            τ θ τ q , ( 17
)
where the minus sign takes into account the adopted conventions for the positive directions of the forces and displacements. Hence, after appropriate substitutions and transformation to the matrix form, the desired stability condition may be written as

1 0 2 FF q TT FF q qq W                         H K H θ θq
HHq [START_REF] Ciblak | Asymmetric Cartesian stiffness for the modeling of compliant robotic systems[END_REF] where  q and θ must satisfy to the geometrical constraints [START_REF] Timoshenko | Theory of elasticity[END_REF].

In order to take into account the relation between  q and θ that is imposed by [START_REF] Timoshenko | Theory of elasticity[END_REF], let us apply the first-order expansion of the function ( , ) g θq that yields the following linear relation

q          θ J J 0 q . ( 19
)
Then, applying the SVD-factorization [START_REF] Strang | Introduction to Linear Algebra[END_REF] of the integrated Jacobian [START_REF] Pigoski | Stiffness mappings employing different frames of reference[END_REF] and extracting from  V , q V the sub-matrices o  V , o q V corresponding to the zero singular values, a relevant null-space of the system (19) may be presented as oo ;

T r qq T q                     V S J J U U V 0
q          θ V μ q V μ ( 21
)
where  μ is the arbitrary vector of the appropriate dimension (equal to the rank-deficiency of the integrated Jacobian).

Hence, the stability condition (18) may be rewritten as an inequality [START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF] that must be satisfied for all arbitrary non-zero  μ . In other words, the considered static equilibrium ( , ) q θ is stable if (and only if) the matrix

o o o o 1 0 2 T FF q T FF q q qq W                                      V H K H V μμ V H H V
o o o o 0 T FF q FF q q qq                               V H K H V V H H V (23)
is negative-definite. It is worth mentioning that the obtained result is in a good agreement with previous studies [START_REF] Sh | Geometrical Approach to The Conservative Congruence Transformation (CCT) for Robotic Stiffness Control[END_REF], where (for manipulators without passive joints) the stiffness properties were defined by the matrix

F    KH
that evidently must be positive-definite for stable configurations. In Section 6 these results are applied for detecting bifurcations and buckling phenomena in typical serial kinematic chains and parallel manipulators.

Stiffness matrix for the loaded mode

The previous section presents technique that generally allows obtaining an exact relation between the elastic deformations and corresponding external force/torque. It is based on sequential computations of loaded equilibriums (and relevant force/torque) for various displacements of the manipulator end-point with respect to its unloaded location. However, though this relation is highly non-linear, common engineering practice operates with the stiffness matrix derived via the linearization.

To compute the desired stiffness matrix, let us consider the neighborhood of the loaded configuration and assume that the external force and the end-effector location are incremented by some small values F , t . Besides, let us assume that a new configuration also satisfies the equilibrium conditions. Hence, it is necessary to consider simultaneously two equilibriums corresponding to the manipulator state variables ( , , , )

Fqθt and (

, , , )         F F q q θ θ t t .
Relevant equations of statics may be written as

  0 ;0 TT q        J F K θ θ J F (24)
and

          0 ; 0 T T                 θ θ θ qq J J F F K θ θ θ J J F F (25)
where ( , ) q Jq θ and ( , )

 Jq θ are the differentials of the Jacobians due to changes in ( , ) q θ . Besides, in the neighborhood of ( , ) q θ , the kinematic equation may be also presented in the linearized form:

( , ) ( , )

q      δ t J q θ δ θ J q θ δ q , (26)
Hence, after neglecting the high-order small terms and expanding the differentials via the Hessians of the function ( , ) T    gq θF (similar to sub-section 4.2), equations ( 24), ( 25) may be rewritten as

( ) ( ) ( ) ( ) ( ) ( ) T F F q T F F q qq q                         
J q , θ F H q , θ q H q , θ θ K θ J q , θ F H q , θ q H q , θ θ 0 [START_REF] Hu | Finite Element Analysis of a Six-Component Force Sensor for the Trans-Femoral Prosthesis[END_REF] and the general relation between the increments of the state variables can be presented as

q T F F q qq q T F F q                                     0 J J F t J H H q 0 J H H K θ0 (28)
The latter gives a straightforward numerical technique for computing the desired stiffness matrix: direct inversion of the matrix in the left-hand side of ( 28) and extracting from it the upper-left sub-matrix of size 66. Similarly, the matrices defining linear relations between the end-effector increment t and the increments of the joint variables θ ,  q can be computed, i.e.:

;;

Cq              F K t θ K t q K t (29)
where

1 ** Cq T F F q q qq q T F F q                         K 0 J J K J H H K J H H K (30)
It worth mentioning that the internal preloading (expressed by the variable 0 θ ) is not included in the latter expression in explicit way, but it directly influences on the variables ( , ) q θ describing the equilibrium configuration and corresponding Jacobians and Hessians, which are elements of [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF]. Besides, in contrast to previous works, here it is possible to obtain supplementary matrices , q  KK that give additional measures of the manipulator stiffness which evaluate sensitivity of the joint coordinates ( , ) q θ with respect to the external loading. In the case when the matrix inverse (30) is computationally hard, the variable θ can be eliminated analytically, using corresponding static equation:

F T F F q              θ k J F k H q , where   1 F      F k K H
. This leads to a reduced system of matrix equations with unknowns F and  q

F T F F q T F F T F F F F q q qq q q                                        q θ J k J J J k H δF δt J H k J H H k H δq 0 . ( 31
)
that may be treated in the similar way, i.e. the desired stiffness matrix is also obtained by direct inversion of the matrix in the left-hand side of ( 31) and extracting from it the upper-left sub-matrix of size 66:

1 F T F F C qq T F F T F F F F q q q qq q q                                    K J k J J J k H K J H k J H H k H (32)
Similar to subsection 4.1, this approach allows reducing the dimension of the inverted matrix from

    66 n m n m      to     66 nn   
, that in the case of Orthoglide corresponds to [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF] It is worth mentioning that the structure of the latter matrix is similar to one obtained for the unloaded manipulator in [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] and differs only by Hessians that take into account the influence of the external load. It should also be noted that, because of the presence of the passive joints, the stiffness matrix of a separate serial kinematic chain is always singular, but aggregation of all the chains for a parallel manipulator produces a non-singular stiffness matrix. Some examples of such aggregation are presented in our previous papers [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF], where the stiffness matrices of separate kinematic chains are expressed with respect to the same end-point and the resulting stiffness matrix for the parallel manipulator is computed as a straightforward sum of these components, i.e.

() i C C i   K K
where the index "(i)" denotes the number of kinematic chain.

The later follows from the superposition principle, since the total external force corresponding to the same end-effector displacement is expressed as the sum of partial forces of the separate chains. Further simplification can be obtained by applying the block matrix factorization technique (Frobenius equation) [START_REF] Gantmacher | Theory of matrices[END_REF] that yields the following expression

                  11 1 1 1 F T F T F F C q q F F F F T F F T F T F F qq q q q q q q T F F T F T qq                                                                  K J k J J k J J J k H H H k H J H k J J k J J J k H J H k J J k J ( 33
)
where the first term exactly corresponds to the classical formula defining stiffness of the kinematic chain without passive joints.

Hence, the presented technique allows computing the Cartesian stiffness matrix in the presence of the external and internal loading. It generalizes our previous results both for serial kinematic chains and for parallel manipulators [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF]. In the following Section, it will be applied to several case-studies that deal with kinematic chains employed in typical parallel manipulators and demonstrate particularities of stiffness analysis of loaded manipulator with passive joints.

Illustrative examples

To demonstrate the efficiency of the proposed technique, let us present two examples that focus the loaded working mode and deal with (i) the stiffness analysis of a serial kinematic chain with passive joints and (ii) the stiffness analysis of a translational parallel manipulator. Within this study, it is assumed that the external wrench is acting upon the manipulator end point and the influence of the gravity is negligible. In previous works [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF], this problem was studied for the unloaded working mode and/or for small deflections in elastic elements, so any non-linear abnormality in the manipulator behavior (buckling, etc.) was non detected. This Section includes only summary of this computational study, while more details are presented in web-appendix [START_REF]Web-appendix[END_REF]

Stiffness analysis of a serial kinematic chain

First, let us consider a serial kinematic chain consisting of three similar links separated by two similar rotating actuated joints. It is assumed that the chain is a part of a parallel manipulator and it is connected to the robot base via a universal passive joint, and the end-platform connection is achieved via a spherical passive joint. For each of these configurations, let us investigate three types of the virtual springs corresponding to different physical assumptions concerning the stiffness properties of the actuators/links. They cover the cases, in which the main flexibility is caused by the torsion in the actuators, by the link bending, and by the combination of elementary deformations of the links.

6.1.1

Geometric model The geometry of the examined kinematic chain (Fig. 3) can be defined as U p R a R a S p where R, U and S denote respectively the rotational, universal and spherical joints, and the subscripts "p" and "a' refer to the passive and active joints respectively. Using the homogeneous matrix transformations, it can be described by the equation L is the length of the links, 0 q is the coordinate vector of the universal passive joint located at the robot base, t q is the coordinate vector corresponding to the passive spherical joint at the end-platform, For presentation convenience, let us also assume that the coordinates 0 q of the universal passive joint are computed to ensure location of the end-effector on the Cartesian axis x. Besides, it is assumed that the external loading is presented by a compressing force applied at the chain end-point and it is directed along x-axis (see Fig. 3).

0 1 1 2 2 3 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) u x s z a x s z a x s s t L q L q L           T R q T T θ R T T θ R T T θ R q (34

6.1.2

Stiffness models In order to investigate possible non-linear effects in the stiffness behavior of this chain, let us consider several cases that differ in stiffness models of the links and actuated joints:

Case of 1D-springs (Model A): the flexible elements are localized in the actuating drives while the links are considered as strictly rigid. It allows, without loss of generality, to reduce the original U p R a R a S p model down to R p R a R a R p and define a single stiffness parameter K  (similar for both actuators) that will be used as a reference value for the further analysis. Besides, it is possible to ignore the end-effector orientation and consider a single passive joint coordinate q (at the base) and two virtual joint coordinates 1  , 2  (at actuators). Case of 2D springs (Model B): the actuators do not include flexible components but the manipulator links are subject to non-negligible elastic deformations in Cartesian xy-plane (bending and compression). Correspondingly, the link flexibility is defined by a 33 matrix that includes elements describing deformation in x-and y-directions and rotational deformation with respect to z-axis.

Case of 3D springs (Model C): the actuators are strictly rigid but the link flexibility is described by a full-scale 3D model that incorporates all deflections along and around x-,y-,z-axes of the three-dimensional Cartesian space. Relevant stiffness matrix of the links has the dimension 66. The kinematics of this model corresponds to the general expression U p R a R a S p , it includes two passive joints  

, s qq incorporating in total five passive coordinates and three virtual-springs with 18 virtual coordinates totally (six for each link).

Stiffness analysis for model A

In this case, the model includes minimum number of flexible elements (two 1D virtual springs in the actuated joints) and may be tackled analytically. However, in spite of its simplicity, it is potentially capable to detect the buckling phenomena at least for S-configuration, because of evident mechanical analogy to straight columns behavior under axial compression. It is also useful to evaluate other initial configurations (-and Z-types), their stability and to compare analytical solutions with numerical results provided by the developed technique.

For this model, consistent solution of geometrical and static equations (5) yields two expressions for "force-deflection" relations that correspond to stable and unstable equilibriums

  sin s K F L       ;   cos( ) 2 cos sin n K qq F L            (35)
where F are the force acting along the x-axis,  is corresponding linear displacement, the subscripts "s" and "n" denote manipulator configuration stability (stable or unstable), K  is actuator stiffness, L is the length of the links, and auxiliary variables , , q

  are expressed via normalized displacement / L   as arccos(1 / 2)      ,       2 arccos 12 6 / 12 4 q         ,   2 arccos 1 3 / 2 / 4       
. It worth mentioning that other possible deformations (along y-axis) are not considered here because, as follows from separate study, in this direction the kinematic chain is kinematically singular due to passive joints and, obviously, the chain does not demonstrate any non-linear phenomena.

Corresponding force-deflections caves are presented in Fig. 4 where there are also shown the bifurcation and analytical expressions for F corresponding to small values of  that are derived from ( 35) by parametric differentiation. The interpretation of this plot is similar to the axial compression of a straight column, which is a classical example in the strength of materials. It should be noted, that the developed numerical algorithm exactly produces the curve corresponding to the stable equilibrium. In this case, the external force

/ F K L  
can not change the manipulator shape, similar to small compressing of straight columns that can not cause lateral deflections. Consequently, for such loading the straight configuration is stable. Further, for

/ 3 / K L F K L  
, the straight configuration may be hypothetically restored but becomes unstable, so any small disturbance will cause sudden reshaping in the direction of the stable trapezoid-type posture. And finally, for

3/ F K L  
, there may exist two types of unstable equilibriums: the trivial straight-type and a more complicated zig-zag one presented in Fig. 4.

Fig. 4. Model A: Force-deflection relations and bifurcations for the initial S-configuration

For other initial shapes (-type and Z-type), the results essentially differ from the above ones (see Fig. 5). In particular, for small deflections, both -and Z-configuration demonstrate rather linear behavior. Moreover, in most of the cases there exist a single stable and a single unstable equilibrium, so the kinematic chain can not suddenly change its shape due to external loading. The only exception is the case of the initial -configuration where there are two stable and two unstable equilibriums, and here there exists a bifurcation of the stable equilibriums corresponding to the cuspidal point of the function () F  where the stiffness reduces sharply. Another conclusion concerns the profile of the force-deflection plots that are highly nonlinear in all cases.

More detailed analysis shows that -configuration demonstrates good analogy with axially compressed imperfect column where the deflection starts from the beginning of the loading and there is no sudden buckling, but the stiffness essentially reduces while the loading increases. However, for Z-configuration that corresponds to the unloaded zig-zag shape, the stiffness behavior demonstrates the buckling that leads to sudden transformation from a symmetrical to a nonsymmetrical posture as shown in Fig. 6. Here, there exist four stable equilibriums (and two unstable ones) that differ in the values of the potential energy.

Hence, in the case of model A, the developed numerical technique provides solutions that correspond to the stable loaded equilibriums and coincide with relevant analytical results. In spite of simplicity of this case study, it allows developing intuition for potentially dangerous manipulator postures with respect to buckling.

6.1.4

Stiffness analysis for Model B In this case, the manipulator stiffness is caused by elasticity in the links while the actuating joints are assumed to be rigid. The elastic deflections (bending and compression) are still restricted by the Cartesian xy-plane but each link includes three virtual joints. Totally, the stiffness model has 11 variables, so it was studied numerically (using the developed technique). The stiffness parameters of the elements were evaluated assuming that the links are the rectangular beams of the length L and the cross-section ab, where 0.02 aL  and 0.05 bL  . For comparison purposes, corresponding stiffness matrices were scaled with respect to the bending coefficient 1/ K  , to keep similarity with model A. The stiffness analysis was performed for above defined S-, and Z-configurations, assuming that the external force is directed along the x-axis. Summary of the modeling results are presented in Fig. 5 and are briefly described below.

For S-configuration, here there is still very strong analogy with the compression of the straight column. In particular, at the beginning of the loading, the links are subject the axial compression and the stiffness is very high. Then, after the buckling, the manipulator changes its shape to become non-symmetrical and the stiffness falls down. The critical force may be also computed using the previous results, as 0

/ F K L  
. For -configuration, the stiffness properties are also qualitatively equivalent to the case of Model A but the stiffness coefficient is slightly lower. For Z-configuration, it was also detected the buckling that occurs if the loading approaches to the critical value 0 1.07 Hence, stiffness analysis of Model B demonstrates qualitative similarity but some quantitative difference compared to Model A. The latter is caused by different arrangement of the elastic elements in the virtual joints that corresponds to other physical assumptions.

/ F K L   .

6.1.5

Stiffness analyses for the model C Finally, for model C, the link elasticity is described in 3D space and corresponding stiffness matrices have dimension 66 (the actuating joints are assumed perfect and rigid, similar to model B). It is also assumed that the links are rectangular beams of the length L with the cross-section ab, where 0.02 aL  , 0.05 bL  and the smaller value a corresponds to zdirection, which was not studied above (in models A, B). To ensure comparability of all examined cases, the link stiffness matrices were parameterized with respect to the bending coefficient of the z-axis K  . In total, the stiffness model includes 23 variables (five for passive joints and 18 for the virtual springs of three links) and it was studied numerically. The stiffness analysis was performed for the same initial configurations (S-, and Z-type) and for the same direction of the external force as for the models A and B.

For S-configuration, the results are qualitatively similar to ones obtained for model B, but here the buckling occurs for essentially lower critical force, 0.16 / KL  , that corresponds to sudden lateral deflection in z-direction 4 . It is worth mentioning that the axial deflection corresponding to the critical force is very low (see [START_REF]Web-appendix[END_REF] for details). In contrast, for configuration, it was detected buckling which does not exist in models A and B. In particular, if the external force exceeds the critical value 0.20 / KL  the stiffness suddenly reduces (it becomes 25 times lower). It corresponds to sudden deflection in z-direction that it was beyond capabilities of previous models. Another interpretation of this buckling phenomenon may be presented as sudden loss of symmetry with respect to xy-plane. For Z-configuration, the results remain qualitatively the same as above, but corresponding critical force is 0.17 / KL  and it is also defined by buckling in z-direction. Therefore, model C yields essentially lower values of critical force compared to models A and B. Besides, here all examined postures demonstrated buckling in z-direction. This presents another source of potential structural instability of kinematic chains that posses the symmetry with respect to a plane.

Hence, the presented results concerning stiffness analysis of serial chains with different types of compliant elements confirm capability of the developed technique to detect (and to evaluate numerically) non-linear effects that may arise due to external loading. Besides, these examples demonstrate meaningful analogy with stiffness behavior of perfect and nonperfect columns under axial compression. It is obvious that, due to potential danger of buckling, classical procedures of manipulator design (including optimization of link cross-sections) must take into account this phenomena and the developed technique contributes to this issue.

Stiffness analysis of translational parallel manipulator

Let us consider now a more sophisticated example that deal with the stiffness analysis of 3-d.o.f. translational manipulator Orthoglide. This manipulator consists of three identical kinematic chains actuated by mutually orthogonal linear drives (see Fig. 7), which are arranged in a special way to ensure almost isotropic workspace kinematic properties and to restrict the end-effector motions in translation only. Such architecture ensures the velocity transmission factors close to 1.0 (similar to the conventional Cartesian-type machines) and provides high stiffness in the non-loaded mode [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF]. However, nonlinear stiffness properties of such manipulator and potential buckling phenomena have not been studied yet.

Each kinematic chain of Orthoglide manipulator includes a foot, a kinematic parallelogram with two bars and two coupling elements, and an end-effector. The elasticity of the manipulator links was evaluated using the FEA-based methodology and special accuracy improvement tools proposed by the authors [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]. As it follows from our study, the rigidity of the parallelogram axel is high compared to the bar and the foot, thus the limb stiffness is described by equivalent 5-d.o.f. virtual spring (see [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]). For the remaining components, the compliance matrices are was taken from [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF].

For the studied manipulator, the stiffness analysis was performed for vertices Q 1 …Q 4 of the dexterous cubical workspace 200200200 mm 3 corresponding to the Orthoglide design specification with the range of the velocity transmission factors [0.5, 2.0]. Besides, similar study was done for the isotropic point Q 0 where the transmission factors are equal to 1.0, similar to conventional Cartesian machine [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF][START_REF] Pashkevich | Design Strategies for the Geometric Synthesis of Orthoglide-type Mechanisms[END_REF]. For all of these points, there were computed non-linear curves of the force-deflection relations assuming that the displacement is directed along the bisecting line of the coordinate system (this force direction is chosen in accordance to the workpiece orientation in milling application which is currently developing for Orthoglide).

As it follows from the obtained results (Fig. 7, Table 2), for all considered points Q 0 …Q 4 , the force-deflection relations are essentially non-linear and include the buckling points where the stiffness reduces radically. Corresponding critical force varies from 2.2 kN for the point Q 2 to 7.8 kN for the point Q 1 , while the critical deflections are in the range 0.5…4.9 mm. The isotropic point Q 0 that is usually used for the benchmarks of such manipulators possesses intermediate values of these indices (critical force 4.6 kN and critical deflection 1.4 mm).

For all cases, the translational stiffness coefficient in the pre-buckling mode is almost constant, but it essentially varies throughout the workspace. In particular, the highest value 16.6 kN/mm is achieved at the Point Q 1 and the lowest one 0.4 kN /mm corresponds to the point Q 2 . However, after the buckling, the stiffness abruptly reduces down to 0.03 0.84 kN/mm  and further, for large deformations, it is about 0.03 0.31 kN /m m  . Hence, the traditional linear stiffness analysis (which ignores the influence of external loading) provides rather reasonable indicators for the prebuckling mode but it does not allow evaluate the range of the loading for which the deflections may be treated as "small" ones and the buckling phenomena do not appear yet. It is worth also mentioning that for the points Q 0 …Q 4 , the physical nature of the buckling is perfectly explained by results from the previous subsection. In these cases, the buckling is caused by sudden change of the configuration of one of the serial chains (from S-or Z-to -configuration, for instance). It worth mentions, that parallelograms introduce some particularities that lead to additional bifurcations and salient points on the force-deflection curves. Besides, in the area of flat singularity (point Q 2 ) it is required some complimentary analysis (for more details see [START_REF]Web-appendix[END_REF]).

Therefore, for this more sophisticated case, the developed technique also allows detecting non-linear phenomena in manipulator stiffness behavior. Moreover, most of these phenomena may be explained and predicted using presented above results for serial chains (potential dangerous configurations, corresponding critical forces and stiffness, etc.). But due to particularity of parallel manipulator geometry, it was detected another type of buckling that is caused by specific spatial location of the kinematic chains with respect to each other.

Discussion

The presented results generalize previous research on the stiffness analysis of the manipulators that usually ignore the effect of the external and internal loading. Moreover, the derived expressions [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF], [START_REF] Long | Optimal structural design of a planar parallel platform for machining[END_REF] for the stiffness matrix are straightforwardly reduced to the known ones corresponding to more simple cases: (i) loaded manipulator without passive joints [START_REF] Chen | Conservative Congruence Transformation for Joint and Cartesian Stiffness Matrices of Robotic Hands and Fingers[END_REF][START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF], and (ii) unloaded manipulator with passive joints [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. Another advantage is that here, in contrast to previous works, the Jacobians and Hessians are computed for the joint variables corresponded to the exact loaded configuration while common approach admits that loading non-essentially influence on the Jacobian and Hessian matrices. Obviously, the proposed approach essentially improves accuracy of the stiffness model and provides it with ability of detecting various buckling phenomena.

However, there is a principal limitation of this method that perfectly takes into account non-linearity of the manipulator geometry but nevertheless uses linear expressions to describe the elasticity of each links (it is a fundamental assumption of all VJM-based models). As follows from our additional study, this limitation is not crucial for robotic manipulators that are design to operate in pre-buckling regime of the structural element where the link deflections are still rather small. Besides, it this research, it is assumed that the external loading is concentrated at the centre point of the end-effector while in general case the loading may be applied to the manipulator joints, or it may be distributed rather than concentrated. Another restriction is related to the manipulator architecture that is assumed to be assembled from strictly parallel serial chains (without additional geometric constraints or cross-linking between the chains). These issues are evidently left for further work.

It is also worth comparing the developed approach with other techniques that are used in neighbouring area, the structural mechanics. At present, there are a number of commercial finite-element software systems that are able to take into account the effects related to the geometric nonlinearity and to evaluate the critical load level and the loaded structure stability through the eigenvalue analysis of the global stiffness matrix. Usually, they implement two types of the numerical procedures: (i) computationally inexpensive "linearized buckling analysis" that assumes that the equilibrium point is very close to the initial configuration (sometimes it includes several iterations in the displacement space); and (ii) numerically cumbersome and straightforward "non-linear buckling analysis" that involves complete incremental tracing of the loaddeflection path. Within both of these techniques, there are two basic versions that assume respectively that the search/tracing is performed (a) in the displacement space for a given fixed load vector, and (b) in the load space for a given fixed displacement. In commercial systems, advantages are usually given to the case (i,b) that provides good numerical stability and efficiency but affords rather moderate accuracy for the critical points.

In the frame of the above classification, the proposed technique corresponds to the case (ii,b). It allows efficiently tracing the complete load-deflection path without visible computational efforts (which is the main difficulty of the FEAbased techniques) and also admits easy incorporation of relevant iterative procedures for assessing the critical force/deflection. This numerical efficiency originates from the VJM-presentation that essentially reduces the matrix size within computations. For comparison, for the prototype of the Orthoglide manipulator created in IRCCyN (Institut de Recherche en Communications et Cybernetique de Nantes), the straightforward finite element analysis leads to the matrix inversion of size about 10 4 10 4 while the developed technique requires several inversions for matrices of size 1010 only.

Hence, the developed technique is perfectly adapted to the stiffness analysis of the robotic manipulators. It is able to take into account influence of the external loading and allowed detecting some uncommon behavior of the robotic manipulators which was not previously reported in robotic literature. These phenomena are directly related to the robot accuracy and must be obviously taken into account during design and analysis.

Conclusions

The stiffness analysis becomes a critical issue in design optimization of robotic manipulators that currently is targeted at achieving high dynamic performances with relatively small link masses and low energy consumption in actuators. This tendency motivates a revision of existing stiffness modeling techniques that must take into account the external and internal loading. To meet this demand, the paper proposes a new systematic method for stiffness modeling of robotic manipulators with passive joints in loaded mode. It is based on a multidimensional lumped-parameter model, which presents the links as pseudo-rigid bodies with 6-d.o.f. virtual springs whose parameters are evaluated via the FEA-modeling and describe both the translational/rotational compliances and the coupling between them. The developed technique allows computing the loaded equilibrium configurations, evaluating their stability, and finding the full-scale "load-deflection" path for any given displacement of the end-effector. It is also proposed the linearization procedure for computing the Cartesian stiffness matrix, which is based on the inversion of the dedicated matrix of larger dimension, composed of the stiffness parameters of the virtual springs and the Jacobians/Hessians of the active and passive joints. These results enable designer to evaluate critical forces that may provoke non-linear behavior of the manipulators, such as sudden failure due to elastic instability (geometrical buckling) which has not been previously studied in robotic literature.

The advantages of the developed technique are illustrated by several examples that deal with parallel manipulator of the Ortholide family and serial kinematic chains employed in such robots. They demonstrate possible non-linear effects that may arise in loaded mode, including essential dependence of the stiffness on the applied force/torque and sudden change of the stiffness if the external wrench exceeds a critical value. Besides, several typical configurations of serial kinematic chains that are potentially dangerous with respect to buckling were detected.

While applied to pure parallel mechanisms with similar kinematic chains and actuators located between the base and foot, the method can be extended to other parallel manipulators to cover different actuator locations and dissimilar chain geometry. So, future work will focus on the stiffness modeling of more sophisticated architectures that include parallel manipulators with the cross-linking between the main kinematic chains and other structures that improve rigidity of the manipulating system. Another prospective research direction is the stiffness analysis of heavy manipulators, for which external loading is caused by gravity and it is distributed within the links. [START_REF] Ciblak | Synthesis of stiffness by springs[END_REF]

FK

  solution of non-linear matrix equation that includes joint stiffnesses and external/internal loadings as parameters c K -Cartesian stiffness at the end-effector ( 66  )  K -aggregated stiffness of the virtual springs ( nn   )  J -Jacobian of the virtual springs ( -stiffness matrix induced by external loading ( nn   ) I K -stiffness matrix induced by constraints (cross-linking) ( nn   )

Fig. 1 .

 1 Fig. 1. Serial kinematic chain with passive joints and its VJM model (Acactuated joint, Pspassive joint, -virtual spring).

Fig. 2 .

 2 Fig. 2. Examples of auxiliary springs in preloaded passive joints.

   φ of the end-platform (or end-effector) in the Cartesian space, the vector 12 ( , , ..

34 

 34

  rotation/translation matrices around/along the z-and x-axes, the spherical joint (incorporating three elementary rotations), 12 , aa qq are the coordinates of the actuated joints,

T

  is the homogenous matrixfunction describing elastic deformations in the links and actuators (they are represented by the virtual coordinates incorporated in the vectors 1 2 3 ,, θ θ θ ). It is obvious that this model can be easily transformed into the form ( , )  t g q θ used in the frame of the developed technique.

Fig. 3 .

 3 Fig. 3. Examined kinematical chain and its typical configurations ( U ppassive universal joint, R a1 , R a2actuated rotating joints, S ppassive spherical joint) To investigate particularities of this kinematic chain with respect to the external loading, let us also consider three typical postures that differ in values of the actuated coordinates:  S-configuration: the links are located along the straight line (Fig. 3a), the actuated coordinates are 12 0 aq qq 

Fig. 5 .Fig. 6 .

 56 Fig. 5. Force-deflection relations for different elasticity models of serial chain with passive joints

Fig. 7 .

 7 Fig. 7. Stiffness of the Orthoglide manipulator in critical points of the cubic workspace

  needs inversion of10 10

				
	matrix only. It should be mentioned that			1

Table 2

 2 Summary of stiffness analysis for Orthoglide manipulator in different workpoints

	Configuration	Stiffness for unloaded mode, [kN/mm]	Critical force, [kN]	Critical deflection , [mm]	Stiffness around the buckling. [kN/mm] cr cr FF  FF 	Stiffness for large deformations, [kN/mm]
	Point Q0	3.23	4.65	1.41	3.36	0.03	0.03
	Point Q1	16.60	7.85	0.48	16.40	0.42	0.31
	Point Q2	0.42	2.24	4.85	0.48	0.06	0.03
	Point Q3	9.95	4.75	0.48	9.60	0.84	0.11
	Point Q4	5.41	2.78	0.52	5.09	0.65	0.06

Here, we distinguish the mechanical stiffness of the manipulator (passive compliance) and the stiffness of the robotic system, which may be adjusted by regulating the virtual elasticity of actuators due to dedicated control algorithm (active compliance).

In general, for non-conservative systems and/or some special parameterisations of end-effector location, in the loaded mode the stiffness matrix may be asymmetrical This issue is discussed in details in[START_REF] Yi | Synthesis of actively adjustable springs by antagonistic redundant actuation[END_REF][START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF][START_REF] Griffis | Global stiffness modeling of a class of simple compliant couplings[END_REF][START_REF] Ciblak | Asymmetric Cartesian stiffness for the modeling of compliant robotic systems[END_REF][START_REF] Ciblak | Synthesis of stiffness by springs[END_REF][START_REF] Pigoski | Stiffness mappings employing different frames of reference[END_REF][START_REF] Howard | On the 6 x 6 Cartesian stiffness matrix for three-dimensional motions[END_REF] and others.

It should be noted that, in the frame of VJM approach, the Cartesian stiffness matrix may be computed for all manipulator postures, including singular ones (singular in usual sense, with respect to actuated joints). This is because the virtual joints produce the Jacobian which is always non-singular, independently of spatial location of the links. Some examples of the stiffness matrices for singular configurations can be found in[START_REF] Pashkevich | Stiffness analysis of 3-d.o.f. overconstrained translational parallel manipulators[END_REF].

For comparison, according to the Euler formula, the local buckling of the links occurs for compressing force of 0.40 K/L, which is 2.5 times higher than the loading that provokes the geometric buckling of the examined mechanism.

Acknowledgements

The work presented in this paper was partially funded by the Region "Pays de la Loire", France and by the EU commission (project NEXT).