actuators Control of a Rehabilitation Robotic Device Driven by Antagonistic Soft Actuators

: Stroke is becoming a widely concerned social problem, and robot-assisted devices have made considerable contributions in the training and treatment of rehabilitation. Due to the compliance and continuous deformation capacity, rehabilitation devices driven by soft actuators are attached to widespread attention. Considering the large output force of pneumatic artiﬁcial muscle (PAM) and the biological musculoskeletal structure, an antagonistic PAM-driven rehabilitation robotic device is developed. To fulﬁll the need for control of the proposed device, a knowledge-guided data-driven modeling approach is used and an adaptive feedforward–feedback control approach is presented to ensure the motion accuracy under large deformation motion with high frequency. Finally, several simulations and experiments are carried out to evaluate the performance of the developed system, and the results show that the developed system with the proposed controller can achieve expected control performance under various operations. Author Contributions: Conceptualization, H.C. and Q.R.; methodology, H.C.; software, H.C.; validation, H.C., H.S. and W.L.; formal analysis, W.L.; investigation, H.S. ; resources, Q.R.; writing—original H.C. H.S.; writing—review and editing, W.L. Q.R.; supervision, H.C.; administration, funding


Introduction
Stroke, a disease caused by the insufficient blood supply to the brain due to blockage or sudden rupture of blood vessels, is becoming a leading cause of mortality and disability. The World Health Organization (WHO) describes stroke as "rapidly developing clinical signs of focal (or global) disturbance of cerebral function, with symptoms lasting 24 h or longer or leading to death, with no apparent cause other than of vascular origin" [1]. According to medical reports, stroke has been the second largest cause of death globally after ischemic heart disease, and up to now, the burden of the stroke remains heavy [2].
As a troubling disease, stroke usually causes a loss of physical strength on one side of the body, paralysis or hemiplegia. More than 2/3 patients affected by stroke suffer from impaired upper limb motor function, which brings great inconvenience to daily life [3]. Medical studies have found that an efficient way to treat stroke is the timely rehabilitation training. Moreover, many trials have shown that task-specific, high-intensity exercises in an active, functional, and highly repetitive manner can enhance motor recovery [4]. However, the traditional exercise therapy performed by the therapists is always tedious and inefficient [5]. To assist patients in accurately, quantitatively and personally completing the therapy, rehabilitation robotic devices are being widely adopted recently.
Most of the rehabilitation robotic devices are actuated by "hard" actuators such as electric motors. These kinds of actuators are able to provide sufficient power and accurate motion. However, limitations such as heavy weight and confined motion direction restrict their usage in rehabilitation. Soft actuators, in contrast, have better biomimetic qualities due to their compliance and flexibility.
Many soft actuators for robotic applications have been investigated, such as pneumatic artificial muscle (PAM) [6][7][8], shape memory actuators (SMA) [9,10], dielectric elastomer ac-tuators (DEA) [11,12] and soft fluidic actuators (SFA) [13][14][15]. Similar to biological muscles, once stimulated, these soft actuators can arise deformation responses such as stretching, bending, tightening, or expanding. Hence, these soft actuators are generally compatible with robotic applications involving human body interaction. Among the actuators explored above, PAMs stand out in rehabilitation robotic applications [16], because PAMs possess high power-to-weight ratios and prominent inherent compliance [17]. Moreover, compared with other aforementioned soft actuators, PAMs are able to generate a greatly larger output force. Thus, many rehabilitation robotic devices actuated by PAMs have been developed.
Although the PAMs adopted in the aforementioned robotic devices exhibited musclelike motion in a rehabilitation process, control of these robotic devices still remains a huge challenge. In the exercise training process, due to the nonlinearity, viscoelasticity, and timedependent properties of soft materials, precision motion control of a PAM-driven robotic device in large deformation motion with high frequency, which is normally necessary for rehabilitation, is very difficult. Previously, a lot of studies have been conducted to improve the robustness control system of PAMs-driven robotic devices, such as a basic proportionalintegral-derivative (PID) control approach proposed in [24]. A modified sliding mode control (SMC) approach [25], an adaptive fuzzy logic control (FLC) approach [26] and an adaptive recurrent neural networks (ARNN) control approach [27] were proposed to handle the nonlinearity problem of the system. However, all these control schemes mentioned above are suitable with low-frequency motion, and few of these approaches can be applied to robotic applications in practice under large deformation with fast response.
Biologically, due to the similarities between PAM and biological muscle, as well as the purpose of assisting patients, an intuitive solution for the design of rehabilitation robots' structure and control system is to mimic the biological system. Most motions of living beings are driven by the muscles in antagonistic pairs, thus one objective of this paper is to investigate the control issues of PAMs in antagonistic pairs. Actually, antagonistic PAMs have been extensively studied, such as the study presented in [28]. However, none of them were designed entirely referring to biology as biceps and triceps. The main purpose of adopting antagonistic pairs is to generate bidirectional actuation and make up for the deficiency of passive recovery characteristics.
Further on, referring to the biological neural control mechanisms, a controller inspired by the cerebellar computing mechanism called cerebellar model articulation controller (CMAC) is adopted in this paper. The CMAC is a kind of adaptive neural network for approximating complex nonlinear functions. It is essentially a technology for nonlinear function mapping, which is a feedforward network virtually. Due to its local learning strategy, CMAC can improve the learning speed while ensuring the approximation performance with the nonlinear functions. Furthermore, CMAC is suitable for continuous input space and is not sensitive to the order of given data as well. Therefore, CMAC is widely used in real-time control of the nonlinear system under complex dynamic conditions [29,30].
This work mainly aims to design a rehabilitation robotic device that can be equipped on the upper limb for stroke patients under an accurate control. The main contributions of this study are listed as follows: (i) Drawing from the biological system, antagonistic pairs in the bio-inspired structure are designed to overcome the uncontrollable passive recovery characteristic, which mimic the coupling relationship of living beings' antagonist muscles. In this structure, the system is able to achieve bidirectional actuation, and the motion can be more compliant with human dynamics. (ii) A knowledge-guided data-driven modeling approach by using the spring-dashpot model is employed to simplify the difficulties in the modeling process, where the model parameters will be determined via system identification. (iii) To achieve good control performance, an adaptive feedforward-feedback controller consisting of a PID controller and a CMAC is proposed, which can perform well in shortening the response time and improving the motion accuracy. Significantly, the effect of CMAC is corresponding to the cerebellum in human beings, and the PID controller is used as error-based feedback control.
The rest of this paper is organized as follows. The design of the robot prototype with the kinematics and dynamics analysis is presented in Section 2. Next, the dynamic models of the PAM, as well as the antagonistic pairs, are established in Section 3. In Section 4, the adaptive feedforward-feedback control scheme is designed and analyzed. Next, both simulations and experiments are carried out to verify and validate the feasibility of the developed robot system in Section 5. Finally, Section 6 gives the conclusions.

