Modeling Inverse Airflow Dynamics Toward Fast Movement Generation Using Pneumatic Artificial Muscle With Long Air Tubes

This study entailed the derivation of an inverse airflow dynamics model that causes pneumatic transmission delay due to the long air tube in pneumatic artificial muscle control. The inverse model of the derived airflow dynamics was used in designing a feedforward pressure control command. We also modulated the feedback command based on the proposed pressure control for precise force tracking performance. The proposed methods were validated by varying the tube lengths in the pressure tracking task and force tracking task. With the proposed method, the pressure tracking task to the sinusoidal desired profile of 2.5 Hz were successfully achieved and the RMSE decreased 63$\%$ under a certain condition.

The air valve is often placed at a considerable distance from the PAM to reduce the inertia of a robot system and end-effector weight [7], [11].In experimental setups thin and soft bendable tubes are practical [10].These tubes cause pressure transmission delays that degrade the control performance [12].
However, how to compensate for the delay in PAM control has not been carefully studied.
This study derives airflow dynamics model to compensate for pneumatic transmission delay due to the long air tube in the PAM system.The inverse model of the derived dynamics is used to design feedforward pressure control commands for a proportional pressure control valve.We also modulated the feedback force commands to achieve precise force tracking performance with long air tubes.The performance of the proposed controller was validated with different lengths of tubes.The contributions of this article are as follows.
1) Derives an inverse airflow dynamics model for feedforward PAM pressure control.2) Develops a pressure feedforward control approach to improve the force control performance.3) Empirically validates the stability of the proposed method.

II. RELATED WORKS
In control of pneumatic systems, it is common to measure the pressure of the actuator and design feedback command to the valve [13], [14], [15].However, many robot applications that need to use long air tube do not allow us to put a pressure sensor on the actuator side.Therefore, in this study, we aim to control the pneumatic control system with long transmission lines in which the feedback based on the measured pressure inside the actuator is not feasible.
Chou et al. attempted to model pressure transmission as a relationship between the inlet and outlet pressures of the tube [16].Tagami et al. controlled the inner pressure of PAM, considering the pressure transmission [12].However, only the pressure loss caused by the friction from the tube's wall was accounted for in the model.In other words, only the simple manually tuned first-order transmission delay was considered.Krichel and Sawadony [17] divided the tube model into multiple finite-volume elements and derived a discrete dynamics formulation.Turkseven and Ueda [18] constructed the force control method of a pneumatic cylinder with a long transmission  line using Krichel's model.However, due to the time-varying relationship between force and pressure in PAM system, the control approach proposed in Turkseven and Ueda's [18] work cannot be directly applied to PAMs.To control PAMs, we need to explicitly take the time-varying contraction rate of PAM into account since the pressure-to-force relationship depends on the contraction rate.In addition, the previous studies have focused on pneumatic cylinder control and a method to handle the difficulty of the time-varying pressure-to-force relationship has not been sufficiently explored.
To cope with these difficulty we use the inner-loop pressure controller rather than the inner-loop force controller adopted in the previous study [18].In the pressure control domain, we derive the inverse airflow model to compensate the delay with considering the contraction rate.In deriving the inverse airflow model, we adopted the airflow dynamics model proposed in previous studies [17], [18].To the best of our knowledge, this is the first study that shows successful experimental results to compensate delay due to the pneumatic transmission line of PAM system.
Proportional flow control valves are more popular than proportional pressure control valves in the pneumatic field.Although this study primarily focuses on proportional pressure control valves due to the safety and usability.The application of the proposed approach to proportional flow control valves and its control results are also presented in Appendix.

III. MODELING PAM SYSTEM WITH LONG TRANSMISSION LINE
PAM system consists of a PAM, air tube, and control valve, as shown in Fig. 1.This section introduces the PAM model and the airflow dynamics model.

