Next Article in Journal
Few-Shot Fine-Grained Image Classification via GNN
Next Article in Special Issue
Myoelectric Control Systems for Upper Limb Wearable Robotic Exoskeletons and Exosuits—A Systematic Review
Previous Article in Journal
A Probabilistic Bayesian Parallel Deep Learning Framework for Wind Turbine Bearing Fault Diagnosis
Previous Article in Special Issue
Soft Wearable Robots: Development Status and Technical Challenges
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Data-Driven Predictive Control of Exoskeleton for Hand Rehabilitation with Subspace Identification

1
Department of Engineering Management & Technology, University of Tennessee at Chattanooga, Chattanooga, TN 37403, USA
2
Department of Mechatronics Engineering, Marmara University, Istanbul 34744, Turkey
*
Author to whom correspondence should be addressed.
Sensors 2022, 22(19), 7645; https://doi.org/10.3390/s22197645
Submission received: 30 August 2022 / Revised: 4 October 2022 / Accepted: 7 October 2022 / Published: 9 October 2022
(This article belongs to the Special Issue Challenges and Future Trends of Wearable Robotics)

Abstract

:
This study proposed a control method, a data-driven predictive control (DDPC), for the hand exoskeleton used for active, passive, and resistive rehabilitation. DDPC is a model-free approach based on past system data. One of the strengths of DDPC is that constraints of states can be added to the controller while performing the controller design. These features of the control algorithm eliminate an essential problem for rehabilitation robots in terms of easy customization and safe repetitive rehabilitation tasks that can be planned within certain constraints. Experiments were carried out with a designed hand rehabilitation system under repetitive and various therapy tasks. Real-time experiment results demonstrate the feasibility and efficiency of the proposed control approach to rehabilitation systems.

1. Introduction

To regain the limb movement ability lost because of any disease or accident, a repetitive and intense rehabilitation process is required. Physiotherapists treat patients during this process. They frequently use rehabilitation robots to help patients perform the right movements with the right intensity while under control inside or outside of the clinic [1]. In recent years, academic or R&D studies were conducted frequently for design and implementation of rehabilitation robotics. Exoskeletons, which allow for direct limb manipulation, or end-effect robots, that support therapy by manipulating the limb’s distal point, are two mechanical structures that are used in rehabilitation robot design. Specially designed hand rehabilitation robots are used to treat post-stroke hand movement limitations. These robots need to be made in accordance with the complex structure of the hand, which has about 20 degrees of freedom; it should also be supported by a robust and adaptive control algorithm so that it can function consistently for each patient [2,3,4,5]. Devices that support active rehabilitation can be used during the patient’s completely lost movement during the passive phase of the repetitive rehabilitation process. The rehabilitation robot should perform fewer movement tasks as the patient recovers more and give them more responsibility. It should assist in this instance with active-assistive rehabilitation procedures. Resistance exercises for muscle strengthening can be done once the patient regains his mobility [6,7]. Robot control design must have an adaptive structure to perform the rehabilitation processes. Active power control is not possible for the direct trajectory control robots described in the literature. Rehabilitation robots are controlled to carry out active, active-assistive, and passive rehabilitation tasks using control algorithms, such as impedance or admittance control [8,9,10].
Systems with nonlinear structured uncertainties can also be controlled using variable structure controllers. These controllers allow for the use of parametric perturbations between lower and upper limits to deal with complexity and noise from the external world [11]. The control of a lower extremity exoskeleton was accomplished in [12] using the anti-disturbance sliding mode controller. Additionally, a robust adaptive sliding mode controller is proposed in [13] to deal with unknown and bounded dynamic uncertainties of the upper limb exoskeleton for passive rehabilitation tasks. A fuzzy approximation-based backstepping control is another approach to the control of an exoskeleton for rehabilitation [14,15].
System identification is crucial for establishing the controller design. The kinematic and dynamic parameters must be accurately determined, and the robot must generate the force required by the patient for it to adapt to the kinematic structure of each patient [16,17].
The limb that the robot works on simultaneously during rehabilitation is breakable and sensitive. While manipulating this limb, the robot must avoid harming it. Both the patient and the robot must cooperate, and the robot can only move the patient’s hand within limits. In this instance, the control algorithm should find the best solution—not just any solution—for the targeted task while taking the constraints into account. Algorithms for model predictive control can be employed to achieve this [2,18].
Model-based control algorithms, such as model predictive control, are frequently used in process control applications because they produce the optimal solutions possible, given the constraints. These algorithms can also be used for rehabilitation procedures because they use past data collected over a specific horizon to determine which estimation model offers the best solution [19,20]. For every patient, a unique model will be developed with unique system parameters and dynamics. A universal model cannot be used in this case. Since MPC algorithms are model-based algorithms, continuous identification data-driven predictive control algorithms can be applied in this situation [21].
The subspace-based parameter estimation approach is used in this study to estimate system parameters. The parameters appropriate for the model defined along a horizon are obtained using the state data and control signal gathered during open loop operation along this horizon [22]. A predictive control algorithm and an optimal control rule are built using the model and parameters acquired here. The control rule is handled independently of the constraints, and experiments are carried out to determine how much each model variable contributes to the control rule’s success.
The rest of the paper is organized as follows: Section 2.1 introduces the suggested control strategy, followed by Section 2.2 and Section 2.3, which describe the subspace identification method and data-driven predictive controller, respectively. The experimental setup is also described in the Section 2. The Section 3 of the paper includes the experimental results. This section examines and presents the impact of all parameters on the success of the DDPC. Finally, every result is discussed and concluded.

2. Materials and Methods

2.1. Exoskeleton for Hand Rehabilitation

