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

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.


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

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 L0, L1, L2 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.

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 modelbased 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 multiuser 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.

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.

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 modelbased 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 multiuser 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.

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).
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.
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.
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.

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 = Ax k + Bu k + Ke k (4) where u k ∈ R m , y k ∈ R l , and x k ∈ R m are the input variables, the output variables, and the state vector variables of the system, respectively; e k ∈ R l is white noise disturbance. The system matrices A ∈ R n×n , B ∈ R n×m , C ∈ R l×n , D ∈ R l×m , and K ∈ R 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 .
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 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: 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].
Γ M ∈ R Ml×n can be described as the extended observability matrix, ∆ d M ∈ R n×Mm as reversed extended controllability matrix (deterministic), and ∆ s M ∈ R n×Ml as the reversed extended controllability matrix (stochastic) [20,28]. Y f can be written with Equations (12)- (14) as below: 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:Ŷ where W p = Y p U p T , U f consists of past inputs and outputs and future inputs, respectively. L w ∈ R Ml×M(l+m) is the subspace matrix corresponding to the past input and output states, and L u ∈ R Ml×Mm 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.
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:

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: 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 or 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: If we solve the cost function, the control rule can be written as follows: 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: With, 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:

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.
Here F m and F l are defined as F m = I m I m . . . I m T ∈ R N c m x m , F l = I l I l . . . I l T with identity matrices I m and I l .
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].
It can be discovered by optimizing the cost function for the DDPC algorithm given in Equation (28) The quadratic programming (QP) algorithm can be used to solve this optimization process. The QP algorithm determines the ideal control signal ∆u Nc while accounting for the constraints.

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.

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 ( ). External force can be used to stimulate a virtual mass-springdamper 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 , can be found.
For active rehabilitation task, the virtual system's position response, , can be used to deviate to the desired reference from the actuator's actual position, (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.
For passive rehabilitation, can be used to create a deviation from the predefined trajectory , as shown in Equation (31), and the controller can make use of this desired reference.

= −
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. 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.

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 ex ). External force can be used to stimulate a virtual massspring-damper system, as shown in Figure 4 and its mechanical parameters can be adjusted depending on the type and degree of rehabilitation required. 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.

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 ( ). External force can be used to stimulate a virtual mass-springdamper 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 , can be found.
For active rehabilitation task, the virtual system's position response, , can be used to deviate to the desired reference from the actuator's actual position, (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.
For passive rehabilitation, can be used to create a deviation from the predefined trajectory , as shown in Equation (31), and the controller can make use of this desired reference.
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. When the measured external force is applied to the virtual system using Equation (29), the virtual system's response x d , can be found.
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.
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 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.

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 (Lu r and Lw r ) are calculated. The Lu and Lw 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 Lu and Lw 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 (Lu r and Lw r ) models and the actual Lu and Lw models are compared.

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 statespace model obtained by the output error method; those applications and results are explained in past studies [23]. By using the model horizon parameters = 30, = 10, the subspace estimation model parameters ( and ) are calculated. The and obtained with the same values as the previous and from the system are then calculated using the input signal to be tested, the and are then calculated using the same values as the previous and from the system ( Figure 5. In the graphics, the results calculated using the reference ( and ) models and the actual and models are compared. 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. 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 sinω 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.

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 statespace model obtained by the output error method; those applications and results are explained in past studies [23]. By using the model horizon parameters = 30, = 10, the subspace estimation model parameters ( and ) are calculated. The and obtained with the same values as the previous and from the system are then calculated using the input signal to be tested, the and are then calculated using the same values as the previous and from the system ( Figure 5. In the graphics, the results calculated using the reference ( and ) models and the actual and models are compared. 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. 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.
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 sinω 0 t + A 1 sinω 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.
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.3 rd/s, where the input signal runs throughout the full stroke in a single alternate, and = 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 ( = + ). The experiment that was performed using input signals with = 0.2 rd/s and = 2, = 0.8 rd/s and = 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. 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. Table 2. Experiments with a two-frequency u input signal (u = A 0 sinω 0 t + A 1 sinω 1 t). The following experiments (13-16) use u input signals that have three distinct frequency components (u = A 0 sinω 0 t + A 1 sinω 1 t + A 2 sinω 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 experiments (13-16) use u input signals that have three distinct frequency components ( = + + ), 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.  Table 3. Experiments with a three-frequency u input signal (u = A 0 sinω 0 t + A 1 sinω 1 t + A 2 sinω 2 t). 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. The following experiments (13-16) use u input signals that have three distinct frequency components ( = + + ), 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.

Experimental Results of Data-Driven Predictive Controller
The predicted model parameters and control parameters Ke and K∆wp 1x2p 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. 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 Nc = 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.

Experimental Results of Data-Driven Predictive Controller
The predicted model parameters and control parameters and Δ 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. 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, = 30, = 10, = 1, = 2, and = 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. Step response according to u signals of different N data lengths.

Analysis of Model Horizon Parameters (p, f)
The system state equations along the chosen horizon are the basis of the subspacebased 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 = 10, = 1, = 2, and = 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. Reference N=3000 N=1500 N=500 Figure 10.
Step response according to u signals of different N data lengths.

Analysis of Model Horizon Parameters (p, f )
The system state equations along the chosen horizon are the basis of the subspacebased 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 Nc = 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. 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 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).  Table 6, each of these values is described in detail.  In the following experiments, the past model horizon p, and other parameters (p = 30, Q = 1, R = 2, and Nc = 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∆wp 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.

Δ
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 Δ parameter calculated in successful results. Figures 13 and 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.  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∆wp parameter calculated in successful results. Figures 13 and 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.

Δ
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 Δ parameter calculated in successful results. Figures 13 and 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.

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

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 Nc , which is the first component of ∆u Nc 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.  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 should be calculated. It is used by the control algorithm to compute the system control signal , which is the first component of Δ along Nc. These tests were conducted with the following parameters: = 30, = 10, = 1, and = 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.  Figure 16. Effect of the Q parameter on step and trajectory response.

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 ′ = . When a maximum counterforce of 7-8 N was applied in the final 10th second, the virtual system response deviates from the reference ′. 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).

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).
. When a maximum counterforce of 7-8 N was applied in the final 10th second, the virtual system response deviates from the reference ′. 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 in the 12th 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 in negative. 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 -6N in the 12th 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. 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.

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.  Institutional Review Board Statement: Ethical review and approval were waived for this study due to no applicable.