A. PAM Pressure and Mass Flow Rate Dynamics Models
The viscous force acts the air in the opposite direction to the flow.As the viscous friction depends on the flow field [19], deriving an analytical solution to the airflow dynamics is challenging.Therefore, an approximated model of the segmented tube elements has been explored in Krichel and Sawadony's [17] work.Fig. 2 shows our two compartment model.Although the PAM volume V pam depending on the contraction rate is nonlinear and difficult to calculate analytically, it is known to be approximated by a quadratic polynomial of [20].V pam and total end volume V total are expressed as V 0 is the PAM volume under the atmospheric pressure.From the natural length l 0 and contracted length Δl of PAM, is defined as (t) = Δl(t)/l 0 .The mass flow rate dynamics through the tube ṁ is expressed as [17] where P P denotes the PAM pressure, and P V is the pressure provided by the proportional pressure control valve.S is the air tube cross-section area, L is the air tube length, and f is the friction to ṁ.The dynamics of P P is expressed as [21] dP where α 1 , α 2 ∈ [1, 1.4] and 1 is used in isothermal change and 1.4 is used in adiabatic change.Appropriate parameters should be selected between the above-mentioned range.We manually tuned the parameters to achieve the best estimation performance and α 1 = α 2 = 1.2 were selected for setup with the proportional pressure control valve.R = 287 J/(kg • K) is specific gas constant and T = 298 K is room temperature.We adopted the approximated friction model as follows [18]: Re (P V , ṁ, P P ) 3/4 η D 2 Re (P V , ṁ, P P ) = ṁD ηSρ (P V , P P , T ) where D is the tube diameter, η = 1.52 × 10 −5 m 2 /s is kinematic viscosity, Re is Reynolds number, and ρ is the air density.Airflow dynamics consists of (2) and (3).A tube is parameterized by L and S. The effect of the tube length L appears in (2) as a coefficient on pressure and in (3) as a part of V total .

B. Pressure-to-Force Model
The relationship between the force F of the PAM and P P and is highly nonlinear and difficult to model accurately.Hildebrandt et al. expressed it well as (5) using polynomial approximation up to the third order of [14].P atm = 0.1013 MPa is the atmospheric pressure.The PAM model accuracy can be improved by adding damping or friction terms including ˙ [22], [23].However, due to vibration caused by ˙ , it is not expected to improve the control performance with a simple implementation.Moreover, as the effect of these terms is not the focus of this study, a static model was chosen.

IV. CONTROLLER DESIGN
This section introduces the outer-loop force controller and the inner-loop pressure controller, as shown in Fig. 3.

A. Force Controller
In our PAM control framework, the reference force F d for the inner-loop pressure controller was derived as follows: where F d is the desired force, K p and K i are the proportional and integral gains for the force feedback control.The force tracking error e F is derived as e F = F − F d .We derived the inverse model of ( 5) to convert the force reference F r to the reference PAM pressure P r as follows:

B. Pressure Controller Based on the Inverse Airflow Dynamics
We use the inverse dynamics of the models introduced in Section III-A to derive the desired pressure P d .From (3), the reference mass flow rate ṁr is derived as P d is then determined as the inverse model of (2) as P V = P d is achieved by the proportional pressure control valve and the control of the valve is described in the Appendix.Feedforward pressure controller designing P d with P r is composed of ( 8) and (9).Feedback force controller designing P r with F d is composed of ( 7), (8), and (9).