The mechanical structure of the exoskeleton that was previously proposed is effective with the simple design. In this structure, the middle and proximal phalanxes serve as links for two 4-bar mechanisms that are sequentially coupled. Additionally, the metacarpal phalangeal (MCP) joint has a range of motion of 55 degrees of flexion, whereas the proximal interphalangeal (PIP) joint has a range of motion of 65 degrees of flexion from a fully extended posture. The linear actuator may move at a maximum speed of 12 mm per second and a maximum force of 45 N. A full hand opening or closing action requires roughly 5 s for a stroke of 50 mm. The mechanical structure and biodynamic fit of the hand make the designed system more practical in terms of usage and productivity than similar ones [2,23,24].
The system consists of L 0 ,   L 1 ,   L 2 linkages. The MCP joint angle ( ϕ 1 ) and PIP joint angle ( ϕ 2 ) rely on the length of the linear actuator ( l 0 ), as shown in Figure 1. Therefore, simultaneous actuation of both joints occurs. Whenever a new user (subject or patient) puts on the exoskeleton on their hand, the exoskeleton and the finger unite to form a single, distinctive 1-DOF system.

2.2. Proposed Control Method

In this study, position tracking was performed with a DDPC designed using subspace identification and model predictive control (MPC) techniques. MPC is a model-based technique that was successfully applied over the years. The difficulties (cost, time, and effort) associated with the identification of a predictive model of the system are major barriers that prevent the widespread adoption of MPC for complex systems.
It is challenging to describe the system model when using rehabilitation robots, orthotics, and prosthetic applications, because the controlled system’s parameters vary depending on the patient or the healing process. Therefore, data-driven algorithms are chosen over model-based algorithms. Model-based control is often a time-consuming and expensive process. Additionally, it is ineffective for a system that will be applied to multi-user processes with various mechanical characteristics. The benefit of the DDPC method is that it contains historical data continuously, making the system sensitive to any changes in mechanical parameters. The need for memory is made clear by the necessity of storing historical data for modeling. This may also be viewed as a drawback.
The input and output signals obtained from open loop data are used as input and output variables in the DDPC. The obtained sub-space matrices can be used to calculate the DDPC algorithm weight parameters. During the specified horizon, the system can keep track of the reference using a control rule based on past input and output data as well as tracking errors. By closely monitoring all activity along the horizon, the rule based on previous data will be able to react quickly. Particularly, the fast and accurate identification of new circumstances will enable quick adaptation. The structure of the proposed control method is shown in Figure 2.

2.3. System Identification Methods

2.3.1. Output Error Method for Identification

The purpose of the output error method is to find the best parametric model according to the given specific criteria. The criteria is the error between the measured noisy output and the simulated model output, as shown in Equation (1).
ϵ = x 2 x ^ 2
Here, x ^ 2 is the estimated state variable of the system calculated using estimated values of unknown dynamic parameters. If the parallel model configuration of the model is similar to 2, we can find the parameter with equations in 3.
x ^ ˙ 2 = a ^ x ^ 2 f ^ s s g n x ^ 2 + b ^ u ,      x ^ 2 0 = ( x ^ 2 ) 0
where a, b, and f s are the unknown parameters to be identified. In this equation, x 2 is the state variable, and u is the input of the system, and both are measured and/or observed.
a ^ ˙ = γ 1 ϵ x 2 ,     b ^ ˙ = γ 2 ϵ u ,    f ^ ˙ s = γ 3 ϵ s g n x 2  
where γ 1 , γ 2   and   γ 3 are the gain of the error effected to parameter. The parameters a, b, and f s , estimated by the output error method, were studied and compared in [2] previous studies.

2.3.2. Subspace Identification Method

Background information on the subspace identification matrices from open-loop data is provided in this subsection. The following section will use these matrices to design a data-driven predictive controller. We will begin by defining the system’s state-space model. As a result, the equations below can be written in state-space form for a linear discrete time-invariant system:
x k + 1 = A x k + B u k + K e k
y k = C x k + D u k + e k
where u k m , y k l , and x k m are the input variables, the output variables, and the state vector variables of the system, respectively; e k l is white noise disturbance. The system matrices A n × n , B n × m , C l × n , D l × m , and K l × l are the state, input, output, feed-through, and Kalman filter gain matrices of the system, respectively.
We assume that the measurements of the inputs u k and the outputs y k for k 1 ,   2 ,   ,   N are available for identification. The input Hankel matrices for u k are represented as U p and U f .
U p u 1 u 2 u N 2 M + 1 u 2 u 3 u N 2 M + 2 u M u M + 1 u N M + 1
U f u M + 1 u M + 2 u N M + 1 u M + 2 u M + 3 u N M + 2 u 2 M u 2 M + 1 u N
where the subscripts ‘p’ and ‘f’ represent the ‘past’ and ‘future’ matrices of the variables. Similarly, Hankel matrices for y k , represented as Y p and Y f defined as (8) and (9), respectively. The dimensions of the matrices are Y p ,   Y f M l × N 2 M + 1 and U p ,   U f M m × N 2 M + 1 , respectively.
Y p y 1 y 2 y N 2 M + 1 y y 3 y N 2 M + 2 y M y M + 1 y N M + 1
Y f y M + 1 y M + 2 y N M + 1 y M + 2 y M + 3 y N M + 2 y 2 M y 2 M + 1 y N  
These data block Hankel matrices are made rectangular in the subspace identification method to reduce the undesirable effects of noise on the identification system. This situation can be achieved via having a large set of data, denoted by the variable N . Moreover, M in Equations (6)–(9) can be understood as the order of the predictor equation. For a successful identification of the system behavior, the order M must be bigger or at least equal to the real system order n as manifested in the dimension of the state matrix A [25]. The system’s past and future state vectors are written as:
X p     x 1     x 2         x N 2 M + 1
X f     x M + 1     x M + 2         x N M + 1
As a result of the derivation of Equations (4) and (5), the equations can be written as below. These equations are known as the subspace matrix input–output equations used for identification [26,27].
Y f = Γ M X f + H M d U f + H N s E f
Y p = Γ M X p + H M d U p + H N s E p
X f = A M X p + Δ M d U p + Δ M s E p
Γ M M l × n can be described as the extended observability matrix, Δ M d n × M m as reversed extended controllability matrix (deterministic), and Δ M s   n × M l as the reversed extended controllability matrix (stochastic) [20,28]. Y f can be written with Equations (12)–(14) as below:
Y f = Γ m A m Γ m A m H m + Δ m Y p U p + H m U f + ( Δ M s A m Γ m H M s ) E p + H M s E f
Since the effect of E f is constant white noise, and cause of the stability of a Kalman filter, Equation (15) can be written to give an optimal prediction expression of the system output Y f as follows:
Y ^ f = L w W p + L u U f
where W p = Y p U p T , U f consists of past inputs and outputs and future inputs, respectively. L w M l × M l + m is the subspace matrix corresponding to the past input and output states, and L u M l × M m is the subspace matrix corresponding to the future inputs. Future output values in Equation (16), as well as the system’s future input, can be formulated as a linear combination of the system’s past input and output states. The system’s behavior will then be described using Equation (16), rather than going back to the identification techniques that yield the traditional transfer function or state-space description of the system.
The following least squares problem is solved to calculate L w and L u from the Hankel matrices.
m i n | | Y f L w L u W p U f | | F 2
This problem can be solved from the orthogonal projection of the row space of Y f into the row space of the matrices W p = Y p U p T . This can be defined by Equation (18) as follows:
Y ^ f = Y f / W p U f

