A Method for Precise Tracking Control of Pneumatic Artificial-Muscle-Driven Exoskeletal Robot

Abstract

Exoskeletal robots are of critical importance in the domain of mechanical boosting. The pneumatic artificial muscle (PAM) is commonly used as a flexible actuator in exoskeletal robots designed for upper limbs due to its high power-to-weight ratio, conformability, and safety. This study establishes a new model based on the existing model to improve its control precision by implementing elastic and frictional forces and empirical coefficients, battling against the time-variant hysteresis that PAM’s output force exhibits. In the meantime, a BP neural network is employed in reverse modeling, followed by the adoption of the least-square-based particle swarm optimization algorithm in order to determine the optimized parameter values. PAM provides the Upper Limb Exoskeletal Robot with appropriate auxiliary power, which can be adjusted to accommodate variations in posture change during the lifting process. PAM is also capable of handling variable loads based on the principle of torque balance, constructing a control system according to the inverse dynamics of exoskeletal robots accompanied by an inverse model of PAM’s output force, and finally, rendering tracking control of the elbow angle during the auxiliary process possible. Finally, the tracking error results are calculated and shown; the maximum angular error in the tracking process is 0.0175 rad, MAE value is 0.0038 rad, RMSE value is 0.0048 rad, and IEAT value is 4.6426 rad. This control method is able to improve the precision of tracking control of the elbow angle of the upper limb–exoskeleton coupled system during the process of lifting goods.

1. Introduction

Exoskeletal robots for upper limbs combine human physical strength with mechanical power, integrating techniques such as bionics, mechanics, transducing, control, and algorithms in themselves [1,2,3], thus becoming intelligent devices that can enhance work efficiency and life quality. In recent years, they have also been intensively researched, both domestically and abroad [4,5,6]. Upper extremity exoskeletons offer several advantages in the field of mechanical boosting, which generally refers to enhancing or augmenting human physical capabilities. By providing mechanical assistance to the upper body, exoskeletons help reduce muscle fatigue. This is valuable in professions that involve repetitive or prolonged use of the upper extremities, such as construction, agriculture, and manufacturing, as it can help prevent overexertion and musculoskeletal injuries. In Europe, three-fifths of workers suffer from work-related diseases, of which 60% are musculoskeletal diseases. In total, 40% of these diseases are upper limb-related. During industrial lifting, workers are required to hold weights above their heads for an extended period. The torque generated by gravity and impounding of the elbows and shoulders of workers weighs 50% of their body weight. Therefore, exoskeletal robots have a significant effect in the domain of industrial lifting and rehabilitation from sports injuries [7,8]. Most exoskeletal devices use electrical, hydraulic, or pneumatic actuators. Safety, conformability, controllability, and power–weight ratio should be taken into consideration when selecting actuators to assist with lifting. Many researchers favor electrical motors for their controllability. For example, Iqbal et al. proposed a novel design of a hand exoskeleton system [9]. The optimization of the exoskeleton device was achieved through a procedure targeting the natural finger workspace and capabilities. Moubarak et al. proposed an upper extremity exoskeleton mounted on a wheelchair [10]. This new device is dedicated to regular and efficient rehabilitation training for weak and injured people without the continuous presence of a therapist. However, its weight is an obvious downside. Hydraulic actuators have an outstanding power–weight ratio, while possible leaking of hydraulic fluid can cause a lot of problems in medical aseptic environments [11]. In contrast thereto, pneumatic actuators possess a favorable power-weight ratio but cause no medical pollution. A pneumatic muscular actuator (PMA) is light and conformable itself. Such conformability softens the counterforce against resistance. In the meantime, a pneumatic muscle is similar to skeletal muscle in terms of biomechanical biconical results [12]. This renders the PMA the best choice for assisting with goods lifting.

The pneumatic muscle assembly consists of three main elements: the outer braided sleeve, the flexible inner bladder, and end-fittings at each of its two ends. This study takes a self-assembled McKibben-style pneumatic muscle as an example, as shown in Figure 1. A PMA operates primarily by inputting pressured air into the bladder, causing it to expand, and applying a radial outward force to the surrounding braided sleeve, which transforms the aforementioned force into axial adduction (or displacement) of the muscle. However, hysteresis, as an inherent characteristic of soft actuators, can seriously affect the accuracy of closed-loop control and practical application [13]. This non-linear characteristic can be observed in the output force through a non-singular and non-convex hysteresis curve [14]. Existing models cannot accurately describe this phenomenon [15], and so far, there is no versatile method to model this hysteresis of soft actuators [16]. However, soft actuators have been more commonly employed on upper limbs in recent years. In 2018, Yu proposed a new 7-DOF upper limb robot exoskeleton [17]. Each arm includes force and torque transducers to detect human motion. In 2020, Chen proposed an Upper Limb Exoskeleton Rehabilitation Robot based on PMA, which can plan the rehabilitation trajectory through curve fitting, and use fuzzy synovial controllers (FSMC) for rehabilitation control [18], but the way to enhance the accuracy of the model is yet to be determined. In 2021, Proietti et al. developed a multi-joint soft exosuit for upper limb assistance and rehabilitation using dynamic GC and joint trajectory tracking controllers [19]. IMU sensors were used to control the pressure for the pneumatic actuators. In 2022, Chiou developed a 4-DOF wearable upper limb rehabilitation assistance exoskeleton system (WURAES) that is suitable for assisting in the rehabilitation of patients with upper limb injuries [20]. a proxy-based output feedback sliding mode control (POFSC) was developed to provide appropriate rehabilitation assistance power for the upper limb exoskeleton and to maintain smooth and safe contact between the WURAES and the patient.

Therefore, the pneumatic muscle output force can affect its control over the motion of the exoskeletal robot. Currently, Upper Limb Exoskeletal Robots’ coupling with humans has reduced their auxiliary performance due to actuator control. Furthermore, they fail to take load variation into consideration. This study takes the perspective of the hysteresis of pneumatic muscle output force, trying to introduce elastic and friction forces into the existing PAM model to enhance it. In the meantime, a BP neural network will be employed in reverse modeling. Furthermore, this study will endeavor to incorporate the principle of torque balancing to handle variable loads, construct a control system according to the inverse dynamics of the exoskeletal robot accompanied by the inverse model of PAM’s output force, and finally, render tracking control of the elbow angle during the auxiliary process possible.

2. Pneumatic-Muscle-Driven Upper Limb Exoskeletal Robot System

In this section, the control system for the pneumatic-muscle-driven Upper Limb Exoskeletal Robot will be introduced. The precision of the output force model of the pneumatic muscle is related to the tracking control precision of the exoskeletal robot. Hence, this section will introduce elastic and frictional forces into the common PAM output force model to establish a new model. Employing collected data from experiments, a BP neural network is built to create an inverse model of output force. This will output the atmospheric pressure value calculated using the output force and its length, enhancing the precision and model itself. Furthermore, this section will elaborate on the control system of this pneumatic-muscle-driven soft Upper Limb Exoskeletal Robot.

