Development of an exoskeleton rehabilitation robot framework for the knee joint based on predictive as-sessment

Common gait deviations in cerebral palsy can be divided into spastic hemiplegia and spastic diplegia gait patterns. These gait deviations require the use of a large amount of data obtained by clinical gait analysis through video, kinematic, kinetic, electromyographic sensors, and plantar pressure data to assess patient gait characteristics. In this paper, we use predictive forward dynamics to investigate the effects of the biarticular hamstrings (HAMS), gastrocnemius (GAS) and biceps femoris short head (BFSH) on the knee joint in gait cycle. We attempted to explore the characteristics for reproducing the patient’s gait by modifying the parameters of model muscles. We applied mild, moderate and severe muscle weakness or contracture to the HAMS, GAS and BFSH muscle groups, respectively, and trained the model to walk at a self-selected speed, showing that with the more severe contracture, the non-swinging phase presented more severe knee hyperflexion and stronger knee torque, and the sensitivity for change is ranked by GAS>BFSH>HAMS. In swing phase, HAMS and GAS contractures aggravate the knee angle, whereas contractures of BFSH have a weak effect on knee angle. Mild HAMS muscle weakness accelerated walking speed, while moderate and severe HAMS muscle weakness hindered walking speed instead. BFSH muscle weakness is more sensitive to knee joint torque. Finally, a variable parameter impedance controller for the lower limb exoskeleton reha-bilitation robot is developed. We apply the knee joint angle and torque parameters optimized by predictive forward dynamics simulation as the expected values for the robot to achieve customized tuning of the motion trajectory for the exoskeleton rehabilitation robot and meet the different rehabilitation stages.


I. INTRODUCTION
Patients with neurological disorders caused by disease or injury, such as stroke and spinal cord injury, commonly present with muscle weakness or contractures which lead to lower extremity dyskinesia [1]. Pathological gait with irregular muscle groups directly causes insufficient force or torque in the lower limb joints [2]. Janet Taylor and Dario Farina [3] found that due to muscle-related changes, a slower preferred walking speed emerged in elderly people as they improved their walking ability. Common gait deviations in cerebral palsy can be divided into gait patterns of spastic hemiplegia (drop foot, horseshoe in different knee positions) and spastic diplegia (true horseshoe, jumping, pronounced horseshoe and squatting) [4].Traditional rehabilitation is not only physically exhausting for the therapist, but 44% of the patients who are rehabilitated by physiotherapy will have future problems [5]. T Mark Campbell and Guy Trudel [6] investigated the associations between knee flexion contracture with range of extension, function, pain, and stiffness of the contralateral knee. G.E. Rose et al. [7] found no significant relationship between time and body mass index for worsening knee flexion range, knee flexion in intermediate stance, peak knee extension in stance and hamstring length, and improvement in mean and maximum hip rotation in children with bilateral cerebral palsy and no history of orthopedic surgery in their gait analysis. By using robotic-assisted rehabilitation devices, it is possible to reduce the intensive work as well as to facilitate the customization of rehabilitation based on the diagnosis of the data obtained [8].
To obtain the reference trajectory, Bin Chen et al. [9] applied a motion capture system to acquisition the 3D kinematic data from the lower body for normal walking of healthy people. Moreover, the joint angles at special timing were derived from the human exoskeleton system (HES) leg geometry constraints. They also modified the designed reference trajectories, such as gait period and amplitude, for different wearers according to their physical characteristics. However, it is a challenge to determine these characteristics, not to mention determining the gait reference trajectory for patients with pathological gait. It was found that inertial or optical motion capture sensors lack accuracy when estimating joint angles during motion capture, which could lead to incorrect data interpretation. Therefore, Halim Tannonus et al. [10] proposed fusion between inertial and visual motion capture sensors to improve the estimation accuracy of joint angles. In another study, Yi Zheng et al. used a human motion capture system: four high-speed cameras, a testbed, and a computer system to study the stability of the exoskeleton robot's climbing posture using an effective data acquisition area in the view of two high-speed cameras [11].
The biartical HAMS and the monoarticular GAS and BFSH muscle groups play the most significant role in knee motion during the walking gait cycle [1][12] [13]. Weakness and contractures of these muscle groups usually occur in conditions such as cerebral palsy, stroke and secondary dysfunction after knee osteoarthritis surgery [14][15] [16]. Carmichael F. Ong et al. [1] applied the predictive dynamics simulation method to investigate how SOL and GAS muscles affect the plantarflexor of the ankle joint and generalized the relationship between them. Kirsten Veerkamp et al. [17] evaluated how gastrocnemius hyperreflexia affects gait kinematics by using predictive simulations.
In this paper, we model pathological gait by predicting forward dynamics, using minimizing the total cost of transport within a self-selected speed while ensuring head stability as a high-level goal. The gait controller utilizes a combination of state machines and low-level control laws to determine the excitation, calculate the optimal motion trajectory to perform a given task. We attempted to explore the characteristics for reproducing the patient's gait just by modifying the parameters of model muscles. Therefore, mild, moderate and severe muscle weakness or contractures were applied to the HAMS, GAS and BFSH muscle groups, respectively, and the models were trained to walk at self-selected speeds. Finally, we developed an impedance control model for the lower limb exoskeleton rehabilitation robot: we adopted the knee joint angle and torque parameters optimized by using predictive forward dynamics simulation as the expected values for the robot in order to achieve customized tuning for the robot motion trajectory. The framework not only realizes progressive rehabilitation training but also reduces the wearer's resistance.