C. Stability Analysis
To analyze the stability of the proposed pressure control method, we evaluated the linearly approximated dynamics of the PAM control system around a state trajectory.We summarize ( 3) and ( 2) as follows: From ( 8), (9), and ( 10), the differential equation of x r = [ ṁr , P r ] T is expressed as By subtracting ( 11) from ( 10) and assuming is obtained.Defining f (x, P V ) as f (x, P V ) = f (x, P V ) ṁ, the Taylor expansion of f (x, P V ) around a certain point x = [ m, PP ] T with the input PV is expressed as Similarly, f (x r , P V ) ṁr is linearized as By substituting ( 13) and ( 14) to (12), the dynamics of e x = x − x r is expressed as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.For (15), the transition matrix Φ A (t, τ ) (t ≥ τ ) can be defined as follows: Then, e x (t) can be derived from e x (τ ) as Here, we discuss the stability of the error dynamics (15) with the transition matrix Φ A , as suggested in Kailath's [24] work.If the condition is satisfied, then the error dynamics ( 15) is considered Lyapunov stable.Furthermore, if the condition is satisfied, then (15) is considered asymptotically stable.

V. EXPERIMENT
We conducted pressure tracking experiments and force tracking experiments to evaluate the proposed method.was mounted on a stage and connected to a proportional pressure control valve (Norgren VP5010B) with an air tube (inner diameter: D = 5 mm).The valve pressure P V was measured by a pressure sensor attached to the exit port of the valve.We also measured P P with a pressure sensor (Honeywell 26PCFFJ6G) connected to the PAM via a thin air tube.The observed value of P P was used only for evaluation and not for control.Δl was measured by a custom-made linear encoder attached at the bottom of a linear slider to derive .F was measured by a load cell (Tec Gihan Co, USL08-H6-3 kN) attached to the PAM tip.The control frequency of the system was 250 Hz.The supply pressure from the compressor (Hitachi SRL-A3.7DW) to the valve was maintained around 0.78 MPa.

B. Controller Settings
For the pressure control task, we compared the following two control methods: 1) Our proposed method using ( 9) and ( 8) to derive P d .2) Baseline method that did not explicitly consider the existence of the long air tube as P d = P r .We derived d ṁr dt by the numerical derivative of ṁr in the proposed method and we applied a first-order low-pass filter 1/(1 + τ s) with a time constant τ = 0.05 to d ṁr dt , to eliminate the high-frequency noise due to the numerical differentiation.
For the force control task, we compared the two control methods: 1) Our proposed method and 2) baseline method.
Feedback gains of each force controller were tuned manually for the sinusoidal desired trajectory with the air tube length of 5 m and frequency ν = 0.5 Hz.The gains of the both method were tuned based on Ziegler-Nichols method [25].Owing to the lack of differential operation, the gains of the baseline method are less prone to oscillations, which allows for higher gains than that proposed as listed in Table I.
The identified parameters used in the controllers are shown in Table II.V 0 was obtained from measurement and the other parameters were identified from the prior experiment with a 4-m tube.For the identification of V 1 and V 2 , the minimization of the square error between measured P P and the estimated value PP was conducted to experimental data when valve pressure was generated as the sinusoidal waves of 1, 2, 3, 4 Hz under the free condition.PP is estimated by ( 2) and (3) with the estimated mass flow rate m.For the identification of d i (i = 0, 1, . .., 4) and b i (i = 1, 2, 3), the minimization of the square error between measured F and the estimated value F was conducted to experimental data, when PAM was extended slowly under isobaric PAM pressure conditions, for each pressure 0.2-0.8MPa in 0.1 MPa increments.

