Design and Control of a Discrete Variable Stiffness Actuator With Instant Stiffness Switch for Safe Human-Robot Interaction

Variable Stiffness Actuators (VSA) have been proposed as an alternative actuation system for manipulators that are utilized for safe physical Human-Robot Interaction (pHRI). However, in the incidents of collision, the need of a fast response in stiffness tuning would rise to ensure safety. In this paper, we present a novel Discrete Variable Stiffness Actuator (DVSA) to be used in a compliant robotic manipulator for safe physical Human-Robot Interaction (pHRI). The novelty of this actuator lies in its design topology which allows the stiffness level to change swiftly among predefined levels without the need of complex stiffness tuning mechanism. Through this topology, three springs in parallel are connected serially between the motor and the link via gear train. The stiffness of the actuator is altered by adding/subtracting the number of involved springs, which can be realized through engagement/disengagement electromagnetic clutches on two of these spring’s shafts. The working principle, and the detailed design of the actuator are illustrated. Moreover, the stiffness model and the dynamic model are presented and discussed thoroughly. In order to validate these mathematical models and achieve optimal control, system identification for the dynamic parameters was performed experimentally on the physical model. Furthermore, the system’s ability of tracking desired trajectory was achieved through the implementation of different control techniques including PID (Proportional-Integral-Derivative), LQR (Linear Quadratic Regulator) and pole placement. The results show the high potential of utilizing the actuator in compliant manipulators. Moreover, DVSA is also characterized for safety in pHRI through Head-Injury Criterion (HIC). Finally, an application of DVSA in human augmentation task (Weight Bearing Task) is presented.


I. INTRODUCTION
A. RELATED WORK Over the past two decades, there has been a growing need for robotic manipulators working closely with human operators. Robot manipulators nowadays are associating, physically interacting or even worn by human operators. Although the general aspects of efficiency and robustness are taken into consideration as being of great importance, safety is The associate editor coordinating the review of this manuscript and approving it for publication was Pedro Neto . also considered to be highly ranked in priority with a growing concern among innovators. Safety should be an intrinsic feature in robots especially in the case of unexpected interactions, or sensor failures [1]. Besides the safety aspect, the interaction between the robot and the operator must show adaptability and force accuracy [2]. The criteria mentioned above have pushed research towards the development of variable impedance actuators (VIA), of which the actuator mechanical properties (inertia, damping, or stiffness) affect the system equilibrium position [1]. These changes will allow the compliance between the robots and the VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ environment or users to provide safe operations and energy efficiency [3]. Compliance in actuators can be classified into two major categories: active and passive. The active compliance is achieved through control by altering the behavior of a high-stiff actuator via software at a limited bandwidth [4], this is well-known as Variable Impedance Control (VIC). The stability of Variable Impedance Controller is subjected to several aspects. One aspect is the variation of the impedance parameters which would affect the stability of the system. An example of the proposed solutions for the stability of the VIC in that regard are proposed in [5], [6]. Furthermore, Variable Impedance Controller has been utilized on compliant actuators and manipulators in [7]- [10], where it shows significant capability in trajectory tracking. On the other hand, this method uses very accurate and expensive force/torque sensors, complex control systems, and has incapability of storing energy and absorbing shocks [11]. Furthermore, the sensors may fail and consequently may cause instability and safety issues during tasks in which the actuated load interacts with the surrounding environment or with humans.
As an alternative solution, passive compliance has been introduced into actuators to lower the control requirements and to improve the system safety in the cases of unexpected interactions or sensor failures. The passive compliance relies on mechanical elastic elements, placed between the internal motors and the actuated load. The internal mechanical compliance not only decouples the inertia of the motors from the load, thereby ensuring safety as an intrinsic feature especially during any kind of human-robot interaction, but can also be used to store energy, especially during tasks in which the kinetic energy can be absorbed during impacts and released when needed. Various methods have been applied to get passive (intrinsic) compliance, including the use of the series elastic actuators (SEA) [12], cable-driven joints [13], pneumatic [14] and hydraulic [15] actuation, etc. The trade-off between the control performance and safety leads the current research to focus and develop joints with variable compliance [16] to be able to cover safe interaction with low stiffness and better control performance with high stiffness. Designs that fall into this category include VSA-I [17], VSA-II [18], AMASC [19], smart MACCEPA [20], Variable Stiffness Actuator based on Gear-Rack Mechanism [21], the biological inspired joint stiffness control mechanism [22] and the beam-based variable stiffness actuator [23]. Moreover, another realization for stiffness altering is achieved through the principle of the lever mechanism. Examples of this type include the AwAS [24], AwAS-II [25], CompAct-VSA [1], the vsaUT [26], the mVSA-UT [27], the vsaUT-II [28] and pVSJ [2]. Although many different novel variable stiffness joints have been developed, many have not been applied successfully in many robot arms except [29] with their complex stiffness tuning mechanisms, bulky sizes, non-ideal stiffness curves, and the relatively slow response in stiffness tuning. The latter disadvantage would negatively impact the safety of such actuator in the presence of sudden changes in load, or incidents of collisions [30].
To overcome the aforementioned disadvantage, a new approach in varying the stiffness is proposed via discretely selecting the level of stiffness by adding/subtracting the number of involved elastic elements. The main advantage of discretizing the stiffness tuning compared to continuous stiffness tuning lies in (1) their capabilities of swiftly change the stiffness, (2) their energy efficiency due to the usage of generally low energy consuming clutching mechanisms and (3) the lower cost in control as stiffness tuning would rely on on/off control for the clutching mechanism [30]- [33]. In this stiffness tuning method, elastic elements can be arranged either in series or in parallel. Examples where elastic elements are arranged in series can be found in the pDVSJ [34] and pDVSJ-II [30]. On the other hand, the first realization where the elastic elements arranged in parallel can be found in the concept of the Series-Parallel Elastic Actuators (SPEA) [35]- [37]. In this paper, we propose a novel design (see, Fig.1) of the actuator inspired by SPEA approach. The additive advantage of this topology will allow more freedom in selecting the level of stiffness without the need of a successive involvement leading to lesser response time in altering the stiffness level. We believe that this novel topology can be applied in discrete variable stiffness actuators for safe and efficient compliant manipulators. Apart from the design issue, another challenging task is the control system design of the compliant actuators due to intrinsic compliance in their structure [38]. The successful implementation of the control involves the accurate modeling of the system and identifying the dynamic parameters correctly [39]. Any uncertainty in the dynamic parameters will eventually affect the performance of the control system. In literature, several algorithms for identifying dynamic parameters have been introduced to enhance the performance of the system [40]- [43].