2.4. Data-Driven Predictive Controller

The data-driven predictive control algorithm uses the linear subspace predictor and the cost function of MPC algorithm.
M and N are lengths of data. Furthermore, y d k + 1 , y k + 1 , u k M are the desired output r , output, and the input, respectively. All I/O data are stored in a database and then can be used again in control.
The MPC algorithm cost function form [29,30] can be written with the prediction and control horizon N p and N c equal to f as follow:
J = k p = 1 N p Y ^ t + k p r t + k p T W Q Y ^ t + k p r t + k p + k c = 1 N c Δ U t + k c T W R Δ U t + k
where W Q and W R are the weight matrices, r t is the reference signal at the current time t , N p and N c are the prediction and control horizon, respectively. N c maybe less than or equal to the prediction horizon N p   N c N p   o r   N c f .
We maintain to improve the basics of DDPC via rewriting the cost function of MPC Equation (19) in quadratic form. Using Equation (16) and the reference signal of r t , we can update the cost function as follows:
J = L w Δ W p + L u N c Δ U N c + Y t r t + 1 T W Q × L w Δ W p + L u N c Δ U N c + Y t r t + 1 + Δ U N c T W R   Δ U N c
If we solve the cost function, the control rule can be written as follows:
Δ U N c = L u N c T W Q L u N c + W R 1 × L u N c T W Q L w Δ W p + Y t r t + 1 = K Δ W p , N c Δ W p K e , N c Y t r t + 1 .
where K Δ W p , N c and K e , N c are the weight of the past data vector and the weight of the tracking error, respectively.
At each time condition, only the first element of Δ U N c is used to calculate the control input U t + 1 , which complies with Δ U t + 1 . Hence, abbreviating the first m rows in Equation (21) gives as below:
Δ U t + 1 = K Δ W p Δ W p K e Y t r t + 1
With,
K Δ W p = [ I m     0 m × M 1 m ] K Δ W p , N c
K e = [ I m     0 m × M 1 m ] K e , N c
where I m is an identity matrix of size m while 0 i × j is a zero matrix with i rows and j columns. Consequently, the control input U t can be written as follows:
U t = U t 1 + Δ U t

2.5. DDPC with Considering Constraints via Quadratic Programming

The ability of MPC and other predictive control algorithms to include constraints in the final control solution is one of their advantages. In this section, a constrained DDPC algorithm is provided, taking the constraints in Equation (26) into account.
F m Δ u m i n Δ u N c F m Δ u m a x F m u m i n u N c F m u m a x F l y m i n Δ y f F l Δ y m a x F l y m i n y f F l y m a x
Here F m   and   F l are defined as F m = I m I m I m T N c m   x   m ,   F l = I l I l I l T with identity matrices I m and I l . Further, u N c = u t + 1 T u t + 2 T u t + N c T T .
The control signal is optimized in the constrained DDPC algorithm while taking into account the constraints placed on the cost function specified in the earlier sections. The inequalities shown in Equation (27) are reached when the constraints defined in Equation (26) are rewritten as a function of [2,20].
I N c m I N c m Γ m Γ m L u N c L u N c Γ l L u N c Γ l L u N c A Q P Δ u N c   F m Δ u m a x F m Δ u m i n F m u m a x F m u t F m u m i n + F m u t F l Δ y m a x L w Δ w p F l Δ y m i n + L w Δ w p F l y m a x F l y t Γ l L w Δ w p F l y m i n + F l y t + Γ l L w Δ w p b Q P
It can be discovered by optimizing the cost function for the DDPC algorithm given in Equation (28) with constraints.
min Δ u N c 1 2 Δ u N c T H Δ u f + Δ u N c T f   s . t .   A Q P Δ u N c b Q P
The quadratic programming (QP) algorithm can be used to solve this optimization process. The QP algorithm determines the ideal control signal Δ u N c while accounting for the constraints.

2.6. Experimental Setup

For the execution of the parameter estimation algorithms, the Matlab/Simulink environment is used. STM32F4107 CPU was used to run control algorithms. The dual-channel H bridge L298 driver IC was used for DC motors and a 9V Li-Po battery was used as the voltage supply for the motor. The Matlab/Simulink environment was also utilized during the studies to gather and store data. Communication is established between the microprocessor and the Matlab/Simulink environment with the UART protocol. In each experiment, control signals or reference position values were sent from the Matlab/Simulink environment to the CPU using the serial communication protocol. The employed experimental setup is shown in Figure 3.
The data-driven predictive control method was validated on a real-time system by different experiments. The system position was directly measured using a linear potentiometer.

2.7. Passive and Active Rehabilitation

