A New Projected Active Set Conjugate Gradient Approach for Taylor-Type Model Predictive Control: Application to Lower Limb Rehabilitation Robots With Passive and Active Rehabilitation

In this paper, a three-order Taylor-type numerical differentiation formula is firstly utilized to linearize and discretize constrained conditions of model predictive control (MPC), which can be generalized from lower limb rehabilitation robots. Meanwhile, a new numerical approach that projected an active set conjugate gradient approach is proposed, analyzed, and investigated to solve MPC. This numerical approach not only incorporates both the active set and conjugate gradient approach but also utilizes a projective operator, which can guarantee that the equality constraints are always satisfied. Furthermore, rigorous proof of feasibility and global convergence also shows that the proposed approach can effectively solve MPC with equality and bound constraints. Finally, an echo state network (ESN) is established in simulations to realize intention recognition for human–machine interactive control and active rehabilitation training of lower-limb rehabilitation robots; simulation results are also reported and analyzed to substantiate that ESN can accurately identify motion intention, and the projected active set conjugate gradient approach is feasible and effective for lower-limb rehabilitation robot of MPC with passive and active rehabilitation training. This approach also ensures computational when disturbed by uncertainties in system.


INTRODUCTION
The number of limb impairment patients who were injured by stroke has increased year by year, and this disease has also been developing in the direction of youth, seriously endangering the health of patients (Zorowitz et al., 2013). Recently, researchers have paid a great deal of attention to robotics to promote the development of scientific and engineering fields (Jin et al., 2017;Jin and Li, 2018;Xie et al., 2019Xie et al., , 2020Zhang et al., 2020). Compared to traditional rehabilitation training methods that see problems like resource consumption, high costs, and long rehabilitation period, lower-limb rehabilitation robots can be deemed as a more effective method for recovering patients' movement function. In virtue of interaction generally existing between the lower-limb rehabilitation robot and the patient, to avoid the a second injury in the patient during rehabilitation training, it is essential to propose a humanmachine interactive control method, which can be utilized to investigate and analyze the lower-limb rehabilitation robot (Fleischer and Hommel, 2008).
Intention recognition is one of the key points for realizing human-machine interactive control methods with lower-limb rehabilitation robots. Generally, motion intentions include joint angles and angular velocities, which can be recognized by decoding bioelectrical signals. Afterwards, the intentions are referenced by the lower-limb rehabilitation robot to complete interaction (Ding et al., 2016;Peng et al., 2018). An appropriate alternative is to establish a relationship between biological signals and movements for the patient. The surface electromyography (sEMG) can be regarded as the biological signals, which is a micro-electrical signal that appears 20-80 ms before the muscle contraction (Fleischer and Hommel, 2008). Involving two approaches to construct the relationship, one of which is a physiological muscle model, such as the Hill muscle model and the Hammerstein muscle model (Hunt et al., 1998;Buchanan et al., 2004), muscle forces and joint motions can be estimated by those models from sEMG signals; meanwhile, the unidentified physiological parameters of those models affect their applications in rehabilitation systems (Han et al., 2015). Another is the regression model, which can be established to connect sEMG signals to indicate intentions in a straightforward manner, rather than considering physiological parameters (Ding et al., 2016). For instance, a BP neural network was exploited to describe the relationship between sEMG signals and motion intentions, which was verified on able-bodied subjects and patients (Zhang et al., 2012); least squares support vector regression was proposed to predict periodic lower-limb motions from multi-channel sEMG signals . Those approaches constitute a connection between human's bioelectrical signals and motion intentions, and intention recognition can thus be realized effectively.
During the rehabilitation training, the motion trajectory, which is a predetermined curve or recognized by sEMG signals, is known to the rehabilitation device, and the patient is assisted by the lower-limb rehabilitation robot to recover. However, it is very important to avoid the risk of second injury for patients in rehabilitation, and a human-machine interactive control method should therefore be considered to increase security and stability of rehabilitation robots. Recently, some classical control methods have been developed and applied to rehabilitation robots, for example, the rehabilitative system was realized by an adaptive control framework and a human-machine interactive method; meanwhile, the potential conflicts between patient and rehabilitation robot were rejected by position-dependent stiffness and predetermined trajectory (Zhang and Cheah, 2015). Pehlivan et al. (2016) presented a minimal assist-as-needed controller for rehabilitation robots, which could provide corresponding assistance for patients during rehabilitation training.
In recent years, model predictive control (MPC) not only considered the constraints of a non-linear system but predicted future states. People have therefore paid more attention to investigating it and applying it to the aerospace, automobile, economics, and robotics fields. As the patient should be assisted by the lower-limb rehabilitation robot within a relatively safe motion range and protected against accidents, the MPC method is suitable for human-machine interactive control of the lowerlimb rehabilitation robot. Generally speaking, the MPC method is designed to optimize multivariable and constrained control systems. On the one hand, a control sequence is created by minimizing an optimization objective function over a finite prediction horizon within state and control constraints. On the other hand, the first optimal solution of non-linear optimization problem feeds back to the non-linear systems, which is utilized to generate the next iteration (Mayne, 2014). A key issue for MPC is that the computational burden of real-time optimization should be reduced through neural networks (Yan and Wang, 2014). To solve this problem is to linearize non-linear systems and discretize the differential term. More and more people have consequently developed some classical methods. For example, a neural network was utilized to identify the unknown nonlinear discrete system, and then one-order Taylor expansion formula was used to linearize the MPC problem (Pan and Wang, 2012). For non-linear continuous systems, in order to discretize the differential term and guarantee the higher precision, a Taylor-type numerical differentiation formula was developed and applied to solve non-linear time-varying optimization problem , non-linear time-varying equations , time-varying matrix inversion (Guo et al., 2017), future dynamic non-linear optimization problem (Wei et al., 2019), time-dependent Sylvester equations (Qi et al., 2020), and so on. As the MPC problem can be seen as an optimization problem, a Taylor-type numerical differentiation formula is thus also suitable for solving the MPC problem online.
Another key issue of MPC to be looked at further is online optimization. In fact, an MPC problem can be converted to a non-linear optimization problem with equality constraints and bound constraints and be solved to obtain the control sequence at each sample time. There are thus numerous algorithms proposed and studied for this non-linear optimization problem. A complex-valued discrete-time neural dynamics is studied by Qi et al. (2019) for solving time-dependent complex quadratic programming (QP), which possesses high accuracy and strong robustness but is only suitable for linear constrained QP problems. A trust region-sequential quadratic programming approach attempts to solve a sequence of QP subproblems of non-linear constraint optimization problems, which is based on trust-region technology and applied by finding the Karush-Kuhn-Tucker (KKT) points. It is noteworthy for the approach that the compatibility of the QP subproblem and the Maratos effect were overcome by adding several linear equations into the traditional trust region-SQP algorithm (Sun et al., 2019). However, the calculation of this algorithm is increased by the added linear equations. In addition, some troubles maybe emerge; for instance, one is the consistency of the coefficient matrix and the other the computational burden. Similarly, conjugate gradient methods can also be regarded as an effective optimization approach that utilizes an iteration point with a steep descent direction to generate conjugate direction and compute a global minimum point instead of solving linear equations of trust region-SQP algorithm (Sun et al., 2019). Some classical conjugate gradient methods include the Hestenes-Stiefel (HS) method (Hestenes and Stiefel, 1952), the Fletcher-Reeves (FR) method (Fletcher and Reeves, 1964), the Polak-Ribiére-Polyak (PRP) method (Polak and Ribière, 1969;Polyak, 1969), the Dai-Yuan (DY) method (Dai and Yuan, 2000a), the Liu-Storey (LS) method (Liu and Storey, 1991), and the conjugate descent (CD) method (Dai and Yuan, 2000b), but those methods were exploited to solve unconstrained optimization problems off-line. The MPC problem usually contains some constraint qualifications. Some modified conjugate gradient methods, which consider the constrained conditions, were thus proposed by modifying the search direction and a projected operator (Dai, 2014;Sun et al., 2018). Besides, a projected gradient method, which projected the gradient into the feasible region, was proposed by Rosen (1960), and some modified conjugate gradient methods were extended by some researchers based on the mentioned methods (Li and Li, 2013;Dai, 2014). Those modified conjugate gradient methods were also applied in optimal robust controllers and robots (Sun et al., 2018). In this paper, a modified conjugate gradient method, which will simultaneously consider equality constraints and bound constraints, will be utilized to solve MPC problem. Furthermore, the proposed algorithm of this paper is further applied to the lower limb rehabilitation robots with passive and active rehabilitation.
There are three significant contributions to be developed in this paper. The primary one is that a new projected active set conjugate gradient approach is developed and investigated, and rigorous proof of feasibility and global convergence is also given. The second is that a relationship between sEMG signals and motion intentions established by an echo state network (ESN) model can identify human motion state. Finally, a numerical simulation about passive and active rehabilitation training of lower-limb rehabilitation robot is illustrated and solved by the proposed method. Surprisingly, the studies on the rehabilitation training of lower-limb rehabilitation robot for MPC problem with projected active set conjugate gradient approach are scarce. This motivates our present study.
The rest of this paper is organized as follows: In section 2, the MPC problem is introduced, and a three-order Taylortype numerical differentiation formula is proposed and utilized to discretize MPC model. In section 3, a new projected active set conjugate gradient algorithm is developed, analyzed, and investigated for the MPC problem, which can be generalized from non-linear systems and lower limb rehabilitation robots. Furthermore, the feasibility and the global convergence of this approach are also proven. The relationship of sEMG signals and motion intentions is established by an ESN model in section 5; meanwhile, passive and active rehabilitation training of lower-limb rehabilitation robot is illustrated and simulated by the proposed method, which is based on sEMG signals with ESN model and MPC problem. Furthermore, the disturbance of dynamic model is also considered through simulation in section 4. Finally, section 6 summarizes the results of lower-limb rehabilitation robot based on the MPC technique and expects future work.