C. Trajectory Tracking Tasks 1) Pressure Tracking Task:
The pressure tracking performance was evaluated based on the sinusoidal trajectory tracking tasks with different bias levels and frequencies: P r (t) = A p sin 2πνt + B p , where the amplitude A p = 0.08 MPa, bias level B p = 0.36, 0.46, 0.56 MPa, and frequencies ν = 0.5, 1.0, 1.5, 2.0, 2.5 Hz.The values of A p and B p were set around the pressure range of the force tracking task in Section V-C2.As in Fig. 5, we conducted the target pressure tracking task in three different conditions as follows.
1) Free condition: one side of the PAM was fixed while the other could freely move on the ball linear guide.2) Fixed condition: both sides of the PAM were fixed.
3) Antagonistic condition: one side of the PAM was fixed, and the other was connected to the antagonistic PAM.The inner pressure of the antagonistic PAM was maintained at a constant value 0.4 MPa.The antagonistic condition emulated the agonistic-antagonistic system composed of two PAMs [26].The pressure tracking tasks at every condition were conducted with a 7 m long tube.In practice, we often used several meter long air tube in our rehabilitation studies [5], [6], [7], [11], [27].To check for length dependence, 3-and 5-m tubes were also used for the experiments under the free condition with B p = 0.46 MPa.We further evaluated the pressure tracking performance to the desired trajectory composed of the two sine functions with a 7-m tube under the free condition: P r = 0.07 sin 2π • 1.3t + 0.05 sin 2π • 2.2t + 0.39 MPa.
2) Force Tracking Task: We then evaluated the force tracking performance with sinusoidal target trajectory: F d = A f sin 2πνt + B f , where the amplitude A f = 100 N , bias levels B f = 400 N and frequency ν = 0.5, 1.0, 1.5, 2.0, 2.5 Hz with 3-, 5-, and 7-m tubes.The force tracking task was conducted under the the antagonistic condition.

A. Pressure Tracking Task
The error between the reference and observed PAM pressure (P r and P P , respectively) was evaluated.Fig. 6 compares the   control performances of the each method with root-mean-square error (RMSE).Fig. 7 shows the results of pressure control at 3-and 5-m tubes.The proposed method improved the control performance under all three conditions.Fig. 8 shows an example of the pressure tracking performances at ν = 2.5 Hz and B p = 0.46 MPa.In this condition, RMSE decreased by 63% from the result of the baseline method.Our proposed method compensates for the delay, while the tracking performance of the baseline approach shows an apparent delay in the reference pressure profile.An example of estimating the value of P P using the proposed airflow dynamics model composed of ( 2) and ( 3) was shown in Fig. 9 for the validation of the proposed model.Fig. 9 shows the behavior of P V and P P along with the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.estimated values of P P during pressure tracking task (7 m tube, ν = 2.5 Hz, and B p = 0.46 MPa).The delay caused by the tube can be observed from the comparison between P V and P P .The proposed model can estimate the delay.Fig. 10 shows the pressure tracking result for the desired trajectory of two sine waves combination.The delay was compensated and even minute desired pressure changes were reproduced by the proposed method.

B. Force Tracking Task
Fig. 11 shows the force tracking performances for a desired sinusoidal trajectory.The control performances of our proposed method were robust against changes of the tube lengths and the frequency, whereas those of the baseline method were not.The superiority of the baseline method over the proposed method was observed due to the higher gain under low frequency and short tube conditions, however the error significantly increased with longer tubes and higher frequency.Fig. 12 shows the force tracking profiles with the 7-m air tube at ν = 2.5 Hz.Our proposed method successfully compensated for the pressure transmission delay.In contrast, an apparent delay was observed in the baseline method result.13.Bode plot of force control results with 7-m air tube.No delay is observed in proposed method, while baseline method was significantly delayed from 1.5 Hz, as shown in phase plot.The bode plots of the force control result using 7 m are shown in Fig. 13.This figure corresponds to the RMSE result with the 7 m in Fig. 11.In the phase plots at 0.5 and 1.0 Hz, there is no significant difference between the values for the both method.However, from 1.5 Hz onwards, a phase delay was observed in the result of the baseline method, and this delay increases with increasing frequency.From the results in Fig. 11, This delay leads to an increase in control errors and the proposed method has compensated this delay.There is no significant difference in the gain plots at all frequencies between the both controllers.

