Flexible Control Strategy for Upper-Limb Rehabilitation Exoskeleton Based on Virtual Spring Damper Hypothesis

: The focus of this work is to design a control strategy with the dynamic characteristics of spring damping to realize the virtual ﬂexibility and softness of a rigid-joint exoskeleton without installing real, physical elastic devices. The basic idea of a “virtual softening control strategy” for a single rigid joint is that a virtual spring damper (VSD) is installed between the motor and the output shaft. By designing the control signal of the motor, the torque output of the joint actuator is softened so that the output has the characteristics of elasticity and variable stiffness. The transfer velocity proﬁle of human limbs reaching from one posture to another always presents as bell-shaped. According to this characteristic, we constructed a trajectory planning method for a point-to-point position-tracking controller based on a normal distribution function, and it was successfully applied to the control of 5-DoF upper-limb rehabilitation exoskeleton. A multi-joint cooperative ﬂexible controller based on the virtual spring damper hypothesis (VSDH) was successfully applied to solve the constrained control problem of the exoskeletons and the self-motion problem caused by redundant degrees of freedom (DoFs). The stability of the closed-loop controlled system is theoretically proven by use of the scalar energy function gradient method and the Riemann metric convergence analysis method.

Existing rehabilitation exoskeletons are mainly driven by rigid motors, which can achieve high control accuracy but low compliance [13]. In the case of excessive movement or sudden abnormal force, it is easy to cause secondary injury to the patient's affected limb [14]. Current exoskeleton robots limit the movement comfort and flexibility of the wearer and have safety problems, due to high-rigidity actuators and low material impact resistance [15]. The movement of rigid exoskeletons presents a mechanical stereotype; the exercise process is uncomfortable, and the effect of the rehabilitation treatment is also affected.
Many practical engineers and academic researchers have shown great interest in studying the application of flexible joints in exoskeleton manipulators [16]. In recent years, research on flexible exoskeleton robots has accelerated, from the laboratory to real life. Flexible exoskeleton robots are characterized by high flexibility, high environmental adaptability, and high compatibility [17].
Flexible joints have been utilized in different kinds of robotic applications, including rehabilitation robots [18,19], walking robots, humanoids [20], and haptics. According to whether the joint involves stiffness regulation, a flexible joint can be divided into two categories: series elastic actuator (SEA) and variable stiffness actuator (VSA) [18]. Regardless of which kind of physical flexible actuator is applied to the upper-limb exoskeleton, it will bring about three problems: First, it will increase the volume, weight, structural complexity, and manufacturing cost of the device. Second, the control of position accuracy will become extremely complex and difficult. Third, it will bring unwanted oscillation.
The main objective of upper-limb rehabilitation exoskeletons is to assist the movement of the affected limb of stroke patients gently and flexibly, move the joints, exercise the muscles, and achieve the purpose of auxiliary rehabilitation treatment [21]. Therefore, the exoskeleton manipulator should not only have a certain degree of rigidity, but also have good flexibility. It should have a rigid skeleton, soft joints, and flexible movements like human limbs to ensure safety and comfort while wearing the exoskeleton [22,23].
To achieve soft joints and flexible movements without installing real physical SEA and VSA in devices, we proposed a flexible control strategy for an upper-limb rehabilitation exoskeleton based on the virtual spring damper hypothesis (VSDH). It mainly includes the VSDH flexible control strategy for single-rigid-joint actuators and the VSDH flexible cooperative control strategy for multi-joint exoskeletons with redundant degrees of freedom (DoF). The innovations and contributions of the proposed methods are fourfold: (1) We propose an innovative concept and method of "virtual softening control strategy" for the rigid joint. It is assumed that a virtual spring damper (VSD) module is installed after the joint motor. By designing the control signal of the motor, the torque output of the joint actuator is softened so that it has the typical characteristics of elasticity and variable stiffness. (2) In this paper, the torque virtual tracking control method based on hyperbolic tangent function (HTF) is applied to deal with the torque control signal so that the output of the motor presents the characteristic of flexibility, avoids mechanical damage to the motor and reducer, and enhances the comfort of human-robot interaction. (3) We propose a trajectory planning method for a point-to-point position-tracking controller based on a normal distribution function. It is successfully applied to the control of the upper-limb exoskeleton. (4) A flexible multi-joint cooperative controller based on VSDH is successfully applied to solve the constrained control problem of the exoskeletons and the self-motion problem caused by redundant DoFs.
The remainder of the paper is organized as follows. Section 2 introduces the flexible control strategy based on VSDH for a single-rigid-joint actuator. The multi-joint cooperative flexible control strategy based on VSDH is proposed, and the feasibility is theoretically proven, in Section 3. Section 4 presents the simulation experiment results. The conclusions are presented in Section 5.

Single-Joint Flexible Based on VSDH
Between the rigid joint motor and the output shaft, it is assumed that a virtual spring damper (VSD) is installed in series to compose a flexible driving module. This module can be regarded as a virtualized flexible joint. The flexible effect is equivalently simulated by a linear spring and a damper in parallel, as shown in Figure 1, where τ m is the output torque of the motor; θ m is the angle of the motor; J m is the equivalent moment of inertia; B m is the total viscous damping coefficient; k s is the stiffness coefficient of virtual torsional spring; B s is the damping coefficient of the virtual damper; τ q is torque at the output end of the flexible VSD module; q is the angle of the output shaft; J q is the effective moment of inertia of the load shaft; B q is damping coefficient of the load shaft; and τ out is external torque.
Actuators 2022, 11, x FOR PEER REVIEW 3 of 24 moment of inertia of the load shaft; is damping coefficient of the load shaft; and is external torque.

Velocity Tracking Control of VSD Flexible Joint
The Laplace transform is performed on both sides of Equation (2): where Θ ( ), ( ), and ( ) represent the Laplace transformation of , , and , respectively. Thus, the transfer function can be obtained as follows: A structural block diagram of the control system of the virtual flexible joint is shown in Figure 2. Assume that the joint motor is driven by a DC servo motor controlled by armature voltage. The dynamic equation of the motor driver is as follows: According to Hooke's law and torque analysis, the output torque of the flexible module can be obtained as follows: According to Newtonian dynamics, the dynamic equation of the flexible joint can be obtained as follows:

Velocity Tracking Control of VSD Flexible Joint
The Laplace transform is performed on both sides of Equation (2): where Θ m (s), Q(s), and T m (s) represent the Laplace transformation of θ m , q, and τ m , respectively. Thus, the transfer function can be obtained as follows: A structural block diagram of the control system of the virtual flexible joint is shown in Figure 2. According to Hooke's law and torque analysis, the output torque of the flexible module can be obtained as follows: According to Newtonian dynamics, the dynamic equation of the flexible joint can be obtained as follows:

Velocity Tracking Control of VSD Flexible Joint
The Laplace transform is performed on both sides of Equation (2): where Θ ( ), ( ), and ( ) represent the Laplace transformation of , , and , respectively. Thus, the transfer function can be obtained as follows: A structural block diagram of the control system of the virtual flexible joint is shown in Figure 2. Assume that the joint motor is driven by a DC servo motor controlled by armature voltage. The dynamic equation of the motor driver is as follows: Assume that the joint motor is driven by a DC servo motor controlled by armature voltage. The dynamic equation of the motor driver is as follows: where u(t) is the control voltage; i(t) is the actual current of motor armature coil; L is the inductance of the armature coil; R is the armature coil resistance; k b is the electromagnetic torque coefficient; . θ m is the angular velocity of the motor shaft; and k m is the motor torque coefficient.  (3) and (5), the system structure diagram of the servo-motordriven virtual flexible joint speed-tracking control can be designed as shown in Figure 3. The q * in the figure represents the reference input.
where ( ) is the control voltage; ( ) is the actual current of motor armature coil; is the inductance of the armature coil; is the armature coil resistance; is the electromagnetic torque coefficient; is the angular velocity of the motor shaft; and is the motor torque coefficient.
Combined with Equations (3) and (5), the system structure diagram of the servo-motor-driven virtual flexible joint speed-tracking control can be designed as shown in Figure  3. The * in the figure represents the reference input. A virtual spring-damper flexible module added in a joint actuator will not only produce time delay and reduce the response speed, but also increase the vibration of the joint, like adding a real spring. An angular-velocity-tracking controller with good performance can be designed by properly selecting the system parameters. In our case, the system parameters are assumed as = = 0.05; = 0.3; and = 0.3. The VSD parameters are set at = 10.0 and = 1.012, based on a large number of simulation experiments. If the input signal is a rectangular wave with an amplitude of 2, the responses of the four state variables , , , are shown in Figure 4. It can be seen from Figure 4 that the motor velocity (red) and the load speed (black) can track the driving torque ( in the figure) very well. The delay time is relatively small, while the response is fast. Most importantly, it can be observed that there is no vibration during the tracking simulation. It not only makes the joint flexible, but also avoids undesired vibration. A virtual spring-damper flexible module added in a joint actuator will not only produce time delay and reduce the response speed, but also increase the vibration of the joint, like adding a real spring. An angular-velocity-tracking controller with good performance can be designed by properly selecting the system parameters. In our case where ( ) is the control voltage; ( ) is the actual current of motor armature coil; is the inductance of the armature coil; is the armature coil resistance; is the electromagnetic torque coefficient; is the angular velocity of the motor shaft; and is the motor torque coefficient.
Combined with Equations (3) and (5), the system structure diagram of the servo-motor-driven virtual flexible joint speed-tracking control can be designed as shown in Figure  3. The * in the figure represents the reference input. A virtual spring-damper flexible module added in a joint actuator will not only produce time delay and reduce the response speed, but also increase the vibration of the joint, like adding a real spring. An angular-velocity-tracking controller with good performance can be designed by properly selecting the system parameters. In our case, the system pa-  It can be seen from Figure 4 that the motor velocity (red) and the load speed (black) can track the driving torque ( in the figure) very well. The delay time is relatively small, while the response is fast. Most importantly, it can be observed that there is no vibration during the tracking simulation. It not only makes the joint flexible, but also avoids undesired vibration. It can be seen from Figure 4 that the motor velocity . θ m (red) and the load speed . q (black) can track the driving torque (τ m in the figure) very well. The delay time is relatively small, while the response is fast. Most importantly, it can be observed that there is no vibration during the tracking simulation. It not only makes the joint flexible, but also avoids undesired vibration.

Torque Tracking Control Based on HTF
When the control torque suddenly increases or disappears, the hyperbolic tangent function (HTF) is used to deal with the torque control signal to avoid mechanical damage to the motor and reducer and to enhance the flexibility of the exoskeleton and the comfort of human-robot interaction.
It is assumed that the output torque is τ m0 at the current moment k, and it needs to be suddenly increased to τ md at the next moment (k + 1). The objective of the torque control signal process is to make the output elastic and flexible like a spring, rather than a sudden change without any delay. There are two signal processing methods: The function f (τ m ) is defined as the map from the torque boost interval [τ m0 ,τ md ] to [−5, 5], which is the independent variable change interval of tan h( * ).
Since the range interval of the function tan h( * ) is [−1, 1], the control torque can be designed as follows:
It is assumed that the control torque τ q of the virtual flexible driving module is stepby-step increased from τ m0 to τ md in N steps, that is, in N signal-sampling period. The control torque can be designed as follows: where k = 0, 1, 2 . . . . . . N; α = N/2; and N can be set to 2, 4, 6, 8, or 10 as required. The value of the independent variable f (k, α) of HTF is shown in Table 1. The driving signal τ m is set as a rectangular wave signal [0→2→−2→2], and N is chosen as 2, 4, 6, 8, and 10; the resulting tracking curve of τ q tracking τ m is shown in Figure 4. Tm denotes the rectangular wave driving signal τ m . Tq(α) denotes the virtual flexible driven torque τ q (α) when α = 1, 2 . . . . . . 5. When α = 5, the number of steps is N = 2α = 10, that is, in 10 control command-sending periods, the output torque of the virtual flexible driving module is increased from τ m0 = 0 to τ md = 2, or decreased from τ m0 = 2 to τ md = −2, or from τ m0 = −2 to τ md = 2. If α = 1, it is only equivalent to inserting an intermediate value in the interval [τ m0 , τ md ], from one-step arrival to two-step arrival. In our case, α can be selected as 4 or 5.
It can be seen from the structure diagram ( Figure 2) that the load module can be simplified to a first-order inertial system, and the transfer function can be simplified as follows: where Q(s) and T q (s) represent the Laplace transformation of q and τ q , respectively. J q and B q are the total effective moment of inertia and the total viscous damping coefficient, respectively. For a joint motor, J q = 5.468 * 10 −4 , B q = 4.932 * 10 −3 , α = 5 is selected, τ q is calculated by (5), the rectangular wave response of angular q and angular velocity . q is as shown in Figure 5, and τ m and τ q (5) in Figure 5 are the same as in Figure 4. q 1 and D q1 represent the angular q and angular velocity . q driven by τ m , respectively. q 2 and D q2 represent q and . q driven by τ q (5), respectively. It can be observed from Figure 6 that the velocity-tracking curve D q2 is smoother with less peak compared to D q1 . A similar situation applies to the angle-tracking curve q 2 compared to q 1 . lows: where ( ) and ( ) represent the Laplace transformation of and , respectively. and are the total effective moment of inertia and the total viscous damping coefficient, respectively. For a joint motor, = 5.468 * 10 , = 4.932 * 10 , α = 5 is selected, is calculated by (5), the rectangular wave response of angular and angular velocity is as shown in Figure 5, and and (5) in Figure 5 are the same as in Figure 4. and represent the angular and angular velocity driven by , respectively. and represent and driven by (5), respectively. It can be observed from Figure 6 that the velocity-tracking curve is smoother with less peak compared to . A similar situation applies to the angle-tracking curve compared to . According to the discussion and simulation results above, it can be concluded that the torque virtual tracking control method based on HTF (4) and (5) can smooth the peak of the response curve, which plays a role as soft buffering, protecting the motor and reducer. By properly selecting the control parameters, the response speed of the joint can be sacrificed as little as possible to ensure the fast tracking of torque, speed, and position. There is no undesired oscillation characteristic, which a real spring flexible joint would produce.
The proposed flexible control strategy for a single rigid joint can not only realize the flexible tracking of the joint, but can also avoid undesired vibration. The simulation results confirm the effectiveness and feasibility of the control method.

Dynamic Equation and Control of Upper-Limb Exoskeleton
A prototype of a 5-DoF upper-limb exoskeleton manipulator was designed in our previous work [24], shown in Figure 7. It can realize movements including shoulder joint abduction/adduction, shoulder joint flexion/extension, shoulder joint internal/external rotation, elbow joint flexion/ extension, forearm supination/pronation, and wrist joint flexion/extension. The exoskeleton was used as the control object.
The length of the upper arm of the exoskeleton manipulator is adjustable, ranging from 22 to 32 cm. The total weight of the upper arm including joint-3 and joint-4 is 4 kg. The length of the forearm can be adjusted from 25 cm to 35 cm. The total weight of the forearm, including joint-5 and wrist, is 3 kg. See Table 2 for the rotation angle range of each joint. According to the discussion and simulation results above, it can be concluded that the torque virtual tracking control method based on HTF (4) and (5) can smooth the peak of the response curve, which plays a role as soft buffering, protecting the motor and reducer. By properly selecting the control parameters, the response speed of the joint can be sacrificed as little as possible to ensure the fast tracking of torque, speed, and position. There is no undesired oscillation characteristic, which a real spring flexible joint would produce.
The proposed flexible control strategy for a single rigid joint can not only realize the flexible tracking of the joint, but can also avoid undesired vibration. The simulation results confirm the effectiveness and feasibility of the control method.

Dynamic Equation and Control of Upper-Limb Exoskeleton
A prototype of a 5-DoF upper-limb exoskeleton manipulator was designed in our previous work [24], shown in Figure 7. It can realize movements including shoulder joint abduction/adduction, shoulder joint flexion/extension, shoulder joint internal/external rotation, elbow joint flexion/extension, forearm supination/pronation, and wrist joint flexion/extension. The exoskeleton was used as the control object. module.

Dynamic Equation and Control of Upper-Limb Exoskeleton
A prototype of a 5-DoF upper-limb exoskeleton manipulator was designed in our previous work [24], shown in Figure 7. It can realize movements including shoulder joint abduction/adduction, shoulder joint flexion/extension, shoulder joint internal/external rotation, elbow joint flexion/ extension, forearm supination/pronation, and wrist joint flexion/extension. The exoskeleton was used as the control object.
The length of the upper arm of the exoskeleton manipulator is adjustable, ranging from 22 to 32 cm. The total weight of the upper arm including joint-3 and joint-4 is 4 kg. The length of the forearm can be adjusted from 25 cm to 35 cm. The total weight of the forearm, including joint-5 and wrist, is 3 kg. See Table 2 for the rotation angle range of each joint.  The length of the upper arm of the exoskeleton manipulator is adjustable, ranging from 22 to 32 cm. The total weight of the upper arm including joint-3 and joint-4 is 4 kg. The length of the forearm can be adjusted from 25 cm to 35 cm. The total weight of the forearm, including joint-5 and wrist, is 3 kg. See Table 2 for the rotation angle range of each joint. Table 2. Exoskeleton joint range of motion.

Joint Action Joint Code Angle Range
Shoulder-1 abduction/adduction Joint-1 0-100 The inertial frame OXYZ is established as shown in Figure 8. The origin of coordinate O is on the center point of joint-1. The OX axis is perpendicular to the coronal plane of the human body and points to the front of the exoskeleton, the OZ axis is outward along the shoulder in the horizontal plane, and the OY axis is perpendicular to the ground and pointing upward. Points A and B represent the center point of the elbow joint and wrist joint, respectively. The link OA (represents the upper arm) rotates an angle q 1 around the axis OX, then rotates q 2 and q 3 around the new axis OZ 1 (not shown in the figure) and axis OA, respectively. The link AB (represents the forearm) rotates q 4 in the plane OAB with respect to the elbow joint-4.
The inertial frame OXYZ is established as shown in Figure 8. The origin of coordinate O is on the center point of joint-1. The OX axis is perpendicular to the coronal plane of the human body and points to the front of the exoskeleton, the OZ axis is outward along the shoulder in the horizontal plane, and the OY axis is perpendicular to the ground and pointing upward. Points A and B represent the center point of the elbow joint and wrist joint, respectively. The link OA (represents the upper arm) rotates an angle around the axis OX, then rotates and around the new axis OZ1 (not shown in the figure) and axis OA, respectively. The link AB (represents the forearm) rotates in the plane OAB with respect to the elbow joint-4.   The Jacobian matrix of the system is derived from Equation (12) and listed in Appendix B.
The Lagrangian dynamic equation of the 5-DoF upper-limb exoskeleton manipulator shown in Figure 7 can be described as follows: where s i = sin(q i ), c i = cos(q i ), i = 1, 2, 3, 4 (the same as below). The Jacobian matrix of the system is derived from Equation (12) and listed in Appendix B.
The Lagrangian dynamic equation of the 5-DoF upper-limb exoskeleton manipulator shown in Figure 7 can be described as follows: where q, . q, .. q ∈ R 4 represent the angular position, velocity, and acceleration. τ ∈ R 4 is the joint torque. S q, . q represents Coriolis/centripetal matrix. G(q) ∈ R 4 is the gravity vector. The generalized inertia matrix M(q) ∈ R 4×4 is symmetric and positive definite. Then, there exist two positive constants h m and h M such that (for any q): The calculation formula of M(q), S q, . q , and G(q) is given in Appendix A. During the movement of the upper-limb exoskeleton manipulator from one position to another, shown in Figure 7, the cooperative rotation of four joints is constrained and limited. See Figure 9. OA and AB represent the upper arm and forearm, separately. B represents the starting position of the exoskeleton endpoint, and B' represents the target position.
The calculation formula of ( ), ( , ), and ( ) is given in Appendix A. During the movement of the upper-limb exoskeleton manipulator from one position to another, shown in Figure 7, the cooperative rotation of four joints is constrained and limited. See Figure 9. OA and AB represent the upper arm and forearm, separately. B represents the starting position of the exoskeleton endpoint, and B' represents the target position.

Constraint 1:
If q3 = 0, then the plane of triangle OAB is in the same plane as BOX, namely, the OX axis is in the plane of triangle OAB.
Constraint 2: The most convenient and comfortable movement path from B to B' is straight line BB', shown in Figure 9. The constraint condition is that the endpoint trajectory B(t) must move in the OBB' plane. If there is a point Q on the arc of the endpoint trajectory, the vertical point to OBB' is P, and QP length is d, then the constraint condition on the control law is to make d = 0.
The initial position coordinates of the endpoints B (xB0, yB0, zB0) and the target position B′ (xB′, yB′, zB′), any point on the endpoint trajectory Q (xB, yB, zB), and the normal vector of the OBB' plane are as follows: where , , represent the unit vector, and , , represents the corresponding coordinate values of three directions. Then, the distance from any point on the endpoint trajectory to the OBB' plane can be obtained as follows: Figure 9. Simplified sketch of exoskeleton manipulator coordinate system.

Constraint 1:
If q 3 = 0, then the plane of triangle OAB is in the same plane as BOX, namely, the OX axis is in the plane of triangle OAB.
Constraint 2: The most convenient and comfortable movement path from B to B' is straight line BB', shown in Figure 9. The constraint condition is that the endpoint trajectory B(t) must move in the OBB' plane. If there is a point Q on the arc BB of the endpoint trajectory, the vertical point to OBB' is P, and QP length is d, then the constraint condition on the control law is to make d = 0.
The initial position coordinates of the endpoints B (x B0 , y B0 , z B0 ) and the target position B (x B , y B , z B ), any point on the endpoint trajectory Q (x B , y B , z B ), and the normal vector of the OBB' plane are as follows: where i, j, k represent the unit vector, and Ψ i,j,k represents the corresponding coordinate values of three directions. Then, the distance from any point on the endpoint trajectory to the OBB' plane can be obtained as follows: The constraint condition is d = 0, that is, the numerator of the formula above is equal to zero. The constraint function is described as follows: The partial derivative (gradient) of φ(q) in q is described as follows: The representation of J φi (I = 1, 2, 3, 4) is shown in Appendix C. The function of the joint driving torque controller is to drive each joint of the exoskeleton to rotate and to make the endpoint move from an initial position B(x 0 , y 0 , z 0 ) to a given target position B (x d , y d , z d ). Based on the virtual spring-damper hypothesis, the flexible multi-joint cooperative controller under constraint is designed as follows: where τ denotes the control torque vector, C is a positive definite diagonal matrix whose diagonal components express virtual damping coefficients, J x (q) = ∂x/∂q is the 3 × 4 Jacobian matrix of x in q (presented in Appendix B), and ∆x = x − x d represents the difference between the real-time measured position x of the end point and the target position x d . K stands for a virtual stiffness parameter, and λ denotes the constraint force coefficient. J φ is the constraint function gradient matrix (presented in Appendix C), and G(q) represents the estimated value of the gravity term, which can be computed in real time based on measured angular data. Let us assume hereĜ(q) = G(q). Substituting (19) into (14), the closed-loop control system equation of the whole upperlimb exoskeleton manipulator can be described in the following form: The controller (19) has the following characteristics: (1) It has the dynamic characteristics of "spring and damper". Equation (19) contains four items. The first item represents the damping moment generated by the virtual damper in joint space, which plays a "braking" role in the system movement. The convergence speed of the system state is changed by designing the value of the damping coefficient matrix C. The second term represents the virtual spring torque. By selecting the stiffness coefficient K, the speed of the endpoint approaching the target position is adjusted. (2) It has solved the "self-motion problem" [25] caused by redundant DoFs. The dynamic Equation (13) defaults that each joint is a cardan joint, resulting in redundant DoFs, and will produce the "self -motion" shown in Figure 8. Two links OA-AB will rotate freely around the OB axis. The third term is the constraint term. It is added to the controller and can limit the arbitrary self-motion of the exoskeleton. (3) To directly solve the control problem from task space to joint space, it is not necessary to solve the pseudo-inverse matrix and inverse kinematics. By selecting the appropriate controller parameters C and K, the state of the system can be driven to converge gradually. That is, when t → ∞ , not only x(t) → x d , but also q(t) → q 0 , . q(t) → 0 , where q 0 represents the most comfortable nominal position of the joint angle when the endpoint reaches the target position in the task space.
The next key problem is to prove that the closed-loop control system Equation (14) is stable, and to verify that the controller is available by use of simulation analysis.

Stability Analysis
The state variables of the closed-loop control system in Equation (13)  q in the joint space, but also ∆x in the task space. Therefore, it is difficult to find the equilibrium points of the system, and it is more difficult to analyze the stability of the equilibrium point by Lyapunov method. Arimoto et al. [26] proposed some concepts related to the reference equilibrium state and its neighborhood, including a six-dimensional manifold in the neighborhood and its stability on the manifold, which solved the problem of stability analysis of the closed-loop control system in Equation (20) well.
The method and idea of stability analysis are as follows: to prove that as t → ∞ , x(t) → x d , then q(t) → q 0 , . q(t) → 0 , just prove that any solution of Equation (20) is stable on a manifold and asymptotically tends to a sub-manifold. Since q, taking an inner product in Equation (20) with . q leads to the following [26]: where E(q, E stands for the total energy of the system, including the quadratic function of . q and ∆x(q) = x(q) − x d . In fact, the scalar function E q, . q is not positive definite with respect to the state vector q, . q ∈ R 8 , because E includes only a quadratic term of threedimensional vector ∆x and a positive definite quadratic function of . q. Therefore, it cannot be used as a Lyapunov function to analyze the stability of the closed-loop control system in Equation (20).
We define the one-dimensional manifold as follows: The constrained 6-dimensional manifold is as follows: We consider a reference equilibrium state posture q 0 , . q 0 = 0 ∈ R 8 with still state (i.e., . q = 0), whose endpoint is located at x d , that is, x q 0 = x d , and hence q 0 , 0 ∈ H 1 . We define a neighborhood N 8 (r) of the reference state q 0 , 0 with radius r in such a manner that: where ∆q = q − q 0 . We define a family of neighborhoods of the reference state q 0 , 0 on the constraint manifold H 6 in such a way that: where E is shown in Equation (23). Through careful analysis and understanding of these manifolds and regions, it can be obtained that any state lying on H 1 ∩ N 8 (r 0 ) is included in N 6 (ε) ∩ N 8 (r 0 ) for any ε ≥ 0. Now, let us consider that the reference equilibrium state q 0 , 0 is stable on the manifold H 1 ∩ N 8 (r 0 ). It is possible to choose a positive parameter r 0 > 0 for Equation (17) and construct a neighborhood N8(r0) such that for any q, . q ∈ N 8 (r 0 ), J x (q) = (∂x B /∂q, ∂y B /∂q, ∂z B /∂q) T is nondegenerate. Then, there exist two positive constants σ 0 and σ 1 such that: where I 3 is the 3 × 3 unit matrix. [27]), and J φ = ∂φ(q)/∂q (see Appendix C). The scalar function E, defined in Equation (19), does not increase as t → ∞ , so it is possible to choose some positive parameter δ 1 such that the following inequality holds: Now, let us consider another extended scalar function [27]: where γ > 0 is a constant that is determined later. The time-derivative of V is as follows: The above formula looks very complex. In fact, the fourth term is the sum of the quadratic function of the angular velocity vector q, that is, they vanish when . q = 0. Therefore, there exists a constant α 0 such that the following formula holds: The third item in Equation (31) meets the following condition: The first three items in Equation (31) can be expressed as follows: .
Substituting Equations (31) and (32) into Equation (30), considering that CP φ C ≤ C 2 , selecting parameters C = cI 4 , K = kI 3 , and noting that h M is the maximum eigenvalue of M(q) over all q, Equation (30) can be reduced to the following: Reasonably selecting parameters 0.004 ≤ c ≤ 0.1, γ = k, and δ 0 ,h M , it is possible to reach the following: We construct an inequality |ab| ≤ (1/8)a 2 + 2b 2 and use it to analyze the inequality of scalar function V from another perspective. and similarly: It follows from Equations (35)-(37) that: .

Simulation Conditions
The upper-limb exoskeleton manipulator introduced in Section 3 and shown in Figure 7 is used as the simulation research object. The parameter matrix model of the closed-loop control system Equation (13) is shown in Appendix A. We note Assumption 1 and 2 and assume that S q,

