Design, modeling, and control of a variable stiffness elbow joint

New technological advances are changing the way robotics are designed for safe and dependable physical human–robot interaction and human-like prosthesis. Outstanding examples are the adoption of soft covers, compliant transmission elements, and motion control laws that allow compliant behavior in the event of collisions while preserving accuracy and performance during motion in free space. In this scenario, there is growing interest in variable stiffness actuators (VSAs). Herein, we present a new design of an anthropomorphic elbow VSA based on an architecture we developed previously. A robust dynamic feedback linearization algorithm is used to achieve simultaneous control of the output link position and stiffness. This actuation system makes use of two compliant transmission elements, characterized by a nonlinear relation between deflection and applied torque. Static feedback control algorithms have been proposed in literature considering purely elastic transmission; however, viscoelasticity is often observed in practice. This phenomenon may harm the performance of static feedback linearization algorithms, particularly in the case of trajectory tracking. To overcome this limitation, we propose a dynamic feedback linearization algorithm that explicitly considers the viscoelasticity of the transmission elements, and validate it through simulations and experimental studies. The results are compared with the static feedback case to showcase the improvement in trajectory tracking, even in the case of parameter uncertainty.


Introduction
Transhumeral amputations (i.e., between the elbow and shoulder) are the second most common type of upper limb amputation, accounting for approximately 16% of all upper limb amputation surgeries [1]. This highlights the importance of developing prostheses concerning not only the hand and wrist, but also the elbow. Most daily tasks require the use of the upper limbs; therefore, the loss of one or both upper limbs can complicate a variety of tasks, which can significantly worsen the quality of life.
If we take the human limb as reference to create a prosthesis as natural as possible, we cannot avoid introducing the concept of variable stiffness, as suggested by Sensinger et al. [2]. In fact, the muscle and bone structure of the human arm allow for stiffness modulation by contracting a pair of agonistic and antagonist muscles. In particular, the alternate muscle contraction leads to the arm motion, whereas the simultaneous contraction leads to stiffness variation [3]. The artificial counterpart of such behavior may be found in variable stiffness actuators (VSAs), which are now extensively studied and employed within the robotic/mechatronic community. In the last decade, several examples of VSAs for robotic manipulators with programmable stiffness have been designed and developed by different research institutes, with the aim of verifying the effectiveness and reliability of the variable stiffness approach (see, e.g., [4][5][6][7]). In particular, Vanderborght et al. [8] comprehensively analyzed VSA uses and possible configurations, starting with an initial distinction between active and passive VSAs. In active VSAs, a motor flanked by a dedicated control system performs the stiffness change based on the signals returned by sensors at the joint. This type of VSA is able to undergo theoretically infinite stiffness changes, and is widely used for industrial [9] and biomedical applications, including in the construction of prostheses [2,10,11].
On the other hand, passive VSAs, also named structurecontrolled systems, generally rely on the change in geometry of some passive spring whose preload is used to achieve stiffness variations [12][13][14]. In terms of prime movers, either traditional electric motors (see, e.g., [15,16]) or intrinsically soft actuators may be employed (e.g., McKibben acutuators [17] or dielectric elastomer actuators [18]). In all cases, owing to the presence of one or more inherently compliant or soft transmission, hereafter referred to as compliant transmission element (CTE), passive VSAs present several advantages over their active counterparts. For example, they protect the joint from unexpected impacts by decoupling the joint and motor through the CTE deformation [19], which makes them safer in the event of collisions with things or people. In addition, they maintain the inherent stiffness of the joint even when the motor is deactivated. Finally, since the CTE may be used as an elastic energy buffer to be released as needed, its inclusion may allow to downsize the motors for a given peak-torque requirement.
Owing to the above advantages, passive VSAs are widely used in collaborative robots (i.e., robots designed to work alongside humans) [20][21][22], exoskeletons for the rehabilitation of upper and lower limbs [23][24][25][26], and prosthetic devices [27]. In the latter reference, Lemerle et al. reported the first example of a variable stiffness elbow joint using an antagonistic motor configuration, which is a common architecture discussed in the literature [14,16,28]. However, major drawbacks of all VSAs, despite being active or passive or employing traditional or soft prime movers, are (1) a rather complex mechanical design due to the need of a couple of prime movers along with additional compliant components to simultaneously adjust both position and stiffness of the joint; and (2) a more complicated control approach, as compared to traditional, single-motor, servosystems, especially when dealing with parameter uncertainties or undesired system behaviors to be compensated (such as, for instance, non-negligible damping [27][28][29]). Within this scenario, it is self-evident that any possible solution aiming at system simplification, either at the mechanical design or control level, is welcome and shall be pursued.
Focusing at first on the mechanical design and, in particular on the CTE implantation, a solution possibly leading to part-count reduction may be based on the compliant mechanism concept [30]. In fact, the growing interest on compliant mechanism with localized [31] or distributed [32] compliance is mainly due to the possibility of making simpler and more compact mechanical systems with completely customized characteristics (such as non-linear springs [33,34], deformable structures employed as rotational/prismatic kinematic pairs [35][36][37], and spatial single-piece mechanisms [38,39]). Naturally, by relying on the deflection of flexible members for their functioning, the performances of a generic compliant mechanism are heavily dependent on the binomial material-morphology. In particular, at present, complex morphologies may be realized with linear elastic materials (e.g., spring steel) or even with low-cost plastics for additive manufacturing (e.g., polylactic acid -PLA or thermoplastic polyurethane -TPU). In the second instance, it is often the case that the compliant element behaves as expected in quasi-static conditions, though it is highly affected by undesired time-dependent phenomena (e.g., damping) during motion. Therefore, in the hypothesis to resort to CTE morphologically conceived as distributed compliant mechanisms but realized with materials affected by viscoelasicity, the mentioned design philosophy aiming at a simplification of the VSA mechanical structure is achieved at the expense of an undesired behavior in dynamic conditions, thus limiting the joint precision during actual tasks. However, such unwanted behavior can be effectively mitigated/compensated by means of appropriate control strategies. In this context, the technical challenge is to design a variable stiffness device comprising intrinsic compliance while compensating for the disadvantages of viscoelasticity.
In light of these considerations, the aim of this study is to develop a novel anthropomorphic variable stiffness elbow joint that can be used as either a human limb prosthesis or a human-like robotic arm, in combination with a robust control strategy to reduce positional error due to CTE hysteresis. This device, shown in Fig. 1, uses two nonlinear CTEs in an antagonistic configuration to achieve a stiffness variation of about 13,500 Nmm/rad. The CTEs, each consisting of an inner and outer ring connected by four slender beams, were designed using a Matlab/ANSYS APDL optimization routine and then prototyped and tested. The experimental results showed good consistency with the finite element analysis (FEA) data, albeit with high hysteresis owing to the internal damping of the material. To overcome this drawback, a mathematical model of the VSA was prepared to develop a robust control strategy. Finally, the joint was fabricated and tested for validation of the model. The rest of the paper is organized as follows. Section 2 reports the static model of the antagonistic VSA; Sect. 3 describes the design, main features, and capabilities of the presented device; Sect. 4 is dedicated to the design and experimental testing of the CTEs that provide the stiffness variation; Sect. 5 discusses the device model and control strategy; and Sect. 6 presents the experimental results of the physical prototype. Finally, Sect. 7 reports the conclusions of the work.