The design of both active and passive rehabilitation tasks requires estimation or measurement of the external force exerting on the exoskeleton’s endpoint. In this study, a micro load cell is placed between end point of the linear actuator shaft and fork joint to measure external force ( F e x ) . External force can be used to stimulate a virtual mass-spring-damper system, as shown in Figure 4 and its mechanical parameters can be adjusted depending on the type and degree of rehabilitation required.
When the measured external force is applied to the virtual system using Equation (29), the virtual system’s response x d , can be found.
m d x ¨ d + b d x ˙ d + k d x d = F e x
For active rehabilitation task, the virtual system’s position response, x d , can be used to deviate to the desired reference from the actuator’s actual position, x 1 (Equation (30)). In this instance, the behavior of the controller is a regulator and keeps the actuator in its actual stroke position. If the external force is greater than zero, the desired reference is different from the actual position of the actuator.
x r = x 1 x d
For passive rehabilitation, x d can be used to create a deviation from the predefined trajectory x r , as shown in Equation (31), and the controller can make use of this desired reference.
x r = x r x d
It is possible to decide how the virtual mass-spring-damper system responds to external force and carries out passive, active, or assistive rehabilitation tasks by setting the parameters within the acceptable range of values.

3. Results and Discussion

3.1. Experimental Results of Subspace Prediction Algorithm

During the tests, the robot was not subjected to any external force and the exoskeleton is tested on a healthy human hand. The predicted model was compared with the state-space model obtained by the output error method; those applications and results are explained in past studies [23]. By using the model horizon parameters p = 30 ,   f = 10 , the subspace estimation model parameters ( L u r and L w r ) are calculated. The L u and L w obtained with the same values as the previous p and f from the system are then calculated using the u t input signal to be tested, the L u and L w are then calculated using the same values as the previous p and f from the system (Figure 5. In the graphics, the results calculated using the reference ( L u r and L w r ) models and the actual L u and L w models are compared.
The output error-based model (OEbM) output is compared with the sub-space model (SPbM) estimation result obtained with the u input signal ( u = A 0 s i n ω 0 t ,     ω 0 = 0.5   rd / sn ve A 0 = 5 ) shown in the Figure 6. The linear actuator’s stroke length, which is considered the system output, is represented by the y-axis on the graph as position.
The percentage error function in Equation (32) was obtained between the OEbM output, which was used as the reference model, and the SPbM output to evaluate the outcomes of the model estimate test, and this number was chosen as the performance criterion.
e = i = 1 n y O E b M y S P b M y O E b M n
This experiment’s error value (e) was calculated as 1.1871%. Table 1 provides the percentage error values of the test results with additional specified input signals.
The range of ω 0 = 0.3   rd / s , where the input signal runs throughout the full stroke in a single alternate, and ω 0 = 2   rad / s , which gives 5% displacement on the motor stroke, is tested using a single frequency component. The first eight experiments revealed a linear correlation between the modeling success and the parameters of the frequency component of the input signal. As a result, success increased and error decreased at higher frequency values, while success increased at lower frequency values. Since the full stroke length could not be studied in the experiments with high frequency u input signals, it can be argued that all the system’s characteristics could not be apprehended in them. This has an impact on modeling success.
The following experiments are performed with a sinusoidal input signal with two separate frequency components ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t ) . The experiment that was performed using input signals with ω 0 = 0.2   rd / s and A 0 = 2 , ω 1 = 0.8   rd / s and A 1 = 3 (Experiment 10) results are shown in Figure 7.
It was noted that in experiments using the u input signal with two frequency components, the average success rate increased. When experiments 9 through 12 in Table 2 are examined, it becomes clear that the error function produces results that are similar across the ranges of experimental parameters. It was noted that the error function is negatively impacted by the separation between the two component frequencies.
The following experiments (13–16) use u input signals that have three distinct frequency components ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t + A 2 s i n ω 2 t ) , as shown in Figure 8. The experiments conducted are more useful than the earlier experiments, as shown in Table 3. Different frequency components are found to increase the success of modeling.
The following tests were conducted using input signals that contained a scanning frequency as shown in Figure 9. These input signals are configurations for input signals that begin with low-frequency components and increase throughout the designated range. The average success is higher in these experiments, as shown in Table 4. This input signal, which has components in several different frequency ranges, helps to clarify the system’s characteristics.

3.2. Experimental Results of Data-Driven Predictive Controller

The predicted model parameters and control parameters K e   and   K Δ w p 1 x 2 p are calculated. The step function and a sinusoidal trajectory response (2 π radians from position 0 to 50%) are examined to analyze the parameters affecting the control algorithm, and the results are discussed.

3.2.1. Effect of Data Length on Control Success

Figure 10 shows the relationship between the calculated control coefficients and the control success as a function of the data size of the input signal (u), which is used as the estimation input signal. In experiments, only the input signal u’s data length is changed, while the other parameters, p = 30 , f = 10 , Q = 1 , R = 2 ,   and   N c = 5 , remain constant. As a result, it was found that as the number of data increases, overshoot decreases, increasing the success of the control. The control coefficients from experiments with the same frequency components but fewer data are seen to negatively affect the success as the number of data decreases. The overshoot increases as N decreases. The quantity of data has no meaningful effect on the rise time.

3.2.2. Analysis of Model Horizon Parameters (p, f)

