A Configurable Architecture for Two Degree-of-Freedom Variable Stiffness Actuators to Match the Compliant Behavior of Human Joints

Living beings modulate the impedance of their joints to interact proficiently, robustly, and safely with the environment. These observations inspired the design of soft articulated robots with the development of Variable Impedance and Variable Stiffness Actuators. However, designing them remains a challenging task due to their mechanical complexity, encumbrance, and weight, but also due to the different specifications that the wide range of applications requires. For instance, as prostheses or parts of humanoid systems, there is currently a need for multi-degree-of-freedom joints that have abilities similar to those of human articulations. Toward this goal, we propose a new compact and configurable design for a two-degree-of-freedom variable stiffness joint that can match the passive behavior of a human wrist and ankle. Using only three motors, this joint can control its equilibrium orientation around two perpendicular axes and its overall stiffness as a one-dimensional parameter, like the co-contraction of human muscles. The kinematic architecture builds upon a state-of-the-art rigid parallel mechanism with the addition of nonlinear elastic elements to allow the control of the stiffness. The mechanical parameters of the proposed system can be optimized to match desired passive compliant behaviors and to fit various applications (e.g., prosthetic wrists or ankles, artificial wrists, etc.). After describing the joint structure, we detail the kinetostatic analysis to derive the compliant behavior as a function of the design parameters and to prove the variable stiffness ability of the system. Besides, we provide sets of design parameters to match the passive compliance of either a human wrist or ankle. Moreover, to show the versatility of the proposed joint architecture and as guidelines for the future designer, we describe the influence of the main design parameters on the system stiffness characteristic and show the potential of the design for more complex applications.