VSA static model
In variable stiffness antagonistic devices [14,15,40], as schematized in Fig. 2, two actuators are employed for each moving joint, which contribute to the position and stiffness of the joint by acting on the two CTEs that are connected between the motors and moving joint. The torques applied on the inner rings of CTEs by the motors, and , are given by (1) = ( ) + J m̈w here = − q and = + q are the angular deflection of the CTEs; J m is the inertia of the motor; = ( ) = (q, ) and = ( ) = (q, ) are the torques applied by the motors to the moving link through the deflection of the CTEs; and are the positions of the motors; and q is the angular position of the moving link. For the static case, ̈ and ̈ both equal zero; therefore, the total torque applied to the joint j is equal to To independently control the position q and stiffness k = j q of the moving joint by guiding the angular positions of the motors ( and ), the torques and must be nonlinear functions of the angular deformations and .
Assuming a quadratic profile, we can write the torque of the motors and joint as follows: where a 1 and a 2 are the first-and second-order coefficients, respectively, of the quadratic torque-deflection profile. The stiffness profile k will then be equal to Equation 6 shows that nonlinear CTEs are essential for achieving stiffness modulation, because linear CTEs would have a zero a 2 coefficient and thus a constant joint stiffness of k = 2a 1 .

VSA design overview
This section describes the design of the variable stiffness elbow. To be used as an upper limb prosthesis, this device must be adaptable to replicate the shape and size of the missing limb. The design of the arm presented in this article derives directly from our previous work [41]. Our first prototype of a passive antagonistic variable stiffness joint, shown in Fig. 3, used a pair of acrylonitrile butadiene styrene (ABS) CTEs and achieved a stiffness variation range of 726-1795 N mm/rad. To transition from a variable stiffness joint to a humanoid elbow joint, two major design challenges must be overcome. The first concerns the geometric constraints of the device owing to the need to replicate the human limb as faithfully as possible. The second concerns the physical capabilities of the joint itself and thus the range of motion and stiffness that can be achieved. According to NASA [42], the human upper arm  The main shaft, motors, IMU sensor, encoder, bearings, motor, and shaft gears were standard parts. All other components were made by additive manufacturing. Specifically, the encoder and forearm transmission gears were fabricated from a thermoset resin (Formlabs Tough 1500) using a Formlabs Form 3 printer; the two CTEs were fabricated from thermoplastic polyurethane (TPU) and all other components were made from polylactic acid (PLA) using an Ultimaker S3 printer.
The final design of the CTE is shown in Fig. 5. It consists of four slender beams with distributed compliance that are connected at either end to the inner and outer rigid rings of the CTE. In this way, when the inner ring rotates, torque is transmitted to the outer ring. Because it would be very difficult to obtain an ideally rigid outer ring in the physical prototype, especially given the dimensions of the beams compared to  Table 1.