C. Empirical Stability Analysis
Fig. 14 shows the matrix norm of Φ A , i.e., ||Φ A (t, τ )||, when we conducted the pressure tracking task at ν = 2.5 Hz and B p = 0.46 MPa using a 7-m tube under the antagonistic   condition.We calculated the norm ||Φ A (t, τ )|| from each starting time; all values of ||Φ A (t, τ )|| were contained in the light orange region.Therefore, all norms converged to zero.This result empirically shows that the PAM system controlled by our method is asymptotically stable.
||Φ A (t, τ )|| shows different behaviors in Fig. 14 for τ = T 1 and τ = T 2 .To investigate their difference, the vector fields of (12) with A (t = t) are plotted in Fig. 15.At t = T 1 , Re(t = T 1 ) is low and the air flow is characterized by laminar flow.f (t = T 1 ) is also low, because of which the eigenvalue of A (t = T 1 ) contains imaginary numbers and the trajectory converges spirally [see Fig. 15(a)].At t = T 2 , Re(t = T 2 ) is high and the air flow is characterized by turbulent flow.f (t = T 2 ) is high, and the corresponding eigenvalue of A (t = T 2 ) contains no imaginary numbers, making the trajectory converge without vibration [see Fig. 15(b)].

VII. DISCUSSION
The proposed method is effective in scenarios where the influence of the airflow dynamics is significant.However, limitations attributed to the low gains resulting can hinder the performance of force control in conditions where the influence of the airflow dynamics is small, such as small pressure changes, short tubes, and small end-volume changes.In fact, in the force control task with 7-m tube under fixed conditions, as shown in Fig. 16 the performance was not better than that of the baseline method; however, this is an unlikely situation in which PAMs are used.

VIII. CONCLUSION
The long tube between the valve and the PAM worsens the control performance of the PAM system.We derived the airflow dynamics using a long transmission line and designed a controller using the inverse dynamics.The proposed method improved the control performance under the significant effect of the pneumatic transmission line such as high-frequency sinusoidal reference profiles or long tubes.
In this study, we focused on improving control performance of a single PAM system.Joint actuators with antagonistic PAMs were studied [26], [28], [29], with the force control term of a single PAM inside the actuator control system.Therefore, improving the control performance of a single PAM can enhance the performance of the antagonistic joint, which can lead to performance improvements in robotic systems with PAM-based joints.In addition, the pressure feedforward controller compensating the pneumatic transmission delay could be applied to other pneumatic actuators [30], [31] and would allow them to operate in fast range.

APPENDIX PROPOSED METHOD WITH FLOW CONTROL VALVE
The proposed method was additionally applied to a proportional flow control valve to evaluate whether our proposed approach is useful even with the different physical input to the system.the PAM system is represented by ( 2), (3), and where ṁV is the mass flow rate provided by the valve, as shown in Fig. 17.Equations ( 2), (3), and (20) are summarized as The mass flow input ṁd to the system was designed as m is the estimated value of ṁ obtained from (2) and (3).
The proposed method for the proportional flow control valve consists of the components represented by ( 8), (9), and (22).K f is the proportional gain for the valve pressure feedback control and tuned as K f = 3 × 10 2 1/s.ṁV = ṁd is achieved by the proportional flow control valve and the control method of the valve is described in the next appendix section.In this article, we used α 1 = 1.4 and α 1 = 1.2 in (3) for the system with the proportional flow control valve (Festo MPYE-5-1/4-010-B), and reidentified V 1 and V 2 as V 1 = 1.09 × 10 −4 m 3 and V 2 = −2.37 × 10 −4 m 3 , respectively.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The baseline method that did not consider the delay due to the long transmission line used the control command in (22) directly using the reference pressure as the desired pressure P d = P r , as explained in Section V-B.
Fig. 18 shows the RMSE result of the pressure tracking task with 7-m tube under the free condition (see also Fig. 5).Our proposed method showed better control performance than the baseline method even with a proportional flow control valve.