B. CONTRIBUTIONS
In this work, we propose the Discrete Variable Stiffness Actuator (DVSA). The purpose of this actuator to be utilized in manipulators for safe physical Human-Robot Interaction tasks. The contributions of this manuscript can be listed as following: • Novelty: The novel design topology of DVSA, which allows the utilization of the fast locking mechanism which would bring instant online selection between significantly large different stiffness levels without conquering the back-driving forces in a continuous way with extra motors. This shows potential advantage of DVSAs for real-time human robot interaction changes. Moreover, this topology provides freedom in the shifting between the stiffness levels without the need of succession involvement of elastic elements.
• Modeling and Prototype Development: The physical prototype of DVSA is developed. The actuator's stiffness model and dynamic model are also developed and illustrated.
• System Parameters Identification: In order to achieve optimal control, the dynamic system parameters were identified using the parameter estimation tool in MATLAB which formulates optimization problem to perform the system identification for parameter estimation. The parameter estimation process involves two steps. In the first step, the dynamical equations of the system with an initial guess and boundaries of the unknown parameters are provided to the tool. In the second step, the actual response of the system is provided to the tool. The tool formulates the optimization problem to minimize the error between the simulated response and actual response of the system by estimating the values of unknown dynamic parameters. The identified dynamic parameters are used during the implementation of the control techniques to follow the desired output trajectory accurately.
• Control System Implementation:In this paper, we have implemented two types of controllers. The first implemented control technique is PID and the second is built on a state feedback tracking system designed with both pole placement and LQR method. For more details on these methods, the readers are encouraged to consult [44], [45]. The robustness of the implemented control techniques is evaluated by testing the system while switching between the different levels of stiffness in real time. In each case, we observed the tracking ability of the controllers to follow the desired trajectory.
• System Evaluation for Safety HRI: in order to utilize the proposed actuator in a HRI application, the actuator has been subjected to Head-Injury-Criterion, which is used to define the maximum operating velocities to ensure safety in HRI applications.
• Application in HRI: The proposed actuator has been utilized in a Human-Robot Cooperation Task (Weight Bearing task).
The rest of the paper is organized as follows: Section 2 presents the proposed DVSA consisting of its working principle, mechanical design, and prototype development.
The detailed model of the actuator comprising stiffness and dynamic modeling is reported in section 3. The dynamic parameter identification through the experimental setup is presented in section 4. In section 5, the control implementation and results are presented. Furthermore, the discussion of the experimental results is introduced to evaluate each controller in different aspects. In section 6, procedure and application in utilizing the proposed actuator in HRI applications are illustrated. Finally, in section 7, the conclusion and future work are outlined.

II. DVSA: DISCRETE VARIABLE STIFFNESS ACTUATOR
In this section, we present the discrete variable stiffness actuator (DVSA) by explaining its working principle, design, and prototype.