CTE design and testing
This section presents the design and characterization of the CTEs used to adjust the stiffness of the device. Referring to Fig. 6, the neutral axis of the beam was determined using a series of six points interpolated over a cubic spline.
The initial and final points ( p in and p fin , respectively) were fixed at the inner and outer rings, respectively, and the positions of the remaining four points ( p 1-4 ) were determined using a Matlab/ANSYS APDL optimization routine, whereby ANSYS was tasked to solve the 1D model of the beam while Matlab managed the optimization procedure in batch.
The optimization of the beam profile to obtain the desired nonlinear behavior is discussed in more detail in our previous article [41]. An overview of the procedures and tools for designing beam-based compliant mechanisms can also be found in ref. [43]. For the sake of brevity, only the optimal values of this particular profile are given here (see Table 2), and only the dimensions of the beam section were changed to tune the torque performance of the device.
Given the nature of the device, we aimed to maximize the torque response of the planar spring while respecting the maximum stress of the material and the maximum dimensions of the robotic joint.
Based on the material used in the previous article (i.e., ABS, Young's modulus E = 1800 MPa, Poisson's ratio = 0.35), an optimization procedure was carried out in which the width B of the beam section (see Fig. 6) was fixed at the maximum value allowed by the overall dimensions (10 mm) and the thickness H was varied cyclically to find the configuration with the maximum response, with the only constraint being that the maximum allowable stress of the material ( ABS, max = 42.5 MPa) was not exceeded. The maximum torque determined in this initial optimization  Fig. 7). Notably, the maximum torque values are about three times greater than those using ABS, for which a thickness of 29 mm would have been required (compared to the threshold thickness of 10 mm allowed based on the dimensions of the elbow) to obtain the same torque as that obtained with TPU. The optimal CTE was then verified using a 3D FEA simulation, as shown in Fig. 8.