The system state equations along the chosen horizon are the basis of the subspace-based estimation model. Using the input signal from the most successful experiment (Experiment 20) from the previous test results, the test results referring to the horizon parameter success characteristic are presented in this section.
As a result of parameter estimation using various historical model horizon values, the control parameters are displayed in Table 5. Experiments are carried out with constant parameters f = 10 , Q = 1 , R = 2 ,   and   N c = 5 by varying the model horizon p between 30 and 60. The effect of the previous model horizon on the system response is investigated in these experiments. Rise time and percent overshoot are used as performance factors in step function experiments. The success factor in trajectory tracking is the mean of the squares of the errors. Table 5 details all these values.
Without considering the experiment taken as p = 50 in the table, it can be said that the other values show that the step response and trajectory tracking success are negatively impacted by an increase in the p value as shown in Figure 11. However, in the experiment where p is set at 50, success increases once more. The control is seen to be totally broken in the following experiment. This assumes that there might be a linear relationship between control success and the past modeling horizon, p. When the response of the system to the sine waveform is analyzed, it is seen that even the worst result (p = 60) follows the trajectory with a certain error (mse = 16.96 from Table 5).
In the following experiments, the past model horizon p, and other parameters ( p = 30 , Q = 1 , R = 2 ,   and   N c = 5 ) are held constant to examine the control success of the future model horizon f. Rise time and percent overshoot are considered performance factors in experiments using the step function. The average of the squares of the errors is used as a performance metric for trajectory tracking success. In Table 6, each of these values is described in detail.
When the tables and graphics are examined, it is observed that there is a nonlinear relationship between the model future horizon f and control success. The overshoot and rise times are close to each other in experiments where the f value of the future horizon is between 20 and 25. The increase in the mean value and standard deviation of Ke and K Δ w p has a negative effect on the rise time. The rise time increases as the f parameter increases, but the average overshoot decreases. The response of trajectory also shows that the overshoot is high in the experiment where f is set to 5 (Figure 12). This overshoot is a result of the Ke value obtained in this experiment being much lower than in the other experiments. Based on gathered data, it is seen that, even if the number of future horizons is chosen as only 5, the steady-state stability is considered as critically stable.
The tables display control coefficients based on the parameters used in the experiments. The graph in Figure 13 shows the control parameters (Ke) for all p and f values in the determined range (5–100). When the Ke value is between 0.2 and 0.3, it is apparent from the experiments, the results of which are evaluated, that the success rate is high. Additionally, the tables display the mean and standard deviation values of the K Δ w p parameter calculated in successful results. Figure 13 and Figure 14 can be used to analyze the suitable numeric selections of p and f parameters for these values. The color scale of the p and f pairs providing successful Ke values is shown in Figure 13. It is clear from all the graphs that the values of p and f for the successful control task can be selected from a range of 5 to 80 and 20 to 60, respectively.

3.3. The Effect of Q and R Parameters on the Control System Succession

It is shown in equation 16 that the parameters Q and R in the general MPC cost function are the weights of the reference error and the differential control signal, respectively. These weights have a direct impact on the success of the control because the data-driven predictive control rule is based on the MPC cost function. The success of trajectory tracking and step response are examined using the experiments in this section with the Q and R parameters. With p = 30, f = 10, and Nc = 5, various Q and R parameters are tested in experiments.
The Q parameter and the rise time have a linear relationship, as can be seen when Table 7 is examined. In addition, the ratio of the Q parameter to R also affects the rise time. It is clearly seen that the R parameter affects the overshoot as seen in Figure 15. The overshoot increased up to a maximum value of 8% in the experiment where R was set to be 50. These tests led to the conclusion that selecting a small Q value would have a positive effect on the rise time as shown in Figure 16. The experiments also confirmed that selecting a low value (Q = 1) and choosing the parameter R within the optimal ranges (10,200) have a positive impact on the response of step and trajectory tracking.
The control horizon Nc in the data-driven predictive control algorithm is described in Equation (16). Unlike f, which is the model horizon, Nc specifies the size of the space from which the system control signal u c should be calculated. It is used by the control algorithm to compute the system control signal u N c , which is the first component of Δ u N c along Nc. These tests were conducted with the following parameters: p = 30 ,   f = 10 ,   Q = 1 ,   and   R = 2 . The impact of Nc on step function response and trajectory tracking response is shown in Table 8 and Figure 17. It is clear from the tables and graphs that the rise time has a linear relationship with Nc’s trajectory tracking response and step function response.

3.4. Passive and Active Rehabilitation

The patient is completely passive during passive rehabilitation tasks. The patient’s hand is completely guided by the exoskeleton. The patient did not exert any force on the exoskeleton for the first 10 s, as shown in Figure 18. The controller in this case uses x r = x d . When a maximum counterforce of 7–8 N was applied in the final 10th second, the virtual system response x d deviates from the reference x r . This makes it possible for the system to react to the patient in accordance with any unexpected issues that might arise on the patient’s finger. The response stiffness of the system to the patient’s hand can be adjusted by appropriately adjusting the parameters of the virtual mechanical system as shown in Equation (26).
The controller functions as a regulator during the active rehabilitation process and tries to keep the current position of the actuator stroke. The reference shifted, as shown in Figure 19, because of the patient applying an external force of about 6 N in the 12 th second. The patient made the flexion movement with a stable force applied to the system. The extension movement was carried out by exerting force in the opposite direction after the 14th second. The system can be operated at greater forces by modifying the virtual mechanical system’s parameters. The system can now perform in resistive mode as a result. The system operates in assistive mode by taking x d in negative.
The study of the robustness of MPC can be approached in a variety of ways. The first focuses on the closed-loop systems’ robustness when created utilizing the nominal system (i.e., neglecting uncertainty). The second makes an effort to achieve robustness within the context of conventional model predictive control by taking into account all feasible realizations of the uncertainty. The third approach addresses this by introducing feedback in the min–max optimal control problem solved online [31]. In this study, according to the results obtained in some experiments, it was evaluated that the system is robust against uncertainties. For example, in the analysis of model horizon parameter experiments, it was observed that the system followed the trajectory even at the worst coefficients (p = 60, f = 5). The model control values, Ke and KΔwp, are determined using this model and calculated as optimum values, showing that the controller is robust enough to handle uncertainties.
The subspace prediction estimation results to be made with the ongoing past data online will guarantee that the controller operates feasibly and continuously during the rehabilitation tasks. Along with model parameters calculated throughout a specific horizon with sub-space prediction, the model includes all model uncertainties associated with the patient’s exoskeleton and measurement noises related to force and position. The use of data-driven estimation methods and model-free control algorithms will improve the effectiveness of studies with patients who can be evaluated in a wide range of spectrums as opposed to computing the model of a biomechanical structure, such as the human hand, with conventional approximate methods.