INTRODUCTION
In the last years, articulated soft robots inspired by the musculoskeletal system of vertebrate animals received increased attention from researchers, since they represent promising solutions to enhance the interactions of robots with unknown and dynamic environments, i.e., the real world (Albu-Schaeffer et al., 2008). In line with this, Variable Impedance Actuators and the subgroup of Variable Stiffness Actuators (VSAs) were widely investigated recently (Vanderborght et al., 2013).
Variable Stiffness Actuators usually rely on two main approaches to modulate the stiffness. In the first one, the stiffness is controlled using software together with a mechanism with fixed impedance properties. In the other one, the stiffness is modulated through a mechanical reconfiguration of the system (Vanderborght et al., 2013). Some literature refers to the first approach as active VSA and the second one as passive VSA (Vanderborght et al., 2013). However, to avoid confusion possibly induced by the notion of passivity, we will refer to the first approach as software-controlled VSA and to the second one as physically compliant VSA. Software-controlled VSAs allow the design of lightweight devices that can in theory simulate any desired stiffness behavior. However, this apparent stiffness emulated by the system relies on an accurate sensing strategy and control computations, and it has been shown that even if the impacts are detected timely, the motors could not be able to react fast enough by solely an impedance control and that the system should therefore be considered as stiff during the impacts (Haddadin et al., 2007). To address these limitations, physically compliant VSAs have been developed with an inherent compliance. They can present several advantages such as shock absorptions, better performances in cyclic or explosive tasks (Albu-Schaeffer et al., 2008;Wolf et al., 2016), and a possible embodiment of specific behavior to improve the control strategies (Visser et al., 2011).
As VSAs are inspired by musculoskeletal systems, their applications as prostheses or part of humanoid devices would seem straightforward. Even though there are some advantages in using VSAs compared to stiff actuators, finding concrete usecases of VSAs is an on-going research topic and requires attention (Sensinger and ff. Weir, 2008;Blank et al., 2014;Stillfried et al., 2018).
To test hypotheses in this objective, we need multi-degree-offreedom (DoF) joints that show abilities similar to those of a human joint, such as the same functional range of motion, a variable stiffness mechanism, the overall shape and mass, and so on. Integrating an inherent compliance inside these systems enables us to test both the variable stiffness ability but also the capabilities due to the compliance, such as the possibility of exploiting the natural dynamic of the system (e.g., the resonance frequency) (Stillfried et al., 2018). Besides, artificial wrists can benefit from a compliant architecture to enhance their manipulation abilities in tight or cluttered spaces (Negrello et al., 2019). Moreover, the availability of various levels of stiffness can enhance the performances in activities of daily living in prosthetic applications (Kanitz et al., 2018). However, the mechanical complexity of such mechanisms increases and constitutes a challenge.
Up to now, most of the proposed physically compliant VSAs have only 1 DoF in position (Vanderborght et al., 2009;Vanderborght et al., 2013;Petit et al., 2015;Lemerle et al., 2019). The classical approach to get multi DoF physically compliant VSAs is to put in series several 1-DoF VSAs, such as the solutions proposed in the DLR Hand Arm System (Grebenstein et al., 2011;Romano et al., 2014;Savin et al., 2019). This approach uses the design of serial manipulators (SMs) that are easier to model and control compared to the parallel manipulators (PMs). The latter, on the other hand, allow more compact designs and have better performances in terms of output torques than SMs for the same motor sizes. Nevertheless, the performances of PMs are very sensitive to their geometric parameters (Siciliano and Khatib, 2007). Alternative works explored the independent setup approach, consisting mainly of a serial arrangement of position motor units with several DoF in position (using either SMs or PMs) and an additional motor unit dedicated solely to the control of the stiffness (Weckx et al., 2014;Barrett et al., 2017;Tan et al., 2017;Rodriguez-Cianca et al., 2019). Finally, other works draw inspiration from the antagonistic arrangements of the musculoskeletal system (Koganezawa, 2017;Stoeffler et al., 2018;Malosio et al., 2019). In the last category, several antagonistic motors are used jointly to control the stiffness and the position of the joint. Theses architectures are solely based on PMs. For more details on the classification and principle of work of VSAs, the reader is invited to refer to (Vanderborght et al., 2013;Wolf et al., 2016).
One limitation of physically compliant VSAs is that one mechanical implementation corresponds to one embedded passive behavior. Although there are devices that are humanlike in terms of shapes and dimensions, and are able to perform activities of daily living (ADL) (Grebenstein et al., 2011;Weckx et al., 2014;Koganezawa, 2017;Stoeffler et al., 2018;Rodriguez-Cianca et al., 2019), there is, for now, no multi-DoF VSA embedding a compliant behavior that matches human joints. It is worth mentioning that there are alternatives to design compact 2 DoF joints with inherent compliance (not necessarily with variable stiffness in the following examples). A first concept is to design an underactuated mechanism (Casini et al., 2017). Its main drawback could be its limited dexterity to do all the desired tasks. Another approach is to have a switchable stiffness to get lightweight and compact devices (Montagnani et al., 2013). This type of solution could be suitable when lightweight systems are needed. However, they do not allow for a continuous and smooth control of the stiffness. Therefore, they lose some advantages of such abilities, e.g., during explosive or cyclic tasks, when a VSA performs better than a fixed compliant mechanism . An additional idea is to use tendon driven mechanisms with compliant pneumatic actuation (Toedtheide and Haddadin, 2018). Then, this system could be lightweight and somehow compact if the remote actuation units are excluded. Yet, this is not fully satisfying in terms of compactness if we consider them, especially when they are pneumatic actuation units. Therefore, all of these approaches do not fully satisfy the criteria of compactness and variable stiffness ability as well as doing various types of tasks. We address this limitation in the following work.
In general, the qualitative representation of the mechanical impedance of human joints is made using stiffness (or compliance) ellipsoids (Mussa-Ivaldi et al., 1985). To modulate the stiffness ellipsoids of the end effector during tasks, it appears that the limb geometry and posture have a large impact on its shape and orientation, whereas the voluntary co-contraction of the participating muscles affects mainly the size of the ellipsoids (Milner, 2002;Perreault et al., 2002). To quantify the mechanical impedance properties of multi-DoF human joints, the traditional approach consists in applying position or force perturbations to the joint of interest and measuring the displacement force or position responses (Mussa-Ivaldi et al., 1985). Yet, this process is complex and long. Investigations are made to simplify and speed up the procedure (Masia et al., 2012;Masia and Squeri, 2015). The measurements usually concern the endpoint impedance (Mussa-Ivaldi et al., 1985;Masia et al., 2012). Yet, there is mechanical impedance data of specific human joints, such as the wrist (Formica et al., 2012;Pando et al., 2014) or the ankle (Kearney et al., 1997;Lee et al., 2014). There is, for now, no data on the mechanical impedance of more complex joints such as the shoulder or hip. Yet, efforts are made toward this direction with the design of low inertia shoulder exoskeleton to measure neuromuscular properties (Hunt and Lee, 2019).
The observations on how human beings modulate the stiffness ellipsoids were successfully exploited to implement a teleimpedance controller to generate human-like motions and desired task space impedance (Ajoudani et al., 2018). Similarly, these observations could be used to reduce the mechanical complexity of human-like VSAs. Indeed, if the human-like compliant behavior is directly implemented in the mechanism, only one additional actuator is ideally required to modulate the overall size of the stiffness ellipsoid, in analogy to the voluntary co-contraction of human muscles. As a matter of fact, from the study of tendon mechanisms (e.g., the Salisbury Hand), we know that to control N DoF in position only N+1 actuators are required. There is then an internal tension that can be used to modulate the stiffness. This was shown for a grasping hand in Mason and Salisbury (1985).
In this work, we propose the concept and modeling of a configurable architecture of a physically compliant 2 DoF VSA that can match the passive behavior of either a human wrist or a human ankle. This joint can control its orientations around 2 perpendicular axes and its overall stiffness as a one-dimensional parameter similar to the co-contraction of human muscles. It is based on a state-of-the-art parallel manipulator that we modified to implement an inherent compliance with a compact architecture. Hence, the variable stiffness (VS) ability of the system relies on the antagonistic arrangement of nonlinear elastic actuators (Vanderborght et al., 2013;Wolf et al., 2016). The compact architecture of the system is based on the PM architecture with the minimum number of required motors to get a VSA for 2 DoF in position (i.e., 3). Thanks to its kinematics, the proposed system could be used as a wrist or ankle joint. It can also be part of more complex joints such as the hip or shoulder joint.
Indeed, these joints are represented with 3 DoF as ball-and-socket joints and one of their three rotations is around a longitudinal axis of the body segment attached to the joint (Wilhelms and Gelder, 2001). Therefore, we can design them with a hybrid architecture of a PM with 2 DoF in series with a 1 DoF rotator, like for instance the humeral rotator proposed in Mazzotti et al., (2015), which can be integrated with a 2 DoF shoulder joint presented in Troncossi et al., (2009). A similar idea could be proposed for a hip joint with a femoral rotator. Therefore, our proposed system could also be used as part of these more complex joints. Figure 1 shows the general principles of the system and some of its potential applications. FIGURE 1 | The proposed 2-DOF variable stiffness joint allows shaping the stiffness (or compliance) ellipsoid at different postures in the design phase, while the volume can be varied in real-time to make the overall behavior stiffer or softer. The kinematics of the 2-DOF VS joint can be used for the design of wrists and ankles of human-like bionic systems (in orange). In principle, this joint can also be used as part of more complex designs for a shoulder or hip human-like joint (in yellow).
Besides the general description of the proposed system in Section 2.1, the main contributions of this paper are: (i) to prove the human-like VS ability of the system with a compact architecture (Section 4.1), (ii) to study the effect of the main design parameters of the system on its passive compliance to provide general guidelines (Section 4.2), and (iii) to provide sets of design parameters to match the passive compliance of a human wrist or ankle based on an optimization process (Section 4.3). Other contributions of this paper include improvements on the kinematic analysis of the rigid structure (used as a basis of our mechanical architecture) done in Sofka et al., (2006a) and the kinetostatic analysis of our proposed compliant structure (Section 3).
The paper is organized as in the following: Section 2 details the design concept of the proposed system and the methodology of its study. It includes, in particular, a description of its general architecture and its main design parameters. Section 3 provides the model and the kinetostatic analysis of the system. Section 4 gives the main results of the current study that are then discussed in Section 5. Finally, Section 6 summarizes the main contributions and concludes the article.

General Architecture of the System
Based on the state of the art of compact 2 DoF in orientation joints based on PMs (refer to (Bajaj et al., 2019) for a quite extended list), we selected the kinematic architecture of the Omni-Wrist III, designed by Rosheim and Sauter (2002) and studied by Sofka et al., (2006a) and Sofka et al., (2006b). This design has several advantages such as 180°hemispherical singularity free movement (Rosheim and Sauter, 2002). In theory, only two legs are needed to actuate the system but the main version has 4 legs. Its general architecture is shown in Figure 2A. In our design, we propose to use a 3-legs configuration with nonlinear elastic elements added in the serial arrangement of each leg. Figure 2B presents its general architecture. The system is a PM composed of a fixed frame (in dark gray on the scheme) and a moving platform (in light gray), i.e., the end effector, linked together with three independent and identical kinematic serial chains, called respectively leg A (in blue), leg B (in green) and leg C (in red). Each leg is composed of the serial arrangement of one motor (M p ), one nonlinear elastic element (S p ), and one kinematic chain (KC p ). The kinematic chains, based on the geometry proposed in Rosheim and Sauter (2002), are composed of 3 linkages each connected together and to the end effector by a succession of 4 noncoplanar revolute joints. Based on the modeling of Sofka et al., these KCs are parametrized by α, which is the angular deviation of the link 2 (in red), and by a characteristic length d (Sofka et al., 2006a). Figure 2C shows the main structure of one KC with its parameters.