2.1. PAM Driver

PAM is a soft actuator. It can overcome the barrier of the conformability of hard actuators, thus enhancing adaptability. However, due to the existence of frictional forces while extending and retracting, PAM exhibits hysteresis, which affects the calculational precision of PAM’s output force.

2.1.1. PAM Output Force Model

The following model was introduced by C.P. Chou [21]. This mathematical model, as shown in equation (1) below, is currently a simplified mechanism of PAM’s output force. Nonetheless, this model is based on the following hypotheses: (a) strands cannot be extended; (b) there are frictional forces between the chamber and the strands enveloping it; (c) the outreach and adduction of muscle are always in an ideal state; no relative motion will take place between the two aforementioned parts; (d) PAM can be approximated as a cylinder by disregarding the structural changes at the muscle end seal. The geometrical relationship of PAM is shown in Figure 2.

$$F(\epsilon ,p)=\pi r_0^{2}p\left[a{(1-\epsilon )}^2-b\right]$$

(1)

In the formula, $a=3/{\tan }^2{a}_0$; $b=1/{\sin }^2{a}_0$; $a_0$ is PAM’s initial braid angle; $r_0$ is the initial radius; $\epsilon =\left(l_0-l\right)/l_0$; $l_0$ is the initial length; $l$ is PAM’s length; and $p$ is PAM’s inner pressure.

In order to accurately model the output force of PAM, the given output force model is obviously insufficient. Elastic forces between weaved strands and the rubber bladder, frictional forces between strands and the rubber bladder, and frictional forces between strands are of necessity when constructing the mathematical model of PAM.

The inner strand of PAM is made from nylon, giving it higher structural strength. This enables the muscle to transform the power provided by air into output force. That is to say, the axial adduction force is actually provided by the strand mesh. Great air pressure is put on the strand during PAM’s axial outreach, causing extortion, and thus, storing energy. This portion of energy can be calculated by measuring the extortion of the strand. As shown in Figure 3, the tractional force on the mesh can be deduced by balancing the axial force. The elastic force generated by the strand can be inferred by calculating the elastic energy’s derivative in the displacement direction, as shown in Formula (2).

$$F_{\mathrm{braid} }=\frac{d{W}_{\mathrm{braid} }}{dL}=N{b}_0{A}_{s_0}\frac{d{V}_{\epsilon }}{dL}=\frac{N^2{b}_0^3{A}_{s_0}^2{E}^2{p}^2\pi nL}{{\left(\pi nNE{A}_{s_0}-pL{b}_0\right)}^3}$$

(2)

In the above formula, $A_{s_0}$ is the nominal cross-sectional area of the strand, $E$ is the strand’s elastic modulus, ${\mu }_i$ represents the property parameter of rubber, N is the number of strands, $W_{\mathrm{braid} }$ is the stored elastic potential energy within, $V_{\epsilon }$ is the strain energy density, and b is the length of a single strand.

The inner expandable chamber of PAM is commonly made of rubber. This is due to the elasticity of the rubber bladder’s inner cavity and the prevention of accurate modeling during the inflation process by the outer enveloping strand of the bladder. The elastic energy of the bladder can be deduced according to the virtual work principle. This energy can be represented indirectly by the elastic energy density W:

$$F_{\mathrm{bladder} }=V_r\frac{dW}{dL}$$

(3)

In this equation, $V_r$ is the bladder volume. Taking $t_0$ and t as the bladder’s inherent thickness and its thickness while inflating, $L_0$ and L are PAM’s length under no external influence and during inflation, respectively. The volume of the bladder is:

$$\begin{array}{c}{V}_r=\pi L_0\left[{\left(\frac{D_0}{2}\right)}^2-{\left(\frac{D_0}{2}-t_0\right)}^2\right]\\ =\pi L_0{t}_0\left(D_0-t_0\right)\end{array}$$

(4)

For rubber material, the elastic energy density W can be represented by the Mooney–Rivlin model:

$$W=C_{10}\left(I_1-3\right)+C_{01}\left(I_2-3\right)$$

(5)

