The system is shown in Figure 1, which includes three degrees of freedom (DOFs) upper limb rehabilitation robot (only two degrees of freedom are adopted in this paper), and an extensive visual interface (Zhang M. M. et al., 2023). The robot mainly consists of three components, the three-axis force and position sensors, the control box and the two handles that can move in three-dimensional space. The real-time sensors’ position and force data can be sampled and recorded. Meanwhile, the real-time position will be displayed on the extensive visual interface in front of subjects to provide subjects with intuitive vision.

The robotic working mode is the human-active mode when the subject dominantly drives the robotic handles’ according to his subjective wishes and the adjustment of the robot handle’s damping coefficient can make both hands in different motion states. This work mode can significantly enhance the participation of subjects to improve motor function (Luo et al., 2019). The bilateral damping coefficients can be adjusted to change subjects’ upper limb movement difficulty.

The damping coefficient k (The healthy side’s k H decides the challenge level during training, and the impaired side’s damping coefficient k I will be adjusted to optimize the mirror-symmetry effects under a challenge level.) are generated according to our proposed supervised method. The max damping coefficient is limited by the maximum speed of the motor.

The dynamic model of the robotic system interacted can be written by M q ¨ + K q ˙ = F h (1) where M = b l o c k d i a g M 1 , M 2 ∈ R 2 n × 2 n with M i = d i a g   m i 1 , . . . , m i n , i = 1 , 2 is the inertial matrix at end-effector, i = 1 and i = 2 separately denote left and right robotic handle, n = 1 to n = 3 respectively represent X-axis to Z-axis. q = q 1 T , q 2 T T ∈ R 2 n × 1 with q i = q i 1 . . . , q i n T ∈ R n × 1 is the position of end-effector in Cartesian space. F h = F h 1 T , F h 2 T T ∈ R 2 n × 1 with F h i ∈ R n × 1 is the interaction force applied on the end-effector by subjects, which can be measured by the force sensor.

In our robotic mirror rehabilitation, the damping coefficients of the handles can be adjusted to minimize the dissimilarity of bilateral kinematic performance. It is assumed that the subject’s unilateral limb movement ability is impaired (an impaired upper limb and a healthy limb). The optimization objective function can be defined as L 1 = ∑ i n 1 2 x i I − x i H 2 (2) where the kinematic performance feature vector X ∈ R n , x i I and x i H refer to the i t h feature of impaired upper limb performance X I , the i t h feature of healthy upper limb performance X H .

The healthy limb, as the master, generates the appropriate kinematic performance X H under the fixed damping coefficient. The damping coefficient of the healthy limb can be modified to change the challenge level according to the training performance, which is friendly to supervisors. The fixed damping coefficient, selected according to the personalized conditions of subjects, determines the upper limit and bottom limit of the reference kinematics performance. A large damping coefficient makes it difficult for the impaired limb to keep up with the healthy limb’s movement, and a small damping coefficient maybe let subjects lose interest in training. When the impaired limb feels uncomfortable during training, the healthy limb can exert more muscle force to improve the performance to reduce the difficulty of the impaired limb under the framework.

The impaired limb, as the slave, can keep up with the performance of the healthy limb with the assistance of the damping adjustment strategy based on proposed model. Figure 2 illustrates the proposed framework as a damping adjustment strategy. The input setting k I of the robot handle on the impaired side leads to the kinematic performance of impaired limb X I , as given in Equation 3. X I = g k I (3) where g * is an unknown function, the input setting k I are the damping coefficient of the impaired limb. The control target is to obtain the optimal input settings k * for which X * = g k * is closest to the kinematic performance of healthy limb X H . In other words, the goal is to find a proper k * to minimize the objective function in Equation 1. In practical, we use incremental form: ∆ X I = g ∆ k I .

The control strategy and conceptual ideas are as follows: the control pipeline is composed of the objective function reflecting the control objective, and the feedback section generates input adjustment ∆ k I . The objective function gives a scalar value reflecting the dissimilarity of kinematic performance of bimanual limbs. Therefore, the function provides a direction of improvement for the impaired feature vector, which is mapped to an input settings adjustment via the feedback section.