A. WORKING PRINCIPLE OF DVSA
The discrete variable stiffness actuator is designed to alter the stiffness by changing the number of involved elastic elements. Fig. 1 shows the working principle of the proposed actuator. A motor with a gear-head is connected to a ring gear which is engaged with three planet gears. The torsional springs are connected in parallel. An inline clutch is serially connected to a torsional spring and then another spur gear (planetary gear) which is then coupled to the main spur gear (ring gear) of the load shaft.
The springs are arranged in parallel to each other and their equivalent stiffness is connected in series between the motor (input) and load (output). One spring will be always connected to avoid total disengagement of the output link, whereas the other two which are attached to clutches from one side will be engaged to change the level of stiffness. When the clutch is engaged, the power transmission goes through the spring to the load. In the case, any clutch is disengaged, the corresponding spring will rotate freely without affecting the stiffness. The use of the clutches facilitates fast switching among the stiffness levels. Table 1 shows a binary representation of the stiffness levels based on the combination of an active/inactive springs. An active/inactive spring is represented by '1' and '0', respectively. The actuator's levels of stiffness can be selected according to the application requirement. Accordingly, the stiffness of each spring would be selected such as the combinations of stiffness for these springs would meet all of the levels as per to Table 1. The design of DVSA shares the topology concept of BpVSJ [31], however, there are two main differences between the two designs; the first difference is that VOLUME 9, 2021 BpVSJ utilizes a sun-planetary gear train, while DVSA utilizes a ring-planetary gear train which significantly affects the compactness of the design. The other difference is the addition of the actuation in DVSA as BpVSJ is a passive device.
Finally, DVSA is categorized under the Series-Parallel Elastic Actuators (SPEA). This category was proposed to minimize energy consumption, or power consumption in robots [33]. In literature, many concepts of SPEAs were proposed. The SPEA [35] and the MACCEPA-Based SPEA [36] are two examples of SPEAs which utilize an intermittent mechanism to recruit the parallel elastic elements in succession. The DVSA benefits from the gear train formation and the clutches to allow the freedom in selecting the level of stiffness without the need of succession involvement of elastic elements.

B. MECHANICAL DESIGN OF DVSA
The Computer-Aided Design (CAD) model of the actuator is reported in Fig. 2. The principle parts of the actuators are: a motor with a gearhead, two stages of ring-planetary gears sets, three springs, two clutches, and two types of ball bearings with different sizes. In particular, sixteen small bearings (single row deep groove ball bearings, OD=16 mm, ID=8 mm, Bore=4 mm) and two big bearings (single row deep groove ball bearings, OD=140 mm, ID=110 mm, bore = 16 mm). The gear ratio of the first ring-planetary gear set is selected to maintain the torque exerted on the clutch lower than the maximum allowable torque of the clutch. The second stage of the ring planetary gear was introduced to achieve again the maximum desired external torque. The two ring gears are (KHK-SI180: 80 teeth). While six planet gears are KHK-MSGA1-20: 20 teeth). These clutches were selected off-the-shelf with maximum torque to size ratio (HUCO SI-19). The dynamic parts are illustrated in Fig. 2B where the output shaft of the motor gearhead is connected with first planetary ring gear set through a ring gear which in turn, rotates the three planet gears set. Each planet gear is attached to a spring holder shaft. The spring holder shaft is firmly holding the one end of the spring while the other end of the spring is attached to the second spring holder shaft. Two of the second spring holder shafts are attached with two electromagnetic friction clutches while the third is directly attached. The output link of the actuator is attached to the ring gear of the second stage planet ring gear set.
The static structure consists of the connecting bars between both stages of the actuator, chassis of both stages, bearing carrier and motor chassis carrier. The bearing carriers encapsulate ball bearings to hold the planet gears and prevent any buckling due to the shear bending caused by the contact forces of the blocked gears.

C. PROTOTYPE DEVELOPMENT OF DVSA
The complete prototype of the actuator is shown in Fig. 3. The selected motor is the brushless flat EC-90 Maxon motor in combination with the reduction gearbox (91:1) along with integrated MILE encoder and Hall Sensor. The selected motor has exceptional features related to power-to-size ratio compactness, robustness, dust, and oil resistant and high precision. The parts of the actuator are fabricated through a computer numerical control (CNC) machine to increase the precision and strength. All the parts are manufactured from Aluminum or Steel to withstand the relatively high torques required of the actuator. Ball bearings are used at different parts to reduce friction and to avoid axis misalignment. The actuator is firmly fixed on a table with customized support.

III. DVSA: MODELING A. STIFFNESS MODELING
To derive the stiffness model, we can use the actuator's kinematics. Consider the schematic diagram shown in Fig. 4, if one end of the actuator is fixed and an input torque  is applied. The input torque will rotate the ring gear. In the scenario, when only one shaft is connected (both clutches are disengaged), the load torque will be transmitted through the corresponding spring (K l ) which deflects and produces a counter torque on the motor shaft. In equilibrium, the input torque (total torque) can be written as follows: where τ l , τ m , τ h are the torque transmitted by each planet gear (low, medium and high respectively) to ring gear. As the low stiffness spring is always connected without a clutch, the torque can be presented as: The other springs which are connected through clutches can be presented as follows: where N 2 is the gear ratio between the planetary gear and the ring gear. θis the angular position at the present time. φ is the angular position at the activation time. ϕ is the backlash angle. t ON ,n is time when the clutch (n) is activated. And K l , K m , K h are stiffness of the springs (low, medium and high respectively). The total torque (τ ) and the equivalent torsional stiffness (K eq ) can be obtained as: B. DYNAMIC MODELING In this subsection, we derive the dynamic equations of our system. The schematic representation of the complete prototype is shown in Fig. 5a. To derive the mathematical model of the system, we consider all the components shown in the orange box of Fig. 5a (except motor and its gearhead). Our approach is to simplify the model without compromising its functionality. Consequently, the effect of the three springs along with their three small gears are transformed to a central virtual axis as shown in Fig. 5b. Accordingly, the inertias of the small gears (planetary gear) are reflected on the virtual center axis is shown in (8). While the inertias of the ring gear and its holder are reflected to the virtual center axis in (9). Similarly, the equivalent inertia for the ring gear (J g2 ) reflecting on the virtual axis is shown in (10). In this work, PTC Creo 5.0 (CAD Software) is used to find the virtual inertia of the springs while the masses of the springs were measured through a weighing scale.
where: J eq1 : The equivalent inertia for the three springs reflecting into the central virtual axis. Accordingly, the system can be further simplified as shown in Fig. 5c. From that figure, the following equations can be driven: where J 1 : The total inertia before the equivalent spring location where all the inertias are reflected on the virtual axis. J 2 : The inertia of the load after the sun gear see Fig. 5b. θ in : The angular position of the input shaft before the sun gear see Fig. 5b.
θ in2 : The angular position of the input side after the sun gear but before the equivalent spring location see Fig. 5b.
θ L : The angular position of the load shaft after the sun gear see Fig. 5.
θ L2 : The angular position of the load side before the sun gear after the equivalent spring location see Fig. 5b.
Consequently, the model becomes as shown in Fig. 5c and the first differential equation is derived see (14).
Substitute (11) into (14) T in N 2 Equation (15) is the differential equation for the input side (motor side) after the model is simplified as shown in the previous steps. In addition, it is necessary to mention the viscous friction (C 1 ) is added to the model to express the friction between the mechanical components. Furthermore, the damping friction term for the motor is neglected due to its very small value of friction coefficient (0.0015) based on the motor catalog. Similarly, the same approach is used for the output side (load side): Equation (16) is the differential equation for the output side (load side) after the model is simplified as shown in the previous steps see Fig. 5c. Finally, to find the state space representation, let the state variables as it follows: Therefore, the state space representation for the system is shown in (17,18) where Based on this, the transfer function G (s) = θ L (s) T in (s) is derived to use it to design a position tracking system based on the PID controller as shown in section V. This is achieved by taking the Laplace transformation as shown below: where:

IV. SYSTEM IDENTIFICATION A. EXPERIMENTAL SETUP
The complete hardware setup and integration of different components of the system is shown in Fig. 6. The actuator is placed firmly in a horizontal orientation to curb the effect of gravity. We use a current mode for sending torque commands to EPOS2 motor controller using a DAQ system. For loading and unloading the springs in the actuator, we use the inertia dynamics clutches which are controlled using relays and the DAQ. The higher-level controller is implemented in MATLAB Simulink (on host PC) which had real-time data exchange and communication using TCP/IP with the target PC (low-level controller).
The control system of the DVSA is developed using MATLAB xPC target technique. The xPC target runs on an independent computer which allows real-time control. Executable code is loaded onto this target PC from a host PC running MATLAB Simulink software. In fact, the clock frequency of the CPU in the target PC that is used to implement the real time controller is 2.4 GHz. Furthermore, the control can be easily implemented with 1 kHz. The TCP/IP is used to establish the communication between the host and the target PCs (offline/online). The communication speed of TCP/IP is equal to 100 Mbps. The real-time data signals, acquired by target PC, are uploaded to the host PC over a direct crossover Ethernet cable (see, Fig. 7). The experimenter can monitor the data signal in the host scope and tune the parameters online in the Simulink model created in the host PC. The DAQ card (PCI 6221M, National Instruments) converts the digital command signal from the target PC into the analog signal and sends it to the EPOS 24/5 motor driver to control the motor. The command signals are generated by the control techniques presented in Section V. The DAQ card is also used to perform analog-digital conversion of the signals from the motor sensors (encoder/hall) and output encoder (E6B2-CWZ3E, 2048 counts per revolution). All the signals are sampled at 1 kHz. The clutches are engaged/disengaged through the 24V relays which are controlled by the target PC using the DAQ card. The supply voltage for the system was given by Siemens power supply 24V/40A.

B. EXPERIMENTAL RESULTS FOR THE ESTIMATED DYNAMIC PARAMETERS
The Simulink block diagram of the mathematical model is shown in Fig. 8. The model depends on the following parameters: C 1 , C 2 , J 1 , J 2 K eq and N 2 . But, N 2 is a gear ratio which is 0.25. Based on that, a Gray box model is used by using a parameter estimation tool in Simulink. It is worth mentioning that C obs is the output matrix for the observer due to θ in and θ L are measured by encoders (see Fig. 8). Which is presented as The parameter estimation tool works by exciting the model with an input signal e.g. chirp signal or a square signal. First, we simulated the dynamic model by using the approximated dynamic parameters from the CAD model as an initial guess for all parameters and later, we excited the model using the same input on the real prototype. In both cases, we recorded the responses elicited by the system. The parameter estimation tool is used to minimize the square of the error between the response of the real prototype and the response of the simulated model. This is achieved by providing the simulated model with an initial guess for the dynamic parameters and their expected boundaries. The estimation process is successful if the error reaches an approximation closest to zero and the corresponding estimated dynamic parameters are obtained. We tested the system under different input signals. Specifically, we initially used a chirp signal with a frequency in the range 0.01-3 Hz and with an amplitude ±35 N·m as shown in Fig. 9 and later, the square signal with a frequency 0.5 Hz and with an amplitude ±38 N·m (see, Fig. 10). Both types of signals are used to estimate the dynamic parameters. As a result, both simulated and measured responses are compared as shown in Fig. 9 and Fig. 10. The example of the estimated dynamic parameters in the case where there is a low stiffness level are listed in Table 2. Please note the significant difference in the values of C 1 and C 2 is due to the presence of a high gear ratio between the motor and input shaft. The same holds true for the inertia values i.e. J 1 and J 2 .

V. CONTROL IMPLEMENTATION AND RESULTS
In this section, we present the design and implementation of two types of controllers to control the system and have reported five case studies along with their results. The objective was to evaluate and compare the performance and robustness of each controller to track the reference trajectory under different conditions e.g. the variations of compliance level in  both cases offline and online (while the system is running). In particular, the implemented controllers are the position PID control and the extended tracking system. Moreover, the extended tracking system is implemented with both pole placement and LQR method.

A. POSITION PID CONTROL
The first implemented controller is a position PID controller to control the angular position of the load shaft θ L . The PID controller is the most common control algorithm used in many applications and various industries because its simplicity and a wide-range usage in operating conditions. In this paper the classical approach is used to design the PID controller.

• Case Study 1
The transfer function in (19) is used to design a position tracking system based on PID controller. The gains of the PID controller are obtained by root locus as a first guess then they are further tuned and adapted with the help of simulation results. The chosen gains for the PID controller were: K p = 100, K d = 12, and K i = 5. As a result, the tracking performance of the implemented controller and input torque are shown in Fig. 11a and Fig. 11b, respectively.

B. DESIGN OF EXTENDED TRACKING SYSTEM
The following equations are used to construct the control system of DVSA. The readers are encouraged to see more details on this method in [45]. Here, we recall the basic equations used to design the implemented controller. e =Âe + Bu e (20) where: where: r (t) = x 3 : The desired responses of θ L . K : The feedback gain matrix see Figure 12. k I : The integral gain see Figure 12. Furthermore, an observer was designed to estimate the unmeasured states. In our experimental setup, θ in and θ L i.e. x 1 and x 3 are the measured states.
In order to computeK , we have used two methods for the tracking system i.e.
(2) Tracking system based on LQR with integration and full state observer based on pole placement method.
The details of each method are following.

1) DESIGN OF A TRACKING SYSTEM BASED ON POLE PLACEMENT
In the pole placement method, the following command in MATLAB is used to find the extended gain matrix as shown in (24).K = place(Â,B, P) (24) • Case Study 2 The desired locations for the eigenvalues for the closed-loop control system are shown in (25) and the reference input is a step input which is then replaced with a sinewave. In this case study, there is no variations of compliance level and the stiffness level in this experiment is K eq = K L . Fig. 13, 14 show the results for the design of an extended tracking system based on the pole placement. In Fig.13a, the graph depicts the experimental results for the comparison between the desired (θ Lref ) and the actual (θ L ) results for the step input command. The motor input torque (T in ) is shown in Fig. 13b. It is noteworthy that there is a bit of latency (0.2 seconds, approximately) which can be the result of the system's backlash or other related hardware issues.
Furthermore, we also used the sine wave with a 0.05 Hz frequency as another type of input into our system in order to observe the robustness of the controller with a different input signal. Fig. 14a shows the comparative results of the desired and actual output (θ L ) of the system for the sinusoidal input while Fig. 14b shows the input torque of the motor.

2) DESIGN OF AN EXTENDED TRACKING SYSTEM BASED ON LQR WITH AN INTEGRATION AND A FULL STATE OBSERVER BASED ON THE POLE PLACEMENT
The second method aims to design an extended tracking system based on LQR. The pole placement method is used to design a full state observer to estimate the unmeasured states. The main task of the tracking system is to track the reference command (θ L ). This controller, in reality, is classified as an optimal controller because it minimizes a quadratic cost function as shown in (26) to calculate the optimal gains (K ) which leads to an optimal response by solving the reduced form for the Riccati equation [33], [34].
Subject toė =Âe +Bu e wherê Q: A semi-positive definite weighting matrix or positive definite weighting matrix.
where: K : Similar to (23). K : The feedback gain matrix see Fig. 15. k I : The integral gain see Fig. 15. S: The co-state matrix. This matrix must be a positive definite matrix in case there is a solution for the reduced form for the Riccati equation.
lamda: The closed loop eigenvalues for the extended tracking system. Fig. 15 shows the extended tracking system with a full state observer which is designed based on (20)- (23) and (26)- (27) to computeK . The closed loop eigenvalues for the extended tracking system based on LQR method is computed also by using the Matlab command in (27). The location of the eigenvalues of LQR are not controlled directly but they can be changed by a trial and error forQ andR matrices. Therefore this operation is executed to get a required performance for the system after that the eigenvalues can be computed as shown in (25).  [45]. k I : is the integral gain, K e : is the observer gains matrix, A, B, C are the matrices from the state space representation, C obs : are the observer C matrix which observes both input and out position states.
However, increasing theQ matrix forces the eigenvalues to move more in the direction of the left side on the s-plane i.e. it increases the relative stability in principle. On another hand, increasing theR matrix also decreases the relative stability to the system because it forces the eigenvalues to move more to the direction of the right side on the s-plane (all the eigenvalues are still on the left plane in the s-domain).

• Case Study 3
Similarly, we did not consider any stiffness variation i.e. the level of stiffness is kept fixed. The selected level of stiffness in this case is also K eq = K L . The values of thê Q andR matrices are obtained by a trial and error method to get an acceptable performance of the system. The closed loop eigenvalues for the extended tracking system are based on LQR are shown in (28). On the other hand, the selected eigenvalues for the closed-loop system for the observer (P o ) is shown in (29). The observer gains matrix (K e ) is computed by using the command in Matlab which is shown in (30).
A comparison between the desired and the actual angle (θ L ) is shown in Fig. 16a for the step input command. While the corresponding input torque T in is shown in Fig. 16b. The performance of the observer can be seen from the Fig's 17. Fig. 17a shows the error between the estimated and the measured outputs (θ in ,θ L ) while Fig's 17b-17c show the comparison between the measured state and the estimated state. As a result, all these figures validate the good performance of the implemented observer In this case, we also tested the system under sinewave with 0.05 HZ as shown in Fig. 18a. The corresponding torque is plotted in Fig. 18b. Fig. 18c reports the dynamic error between the estimated and the measured states to evaluate the performance of the observer which shows that the observer is working well.
In the following sections, we test the robustness of the implemented controllers using two case studies. The case FIGURE 17. Results of observe performance testing experiment, (a) Dynamic errors between the actual responses (y act ) and the estimated responses (y est ) for a step input. (b) A comparison between the measured and the actual θ in in rad/s for a step input. (c) A comparison between the estimated and the actual θ L in rad/s for a step input. studies are performed after selecting different levels of stiffness while the system is running and with a different reference input each time. This is achieved by performing a schedule for changing the stiffness level within the system at specific intervals which will be described in depth later. We tested all the levels of stiffness without changing the controllers' parameters.

• Case Study 4
In this case study, we consider the scenario of changing the stiffness online through the activation/deactivation of the clutches while moving in a step function as shown in Fig. 19a.
• During the first five seconds, both the medium and high stiffness level springs are disengaged, thus K eq = K l .
• During the following period [5-10] seconds, the medium stiffness level spring is engaged only, therefore K eq = K l + K m . • During the following period [10][11][12][13][14][15] seconds, the medium stiffness level spring is disengaged, while the high stiffness level spring is engaged, based on that K eq = K h + K l .
• Finally, in the last five seconds, both the medium and high stiffness level springs are engaged, therefore K eq = K h + K m + K l . Under these conditions, all the control techniques are applied to evaluate the performance of each in tracking the desired trajectory. The output (load) angle θ L was recorded for each case and the comparative results are shown in Fig. 19c. As a result, all the controllers were able to track the reference trajectory with reasonable accuracy. The PID and state feedback with pole placement have oscillation at the time of switching the stiffness level. While LQR controller resulted to be more robust even during the stiffness level switching. Moreover, LQR showed improved performance in terms of steady state error and oscillation. The schedule of the clutch engagement/disengagement (a) and corresponding stiffness level (b). The stiffness level is changed online (while the system is running) under the step input as reference. (Note, CM: Clutch associated with spring with middle level stiffness, CH: Clutch associated with spring with high level of stiffness) (c) A comparison between the desired and the actual (θ L ) for a step input to evaluate the performance of all the presented control schemes under the different level of stiffness obtained by changing the status of the clutches.

• Case Study 5
In this case, we test the system performance using the reference input as a sinewave having a frequency of 0.05 Hz. We tested the system by changing the level of stiffness through clutches while the system is running. The engagement or disengagement schedule of clutches is shown in Fig. 20a which can be summarized Fig. 20b as follows.
• During the first twenty-five seconds, there is only low level of stiffness i.e. both the medium and high stiffness level springs are disengaged, thus K eq = K l .
• During the following period  seconds, the low and medium stiffness level springs are engaged, therefore K eq = K l + K m .
• During the following period  seconds, the medium stiffness level spring is disengaged, while the low and high stiffness level springs are engaged, based on that K eq = K h + K l .
• Finally, in the last twenty-five seconds, all the springs are engaged i.e. both the medium and high stiffness level springs are engaged along with low level of stiffness spring, therefore K eq = K h + K m + K l .
Under these conditions, all the control techniques are applied to evaluate the performance of each in tracking the desired trajectory. The output load angle (θ L ) was recorded for each case and the comparative results are shown in Fig. 20c. As it can be seen that LQR showed consistent performance in terms of robustness as explained above even with the different reference input (i.e. sinewave).

C. RESULTS DISCUSSION
This section aims to discuss the presented results in the previous five case studies. The first three cases tested the tracking system which utilizes the PID controller, the pole placement, and the LQR respectively. These three cases have the same conditions e.g. the variation of compliance level is not taken into consideration. Moreover, all the controllers are designed based on the model when the stiffness level is low. The stiffness level is low in the first three cases as the medium and high springs are disengaged. It is important to mention a complicated, nonlinear, and dynamical friction such as backlash friction which is not taken into consideration because of the simplicity of the model where it may cause a steady state error to occur. Table 3 shows a comparison between the three controllers based on the settling time, the percentage overshoot, and the steady state error. In principle, all these controllers can be used in several applications. However, the authors tried to find the best controller from the vantage of performance, the robustness, and the energy saving. As stated in Table 3, the PID controller has the minimum settling time but it suffers from high values with the steady state error and the percentage overshoot in comparison to others. On the contrary, LQR has the minimum value for the steady state error with almost the same percentage overshoot in comparison to the pole placement. In conclusion, LQR is the best controller among them due to the performance as shown in Table 3 and the aspect of energy consumption shown later. This is achieved by observing the energy consumption for the LQR (see Fig. 16b) which is less than the energy consumption for the pole placement as shown in Fig. 13b. However, this point would be more interesting when it is observed and once the reference signal in this case, the pole placement, is equal to 0.5 rad (see Fig. 13a) while in case of LQR is equal to 2 rad as shown in Fig. 16a. This means the LQR consumes less energy even though it moves more. Finally, the fourth and the fifth case studies are designed to test the robustness of each controller determined by the variation of compliance level i.e. changing the stiffness level has led it to an increase in the uncertainty of the model especially, if we are aware of all the control parameters and know they have not changed at all. Furthermore, the control parameters are computed based on the dynamical model for the system when the equivalent stiffness is low (K eq = K l ).
It is observed in Fig's. [19][20], that the best controller is LQR according to its robustness during its performance which is almost unchanged. It is also worth noting, that the VOLUME 9, 2021 LQR has almost the same values for the overshoot and the steady state error even though the stiffness level has changed. The e performance for the pole placement is acceptable throughout all the interval periods except when the stiffness level has changed. During that time, the performance was poor because of the increase in the overshoot percentage. The PID controller had the worst performance and can be observed the whole time when changing the status of the clutches as shown in the figures.

VI. DVSA IN HUMAN-ROBOT INTERACTION APPLICATION
In this section, DVSA is characterized for safety in HRI. The Head-Injury Criterion (HIC) is used to determine the maximum operating velocities of a robotic link for a specific weight and joint stiffness [46]. Later, an example to illustrate the utilization of DVSA in human augmentation task (Weight Bearing Task) is presented.

A. HEAD INJURY CRITERION (HIC)
In the human robot cooperation field, a robotic system must not injure people during normal operations, nor during the incidents of operational error or mechanical failure. To quantitatively assess this risk, several standards have been developed. Head Injury Criterion (HIC) is one of the most commonly used criteria in many researches [47]- [49]. A HIC value of 100 m 5/2 s −2 is considered as the threshold in human robot cooperation tasks. It is safe for human to physically interact with the robot if HIC is less than 100 m 5/2 s −2 . For variable stiffness actuator (VSA), in [50], the formula for approximating the HIC based on the mass-spring-mass model was developed as a design criterion to minimize the risk of injury.
where K cov is the contact interface stiffness, M oper is the mass of head, v is the relative velocity between robot and head, M rob is the effective mass of the robot during the impact, and it can be defined as: where M link is the link mass, M rot is the rotor mass, K transm is the transmission stiffness which decouples the link mass M link and rotor mass M rot , γ is the rigid stiffness which is usually a large enough constant. There are 4 levels of transmission stiffness in our DVSA, which are 6.2335 N · m/rad, 29.0335 N · m/rad, 34.5335 N · m/rad and 57.3335 N · m/rad. The other parameters are: =3000 N · m/rad, M Oper = 4 kg, K cov = 25 kN/m, M rot = 1.35 kg, and using the link which mass is M link = 0.463 kg we can get the relation between the HIC and the relative velocity which is shown as Fig. 21.
Under the limitation of HIC=100 m 5/2 s −2 , the maximum velocity for minimum stiffness (level 1) is 2.2166 m/s; for level 2 stiffness, the maximum velocity is 2.1869 m/s; the maximum velocity for level 3 stiffness is 2.1800 m/s; for the highest stiffness (level 4), the maximum velocity is 2.1521 m/s. These maximum velocities yielded different stiffness levels and will be used in the future for human robot cooperation application.