where $I_1={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$, $I_2={\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2+{\lambda }_1^2{\lambda }_3^2$, C_ij is the material constant, and I₁ and I₂ are the Green strain invariant. ${\lambda }_i$ is the main elongation. Because the deformation process of rubber material is considered to be isotropic and incompressible, ${\lambda }_1{\lambda }_2{\lambda }_3=1$, and

$$I_2=\frac{1}{{\lambda }_3^2}+\frac{1}{{\lambda }_1^2}+\frac{1}{{\lambda }_2^2}$$

(6)

Usually, ${\lambda }_1=\frac{L}{L_0},{\lambda }_2=\frac{D}{D_0},{\lambda }_3=\frac{1}{{\lambda }_1{\lambda }_2}$. This paper takes into consideration the wall thickness of muscles, and chooses the following representation method for the primary elongation rate:

$${\lambda }_1=\frac{L}{L_0},{\lambda }_2=\frac{D-t}{D_0-t_0},{\lambda }_3=\frac{t}{t_0}$$

(7)

where t₀ represents the initial thickness of the rubber, while t represents the current thickness of the rubber. Therefore, the stored elastic energy within the rubber bladder can be deduced as follows:

$$\begin{array}{c}{F}_{\mathrm{bladder}}=V_r\left\{2C_{10}\left[{\lambda }_1\frac{d{\lambda }_1}{dL}+{\lambda }_2\frac{d{\lambda }_2}{dL}+{\lambda }_2\frac{d{\lambda }_2}{dL}\right]\right.\\ \left.+2C_{01}\left[{\lambda }_1\left({\lambda }_2^2+{\lambda }_3^2\right)\frac{d{\lambda }_1}{dL}+{\lambda }_2\left({\lambda }_1^2+{\lambda }_3^2\right)\frac{d{\lambda }_2}{dL}+{\lambda }_3\left({\lambda }_1^2+{\lambda }_2^2\right)\frac{d{\lambda }_3}{dL}\right]\right\}\end{array}$$

(8)

In the self-assembled PAM, the bladder is mixed with strands, leading to full contact between its cylindrical surface and the strands. The stress between the bladder and outer interwoven strands equals the air pressure in the cylinder. Thus, the frictional force between the bladder and strand should be:

$$F_{s-r \mathrm{friction} }=f^{\prime }\cdot P\cdot \pi DL$$

(9)

During the adduction of the pneumatic muscle, relative displacement will occur between strands. The frictional force caused by it will directly hinder the accuracy of such a model. Davis and Caldwell, based on the Hertzian Contact Theory, proposed the formula for frictional force between PAM strands [22].

$$F_{\mathrm{strand}\text{-}\mathrm{strand} }=f_{\mathrm{strand}\text{-}\mathrm{strand} }\frac{S_{\mathrm{contact} }}{S_{\mathrm{scale} }}p$$

(10)

In this formula, $f_{\mathrm{strand}\text{-}\mathrm{strand} }$ is the friction coefficient between interwoven strands, $S_{\mathrm{contact} }$ is the contact area between strands, $S_{\mathrm{scale} }$ is adjustment coefficient for the contact area, $F_{\mathrm{comp} }$ is the extrusion force between strands, ${\theta }_{\min }$ is the minimum braiding angle, and v is the strands’ Poisson rate.

In the existing relevant model, mechanism modeling based on PAM does not consider the impact of elastic force and friction in the process of PAM movement, so it is more considered as empirical modeling and is not based on mechanism modeling, resulting in complex expressions and too many parameters. This study, based on a common model of output force, implements the elasticity of the strands and rubber bladder, the friction between the two parts, and the friction between strands into the PAM mathematical model, which will enable the model to better simulate muscle performance. In the meantime, the empirical coefficient is introduced as a constant, and the coefficient in question will also be determined using the Particle Swarm Algorithm. This will not only lower the difficulties in obtaining empirical coefficients, but also raise their precision to a higher level. The final optimized model is shown below:

$$\left\{\begin{array}{c}F=k_1\pi r_0^{2}p\left[a{(1-k_2\epsilon )}^2-b\right]+k_3-F_{\mathrm{bladder} }-F_{\mathrm{braid} }-\left(F_{s\text{-}\mathrm{sfriction} }+F_{s\text{-}\mathrm{rfriction} }\right)inflation\\ F=k_1\pi r_0^{2}p\left[a{(1-k_2\epsilon )}^2-b\right]+k_3-F_{\mathrm{bladder} }-F_{\mathrm{braid} }+\left(F_{s\text{-}\mathrm{sfriction} }+F_{s\text{-}r \mathrm{friction} }\right)deflation\end{array}\right.$$

(11)

The purpose of system identification is to determine the best parameter scheme according to the provided input and output data. Aiming at producing an output force model, this study applies the least-square method as the function of adaptability, stops when a suitable value is reached, and then, obtains the best parameter using the Particle Swarm Algorithm. According to Hill’s tension–velocity curve, the maximum peak power is about 0.1F_max × v_max. The termination criterion of PSO in this paper is that the fitness function is less than or equal to 0.1 N, or the maximum peak power of PAM is more than 0.1F_max × v_max.

$$f(F_{\mathrm{i}})=\frac{1}{H}{\displaystyle \sum _{i=1}^H{\left(F_i-{\hat{F}}_i\right)}^2}$$

(12)

Here, $f(F_{\mathrm{i}})$ is the objective function of the particle swarm optimization algorithm, F_i is the actual output force of the PAM system, ${\hat{F}}_i$ is the output value in the FN-QUPI model, and H is the number of relevant samples during the identification process.

2.1.2. Inverse Model of PAM Output Force

According to the output force formula, output force is determined based on load and atmospheric pressure, while the inverse model is responsible for calculating atmospheric pressure from output force, load, and length. Due to a strong non-linear relationship, great challenges have been posed for determining the inverse model. Therefore, this study uses a BP neural network as a method to build an inverse model of PAM’s output force.

A BP neural network consists of an input layer, a hidden layer, and an output layer. Layers are all connected to each other; however, for units in the same layer, there is no connection. Information is transmitted upward, from the input layer to the output layer. The BP neural network algorithm adjusts and modifies the connection weighting factor through the inverse transmission of the network output error, thus minimizing its error and establishing mapping from the input to output among a given number of samples. The three-layer BP neural network has a topological structure, as described below: For the input layer, there are n nodes; there are q nodes on the hidden layer; and there are m nodes on the output layer. The weighting factor for connections between input layer and hidden layer is v_ki; for connections between the hidden layer and output layer, this factor is w_jk. The BP neural network is trained as follows:

• Initialize the weighting factor: randomly assign factors in the network v_ki (0) and w_jk (0)
• Output of the calculated sample data is given, employing the initialized values determined in step 1. The neural network will then give outputs for every neuron. The output of the hidden Layer is:

$$h_k=f_h\left({\displaystyle \sum _{i=1}^n{v}_{ki}}{x}_i-b_k\right)$$

(13)

Here, k = 1, 2, …, q:x are sample input data, b_k is neuron bias, and f_h is the transfer function of the hidden layer neuron.

The nodes on the output layer give outputs, as shown below:

$$y_j=f_o\left({\displaystyle \sum _{k=1}^q{w}_{jk}}{h}_k-b_j\right)$$

(14)

Here, j = 1, 2, …, m, and f_o is the transfer function for output layer.

3.

Then, errors are calculated reversely: According to the expected output and the actual output of all neurons from the given sample data, the errors of all neurons can be calculated in reverse with the result. The input error is:

$$E=\frac{1}{2}{\displaystyle \sum _{j=1}^m{\left(g_j-y_j\right)}^2}$$

(15)

4.

The correction coefficient: According to the principle of the gradient drop of the error, the weighting factor, previously assigned in step 1, will be corrected layer by layer; the correction coefficient for the weighting factor is:

$$w_{jk}^{N+1}=w_{jk}^N+\eta {\delta }_j^o{h}_k$$

(16)

$${\delta }_j^o=\left(g_j-y_j\right)y_j\left(1-y_j\right)$$

(17)

The correction coefficient for the input layer’s weighting factor is:

$$v_{ki}^{N+1}=v_{ki}^N+\eta {\delta }_k^h{x}_i$$

(18)

$${\delta }_k^h=\left[{\displaystyle \sum _{k=1}^q\left(g_j-y_j\right)}{y}_j\left(1-y_j\right)w_{jk}\right]h_k\left(1-h_k\right)$$

(19)

5.

Loop condition and Termination: If the errors do not fulfill the expectation, then the calculation will enter a loop by jumping to step 2. If expectations are met, then the iteration will be terminated, and the final weighting factor and errors during training are recorded.

In order to collect PAM’s variation in displacement, air pressure, load, and output force, The experimental setup as shown in Figure 4 is assembled. A pneumatic cylinder (RIGHT PNEUMATIC, Ningbo, China) is used to represent the load, thus reducing random errors caused by weight suspension. In this experiment, PAM is horizontally placed. One of its ends is anchored to the experiment surface through a force transducer (MEIKONG, Hangzhou, China); the other end is connected to a baffle equipped with a laser displacement transducer (Panasonic INDUSTRY, Tangshan, China); this is further connected to a pneumatic cylinder. The inner and outer air pressures of PAM are regulated by a proportional valve (NUONENGTAI, Guangdong, China). As shown by Figure 5, DSPACE (dSPACE, Paderborn, Germany) on a computer can transmit an analog signal to the proportional valve, which will, with the aid of an air compressor built into the air source, convert the voltage signal into a pressure signal and either inflate or deflate the pneumatic cylinder and PAM. The laser displacement sensor and tension–compression sensor will output the displacement and tension in the form of voltage signals and send them back to the computer. The recorded voltage signals sent by the transducers are then translated into displacement and tension through Simulink. Here, the input voltage of PAM and the pneumatic cylinder, the output force of the tension–compression transducer, and the displacement of the laser displacement transducer are given. After adjusting the parameters in MATLAB 2022b, the discrete data regarding air pressure for adduction and outreach, displacement, and output tension are implemented into the neural network and trained.

2.2. Flexible Upper Limb Exoskeletal Robot

The structure of the Flexible Upper Limb Exoskeleton Robot designed in this study is indicated in Figure 6. To enhance the strength of human upper limbs and to make the device easy to wear so that the wearer can carry goods without hinderance, the robot has the following features in terms of structural design: (1) Adjustable size to adapt to a larger group of people. The rotational centers of exoskeleton joints are in conformity with those of human joints, so there is no conflict between the exoskeleton joints and real human joints. (2) Flexible material is chosen to lower the overall weight. With the pneumatic muscle chosen as the driving method, flexible textiles are used to bond the person and the exoskeleton, so there will not be any limitations like those that occur with rigid materials. (3) The overall assistive function is perfected in strength. To reduce damage to the shoulder and back joints of the human body during the goods-carrying process, supporting structures for the shoulder and back are added. At the shoulder joint, bearing hinges are used to improve the angle limitations of the robot during the adduction and outreach processes. The redundant structure of the holder at the back is reduced to lower the overall weight of the exoskeleton. In this case, with the minimum strength that human muscles can bear, more weight can be loaded. To enhance the adaptability of the exoskeleton robot, the wrist joint and the forearm structure are simplified.

In terms of joint design, due to the complicated physiological structure of the human body, human joints do not simply rotate on a plane surface. In this study, the joint rotation limit is simplified on the premise that the key joint motion during goods carrying is not affected. For arms, the angle limitation of shoulder joints is simplified to three levels, and the elbow and wrist joints are expressed by buckling and stretching freedom on the sagittal plane to satisfy movement needs. During the movement of shoulder joints, the Bowden cables are passively driven to transmit the force at a long distance within a small space, which also enhances the flexibility of the joints. When operators are carrying goods while wearing exoskeleton robots, since the shoulder joint is experiencing a certain weight, the upper arm of exoskeleton robots may lean backward. To guarantee the stability of the shoulder joints, one end of the Bowden cables is fixed on the back holder of the exoskeleton, and the other end is fixed on the guide pulleys on the shoulder. When carrying heavy goods, the elbow is bent at 90 deg, at which time the Bowden cables are tight and can provide some assistive force to stabilize and relieve the shoulder from fatigue.

To enhance the auxiliary force at the back without adding more weight to the exoskeleton, new textiles and flexible materials are employed. Thus, the elastic bands will store more elastic potential energy after expanding. When bending over, the flexible material will expand, so when the operator resumes an upright posture, the device will provide auxiliary force to the upper body.

To reduce the weight of the robot, its battery is separated from the main body, and the driving components, magnetic valve, and other hardware are integrated into an independent control cabinet, which can effectively reduce the weight of the robot so that the human body will not bear so much weight. The condition of the joint inversion was designed to averted, and the operator’s physical health would be protected.

3. The Tracking Control System of the Flexible Upper Limb Exoskeletal Robot during the Lifting Process

This section is dedicated to the introduction of the tracking control system of the Flexible Upper Limb Exoskeletal Robot. Quintic polynomials were employed while designing the target angle, angular speed, and angular acceleration for the elbow joint. While lifting goods with the assistance of the Flexible Upper Limb Exoskeleton, the load will adapt itself to the changing posture. The exoskeletal robot can deal with variable loads by applying the principle of torque equivalence, performing an analysis on the system from a kinematic and dynamic perspective while taking the differences between the expected and actual angle, angular speed, and acceleration into consideration, and finally, calculating the torque through robotic inverse dynamics. The difference between the obtained torque and the passive torque provided by the load is called the control torque. This torque is applied to PAM. Then, the corresponding air pressure can be calculated according to the principle of moment equivalence and the inverse model of PAM’s output force. Thus, it helps in designing and developing the controller with less effort and at cheaper cost [23].

3.1. Handling of Variable Load by Upper Limb Exoskeleton

First of all, a body-fixed coordination system is built based on the left and right elbow joints, designating the center point of both elbow joints as its point of origin. While lifting with the Upper Limb Exoskeleton, the angle of lifting varies, and the interactive forces between cargo and robot and the vertical component force of gravity will also change accordingly. Thus, the lifting process is a problem of variable load.

Taking the left arm, the rotation angle of the elbow joint during lifting is q, as shown in Figure 7. a represents the distance from the attachment point of the pneumatic muscle to the elbow joint on the upper arm. b represents the distance from the attachment point of the pneumatic muscle to the elbow joint on the forearm. r represents the vertical distance from the elbow to the pneumatic muscle. The main parameters of the upper extremity exoskeleton’s system designed in this paper are shown in

Table 1. Main parameters of the upper extremity exoskeleton.

Parameter	a (m)	b (m)	m (kg)	ZZ (kg·m2)
Value	0.165	0.21	5	0.05

: a is the distance between the pneumatic muscle and the upper arm connection point and the elbow joint; b is the distance between the pneumatic muscle and the forearm connection point and the elbow joint; m is the mass of the load being carried; the upper extremity exoskeleton’s handling process is represented by rotation around the Z axis, and its moment of inertia regarding the Z axis rotation is ZZ.

According to the second law of Newton, the following analysis can be conducted regarding the lifted object:

$$T-mg\cos \left(q\right)=mr\ddot{q}$$

(20)

$$T=mr\ddot{q}+mg\cos \left(q\right)$$

(21)

$$f=mg\sin \left(q\right)$$

(22)

Regarding the Upper Limb Exoskeletal Robot, such an analysis of interaction forces can be conducted, from which conclusions can be drawn that the total pressure of the goods carried can be determined by multiplying the pressure in the center-of-mass coordinate system and the direction vector in the world coordination system. Therefore, the total pressure put on the robotic arm can be calculated as follows:

$$T^{\prime }=T\left(\begin{array}{c}-\sin \left(q\right)\\ \cos \left(q\right)\\ 0\end{array}\right)$$

(23)

The frictional force between the goods and the robotic arm is:

$$f^{\prime }=mg\sin \left(q\right)\left(\begin{array}{c}-\cos \left(q\right)\\ -\sin \left(q\right)\\ 0\end{array}\right)$$

(24)

Therefore, under the world coordination system, the joint force put upon the robotic arm is:

$$\begin{array}{cc}F& =T\left(\begin{array}{c}\sin \left(q\right)\\ -\cos \left(q\right)\\ 0\end{array}\right)+mg\sin \left(q\right)\left(\begin{array}{c}-\cos \left(q\right)\\ -\sin \left(q\right)\\ 0\end{array}\right)\\ & =\left(mr\ddot{q}+mg\cos \left(q\right)\right)\left(\begin{array}{c}\sin \left(q\right)\\ -\cos \left(q\right)\\ 0\end{array}\right)+mg\sin \left(q\right)\left(\begin{array}{c}-\cos \left(q\right)\\ -\sin \left(q\right)\\ 0\end{array}\right)\\ & =\left(\begin{array}{c}-mr\ddot{q}\sin \left(q\right)\\ mr\ddot{q}\cos \left(q\right)-mg\\ 0\end{array}\right)\end{array}$$

(25)

The passive torque generated by the load on the Upper Limb Exoskeleton robot can be expressed as Equation (24) in this study. By utilizing the principle of symmetry, it can be deduced that this torque is equally divided between the left and right mechanical arms. Thus, the torque experienced specifically by the left arm can be calculated using Equation (25).

$$p_{\tau }=\left(\begin{array}{c}r\cos \left(q\right)\\ r\sin \left(q\right)\\ 0\end{array}\right)\times F$$

(26)

$${}_l{}^p\tau =0.5\times p_{\tau }$$

(27)

The torque generated by the load causes a downward movement tendency in the Upper Limb Exoskeleton’s mechanical arm. To maintain balance, the pneumatic muscles need to provide equal and opposite torque. The control variable in the system is the air pressure. In order to accomplish the task of assisting in lifting, it is necessary to calculate the corresponding air pressure for the pneumatic muscles at this time. Taking the pneumatic muscle on the left mechanical arm as an example, the length of the pneumatic muscle at this time can be expressed using Equation (27). By using the inverse model of the muscle output force of the pneumatic muscle, the corresponding air pressure for the pneumatic muscle at this time can be solved, as shown in Equation (29).

$$r=\frac{ab\sin \left(\frac{\pi }{2}-q\right)}{{}^{l}l}$$

(28)

$${}^{l}l=\sqrt{a^2+b^2-2ab\cos \left(0.5\times \pi -{}^{l}q\right)}$$

(29)

$${}_l{}^{p}f=\frac{{}_l{}^p\tau }{{}^{l}r}$$

(30)

$${}_l{}^p\mathrm{p}={}^{l}F{}^{-1}\left({}^{l}l,{}_l{}^{p}f\right)$$

(31)

3.2. Kinematic and Dynamic Analysis of Upper Limb Exoskeleton during Lifting Process

In this paper, the Lagrangian formulation is adopted for dynamic modeling, providing a description of the coupling between the torque and structural motion of the elbow joint and pneumatic muscles in the Upper Limb Exoskeleton. The main principle involves establishing a center-of-mass coordinate system C at the position of the load’s center of mass and a body-fixed coordinate system J at the axis of the left and right elbow joints. This enables the coordinate transformation of the center-of-mass coordinate system relative to the world coordinate system W, as shown in Figure 8. After the elbow joint rotates by an angle $\theta$, the new coordinate axis is as shown by the red coordinate axis. The left and right arms are represented by 1 and 2, respectively. The position of the origin of J1 relative to W is denoted as ${\left(\begin{array}{ccc}0& 0& -L\end{array}\right)}^T$, where $q_1,{\dot{q}}_1,{\ddot{q}}_1$ represent the rotational angle, angular velocity, and angular acceleration of J1 relative to W, respectively. The position of the origin of C1 relative to J1 is denoted as ${\left(\begin{array}{ccc}{}^{1}C_x& {}^{1}C_y& {}^{1}C_z\end{array}\right)}^T$.

Taking the left arm of the exoskeleton robot with the upper limb as an example, when the joint rotation angle is q1, the orientation matrix ${}_{J1}{}^{W}T$ of the conjoined body coordinate system J1 is relative to the world coordinate system W, the orientation matrix ${}_{C1}{}^{J1}T$ of the centroid coordinate system C1 is relative to the conjoined body coordinate system J1, and the orientation matrix ${}_{C1}{}^{W}T$ and rotation matrix ${}_{C1}{}^{W}R$ of the centroid coordinate system C1 are relative to the world coordinate system W.

$${}_{J1}{}^{W}T=\left(\begin{array}{cccc}\cos \left(q_1\right)& -\sin \left(q_1\right)& 0& 0\\ \sin \left(q_1\right)& \cos \left(q_1\right)& 0& 0\\ 0& 0& 1& -L\\ 0& 0& 0& 1\end{array}\right)$$

(32)

$${}_{C1}{}^{J1}T=\left(\begin{array}{cccc}1& 0& 0& {}^{1}C_x\\ 0& 1& 0& {}^{1}C_y\\ 0& 0& 1& {}^{1}C_z\\ 0& 0& 0& 1\end{array}\right)$$

(33)

$${}_{C1}{}^{W}T=\left({}_{J1}{}^{W}T\right)\left({}_{C1}{}^{J1}T\right)$$

(34)

$${}_{C1}{}^{W}R={}_{C1}{}^{W}T\left(1:3,1:3\right)$$

(35)

Based on the pose matrix ${}_{C1}{}^{W}T$, the Roll–Pitch–Yaw (RPY) Euler angles ${\alpha }_1,{\beta }_1,{\gamma }_1$ of the centroid coordinate system C1 with respect to W can be obtained.

$${\alpha }_1=\arctan \left(\frac{{}_{C1}{}^{W}T\left(3,2\right)}{{}_{C1}{}^{W}T\left(3,3\right)}\right),{\beta }_1=\arcsin \left(-{}_{C1}{}^{W}T\left(3,1\right)\right),{\gamma }_1=\arctan \left(\frac{{}_{C1}{}^{W}T\left(2,1\right)}{{}_{C1}{}^{W}T\left(1,1\right)}\right)$$

(36)

The Jacobian matrix ${}_{C1}{}^{W}J{}^O$ of the centroid coordinate system C1 with respect to the angular velocity of W is as follows:

$${}_{C1}{}^{W}J{}^O=\left(\begin{array}{cc}\begin{array}{c}\frac{\partial {\alpha }_1}{\partial q_1}\\ \frac{\partial {\beta }_1}{\partial q_1}\\ \frac{\partial {\gamma }_1}{\partial q_1}\end{array}& \begin{array}{c}\frac{\partial {\alpha }_1}{\partial q_2}\\ \frac{\partial {\beta }_1}{\partial q_2}\\ \frac{\partial {\gamma }_1}{\partial q_2}\end{array}\end{array}\right)$$

(37)

Therefore, the angular velocity ${}_{C1}{}^W\omega$ of the centroid coordinate system C1 with respect to W is:

$${}_{C1}{}^W\omega =\left({}_{C1}{}^{W}J{}^O\right) \dot{\mathrm{q}}$$

(38)

Based on the pose matrix ${}_{C1}{}^{W}T$, the position (x₁, y₁, z₁) of the centroid coordinate system C1 with respect to W is:

$$x_1={}_{C1}{}^{W}T\left(1,4\right),y_1={}_{C1}{}^{W}T\left(2,4\right),z_1={}_{C1}{}^{W}T\left(3,4\right)$$

(39)

The Jacobian matrix ${}_{C1}{}^{W}J{}^P$ of the centroid coordinate system C1 with respect to the linear velocity of W is:

$${}_{C1}{}^{W}J{}^P=\left(\begin{array}{cc}\frac{\partial x_1}{\partial q_1}& \frac{\partial x_1}{\partial q_2}\\ \frac{\partial y_1}{\partial q_1}& \frac{\partial y_1}{\partial q_2}\\ \frac{\partial z_1}{\partial q_1}& \frac{\partial z_1}{\partial q_2}\end{array}\right)$$

(40)

Therefore, the angular velocity ${}_{C1}{}^{W}v$ of the centroid coordinate system C1 with respect to W is:

$${}_{C1}{}^{W}v=\left({}_{C1}{}^{W}J{}^P\right) \dot{\mathrm{q}}$$

(41)

The inertia matrix ${}_{C1}{}^{C1}I$ of the load with respect to the centroid coordinate system C1 is:

$${}_{C1}{}^{C1}I=\left(\begin{array}{ccc}{}^{1}I_{XX}& -{}^{1}I_{XY}& -{}^{1}I_{XZ}\\ -{}^{1}I_{XY}& {}^{1}I_{YY}& -{}^{1}I_{YZ}\\ -{}^{1}I_{XZ}& -{}^{1}I_{YZ}& {}^{1}I_{ZZ}\end{array}\right)$$

(42)

The inertia matrix ${}_{C1}{}^{W}I$ of the center-of-mass coordinate system C1 with respect to W is:

$${}_{C1}{}^{W}I=\left({}_{C1}{}^{W}R\right)\left({}_{C1}{}^{C1}I\right){\left({}_{C1}{}^{W}R\right)}^T$$

(43)

The kinetic energy ${}^{W}T_1$ generated by the movement and rotation during the lifting process of the left arm of the exoskeleton robot is:

$${}^{W}T_1=0.5{\left({}_{C1}{}^W\omega \right)}^T\left({}_{C1}{}^{W}I\right)\left({}_{C1}{}^W\omega \right)+0.5M_1{\left({}_{C1}{}^{W}v\right)}^T\left({}_{C1}{}^{W}v\right)$$

(44)

The gravitational potential energy ${}^{W}U_1$ is:

$${}^{W}U{}_1=M_{1}g{y}_1$$

(45)

Therefore, the total kinetic energy T of the exoskeleton robot system is:

$$T={}^{W}T{}_1+{}^{W}T{}_2$$

(46)

The total potential energy U is:

$$U={}^{W}U{}_1+{}^{W}U{}_2$$

(47)

Based on the Lagrange function, the total energy L is:

$$L=T-U$$

(48)

The Lagrange equation is expressed as follows:

$$\begin{array}{l}\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)-\frac{\partial L}{\partial q}={}^a\tau +{}^p\tau \\ \frac{d}{dt}\left(\frac{\partial \left(T-U\right)}{\partial \dot{q}}\right)-\frac{\partial \left(T-U\right)}{\partial q}={}^a\tau +{}^p\tau \\ \frac{\partial }{\partial {\dot{q}}^T}\left(\frac{\partial T}{\partial \dot{q}}\right)\ddot{q}+\frac{\partial }{\partial q^T}\left(\frac{\partial T}{\partial \dot{q}}\right)\dot{q}-\frac{\partial T}{\partial q}+\frac{\partial U}{\partial q}={}^a\tau +{}^p\tau \end{array}$$

(49)

Therefore, it can be concluded that:

$$M\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)+G\left(q\right)={}^a\tau +p_{\tau }$$

(50)

where $M\left(q\right),C\left(q,\dot{q}\right),G\left(q\right)$ is:

$$\begin{array}{l}M\left(q\right)=\frac{\partial }{\partial {\dot{q}}^T}\left(\frac{\partial T}{\partial \dot{q}}\right)=\left(\begin{array}{cc}\frac{{\partial }^{2}T}{\partial {\dot{q}}_1\partial {\dot{q}}_1}& \frac{{\partial }^{2}T}{\partial {\dot{q}}_2\partial {\dot{q}}_1}\\ \frac{{\partial }^{2}T}{\partial {\dot{q}}_1\partial {\dot{q}}_2}& \frac{{\partial }^{2}T}{\partial {\dot{q}}_2\partial {\dot{q}}_2}\end{array}\right)\\ C\left(q,\dot{q}\right)=\frac{\partial }{\partial q^T}\left(\frac{\partial T}{\partial \dot{q}}\right)\dot{q}-\frac{\partial T}{\partial q}=\left(\begin{array}{cc}\frac{{\partial }^{2}T}{\partial q_1\partial {\dot{q}}_1}& \frac{{\partial }^{2}T}{\partial q_2\partial {\dot{q}}_1}\\ \frac{{\partial }^{2}T}{\partial q_1\partial {\dot{q}}_2}& \frac{{\partial }^{2}T}{\partial q_2\partial {\dot{q}}_2}\end{array}\right)\left(\begin{array}{c}{\dot{q}}_1\\ {\dot{q}}_2\end{array}\right)-\left(\begin{array}{c}\frac{\partial T}{\partial q_1}\\ \frac{\partial T}{\partial q_2}\end{array}\right)\\ G\left(q\right)=\frac{\partial U}{\partial q}=\left(\begin{array}{c}\frac{\partial U}{\partial q_1}\\ \frac{\partial U}{\partial q_2}\end{array}\right)\end{array}$$

(51)

3.3. Establishment of Control Law for Upper Limb Exoskeleton

In order to complete the power-assisted handling task of the Upper Limb Exoskeleton, the elbow joint angle, angular velocity and angular acceleration are subject to target trajectory planning as the reference input of the control system. In this paper, the target angle process is constructed by the fifth-order polynomial, the total time is 5 s, the control period is 0.01 s, and the lifting angle changes from 0 to 0.3 rad. Taking the left arm as an example, the target angles of the left and right arms are shown in Figure 9. The fifth polynomial is used as the target movement trajectory for the following reasons: (1) The time and angle of conformation to human-body handling process habits. (2) The fifth polynomial can ensure a continuous change in the trajectory angle process, and the corresponding angular velocity and angular acceleration are also continuous changes, which can ensure safety and comfort. (3) It can be set to specify the initial and end angle, angular speed, and angular acceleration. Therefore, the fifth polynomial is a good target trajectory function for complex human dynamics.

In the process of the Upper Limb Exoskeleton carrying objects, the pressure value of the pneumatic muscle is controlled according to the change in the load and the angle error value to maintain the tracking of the elbow angle. Therefore, it is necessary to establish a control rate according to the difference between the actual angle, angular velocity, and angular acceleration and the target value during the rotation process. In order to realize the tracking process of the angle of the handling process, as long as $\ddot{q}\to {}^g\ddot{q}$, there must be $q\to {}^{g}q,\dot{q}\to {}^g\dot{q}$; according to the previous dynamic model, it can be determined that:

$${}^a\tau =M\left(q\right)\times \ddot{q}+C\left(q,\dot{q}\right)+G\left(q\right)-p_{\tau },\ddot{q}\to {}^g\ddot{q}$$

(52)

The error function is constructed as follows:

$$e={}^{g}q-q,\dot{e}={}^g\dot{q}-\dot{q},\ddot{e}={}^g\ddot{q}-\ddot{q}$$

(53)

When $K_d,K_p$ is a positive diagonal matrix, the equation $\ddot{e}+K_d\dot{e}+K_{p}e=0$ is asymptotically stable in a large range. Therefore, when $K_d,K_p$ is a positive diagonal matrix, $\ddot{e}\to 0,\dot{e}\to 0,e\to 0$. The following formula can be obtained by substituting $\ddot{e}={}^g\ddot{q}-\ddot{q}$ into the dynamic model. It is also proven that the control method can realize tracking control and the tracking process is convergent.

$${}^a\tau =M\left(q\right)\left({}^g\ddot{q}+K_d\dot{e}+K_{p}e\right)\ddot{q}+C\left(q,\dot{q}\right)+G\left(q\right)-{}^p\tau$$

(54)

The control torque of pneumatic muscle obtained via discretization is:

$${}^a\tau {}_k=M\left(q_{k-1}\right)\left({}^g\ddot{q}{}_k+K_d{\dot{e}}_{k-1}+K_p{e}_{k-1}\right)+C\left(q_{\mathrm{k}-1},{\dot{q}}_{k-1}\right)+G\left(q_{k-1}\right)-{}^p\tau _{k-1}$$

(55)

According to the principle of torque equivalence, the torque required by the pneumatic muscle can be converted into the output force of the pneumatic muscle. Combined with the length of the pneumatic muscle, the corresponding air pressure value of PAM can be obtained according to the inverse model of the output force of the pneumatic muscle.

3.4. Control System of Upper Limb Exoskeleton Handling Process

The control system in this paper is a rigid flexible coupling system of the pneumatic muscle and Upper Limb Exoskeleton. The air pressure value obtained by the control law is input to the pneumatic muscle, and the length of the pneumatic muscle is calculated using the actual value of the elbow angle. The output force of the pneumatic muscle can be obtained from the length and air pressure of the pneumatic muscle according to the output force formula of the pneumatic muscle. The torque at this time is obtained using the product of the pneumatic muscle output force and vector diameter. This torque acts on the Upper Limb Exoskeleton system to drive the object to move.

When the object is moving, the torque value provided by the load can be obtained according to the principle of torque equivalence. The torque of the elbow joint can be obtained from the difference between the actual value and the target value of the elbow joint angle, angular velocity, and angular acceleration through inverse dynamics calculation of the robot. The difference between the torque and the torque provided by the load is the compensation torque acting on the pneumatic muscle, that is, the output value of the control law. Then, according to the torque equivalent principle of force, the output force acting on the pneumatic muscle is obtained. According to the inverse model of the output force of the pneumatic muscle built using the BP neural network, the control air pressure of the pneumatic muscle can be obtained, which acts on the pneumatic muscle. A control system flow chart of the Upper Limb Exoskeleton handling process is shown in Figure 10.

The physical control process of assisting the Upper Limb Exoskeleton is shown in Figure 11. Firstly, the elbow joint angle signal, angular velocity signal, and angular acceleration signal are collected through the Arduino controller. Three signals are input into the controller, and the difference between the actual signal and the target signal is calculated using the control rate to obtain the control air pressure acting on the pneumatic muscle. The control air pressure combined with the length of the pneumatic muscle acts on the control system, and the force output of the pneumatic muscle is calculated to drive the assistance of the Upper Limb Exoskeleton in transportation. In order to provide better power to the system, a 24 V DC power supply is used. The voltage of the Arduino controller is 5 V, and the voltage is set to a safe and suitable level through the step-down module and voltage stabilization module, providing power support for the system components and driving circuits. The solenoid valve is a normally closed MATRIX pneumatic solenoid valve, which adjusts the charging and discharging speed of the pneumatic muscle by adjusting the switching frequency of the charging and discharging ports of the on/off valve. The maximum frequency of operation of the MATRIX pneumatic solenoid valve is 300 Hz; the driving current is 0.7 A, and the supply voltage is 12 V. The air pump provides an air source to the pneumatic muscle, and three states, including the inflation, deflation, and pressure maintenance of the pneumatic muscle, can be achieved through energization or de-energization of the valve body of the electromagnetic valve.

4. Results and Discussion

In this section, firstly, a process analysis of the change results of the intermediate variables, such as the output torque of the pneumatic muscle and the torque of the elbow joint of the Upper Limb Exoskeletal Robot during the handling process, is carried out. Finally, the tracking results of the actual value and target value of the angle, angular velocity, and angular acceleration of the Upper Limb Exoskeleton handling process are analyzed and discussed. Finally, the error value proves that the control system designed in this paper is effective in the handling process, and the accuracy of the angular velocity and angular acceleration tracking process is good.

4.1. Model Accuracy Evaluation Method

In order to evaluate the accuracy of the pneumatic muscle hysteresis model and the position tracking control of the Upper Limb Exoskeletal Robot, three evaluation criteria are introduced for quantitative analysis, which are the maximum error, the mean absolute error (MAE), the root mean square error (RMSE), and the integral of time and absolute error (ITAE), and the calculation formulas are as follows:

$$\mathrm{MAE}=\frac{1}{H}{\displaystyle \sum _{j=1}^H\left(l_j-{\tilde{l}}_j\right)}$$

(56)

$$\mathrm{RMSE}=\sqrt{\frac{1}{H}{\displaystyle \sum _{j=1}^H{\left(\left(l_j-{\tilde{l}}_j\right)\right)}^2}}$$

(57)

$$\mathrm{ITAE}={\displaystyle {\int }_0^{t}t}|e(t)|dt$$

(58)

where H is the total number of data, l_j $l_j$ is the actual value, and ${\tilde{l}}_{\mathrm{j}}$ is the value obtained by the FN-QUPI model. e(t) is the error at t seconds.

4.2. Discussion of Intermediate Variables in the Process of Upper Limb Exoskeleton Transportation

(1) Load Variation Analysis.

In the process of carrying objects using the Upper Limb Exoskeleton, the interaction force and friction force between the load and the Upper Limb Exoskeleton will change with a change in the elbow joint angle. In order to maintain the balance of objects, the torque provided by the elbow joint will also change. According to the right-hand rule, the torque is negative, as shown in Figure 12a. At the initial stage of handling, the torque provided by the elbow joint fluctuates slightly and changes gently. After 2 s, the torque provided by the elbow joint increases approximately in direct proportion with the increase in the angle. After 4 s, the torque changes tend to be gentle and tend toward a certain value, and the motion curve is smooth, indicating that the torque fluctuation during the movement of the elbow joint is small and has no excessive impact on the arm. According to the principle of torque equivalence and the inverse model of pneumatic muscle output force, the torque of the elbow joint is converted into the pneumatic muscle pressure value. As shown in Figure 12b, the pneumatic muscle pressure value and the torque of the elbow joint show similar changes.

(2) Kinematic and Dynamic Simulation Analysis.

In the process of carrying the Upper Limb Exoskeleton, the control object is the pneumatic muscle pressure. The change in pneumatic muscle pressure is transformed into the length change of the pneumatic muscle, and then, the angle of the elbow joint is changed to drive the movement of the elbow joint. The change in the length of the pneumatic muscle during the assistance process is shown in Figure 13a. At the beginning, the pneumatic muscle is at its original length, so the change is 0. Then, with the process of handling, the pneumatic muscle begins to compress, and the compression changes gently at the beginning. After 0.5 s, the change speed increases, and the speed reaches the maximum at 2.5 s. After that, the approximate symmetrical speed change table is small, and tends to be flat after 4.5 s. It is proven that when the elbow joint angle is changed by the quintic polynomial, the aerodynamic muscle shape variable also changes similarly. The length change process of the pneumatic muscle can be obtained by calculating the sum of the shape variable and the original length of the pneumatic muscle, as shown in Figure 13b. The change speed and acceleration of pneumatic muscle length are shown in Figure 13c,d.

The output of the control system in this paper is the control torque that needs to be compensated by the elbow joint in the process of carrying the Upper Limb Exoskeleton, and its change within the control time is shown in Figure 14. It can be seen from the control rate that the elbow joint torque needs to be compensated according to the difference between the target angle and the actual angle. The image is a smooth curve, which proves the compensation performance of the control system. The control torque of the elbow joint changes gently at the initial time, because the load is predicted at the initial time, the elbow joint provides the corresponding torque, and the control torque also changes gently due to the gentle change in the elbow joint angle, which changes rapidly after 1.5 s and tends to be constant after 4.5 s. The corresponding pneumatic muscle control air pressure process according to the principle of torque equivalence and the inverse model of pneumatic muscle output force is shown in the figure. At the beginning, because of the load prediction, the elbow joint is offset by the load force brought about by the power, so the initial control air pressure value is 0. After that, it changes gently, changes rapidly after 1.5 s, and tends to be constant after 4.5 s.

4.3. Analysis of Trajectory Tracking Results

The fitting results of the actual elbow angle, angular velocity, and angular acceleration with the target value and the corresponding time-varying error are shown in Figure 15.

Figure 15a shows the fitting curve of the actual angle of the elbow joint of the Upper Limb Exoskeletal Robot and the target angle. From the trajectory curve, it can be seen that the elbow joint movement is relatively stable and smooth, and the corresponding real-time error diagram is shown in Figure 15b. Figure 15c shows the fitting curve between the actual angular velocity and the target angular velocity of the elbow joint of the Upper Limb Exoskeletal Robot. The initial velocity and end velocity of the elbow joint are both 0. The angular velocity function is continuous. The robot is stable in the process of movement without rigid impact. The corresponding real-time error diagram is shown in Figure 15d. Figure 15e shows the fitting curve between the actual angular acceleration and the target angular acceleration of the elbow joint of the Upper Limb Exoskeletal Robot. The angular acceleration at the initial position and the end position of the joint is 0, and there is no flexible impact. The corresponding real-time error diagram is shown in Figure 15f. During the trajectory tracking of the arm lifting action of the Upper Limb Exoskeleton, the process runs smoothly, with continuous speed and acceleration, and no rigid and flexible impact. At the same time, it can be seen from

Table 2. Error value of tracking control.

	MAX	MAE	RMSE	ITAE
Angle	0.0175	0.0038	0.0048	4.6426
Angular velocity	0.0031	0.0012	0.0013	1.4645
Angular acceleration	0.0049	0.0018	0.0021	2.2687

that the maximum angular error in the tracking process is 0.0175 rad, MAE value is 0.0038 rad, RMSE value is 0.0048 rad, and ITAE value is 4.6426 rad. At the same time, combined with the error value of acceleration and angular acceleration, it can be seen that the control method in this paper has the ability to complete Upper Limb Exoskeleton handling assistance, and the tracking control process has high accuracy.

5. Conclusions

The hysteresis of the output force of the pneumatic muscle affects the tracking control accuracy of the Upper Limb Exoskeletal Robot driven by it. In this paper, based on the traditional pneumatic muscle output force, the influence of elastic force and friction force on PAM’s output force and the empirical coefficient are added to PAM’s pressure, and PAM’s expansion and load, to establish a new output force model, and an inverse model of output force is established using the BP neural network. Finally, the particle swarm optimization algorithm is used to determine the optimal parameters in the model. With the change in posture in the process of handling, the variable load problem in the process of handling is handled according to the principle of torque equivalence. The control system is built based on the inverse model of the exoskeleton robot’s inverse dynamics combined with the output force of the pneumatic muscle, and elbow angle tracking control in the process of Upper Limb Exoskeleton handling is realized. The results of the mean absolute error and root mean square error show that the control method can achieve high accuracy in elbow angle tracking control of the Upper Limb Exoskeleton system during the handling process. This research result is not only of great help to PAM in the field of assistance and rehabilitation, but also of reference significance regarding the control of flexible actuator exoskeleton robots. However, the accuracy of the proposed control method is largely dependent on the precision of the inverse model of the pneumatic muscle output force, and the calculation of the load depends on the elbow angle and the variation in the angle, which requires a more accurate sensor. In this paper, a BP neural network is used to calculate the inverse output force model, which has no generalization. Therefore, future work will also focus on these limitations.