VALVE CONTROLLER SETTINGS
The model of proportional pressure control valve (Norgren VP5010B) and the controller designed accordingly are expressed, respectively, as [27] where τ v is time constant, u bias is the voltage bias, and c p = 1 × 10 −1 MPa/V is the conversion parameter from the voltage to the pressure.The parameters were identified as τ v = 6.15 × 10 −2 1/s, u bias = −3.00× 10 −1 V.For the model of the proportional flow control valve (Festo MPYE-5-1/4-010-B), the mass flow rate is defined as [32] ṁV = (L c + max{(u − U c )b − L c }) φ(P c , P V ) where φ(P a , P b ) is the thin-plate flow function between the pressure P a and P b (for more detail see [32] and [33]).Both valves are controlled by the multifunction board, where the board communicates with real-time OS Xenomai through real-time Ethernet as [7] and [11].

Fig. 2 .
Fig. 2. Two compartment model of transmission line.Air tube is divided into two sections.

Fig. 3 .
Fig. 3. Block diagram of proposed method comprising force feedback and pressure feedforward controllers.

Fig. 4 .
Fig. 4. PAM is mounted on stage and connected to proportional pressure control valve with air tube.Control system is equipped with load cell, pressure sensor, and linear encoder.(a) Experimental setup.(b) Sensor allocations.

Fig. 4
Fig. 4 shows our experimental setup and sensor allocations.The PAM bladder (cut from Festo DMSP-20, diameter: 20 mm, length: 325 mm) plugged with custom-made PAM both ends

Fig. 6 .
Fig. 6.Pressure control performances: RMSE of each pressure control experiment with 7 m length tube.Proposed method improves control accuracy in all conditions.

Fig. 7 .
Fig. 7. Performances of each pressure control at the experiment with 3 and 5 m length tube.Proposed method improves control accuracy even with 3 and 5 m tubes.

Fig. 8 .
Fig. 8. Pressure tracking profiles with 7-m air tube for desired sinusoidal trajectory at ν = 2.5 Hz and B p = 0.46 MPa.Delay caused by long tube is compensated by proposed method.

Fig. 9 .
Fig. 9. P V , P P , and estimated P P by proposed model.

Fig. 10 .
Fig. 10.Pressure tracking profiles for desired trajectory of two sine waves combination.Even minute pressure changes are reproduced by proposed method.

Fig. 11 .
Fig. 11.Force tracking performances for sinusoidal trajectory with three different tube lengths.Proposed method was robust against changes in tube lengths and frequency.

Fig. 12 .
Fig.12.Force tracking profiles with 7-m air tube for desired sinusoidal trajectory at ν = 2.5 Hz.Proposed method can track without delay, while baseline method causes significant delays.
Fig. 15.Vector fields for each A (t). Trajectory converges in spiral at t = T 1 .Trajectory converges without vibration at t = T 2 .A (t) have different dynamic properties due to Re(t).(a)A (t = T 1 ).(t = T 2 ).

Fig. 16 .
Fig. 16.Tracking performances for desired sinusoidal force trajectory with frequency with 7 m under fixed condition.

Fig.
Fig. Input of system is mass flow rate generated by valve.

Fig. 18 .
Fig.18.Pressure control performances with proportional flow control valve.Proposed method showed better control performance than baseline method.
L c and L r are the minimum aperture areas of the supply and exhaust ports.U c and U r are the voltage values when the respective ports are sealed, and b is the conversion factor from voltage to aperture area.By using the inverse model of (25) the feedforward controller of the flow control valve is designed asu = ⎧ ⎨ ⎩ U r + ṁd +L r φ(P V ,P room ) bφ(P V ,P room ) ( ṁd < 0) U c − ṁd +L c φ(P c ,P V ) bφ(P c ,P V ) ( ṁd ≥ 0).(26)We adopted a first-order low-pass filter with a time constant τ = 0.05 to derive ṁd to prevent oscillation.The parameters were identified as U c = 5.51 V, L c = 2.14 × 10 −7 m 2 , U r = 4.65 V, L r = 3.65 × 10 −7 m 2 , b = 4.76 × 10 −6 m 2 /V.