4. Conclusions

In this work, a data-driven predictive control is developed for a hand exoskeleton robot used for rehabilitation. The designed control rule was used in a set of experiments, and the results were presented. The experiments are intended to examine how the parameters affecting the suggested control algorithm influence the success of the control. A data-driven predictive control algorithm is optimization-based and certain constraints are added to the problem during the optimization process; it then suggests the best solution within the boundaries set by those constraints. The rehabilitation process aims to regain the patient’s lost mobility by having them perform repeated exercises that are suited to their situation.
In our experiments to evaluate the performance of the proposed control algorithm, data length is investigated for subspace prediction, and it is expressed that, within a certain range, data length was linearly related to modeling success. On the success of controlling for the DDPC rule, the effects of the p, f, Nc, Q, and R parameters are examined and discussed separately. When all the test results are evaluated, we can conclude that the suggested solution is suitable for rehabilitation processes because it offers the best solutions, while still considering some limitations.
It was shown that the exoskeleton controller operates in passive, active, and assistive modes with the benefit of an auxiliary reference created using the measured external force.

Author Contributions

Conceptualization, E.K. and G.A.; methodology, G.A.; software, G.A.; validation, E.K. and G.A.; formal analysis, E.K.; investigation, G.A.; resources, E.K.; data curation, E.K.; writing—original draft preparation, G.A.; writing—review and editing, E.K.; visualization, G.A.; supervision, E.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Ethical review and approval were waived for this study due to no applicable.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Yue, Z.; Zhang, X.; Wang, J. Hand rehabilitation robotics on poststroke motor recovery. Behav. Neurol. 2017, 2017, 3908135. [Google Scholar] [CrossRef] [PubMed]
  2. Akgun, G.; Cetin, A.E.; Kaplanoglu, E. Exoskeleton design and adaptive compliance control for hand rehabilitation. Trans. Inst. Meas. Control 2020, 42, 493–502. [Google Scholar] [CrossRef]
  3. Du Plessis, T.; Djouani, K.; Oosthuizen, C. A review of active hand exoskeletons for rehabilitation and assistance. Robotics 2021, 10, 40. [Google Scholar] [CrossRef]
  4. Reitan, I.; Dahlin, L.B.; Rosberg, H.E. Patient-reported quality of life and hand disability in elderly patients after a traumatic hand injury—A retrospective study. Health Qual. Life Outcomes 2019, 17, 148. [Google Scholar] [CrossRef] [Green Version]
  5. Dovat, L.; Lambercy, O.; Gassert, R.; Maeder, T.; Milner, T.; Leong, T.C.; Burdet, E. HandCARE: A cable-actuated rehabilitation system to train hand function after stroke. IEEE Trans. Neural Syst. Rehabil. Eng. 2008, 16, 582–591. [Google Scholar] [CrossRef]
  6. Kabir, R.; Sunny, M.S.H.; Ahmed, H.U.; Rahman, M.H. Hand Rehabilitation Devices: A Comprehensive Systematic Review. Micromachines 2022, 13, 1033. [Google Scholar] [CrossRef]
  7. Perry, J.C.; Trimble, S.; Machado, L.G.C.; Schroeder, J.S.; Belloso, A.; Rodriguez-de-Pablo, C.; Keller, T. Design of a spring-assisted exoskeleton module for wrist and hand rehabilitation. In Proceedings of the 2016 38th Annual International Conference, Orlando, FL, USA, 16–20 August 2016. [Google Scholar]
  8. Tsoi, Y.H.; Xie, S.Q. Impedance control of ankle rehabilitation robot. In Proceedings of the 2008 IEEE International Conference on Robotics and Biomimetics, Bangkok, Thailand, 22–25 February 2009; pp. 840–845. [Google Scholar]
  9. Li, X.; Liu, Y.H.; Yu, H. Iterative learning impedance control for rehabilitation robots driven by series elastic actuators. Automatica 2018, 90, 1–7. [Google Scholar] [CrossRef]
  10. Culmer, P.R.; Jackson, A.E.; Makower, S.; Richardson, R.; Cozens, J.A.; Levesley, M.C.; Bhakta, B.B. A control strategy for upper limb robotic rehabilitation with a dual robot system. IEEE/ASME Trans. Mechatron. 2009, 15, 575–585. [Google Scholar] [CrossRef]
  11. Jafarov, E.M.; Parlakçi, M.N.A.; Istefanopulos, Y. A new variable structure PID-controller design for robot manipulators. IEEE Trans. Control Syst. Technol. 2004, 13, 122–130. [Google Scholar] [CrossRef]
  12. Mo, L.; Feng, P.; Shao, Y.; Shi, D.; Ju, L.; Zhang, W.; Ding, X. Anti-Disturbance Sliding Mode Control of a Novel Variable Stiffness Actuator for the Rehabilitation of Neurologically Disabled Patients. Front. Robot. AI 2022, 9, 9. [Google Scholar] [CrossRef]
  13. Riani, A.; Madani, T.; Benallegue, A.; Djouani, K. Adaptive integral terminal sliding mode control for upper-limb rehabilitation exoskeleton. Control Eng. Pract. 2018, 75, 108–117. [Google Scholar] [CrossRef]
  14. Chen, Z.; Li, Z.; Chen, C.P. Disturbance observer-based fuzzy control of uncertain MIMO mechanical systems with input nonlinearities and its application to robotic exoskeleton. IEEE Trans. Cybern. 2016, 47, 984–994. [Google Scholar] [CrossRef]
  15. Rahmani, M.; Rahman, M.H. An upper-limb exoskeleton robot control using a novel fast fuzzy sliding mode control. J. Intell. Fuzzy Syst. 2019, 36, 2581–2592. [Google Scholar] [CrossRef]
  16. Yu, H.; Huang, S.; Chen, G.; Thakor, N. Control design of a novel compliant actuator for rehabilitation robots. Mechatronics 2013, 23, 1072–1083. [Google Scholar] [CrossRef]
  17. Chen, J.; Huang, Y.; Guo, X.; Zhou, S.; Jia, L. Parameter identification and adaptive compliant control of rehabilitation exoskeleton based on multiple sensors. Measurement 2020, 159, 107765. [Google Scholar] [CrossRef]
  18. Ulkir, O.; Akgun, G.; Nasab, A.; Kaplanoglu, E. Data-Driven Predictive Control of a Pneumatic Ankle Foot Orthosis. Adv. Electr. Comput. Eng. 2021, 21, 65–74. [Google Scholar] [CrossRef]
  19. Maciejowski, J.M. Predictive Control: With Constraints; Pearson Education: London, UK, 2002. [Google Scholar]
  20. Mardi, N.A. Data-Driven Subspace-Based Model Predictive Control. Doctoral Dissertation, RMIT University, Melbourne, Australia, 2010. [Google Scholar]
  21. Ulkir, O.; Akgun, G.A.Z.İ.; Kaplanoglu, E. Real-time implementation of data-driven predictive controller for an artificial muscle. Stud. Inform. Control 2019, 28, 189–200. [Google Scholar] [CrossRef] [Green Version]
  22. Rossiter, J.A. Model-Based Predictive Control: A Practical Approach; CRC Press: Boca Raton, FL, USA, 2017. [Google Scholar]
  23. Akgün, G. Data Driven Predictive Control of Exoskeleton for Rehabilitation. Doctoral Dissertation, Marmara Universitesi, Istanbul, Turkey, 2019. (In Turkish). [Google Scholar]
  24. Akgun, G.; Kaplanoglu, E.; Cetin, A.E.; Ulkir, O. Mechanical design of exoskeleton for hand therapeutic rehabilitation. Quest J. J. Res. Mech. Eng. 2018, 4, 9–17. [Google Scholar]
  25. Van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems: Theory—Implementation—Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  26. Kadali, R.; Huang, B.; Rossiter, A. A data driven subspace approach to predictive controller design. Control Eng. Pract. 2003, 11, 261–278. [Google Scholar] [CrossRef]
  27. González, A.H.; Ferramosca, A.; Bustos, G.A.; Marchetti, J.L.; Fiacchini, M.; Odloak, D. Model predictive control suitable for closed-loop re-identification. Syst. Control Lett. 2014, 69, 23–33. [Google Scholar] [CrossRef]
  28. Meidanshahi, V.; Corbett, B.; Adams, T.A., II; Mhaskar, P. Subspace model identification and model predictive control based cost analysis of a semicontinuous distillation process. Comput. Chem. Eng. 2017, 103, 39–57. [Google Scholar] [CrossRef]
  29. Venkat, A.N.; Hiskens, I.A.; Rawlings, J.B.; Wright, S.J. Distributed MPC strategies with application to power system automatic generation control. IEEE Trans. Control Syst. Technol. 2008, 16, 1192–1206. [Google Scholar] [CrossRef]
  30. Kaplanoğlu, E.; Arsan, T.; Varol, H.S. Predictive control of a constrained pressure and level system. Turk. J. Electr. Eng. Comput. Sci. 2015, 23, 641–653. [Google Scholar] [CrossRef]
  31. Mayne, D.Q.; Rawlings, J.B.; Rao, C.V.; Scokaert, P.O. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789–814. [Google Scholar] [CrossRef]