Problem Description
In this subsection, consider the following non-linear control system: where x k ∈ R n is a system state variable, u k ∈ R m is a control input signal, A(x k ) ∈ R n×n , B(x k ) ∈ R n×m , C(x k ) ∈ R n are the state matrix, control input matrix and constant matrix, respectively. y k ∈ R n denotes the system output, and h(·) is a non-linear function. MPC for the non-linear control system is devoted to generating a sequence of control signals by minimizing an objective function repeatedly over a finite moving prediction horizon with system state and input constraints satisfied simultaneously. For the non-linear control system (1), the MPC problem can be described as a non-linear discrete-time optimal control problem within input and state constraints: (2) where r(k + i|k) and y(k + i|k) are the desired output and the predicted output for ith step ahead from kth sampling instant; u(k + j|k) = u(k + j|k) − u(k + j − 1|k) denotes the control increment; N and N u are the prediction horizon and control horizon, respectively; X k = {x(k+1|k), x(k+2|k), . . . , x(k+N|k)}, U k = {u(k|k), u(k + 1|k), . . . , u(k + N u − 1|k)}; Q ∈ R n×n and R ∈ R m×m are real positive definite matrix, · Q and · R are Euclidean norms defined as

Three-Order Taylor-Type Discretization for MPC
In this subsection, a three-order Taylor-type numerical differentiation formula with truncation error of O(h 2 ) is constructed for the first-order derivative approximation and is exploited to discretize MPC problem (Jin and Zhang, 2014). To obtain the higher-order truncation error, Guo et al. (2017) proposed a novel Taylor-type numerical differentiation formula for the time-varying matrix inversion.
, h denotes the sampling gap. Subsequently, a three-order Taylor-type numerical differentiation formula can be obtained as follows: with a truncation error of O(h 3 ).
Proof. See Guo et al. (2017). According to the Lemma 1, the non-linear control system (1) can be discreted as Therefore, the MPC problem (2) can be rewritten in the following form: where x i = x(k + i|k), u j = u(k + j|k), and y i = y(k + i|k), and another symbols see Appendix.
In general, the MPC problem can be solved by the conjugate gradient approach (Šantin and Havlena, 2011), and a sequence can be regarded as an initial control input sequence of the MPC problem at Step i. In addition, the first element u * 0 of optimal solution u * can be seen as a feedback control law for non-linear control system (1).
For simplicity, the MPC problem (5) is equivalent to the following non-linear optimization problem with linear equality constrain and bound constrain: where x = (x, u), Ŵ(·) is a continuously differentiable function, is a constant matrix, b is a constant vector, and

PROJECTED ACTIVE SET CONJUGATE GRADIENT ALGORITHM FOR NON-LINEAR CONSTRAINED OPTIMIZATION
In this section, a projected active set conjugate gradient algorithm is proposed to solve the following non-linear optimization problem: and the convergence of the proposed approach is developed, investigated, and analyzed as follows.

Projected Active Set HS-Type Conjugate Gradient Algorithm
To further analyze projected active set conjugate gradient algorithm, some basic definitions and notions should be revisited in this subsection. Let x * be a stationary point of (7), and consider the following active set: Furthermore, define as a set of free variables, where L * is the complement of H * . Therefore, the KKT conditions for the problem (7) can be converted as follows: where ∇Ŵ i (x) is the ith element of the gradient for Ŵ at x. According to the literature (Kanzow and Klug, 2006), H(x) and L(x), which approximate the active set and the free variables set can be defined as follows: is a projection function defined as P (x) = argmin ω∈ x − ω 2 . ψ 0 is a positive scalar, which should be sufficiently small. Furthermore, it should satisfy the following inequality: In what follows, let x k be the point of iteration k, and for simplicity, we abbreviate H(x k ) and L(x k ) as H k and L k .
Frontiers in Neurorobotics | www.frontiersin.org According to the literature (Cheng et al., 2014), the active set H k will be divided into the following three parts: . It is demonstrated that the active set H k 1 can be seen as the equality constraints of (7), and according to Rosen's gradient projection method (Rosen, 1960;Dai, 2014), an active set projection matrix is given as follows: where and γ 0 is a positive constant. In what follows, the search direction d k can be rigorously proved as a feasible descent direction of Ŵ at x k for non-linear optimization problem (7). Theorem 1. Suppose that x k ∈ {x| x = b, x ∈ } holds, and x k is not a stationary point of (7), d k is defined by (10), then the search direction d k is a feasible descent direction of Ŵ at x k for non-linear optimization problem (7).
Proof. According to (8) and (10), and as per definition of the projection matrix, the following two cases can be generalized.
Case 1. If k = 0 or ∃i ∈ H k 2 ∪ H k 3 , then the inequality can be directly computed as follows: Case 2. If k ≥ 1 and ∀i ∈ L k , then the following inequality can be directly obtained: Now, the proof of the feasibility for d k is shown as follows, and it is further inferred that If E k is not an empty set, it can be seen that If E k is an empty set, then M k = . Owing to Equation (10), then the following equation can be generalized as Owing to Equations (13)-(15), then we have M k d k = 0 for all k ≥ 0. It is also inferred that d k is a feasible descent direction for non-linear optimization problem (7).
2 According to the above analysis and investigation, the projected active set HS-type conjugate gradient algorithm (PASHS) is developed, analyzed, and verified for non-linear optimization problem (7).
Frontiers in Neurorobotics | www.frontiersin.org Remark 1. A projected matrix P k is computed by an active set, which ensures iteration point satisfying equality and bounded constraints of non-linear optimization problem (7). Assume that if the component the active set H k is not contained in the previous iteration point x k , the search direction d k is updated by the second formula of (10); otherwise, the search direction d k is generated by the projected gradient method, which is the first formula of (10). Furthermore, combining with Armijio-type line search, it can be proved that the proposed PASHS algorithm guarantees the feasibility and global convergence for non-linear optimization problem (7).

Convergence Analysis
In this subsection, to further investigate the convergence of the Algorithm 1 (PASHS) for non-linear optimization problem (7), some basic assumptions should be revisited and introduced in this subsection.
Assumption 1. The level set is bounded. Assumption 2. Given that the objective function Ŵ : R n → R is continuously differentiable on an open set N ⊆ D and its gradient is Lipschitz continuous, there exists a positive constant W > 0 that satisfies the following inequality: As {Ŵ(x k )} is a descending sequence, it is clear that the sequence {x k } generated by Algorithm 1 (PASHS) is contained in D. In addition, according to Assumption 1, it is inferred that the gradient of Ŵ is bounded, i.e., there exists a positive constant γ > 0 such that Since the matrix P is a projected matrix, suppose that there exists a positive constant C > 0, and the following inequality can be obtained: Lemma 2. Assume that the iterative sequence {x k } generated by Algorithm 1 (PASHS). The step size η k is obtained via the Armijo line search rule (16), and then there exists a positive constant c 0 > 0 such that the following inequality holds for sufficiently large k.
Proof. According to Armijio-type linear search rule (16), the following inequality can be obtained: Combined (11) and (12), and the properties of the projected matrix P k , the inequality can be generalized as follows Now, the following two cases can be taken into account and be utilized to prove (21). Case 1: If η k = 1, according to Equations (11), (12), and (23), and by applying the Cauchy-Schwarz inequality, we can derive that Thus, the inequality (21) holds.
Case 2: If η k < 1, assume that the Armijio-type line search rule is not true, there thus exists a positive constant ρ −1 η k such that the following inequality holds true: Using the mean-value theorem and Assumption 1, there exists a positive constant ξ k ∈ (0, 1) such that x k + ξ k ρ −1 η k d k ∈ D and (25) Combining inequality (24) and Theorem 1, the following inequality can be directly computed as Let c 0 = min{1, (1 − δ)ρ/W}, the conclusion is true.
Proof. According to the Assumption 2 and Algorithm 1 (PASHS), the following inequality can be directly obtained: In addition, in term of the search direction (10) and Theorem 1, the following inequality can be derived as (28) where λ max (P k ) is the maximum eigenvalue of the projected matrix P k .
2 Theorem 2. Suppose that Assumption 1 holds. The iterative sequence {x k } is generated by Algorithm 1 (PASHS), and then Proof. According to (21), there exists a positive constant c 0 such that Combining Algorithm 1 with Lemma 3, it implies that lim k→∞ η k d k 2 = 0.
The following inequality can be generalized as k → ∞, Hence lim k→∞ inf P k ∇Ŵ k = 0.
2 Remark. Owing to Theorem 2, it can be seen that the Algorithm 1 (PASHS) is globally convergent for the non-linear optimization problem (7). Combing the ESN learning algorithm and Algorithm 1 (PASHS), the optimal controller of the MPC problem can thus be solved rapidly; this is used for the patients through a lower-limb rehabilitation robot with passive and active rehabilitation training.

SIMULATIONS AND RESULTS
In this section, the proposed PASHS algorithm with MPC technique is applied to the passive rehabilitation training of the two-link lower-limb rehabilitation robot. Moreover, combining the ESN model and intention recognition, the MPC and PASHS algorithm also are utilized to active rehabilitation training.

Two-Link Lower-Limb Rehabilitation Robot With MPC
The general dynamic model of two-link lower-limb rehabilitation robot is shown as follows (He et al., 2015): where q,q,q ∈ R 2 are angle, angular velocity and angular acceleration of hip and knee, respectively; τ ∈ R 2 is a torque for the rehabilitation robot, which represents admissible control inputs; D(q) ∈ R 2×2 is a positive-definite inertia matrix; C(q,q) ∈ R 2×2 is a centrifugal and Coriolis term; and G(q) ∈ R 2 is related to gravity term. The state space expression of (32) can be described as where y ∈ R 2 is the end-effector position coordinates, h(·) is a function mapping angles of the rehabilitation robot to the position coordinates. The schematic of a two-link lower-limb rehabilitation robot is shown in Figure 1. Figure 3, θ 1 = q 1 , θ 2 = q 2 , C and r represent the the hip joint angle, the knee joint angle, center, and radius of the reference trajectory (which can be defined as a circle), respectively. The lengths of the links are l 1 = 0.35 m and l 2 = 0.32 m; the mass and inertia of two links are m 1 = 1.8 kg, m 2 = 1.65 kg and I 1 = 1 4 m 1 l 2 1 kg·m, I 2 = 1 4 m 2 l 2 2 kg·m, respectively; the gravity constant is g = 9.801 m/s 2 . In addition, the center and radius are C = (0.5, 0) and r = 0.1 m, respectively.

As shown in
The parameters of the algorithm 1 (PASHS) are chosen as follows: ε = 10 −6 , ψ 0 = 10 −5 , γ 0 = 10 −3 , δ = 0.0001, ρ = 0.5, M 0 = , and the initial step size is selected as (Dai, 2011) .  Due to the proposed MPC model in section 3, the parameters of the MPC will be defined as follows: the prediction horizon is N = 5 and the control horizon is N u = 5; the sampling time is h = 0.01s. The initial state of end-effector position coordinates is y 0 = (0.6, 0), and the motion-task duration is 25 s with 5 s per cycle. The following experiments are conducted under different constraints of torque: while the hip joint angle and the angular velocity are constrained in interval − 2 3 π, 2 3 π , and [−π, π], the knee joint angle and the angular velocity are constrained in interval − 4 3 π, 0 , and [−π, π] (Jin and Zhang, 2011).

The Passive Rehabilitation Training With Different Torques Constraints
Example 1: Consider the following control torques constraints where τ 1 and τ 2 correspond to the torques of the hip joint and the knee joint of a two-link lower-limb rehabilitation robot, respectively. The numerical results of this situation are shown in Figures 2-4. Figure 2 represents the curves of angle and angular velocities of the hip and knee for the two-link lower-limb rehabilitation robot, and Figure 3 is the limit cycles of angle and angular velocities. From Figures 2, 3, it is further inferred that the angles and angular velocities of the lower-limb rehabilitation robot present the periodic properties; it also verifies that the proposed approach is feasible and effective. Figure 4A denotes the control torque vs. time, and the hip joint and knee joint of the rehabilitation robot can be controlled by torques τ 1 and τ 2 , respectively. As shown in Figure 4A, it can be seen that the control torques change periodically for two-link lowerlimb rehabilitation robot, which can help the injured patients to do rehabilitation training stably via Algorithm 1 (PASHS) with MPC technique. Figure 4B represents the tracking errors of the real position of end-effector and desired trajectory, while e x and e y are the tracking errors of horizontal ordinate and longitudinal coordinates. As shown in Figure 4B, the absolute value of the tracking errors is also smaller than 0.002 m, which also infers that the lower-limb rehabilitation robot could implement the passive rehabilitation training efficiently by the desired trajectory and MPC technique. It thus further  demonstrates that the theoretical analyses are feasible and reliable. Besides, it is very important to take into consideration the energy consumption of rehabilitation training in real-world rehabilitation implementations, therefore, control torques should be constrained within relatively reasonable bounds. As shown in Figure 4A, the optimal control input is obtained by online solving of the MPC problem via Algorithm 1 (PASHS). However, it can be seen from Figure 4A, that the bounded condition is too large for non-linear optimization problem with MPC model, that is, the real-time control torque τ 1 is around 13 N·m. In other words, to achieve the low-energy consumption, the boundary constraints can be reduced to around 10 N·m, which are utilized to control the two-link lower limb rehabilitation robot to realize the rehabilitation cycle movement of the injured lower limb.
The numerical results of this situation are shown in Figures 5-7. Figure 5 represents the curves of angle and angular velocities for the two-link lower-limb rehabilitation robot, and Figure 6 shows the limit cycles of angle and angular velocities, respectively. As can be seen from Figures 5, 6, the angular velocities of the lower limb rehabilitation robot are affected by the control torques with constraint conditions limited to −10-10 N·m. However, the injured limb can stably complete rehabilitation training activities via the algorithm 1 (PASHS) with MPC technique. Figure 7A shows the control torque vs. time, and it can be seen that the control input τ 1 could be constrained in between −10 and 10 N·m, and the smoothness of angular velocities maybe influenced by the constraint conditions. However, it further infers that the stability can be maintained for the two-link lower limb rehabilitation robot. Figure 7B plots the tracking errors of the real trajectories of end-effector and desired trajectories, while e x and e y are the same with the definition of Example 1. The proposed approach is therefore suitable for passive rehabilitation training of lower-limb rehabilitation robot.
Example 3: This example shows a comparison of different algorithms with the MPC solution.
In order to compare the advantages of PASHS algorithm, sequential quadratic programming (SQP) is selected to compare with our algorithm (Sun et al., 2020). The simulation problem is chosen as example 1, and the results are as follows: Figure 8A represents the horizontal ordinate tracking errors e x of lower-limb rehabilitation robots with passive rehabilitation training, and Figure 8B means the longitudinal ordinate tracking errors e y . From those two figures, it can be seen that the   tracking errors of SQP are almost > 0.1 m at every iteration, and sometimes the errors were closed to 1m. This is because an optimal solution for SQP may not be in a feasible region. However, the tracking errors of PASHS are always smaller than 0.01 at every time. The optimal solution of the PASHS algorithm was satisfied by the constraint conditions because of the projective matrix and active set. PASHS was therefore more suitable for the MPC of lower-limb rehabilitation robots. Figure 8C is the running time of two algorithms at every optimization; according to this figure, the running time of PASHS  was nearly always smaller than 0.2 s, and the SQP running time was around of 0.4 s. This can account for the real-time computing capability of PASHS algorithm with the MPC technique.
Example 4: The example with parameters' perturbation 1.5 times.
This example is reported the influence of uncertainties in system, and the 1.5 times parameters' perturbation is introduced into lower-limb rehabilitation robots. During the simulation, the parameters are selected as l 1 = 0.35 m, l 2 = 0.32 m, m 1 = 2.7 kg, m 2 = 2.475 kg, I 1 = 0.0827 kg·m, and I 2 = 0.0634 kg·m. Other conditions are the same as Example 1. The results of this example are shown as follows. Figures 9A,B represent the angles and angular velocities of lower-limb rehabilitation robot with 1.5 times parameters' perturbation, respectively. As can be seen from these figures, the angles and angular velocities are changed periodically and stably, although the model is disturbed. Figure 10A is the torques of lower-limb rehabilitation robot which is subjected to 1.5 times parameter perturbation, and Figure 10B represents the horizontal ordinate tracking errors e x and the longitudinal ordinate tracking errors e y of lower-limb rehabilitation robots with 1.5 times parameters' perturbation. From Figure 10A, we can find that torques are limited between −15 and 15 N·m because the mass of lower-limb rehabilitation robot is added. However, the absolute value of tracking errors are smaller than 0.0015 m according to Figure 10B. The PASHS algorithm could thus solve MPC problems of the lower-limb rehabilitation robot with uncertainties in the model, and high accuracy could also be guaranteed.

sEMG-Based Active Rehabilitation Training
In this subsection, the active intention of injured patients is regarded as one of the most important rehabilitation steps. Furthermore, the joint trajectories of injured lower limb can be identified via the mentioned ESN model based on the active motion intention, which can be seen as the desired trajectories of lower limb rehabilitation robots. A numerical simulation is illustrated and analyzed for two-link lower limb rehabilitation robot with ESN model and MPC technique. The technical diagram of sEMG-based active rehabilitation training and intention recognition is shown in Figure 11.
The active rehabilitation training consists of two parts. The first one is intention recognition, which collects and preprocesses raw sEMG signals, and then motion intention is identified by the ESN model. The details of first part is described in the following description.
During the data acquisition stage, a subject sits on a chair and swings the shank periodically. The sEMG signals of seven muscles of leg, which include the vastus rectus muscle (VR), semitendinosus muscle (SM), tibialis anterior muscle (TA), gastrocnemius muscle (GM), vastus lateralis muscle (VL), biceps muscle of thigh (BM), and extensor pollicis longus (EP), need to be recorded through data acquisition unit, respectively (Tong et al., 2014). The acquisition device is BIOPAC MP160, which can simultaneously capture eight channels of sEMG signals at the default 2 kHz sample rate. Angles and angular velocities of knee and ankle are recorded by inertial measurement unit (IMU), which selects 100 Hz as the sample rate. Due to the sample rate of sEMG signals is higher than IMU, the sub-sampling technology should be implemented by (Zhang et al., 2012) where sEMG pre (k) represents the sEMG signals after subsampling at k times. The position of the electrodes and raw sEMG signals are shown in Figure 12.
We then use neural network technology to establish the relationship between sEMG signals and motion state. This is due to the original sEMG signals being contaminated by different measurement noises, such as direct current bias and baseline noise (Law et al., 2011). The raw sEMG signals need to be preprocessed, which includes a high-pass filter with 50 Hz high cut-off frequency, full-wave rectification technology, low-pass filter with 5 Hz low cut-off frequency and normalized technology (Han et al., 2015). The sEMG signals can be seen as the input signals of neural network when the noise of raw sEMG signals is eliminated by the mentioned methods.
The ESN model is a kind of recurrent neural network that is composed of an input layer, a hidden layer where neurons interconnect randomly, and an output layer. The network architecture is depicted in Figure 13.
The mathematical model of the ESN model can be obtained as where F (·) is an activation function and commonly generates from F (x) = tanh(x); M X ∈ R l×l , M U ∈ R l×n , M F ∈ R l×m , M Y ∈ R m×l are the internal connection weight of the hidden layer, the input layer to the hidden layer connection weight matrix, the output layer to the hidden layer feedback weight matrix, and the hidden layer to the output layer connection weight matrix; X and Y are the echo state and output vectors of the ESN model, respectively. Assume that M X , M U , and M F are unmodifiable during the ESN model training (Pan and Wang, 2012). The ESN algorithm with off-line learning is summarized as follows.