Experimental evaluation
A prototype of the optimized TPU CTE was fabricated using a fused deposition modeling (FDM) Ultimaker S3 printer and tested with a customized experimental setup comprising a Kollmorgen DC motor, six-axis ATI load cell, and pair of connecting ABS flanges (see Fig. 9).
The inner ring of the CTE was connected to the motor flange (black component), while the outer ring was connected to the load cell flange (red component). To control the motor drive and read the torque values captured by the load cell, a real-time LabVIEW script was used with an embedded NI-cRIO controller. A sinusoidal motion law was then imposed on the motor ( 15 • amplitude, 1 Hz frequency) and the corresponding torque/rotation law of the CTE was obtained.
As shown in Fig. 10, the experimental results were in agreement with those of the FEA analyses, minus possible errors due to the 3D printing process. The torque-deflection relationship curve of the TPU CTE, as mentioned in Sect. 3, has a similar shape to that obtained for the ABS CTE [41]. It is important to note that this curve must be asymmetric to avoid zero stiffness at the minimum point and thus a singularity point in the control. Instead, the symmetry of the system results from the fact that the two CTEs are mounted in an antagonistic configuration.
Dynamic characterization was performed to visualize the effects of the internal damping of the material. There was an evident hysteresis effect in the loading and unloading phases of the spring, especially when compared to the ABS CTE used in the previous work [41]. The damping coefficients of the ABS and TPU CTEs were estimated using Recurdyn software, which is able to reproduce the dynamics of the CTEs under the same conditions as those used in the experimental tests. Damping coefficients of 0.0025 and 0.04 N s/mm were achieved for the ABS and TPU CTEs, respectively. This corresponds to a 16-fold difference; therefore, the effect of damping cannot be neglected for the TPU CTE, as it will have a significant effect on the control of the device (Fig. 11). For this reason, in the following sections, the model of the device is developed with consideration to damping, and a control strategy is proposed to compensate for its effect.

Modeling and control
Referring to Fig. 2, the robot dynamic model is composed by the dynamics of 3 rigid bodies (the link and 2 actuators), coupled through elastic transmission elements. The dynamic model of the system is as follows: whole manipulator can be written by grouping the actuators in two sets, denominated here as and , described by two equations, one for each actuators set: where q is joint positions, J q is the link inertia matrix, b q is the joint damping coefficient, and are the actuator positions, and J m is the diagonal inertia matrix of the actuators, while and are the torques applied by the actuators and e is the load torque applied to the link (e.g., the gravity torque acting on the output link).
In case of quadratic transmission elements, such as in [14,40], the coupling torque between the actuators and the moving link can be expressed as where is the the spring displacement = − q and = + q for the two actuators, respectively, a 2 and a 1 are the quadratic spring coefficients, and b 1 is the transmission elements' damping factor.

Feedback linearization
Since we are interested to control both the position and the stiffness of the joint, it is useful to define a new output vector that contains directly this information: where k is the mechanical stiffness of the joint and k 0 = 2 a 1 is the stiffness value in the equilibrium condition. By assuming = ( − )∕2 and s = + , it is possible to write where = − and k = + , k(s) > 0 is a strictly positive function representing the generalized joint stiffness, (q − ) is a odd strictly monotonically increasing functions representing the generalized joint displacement, while (q − , s) is a function such that (0, 0) = 0.
In case of transmission elements characterized by the torque displacement relation (10), the following relations hold: By dropping the dependance of k from s and inverting its relation, we obtain It can be also noted that k = 2 a 2s .
Since the joint stiffness k is a function of s, that is in turn a function of the actuators configuration, k is in general a function of time, i.e., k = k(t) . Moreover, it is k(t) > 0 for all t, since it has no physical meaning to consider negative stiffness while if the stiffness drops to zero the joint/transmission would lead to an unactuated system.
Taking into account for model uncertainties, the system dynamics can be rewritten as Moreover, it is assumed that all the effects due to frictions, dead-zones, non-modeled dynamics, parameters variability, etc., not considered in the other terms of the dynamics, can be collected by introducing the additive functions of time Therefore, the controllable input is u = [ k ] T , the non-controllable input is w = e , and the robot state is x = [q̇q̇k̇k] T . Considering the new output vector y c , the input u, and the state vector x, the dynamic model of the VSA joint in the nominal case (i.e., without model uncertainties) can be rewritten in the state space form: where x ∈ ℝ 6 , u, y c ∈ ℝ 2 and w ∈ ℝ . In particular, for the model of the antagonistic actuated joint: Considering the control input u and the output information y c , Eqs. (11)-(15) define a square nonlinear system that can be linearized via output feedback if suitable conditions are satisfied [40,44,45].
Starting from the analysis on the feedback linearization of antagonistic actuated robots carried out in the previous chapter, we have to check, first of all, some conditions on the relative degrees of the output information. In this preliminary analysis, only the controllable input u of the system is taken into account while the external torque e applied to the joint is neglected. In this case, for the output q it is possible to write: while for the output S: where L i f h q (x) stays for the i-th Lie derivative of h(x) with respect to the generic function f(x) and L g and L g k denote the restriction of the Lie derivative respectively to the and k component of the input vector while L g ( ,k) denotes both the cases above mentioned.
If the disturbance input, the external torque applied to the joint, is considered, the relative degrees of the output can be defined as: while for the output k: On the basis of this consideration, it is possible to state that the joint position cannot be decoupled from the external torque, while the joint stiffness is not affected by the external torque, since the relative degree of the stiffness w.r.t the control input is lower that w.r.t. the disturbance. The external torque e will be neglected in the sequel to simplify the analysis. The problem of feedback linearization under the effects of an external torque will be an object of future research.