Robotic Device
In this section, the mechanical structure and its working principle are explained, and then its kinematics and dynamics analysis is given.

Rehabilitation-Arm Design
Usually, muscles in biological systems are placed in antagonistic pairs, which are more flexible and faster in accomplishing assigned tasks compared with a single muscle. Figure 1a shows three common antagonistic configurations, which are pulley type, lever type and pendulum type, respectively. Each type of configuration has its own characteristics, and they have been widely used in engineering practice. However, when looking for the optimal solution from biological systems, a different configuration illustrated in Figure 1b is found. Inspired by this configuration, the structure for the rehabilitation pneumatic robotic arm is designed. Thus, a robot manipulator shown in Figure 1c is developed in this paper. The motivations of adopting this bio-inspired structure are: (i) The robot manipulator is designed for stroke patients' rehabilitation, this structure is similar to physiological structure of human arms. It could mobilize the muscle groups of the arm to a greater extent in principle, and thus it is able to achieve better rehabilitation results. (ii) Unlike conventional antagonistic pairs, in this paper, the soft and deformable actuators can achieve antagonistic effect even if there is a distortion during transmission, which is able to obtain good operational feasibility in practice.  For the sake of illustration, the upper arm is defined as B1 while the forearm as B2. As mentioned above, the McKibben PAM P1 represents the biceps and P2 represents the triceps. Then, P1 and B1 are connected at N 1 , P1 and B2 are connected at N 2 , P2 and B1 are connected at N 4 , and P2 and B2 are connected at N 3 . The distances between the endpoints of P1/P2 and the corresponding connection points are defined as d 1 , d 2 , d 3 , d 4 .

Kinematics and Dynamics Analysis
In the initial state shown in Figure 2b, P1 is divided into two sections with the center of the elbow joint (rotation axis) as the reference, the lengths of these two sections are noted as d 5 , d 6 . Similarly, d 7 , d 8 could be obtained with P2. Define the length of P1 as L 1 and the length of P2 as L 2 , then the initial lengths of PAMs can be expressed as When the robot manipulator rotates to an angle of θ as shown in Figure 2c, it can be seen obviously that both PAMs are deformed to different angles. According to geometric principles, the kinematics can be obtained. Firstly, the relationship between the length of P2 and the rotation angle θ can be simply expressed as The length of P1 can be obtained through the Pythagorean theorem: Since the human arm is naturally drooping, the rotational moment of the forearm is considered as where F 1 and F 2 represent the force driven by P1 and P2, respectively, l F 1 and l F 2 are the corresponding moment arm. M G is the representation of the gravitational moment, and l F is the vertical distance between the force direction F 1 and the forearm B2. l O G is the distance between O and O G shown in Figure 2b. According to the Law of Cosines, l F can be obtained from the following equation In addition, the dynamic model derived via Euler-Lagrange equation for this robot manipulator is where J is the inertia coefficient, C is the damping coefficient and G(θ) is the gravitational torque. Thus, the following equation can be derived where J(θ,θ) includes the torque caused by gravity and other external disturbances.

System Modeling
Modeling is important for the control scheme design [31,32]. In this section, the system modeling of the PAM and the robot manipulator are presented.

Modularization Design
Before modeling, the selection of modeling objects should be determined. In this paper, there are two objects: one is the whole robotic system, and the other is the single McKibben PAM.
If the former object is chosen, the input-output relationship can be established directly through a complete system model, which is really intuitive and would reduce the complexity of the internal mechanism analysis of the system to a certain extent. There is no denying that this approach is well suited to the data-driven modeling method.
For the latter one, only the actuator's input-output relationship can be obtained. In this way, a second-step modeling process should be performed by integrating the kinematics and dynamics model of the antagonistic structure. Through this method, easier controller design can be achieved, and thus at least the real-time state of each actuator is predictable, which makes the interaction process safer.
In this paper, the modeling on each actuator is adopted. On the one hand, medical devices naturally require higher accuracy and safety guarantee, so a kind of completely controlled medical device can better meet patients' needs for rehabilitation and provide more personalized and appropriate treatment schemes. As for the extension-flexion elbow motion mentioned above, the degrees of function loss in biceps and triceps in different patients are different in the general sense. Therefore, when carrying out rehabilitation training movements, it is necessary to regulate the auxiliary strength of muscles according to the principle of "assist as needed" [33], rather than taking the motion trajectory as the only evaluation index. On the other hand, the single actuator model could be treated as a basic module, and in this way, the modular platform construction can be realized. By coupling a lot of basic modules through antagonism or other connection modes, diverse platforms can be built, not limited to the robotic elbow joint manipulator described in this paper, which is a great foundation for future work. In addition, the core value of modular robots lies in their various structural forms, which enables the robots to generate many types of movements and adapt to different environments.