Figure 1. Designed exoskeleton for hand rehabilitation.
Figure 1. Designed exoskeleton for hand rehabilitation.
Sensors 22 07645 g001
Figure 2. Structure of proposed data-driven predictive controller block diagram.
Figure 2. Structure of proposed data-driven predictive controller block diagram.
Sensors 22 07645 g002
Figure 3. Experimental setup.
Figure 3. Experimental setup.
Sensors 22 07645 g003
Figure 4. Virtual mass-spring-damper system.
Figure 4. Virtual mass-spring-damper system.
Sensors 22 07645 g004
Figure 5. Subspace prediction validation procedure.
Figure 5. Subspace prediction validation procedure.
Sensors 22 07645 g005
Figure 6. Selection of the sinus input signal ω 0 = 0.5   rd / s and A 0 = 5 , and comparison of results. (a) Single sinus signal, (b) Comparison of results.
Figure 6. Selection of the sinus input signal ω 0 = 0.5   rd / s and A 0 = 5 , and comparison of results. (a) Single sinus signal, (b) Comparison of results.
Sensors 22 07645 g006
Figure 7. Results of experiment 10: (a) input and output signal, (b) results.
Figure 7. Results of experiment 10: (a) input and output signal, (b) results.
Sensors 22 07645 g007
Figure 8. Results of experiment 10: (a) input and output signals, (b) results.
Figure 8. Results of experiment 10: (a) input and output signals, (b) results.
Sensors 22 07645 g008
Figure 9. Results of experiment 20: (a) input and output signal, (b) results.
Figure 9. Results of experiment 20: (a) input and output signal, (b) results.
Sensors 22 07645 g009
Figure 10. Step response according to u signals of different N data lengths.
Figure 10. Step response according to u signals of different N data lengths.
Sensors 22 07645 g010
Figure 11. Success of different past model horizon (p) values in tracking trajectories and success of response of step function.
Figure 11. Success of different past model horizon (p) values in tracking trajectories and success of response of step function.
Sensors 22 07645 g011
Figure 12. Effect of different future model horizon (f) values on the success of tracking the step response and trajectory.
Figure 12. Effect of different future model horizon (f) values on the success of tracking the step response and trajectory.
Sensors 22 07645 g012
Figure 13. Ke and K Δ w p calculated for the 5–100 range of p and f values.
Figure 13. Ke and K Δ w p calculated for the 5–100 range of p and f values.
Sensors 22 07645 g013
Figure 14. The standard deviation value of K Δ w p calculated for the 5–100 range of p and f values.
Figure 14. The standard deviation value of K Δ w p calculated for the 5–100 range of p and f values.
Sensors 22 07645 g014
Figure 15. Effect of the R parameter on step and trajectory response.
Figure 15. Effect of the R parameter on step and trajectory response.
Sensors 22 07645 g015
Figure 16. Effect of the Q parameter on step and trajectory response.
Figure 16. Effect of the Q parameter on step and trajectory response.
Sensors 22 07645 g016
Figure 17. According to Nc, the response of step function and trajectory.
Figure 17. According to Nc, the response of step function and trajectory.
Sensors 22 07645 g017
Figure 18. Passive rehabilitation task.
Figure 18. Passive rehabilitation task.
Sensors 22 07645 g018
Figure 19. Active rehabilitation task.
Figure 19. Active rehabilitation task.
Sensors 22 07645 g019
Table 1. Experiments with a single frequency input signal u = A 0 s i n ω 0 t .
Table 1. Experiments with a single frequency input signal u = A 0 s i n ω 0 t .
Experiments ω 0   rd / s A 0 Error (e)
10.350.0187
20.450.0266
30.550.0337
40.650.0395
50.750.0727
60.850.0825
7150.1431
8250.9273
Table 2. Experiments with a two-frequency u input signal ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t ) .
Table 2. Experiments with a two-frequency u input signal ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t ) .
Experiment ω 0   rd / s ω 1   rd / s A 0 A 1 Error (e)
90.30.4230.0127
100.20.8230.0267
110.20.8140.0360
120.50.1220.0235
Table 3. Experiments with a three-frequency u input signal ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t + A 2 s i n ω 2 t ) .
Table 3. Experiments with a three-frequency u input signal ( u = A 0 s i n ω 0 t + A 1 s i n ω 1 t + A 2 s i n ω 2 t ) .
Experiment ω 0   rd / s ω 1   rd / s ω 2   rd / s   A 0 A 1 A 2 Error
130.30.40.72220.0088
140.30.40.12220.0239
150.30.40.12250.0251
160.50.40.12250.0211
Table 4. Experiments with input signals that configured by scanning frequency.
Table 4. Experiments with input signals that configured by scanning frequency.
Experiment k 0 A 0 e (Error)
170.00340.0149
180.00550.0200
190.0150.0083
200.0130.0048
Table 5. Control parameters calculated using p.
Table 5. Control parameters calculated using p.
Model Horizon (p)Ke Mean   of   K Δ w p 1 x 2 p Standard   Deviation   of   K Δ w p 1 x 2 p Response of Step FunctionTracking Performance
Rising Time (s)Overshoot %mse
300.19080.07300.24993.880.174.3138
400.28530.08520.24935.770.2018.4100
500.37460.09040.23993.840.181.4763
600.37800.10450.2202* None* None16.9600
* In this experiment, no rise or overshoot was observed during the experiment.
Table 6. Control parameters defined using f.
Table 6. Control parameters defined using f.
Future Horizon (f)Ke Mean   of   K Δ w p 1 x 2 p Standard   Deviation   of   K Δ w p 1 x 2 p Response of Step FunctionTracking Performance
Rising Time (s)Overshoot %mse
50.04760,00660.02223.380.09629.0828
100.19080.07330.24003.850.03224.3138
150.21910.09960.29463.940.02201.2027
200.23850.15800.39424.360.00183.9737
250.16880.11880.29804.480.025410.5345
500.19510.18970.46226.910.052460.1988
Table 7. Effect of Q and R parameters on control response.
Table 7. Effect of Q and R parameters on control response.
QRKe Mean   of   K Δ w p 1 x 2 p Standard   Deviation   of   K Δ w p 1 x 2 p Response of Step FunctionTracking Performance
Rising Time (s)Overshoot %mse
520.76940.19690.8251* None* None288.1818
420.56930.15130.61845.6540.038543.5554
320.49750.13540.55126.6310.008465.2003
220.43290.12020.47763.0930.068833.0765
120.25890.08060.29813.0310.005011.8783
130.18970.06480.22693.5280.01186.9044
150.12930.05100.16493.1890.02863.7065
1100.08070.03980.11523.4930.04902.7420
1200.05490.03370.08893.9020.06262.4434
* In this experiment, no rise or overshoot was observed during the experiment.
Table 8. Effect of Nc parameter on control response.
Table 8. Effect of Nc parameter on control response.
N c Ke Mean   of   K Δ w p 1 x 2 p Standard   Deviation   of   K Δ w p 1 x 2 p Response of Step FunctionTracking Performance
Rising Time (s)Overshoot %mse
50.16880.11880.29803.2320.025410.5345
100.17990.07980.24293.2240.01528.5312
150.22200.07680.26693.2100.00843.9609
200.25890.08060.29813.1800.008411.8783
250.34230.08700.35623.1920.068417.3509
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kaplanoglu, E.; Akgun, G. Data-Driven Predictive Control of Exoskeleton for Hand Rehabilitation with Subspace Identification. Sensors 2022, 22, 7645. https://doi.org/10.3390/s22197645

AMA Style

Kaplanoglu E, Akgun G. Data-Driven Predictive Control of Exoskeleton for Hand Rehabilitation with Subspace Identification. Sensors. 2022; 22(19):7645. https://doi.org/10.3390/s22197645

Chicago/Turabian Style

Kaplanoglu, Erkan, and Gazi Akgun. 2022. "Data-Driven Predictive Control of Exoskeleton for Hand Rehabilitation with Subspace Identification" Sensors 22, no. 19: 7645. https://doi.org/10.3390/s22197645

APA Style

Kaplanoglu, E., & Akgun, G. (2022). Data-Driven Predictive Control of Exoskeleton for Hand Rehabilitation with Subspace Identification. Sensors, 22(19), 7645. https://doi.org/10.3390/s22197645

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop