The Design and Adaptive Control of a Parallel Chambered Pneumatic Muscle-Driven Soft Hand Robot for Grasping Rehabilitation

Abstract

The widespread application of exoskeletons driven by soft actuators in motion assistance and medical rehabilitation has proven effective for patients who struggle with precise object grasping and suffer from insufficient hand strength due to strokes or other conditions. Repetitive passive flexion/extension exercises and active grasp training are known to aid in the restoration of motor nerve function. However, conventional pneumatic artificial muscles (PAMs) used for hand rehabilitation typically allow for bending in only one direction, thereby limiting multi-degree-of-freedom movements. Moreover, establishing precise models for PAMs is challenging, making accurate control difficult to achieve. To address these challenges, we explored the design and fabrication of a bidirectionally bending PAM. The design parameters were optimized based on actual rehabilitation needs and a finite element analysis. Additionally, a dynamic model for the PAM was established using elastic strain energy and the Lagrange equation. Building on this, an adaptive position control method employing a radial basis function neural network, optimized for parameters and hidden layer nodes, was developed to enhance the accuracy of these soft PAMs in assisting patients with hand grasping. Finally, a wearable soft hand rehabilitation exoskeleton was designed, offering two modes, passive training and active grasp, aimed at helping patients regain their grasp ability.

1. Introduction

There has been an increase in patients with hand dysfunction caused by diseases such as hemiplegia and strokes [1]. According to statistics, around 13.68 million people suffer new strokes globally each year, with about 70% of survivors experiencing some degree of upper limb and hand motor function impairment [2]. About 67% of stroke complication patients are still unable to use their hands even four years after onset. Key goals of hand rehabilitation include grasping objects, transferring, manipulating, coordinating finger movements, and maintaining flexibility [3]. Repetitive hand flexion and extension training is required to restore motor nerve function [4,5]. In clinical practice, rehabilitation is primarily conducted by physicians, and while effective, this approach has low efficiency. Consequently, researchers globally have extensively investigated the design and control of hand rehabilitation robots. Rigid structure hand exoskeletons [6,7], mainly driven by linkages, motors, and ropes, and underdriven structure hand rehabilitation robots [8] using spring drives offer precise control but have poor wearability. These rigid or semi-rigid devices can cause secondary injury during improper operation.

Soft robots driven by pneumatic artificial muscles (PAMs), made from soft materials such as silicone rubber and fabric, produce movements like bending, elongation, and twisting through changes in air pressure. This significantly enhances wearability and safety compared to rigid robots [9]. Soft rehabilitation robots offer unparalleled advantages in human–robot interaction, rehabilitation training, and assisted grasping.

Research has shown that PAM actuators, using pneumatic pressure as a power source, have many advantages. They are a central focus and challenge in the design and control of soft rehabilitation robots and represent an important trend for future development [10]. In recent years, soft PAMs that fit the body have been applied in soft rehabilitation exoskeletons. A typical bending PAM is made of silicone rubber airbags. Asymmetric airbags or limiting layers with fibers result in varying degrees of elongation on both sides of the PAM, enabling unidirectional bending [11]. Han et al. [12] designed a fiber-reinforced soft pneumatic bionic actuator that bends by utilizing the strain difference between the chamber and the limiting layer. Wang Tianyu et al. [13] from Shanghai Jiao Tong University proposed a pneumatic grid torsion actuator with a programmable grid shape design. Liu Chih-Hsing et al. [14] developed a soft pneumatic bending actuator using a grid-type structure. Their hand claw comprises three grid-type bending actuators connected in parallel by end actuators. Panagiotis P. et al. [15] utilized the enhancement effect of fiber threads in different directions to design a new bending PAM. Integrating bending PAMs with soft fabric gloves to create a hand rehabilitation exoskeleton minimizes damage during repetitive rehabilitation training. This approach is crucial for improving rehabilitation efficiency in clinical settings, enabling remote rehabilitation, and reducing the repetitive workload of physicians [16].

However, existing research on soft hand rehabilitation devices faces several problems, such as a lack of soft multi-degree-of-freedom PAMs and an inability to create lightweight PAMs that adapt to the natural range of motion of fingers of different lengths. Consequently, many scholars have turned their attention to designing novel cavity and braiding structures, optimizing PAM structural parameters, and modeling PAM kinematics through finite element simulation.

Due to the characteristics of soft materials, such as their nonlinearity and complex fiber-constrained structures, realizing precise position control of hand rehabilitation exoskeletons driven by PAMs is a significant research direction [17]. Many researchers use simple PID control methods, but the control accuracy and speed are insufficient [18]. Some researchers delve into the dynamic characteristics of PAM to establish its dynamic model [19]. For the modeling of soft PAM actuators, empirical model building, geometric model analysis, and finite element analysis methods are commonly used [20]. Wang Boran et al. [21] proposed a way to establish the relationship between the compilation angle and output force using finite element simulation to analyze the PAM. For the very complex nature of the dynamics of fluid soft actuators, Wang Tao et al. [22] used a second-order transfer function to describe the kinetic behavior of the motion from the driving air pressure to the bending angle of the soft actuator. Given the advantages of neural networks in handling nonlinear problems, applying neural networks to PAM control can effectively reduce tracking errors and achieve better position control and trajectory tracking. Elgeneidy. K et al. [22] proposed a purely data-driven PAM modeling method for controlling the bending angle of the PAM actuator, and a data-driven modeling technique using a regression analysis and artificial neural networks was used to derive an empirical model based on the generated experimental data.

In this paper, we design a soft wearable hand rehabilitation exoskeleton for stroke patients. We develop parallel chamber bidirectional bendable PAM actuators, analyze the static bending characteristics using a finite element analysis, and optimize their structural parameters. The reduced order forms of aerodynamic muscle strain energy, generalized force, kinetic energy, and gravitational potential energy were substituted into the Lagrangian equation to obtain the dynamic equation of PAM, and considering the need for precise position control of PAMs during rehabilitation grasp training, an optimized radial basis function neural network adaptive control algorithm (RBFNNO) based on parameter and hidden layer node optimization is designed. Finally, we design a hand rehabilitation exoskeleton that uses PAMs to drive the patient’s fingers for grasp training. The rest of the paper is organized as follows: Section 2 describes the design and fabrication of the parallel chamber PAMs and soft exoskeleton. Section 3 presents the PAM dynamic modeling and control methods. Section 4 details the hand rehabilitation experiments. Section 5 contains the discussion, followed by the conclusion in Section 6.

2. Bio-Inspired Mechanical Design

Strokes, hand trauma, and nerve injuries lead to hand dysfunction, significantly impacting patients’ daily lives. The hand is a highly complex mechanism with over 20 degrees of freedom (DOFs). The DOFs of the hand joints are illustrated in Figure 1. In patients with hand dysfunction, muscle contraction and impaired joint movement reduce the DOFs and range of motion.

Traditional hand rehabilitation primarily relies on one-on-one repetitive training by physicians, but the low efficiency and insufficient training volume can result in suboptimal rehabilitation outcomes. Hand rehabilitation exoskeletons are clearly beneficial in reducing repetitive labor and enhancing rehabilitation quality. The hand exoskeleton is designed to drive the fingers to bend at specific angles and transmit appropriate forces to the fingers. To design a hand exoskeleton with a high fit, flexibility, and dexterity, the bidirectional bendable parallel chamber PAM is selected as the actuator to drive finger movement through muscle flexion.

Considering the requirements for daily hand rehabilitation, this paper designed a parallel chamber PAM actuator by comparing the fiber-reinforced actuators developed by current researchers. Inspired by the motion of hand tendons, a parallel chamber structure was used in the design of the bending PAM to accommodate the needs of finger flexion and extension movements, as shown in Figure 2. The overall design of the actuator is cylindrical to ensure high power density and fast response. To meet the requirements for effective hand rehabilitation, the actual bending angle of the PAM should exceed 100°, and the fingertip force should be greater than 0.8 N.

The soft PAM consists of a silicone rubber shell and two independent asymmetric chambers. To limit the radial expansion and torsional bending of the PAM, Kevlar fiber lines are wrapped around the surface in a double-helix pattern, with the intersection of the helix ideally aligned along the axis of the chamber. The number of turns of the fiber lines affects the actual bending performance of the PAM. Too many turns increase stiffness, requiring higher inflation pressure to achieve the same bending effect, which reduces the PAM’s power density. Conversely, too few turns do not adequately constrain the PAM’s expansion, leading to excessive radial expansion that severely impacts bending performance.

Since the angles required for finger flexion and extension are different, Chambers A and B were designed with different radii. Chamber A is larger than Chamber B. When Chamber A is inflated, it produces significant bending to assist finger flexion, while Chamber B generates less bending to aid finger extension. Compared to single-chamber PAMs, our design of a double-chamber PAM enhances the extension function, providing better assistance for patients with muscle weakness.

In order to verify the actual bending effect, the simulation analysis of the PAM bending angle was carried out using ABAQUS finite element software (Version 6.14). The analysis of the soft body drive model using ABAQUS software mainly includes creating parts, setting material properties, and dividing the mesh, assembly, boundary conditions, and loads.

The Yeoh model can be simply used for the deformation behavior of nonlinear materials, and the model parameters can be well fitted by uniaxial tensile experiments, and the Yeoh model strain density function model is as follows.

$$W={\sum }_{i=1}^N{C}_{10}(I_1-3{)}^i+{\sum }_{k=1}^N{\displaystyle \frac{1}{d_k}}(J-1{)}^{2k}$$

(1)

$$W=C_{10}(I_1-3)+C_{20}(I_1-3{)}^2$$

(2)

where N, $C_{i0}$, and $d_k$ are material constants determined by uniaxial tensile tests. The initial shear modulus $\mu =2C_{10}$.

The control variable method was used to simulate and compare various factors affecting the bending performance of PAMs in finite element simulations, as shown in Figure 3. This approach helps in selecting optimal parameters that meet finger rehabilitation requirements. The specific simulation results are detailed in Section 4.2. The optimized design parameters are listed in

Table 1. PAM parameters.

Parameters	Numerical Value
Drive length/mm	100/120
Drive radius/mm	8
Inner diameter of cavity A/mm	3.5
Inner diameter of cavity B/mm	3.0
Distance of cavity A center deviation from drive center/mm	3.5
Distance of cavity B center deviation from drive center/mm	3.0
Number of spiral wire windings N	34/42
Winding pitch/mm	3
Drive mass/g	25/31

. A double-helix wire with a pitch of 3 mm was chosen for winding the PAM body, which has lengths of 100 mm and 120 mm. To prevent the fiber-wound wire from being exposed outside the silicone rubber housing, an additional 1 mm silicone rubber sleeve was added to cover the exterior of the PAM.

We employ silicone casting to produce PAMs, utilizing a 1:1 mixture of Dragon-skin 30 silicone. The mold was designed using SOLIDWORKS based on specific design parameters and fabricated with a Stratasys F170 3D printer. The manufacturing process for the PAM, as illustrated in Figure 4, involves several key steps: mold design and production, mixing the silicone solution, pouring the matrix with cavities, attaching fiber threads, pouring the outer silicone, and assembling the end caps.

Velcro is sewn onto the fabric glove. Velcro is used to secure the PAMs. The PAMs were secured to the distal interphalangeal (DIP) and metacarpophalangeal (MCP) joints with straps, fitting well with the patient’s fingers, as shown in Figure 5. To guarantee the durability of the PAM, each actuator is tested multiple times to confirm proper operation within a specific pressure range.

The bending angle range of the five fingers of the hand exoskeleton is determined by the PAMs. The thumb and little finger use 100 mm PAMs, while the index, middle, and ring fingers use 120 mm PAMs. As shown in

Table 2. Actual range of motion of PAMs.

Finger	Flexion Angle	Extension Angle
Thumb	0–89°	0–12°
Index finger	0–118°	0–21°
Middle finger	0–120°	0–18°
Ring finger	0–119°	0–20°
Little finger	0–91°	0–11°

, the five PAMs made for exoskeletons have different bending ranges. Patients with hand injuries cannot bend their fingers normally, making it difficult to reach the maximum bending angle of healthy fingers. It is also difficult to stretch after bending the fingers. The PAMs drive the movement of the patients’ fingers, enabling them to complete rehabilitation training.

The hand exoskeleton assisted the patient in performing various rehabilitation exercises. At different rehabilitation stages, the hand exoskeleton offers various modes to optimize rehabilitation outcomes, aiding in nerve recovery and enhancing hand movement abilities. In the early stage of rehabilitation, a passive training mode involving repetitive flexion and extension was designed for patients. The hand exoskeleton drove the patient’s fingers to complete these passive training tasks. As the patient’s hand strength and nervous system had recovered to some extent, an active grasping rehabilitation training mode based on the mapping of the healthy and affected side was implemented. The angular signal of the patient’s healthy hand during object grasping was collected by a motion capture system and used to control the exoskeleton to assist in grasping objects.

3. Methods

3.1. Pneumatic Muscle Dynamic Modeling

The PAM is primarily composed of silicone rubber and fiber threads. The relationship between the bending angle and the input air pressure can be established by analyzing a geometric model. The model is based on the following assumptions [23]: the silicone rubber material is incompressible when charged with air pressure, meaning the total volume of the silicone rubber remains constant; the actuator remains cylindrical during inflation and bending [24]; the actuator maintains a circular shape during bending and exhibits an overall constant curvature change; and the fiber helix restricts radial extension and is always in contact with the outer surface of the PAM.

The dynamic model of PAM is shown in Figure 6. The nonlinear terms in the strain energy density function and the sine and cosine terms in the robot dynamics are approximated using Taylor series. Lagrangian dynamics are an element of analytical mechanics that combines work, energy, and generalized coordinates to reformulate Newtonian mechanics. Unlike Newtonian mechanics, Lagrangian dynamics do not rely on a spatial coordinate system and do not require analysis of the binding forces within the system.

The cross-sectional and side views of the PAMs are shown in Figure 7. The bending radius R, eccentricity of cavity A $e_1$, eccentricity of cavity B $e_1$, and bending angle $\theta$ are selected to describe the bending degree of the actuator under varying air pressures. The radius $R_0$ represents the original radius of the pneumatic muscle, while $r_{10}$ and $r_{20}$ denote the original radii of cavities A and B, respectively. $L_1$ and $L_2$ represent the lengths of the central axes of cavities A and B after the pneumatic muscle inflates and deforms.

The strain energy density function of the Yeoh model is used to calculate the elastic strain energy, and the strain energy density function can be expressed as follows:

$$W=C_{10}\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2-3\right)$$

(3)

In the above equation, ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are the axial, radial, and circumferential principal stretch ratios of the PAM, respectively. $C_{10}$ is the Yeoh model material parameter.

When a certain air pressure is filled in cavity A, the dimensional variation of the driver cavity in three directions—axial, radial, and circumferential—can be obtained according to the driver wall thickness parameters [25]:

$$\left\{\begin{array}{c}{\lambda }_1={\displaystyle \frac{l_1}{L_0}}\\ {\lambda }_2={\displaystyle \frac{r_1}{R_0}}\\ {\lambda }_3=({\displaystyle \frac{t_1}{t_{10}}}+{\displaystyle \frac{t_2}{t_{20}}}+{\displaystyle \frac{t_3}{t_{30}}}+{\displaystyle \frac{t_4}{t_{40}}})/4\end{array}\right.$$

(4)

where $L_0$ represents the initial length of the central axis and $R_0$ denotes the initial radius of the PAM when the PAM is not deformed.

The thickness of the front and rear sides of cavity A are $t_{10}$ and $t_{30}$. The thickness of the left and right sides of cavity A are $t_{20}$ and $t_{40}$. After inflating cavity A, the thicknesses $t_{10}$/$t_{20}$/$t_{30}$/$t_{40}$ become $t_1$/$t_2$/$t_3$/$t_4$.

The PAM inflation deformation resembles a circular arc with the center of circle O, which can be obtained by geometric structure analysis:

$$\left\{\begin{array}{c}{t}_1=R-e_1-r_1\\ t_2=\sqrt{R^2-{e_1}^2}-r_1\\ t_3=R+e_1-r_1\\ t_4=\sqrt{R^2-{e_1}^2}-r_1\end{array}\right.$$

(5)

By the binding of the radial fiber threads, the driver is barely deformed in the radial direction. It can be assumed that the radial stretch ratio ${\lambda }_2$ is 1. The axial stretch of the silicone rubber varies with the height in the cross section, and the axial stretch ratio ${\lambda }_1$ can be set as shown in the following equation, where $\beta$ is the variation of the thickness along the bottom to the top.

$${\lambda }_1=1+{\displaystyle \frac{\theta }{L}}\beta$$

(6)

$${\lambda }_3={{\lambda }_1}^{-1}$$

(7)

The strain energy density is therefore rewritten as the following equation, where G is the shear modulus of the silicone rubber material.

$${\rho }_w={\displaystyle \frac{G}{2}}\left({{\lambda }_1}^2+{{\lambda }_1}^{-2}-2\right)$$

(8)

The strain energy $E_w$ is expressed as integrating the entire silicone rubber, ignoring the uninflated cavity B.

$$E_w=2{\int }_{-{\displaystyle \frac{\pi }{2}}}^{{\displaystyle \frac{\pi }{2}}}d\phi {\int }_0^{R_0}{E}_{w}Ltdt-2{\int }_0^{{\displaystyle \frac{\pi }{2}}}d\phi {\int }_0^{2e_{1}sin\phi }{E}_{w}Ltdt$$

(9)

The gravitational potential energy $E_g$ is mainly determined by its average relative height in the vertical direction, and the gravitational potential energy can be expressed as follows:

$$E_g=-mg{\displaystyle \frac{{\int }_0^{\theta }R\mathit{sin}\psi d\psi }{\theta }}=-{\displaystyle \frac{mgL}{{\theta }^2}}\left(1-\mathit{cos}\theta \right)$$

(10)

where $m$ is the PAM mass and $g$ is the gravitational acceleration. The analysis of kinetic energy reduces the drive to a curve, and the coordinates of each point on the drive can be expressed in the following equations:

$$x(s)={\displaystyle \frac{L}{\theta }}(1-\mathit{cos}{\displaystyle \frac{s\theta }{L}}))$$

(11)

$$y\left(s\right)={\displaystyle \frac{L}{\theta }}\mathit{sin}{\displaystyle \frac{s\theta }{L}}$$

(12)

where $s$ is the arc length between the point on the PAM and the origin.

The total kinetic energy $E_k$ can be calculated by the following equation:

$$E_k={\displaystyle \frac{m}{2L}}{\int }_0^L\left[{\left({\displaystyle \frac{dx}{dt}}\right)}^2+{\left({\displaystyle \frac{dy}{dt}}\right)}^2\right]ds=m{L}^2{\theta }^2\left[{\displaystyle \frac{1}{{6\theta }^2}}+{\displaystyle \frac{1}{{\theta }^4}}\left(1+\cos \theta \right)-{\displaystyle \frac{2}{{\theta }^5}}\sin \theta \right]$$

(13)

According to the principle of virtual work, the generalized force caused by the inflation pressure can be expressed as follows:

$$f_p=u{\displaystyle \frac{\partial V_h\left(\theta \right)}{\partial \theta }}=\mathsf{\pi }\left[2r_1\left(L_0+e_1\theta \right)·{\displaystyle \frac{d{r}_1}{d\theta }}+r_1^2·e_1\right]u=k_{p}u$$

(14)

where $u$ is the air pressure, $V_h\left(\theta \right)$ is the volume of the cavity A after inflation and deformation, and $k_p$ is the coefficient used to simplify the expression.

Based on the potential and kinetic energies calculated above, the Lagrangian function is as follows:

$$L_f=E_k-E_g-E_w$$

(15)

The Lagrangian kinetic equations are as follows:

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial L_f}{\partial \dot{\theta }}}-{\displaystyle \frac{\partial L_f}{\partial \theta }}=f_p-\xi \dot{\theta }$$

(16)

where $\xi$ is the coefficient used to account for the damping characteristics and viscous resistance.

Using series expansion ${(1+k\theta )}^{-2}=1-2k\theta +3{\left(k\theta \right)}^2-4{\left(k\theta \right)}^3$ for Equation (8), we can obtain the following:

$${\rho }_w={\displaystyle \frac{G}{2}}[1+4{\displaystyle \frac{{\theta }^2}{L^2}}{\beta }^2-4{\displaystyle \frac{{\theta }^3}{L^3}}{\beta }^3]$$

(17)

So, Equation (9) can be simplified as

$$\begin{array}{c}{E}_w=G{\int }_{-{\displaystyle \frac{\pi }{2}}}^{{\displaystyle \frac{\pi }{2}}}d\phi {\int }_0^{R_0}[1+4{\displaystyle \frac{{\theta }^2}{L}}{t}^3{sin}^{2}t-4{\displaystyle \frac{{\theta }^3}{L^2}}{t}^4{sin}^{3}t]dt\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\int }_0^{{\displaystyle \frac{\pi }{2}}}d\phi {\int }_0^{2e_{1}sin\phi }{\int }_0^{R_0}[L+4{\displaystyle \frac{{\theta }^2}{L}}{t}^3{sin}^2\phi -4{\displaystyle \frac{{\theta }^3}{L^2}}{t}^4{sin}^3\phi ]dt=k_1{\theta }^2-k_2{\theta }^2-k_3{\theta }^3\hfill \end{array}$$

(18)

The sine and cosine functions in Equations (10) and (13) use the series expansion and ignore the terms above the sixth order in the cosine function and above the seventh order in the sine function as follows:

$$E_g=mgL\left({\displaystyle \frac{{\theta }^2}{24}}-{\displaystyle \frac{1}{2}}\right)$$

(19)

$$E_k=\left({\displaystyle \frac{1}{40}}-{\displaystyle \frac{{\theta }^2}{1008}}\right)m{L}^2{\dot{\theta }}^2$$

(20)

By substituting Equations (14) and (18)–(20) into Equation (16), the dynamic model can be obtained as follows:

$$\left[\left({\displaystyle \frac{1}{20}}-{\displaystyle \frac{{\theta }^2}{504}}\right)m{L}^2\right]\ddot{\theta }+\left[\zeta -{\displaystyle \frac{\theta }{252}}m{L}^2\dot{\theta }\right]\dot{\theta }+\left[2\left(k_1-k_2\right)-3k_3\theta +{\displaystyle \frac{mgL}{12}}\right]\theta =k_{p}u$$

(21)

When the derivative term is neglected, the equation simplifies to a static model.

3.2. RBFNN Adaptive Control

3.2.1. Optimization Methods for RBFNN Nodes and Parameters

The three-layer structure of the RBFNN is shown in Figure 8, including the input layer, hidden layer, and output layer. The input is the data source node. The transformation from the input layer to the hidden layer is achieved using a nonlinear Gaussian basis function, while the transformation from the hidden layer to the output layer is linear. The output of the hidden layer is represented as a quadratic term through the weights.

The choice of the radial basis function in an RBFNN determines its mapping ability. The Gaussian kernel function is the most used radial basis function, as shown in the following equation:

$$h_j(t)=\mathit{exp}(-{\displaystyle \frac{{‖x\left(t\right)-c_j\left(t\right)‖}^2}{2b_j^2}}), j=1,\dots ,m$$

(22)

In Equation (20), $c_j\left(t\right)$ and $b_j$ represent the center and width of the kernel function, respectively, and are calculated using an optimization method based on an improved k-means clustering algorithm.

The desired trajectory angle $\theta$ and angular velocity $\dot{\theta }$ of the PAM are sampled to form a sample set $D=\left\{x_1,x_2,\cdots ,x_m\right\}$. A random sample from $D$ is selected as the centroid ${\mu }_1$. Additional centroids are selected based on the following criteria. When selecting a new centroid ${\mu }_{r+1}$, the minimum distance between the samples ($x_j)$ and the selected centroids (${\mu }_1,{\mu }_2\cdots {\mu }_r$) is calculated as follows:

$$D\left(x_j\right)=\arg min{‖x_j-{\mu }_{1,2\cdots ,r}‖}^2$$

(23)

The larger the value of $D\left(x_i\right)$ is, the higher the probability that the sample will be selected until k centroids are selected is. Next, the distance $d_{ji}={‖x_j-{\mu }_i‖}_2$ is calculated from each sample $x_j$ to each centroid.

$$d_{ji}={‖x_j-{\mu }_i‖}_2$$

(24)

The variable $x_j$ is assigned to cluster $C_i$, corresponding to the nearest centroid. Finally, the centroid of each cluster$C_i$ is computed to update the original centroids:

$${\mu }_i=\left({\displaystyle \frac{1}{\left|C_i\right|}}\right)·\sum x,x\in C_i$$

(25)

The final centroids ${\mu }_{1,\cdots ,k}$ are the centers $c_j\left(t\right)$ of the kernel functions.

Next, the average distance between each centroid and the other centroids is calculated as the basis for width selection:

$$meanD\left({\mu }_j\right)={\displaystyle \frac{\sum _{i\ne j}dis\left({{\mu }_i,\mu }_j\right)}{k-1}}$$

(26)

For the cluster $C_j$ belonging to the same cluster center ${\mu }_j$, the variance and scaling factor are calculated as follows:

$$S_j={\displaystyle \frac{1}{size\left(C_j\right)}}\sum _{x_i\in C_j}dis{\left(x_i,{\mu }_j\right)}^2$$

(27)

$${\epsilon }_j={\displaystyle \frac{S_j}{\frac{1}{k}\sum _{j=1}^k{S}_j}}$$

(28)

The width corresponding to each center is then given by the following:

$$b_j={\epsilon }_j·meanD\left({\mu }_j\right)$$

(29)

For a specific approximation error, the distribution of hidden nodes optimized by the improved k-means algorithm and the variance-based width selection scheme will significantly reduce the number of hidden nodes and the complexity of the RBF network structure. The advantages of the optimized allocation scheme for hidden layer nodes can be understood in the following ways: all hidden layer nodes are effective for the approximation task, and unnecessary hidden nodes disappear, thus significantly reducing the number of hidden nodes, especially for high-dimensional input RBFs; each hidden layer node works at its maximum approximation capacity.

3.2.2. Control Law

Equation (21) can be simplified as a static model of PAMs as follows:

$$M\left(\theta \right)\ddot{\theta }+C\left(\theta ,\dot{\theta }\right)\dot{\theta }+G\left(\theta \right)=k_{p}u$$

(30)

where $\theta$,$\dot{\theta }$, and $\ddot{\theta }$ are the vectors of PAM bending angle, angular velocity, and angular acceleration, respectively, and $u$ is the input air pressure vector. This equation is rewritten as a second-order nonlinear system of bending PAMs.

$$\ddot{\theta }=-{\displaystyle \frac{1}{M\left(\theta \right)}}\left[C\left(\theta ,\dot{\theta }\right)\dot{\theta }+G\left(\theta \right)\right]+{\displaystyle \frac{k_p}{M\left(\theta \right)}}u$$

(31)

State variables are used to represent $x_1=\theta$, $x_2={\dot{x}}_1$ and $y=x_1$.

$$\left\{\begin{array}{c}{x}_2={\dot{x}}_1\hfill \\ {\dot{x}}_2=f(x_1,x_2)+g(x_1,x_2)u\hfill \\ y=x_1\hfill \end{array}\right.$$

(32)

where $f\left(\theta ,\dot{\theta }\right)$ and $g\left(\theta ,\dot{\theta }\right)$ are nonlinear functions and $u\in R^n$ and $y\in R^n$ are the control inputs and outputs of the system, respectively.

The desired tracking trajectory is $y_d$; then, the error can be expressed as follows:

$$e=y_d-y=y_d-x_1,E=\left(e,\dot{e}\right)$$

(33)

Substituting the above equations into the second-order nonlinear system yields the error system:

$$\ddot{e}+k_{p}e+k_d=0$$

(34)

Let $K={\left(k_P,k_d\right)}^T$ so that the roots of polynomial $s^2+k_{d}s+k_p=0$ are all in the left half of the complex plane. Then when $t\to \infty$, $e\left(t\right)\to \infty$ and $\dot{e}\left(t\right)\to \infty$.

The RBFNN is employed to approximate the nonlinear function:

$$f=W^{T}h\left(x\right)+\epsilon$$

(35)

Replacing the nonlinear function in Equation (29) with the RBFNN output gives the control law:

$$\left\{\begin{array}{c}u={\displaystyle \frac{1}{g\left(x\right)}}\left[-\stackrel{̑}{f}\left(x\right)+{\ddot{y}}_d+K^{T}E\right]\\ \stackrel{̑}{f}\left(x\right)={\stackrel{̑}{W}}^{T}h\left(x\right)\end{array}\right.$$

(36)

The constant $\gamma$ and constant matrix $Q$ are the parameters that need to be set in the experiment. The designed adaptive law for weight updating can be expressed as follows:

$$\dot{\widehat{W}}=-\gamma E^{T}PBh\left(x\right)$$

(37)

The RBFNN adaptive control is illustrated in Figure 9. The center $c_j$ and width $b_j$ of the Gaussian function in the neural network are determined using an improved k-means algorithm, while the weights $\widehat{W}$ are updated in real time according to an adaptive law. This allows the neural network to continuously approximate the unknown nonlinear components $\stackrel{̑}{f}\left(x\right)$ of the PAM dynamic model, ultimately computing the output $u$ based on the control law. Stability analysis can be found in Appendix A.

4. Experiments and Results

4.1. Experiment Setup

The hand exoskeleton system, as shown in Figure 10, consisted of an upper computer, a data acquisition module, an air pump, proportional valves, bending angle sensors, and the hand exoskeleton. To accommodate varying finger lengths, flexible angle sensors, such as the customizable version of the Sparkfun FLEX4.5, were used. The data acquisition module (Product Model: USB-6210; Company: National Instruments) provided real-time feedback of bending angle signals to the upper computer. The upper computer processed these data and output corresponding voltage signals through the data acquisition module to control the proportional valve. The proportional valve adjusted the air pressure in the PAM cavity according to changes in input voltage, resulting in the PAMs bending to different angles.

We conducted three sets of experiments. The first set involved using finite element simulation software, Abaqus, to analyze the impact of various factors such as wall thickness, length, radius, and cavity shape on the performance of PAMs. Based on the optimal parameters identified from these simulations, PAMs were fabricated, and their bending characteristics were compared between physical experiments and simulations. The second set of experiments focused on controlling a single PAM, comparing the proposed RBFNNO control algorithm with other algorithms to demonstrate its effectiveness. The third set of experiments tested the hand exoskeleton in both passive training mode and active grasping mode. In passive training mode, the PAMs assisted the fingers in moving along a specified trajectory. In active grasping mode, the motion trajectory of the healthy hand was captured using a motion capture system, and the exoskeleton worn on the affected hand was controlled to perform the corresponding grasping task.

To prevent data loss due to occlusion during movement, we selected reflective markers with a diameter of 3 mm, distributing them evenly across the joints of the hand. The placement of the hand markers was shown in Figure 11. Each of the five fingers had three markers attached at the MCP, proximal interphalangeal (PIP), and DIP joints to calculate finger bending angles. A reflective marker was placed on the wrist to serve as the origin of the spatial coordinate system. MATLAB was used to process the coordinate changes of each marked point and convert them into joint angles

4.2. Simulation and Physical Testing of PAM

To simulate the flexion and extension movements of fingers, PAMs operate in two modes: flexion and extension. Several factors, including the length of the parallel chamber actuator, cavity shape, eccentric structure of the pneumatic chamber, cavity radius, and helix line, influence the bending deformation effect of these pneumatic muscles. This section explores the factors affecting the bending and extension deformations of PAMs through a simulation analysis.

As shown in Figure 12a,b, the best performance is achieved with a semicircular cross-section. However, it is challenging to design a suitable parallel chamber structure, so we opted for a circular cavity. Figure 12c,d indicated that longer PAMs result in greater bending and extension angles. Considering the length of the fingers, we chose 120 mm long PAMs to assist the index, middle, and ring fingers, and 100 mm long PAMs for the other fingers. Figure 12e,f showed that larger cavity radii can enhance bending ability under high pressure. To maintain a sufficient wall thickness and prevent rupture during inflation, we set the radius to 3–3.5 mm based on finger size. Figure 12g,h demonstrated that both larger and smaller spiral angles of the fiber lines reduce bending performance, especially at higher pressures. Therefore, we selected a spiral angle of 2° for the bi-directional winding of the fiber lines.

Based on the analysis of the simulation results, the parameters affecting the bending of the PAM, ranked from largest to smallest influence, are as follows: wall thickness, length, chamber diameter, cavity shape, and center diameter. To achieve an optimally bending PAM, geometric parameters such as longer length, larger chamber diameter, smaller center diameter, and thinner wall thickness are preferred. However, considering actual rehabilitation needs, manufacturing difficulty, durability, large pressure tolerance, and power density, we selected two PAMs with a length of 100/120 mm, a chamber diameter of 3.5/3.0 mm, a center diameter of 8 mm, and a wall thickness of 1.5 mm for physical fabrication.

To verify the effectiveness, we performed a finite element analysis on the flexion and extension movements of the PAMs. The simulation applied a load to one of the parallel chambers, using pressure settings from 0 to 180 kPa in increments of 20 kPa. The displacement contour plots of the PAM in ABAQUS for this state were shown in Figure 13. The first seven images illustrated the application of air pressure from 0 to 180 kPa in Chamber A, corresponding to finger flexion. The last two images showed air pressure applied from 30 kPa to 60 kPa in Chamber B, corresponding to finger extension.

Next, we manufactured the PAM based on the selected parameters. Inflation pressure was incrementally increased by 20 kPa, ranging from 0 kPa to 180 kPa. Once the PAM reached the maximum inflation pressure and stabilized, the steady-state bending angle and fingertip force were measured. Ten steady-state tests were performed at each pressure level. To ensure accuracy, the data were averaged after excluding the maximum and minimum values.

The experimental test results shown in Figure 14 closely match the static model simulations of the PAM. Some errors in the bending angle curve during the actual test were due to the exclusion of the gravity term in the simulation and modeling. The maximum forward bending angle reached 125°, while the reverse bending angle reached 50°. The tests confirmed that the PAM’s bending angles exceeded the normal range for finger bending and extension. Therefore, PAMs can be integrated with soft fabric gloves to assist patients in repetitive grasp training.

4.3. PAM Angle Tracking Control

To validate the effectiveness of the proposed control algorithm, both simulations and experimental tests were conducted on the control of a single PAM. The desired trajectory was set to $y_d=0.5\mathit{sin}\pi t+0.5$. The initial state was $x_1=0$, $x_2=0$. Parameters were set as $Q=\left[\begin{array}{cc}500& 0\\ 0& 500\end{array}\right]$, $k_d=50$, $k_p=30$, and the adaptive parameter $\gamma$= 1200. A duration of 2000 s was selected for trajectory tracking. The first 20 s (initial state) and the last 20 s (stable tracking) were analyzed separately, with a sampling interval of 0.01 s.

The distribution of RBFNN hidden layer nodes was optimized using the improved k-means algorithm. The number of clustering centers, which corresponded to the hidden layer nodes of the RBFNN, was set to 5, 10, 15, and 20. For each cluster center, a variance-based width optimization algorithm was used to calculate the width of each radial basis function. The initial weights of the network were set to zero and were iteratively updated using an adaptive rate. The parameters $c_i$ were determined based on the clustering centers selected by the algorithm.

As shown in Figure 15, when tracking the desired trajectory, the RBFNN becomes more complex with an increasing number of hidden layer nodes, resulting in improved tracking of the PAM trajectory. The PAMs start at rest, leading to a significant initial tracking error. Due to the nonlinearity and hysteresis, the tracking error is most pronounced in the inflation deflation transition stage.

The more hidden layer nodes there are in the initial state, the better the trajectory tracking effect is, and the number of hidden layer nodes accelerates the initial tracking effect. However, as shown in Figure 16, after achieving stable tracking, the error of the RBFNN with 20 hidden layer nodes increases and becomes larger than that of the RBFNN with 10 and 15 hidden layer nodes. The results indicate that reducing the number of hidden layer nodes not only reduces the complexity of the neural network but also improves stability.

To verify the advantages of the RBFNNO method, fuzzy PID control and BP neural network control were selected for a comparison analysis. The PAM trajectory tracking is shown in Figure 17. It can be seen that the error convergence is fast, and the control accuracy is high when using the RBFNNO method.

4.4. Hand Exoskeleton Control

To validate the role of the hand exoskeleton in aiding rehabilitation training for patients, passive rehabilitation training for each finger was conducted after the patient was fitted with the exoskeleton. Each finger was controlled by individual PAMs with minimal coupling between fingers. The desired trajectory was tracked using the RBFNNO methods.

Based on the rehabilitation needs, the maximum bending angle was designed to be 1.2 radians, with one inflation and deflation cycle lasting 25 s for slow rehabilitation. When the angle signal is tracked stably, the angle tracking trajectory and angle trajectory tracking error results for the five fingers are shown in Figure 18.

The maximum error of the five fingers is primarily due to the hysteresis effect and jitter. As shown in

Table 3. Trajectory tracking error for five fingers.

Finger	Maximum Error/Rad	Mean Error/Rad	Mean Square Error/Rad
Thumb	0.3308	0.0651	0.1219
Index finger	0.3596	0.0788	0.1141
Middle finger	0.3291	0.0815	0.1183
Ring finger	0.3194	0.0803	0.1179
Little finger	0.3517	0.0705	0.1217

, the mean squared deviation of the thumb is greater than that of the index finger despite having relatively small mean error values. From the Figure 18, it is evident that the hysteresis effect of the thumb is slightly smaller than that of the index finger. For the shorter PAMs, it is generally easier to return to the initial position during the deflation phase, although their jitter and vibration are more pronounced.

In the active grasping mode, the trajectories of the five healthy fingers during object grasping were first collected to serve as the target trajectory for the hand exoskeleton. The subjects were instructed to sit in a relaxed state at the center of the motion capture space and perform full grasping movements using their healthy hand. The hand movement coordinates were recorded from the start of grasping, through object manipulation, to the completion of the grasping action.

Figure 19a showed the angle changes of the thumb and little finger during gripping actions in a healthy hand. Figure 19b displayed the angle changes of the remaining fingers during the same action. The extracted finger trajectories from the healthy hand served as the desired trajectories for the exoskeleton.

The exoskeleton assisted the affected hand in completing grasping movements according to the trajectories recorded from the healthy hand. Figure 20 showed the gripping action of a healthy hand and the affected hand wearing the exoskeleton while grasping the same object. Figure 21 illustrated the trajectory changes of the affected hand’s index finger during the gripping process along with the error compared to the healthy hand’s index finger trajectory. As demonstrated in Figure 20 and Figure 21, the hand rehabilitation exoskeleton effectively aids patients with muscle weakness in grasping different objects.

5. Discussion

By analyzing the mechanism of finger joint movement and the rehabilitation needs of the hand, a parallel chamber PAM with bidirectional bending capabilities was designed to assist patients in finger flexion and extension movements. The bending performance of the PAM with different structural parameters was compared using a finite element analysis, and the structural parameters were optimized based on actual rehabilitation needs. The mold design was carried out in SOLIDWORKS, and the mold was fabricated using 3D printing technology.

Compared to rigid hand rehabilitation exoskeletons [26] and motor-driven soft hand rehabilitation exoskeletons [27], the designed actuator is lighter and soft, with a single PAM weighing only 25 g, making it more comfortable to wear. Compared to existing soft hand exoskeletons driven by PAMs, our developed exoskeleton features an improved structure and algorithm, resulting in enhanced performance and higher precision. The bending angle of the PAM designed by Li et al. [28] for hand rehabilitation is only around 0 to 70°, which is much smaller than the −50° to 125° range of our developed PAM. The fingertip force of the rehabilitation glove developed by Han et al. [12] is approximately 0.40 N to 0.50 N, which is lower than the fingertip force of 1.25 N of the exoskeleton we developed. Additionally, many PAM-driven hand rehabilitation exoskeletons [15] only assist with finger bending and cannot help patients with muscle weakness extend their fingers for repeated grasping. The bidirectional bending capability of our PAM better supports finger extension and enables patients to complete multiple grasping tasks in daily life.

We derived a dynamic model of the PAM by substituting the reduced-order forms of strain energy, generalized force, kinetic energy, and gravitational potential energy into the Lagrange equation. Given the requirement for precise position control of PAMs during rehabilitation training, an optimized radial basis function neural network adaptive control algorithm was developed, focusing on parameter and hidden layer node optimization. The algorithm utilizes a k-means clustering algorithm to determine the centers of the radial basis functions and optimizes their widths based on sample distribution, thereby accelerating network convergence and enhancing the position control accuracy. This approach offers higher accuracy and greater stability compared to control methods such as PID, fuzzy logic, and unoptimized neural networks used in other hand rehabilitation exoskeletons. For example, the PID control method used by Tejada et al. [29] can only bend the pneumatic muscles to 90% of the maximum expected angle. Haghshenas et al. [30] reported an average error of 12.4° in the movement of the index finger during bilateral therapy.

We designed both a passive training mode and an active grasping mode for the hand exoskeleton. In the early stages of rehabilitation, patients undergo passive flexion and extension training for individual fingers or the entire hand. In the later stages, grasping rehabilitation training is implemented using joint trajectory mapping between the healthy and affected hands to enhance hand muscle strength during active rehabilitation. The combination of these modes provides more effective support for patients in completing rehabilitation tasks and daily activities.

Our research currently has some shortcomings. Single fingers can only flex and extend, lacking the ability for independent control of multiple joints in one finger or for abduction/adduction. The thumb’s multi-degree-of-freedom motion requires more advanced PAMs. We plan to improve the design of PAMs and the exoskeleton to better match hand movement characteristics. Additionally, using motion capture systems or multiple bending sensors may make the exoskeleton cumbersome and difficult to use, so we aim to investigate PAMs with embedded sensors to enhance wearability. Since hand sizes vary among users, the required sizes for the exoskeleton and PAMs will differ as well. To improve usability and simplify manufacturing, we need to develop an automated design tool that customizes PAM and exoskeleton parameters based on individual hand sizes.

In terms of control, we also need to reduce response time and enhance human–machine interaction. The current control method has significant errors and delays during the inflation and deflation phases, which require improvement. In the future, we could maintain a certain pressure in one chamber while inflating the other. This control method, resembling antagonistic muscle pairs, would allow us to adjust the stiffness of the PAM. Additionally, we hope to incorporate physiological signals, such as electromyographic (EMG) signals, to monitor the patient’s movement capabilities in real time. This would enable the exoskeleton to automatically switch between active and passive grasping training.

6. Conclusions

This paper explores the design, modeling, and control methods of bidirectional bendable PAMs for the soft hand rehabilitation exoskeleton. A PAM dynamic model based on the Lagrangian method is established, and an adaptive control algorithm using an RBFNN with parameter and node optimization is designed to achieve accurate position control. The integration of PAMs with fabric gloves results in a softer and more wearable hand rehabilitation exoskeleton compared to motor-driven hand exoskeletons, better aligning with the rehabilitation needs of stroke patients. We have developed both passive training and active grasping modes to enhance patients’ grasping ability and increase their motivation for rehabilitation. Next, we plan to enhance the functionality of exoskeletons by improving the performance of PAM and integrating physiological signals such as EMG. We plan to collaborate with hospitals in the future to conduct more experiments involving patient participation. This involvement will better validate the effectiveness of the rehabilitation exoskeleton across various rehabilitation modes.