Knowledge-Guided Modeling
In the process of modeling the single actuator model, a simplified spring-dashpot model mentioned in [12] is employed to describe the single-dimensional physical properties of PAM because of the similarity with biological muscle in viscoelasticity. The modeling method here is similar to our previous work [34,35]. This method has been verified to be a feasible solution, and this work aims to verify that this kind of modeling method has generality through the modeling work of different objects whose properties are similar to biological muscles.
As shown in Figure 3, the decisive factors in the modeling process are the model parameters (i.e., linear spring stiffness k, spring stiffness k h , viscous friction coefficient c h , and system order N). These parameters are eventually determined through the system identification. F active refers to the active force produced by the system. In this paper, it is the contraction force generated by the PAM driven with high pressure. F passive refers to the passive damping force, including resistance force and other environmental stresses during the motion. Since the fixed environmental force can be directly included in the model parameters for simplification, only resistance force ζ would be concerned when considering the external effects. The model uncertainties are mainly due to the equivalent massm of the robotic system and the resistance force ζ during the motion. Hence, the dynamic equations for this system are given bỹ In this paper, the best fitting effect arises when N = 4. Hence, the state-space model can be expressed asẊ where x is the displacement of the linear spring, x h represents the spring deformation in the h-th spring-dashpot module (h = 1, 2, · · · , N, N is the system order), k h represents its spring stiffness, and c h shows its viscous friction coefficient, F active is the active driving force, ζ illustrates the resistance force during the motion, and F passive contains all the environmental stresses,m is the equivalent mass of the whole system.

System Identification
According to (8) shown above, the system input is the resultant force F active , and the output is the displacement x. However, in the actual operation process, it is difficult to provide the force signal directly. In the actuation with PAM, the actual input signal is pressure, and its output includes contraction force and deformation. Thus, the conversion relationship between force F active and pressure P should be established before system identification.
As can be seen in the actuation process of PAM, the contraction force increases as the input pressure increases, and it is directly associated with the end point's movement under the unstrained case. In this paper, the end point's movement equals to the system's displacement, so it can be inferred that the contraction force is affected by both the pressure and the current displacement.
The force-pressure relationship of the PAM actuator is measured under isometric constraint, the experimental results are shown in Figure 4. Obviously, in the case of consistent displacement, the force-pressure curve has a strong linear characteristic, so (18) can be rewritten as The specific expression of the force-pressure curve can be roughly determined through curve fitting method. The results are shown in Figure 5, whose coefficients of determination are R 2 = 0.9878.
Thus, the final expression of the force-pressure relationship is After the definition of the F active − P relationship, the system identification can be carried out. In this part, the data collected are mainly the input pressure and the output displacement. Here, the collected pressure data would be converted into actual contraction force data, which will be directly substituted into (8) for identification. In the experiment for system identification, a signal presented as a sinusoidal sweep force signal with the frequency of 0.2 to 1 Hz and the amplitude of 280 N is firstly applied to the PAM. To further improve the range of identification as well as the applicability, the PAM is also driven by another identification sinusoidal signal with the amplitude of 30 to 300 N and the frequency of 0.5 Hz. Following that, the System Identification Toolbox in MATLAB is used to estimate the parameters in (8) and (9). Figure 6a shows the identification result. Besides the identification, several regular signals (such as square wave and triangle wave) with a low frequency and a high frequency are also used to validate the identified model, the results of which are shown in Figure 6b,c.  Generally, the identified parameters are listed in Table 1. The identification accuracies of each curve are shown in Table 2.   Since the displacement x has been taken as the actual output, the final output of the two PAMs can be expressed as As can be known from (7), the model of each PAM in the antagonistic pairs can be written as where sign(ẋ k ) returns 1 or −1 depending on whetherẋ k is positive or negative, and Thus, wherex k = [x k1 x k2 x k3 x k4 ] T , and P k is the current inflation pressure of the k-th PAM.
Considering model uncertainties and disturbances, the following model is obtained wheref k (x k ,ẋ k ),h k (ẋ k ) andf ka (x k ) are the nominal models denoted by system identification, δ f k (x k ,ẋ k ), δh k (ẋ k ) and δ f ka (x k ) are the model uncertainties, and d k (x k ) is the external disturbance. Lumping the model uncertainties, the external disturbances and the coupling effect together as D k (x k ,ẋ k ) ≡ δ f k (x k ,ẋ k ) + δh k (ẋ k ) + δ f ka (x k )P k + d k (x k ) + δ k , which is unknown but bounded. Specially, the bounded condition can be due to the mechanical constraints and the actuation limit.
Remarkably, the system can be considered as a bilinear system, which is an important subclass of nonlinear system [36]. Furthermore, the nonlinearity can come from the model uncertainties due to the system nonlinearity and the nonlinear components in the distur-bance. Thus, each PAM in the antagonistic pairs can be treated as a nonlinear system with unknown but bounded uncertainties and disturbance.

Control Scheme
In this section, the control scheme for the developed robotic system is proposed and designed in detail.

Feedforward-Feedback Control Scheme
To ensure the control accuracy as well as a quick response, a feedforward-feedback control scheme is proposed in this paper, the block diagram of which is illustrated in Figure 7. Here, the feedback controller mainly refers to the error-based controller. It is worth noting that both PAMs are in the same control structure and therefore their control laws are significantly similar. From (25), it shows that the ideal system input can be expressed as where k p and k d are two positive constants,C(e k ,ė k ) is the function of the motion error. As shown in the equation, the ideal controller can be regarded as a controller composed of feedforward control and feedback control essentially. u f eed f orward represents the feedforward control for desired motion, u f eedback represents the feedback control for position error, and δD denotes the compensation signal for uncertainties and disturbances.
Thus, a feedforward-feedback control scheme is suitable for the system developed in this paper. Moreover, a basic PID controller is adopted as the feedback controller due to its advantages of simple structure and easy implementation [37], the detailed design of which is shown below. u pk = K k,P e k + K k,I e k dt + K k,Dėk , where K k,P , K k,I , K k,D are the proportional gain, the integral gain, the derivative gain, respectively, and e k is the motion error of the k-th actuator.

Cerebellar Model Articulation Controller
In this paper, the selection of the PAM actuator is due to its similarities to biological muscles, and the modified antagonistic structure is inspired from biological system. Therefore, developing a controller with bioinspiration can be a way to control this bio-inspired robot manipulator. In the spirit of this bio-inspired actuator, an artificial neural network controller based on the model of the mammalian cerebellum, namely, CMAC can be a good candidate [38].
Significantly, the functions of CMAC in the control scheme are mainly to approximate the feedforward controller mentioned above as well as to compensate for the disturbances because this neural network controller can offer good learning ability to approximate the system and even the nonlinear part. Meanwhile, the use of PID controller along with the CMAC can help the control system to stabilize the robot in the early phase which is good to improve the control system performance [39].
As observed from Figure 8, the CMAC is composed of three layers (input layer, association cell layer and response output layer) with two mappings (concept mapping and actual mapping). Then, the CAMC could be regarded as a pair of mappings given by where U denotes the input vectors, AC and AP represent Access Concept Memory and Access Physical Memory, respectively, and y denotes the output vector. For the PAM, the control scheme with CMAC is illustrated in Figure 9. There are two inputs used in U: the displacement feedback information x and the desired motion reference x d . Both signals are processed through the mappings and then the synaptic weights are located. Finally, the sum of the activated weights is sent as the control output to the system. Meanwhile, the weights are adjusted based on a rule with the objective of minimizing the motion error.  According to the cerebellar model, the detailed design of the CMAC is shown as follows.
(1) Input Quantization. To convert the continuous inputs into finite and discrete sets, the input vectors must be quantized. The uniform quantization interval is computed by where s = [s 1 s 2 ] T = [x x d ] T represents the input vector, s i,max = max(s i ) and s i,min = min(s i ) are the maximum and minimum of s i , respectively, r i ∈ N + determines the quantization resolution, and q i is the number of the quantization interval. In this step, the input space is divided into different regions. Each region is corresponding to a unique weight. Upon to the activation of the region, its corresponding weight will be activated.
(2) Association. To find out the regions needed to be activated corresponding to the current input vector, the Gaussian function is used as the basis function in this paper, which is expressed as where ϕ ij i (s i ) represents the quantization distance of the i-th input s i to the j i -th region, µ ij i represents the end points of the j i -th region and σ ij i represents the width of the j i -th region. The number of the regions for the input s i could be defined as n i = q i + 1. In this paper, the µ ij i ) and σ ij i are set according to the quantization interval q i , where σ ij i is set equally in each region, which means σ ij i = σ i . Evidently, a larger n i can cause a higher resolution, but a more complex computation. The basis function used in this paper is defined as ), j = 1, 2, · · · , n r and j i = 1, 2, · · · , n i , (31) where n r = n 1 × n 2 denotes the number of the respective fields. In another way, the mapping vector can be expressed by the form of (3) Output. The output of the CMAC is the sum of the activated weights in the weight memory w: where u c = u CMAC is the output of controller to the system,w = [w 1 w 2 · · · w n r ] T ∈ R n denotes the estimated weight vector. (4) Weight Update. To adjust the weights in each cycle, the δ learning rule [40] is employed, then the weight adjustment index is defined by where ε =ë + k dė + k p e is the filter error and e = x d − x is the motion error. The objective of CMAC is to minimize the adjustment evaluation function Q(w). By applying the gradient descent algorithm, the learning law is given bỹ where p denotes the p-th iteration and η is the learning rate. Substitute u c into (25), the error model of the system can be expressed aṡ where e u =ū − u c is the approximation error between the desired control outputū and the actual CMAC control output u c , E = [eė] T represents the error state. Thus, we havē Hence, ∂Q ∂w Then, the learning law is given bỹ In the designed CMAC, because the learning law is related to the filter error, which is determined by the motion error, an adaptive control is formed.
As the two actuators have a similar control law, the CMAC controller could be given by where k = 1, 2 represents the first actuator PAM1 and the second actuator PAM2, respectively. Hence, the whole controller, combining the adaptive CMAC and the PID controller, is shown as follows. u k = u ck + u pk =w T k φ k + K k,P e k + K k,I e k dt + K k,Dėk . (41)

Simulation and Experimental Results
In this section, both simulations and experiments are carried out to evaluate the control performance of the proposed feedforward-feedback control scheme.

Simulation
The parameters used in the simulation are given in Table 1. According to the demand in medical rehabilitation of stroke patients, a proper trajectory is the periodic repeatable motion between the target position and the initial relaxation position (θ = 0).
For the purpose of verifying the control performance of the proposed system, the following three simulations are conducted to evaluate the control performance.
As mentioned in Section 4, the feedforward-feedback control scheme is composed of CMAC and PID controllers. A sinusoidal trajectory with an amplitude of 0.08 m is used in the simulation.
For the comparison purpose, the control schemes and controllers are applied to the system: (i) PID-only, (ii) CMAC-only, and (iii) feedforward-feedback. Figures 10-12 shows the control results with the PID-only, CMAC-only, and feedforward-feedback control schemes, respectively.  Figure 12. Sinusoidal wave trajectory tracking performance with feedforward-feedback control scheme.
As shown in Figure 12, the PID-only control scheme can achieve a good tracking performance, but there is a certain response lag in each cycle. By comparing the PID-only control scheme with the CMAC-only control scheme, the CMAC-only control scheme can also help the system to track the desired trajectory, but the trajectory tracking error is smaller while using the PID-only control scheme.
After applying the feedforward-feedback control scheme, combining the CMAC and PID controller, the tracking performance is shown in Figure 12, the improvement of the response is apparent. From the results, both response speed and accuracy are satisfactory. This is because the system nonlinearity can be well handled by the CMAC and the feedforwad control can help to shorten the response time. This shows the advantages of applying the CMAC in the proposed control scheme.
Considering the practical application (repetitive rehabilitation), the saturated sinusoidal trajectory is used as the desired motion trajectory. As shown in Figure 13, the displacement (length) of PAM switches periodically between 0.01 m and 0.045 m, mimicking the deformation process of biological muscle in a repetition motion. The result shown in Figure 13 is consistent with the previous simulation results. It can achieve good performance in terms of response speed and accuracy in each cycle, but some errors still exist at the peaks of the trajectories. (2) PAMs in antagonistic pairs control simulation.
After the simulation of a single PAM, it is necessary to evaluate the control performance of the manipulator with antagonistic pairs. Similarly, with the consideration of practical application, the desired trajectories are set as the saturated sinusoidal trajectory. Furthermore, with the concern of the special case, both sides of rotation should be tested. Here, the positive rotation ranges from 2.5 • to 60 • , and the negative rotation from 2.5 • to −10 • (the manipulator that hangs naturally has a slight angle instead of 0 • ) are tested. Figure 14 shows the simulation results on the trajectory tracking with the manipulator using the proposed control scheme. As can be seen, both desired trajectories in the bidirectional rotation can be well tracked. By comparing these results to the control of the single PAM, it can be found out that the response speed and accuracy can be improved with the assistance of the antagonistic pairs.

(3) PAMs in antagonistic pairs control simulation with uncertainty and disturbance.
To evaluate the robustness of the proposed system, simulation while the uncertainty and disturbance are added to the system is carried out, respectively. Two situations are tested in this simulation, the first one is to increase the equivalent mass of the whole system to be two times of the identified equivalent mass, which represents the model uncertainty, and the second one is to introduce an input disturbance into the system after about five seconds. Figure 15 shows the control system tracking performance under these two situations. As can be observed from the figure, the proposed control scheme is still able to help the system to track the desired trajectories well no matter there is an uncertainty or a disturbance exist in the system. Therefore, it can be concluded that the proposed control scheme can achieve guaranteed robustness. In all, the simulation results indicate that the proposed feedforward-feedback control scheme performs well on the motion trajectory tracking of the rehabilitation pneumatic robotic device.

Experiments
To further validate the effectiveness and feasibility of the proposed control scheme, real-time experiments of motion trajectory tracking have been conducted in a practical experimental system setup shown in Figure 16. In this setup, two PAMs (FESTO, DMSP-10-200N-RM-CM) are used in antagonistic pairs to construct the robot manipulator. Two solenoid switching valves (MHE2-MS1H-3/2G-QS4-K) are used to control inflation/deflation of the PAMs and four solenoid valves (MHJ10-S-2.5-QS-4-MF) are used to adjust the airflow to the PAMs. Furthermore, two pressure sensors (SPTE-P10R-Q4-B-2.5K) are used to measure the pressure inside the two PAMs, respectively. Moreover, an OptiTrack motion capture system (with five cameras) is used to measure the displacements of the PMAs and the rotation angle of the manipulator. Lastly, a computer is used to receive the sensing information, implement the proposed control scheme, send out the control input signals to the robot manipulator. Considering that keeping the arm in the elbow flexion state for a period of time is more effective in the arm rehabilitation training, the sinusoidal wave is replaced by four more appropriate trajectories in the real-time experiment. Four different experiments are carried out to explore the effectiveness of the proposed rehabilitation pneumatic robotic device with antagonistic pairs. The results are evaluated by the displacement accuracy of the PAMs and the rotation angle error of the manipulator.
In this paper, the antagonistic pairs are prepared for the large displacement and/or high-speed operations, thus this property should be tested. As shown in Figure 17a, a trapezoidal wave with a high frequency of 2.6 Hz is set as the desired trajectory, ensuring that the manipulator rotates from 3.2 • to 25 • periodically.
From the results, the actual motion tracks well with the desired trajectory in spite of little lag. Due to the high frequency, the response is not timely in the rising edge stage, but it still can achieve the target position finally at the steady state.
This experiment is designed to verify the feasibility in various large displacement situation. In this part, a square signal with changing amplitude of around 20 • is applied. As shown in Figure 17b, the manipulator performs well with the low frequency and large amplitude desired trajectories, but it is still a bit slow at the rise time. Considering that the overshoot and the steady state error are more important for a rehabilitation device, the tracking performance on this square trajectory is satisfactory for the rehabilitation application. Besides that, the gradually changing amplitude can assist the patients in increasing their motion range gradually rather than abruptly.