A Velocity Driving Control Method Base on SNDF
In order to compare the simulation results of different control strategies, we propose an innovative velocity drive control method based on standard normal distribution functions (SNDF). This idea comes from the research conclusion in [26], that is, human skilled multijoint reaching is characterized as follows: (a) the profile of the endpoint trajectory in task space becomes closely rectilinear, (b) the velocity profile of it becomes symmetric and bell-shaped, and (c) the acceleration profile has double peaks. As we all know, the SNDF curve is a bell shape. Therefore, we construct the driving velocity control signal of the closed-loop control system Equation (14) by use of SNDF. The design steps are as follows: Step 1: Determine the initial posture q(0) and target posture q(∞); Step 2: Map the interval [q(0), q(∞)] to the independent variable interval [-4, +4] of SNDF, and construct the angular velocity . q(t) as follows: .
Step 3: The angular velocity Equation (47) is used as the driving angular velocity of the joint motor. According to the properties of the SNDF, the area enclosed by the function curve and the independent variable coordinate axis is equal to 1. Thus, we can get: Then, it can be found that: Function Ω = ∆q . q(t) can be used to construct the driving control torque of joints, as follows:

Simulation Results
Firstly, the simulation results under the action of Equation (50) are shown in Figure 10. Figure 10 shows the four constructed angular velocities Ω i , i = 1, 2, 3, 4. Those four curves are obviously applied to standard normal distribution. Figure 11 shows the change of posture and the endpoint trajectories in Cartesian space. The initial positions of the upper arm and forearm are located at OA and AB, respectively, while the end positions are at OA' and A'B'. x A (t) represents the trajectory of point A, while x B (t) represents the trajectory of endpoint B. From the simulation results, the target position is accurately reached, and the control task is fulfilled.

Simulation Results
Firstly, the simulation results under the action of Equation (50) are shown in Figure  10. Figure 10 shows the four constructed angular velocities , = 1, 2, 3, 4. Those four curves are obviously applied to standard normal distribution. Figure 11 shows the change of posture and the endpoint trajectories in Cartesian space. The initial positions of the upper arm and forearm are located at OA and AB, respectively, while the end positions are at OA' and A'B'. ( ) represents the trajectory of point A, while ( ) represents the trajectory of endpoint B. From the simulation results, the target position is accurately reached, and the control task is fulfilled.  Secondly, let us simulate the control effect of Equation (19) without constraint item , and assuming that the gravity term is fully compensated, that is:

Simulation Results
Firstly, the simulation results under the action of Equation (50) are shown in Figur 10. Figure 10 shows the four constructed angular velocities , = 1, 2, 3, 4. Those fou curves are obviously applied to standard normal distribution. Figure 11 shows the chang of posture and the endpoint trajectories in Cartesian space. The initial positions of th upper arm and forearm are located at OA and AB, respectively, while the end position are at OA' and A'B'. ( ) represents the trajectory of point A, while ( ) represent the trajectory of endpoint B. From the simulation results, the target position is accuratel reached, and the control task is fulfilled.  Secondly, let us simulate the control effect of Equation (19) without constraint item , and assuming that the gravity term is fully compensated, that is:  Secondly, let us simulate the control effect of Equation (19) without constraint item λJ T φ , and assuming that the gravity term is fully compensated, that is:  Figure 12, denoted by Dq1, Dq2, Dq3, and Dq4, respectively. Unlike the velocity distribution shown in Figure 10, the angular velocity increases rapidly at the beginning, reaching the maximum peak in 0.5 s. The posture and endpoint trajectories under the action of Equation (50) (green curve) and Equation (51) (blue curve) are shown in Figure 13. Obviously, despite starting from the same initial point, the two endpoint trajectories are completely different. Figure 14 shows the change process of the four joint angles with two different control methods: the curve qi-1 indicates (50), while the curve qi-2 indicates (51). As can be seen from Figure 14, the angles change processes of solid (50) and dotted lines (51) of the same color are completely different, despite starting from the same initial angle (−30 • 20 • , 0 • , 60 • ). In particular, the change trend of angle q 1 is completely opposite: one positive (q 1-1 ) turn, and the other reverse (q 1-2 ). The redundant DoFs in the system lead to the rotation of the angle from another direction in the process of tracking the change of spatial position. Figure 15 shows the change process of the endpoint position coordinates x b (t), y b (t), and z b (t). The solid line represents method in Equation (50), and the dashed line represents method in Equation (51). Under the action of the two control methods, the movement trends of the two endpoints are consistent, although one speed-up is earlier, and the other is later. shown in Figure 12, denoted by 1, 2, 3, and 4, respectively. Unlike the velocity distribution shown in Figure 10, the angular velocity increases rapidly at the beginning, reaching the maximum peak in 0.5 s. The posture and endpoint trajectories under the action of Equation (50) (green curve) and Equation (51) (blue curve) are shown in Figure 13. Obviously, despite starting from the same initial point, the two endpoint trajectories are completely different. Figure 14 shows the change process of the four joint angles with two different control methods: the curve i-1 indicates (50), while the curve i-2 indicates (51). As can be seen from Figure 14, the angles change processes of solid (50) and dotted lines (51) of the same color are completely different, despite starting from the same initial angle (−30° 20°, 0°, 60°). In particular, the change trend of angle q1 is completely opposite: one positive ( 1-1) turn, and the other reverse ( 1-2). The redundant DoFs in the system lead to the rotation of the angle from another direction in the process of tracking the change of spatial position. Figure 15 shows the change process of the endpoint position coordinates    shown in Figure 12, denoted by 1, 2, 3, and 4, respectively. Unlike the veloc ity distribution shown in Figure 10, the angular velocity increases rapidly at the beginnin reaching the maximum peak in 0.5 s. The posture and endpoint trajectories under the a tion of Equation (50) (green curve) and Equation (51) (blue curve) are shown in Figure 1 Obviously, despite starting from the same initial point, the two endpoint trajectories ar completely different. Figure 14 shows the change process of the four joint angles with tw different control methods: the curve i-1 indicates (50), while the curve i-2 indicates (51 As can be seen from Figure 14, the angles change processes of solid (50) and dotted line (51) of the same color are completely different, despite starting from the same initial ang (−30° 20°, 0°, 60°). In particular, the change trend of angle q1 is completely opposite: on positive ( 1-1) turn, and the other reverse ( 1-2). The redundant DoFs in the system lea to the rotation of the angle from another direction in the process of tracking the change o spatial position. Figure 15 shows the change process of the endpoint position coordinate     Thirdly, let us simulate the control effect of Equation (19). In the control method in Equation (19), the constraint term = ( / ) is added to the control law to ensure that the motion trajectory of the endpoint is kept within the plane formed by the initial position (0) and the target position (∞) as much as possible. It could be found that the system performs better when the elements of diagonal matrix C and K are selected as ( , , , ) = (8,10,6,2) , ( , , , ) = (25,35,30,20), and = 50 . The simulation results are shown in Figures 16 and 17. The variation tendencies of the four joint angles for the control method in Equation (51) and Equation (19) are shown in Figure  16. The dashed curves are the same as in Figure 14 and represent the results of control method in Equation (51). The solid curves in the Figure 16 represent the results of control method in Equation (19). Figure 17 shows the joint posture changes and endpoint trajectories under three control methods. The curves (by Equation (50)), (by Equation (51)), and (by Equation (19)) represent the three trajectories of endpoints under three different control methods, respectively.  Thirdly, let us simulate the control effect of Equation (19). In the control method in Equation (19), the constraint term = ( / ) is added to the control law to ensure that the motion trajectory of the endpoint is kept within the plane formed by the initial position (0) and the target position (∞) as much as possible. It could be found that the system performs better when the elements of diagonal matrix C and K are selected as ( , , , ) = (8,10,6,2) , ( , , , ) = (25,35,30,20), and = 50 . The simulation results are shown in Figures 16 and 17. The variation tendencies of the four joint angles for the control method in Equation (51) and Equation (19) are shown in Figure  16. The dashed curves are the same as in Figure 14 and represent the results of control method in Equation (51). The solid curves in the Figure 16 represent the results of control method in Equation (19). Figure 17 shows the joint posture changes and endpoint trajectories under three control methods. The curves (by Equation (50)), (by Equation (51)), and (by Equation (19)) represent the three trajectories of endpoints under three different control methods, respectively. Thirdly, let us simulate the control effect of Equation (19). In the control method in Equation (19), the constraint term J T φ λ d = ∂ φ /∂ q T λ d is added to the control law to ensure that the motion trajectory of the endpoint is kept within the plane formed by the initial position x(0) and the target position x(∞) as much as possible. It could be found that the system performs better when the elements of diagonal matrix C and K are selected as (c 11 , c 22 , c 33 , c 44 ) = (8, 10, 6, 2), (k 11 , k 22 , k 33 , k 44 ) = (25,35,30,20), and λ d = 50. The simulation results are shown in Figures 16 and 17. The variation tendencies of the four joint angles for the control method in Equation (51) and Equation (19) are shown in Figure 16. The dashed curves q i−2 are the same as in Figure 14 and represent the results of control method in Equation (51). The solid curves q i−1 in the Figure 16 represent the results of control method in Equation (19). Figure 17 shows the joint posture changes and endpoint trajectories under three control methods. The curves BB (by Equation (50)), BB (by Equation (51)), and BB (by Equation (19)) represent the three trajectories of endpoints under three different control methods, respectively.

Comparative Analysis
By analyzing the above simulation results, the following conclusions can be drawn: (1) The system in Equation (20) (19) are completely correct; this trend is the same as tha of Equation (50) and completely different from that of Equation (51). Comparative analysis of the three endpoint trajectories is shown in Figure 17. The trajectory o ( ) controlled by Equation (19) is closely rectilinear, and the accuracy is the highest

Comparative Analysis
By analyzing the above simulation results, the following conclusions can be drawn: (1) The system in Equation (20) Figure 13, it can be seen that the trajectories are completely different. (3) The constraint term added to the controller limits the arbitrary self-motion o the exoskeleton manipulator. It can be seen from Figure 16 that the posture paths o ( ) controlled by Equation (19) are completely correct; this trend is the same as tha of Equation (50) and completely different from that of Equation (51). Comparativ analysis of the three endpoint trajectories is shown in Figure 17. The trajectory o ( ) controlled by Equation (19) is closely rectilinear, and the accuracy is the highest

Comparative Analysis
By analyzing the above simulation results, the following conclusions can be drawn: (1) The system in Equation (20) q(t) → 0 . That is to say that q(t) and . q(t) can both be controlled by ∆x B . (2) The "self-motion problem" caused by redundant DoFs will lead to the uncertainty of the posture paths. Two control methods in Equations (50) and (51), achieved completely different posture change paths and endpoint trajectories for the same control target ∆x B → 0 . It can be seen from Figure 14 that the trend of solid and dotted lines of the same color is completely different. The joint angle changes monotonically from q(0) to q(∞) under the control of Equation (50). However, the joint angle under the control of Equation (51) varies according to an uncertain law and finally reaches an uncertain angle position, although the control objectives of ∆x B ≈ 0 and . q(t) ≈ 0 have been achieved already. Compared with the two trajectoriesB B andB B of endpoint B in Figure 13, it can be seen that the trajectories are completely different.
(3) The constraint term J T φ λ d added to the controller limits the arbitrary self-motion of the exoskeleton manipulator. It can be seen from Figure 16 that the posture paths of q(t) controlled by Equation (19) are completely correct; this trend is the same as that of Equation (50) and completely different from that of Equation (51). Comparative analysis of the three endpoint trajectories is shown in Figure 17. The trajectory BB of q(t) controlled by Equation (19) is closely rectilinear, and the accuracy is the highest.

Experimental Verification
We use the prototype of an upper-limb rehabilitation exoskeleton designed by our team, as shown in Figure 7, to verify the feasibility and the reliability of the VSDH algorithm (19). A 24-year-old male volunteered to conduct the experiment. The initial posture of the exoskeleton is vertical down, as shown in Figure 18a, and the end position coordinates are (44, −30, 17) (cm). The target position is above the front of the body, shown in Figure 18b, and the end position coordinates are (50, 10, −30) (cm). In one experiment, an action of lifting and falling was performed. The initial parameters of the controller are set as (c11, c22, c33, c44) = (12.85, 16.6, 3, 8.7), (k11, k22, k33, k44) = (20,30,40,30), and λ = 0.80. Collect the voltage signal data output by the encoder that comes with the servo motor (sampling period is 0.05 s) and convert it into joint angle data. A set of four joint angle data sets is plotted, as shown in Figure 19. The dotted lines sq1, sq2, sq3, and sq4 in Figure 19 represent the digital simulation data curves of the four joint angles, respectively, the tortuous solid lines q1, q2, q3, and q4 represent the measured data curves of the four joint angles, respectively, and the smoother solid lines represent the data curve after filtering. It can be seen from Figure 19 that the experimental data are highly close to the simulation results. The error between them is caused by the motion uncertainty of the human upper limb.

Experimental Verification
We use the prototype of an upper-limb rehabilitation exoskeleton designed by ou team, as shown in Figure 7, to verify the feasibility and the reliability of the VSDH algo rithm (19). A 24-year-old male volunteered to conduct the experiment. The initial postur of the exoskeleton is vertical down, as shown in Figure 18a, and the end position coord nates are (44, −30, 17) (cm). The target position is above the front of the body, shown i Figure 18b, and the end position coordinates are (50, 10, −30) (cm). In one experiment, a action of lifting and falling was performed. The initial parameters of the controller are se as (c11, c22, c33, c44) = (12.85, 16.6, 3, 8.7), (k11, k22, k33, k44) = (20,30,40,30), and λ 0.80. Collect the voltage signal data output by the encoder that comes with the servo moto (sampling period is 0.05 s) and convert it into joint angle data. A set of four joint angl data sets is plotted, as shown in Figure 19. The dotted lines 1, 2, 3, and 4 i Figure 19 represent the digital simulation data curves of the four joint angles, respectively the tortuous solid lines q1, q2, q3, and q4 represent the measured data curves of the fou joint angles, respectively, and the smoother solid lines represent the data curve after fi tering. It can be seen from Figure 19 that the experimental data are highly close to th simulation results. The error between them is caused by the motion uncertainty of th human upper limb.
The overall feeling of the experimenter is that the exoskeleton runs smoothly, th acceleration and stop are gentle and natural, and there is no discomfort of sudden change Repeated experiments verified that the VSDH algorithm in Equation (19) has good feas bility and reliability.

Experimental Verification
We use the prototype of an upper-limb rehabilitation exoskeleton designed by our team, as shown in Figure 7, to verify the feasibility and the reliability of the VSDH algorithm (19). A 24-year-old male volunteered to conduct the experiment. The initial posture of the exoskeleton is vertical down, as shown in Figure 18a, and the end position coordinates are (44, −30, 17) (cm). The target position is above the front of the body, shown in Figure 18b, and the end position coordinates are (50, 10, −30) (cm). In one experiment, an action of lifting and falling was performed. The initial parameters of the controller are set as (c11, c22, c33, c44) = (12.85, 16.6, 3, 8.7), (k11, k22, k33, k44) = (20,30,40,30), and λ = 0.80. Collect the voltage signal data output by the encoder that comes with the servo motor (sampling period is 0.05 s) and convert it into joint angle data. A set of four joint angle data sets is plotted, as shown in Figure 19. The dotted lines 1, 2, 3, and 4 in Figure 19 represent the digital simulation data curves of the four joint angles, respectively, the tortuous solid lines q1, q2, q3, and q4 represent the measured data curves of the four joint angles, respectively, and the smoother solid lines represent the data curve after filtering. It can be seen from Figure 19 that the experimental data are highly close to the simulation results. The error between them is caused by the motion uncertainty of the human upper limb.
The overall feeling of the experimenter is that the exoskeleton runs smoothly, the acceleration and stop are gentle and natural, and there is no discomfort of sudden change. Repeated experiments verified that the VSDH algorithm in Equation (19) has good feasibility and reliability.    The overall feeling of the experimenter is that the exoskeleton runs smoothly, the acceleration and stop are gentle and natural, and there is no discomfort of sudden change. Repeated experiments verified that the VSDH algorithm in Equation (19) has good feasibility and reliability.

Conclusions
An upper-limb rehabilitation exoskeleton should have a rigid skeleton, soft joints, and flexible movements like human limbs so that it can meet ever-increasing demands for high flexibility, high environmental adaptability, and good comfort. Flexible joints have been successfully utilized in rehabilitation robots and shown some benefits. The flexible control strategy based on VSDH proposed in this paper was successfully applied to the controller design for an upper-limb rehabilitation exoskeleton to solve the virtual softening problem of a single-rigid-joint actuator and the constrained control problem of a multi-joint exoskeleton with redundant DoFs. Theoretical design and simulation results confirm that it is feasible to soften the torque output of the single-joint actuator by designing the control signal of the motor according to VSDH. The multi-joint cooperative flexible controller in Equation (19) presented the dynamic characteristics of spring damping and has solved the self-motion problem caused by redundant DoFs.
The control methods proposed in this paper are mainly suitable for passive rehabilitation training with an upper-limb exoskeleton. The initial pose and target position are set artificially. In future work, the focus of this project is to carry out active rehabilitation training according to the patient's active exercise intention. Apart from this issue, the unmeasurable noise caused by system vibration would also influence the precision of the controller; thus, we will be focusing on solving that problem too.
Author Contributions: All authors have made substantial contributions to the research design. Material preparation, data collection, and analysis were performed by D.K. Design of the prototype was performed by W.W. and Y.S. The first draft of the manuscript was written by D.K. and L.K. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
We thank the participants in this study for their willingness to devote energy and sweat to completing the research work of this project.

Conflicts of Interest:
The authors declare that they have no competing interests.