Variable Stiffness Ability
As it is common in VSA with antagonistic architecture, to get the variable stiffness ability, the stiffness characteristic of the elastic element has to be nonlinear. There are many ways to design nonlinear springs, to cite a few, the physical implementation can rely on tendons and linear springs (Catalano et al., 2011), on a combination of belts and linear springs (Lemerle et al., 2019), on a combination of torsional spring with guide shafts with varying radius (Koganezawa et al., 2005;Koganezawa, 2017), or on hydraulic or pneumatic solutions (Stoeffler et al., 2018).
Nevertheless, nonlinear stiffness alone is not sufficient to grant the variable stiffness control of the system, which is equivalent to prove that we can control independently the equilibrium orientations of the system and its overall stiffness along one dimension. To prove this analytically, the kinetostatic analysis of the system is first presented in Section 3. Then, the final arguments to show the VS control of the proposed system are given in Section 4.1.  (Rosheim and Sauter, 2002) is composed of 4 legs linking a fixed frame and a moving platform. Each actuated leg is composed of the serial arrangement of one motor (M * ) and one kinematic chain (KC * ). Each non-actuated leg is composed of one kinematic chain (KC * ) (B) Represents the architecture of the proposed system. We use a configuration with 3 legs only. In each leg, we add elastic elements (S * ) between the motor (M * ) and the kinematic chain (KC * ). (C) shows the details of one kinematic chain (each KCs are identical). It is composed of 3 linkages (1, 2 and 3) and 4 noncoplanar revolute joints (I, II, III, IV). The joint IV is linked to the end effector and the joint I is linked either to the output of the motor in the Omni-Wrist III or to the output of the elastic element in our proposed architecture.
Frontiers in Robotics and AI | www.frontiersin.org March 2021 | Volume 8 | Article 614145 In addition, we would like to show the effects of the modulation of the stiffness on the overall compliance of the system. To represent them graphically, we need to derive a model of the system (done in Section 3) and to model the elastic elements. In this study, we make the deliberate choice of adopting a hyperbolic sine spring, whose mechanical characteristics have the form where f S* stands for the force exerting by the spring, K * for the spring constant, δ * the deflection of the spring, and δ 0* a characteristic deflection of the spring. This choice is motivated by the fact that springs with a similar characteristic have been successfully used in other actuators (Catalano et al., 2011;Lemerle et al., 2019), and recent works suggested that this type of elastic behavior replicate several attributes of an antagonistic pair of muscles driving a joint (Garabini et al., 2017). This restricts the conclusions drawn in the next sections to be exact only when dealing with this specific kind of springs. Nevertheless, it is expected to observe similar trends in all springs with similar stress-stiffening characteristics. For the sake of clarity, the detailed graphical representation methodology is described after the complete analysis of the system in Section 3.8. The effects of the modulation of the stiffness are shown in Section 4.1.

General Design Guidelines
As previously described, there are various ways to implement the nonlinear elastic elements. In a first approach, we leave this choice to the designer to focus on the design of the KCs and their arrangement. Therefore, as general design guidelines, we study the effects of the purely geometric parameters (simply referred to as geometric parameters) on the compliance of the proposed system. To do this study, the compliance of the system should be modeled as a function of these geometric parameters. This is done in Section 3. The first step in the modeling is to define the geometric parameters. They correspond to the parameters of the KCs α and d (Sofka et al., 2006a), and the angular arrangement of the first joint of each leg around the basis, namely η A , η B , η C (refer to the parametrization of the system in Section 3.1 for the exact definitions of these parameters). The angular positions of each leg are defined according to the scheme shown in Figure 3.
In the case of a fixed distribution of the legs around the basis, α and d are not enough to modulate the compliant behavior of the system (refer to Section 4.2 for more details). That is why we need to extend the modeling done in Sofka et al. (2006a) by adding the angular positions of the legs.
Moreover, to study separately the effects on the compliance due to the kinematic architecture from the effects due to the internal variation of the compliance of the elastic elements, the compliance characteristic of each leg is modeled as a constant. Concretely, to study the impact of the angular position of the legs, two types of configuration are simulated. Firstly the leg A is considered as fixed with η A 0 and the couple (η B , η C ) is varying symmetrically, with η B in {π/4, π/3, π/2, 2π/3, 3π/4} and η C −η B . Secondly, the legs A and C are fixed, respectively equal to π/4 and −π/4 and η B is varying in the set {−3π/4, −2π/3, −π/2, −π/3, π/3, π/2, 2π/3, 3π/4, π}. To study the impact of the angular deviation of the middle link of each leg α, this angle is varying in the set {π/6, π/4, π/3, 5π/12}. Finally, to study the impact of the characteristic length d of each leg, this length is varying in the set {25 mm, 50 mm, 75 mm, and 100 mm}. All the simulations and computations are performed using MATLAB. The results are shown in Section 4.2.

Design Parameters to Match Compliance of Human Wrist and Ankle
One of the main goals of this study is to show that we can implement passive compliant behaviors of human joints. In this work, we focus on two different human joints: the wrist and the ankle. To find the appropriate sets of design parameters to match their passive compliant behavior, we derived an optimization process. In the first step in Section 2.4.1, we need to define the behaviors we would like to embed based on the mechanical impedance data of human joints. Then, in Section 2.4.2 we define the design parameters of the proposed architecture. Finally, to derive the cost function of the minimization algorithm, we need to characterize the compliance of the system as a function of the design parameters. This is done in Section 3 and the cost function is defined in Eq. 58. Results of the optimization process are shown in Section 4.3.

Parameters of Human Wrist and Human Ankle
Based on the state-of-the-art, we selected the biomechanical data of a human wrist presented in Pando et al. (2014) and of a human ankle presented in Lee et al. (2014).
In Pando et al. (2014), Pando et al. studied, among others, the passive wrist stiffness in the two dimensions that are the flexion/ extension and the radial/ulnar deviation. They extracted from their data the stiffness ellipses around the neutral position of the wrist in case of no voluntary co-contraction of the subjects to evaluate thus the passive stiffness of the articulation. They found that the human wrist stiffness is anisotropic and the direction of the greater magnitude of the average stiffness ellipse is around 12.1°from the pure radial deviation in the counter-clockwise direction toward a pure flexion. The ratio between the magnitude of the two axes of the stiffness ellipse is about 2.69 and using the area of the ellipse [5.61 (Nm/rad) 2 ], we can extract the lengths of the two semi-axes of the ellipse, which are then 2.19 and 0.81 Nm/rad.
In Lee et al. (2014), Lee et al. studied the mechanical impedance of a human ankle with relaxed muscles, i.e., in passive conditions, both in seated and standing positions. They found that for low-frequency motion the stiffness dominates the mechanical impedance of the ankle, which is our application case as we consider the static stiffness in our design. Its characteristics in standing positions are: a tilt angle of 10°, a magnitude in the major principal axis direction of 23.86 Nm/rad, and a magnitude in the orthogonal direction of 8.14 Nm/rad (the ratio is 2.99).
Considering these data, we can derive the average Cartesian compliance ellipse and the associated compliance matrix for both a human wrist and ankle, using similar parametrization as the one proposed in Section 3.7.

Design Parameters Definition
In a second step, we should define all the design parameters of the proposed architecture. They can be divided into geometric parameters (previously described), elastic parameters, and internal modulation parameters. As a reminder, there are 5 geometric parameters: α, d, η A , η B and η C .
Each nonlinear spring can be characterized by nonlinear functions. Their explicit formulations depend on mechanical implementation. They can have various forms such as polynomial functions of degree higher than 2 (e.g., quadratic springs), trigonometric functions, exponential or hyperbolic functions to name but a few (English and Russell, 1999;Koganezawa et al., 2005;Vanderborght et al., 2009;Catalano et al., 2011;Petit et al., 2015;Lemerle et al., 2019). Therefore the number of their associated parameters depends on the complexity of the model. These parameters are referred to as elastic parameters. Noting n elastic the number of parameters for one nonlinear spring, there are in total 3 × n elastic elastic parameters in the proposed system. They can be highly variable depending on the mechanical implementation. Moreover, it is not easy to derive realistic values of them to fix the bound constraints for the minimization problem. Therefore, it is better if the solution provided by the algorithm is robust or independent with the values of these parameters.
The internal modulation parameter, called λ, stands for the onedimensional quantity we can use to control the stiffness in the proposed system. Unlike the previous parameters, λ is not fixed to one mechanical implementation and can therefore be different if we try to match several compliant behaviors. Therefore, there is one additional parameter for each desired behavior.
Therefore, the number of design parameters, noted n param , to match the specified targeted behaviors is defined as where n elastic stands for the elastic parameters of one leg and n targets the number of compliant targeted behaviors. For this study, we model the nonlinear elastic elements as hyperbolic sine springs for similar reasons as previously described in Section 2.2. Their mechanical characteristics have the form defined in Eq. 1. Therefore, there are 2 elastic parameters for each leg: (i) the spring constant K * and (ii) the free length δ 0* . All the design parameters are summarized in Table 1.

MODEL AND ANALYSIS
In this section, we describe the kinetostatic analysis of the proposed system with two goals in mind: (i) to derive the compliance matrix as a function of the design parameters and (ii) to show the variable stiffness ability of the system.
Firstly the kinematic analysis done by Sofka et al. (2006a) needs to be extended to include the angular arrangement of the legs. Then, we describe the inverse and differential kinematics of the system in a closed form, showing that it is a 2 DoF joint. We provide an explicit and complete formulation of the general mapping of the generalized pose of the end effector as a function of a minimal parametrization. This result is then used to analyze the static equilibrium of the system and describe its static stiffness and Cartesian compliance, required for the optimization algorithm. A list of the main symbols used in the analysis can be found in Table 2.

Parametrization of the System
Define the reference coordinate frame attached to the fixed base of the mechanism where O 0 marks its origin and (X 0 , Y 0 , Z 0 ) form an orthonormal basis, noted also B 0 . Define also, with analogous notation, the reference frame attached to the end effector of the mechanism C ee (X ee , Y ee , Z ee , O ee ).
A possible choice of parameters that easily parametrizes the generalized pose of the end effector is where the first three components stand for the Euler angles of the end effector, defining, therefore, its angular orientation about respectively the axis X 0 -Y 0 -Z 0 , namely (roll, yaw, pitch) and the last three components are the position of its center O ee in B 0 . Denoting with Rot(z, θ) and Trans(z, d) the homogeneous matrices associated to the transformations of a simple rotation of θ around the axis z → , and of a simple translation of d along the axis z → , respectively, it is possible to write i.e., 0 T ee ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ c z c y c z s x s y − s z c x c z s y c x + s z s x x ee s z c y s z s x s y + c z c x s z s y c x − c z s x y ee −s y s x c y c x c y z ee where, c k and s k stand respectively for cos(α k ) and sin(α k ), for k in {x, y, z}.
Considering now each leg (A, B, and C) of the mechanism, it is possible to parametrize the configuration of C ee also by following the associated serial kinematic chain, e.g., using the Denavit-Hartenberg (DH) convention (Siciliano and Khatib, 2007). Note first that the three legs share the same kinematic structure, which is shown in Figure 4. As such, they can all be parametrized with the same DH table given in Table 3. Note that, by the construction of the system, linkages 1 and 3 are identical, and d and α are fixed parameters that define the shared legs geometry, while η is different for each leg.
Calling (q 1 , q 2 , q 3 , q 4 ) and η the joint angles and the base angle of the generic leg, the global transformation matrix going from the system base to the end effector through the leg 0 T l ee is defined as where all the intermediate matrices are defined based on the parameters in Table 3. After calculations, we obtain: where all the non-trivial terms are defined in Supplementary Appendix A.

Closing the Kinematic Chain
When the joint is mounted, i.e., the legs are assembled, the assembly constrains the final frame associated to each leg (and the generic frame 0 T ee ) to be equal, i.e., Orthonormal basis of the frame C ee x Generalized pose of the end effector α x α y α z T Angular orientations of the end effector x ee y ee z ee T Cartesian position of O ee i T j Transformation matrix from C i to C j i T lj Transformation matrix from C i to C j through leg l ∀i ∈ {1, 2, 3, 4}, q i Angular position of i th joint for a generic leg α Angular deviation of the middle link of each leg d Characteristic length of each leg η Generic angular orientation of one leg η l Specific angular orientation of leg l (g i ) i∈{1,...,7} and (h i ) i∈{1,2,3} families of functions defined respectively in Eqs. (11) and (12)

Positions of the actuated joints τ act
Torques of the actuated joints f IK : Torque generated by the system ∀* As it is well-known, the solution of Eq. 8, is fundamental to yield both the direct and inverse kinematics of the parallel system. An important aid to facilitate this comes from Sofka et al., (2006a), where it was found that a feasible mounting of each leg holds for Therefore, after simplifications, Eq. 7 can be written as where is a family of functions defined as and (h i ) i ∈ {1,2,3} is a family of functions defined as The families (g i ) i ∈ {1,...,7} and (h i ) i ∈ {1,2,3} are introduced to improve the readability of the matrices and the calculations.
As expected from Sofka et al. (2006a), and as we will confirm with our derivation, the assembled system has 2 DoF. Thus, it is possible to parametrize the position and orientation of the reachable end effector configurations with only two variables. For convenience we choose To fully characterize the set of feasible end effector configurations it is necessary to specify the other 4 variables in Eq. 3.
To derive the angular orientation α x as a function of u, a key point relies on noting that 0 T lee (3, 2) 0 T lee (2, 3), and then, by using Eq. 5 and rearranging the terms, it is possible to extract α x as α x arctan s y s z c y + c z , a result in agreement with (Sofka et al., 2006a). Then, the analysis can be extended by deriving an explicit formulation of the position of the center of the end effector (x ee , y ee , z ee ) as a function of u. Only the key results are given in the following part and the reader is invited to refer to Supplementary Appendix B for more detailed proofs of the different assumed expressions.Based on Eqs 5 and 8, we have FIGURE 4 | Kinematic diagram of a generic leg. It is composed of 3 linkages (1, 2, and 3), and 4 noncoplanar joints. It is parameterized by its angular deviation α, its characteristic length d and its angular orientation η. The coordinates frames are also represented. Its associated DH table is given in Table 3.
Frontiers in Robotics and AI | www.frontiersin.org March 2021 | Volume 8 | Article 614145 x ee 0 T lee (1, 4) (16a) y ee 0 T lee (2, 4) (16b) z ee 0 T lee (3, 4), (16c) from which it can be shown [based on (10)] that which gives the distance square of the point O ee to the point O 0 . This distance is thus independent of the positions of the legs as it depends only on d and α which are two parameters inherent of the mechanical design of the system and could also be found based on a geometric closed form of the system. Defining such distance as and using the previous equations, it can be shown that According to the convention defined previously, we have ∀ α y , α z ∈ ( − π/2; π/2), x ee ≥ 0, therefore, using Eq. 19a, we obtain Then, the system defined by Eqs. 19b and 19c has 4 pairs (y ee , z ee ) of solutions. Yet, according to the convention previously defined, we know that y ee should increase when α z is increasing and z ee should decrease when α y is increasing. Therefore, there is only one suitable solution for the pair (y ee , z ee ), which is To conclude, the proposed system has 2 DoF in orientation and the explicit formulation of the general mapping T is as This mapping is always defined for all (α y , α z ) in (−π/2; π/2).

Inverse Kinematics Model
This subsection aims to find the position of all the joints of each leg of the system as a function of the pose of the end effector, x. Based on the chosen parametrization, calculations for each leg are very similar. Therefore, only generic calculations are presented in the following section. Moreover, using the reduction of parametrization, presented in Eq. 9, the aim is equivalent to find only q 1 and q 2 as functions of x. The reader is invited to refer to Supplementary Appendix C for a more detailed demonstration of the following equations.
Firstly, to obtain q 2 , the idea is to isolate h 3 (q 2 , α) defined in Eq. 12. It can be shown that c 2 is defined by the equation From Eq. 23, it is possible to extract q 2 as a function of x and the parameters of the system, as follows: Besides that, it can be shown that c 1 and s 1 should satisfy the following system whered is defined as As q 2 is known in Eq. 24, it is possible to extract from Eq. 25 q 1 as a function of only x and the parameters of the system, as follows: The complete solution for a generic leg is then It can be easily adapted for each leg, by replacing the generic parameters and variables with the proper ones. For instance, by replacing in the previous equations, (q 1 , q 2 , q 3 , q 4 ) with (q A1 , q A2 , q A3 , q A4 ) and η with η A , it is possible to obtain the Frontiers in Robotics and AI | www.frontiersin.org March 2021 | Volume 8 | Article 614145 inverse kinematics of leg A, namely T A ee , and analogously for legs B and C.

Differential Kinematics
The explicit inverse kinematics model is given as where q act represents the vector of the actuated joints of the system, defined as the function f IK can be written as and u stands for the minimal parametrization of the system, defined in Eq. 13.By differentiation of Eq. 29, we obtain where the matrix A u is defined as where for any function f and variable x i , z x i f stands for the partial derivative of f with respect to x i . An explicit formulation of A u can be computed based on the results presented in Section 3.3 and is reported in Supplementary Appendix D. This formulation depends explicitly on the pose of the end effector.

Static Equilibrium
Using the kinetostatic duality, we obtain where τ u in R 2 represents the torques of the end effector associated with u and τ act stands for the torques of the actuated joints, i.e., it could be denoted as For a given external wrench f, acting on the end effector, the static equilibrium gives us, using Eq. 34, the relationship between the external wrench f and the torque of the actuated joints as As we assume to be working outside of the singularities, the rank of the matrix A u is 2, so it has a pseudo-inverse. So, it is possible to derive from Eq. 36, the torque τ * act required to balance a given external wrench f * acting on the end effector, as where λ is a scalar real and N 0 represents a basis of the solutions of The reader is invited to refer to Section 4.1 for details on the calculation of N 0 .

Static Stiffness
To define the static stiffness, we first need to define τ act as a function of the torque of the elastic elements. Define f S A , f S B f S C the torque functions of the nonlinear elastic elements placed on each leg of the system. Letting * ∈ {A, B, C}, f S * represents thus the torque of the actuated joint on the leg * and is a function of the deflection δ * of the elastic mechanism defined as where θ M* denotes the position of each motor M * . θ MA , θ MB and θ MC are the controlled input variables of the proposed system. Therefore, τ act can be written as where δ stands for the vector of the deflections of the elastic elements defined as Let's define σ ee , the static stiffness of the end effector as Using the chain rule we have that σ ee zτ u zτ act zτ act zq act zq act zu .
Based on Eq. 34, we have and based on Eqs. 29 and 33, we have Moreover, we have where K(δ) is a diagonal matrix in R 3×3 with the stiffness of each leg on its diagonal, defined as Therefore, we obtain

Cartesian Compliance and Cost Function
To represent the stiffness characteristics of the system, we decide to represent in the space its Cartesian compliance, the quantity denoted by c in R 3×3 such that dp cdf , where p in R 3 is a point attached to the end effector and f in R 3 is the Cartesian external forces applied at this point, i.e., f represents only the forces applied on p in the directions x, y and z. This decision is motivated by the fact that, as shown in the previous section, the end-effector of the system has only 2 Dof, which means that the system will not allow motions in one of the three Cartesian directions of space, thus in that same direction the stiffness would be infinite. Therefore its inverse (i.e., the compliance) will be more manageable mathematically, and easier to represent graphically. According to the minimal parametrization defined in Eq. 13, p is a function of u. So, we have dp Jdu, where the matrix J is in R 3×2 and denotes the Jacobian matrix of the function p with respect to the variables of u.
where τ u denotes the forces associated with u and f stands for the Cartesian forces of p, from which we obtain in the static equilibrium Moreover, using the definition of the static stiffness of the end effector in Eq. 42, we obtain dτ u σ ee du. (53) The matrix K is always invertible and as previously said, the rank of the matrix A u is 2, as we assume to be working outside of the singularities, so it has a pseudo-inverse. Therefore, based on Eq. 48, σ ee is invertible. So, So, by injecting Eq. 52 in Eq. 54, we obtain which gives us, based on Eq. 50, the expression of the Cartesian compliance c such as dp cdf with c Jσ −1 ee J T .
Knowing the compliance matrix, it is possible to derive the cost function used in the optimization algorithm such as where c ref stands for the desired compliance matrix based on the definition of the desired behavior and n targets the number of desired compliant behaviors. As a note, the sum is only done on the lower-diagonal and diagonal components (and non-null) of the compliances matrices as they are symmetric. The minimization problem is done using MATLAB, based on the sqp-algorithm.

Graphical Representation Methodology
To get a visual representation of the Cartesian compliance, we plot the ellipsoid associated with the matrix c. This ellipsoid is defined using the singular value decomposition of c. Its axes are directed by the eigenvectors of the matrix and their half lengths are defined by their respective associated eigenvalue. The ellipsoid is centered in p. Hence, the ellipsoid represents the possible deflection of the system in the space for a normalized external wrench, and thus the greater the ellipsoid is, the more compliant is the system. Figure 5 shows an example of the various representations and views of the compliance ellipsoids for a specific configuration of the system. One way to represent the compliance ellipsoids is to draw them in 3 dimensions. They can be viewed either from a 3D perspective, such as in Figure 5A, or from a top view of the system [i.e., on a projection on the plane (y, z)] such as in Figure 5B. A second option is a representation in 2D, where the same ellipsoids of the Cartesian compliance are plotted in a plane (L y ,L z ), such as in Figure 5C. This plane corresponds to an arbitrary scaled distance between the points while following the surface of all of the reachable points. It represents therefore the ellipsoids as if they were always seen from the normal of every point. Note that these ellipsoids are not defined in the (L y ,L z ) plane, so not to be affected by the representation singularities associated with the selected parametrization inL z ± π/2. On the contrary, they are defined in the tangent plane to each point of the curved surface in Figure 5A. Note also that forL z ± π/2 the ellipsoids are arbitrarily oriented choosingL x 0. Note that we will preferentially use this representation in the rest of the paper, referring to it as flattened representation.
The analogous of meridians and parallels have been plotted also to help the understanding of the corresponding ellipses in the example figures. In each representation, the locus of all the possible points reachable of p is also plotted as a colored area even if the ellipsoids are displayed only in some discrete points. This area is referred to as the workspace of the system. The flattened representation has the advantage not to distort the ellipsoids and not to be sensitive to the Cartesian workspace which may be affected by the design parameters of the mechanism. So, this representation will mainly be used thereafter. Finally, the angular orientation of the axes of the first joint of each leg is plotted as small black lines in the different graphs.
Besides, one thing to notice is that the ellipsoids are only ellipses, tangent to the workspace. This is in agreement with the consideration made at the beginning of this subsection, which motivated the choice of representing the system compliance instead of its stiffness.

Variable Stiffness Ability of the System
In this section, we prove that it is possible to have a decoupled control of the position and the stiffness of the system, i.e., that it is possible to modify the spring load without moving the end effector.
Let us consider first the static stiffness given in Eq. 48 and show that it can be modulated. For a given position u, the matrix A u is decided. Moreover, in case of the elastic elements are nonlinear, the matrix K depends on the deflection δ which can be controlled through the position of the motors θ M* . Therefore, the matrix K and thus the static stiffness of the system σ ee can be controlled using nonlinear elastic elements in series of each kinematic leg. This result agrees with the one in case of the design of one DoF physically compliant VSA using the antagonistic approach, as reported in the literature (Vanderborght et al., 2013).
Similarly, the Cartesian compliance of the system defined in Eq. 57 can be controlled if the elastic elements are nonlinear by modifying the positions of the motors θ M * .
Let us prove now that the compliance can be modulated at a fixed position. The first step is to derive a possible expression of N 0 , a basis of the space of the solutions of Eq. 38. The matrix A u , defined in Eq. 33 can be written as where Q y and Q z stand respectively for the first and second columns of matrix A u . Therefore, X in R 3 is a solution of Eq. 38 if and only if This system of equations has a space of solution of dimension 1 if and only if Q y and Q z , i.e., the two columns of A u , are noncollinear which can be written as det(A T u A u ) ≠ 0. This is always true as we assume to be working outside of the mechanical singularities, i.e., the rank of A u is equal to 2.
Therefore, the Solutions of Eq. 38 are described as A possible definition of N 0 is thus where . × . and . 2 denote respectively the cross product and the euclidean norm of vectors in R 3 . Therefore, along the space described by N 0 , we have a one dimensional space where τ act is varying without affecting the equilibrium position of the end effector. The matrix A u is thus constant as the pose of the end effector is decided and therefore, q act is fixed. Assuming the non-linearity of the elastic elements, there is a space to control the compliance of the end effector without modifying the pose of the end effector. This is a one-dimensional space and can be therefore used to control a one-dimensional variable linked to the stiffness, such as the sum of the potential elastic energies stored in each elastic elements.
The effect of the internal modulation is studied while varying the parameter λ defined in Eq. 37, with the model of the elastic elements defined in Eq. 1. It has to be noticed that with this definition of f S* , the stiffness of each elastic element is an even function in λ, as f S* is an odd function in λ, using Eq. 37. Figure 6 represents the Cartesian compliance ellipses for various values 1 of λ. We can see how, acting on λ, the ellipses scale homogeneously in all directions, without appreciable changes in shape and orientation like the cocontraction of human muscles (Milner, 2002;Perreault et al., 2002).
Moreover, Figure 7 shows the evolution of the legs torque and stiffness, and the output force at three different positions of the end effector, highlighting how the actuation of λ changes the stiffness but does not change the external equilibrium. Indeed, the external wrench balanced at the end effector (f ay , f az ) is not modified whereas the stiffness of each leg is affected and thus the stiffness of the end effector.
As a final note, it is possible to show that the stiffness of the system is increasing when the external wrench is increasing. In addition, although the stiffness increases slightly in all directions, the increase is larger in the direction of the applied external force. For sake of space, the associated figures are not displayed here as this is an expected behavior for a VSA (Vanderborght et al., 2013).