(3) Low-frequency convex trajectory:
To achieve a better auxiliary effect, it is necessary to design a trajectory according to the biological behavior and rehabilitation goals. The purpose of the previous two experiments is to verify the effectiveness of the control scheme and the manipulator, and its evaluation is determined by the response speed and accuracy. However, it is not sufficient for practical medical applications. For example, the speed of motion should be limited to ensure the smoothness of the trajectory and avoid the cause of secondary injury during the recovery process. The rehabilitation course should also be gradual, which means that each course of treatment should be improved on the basis of the last course of treatment, it cannot be abrupt.
Referring to the above considerations on the practical application, a convex trajectory was designed and the rotation rate was limited (a speed controller as the main circuit is added, forming a kind of cascade control). As shown in Figure 17c, the system has a fast response and is able to reach the desired position eventually. In the tracking process, the actual angle of the manipulator increases gently, which puts little load on the arm and provides a period of adaptive phase. It should be pointed out that the rehabilitation device is essentially a passive device, and thus the human-robot interaction should be emphasized when assisting body movements. Otherwise, the therapeutic effect can hardly be improved, and it might cause even worse malignant effects.
According to the results of this experiment, it can be determined that the proposed system could meet the design expectation for actual application with the low-frequency convex trajectory. To be specific, the manipulator will enlarge its motion range. In this experiment, assuming that the patients have already adapted to the motion range from 0.2 • to 16 • , then extending the upper range of the motion to 19 • , and thus the motion range from 0.2 • to 19 • , is set as the new goal of rehabilitation. Considering the exercise load and rehabilitation efficiency, the motion phase is divided into two steps. In the first step, the angle rotates from 0.2 • to 16 • . Then in the second step, the angle rotates from 16 • to 19 • . In this way, the patients only need to concentrate on the rehabilitation training with the motion range from 16 • to 19 • . After adapting to the new extended range, the integrated square-like trajectory will be applied for a new training. Repeat the above process, the patients would increase their exercisable motion range step by step.

(4) High-frequency convex trajectory.
Based on the third experiment, a repeated experiment is carried out with a highfrequency convex trajectory in order to further demonstrate the feasibility of rehabilitation. In the third experiment, the reason to choose convex trajectory has been explained, and the design of the rehabilitation project has also been introduced. However, a crucial transition phase should be added in the above rehabilitation process, i.e., the high-frequency motion stage. Apparently, the low-frequency movements occupy a major part in the rehabilitation to allow the patients to adapt to various motions. However, the accomplishment of high-frequency movements can help to reflect the recovery effect. Thus, a gradual training process from low frequency to high frequency is indispensable for any motion in rehabilitation.
Hence, in this experiment, another convex trajectory with a higher frequency but the same amplitude is applied, the results of which are shown in Figure 17d. From the figure, the actual motion can achieve the desired target at the end of the response. There is a difference between the desired trajectory and the actual trajectory due to the added design with the motion speed limitation. In spite of this, the system performs well on the motion trajectory tracking.
In summary, although the rise time is a bit slow, the steady state is very accurate where the steady state errors under the high-frequency trajectory and the low-frequency trajectory are within 0.8 • and 0.5 • , respectively. According to the results and analysis shown above, it can be revealed that the proposed control scheme is effective for the motion trajectory tracking with the developed manipulator, which helps the system achieve good tracking performance for the medical application.

Conclusions
In this paper, the control issues of a rehabilitation robotic device driven by soft PAMs for stroke patients has been investigated. The main challenge of the control issues lies in the difficulty of dealing with the nonlinearity and the passive recovery characteristic. To overcome the challenges, an antagonistic coupling structure and a feedforward-feedback control scheme are proposed. Moreover, for medical applications, a rehabilitation robotic device driven by two PAMs in antagonistic pairs is designed.
To verify the feasibility of the proposed medical rehabilitation robotic device, both theoretical analysis as well as simulations and experiments are conducted. The results show that the designed robotic system can achieve a good control performance with the proposed control scheme. Furthermore, the applicability with a more complex and flexible medical rehabilitation robotic device will be studied in future work.

Conflicts of Interest:
The authors declare no conflict of interest.