II. METHODS
Our control framework consists of a predictive forward simulation framework and a lower limb exoskeleton rehabilitation robot based on a variable impedance controller (see Fig. 1). The neural model used in the predictive forward simulation platform is a closed-loop controller. The excitation is based on the following information from the sensor feedback: Joint position/velocity with PD control, the muscle spindles model adopts the Hill-type [16] [18], and the neural delay. Our musculoskeletal model uses Carmichael F. Ong's [1] 18 Hill-type muscle-tendon units (MTU), nine for each leg. We implement the musculoskeletal model in OpenSim 3.3 (NCSRR, USA) and the predictive forward simulation in Simulated Controller OptimizatioN environment (SCONE, https://scone.software). The kinematics and kinetics modeling for lower limb exoskeleton rehabilitation robot, and the design of impedance controller are implemented using MATLAB (MathWorks, Natick, Massachusetts, USA).

A. Overview of the framework
The predictive forward simulation framework has three main data files: Model.osim, InitialFile.par and Measures.scone (see Fig. 1). The model files are the muscle-tendon geometry and the muscle-tendon actuated model, and the initialization files defined the initial values and thresholds of parameters. For instance, by changing the pelvis velocity along the x-axis direction, it is possible to have the model walk at different speeds. The lower limb exoskeleton rehabilitation robot based on variable parameter impedance controller is mainly composed of robot kinematics and dynamics modeling, impedance controller design and parameter tuning. The yellow blocks (see Fig. 1) indicate the stored files of the model. The blue lines represent the need for human intervention and manipulation. In this paper, one of the main tasks is to establish the modification rules for the knee joint in gait cycle. The predictive simulation of knee joint motion for different pathological gait patterns can be thus reproduced, so that customized pattern adjustment of the exoskeleton robot is achieved. The physiotherapist increases the bias value based on this simulation pattern, allowing the rehabilitation robot to gradually train the affected limb and avoid secondary damage to the patient.

B. Modeling
The predictive forward dynamics simulation framework relie s on a musculoskeletal model. We use OpenSim 3.3 to imple ment the development and parameterization for the musculos keletal model. To model the human lower limb motion more realistically, our exoskeleton robot designed the center of gra vity of the rod in [19] at the center of the rod instead of at the joint (see Fig. 4).

1) MUSCULOSKELETAL MODEL
The musculoskeletal model used for the predictive dynamic simulation is based on an adult, 1.8m height and 75.2kg weight, used to simulate lower limb gait. The musculoskeletal geometry uses a trunk and two three-segment legs (see Fig. 2) to reduiyu present the human body which is a muscle reflex model proposed by Hartmut Geyer and Hugh Herr [20].
Considering the peak isometric forces, the muscle groups with similar functions in the lower limbs were combined into one MTU, so we obtained 9 MTU representing each leg: gluteus maximus (GMAX), biarticular hamstrings (HAMS), iliopsoas (ILPSO), rectus femoris (RF), vastus (VAS), biceps femoris short head (BFSH), gastrocnemius (GAS), soleus (SOL), and tibialis anterior (TA) [21][22] [23]. The tendon slack length for each MTU was calculated using experimental data [23]. Previous studies in [24] [25] found that representing muscle paths as a single line tends to overestimate length changes, so we set the maximum muscle fiber contraction velocity to 15 optimal fiber lengths per second ( 0 / ). The tendons are modeled as nonlinear springs that generate torque when the joint is hyperflexed or hyperextended. Ligaments generate torque when the hip is flexed over 120° or extended over 30°, the knee is flexed over 140° or extended over 0°, the ankle is dorsiflexed over 20° or plantarflexed over 40° [1].

FIGURE 2. Predictive forward framework for dynamic optimization. Relationship between gait controller, muscle model, contact model, measurement feedback and forward simulation. We implement musculoskeletal model in OpenSim 3.3, and use SCONE for forward prediction simulation.
We apply the Hill-type model [15] with three elements in each MTU model: contractile element (CE), parallel-elastic element (PE), and serial-elastic element (SE) (see Fig. 3). The following relations hold true for our model: where , and are complete musculotendinous actua tor, tendon force and muscle force, respectively. , , are the complete length of the model unit, muscle length and tendon slack length.
denotes the tendon stiffness propertie s. Muscle excitation, , represents the neural signals from the central nervous system and is a value between 0 and 1, repre senting the discharge rate of neurons. The excitation-activatio n model is represented by the first-order delay of equation (4). The model used the time delay ( ) parameters of [1].
For all positive feedback and PD control laws, the paramet er for time delay, is set for each muscle depending on the mos t proximal joint over which the muscle crosses based on valu es from previous simulation work using reflex-based controls [20][26]: time delay is 5 for the hip, 10 for the knee, a nd 20 for the ankle. The time delays of 20 , 10 , and 5 represent long, medium and short neural signal delays. They are not tuned but estimated from the time gaps between M-wave and H-wave of H-reflex experiments (for details see [27]). When considered with a first-order activation time cons tant of 10 , this delay better represents the short-latency TA suppression (56 to 74 ) observed in experiments [28]. Mus cle activations are calculated from the muscle excitations usin g a first-order dynamic model, with activation and deactivatio n time constants of 10 and 40 , respectively [29]. Muscl e activation, , and muscle excitation are represented by the d ifferential equation (4), where ℎ is the step size (1 2400 ⁄ ), and denote the muscle activation and excitation values at the -th timestep.
( ) is passive-element force-length curve [30] for inputs that are not normalized. The relationship between muscle force-muscle fiber length, ( ) , and mus cle force-muscle fiber length velocity, (̇), is based on th e relevant conclusions of [30] [31].
The architecture of the adopted Hill-type muscle-tendon actuator model. PE is the passive elastic element, CE is the muscle contractile element, and SE is the serial-elastic element that represents the tendon. is the complete musculotendinous actuator. , , are the complete length of the model unit and muscle length.

2) KINEMATIC AND KINETIC MODELING FOR LOWER LIMB EXOSKELETON REHABILITATION ROBOT
The lower limb exoskeleton rehabilitation robot system is a c omplex nonlinear highly coupled dynamic system [32]. The m odel of the robot has three parts: robot kinematics model, dyn amics model and variable parameter impedance controller.  Planar mechanical structure with two rotating joints (see Fi g. 4). The notation is solved by: for = 1,2, denotes the jo int turning angle, which is also used as a generalized coordina te, denotes the mass of link , denotes the length of link , denotes the distance between the previous joint and the c enter of mass of link , and denotes the rotational inertia of link about the axis that passes through its center of mass and points out of paper. We use the Denavit-Hartenberg joint variables as generaliz ed coordinates and thus are able to efficiently derive the Jaco bian matrix expressions to calculate the kinetic energy as with: The inertia matrix is equation (8) combined with the trigonometric constant equati on (9) is obtained as Calculating the Christoffel symbols, we get ] . (11) The gravity moment vector ( ) is given by ] . (12) The dynamic equation of the lower limb exoskeleton rehab ilitation robot is VOLUME XX, 2017 1 Where ∈ ℝ 2 is the angular displacement of the joint, ( ) ∈ ℝ 2×2 is the inertia matrix of the robot, ( ,̇) ∈ ℝ 2 denotes the Centrifugal and Coriolis forces, ( ) ∈ ℝ 2 is the gravity term, ∈ ℝ 2 is the control torque and ∈ ℝ 2 is the applied perturbation. Simulation parameters for the robot (see Table 1).

1) GAIT CONTROLLER FOR PREDICTIVE FORWARD SIMULATION FRAMEWORK
We trained a planar musculoskeletal model actuated by 18 MTU actuators to walk by optimizing the parameters of a gait controller. The parameters are evaluated using a combination of high-level state machines and low-level control laws to calculate muscle excitation. The high-level state machine has two states in stance: early stance (ES), mild stance (MS), and three states in swing: pre-swing (PS), swing (S), and landing preparation (LP). Five transitions are generated between the five states, where the transitions associated with landing and standing are determined by comparing the ground reaction force of the ipsilateral foot to the threshold. In contrast, the transitions related to swing is determined by comparing the horizontal distance of the ipsilateral foot from the pelvis to the threshold. The five transitions are ES to MS, where the horizontal distance between the ipsilateral foot and the pelvis is less than the threshold; PS to S, where the ground reaction force of the ipsilateral foot is below the threshold; S to LP, where the horizontal distance between the ipsilateral foot and the pelvis is greater than the threshold; LP to ES, the ground reaction force of the ipsilateral foot is greater than the threshold; MS to PS, which is not controlled by the free parameter, occurred when the contralateral leg entered the ES state. The transition for each state is determined by the activation of the low-level control law. The low-level control laws include signal constants, feedback based on muscle length, muscle velocity and muscle force, and PD control based on pelvic tilt angle. The positive and negative feedback is denoted by (+) and (-), respectively (see Fig. 5). In the low-level control laws, given the MTU model, the positive force feedback law, positive length feedback law, positive length velocity feedback law and Muscle-driven PD control law are defined as Where ̃( − ) , ̃( − ) and ̇( − ) are the model force normalized by , and ̇ with a time-dela y of , respectively. , , , are proportional coefficients, differential coefficients and joint angles of the PD controller, respectively. PD controller is to ensure the stability of the con trolled joint movement. The objective function, , quantified high-level tasks of walking: The goal is to minimize the gross cost of transport ( ) withi n the specified speed ( ), while ensuring head ( ℎ ) stabi lity. To balance the competitive objectives, the weights w ere manually adjusted to the following values[1]: = 1 / , = 10000 −1 , ℎ = 0.25 3 / 2 These weights determine the priority of the solution, that is, the contribution of the ℎ term is greater than the contribu tion of , while ensuring that is minimal The forward dynamics model for gait controller consists of 70 free optimization parameters, 16 joint offset, four range th resholds and load thresholds for swing and stance, respectivel y [33]. The gait controller is implemented in SCONE, includi ng leg states update, target features update and computation o f the excitation signal (simulating the central nervous system), with output as ( − ). The gait controller computes musc le excitations, ( ), for a musculoskeletal model to generate a forward simulation Sensory feedback, based on the model' s muscle and joint states, is used in a feedback loop with the g ait controller. The objective function quantifies the performan ce of each simulation, and a Covariance Matrix Adaptation E volutionary Strategy (CMA-ES) optimization method update s the values of the variables in the optimization (see Fig. 2). VOLUME XX, 2017 9

2) VARIABLE PARAMETER IMPEDANCE CONTROLLER FOR ROBOT
The lower limb exoskeleton rehabilitation robot system is a complex nonlinear highly coupled dynamic system [32]. The wearable lower extremity rehabilitation exoskeleton is a human-machine interactive robotic system, thus requiring not only the wearer to be comfortable, but the physiotherapist can adjust the robot's damping force and stiffness force according to the patient's rehabilitation level. The impedance control relationship is where, represents the contact force with the environment, = − . , , and denotes the inertia matrix of the impe dance model, damping matrix and stiffness matrix of the robo t, respectively [34]. denotes the error between actual positio n and desire position ( = − ). Considering that the reha bilitation robot we developed is mainly used in the passive tra ining stage of the primary training for the affected limb. Owi ng to the relatively gentle and comparatively small joint angu lar accelerations during this training stage, the inertia matrix ( ) of the robot kinematics and the inertia matrix of the i mpedance model are considered equal. Thus, when equation (13) is applied to the external environmental forces to obtain weakly coupled impedance control in joint space, as shown in the following nonlinear equation: where, ∈ ℝ 2 is the joint angular displacement of the joint. ∈ ℝ 2 is the desired value of joint angular displacement an d. ( ) ∈ ℝ 2 is the gravity term.
∈ ℝ 2 is the desired dam ping matrix.
∈ ℝ 2 is the desired stiffness matrix. and are generated by the variable impedance controller (see Fi g. 6). Equation (20) is a simplified spring-damped first-order dynamic system via a second-order system consisting of sprin g-damped-mass. This processing eliminates the torque feedba ck term and holds the possibility for subsequent simplificatio n of the robot hardware. The red dotted box (see Fig. 6) shows the exoskeleton robot ankle and knee joints both applying the impedance control of equation (20). The variable impedance controller in the black dotted box (see Fig. 6) allows the physiotherapist to adjust the impedance parameters and according to the rehabilitation training needs. For example, the physiotherapist can increase the damping coefficient when active training of the affected limb is required, and similarly, the elasticity coefficient can be increased when the robot is required to drive the affected limb for passive training. The rules for determining the impedance parameters will be investigated in our next work.

1) VALIDATING THE MODEL'S GAIT OVER A RANGE OF SPEEDS
We validated the ability of the predictive forward model and applied the optimization framework to capture the walking tr ends at four different speeds: three prescribed speeds 0.6 m/s, 1 m/s and 1.4 m/s, and a self-selected speed. Individual comp arisons of each speed with the experimental data of Schwartz et al [36] are (see Fig. 7).
The simulated kinematic and kinetic adaptations, joint angl es matched the trends observed in the experimental data [36] [37]. The optimized gait generated by walking at a constant s peed of 0.6 / had significantly larger hip (see Fig. 7a) and knee angles (see Fig. 7b) in the late swing phase. At the pres cribed speed of 1 / , the hip and knee angles were essential ly the same as the self-selected speed in gait cycle, but in the pre-swing phase, the knee angles (see Fig. 7b) were smaller c ompared to the experimental data. With the prescribed speed of 1.4 / , the hip joint angle (see Fig. 7a) was out of range i n the pre-swing phase and the knee joint (see Fig. 7b) entered flexion earlier, indicating that the optimized gait speed was i ncreased.
Overall, the generated optimized gait for walking at a constant speed of 0.6 / exhibited greater angles at the hip, knee, and ankle joints than the other three speed patterns in gait cycle. The self-selected speed was the optimized gait generated by setting the initial speed in the range of (0.5-1.5 / ). The trend for the self-selected gait was most similar to the prescribed speed of 1 / . In addition, the ankle joint (see Fig. 7c) of our model showed a significant early entry into dorsiflexion compared to the experimental data. However, this paper mainly focuses on finding the impacts of HAMS, GAS and BFSH on the knee joint. Therefore, it has no influence on the subsequent study.
This indicates that the magnitude of the prescribed speed is not necessarily related to the change of each joint angle in the optimized gait. The self-selected speed could find a solution to the optimization framework and be insensitive to the initial guesses. So, we subsequently used self-selected speed to validate the relationship between HAMS, GAS and BFSH and knee joint.    The previous results gave us confidence in the forward predictive model. We then used the framework to study how the HAMS, GAS and BFSH muscle groups act on the knee joint in gait cycles under contracture or muscle weakness. We modified the optimal fiber lengths of the HAMS, GAS, and BFSH muscles in the musculoskeletal model to 85%, 70%, and 55% of their original value, respectively, represented by mild, moderate, and severe degrees of contracture. Similarly, we modified the maximum isometric forces of these three muscle groups to 25%, 12.5%, and 6.25% of their original value, respectively, and expressed them as mild, moderate, and severe muscle weakness, respectively. Angle simulation with muscle contracture demonstrates that when different degrees of contracture are applied to HAMS (see Fig. 8a), GAS (see Fig. 8b) and BFSH (see Fig. 8c), the knee joint exhibits hyperflexion during the contact phase with the ground, and the degree of sensitivity from strong to weak is GAS>BFSH>HAMS. Contractures of HAMS (see Fig. 8a) and GAS (see Fig. 8b) results in knee hyperflexion during the swing phase, but it is noted that when the HAMS contracture is severe (see Fig. 8a), the knee flexion during the swing phase is instead significantly relieved. The more severe GAS contracture (see Fig. 8b), the more severe knee flexion. The moment simulation results indicate (see Figs. 8(d)-(f)) that all three muscle contractures cause the knee torque to increase, and as the contracture becomes more severe, the knee torque becomes greater. Muscle weakness in the GAS (see Fig. 9b) and BFSH (see Fig. 9c) muscle groups had little effect on knee flexion, but mild muscle weakness in the HAMS (see Fig. 9a) had an accelerating effect on walking speed, while moderate and severe HAMS muscle weakness hindered walking speed instead. The effect of muscle weakness in BFSH is more sensitive to the knee moment (see Fig. 9f ): the more severe the muscle weakness of BFSH, the weaker the knee moment in the whole gait cycle.
On the basis of these results, we assessed the level of influence by HAMS, GAS and BFSH muscle groups on knee joints in gait cycle (see Table 2). (+) represents the enhancement effect, (-) represents obstruction. The number indicates the intensity. The larger the number, the heavier the impact. The assessment reference value reflects the extent of impact by HAMS, GAS and BFSH muscle groups on knee joint in gait cycle.
The physical therapist is able to modify the MTU paramete rs in the predictive forward simulation framework referencing Table 2. The contracture effect of the MTU actuators is obtai ned by modifying the optimal fiber length of the muscle mod el, similarly, the muscle weakness effect of the MTU actuator s is obtained by modifying the maximum isometric force of t he muscle model. It is experimentally demonstrated that each change in MTU parameters has a different effect on the opti mization of knee gait. Our aim is to reproduce the gait simula tion of the affected limb using the forward predictive simulati on framework in order to obtain real-time data on knee joint a ngle and knee joint moment in gait cycle. The parameters mo dified by the physiotherapist are used as the desired values fo r the lower limb exoskeleton rehabilitation robot controller to realize the robot's spatial motion to track the movement of th e affected limb offline. In addition, physiotherapists can adjus t the damping and stiffness parameters offline through the pat ient's condition, so that the robot can not only quickly follow the affected limb trajectory, but also adjust the damping force and stiffness force to realize active rehabilitation training or passive rehabilitation training.

2) EXPERIMENTAL VALIDATION OF MUSCLE ACTIVA-TION
In the gait cycle, muscle actuators are directly reflected by m uscle activation on its surface [26]. Therefore, we collected e xperimental data on muscle activation in the HAMS, GAS an d BFSH muscle groups of the normal human right leg in gait cycle. The height and weight of the testee are approximately 175 and 80 , respectively. Electromyograms (EMG) se nsor use DTing XS01 type dry electrode biosensor (see Fig. 1 0) and the sampling frequency is 1000 , the signal-to-nois e ratio is 59. 5 , the sampling accuracy is 12 , the amplif ication is 700 times. DTing XS01 is easy to measure but not e asily fixed to the skin. Therefore, we fix the sensor with nylo n straps and tape to avoid relative displacement to the skin. T he EMG sensor communicates with the handheld device via Bluetooth in real-time to record and generate image informati on.
Considering that the intensity of muscle activation increases with exercise time, the surface EMG signal is enhanced. Therefore, we intercepted the EMG signals of these three muscle groups in the latter part of the experiment for analysis, for it was found that the cycle changes of EMG signals were relatively stable in the latter part of the experiment (see Fig.  10). We then used the muscle on-off time to indicate the muscle activation, and compared the binarized processed EMG signals with simulated muscle activation (normal gait, self-selected speed) (see Fig. 11).
The simulated muscle activity presented many of the significant features observed in the experiments. During early stance and landing preparation, it was the increase in body weight that made the HAMS very active. In mid-stance, when the knee joint needs to be driven to start flexion, the GAS behaves quite actively. While the swinging phase, the BFSH was active, allowing maximum knee flexion. There were some differences between our simulated data and experimental EMG, in which the activity of HAMS and BFSH during landing preparation was lower than the experimental data. The activation sequence of the three muscle groups also ensured a continuous flexion-extension of the knee joint in gait cycle.