We define the control strategy for input adjustment ∆ k I as in Equation 4. ∆ k I = g − 1 ∆ X I − k I (4) where g − 1 : R n →   R 1 is established by the broad learning system method (the inverse mapping of g * in (3)). ∆ k I denotes the change of damping coefficient. In addition, ∆ X I refers to expected kinematic performance changes, denoted by the negative derivative of Equation 2: − ∇ L 1 X I = X H − X I . In other words, the negative gradient of the target equation provides an improved direction for the feature vector X I . There the control law has the form in Equation 5. ∆ k I = f X H − X I − k I (5) where f : R n →   R 1 needs to be learned, play a role in making damping adjustments based on the state performance.

The BLS, which is developed from a random vector functional-link neural network (Pao et al., 1994), provides a fast and efficient, and incremental learning way to learn and remodel in broad expansion without retraining the whole model.

In this learning system, the original input vectors are first mapped into random features through feature nodes by some feature mappings. Then, these features are sent into the “enhancement nodes” for nonlinear transformation to generate enhancement features. Finally, both the mapped features and enhancement features are fully connected to the output layer, and the full connection weights are the parameters that can be calculated rapidly by ridge regression of the pseudoinverse or can be optimized by gradient descent.

Given the training datasets X , Y and n feature mappings ϕ i , then the i t h feature matrix can be described as in Equation 6. Z i = ϕ i X W e i + ϑ e i , i = 1 , 2 , … , n (6) where X ∈ R N × M denotes the kinematic performance features (details are shown in Table 1 in Experiments Configuration Section). N is the number of input data, M is the number of features each sample, C is the dimension of outputs, weights W e i and bias term ϑ e i are randomly initialized and need to learn. ϕ i is the linear transformation function. In this work, the bias term is set to 0.

We denote Z n ≜ Z 1 , Z 2 , … , Z n as the sets of n groups of mapped feature nodes. The sparse autoencoder is utilized to generate the sparse features Z n of input. Then, Z n is fully connected to the enhancement nodes for nonlinear activation. Similarly, we obtain the outputs of the j t h group of enhancement nodes by Equation 7. H j = ξ j Z n W h j + ϑ h j , i = 1 , 2 , … , m (7) where ξ j = tanh   x is the nonlinear activation function, weight W h j and bias ϑ h j are randomly initialized. We denote the output matrix of the enhancement layer by H m ≜ H 1 , H 2 , … , H m .

Therefore, the output of a BLS, the damping coefficient adjustment Y ^ , can be denoted as in Equation 8. Y ^ = Z 1 , Z 2 , … , Z n , H 1 , H 2 , … , H m = W m Z n , H m W m (8) where W m are the weights connecting the layer of mapped feature nodes and the layer of enhancement nodes to the output layer, and it can be easily obtained by the ridge regression approximation of pseudoinverse Z n , H m + as in Equation 9. W m ≜ Z n , H m + Y (9)

In the rehabilitation scenarios, the new training data is generated from the system after each training trial and we need to build a model that is easily adaptive to the new data or different subjects’ data. Retraining the model again with the whole training data is time-consuming. However, only by updating some weights can make the BLS be adapted to the new data, which is very practical.

Given that X a , Y a refers to the new training data for a BLS, where X a and Y a separately denote new input data and corresponding target output. The mapped feature nodes for X a can be written as in Equation 10. Z a n = ϕ X W e 1 + ϑ e 1 , … , ϕ X W e n + ϑ e n (10) where Z a n denotes the feature nodes derived from the input data X a using a nonlinear mapping function ϕ . Then, the model’s output after processing the mapped feature nodes through the enhancement layer is denoted as in Equation 11. A x ≜ Z a n , ξ Z x n W h 1 + ϑ h 1 , … , ξ Z x n W h m + ϑ h m (11) where Y a T denotes the true output, and A x T W m represents the predicted output. B is an adjustment factor. Final, the incremental weights are denoted as in Equations 12, 13. W a m = W m + Y a T − A x T W m B (12) where B T = C + , i f   C ≠ 0 1 + D T D − 1 A m + D , i f   C = 0   C = A x T − D T A m   D T = A x T A m + (13)