Algorithm 2. (ESN learning algorithm)
Step 0. Initialize echo state X 0 and randomly obtain a matrix M. NormalizeM = M/|λ max | and generate M X = α XM , M U ∈ R l×n , M F ∈ R l×m , where α X < 1 is a spectral radiuses of M X and |λ max | is a spectral radius of M.
Step 1. Compute the echo state by for k = 0, 1, . . . , N, whereÛ k andŶ k are the kth input and output reference data from the training dataset.
Step 2. Collect the reservoir state X and target state Y as follows: Step 3. Off-line compute matrix M Y = (X + Y) T , where X + is the pseudo-inverse of X.  During the training process, the preprocessed sEMG signals are recorded by the Biopac MP160 system with seven channels from 0 to 46 s. Furthermore, the joint trajectories of injured lower limb are recorded by IMU, the data set from 0 to 32 s is then collected as a training set, and the remainder of the data set is regarded as a test set. The training set is utilized to train the ESN model, and the testing set is exploited to simulate lower limb rehabilitation robots with active rehabilitation training. The real angles and angular velocities are recorded by IMU, which aims at demonstrating the accuracy of proposed method. The ESN model has seven input neurons and four output neurons, the numbers of ESN hidden layer neurons are 100, and the parameters are selected as follows: α X = 0.5. The results of ESN training and testing are shown in Figures 14, 15. Figure 14 represents the intention recognition results of injured lower limb via ESN model, where θ 1 and θ 2 mean knee and ankle angles, respectively. Red solid lines denote real joint trajectories of knee and ankle angles, and blue solid lines represent training and testing results of ESN learning algorithm. Figure 15 shows the training and testing results of knee and ankle angular velocities through ESN learning algorithm.θ 1 anḋ θ 2 mean the angular velocities for knee and ankle, which are shown by blue solid lines. Red solid lines represent real angular velocities, which are recorded by IMU. As can be seen from Figure 14, the injured lower limb swings the calf at the 5th second, and then the knee periodically stretches and flexes about at about 41 s. When the knee joint angle reaches 0 rad, the knee joint swings to its maximum position. During the swing phase, the knee joint can flex more than −1.8 rad. The angle of the ankle joint is between 1.8 and 2.2 rad, which is always plantar flexion. Figure 15 shows that the angular velocity of the knee joint alternates between −2 and 2 rad/s with a cycle of about 4 s, while the angular velocity of the ankle joint varies slightly between −0.5 and 0.5 rad/s. The motion intention of injured lower limbs can be identified by the ESN learning algorithm with multichannel sEMG signals from 32 to 45.3 s. Meanwhile, it is also inferred that the proposed method shows superior performance for intention identification of injured lower limb.
The second one is an MPC problem; the joint trajectories of injured lower limb can be identified via an ESN model based on active motion intention, which is viewed as desired trajectories of the two-link lower limb rehabilitation robot. The control law  generated by the Algorithm 1 (PASHS) is transmitted to two-link lower limb rehabilitation robot, which aims at assisting patient to do rehabilitation training. The next predictive state also can be computed by the optimization results of MPC problem (2), which feeds back to the two-link lower limb rehabilitation robot system. Meanwhile, the MPC problem can be seen as follows: where k i,j k + i, j|k, and θ k ,θ k , and τ k represent the angle vector, angular velocity vector, and torque vector of two-link lower limb rehabilitation robot at prediction horizon and control horizon, respectively. θ d (k + i|k) means the desired trajectory at k + i time, and τ (k + j|k) is the control input increment. Q = 5I and R = I, where I is the identity matrix and the index N = 3 and N τ = 3. The trained ESN model can effectively identify the joint angle and angular velocity of injured lower limb from the multichannel sEMG signals. The lower limb rehabilitation robot takes the results of recognition as the desired trajectories. Combining the ESN model and MPC technique, the human-machine interactive control method is developed, investigated, and analyzed for lower limb rehabilitation robot and injured lower limb in this paper. Besides, to design the human-machine interactive controller, the joint angle and angular velocity that can be regarded as desired trajectories are identified by ESN learning algorithm from 32 to 46 s. The numerical results are shown in Figure 16.
Figures 16A-C represent the numerical results of angles, angular velocities, and torques of lower-limb rehabilitation robot during rehabilitation training process. As shown in Figure 16, the blue solid lines mean the results of knee and the red solid lines represent the results of ankle, respectively. In light of Figure 16, it can be seen that the patients can be stably driven to do rehabilitation motion by lower limb rehabilitation robot, which be oriented by a human's active intention. It is also inferred that the rehabilitation robot appropriately generates torques, which assist patients to do rehabilitation training and avoid the second injury. It is thus further demonstrated that it is very practical to train the injured lower limb through a human-machine interactive control method with multichannel sEMG signals. In other words, it is also verified that the ESN learning algorithm and Algorithm 1 (PASHS) are feasible and effective for the rehabilitation training of injured lower limb.

CONCLUSIONS
In this paper, to obtain an optimal controller of a nonlinear system, an MPC problem firstly solved by a new PASHS algorithm has been proposed and analyzed by exploiting the three-order Taylor discretization formula to linearize and discretize the constraint conditions. Furthermore, the PASHS approach not only takes advantage of a projected operator, but it also integrates the active set into HS conjugate gradient methods; the optimal controller can thus be rapidly solved for a non-linear optimization problem. Moreover, the feasibility and global convergence have been rigorously proved in this paper. Some numerical results have been presented and analyzed to substantiate the feasibility, effectiveness, and superiority of the developed human-machine interactive control method for passive/active rehabilitation training. The ESN model with multichannel sEMG signals also has been proposed for intention recognition, which could identify the joint angles and angular velocities of the injured lower limb to realize active rehabilitation training. In other words, passive rehabilitation makes patients train through fixed-based trajectories of injured lower limb; however, the desired trajectories of active rehabilitation training are identified by ESN learning algorithm with multichannel signals. Besides, combining with MPC technology and Algorithm 1 (PASHS), human-machine interactive control has been developed, investigated, and analyzed for two-link lower limb rehabilitation robot. The numerical results have inferred that the proposed method could be effectively applied to passive/active rehabilitation training. The proposed method has also solved a problem that creates uncertainty in the model. In future work, more effective and real-time methods will be developed and investigated in the solution of MPC problem and applied to the rehabilitation of patients, such as upper limb rehabilitation training, assisting patients to walk on the plane, or up and down stairs.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
TS: methodology, software, and review. YT: investigation and writing original draft. ZS: conceptualization and methodology. BZ: methodology. ZP: algorithm. JY: methodology and review. XZ: editing. All authors contributed to the article and approved the submitted version.