3) SIMULATION OF ROBOT KNEE JOINT TRACKING
The approximate clinical gait knee flexion angle and torque parameters obtained using the predictive forward simulation framework are used as inputs to the lower limb exoskeleton rehabilitation robot controller with variable impedance parameters. The input values and impedance parameters allowed for adjustment by the physiotherapist depending on the rehabilitation situation. We performed simulations for the variable parameter impedance controller with gravity compensation in MATLAB (MathWorks, Natick, Massachusetts, USA). Three sets of impedance parameters ( = 10 / , = 5000 / , = 50 / , = 500 / and = 50 / , = 5000 / ) were set and external sinusoidal perturbations disturbance was applied during 35%-40% of the gait cycle. The robot knee joint angle trajectory tracking and torque tracking in gait cycle were obtained separately (see Fig. 12).
The simulation results show that the knee joint angle trajectory tracking and torque tracking of our established robot controller model are effective. When external perturbations are applied, both robot knee joint angle tracking and torque tracking vary greatly with different impedance parameters. Specifically, the large stiffness of the impedance controller makes the robot knee joint angle tracking excellent in realtime, and the large damping of the impedance controller renders the robot knee joint show more obvious elasticity when it receives external resistance.
The knee moment for the optimized gait is consistent with the overall trend of the experimental data, however, between 35% and 45% of the gait cycle (from the right foot leaving the ground to the start of swing). It might be caused by the error between the model center of gravity and the normal human model center of gravity position during the motion. During 80% to 100% of the gait cycle (from right heel touching the ground to left toe off the ground), the simulation deviated from the experiment, and the model trunk tilted forward prematurely. The knee joint angle in the optimized gait of a normal human is smaller than the experimental data in the stance phase. The optimized model exhibited leg stiffness in the stance phase, probably owing to the inadequate design of the model muscle actuators, which could not mimic completely the effect of human lower limb muscle groups on the knee joint.

Ⅳ. CONCLUSIONS
In this paper, a lower limb exoskeleton rehabilitation robot fr amework based on predictive forward simulation was propos ed, which mainly consists of a clinical gait assessment frame work based on predictive simulation and a variable parameter impedance robot control framework including gravity compe nsation. The paper established a musculoskeletal model for pr edictive forward simulation and opted for self-selected speed to investigate the impact relationships of the model's HAMS, GAS and BFSH muscle groups on the knee joint in gait cycle, respectively. On this basis, it was proposed that the paramete rs of the predictive simulation muscle model were appropriat ely modified according to the affected limb diagnosis combin ed with the impact relationships of the model muscle groups on the knee joint to obtain the simulation model and the optim ized parameters of each joint to obtain the maximum matchin g of the clinical gait. Real-time parameters of knee joint angle and torque in gait cycle were applied as the expected values of impedance controller for exoskeleton rehabilitation robot i n this paper, and different angle trajectory tracking and torque tracking effects of the robot were obtained by adjusting dam ping matrix and stiffness matrix parameters of impedance con troller.
We developed an impedance controller with variable parameters for the lower limb exoskeleton rehabilitation robot. Although gravity compensation is taken into account, we equate the inertia matrix of the impedance model and the inertia matrix of the robot model, however, in the later stage of rehabilitation training, the acceleration of the rehabilitation robot will increase with the strength enhancement of the affected limb, and the performance of our impedance controller will decline. To achieve the ideal impedance control, it is necessary to meet the requirements of obtaining real-time feedback of the robot joint angle and joint angular velocity, real-time feedback of the robot joint torque, and torque control of the robot. Therefore, our impedance controller is only suitable for the early stage of rehabilitation when the rehabilitation movement is slow. Actually, for the actual robots, the selection of impedance parameters varies greatly due to their different kinetic parameters. Our team is trying to solve this problem. First, we intend to obtain the mechanical impedance parameters based on online recognition. Then we apply fuzzy neural networks to dynamically tune the target impedance control parameters to adapt to the changes of the affected limb. In addition, we also plan to design robust adaptive PD controllers to achieve smooth trajectory tracking for passive training modes in which the upper bound of the perturbed signal is known. To realize the ideal impedance control, we need to satisfy the requirements of obtaining realtime feedback of the robot's joint angle and joint angular velocity, real-time feedback of the robot's joint torque, and the robot's torque control. These are also the main contents of our next work.