General Design Guidelines
In this section, we give the results of the effects of the geometric parameters on the compliance of the system according to the methodology described in Section 2.3. To show their effects, the Cartesian compliance ellipsoids are plotted at the center of the end effector (i.e., p in the previous analysis is O ee ) for several positions of the end effector [parametrized by u α y α z T ] and various geometries defined by the design parameters. When it is relevant, the workspace of the system is also plotted as a colored area even if the ellipsoids are displayed only in some discrete points. Figure 8 represents the different figures of the simulations when the various design parameters are varying. Figures 8A,B show the results when the angular positions of the legs are varying in a symmetric or asymmetric way. We can see that in both cases, depending on the angular position, it is possible to stretch the ellipses in one of its two main directions by modifying the value of the angle and to rotate them. Figures 8C-E show the effects of α on the Cartesian compliance ellipses and the workspace. We can notice that the compliance ellipses are globally increasing when α is increasing. Figures 8F-H show the effects of d on the Cartesian compliance ellipses and the workspace. We can notice that the compliance ellipses are globally increasing when d is increasing.
As it is shown in the different figures, all the design parameters do not have the same effect nor the same intensity on the compliance characteristic of the system. For instance, the effects of α are more important in terms of FIGURE 6 | System Compliance in flattened representation, as a consequence of the internal modulation while λ is varying in the set {0 Nm, 2 Nm, 5 Nm, 10 Nm}. Simulations are made at the center of the end effector, with η A , η B and η C evenly distributed around the basis, α π/4, d 50 mm, and all the elastic elements identical (K * 1 Nm/rad and δ 0* 1 rad).
stretching the ellipses, especially at the borders of the workspace, whereas d has more a global influence everywhere that can be compensating by using different elastic elements. Moreover, it is important to notice that both previous results go together with a stretch of the surface of the feasible positions of the end effector, shown in Figures 8C,D,F,G. This explains that the increase in compliance is also an effect of the increase of the lever arm of the forces. Combining the previous two effects, it is conceivable to change the shape of the ellipses while maintaining the feasible surface and ellipses area roughly constant by changing alpha and d simultaneously in opposite directions. However, it appears that α and d are not enough to modulate the ratio of the ellipse at (0,0) position where we observe always a circle for an even distribution of the legs as shown in Figures 8E,H. This shows a limitation of using only these two parameters to modify the embedded behavior in the system and justify the use of the angular positions of the legs as geometric parameters.

Design Parameters for Human-Like Joints
In this section, we derive sets of parameters to match the passive compliant behavior of a human wrist (Pando et al., 2014) and a human ankle . In both cases, the idea is to find the best set of parameters to match the targeted ellipse at the neutral position, defined in Section 2.4.1.
Regarding the input of the algorithm, we would like to match 3 characteristics. Therefore, with 12 design parameters the solutions are not unique. We can then set additional constraints besides the realistic values of the design parameters. As we are considering passive behavior, we set λ equal to 0. This means that in absence of external load, no motor torque is required to get the passive behavior of the FIGURE 7 | Evolution of the torque of each leg τ leg , their respective stiffness σ leg and the force balanced by the end effector (f ay , f az ) (in the directions of α y and α z ) for λ varying between -10 Nm and 10 Nm and for various equilibrium positions u eq and external forces f eq (f ay ext , f az ext ). (A) corresponds to u eq (0, 0) deg − f eq (0, 0) Nm, (B) to u eq (90, 0) deg − f eq (5, 0) Nm and (C) to u eq (60,30) deg − f eq (3, −5) Nm.
system. Moreover, we consider an even distribution of the legs around the basis. Figure 9 represents solutions for a wrist and an ankle considering theses constraints. Their associated design parameters are reported in Table 4. Therefore, it is possible to get feasible systems based on the proposed architecture with a passive compliant behavior similar to the human wrist or a human ankle.