The modeling ideas and methods have been given before, and now we will begin to establish the control strategy. There are a series of data collection and processing tasks before training the model.

Generally speaking, there is a positive correlation between the damping coefficient and kinematic performance. Therefore, the dissimilarity of the kinematic performance of bilateral limbs has a trend as shown in Figure 3B with the impaired damping coefficient increasing from min to max, and there is an optimal k * (blue point) where the dissimilarity is minimum. The optimal k * converges to a certain value greater than the damping coefficient of the healthy side.

In the pre-train experiment, the kinematic performance of the impaired side changes with the increasing k I as shown in Figure 3C. The blue horizontal line represents the approximate kinematic performance of the healthy side under a fixed damping coefficient and the orange upward trend curve represents the kinematic performance of the impaired under the uniform increase of the damping coefficient k I . The intersection (blue point) of the two lines is the point with the smallest bilateral dissimilarity, which corresponds to the optimal k * . X 0 I ( k 0 I ) represents the performance X 0 I under the damping coefficient k 0 I . The input data is the performance difference X i H − X i I and the label is k * − k i I .

Algorithm 1.

Train and Application of Feedback model.

Remark 1: The selection of hyperparameters in the BLS model can use the grid search method, and at the same time, select a smaller number of neurons to prevent the model from overfitting when the accuracy is high.

Obtained the BLS model, we can achieve a whole closed loop control pipeline. We feed the input vector into the BLS model to get the damping adjustment, which is sent to the robot device in real-time via the UDP protocol after fine-tuning.

In actual training work, the damping adjustment predicted by the BLS model mainly considers how to quickly make the bilateral limbs achieve the same kinematic performance. Large damping adjustment will give rise to the unsmooth movement (when the damping adjustment exceeds the subject’s minimum sensitivity value, subjects will feel the sudden changes in damping coefficient and unsmooth movement). To solve the above deficiency, three tricks in the following, are integrated and applied to the output results of the BLS model to achieve fast and smooth damping adjustment.

Trick 1: Introduce a saturation function to limit the maximum value of each damping adjustment output ∆ k I from the BLS model. The function is as in Equation 14. ∆ k = ∆ k , ∆ k <   ∆ k l i m i t + ∆ k l i m i t , ∆ k   >   + ∆ k l i m i t − ∆ k l i m i t , ∆ k   <   − ∆ k l i m i t (14) where ∆ k l i m i t is the limit of damping coefficient adjustment.

Trick 2: Perform n -time damping interpolation within the time sliding window T through logarithmic interpolation in Figure 4. The damping coefficients after interpolation are sent to the robot via the UDP package in turn. The interpolation calculation formula can be written as in Equations 15, 16. s t r i d e = e 1000 * k i + ∆ k − e 1000 * k i n u m (15) k i + j = ln e k i + j * s t r i d e \ 1000 j ∈ 1 , n u m (16) where s t r i d e is interpolation’s step size, n u m is the number of interpolations. k i is the i t h time window’s damping coefficient.

Trick 3: To increase the stability of the entire system, when the loss S is lower than the particular threshold, we reduce the damping adjustment by setting the decay factor δ ∈ 0 , 1 in Equation 17. ∆ k I = ∆ k , S ≥ T h r e s h o l d   δ ∆ k , S < T h r e s h o l d   (17)

The fixed-step model is a method of fixing each adjustment’s step value without offline training compared with the supervised learning methods.

Fixed-step model: We also utilized the fixed-step model for comparison. There are only two options for adjusting the damping coefficient: increasing and decreasing + 1 , − 1 . We set the damping coefficient adjustment with a fixed step value, which is the maximum step value at which the subjects cannot feel the damping coefficient’s sudden change. The fixed-step control model focuses on the kinematic features to vote for deciding the sign of the damping coefficient. The formula can be written as in Equations 18, 19. s i g n = ∑ i = 1 n S I G N x i I − x i H (18) ∆ k = s i g n * ∆ k ′ (19) where S I G N is a sign function; sign’s value is −1 or +1; x i I and x i H are the i t h feature of the impaired and healthy upper limbs.