Static feedback linearization
From Eqs. (16)-(21), we can state that, considering the controllable input u, the vector relative degree of q is 3 while the one of k is 2. The sum of the vector relative degrees of the output is then not equal to the dimension of the state of the system, so the conditions for the solution of the full feedback linearization problem are not satisfied [44]. Anyway, the system can be partially linearized via static feedback if it is possible to define a nonsingular coordinates transformation from the space of the original state variables to the state space of the partially linearized one [46]. To solve this problem, since the additional entries of the coordinates transformation (the ones that must be added to make the coordinates transformation nonsingular) can be arbitrary chosen, it is convenient, if possible, to make the nonlinear part of the new system independent from the inputs.
The state of the partially linearized system z can be split into the linear part: and its nonlinear part: where z 6 will be chosen in suitable way as shown in the sequel. Hence, the overall system state in the new coordinate system is: The dynamics of the partially linearized system can be then written as: in which Q(x) is the so-called decoupling matrix and the matrices A, B, and C are in the Brunowsky canonical form where Eq. (27) represents the nonlinear part of the system and z = Φ(x) represents the coordinates transformation between the original and the partially linearized system. Obviously the term L 2 f h q (x) =q can be expressed in terms of the state vector x by using the second element of Eq. (13), while Φ 6 (x) must be chosen in such a way to make the coordinate transformation Φ(x) nonsingular.
It is important to note that the nonlinear residual part of the system (z l , z n , u l ) is unobservable from the outputs, as can be seen looking at Eq. (28). In particular, the stability of the zero dynamics is a necessary condition for the controllability with bounded internal state of the overall system [44]. Looking at the first five rows of dΦ(x) dx : (26) z l = Az l + Bu l (27) z n = (z l , z n , u l ) it is clear that the third and the fourth columns are linearly dependent. Therefore, Φ 6 (x) must be choosen to make this matrix full rank. By choosing z 6 =ṫ he state space transformation Φ(x) is nonsingular in the whole state space, and in particular in the origin (since k is always positive by assumption). The state space transformation Φ(x) is then: The state space transformation form the original to the linearized system can be then expressed only by means of the state space information of the original system. The dynamic equations that describe the nonlinear part of the system is then: Even if the nonlinear residual part of the system depends on the input, it is important to note that the zero dynamics is asymptotically stable in the whole state space, so the system can be controlled with bounded internal state.
To define the input of the linearized system, see Eq. (29), the nonlinear terms L 3 f h q (x) and L 2 f h k (x) and the decoupling matrix Q(x) must be defined first: By considering Eqs. (17)- (21), the decoupling matrix Q(x) is defined as that is always nonsingular. The control law that allows to obtain the linearized system is then: where v = v q v k T is the new input of the system. Note that this control law is static, since it can be be computed directly from the system state information x and from the new reference input v. The dynamics of the system Eqs. (26)- (29) is, in the z coordinates: By recalling the definition of the matrices A, B, and C, it is possible to write: An important point to note is that if the damping coefficient of the transmission elements b tends to zero (considering also the case in which b can be neglected), the control law Eq. (34) becomes ill conditioned and the zero dynamics Eq. (32) is not exponentially stable anymore. Due to the structure of the matrices A, B, and C, the system described by Eqs. (35) and (37) is completely controllable and observable, while Eq. (36) represents its unobservable part. This implies that the state z l can be estimated by means of an asymptotic state observer and of the change of coordinates Φ(x) or, since the position of each rigid body is directly measurable, the velocities can be estimated in many ways, e.g., by means of state variable filters or adaptive windowing algorithms [47].
For the control of the partially linearized system, an outer control loop can be used to solve the problem of simultaneous decoupled stiffness and position trajectory tracking. From Eq. (38), applying the control laws: where K 2 q , K 1 q , K 0 q , K 1 k , and K 0 k are such that are Hurwitz polynomials, the convergence to zero of the tracking error is ensured. If the desired trajectories q d is continuous together with its derivatives up to 3th and k d is continuous together with its fist derivative, also the asymptotic trajectory tracking is achieved. The control law in Eqs. (39) and (40) is equivalent to a static state feedback plus feedforward in the state space of the linearized system: and Φ (q,k) (x) stays for the restriction of the coordinate transformation to the state component z l : The matrix K can be also obtained via direct eigenvalues assignment or through the solution of the CARE equation with an opportune choice of the weight matrices.

Dynamic feedback linearization
Since the sum of the relative degrees of the inputs is not equal to the dimension of the state space of the system, to achieve full state linearization and to solve the problem of a vanishing damping factor b , a dynamic feedback compensator must be designed.
By choosing, for the linearized system, the new state vector: the dynamics of the linearized system can be then written as: [3] k − k 0k T (48) z = Az + Bu l where while the matrices A, B and C change in: The expression of L 4 f h q (x) is not reported for brevity, since it can be derived by computing the Lie derivative of L 3 f h q (x) with respect to f(x), but it is important to note that, also in this case, it is an algebraic function of the system state x. Moreover, note that if the coefficient b vanishes in Eq. (50), it is possible to show that the case of static feedback linearization with no damping in the transmission elements is recovered. The coordinates transformation between the original and the linearized system is then: Now, by defining the dynamic compensation law: we obtain the fully linearized system that, together with Eq. (54), give the whole dynamics of VSA joint plus the dynamic controller. By recalling the definition of the matrices A, B, and C in Eq. (51), it is possible to write:

Elbow joint physical prototyping
A physical prototype of the device was produced by additive manufacturing and tested to verify its functionality. Figure 12 shows the actual change in stiffness of the joint with a 1-kg weight attached to the end of the forearm section. Stiffness values of 15 830 and 2417 N mm/rad were achieved in the maximum and minimum stiffness configurations, respectively.
To validate the control algorithm presented in Sect. 5, the physical prototype used in the previous work [41] (see Fig. 3) was subjected to further experimental tests to compare the results. In particular, the errors related to the position and stiffness of the joint were compared in the following cases: (a) static feedback control without damping; (b) static feedback control with damping; and (c) dynamic feedback control with damping. All three cases were considered by providing the joint with both a sinusoidal and stepwise position setpoint (Fig. 13) and a stepwise stiffness adjustment (Fig. 14) evaluated both individually (Fig. 15) and during sine motion (Fig. 16). Figures 17 and 18 reveal that the positional error was considerably reduced when moving from case a (static feedback without damping) to cases b and c (static and dynamic feedback with damping, respectively), which take the damping of the material into consideration. This demonstrates the need for a dynamic feedback linearization algorithm that explicitly considers the viscoelasticity of the CTEs.

Conclusion
This article reports the design of a variable stiffness elbow joint using a pair of motors and two nonlinear CTEs with an antagonistic configuration. The design of the device, which features a reduced number of components and is suitable for 3D printing prototyping, matches the appearance of a human upper limb so that it can be used as a prosthesis or as a component of a humanoid robot. The CTEs, consisting of two rigid rings connected by four beams with distributed compliance, were designed to have nonlinear characteristics and provide maximum variation in stiffness while remaining within the defined maximum dimensions of the device. Subsequent FEA and experimental tests confirmed the design values. Non-negligible hysteresis was present, which was counteracted by a robust control strategy that was experimentally verified using a simplified setup. Future work will extend the experimental testing of the final device both in the laboratory and in situations of normal daily use. In addition, a control system using electromyography signals will be implemented for use by amputees.