DISCUSSION
In this work, we proposed and studied a new concept to design a 2 DoF VSA. To reduce the complexity of the mechanical implementation and increase the compactness of the system, we use only 3 motors. This is the minimum number of motors to be able to have a 2 DoF mechanism fully FIGURE 8 | Effects of the geometric parameters on the system compliance. (A) represents the system compliance in flattened representation, as the angular positions of the two legs η B and η C vary in a symmetric way. η A is fixed and equal to 0, η C −η B and η B is varying in the set {π/4, π/3, π/2, 2π/3, 3π/4}. Simulations are made at the center of the end effector, with α π/4, d 50 mm and c (1 rad/Nm 1 rad/Nm 1 rad/Nm). (B) represents the system compliance in flattened representation, when the angular position of one leg (η B ) varies in an asymmetric way. η A and η C are fixed and respectively equal to π/4 and −π/4 and η B is varying in the set {−3π/4, −2π/3, −π/2, −π/3, π/3, π/2, 2π/3, 3π/4, π}. Simulations are made at the center of the end effector, with α π/4, d 50 mm and c (1 rad/Nm, 1 rad/Nm, 1 rad/Nm). (C) and (D) represent the workspace of the system in 3D view and side view when α is varying in the set {π/6, π/4, π/3, 5π/12}. (E) is the associated flattened representation of the compliance ellipsoids. Simulations are made at the center of the end effector, with η A , η B and η C evenly distributed around the basis, d 50 mm and c (1 rad/Nm, 1 rad/Nm, 1 rad/Nm). (F) and (G) represent the workspace of the system in 3D view and side view when d is varying in the set {25 mm, 50 mm, 75 mm, and 100 mm}. (H) is the associated flattened representation of the compliance ellipsoids. Simulations are made at the center of the end effector, with η A , η B and η C evenly distributed around the basis, α π/4 and c (1 rad/Nm, 1 rad/Nm, 1 rad/Nm).
Frontiers in Robotics and AI | www.frontiersin.org March 2021 | Volume 8 | Article 614145 actuated in positions and partially actuated in stiffness. Therefore, our solution is more compact than most of the solutions existing in the current state-of-the-art (refer for instance to the works of (Grebenstein et al., 2011;Romano et al., 2014;Weckx et al., 2014;Koganezawa, 2017;Tan et al., 2017)). To the best knowledge of the authors, only the work in (Stoeffler et al., 2018) follows the same approach of using only 3 motors. Obviously, a limitation of our proposed solution relies in the loss of DoF in the control of stiffness of the system with respect to the existing 2 DoF VSAs. Indeed, our system can control the overall size of its endpoint stiffness ellipse but not its stiffness independently in the two directions of motion. Yet this approach is motivated by the fact that human beings modulate their stiffness through both their posture (exploiting the geometry and properties of their limbs) and co-contraction (Milner, 2002;Perreault et al., 2002). More interestingly, their effects are quite distinct as the limb geometry and posture have a large impact on the shape and orientation of the ellipse and the co-contraction on the overall size (Perreault et al., 2002). Therefore, by implementing physically the human joint geometry and allowing the control of the overall size of the stiffness, we can design human-like joints with the potential of natural control of both the position and stiffness. To the best knowledge of the authors, there is no other work following this approach.
However, the challenge of the mechanical implementation of the passive compliance of human joints remains. That is why we derived a model of our system and studied the effect of its main design parameters on the characteristics of the compliance ellipses. Additionally to the guidelines given in Section 4.2, it appears that the influences of the angular positions of the leg are interesting as they allow us to stretch or rotate the compliance ellipses in certain directions. Combining all of the effects, it will be possible to embed some specific behavior in the design for a specified application. And indeed, we derived sets of design parameters to match the passive behavior of a human wrist and a human ankle.
For now, biomechanical data exists only at the neutral position of the joints. However, in the future, additional data could be available on various positions of the workspace of the joints. In this case, we may want to match several ellipses. In the proposed system, it is possible to increase the number of desired ellipses. Figures 10A,B present solutions for 2 desired ellipses. The case 1.1 corresponds to a solution with an additional constraint of the vector λ equal to 0. As it is shown, this solution may be insufficient. Case 1.2 stands for a solution without any FIGURE 9 | (A) Compliance ellipses at the neutral position in absence of external load for a wrist joint. The dashed black line stands for the targeted ellipse extracted from human biomechanical data in Pando et al. (2014). The blue line represents compliant behavior of the proposed system with an even distribution of the legs around the basis. (B) Compliance ellipses at the neutral position in absence of external load for an ankle joint. The dashed black line stands for the targeted ellipse extracted from human biomechanical data in Lee et al. (2014)). The blue line represents compliant behaviour of the proposed system with an even distribution of the legs around the basis. additional constraints. In this case, the two ellipses are well matched. It shows that the system could match additional passive behaviors, for instance by relaxing the constraint of passivity, previously explained (λ equal to 0). The associated design parameters of these 2 cases are given in Table 5.
To go further, we could also apply the proposed concept to more complex joints (such as the shoulder or the hip). Unfortunately, such biomechanical data are not yet available. We are more than aware that this is not an easy task and we can only encourage the biomechanical scientific community to pursue their current efforts in this direction, such as the work proposed in Hunt and Lee (2019).
As another further step, we could imagine having a lot of desired ellipses. The case 2.1 for instance corresponds to an isotropic FIGURE 10 | (A) (case 1.1) represents the compliance ellipses at the neutral position in absence of external load with λ null. (B) (case 1.2) represents the compliance ellipses at the neutral position in absence of external load without any additional constraints. (C) (case 2.1) represents the compliance ellipses in a large portion of the workspace. In all figures the dashed black lines stand for the targeted ellipses and the blue lines stand for the solutions. compliant behavior of the system in a large portion of the workspace. Such a system could be conceived for applications using teleoperated robots where we would like to have the same behavior over the workspace. Concretely, the required behavior is now a compliance circle over 25 positions over a grid of 5 × 5 positions between −60°and 60°in both directions. In this case, we have thus 75 characteristics to match. Similar results are obtained if the vector 2 λ is free and if λ is set to 0, ensuring thus the passivity of the compliant behavior. Therefore, only the results when λ is null, noted 0 25 , are displayed here. Figure 10C represents the result of the optimization algorithm and the parameters are given in Table 5. As it is shown in Figure 10C, the behavior is not exactly matched but it seems acceptable over a smaller range of the workspace.
Obviously, in the case of several desired ellipses, there are some mechanical limitations but we can refer to Section 4.2 to get the general trends and to design useful desired behavior. And it seems possible to design artificial joints that match the passive compliant behavior of human articulations in various positions of the workspace, to investigate the potential of such features in robotic systems.
Another limitation of this study is that we did not explore the space of all the possible spring characteristics. Indeed, we used a model (referred to as hyperbolic sine springs) that we already implemented in previous works (Catalano et al., 2011;Lemerle et al., 2019), as similar elastic behaviors can reproduce some attributes of an antagonistic pair of muscles driving a joint (Garabini et al., 2017). However, to improve the human-likeliness of the output stiffness characteristic, we could study more extensively on the best nonlinear spring characteristic that should be implemented. This was out of the intended scope of this work, but future studies could investigate further on this topic using additional biomechanical data to implement human-like stiffness behaviors.

CONCLUSION
This paper presents the concept of a new configurable 2 DoF variable stiffness joint. The system is based on redundant parallel nonlinear elastic actuation, using the antagonistic approach inspired by the human musculoskeletal system. The kinematic structure of the system is based on the architecture of the Omni-Wrist III, designed by Rosheim et al. and we propose to add in series nonlinear elastic systems to obtain a VSA. To get a compact design, only three motors are used and we prove that this allows us to control both the position and the overall stiffness of system. A one-dimensional quantity linked to the stiffness can thus be controlled like the voluntary co-contraction of human muscles. Moreover, we outline the impact of the main design parameters on the compliant behavior of the system as general guidelines. We describe an optimization methodology to fine-tune the mechanical implementation of this type of system to match the specific passive behavior of a human wrist and ankle.
We believe that the proposed design could find many applications, including industrial manipulators and prosthetic devices for various joints, due to its versatility. The sets of parameters obtained to match the passive behavior of a human wrist are particularly relevant as they open the possibility to design a compact system with an embedded passive behavior close to the human one. A new prototype of a compact physically compliant variable stiffness wrist is currently under development.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.