In this part, we show the design of the experiment for demonstrating the effectiveness of the method proposed in Part D. The experimental setup is composed of the bilateral robot device and two subjects with different upper limb strength imitating the impaired limb and the healthy limb of the patient with partial unilateral movement disorders. The user datagram protocol (UDP) was utilized to send data from the robot control box to control the program and send the damping coefficient from the control program to the robot control box. The communication frequency is 200 Hz.

The kinematics dissimilarity, which is the optimization objective of the proposed framework, reflects the subjects’ movement. Figure 6 shows the evolution of the mean kinematics dissimilarity for six subjects under the three conditions in the 60s. As can be seen, the dissimilarity of bilateral kinematics performance without assistance remains at a high value. Under the adjustment policy, the bilateral dissimilarity gradually decreases to a low stable value. The adjustment time delay is defined as the duration from the beginning of adjustment to the state that dissimilarity reaches around the mean value. The results are given in Figures 6, 7 show that the adjustment time delay of the BLS method is smaller than the Fix-step method.

From the adjustment delay perspective, the BLS method can rapidly assist participants in completing mirror-symmetry movement on average. When patients feel that the impaired limb is uncomfortable or painful, he/she can adjust the kinematic performance of healthy limb to help the impaired limb relieve the trouble instantaneously. (Such as reducing movement speed).

The evolution of the left damping coefficients ( K I ) in 60s under three different conditions are shown in Figure 6. The results show that the supervised learning modeling method based on the participant training data is effective.

Three optimization indicators, including mean velocity error number of peak velocity error, and maximum peak velocity error, are used to analyze the effect of kinematics equivalent mirror-symmetry movement. As is seen in Figure 8, the Fixed-step model and BLS model can significantly reduce the difference in bilateral upper limbs’ kinematic performance error, especially the mean velocity error and the maximum peak velocity error. The “number of peak velocity error” did not increase significantly, indicating that the adjustment did not cause the subject’s movement to be unsmooth.

According to the optimization and comprehensive indicators’ evaluation results, it can be inferred that the proposed framework can achieve a great mirror symmetry effect without causing unsmooth movement. Further verification results are given in Section 3-B.

Meanwhile, we analyze variance on the dissimilarity of six participants under three different conditions in Figure 8. The test results show that both the Fixed-step method and BLS method can significantly reduce the difference of bilateral upper limbs’ performance error.

To verify the effect of the proposed framework, the mirror position error (mirror the real-time position of one side to the opposite side and calculate the Euclidean distance) and the difference of bilateral force are adopted. The error of bilateral position under three different conditions is shown in Figure 9. As can be seen, the position error has converged to low values under these adjustment models. It can be seen from the verification results that the Fixed-step method and BLS method can assist subjects with weak strength in achieving great mirror-symmetry effects.

The difference of bilateral force under three different conditions is shown in Figure 10. Under the three conditions, the average force levels of the weak subjects were similar. Combining the results of the two verification indicators in Figures 9, 10, the proposed method has the potential to help patients with unilateral motor dysfunction complete mirror rehabilitation training.

Patients with unilateral limb dyskinesia are accompanied by rhythmic and discrete movements (Montes et al., 2014; Leconte and Ronsse, 2016b). Thus, there is no doubt that movement smoothness is also an indispensable consideration to promote the recovery of a patient’s motor function (Krebs et al., 1998; Balasubramanian et al., 2012). There is a tradeoff between the damping coefficient adjustments and movement smoothness. To observe the influence of the damping coefficient’s adjustment on the smoothness of the participants’ bilateral upper limbs under the two different conditions, average acceleration, and the SPARC are compared in Figure 11. With the help of the damping adjustment strategy, the subjects’ average speed variability is closer to those of the strong group. The auxiliary strategy based on the BLS model can make the impaired side’s volatility of speed closer to that of the healthy side compared to the adjustment method based on the fixed-step model. From the SPARC indicator, the adjustment of damping can make the weaker subjects move more smoothly. As can be seen, the proposed framework allows the weaker subjects to move more naturally and smoothly like the stronger ones.

Although the proposed model has not been applied to clinical experiments, the experimental results show this model’s good performance in ordinary people’s simulation experiments. Our future work will focus on conducting clinical trials to evaluate the effectiveness of the proposed robotic mirror rehabilitation.