B. HUMAN-DVSA COOPERATION IN WEIGHT BEARING TASK
In this section, we illustrate an example of one application of DVSA in the human-robot cooperation task: weight bearing. The experimental setup has been shown as Fig. 22. One link is installed to the output shaft of DVSA, and the load can be added at the end of this link. The distance between the output shaft and the handle is 300mm, and the weight of the whole link is 0.864kg. In this weight bearing experiment [51], the link is set horizontally with a (2.5kg) load at the end. The objective of the collaborative task is to maintain the horizontal stance throughout the subject and the DVSA. The subject will collaborate with the DVSA at its 4 levels of stiffness. The stiffness of the DVSA will be changed from the highest to the lowest. In the meantime, the muscle surface electromyography (EMG) [52] sensor (Delsys Trigno Avanti Platform) is used to record the activity of the bicep muscle. To estimate the metabolic power expended by muscles, the EMG raw data is processed with amplitude analysis.
Two subjects (Male: 24 years, healthy. And Female: 35 years, healthy) volunteered to participate in the experiment. The experiment consists of 5 time intervals, with 5 seconds for each interval. In the first interval (t∈[0,5]s), the stiffness of the DVSA was set to its maximum, the EMG signals of both the subjects were recorded and the average of their muscle activity is plotted. It is observable in Fig. 23 that the subjects have not utilized much of their effort as the DVSA was holding most of the load during the first interval (t∈[0,5]s). During the second time interval (t∈ [5,10]s), the stiffness of DVSA was lowered to the third level, and the EMG sensor recorded a slight change in the muscle activity at (t∈ [5,6]s), however, it maintained its minimum level indicating that the DVSA is still holding most of the weight. During the third time interval (t∈ [10,15]s), the stiffness level was lowered to the second level. The EMG sensor recorded a significant change in muscle activity than the first two intervals, which then faded again indicating that the load had been carried by the motor. During the fourth interval (t∈ [15,20]s), the lowest stiffness level had been set. The EMG sensors recorded significant muscle activity within the first 2 seconds where it faded at a slower pace compared to the previous intervals. This indicates that at the lowest stiffness, more muscle activity is needed to stabilize the load. During the last time interval, the DVSA had been switched off. The EMG recorded continuous muscle activity, indicating that the subjects were only utilizing their muscles to keep holding the weight. From the previous experiment, it can be seen that the DVSA can be utilized in Human-Robot Collaborative tasks where such tasks can be extended into the realm of rehabilitation and sports training purposes.

VII. CONCLUSION AND FUTURE WORK
In this paper, a novel discrete variable stiffness actuator is presented. This actuator will be used in a compliant manipulator that will be applied for safe Human-Robot-Interaction (HRI).
The working principle of this actuator falls under the category of Series-Parallel-Elastic-Actuators. However, it franchises itself through its novel topology of changing the stiffness level. Firstly, the detailed design and development of the actuator are presented. Later, the mathematical model consisting of stiffness and dynamic model is detailed. Moreover, the model parameters were estimated through parameter estimation tool in SIMULINK, and they are validated throughout the experiments on the real system. Furthermore, we implemented different control techniques on the prototype of the actuator to evaluate their ability to track the desired reference trajectory under different conditions associated with the change in stiffness level and different reference signals. In particular, the position tracking system based on PID (Proportional-Integral-Derivative) and extended tracking system are implemented. The extended tracking system is implemented both with the pole placement and with the Linear Quadratic Regulator. The results show that the actuator has a high potential to be used as compliant joint in manipulators. Moreover, we also characterized the DVSA for safety in HRI using Head-Injury Criterion (HIC) and an application of the DVSA in human augmentation tasks (Weight Bearing Task) is presented. The current design is significantly more compact when compared to its counterpart presented in [31]. However, the bulkiness of the current version is referred to the fact that all the used parts were off-the-shelf. Therefore, the room for improvement to achieve more compactness in future can be summarized by designing and developing customized parts like stronger and more compact clutches. Moreover, utilizing newly developed flat springs may help in realizing a more compact version of DVSA. Furthermore, we are aiming to develop the complete compliant manipulator using multiple joints proposed by the DVSA in the future.