Generalized Multiphase Dynamic Modeling and Precision Interaction Force Control of a Walking Lower Limb Hydraulic Exoskeleton

Wearable lower limb hydraulic exoskeletons can be used to augment the human performance in heavy load transportation. Nonlinear and walking phase-dependent dynamics make the lower limb hydraulic exoskeleton become difficult to be modeled. This paper presents a generalized multiphase dynamic modeling method in which the dynamic model of each walking phase can all be solved based on a general higher dimensional dynamic model and different holonomic constraints. Compared to traditional lower limb exoskeleton modeling methods where the modeling of each walking phase is done independently, the proposed method is simple and applicable to arbitrary walking phases, especially for double leg support phase (closed-chain dynamics). Based on the established dynamic models, MIMO adaptive robust cascade force controllers (ARCFC) are designed both for double leg support phase and single leg support phase to effectively address high-order nonlinearities and various modeling uncertainties in hydraulic exoskeletons. An additional torque allocation method is proposed to deal with the overactuated characteristic in double leg support. Comparative simulations are conducted to verify the excellent human-machine interaction force control performance of the proposed scheme.


Introduction
Wearable lower limb exoskeletons are intelligent humanmachine integrated devices used to augment the performance of the wearer in heavy load transportation [1,2], such as soldier marching, earthquake rescue, and construction site. Thanks to the large ratio of power-to-weight, it is suitable to adopt hydraulic actuators in developing such systems which need to be a small size while providing large force. In hydraulic lower limb exoskeleton, the wearer provides motion commands, while hydraulic actuators supply enough actuation force to support heavy loads. When the exoskeleton tracks the human motion precisely, little load force can be felt by the wearer. Since the human-machine interaction force is closely related to the trajectory tracking error between the wearer and exoskeleton, the control objective can finally be transformed into minimizing the humanmachine interaction force.
As for the interaction force controller design of a walking hydraulic lower limb exoskeleton system, a number of challenging issues exist. First, different from 1-DOF or single leg exoskeleton, there exist multiple walking phases in the walking lower limb exoskeleton, such single leg support phase, double leg support phase, and subphases during double leg support (like toe off and heel strike). Moreover, for double leg support phase and its subphases, the two feet are landing the ground at the same time leading to a closedchain mechanism. Also overactuated conditions exist in closed chain dynamics. All these make the dynamic modeling and controller design of lower limb exoskeleton become more difficult. Besides, as a nondesired force output source, high-order nonlinearities and various model uncertainties exist in hydraulic actuators which make the controller design of hydraulic lower limb exoskeleton more challenging than that of motor-driven system.
As for the multiphase dynamic modeling of lower limb exoskeletons, usually the dynamic model of each walking phase is established separately. For single leg support phase, since it is a serial-chain mechanism, the dynamic model is usually obtained straightforwardly using Lagrangian equations. In [3], the lower limb exoskeleton in single leg support is described as serial chain of 7 segments. Using Lagrangian equations of the second kind, the planar model of single leg support is obtained. In [4], the Lagrange equations of the second kind is adopted to model the single leg support phase of 6-link biped robot. For double leg support phase or its subphases (such as toe off and heel strike), since it is a closed-chain mechanism with overactuation, the Lagrangian of the second kind cannot be applied directly. In [3], the lower limb exoskeleton in double leg support is partitioned into two 3-DOF serial manipulators. Lagrangian equations of the second kind are used to model the dynamics of each 3-DOF manipulator, while Newtonian mechanics is adopted to describe the relationship of the two parts. In [4], the Lagrange formulation of the first kind using Lagrange multipliers is adopted for dynamic modeling of biped robot in double leg support. However, all these methods focus on the dynamic modeling of each walking phase independently in which different general coordinates needed to be established and different Lagrangian modeling process need to be conducted for each walking phase. This obviously leads to a complicated modeling process when there exist various walking phases. Also, the joint positions may become discontinuous due to the role switching of the swing and support leg [5].
As for the human-machine interaction force controller design, various control schemes are proposed [6]. One method is to design the exoskeleton controller without measuring the human-machine interaction force. In [7], an inverse dynamics-based positive feedback controller is designed for the Berkeley lower extremity exoskeleton so that the sensitivity to the human force can be increased. In [8], first the inverse dynamics of the exoskeleton is adopted to estimate the joint torque, and then, a fictitious gain is designed to increase the sensitivity of the human body. However, these methods ignore the model uncertainties in computing the inverse dynamics leading to a poor robust performances. Another line of thought is to design a cascade force controller based on the measured human-machine interaction force. Model-free PID controller [9], admittance controller [10], and human-machine cooperation controller [9] are often used in the high-level control algorithm design to generate the desired motion command, while PID [11] or dynamic model compensation are often used in the low-level controller design to achieve the motion tracking. In these cascade force control schemes, the human motion intent is inferred from the fixed dynamic model or human data which cannot be generalized to different wearers. Also, the proposed low-level motion tracking algorithm cannot guarantee the fast response and accurate tracking of human motion in the presence of strongly coupled nonlinearities and various model uncertainties. Other control methods have considered the model uncertainties in the controller design, such as the adaptive impedance control [12], neural network control [13], sliding mode controller [14,15], and robust output feedback-assistive control [16]. However, these methods are mainly focused on the controller design of motor-driven exoskeletons which cannot be easily extended to the control of hydraulic exoskeletons. The reason is that the order of hydraulic exoskeleton dynamics is higher (it is a three-order dynamics for hydraulic exoskeleton in this paper, while it is usually a two-order dynamics for motor-driven exoskeleton systems). Also, various model uncertainties exist in the hydraulic actuators, like the load change, hydraulic parameters variation (e.g., bulk modulus), external disturbances, and leakage.
Recently, an adaptive robust cascade force control algorithm is proposed for 1-DOF and 3-DOF hydraulic exoskeleton system [17,18]. Robust interaction force control performance to various model uncertainties is guaranteed. In this paper, the problem is extended to a walking lower limb exoskeleton. Compared to the 3-DOF single leg exoskeleton, multiple walking phases exist when two legs move, which makes the human motion intent inference method of lower limb exoskeleton different from that of single leg exoskeleton. Moreover, the closed-chain dynamics and overactuated characteristic in double leg support and its subphases make the interaction force controller design become more challenged.
This paper presents a generalized multiphase dynamic modeling method and a robust interaction force control scheme for hydraulic lower limb exoskeleton. The principal contributions can be summarized as follows: (1) A generalized multiphase dynamic modeling method is proposed for lower limb exoskeleton which is applicable to arbitrary walking phase, especially for double leg support phase (closed-chain dynamics).
Since the generalized coordinates are the same and the Lagrangian modeling process will only be done once, the multiphase dynamics modeling process of a walking lower limb exoskeleton becomes simple (2) Multiphase adaptive robust interaction force controllers are designed to effectively deal with overactuated characteristic, negative effect of high-order nonlinearities of hydraulic system, various parameter uncertainties, and modeling errors. Robust interaction force control performance is achieved

Physical Modeling
Since human dynamics is very complicated and may not be conveniently applied for exoskeleton controller design, in  Figure 1.
For a floating exoskeleton in which positions of the two exoskeleton feet can be changed, as shown in Figure 2, we can define nine generalized coordinates to completely describe the positions of such system. Denote where x l , y l , and q l represent the left foot positions in the Cartesian coordinate frame, q 1 ∼ q 6 represent the joint positions in the left leg and the right leg, τ 1 ∼ τ 6 represent the joint torque from the actuators, and F L and F R represent ground contact force at the left and right foot. Using Lagrangian equations, the dynamic model of the above floating exoskeleton can be described as where MðqÞ, Cðq, _ qÞ, and GðqÞ represent the system matrices and gravity force, B a represents the joint-torque projection matrix, and J L and J R represent the Jacobian matrix in the left and right foot.

Dynamic
Model for Each Walking Phase. In different walking phases, the contact condition between the exoskeleton foot and the ground is different which leads to different holonomic constraints. Combining the same general dynamic model (Equation (2)) with different holonomic constraints, we can finally solve the detailed dynamic model for each walking phase through solving amount of linear equations. It should be noted that the generalized coordinates and the general high dimensional dynamic model are the same for the modeling of each walking phase which are quite different from the traditional multiphase dynamic modeling method where the generalized coordinates will be redefined and different Lagrangian modeling process will be conducted in each walking phase.
Define the following terms: (1) Left Leg Support Dynamics. For left leg support, as shown in Figure 1(a), the positions and the rotation angle of the left foot are fixed. Thus, the following holonomic constraints exist where C xl , C yl , and C ql are constant values. Besides, since the right leg is the swing one, there is no ground contact force in the right foot, namely 3 Applied Bionics and Biomechanics Differentiating (4) while noting (2) and (5), we can obtain Here, there are 15 unknown variables (€ q, F L , and F R ) and 15 equations (Equation (6)). Using T act as input, € q, F L , and F R can be solved. Finally, the dynamics can be obtained as where M Lsp , C Lsp , and G Lsp are matrices computed from M, C, and G, as follows: (2) Right Leg Support Dynamics. For right leg support, the positions of the right foot are fixed. Thus, the following holonomic constraints exist: where C xr , C yr , and C qr are constant values. Besides, since the left leg is the swing one, there is no ground contact force in the left foot, namely Since _ Differentiating (9) while noting (11) and (2), we can obtain Here, there are 15 unknown variables (€ q, F L , and F R ) and 15 equations (Equation (12)). Treating T act as input, € q, F L , and F R can be computed. Finally, the dynamics can be obtained as where M Rsp , C Rsp , G Rsp , M Rsp2 , C Rsp2 , and G Rsp2 are matrices computed from M, C, and G, as follows:   Applied Bionics and Biomechanics (3) Double Leg Support Dynamics. For double leg support, both the positions of the right foot and the left foot are fixed. Thus, the following holonomic constraints exist: x l = C xl , x r = C xr , Differentiating (15) while noting (11) and (2), we can obtain Here, we have 15 unknown variables (€ q, F L , and F R ) and 15 equations (Equation (16)). Treating T act as input, € q, F L , and F R can be computed.
Finally, the dynamics can be obtained as where all the new matrices in (17) can be computed from M, C, and G. Define the following terms: M Dsp , C Dsp , G Dsp can be computed as follows: (4) Right Heel Strike Dynamics. For right heel strike walking phase, as shown in Figure 1(b), the positions and the rotation angle of the left foot are fixed. Also, the positions of the right foot heel are fixed. Thus, the following holonomic constraints exist: x l = C xl , x r = C xr , where C xr and C yr are constant values. Besides, since the right leg contact the ground on a point, there is no ground contact torque in the right foot, namely Differentiating (20) while noting (2) and (11), we can obtain Here, there are 15 unknown variables (€ q, F L , and F R ) and 15 equations (Equation (22)). Treating T act as input, € q, F L , and F R can be solved. Denote Finally, the dynamics can be obtained as where all the new matrices in (24) can be computed from M, C, and G. Define the following terms: M Hs , C Hs , and G Hs can be computed as follows:  Figure 3(d), the positions and the rotation angle of the right foot are fixed. Also, the positions of the left foot toe are fixed. The holonomic constraints are similar to the right heel strike walking phase. Following the same computation process as that of right heel strike walking phase, the dynamic model of toe off walking phase can be obtained. The detailed computation process is omitted here for simplicity.
2.2. Hydraulic Actuator Dynamics. The pressure and flow rate dynamics of hydraulic cylinders can be described as [18] where the definition of all the terms can be seen in the notation list at the Appendix.

Human-Machine
Interface Dynamics. The model of human-machine interaction force may be complex; for example, the model may be uncertain and/or varies. Thus, it is difficult to obtain a precise model that can describe the properties of actual human/robot attachment precisely. Also, a much precise human/robot interface dynamic model may be much complex to be used for designing the controller. Thus, in the paper, only the main property of the interface is considered. Since a belt is often adopted in connecting the robotic leg with the human leg. Thus, a spring model with unknown stiffness can describe the main compliant property of the interface. As for other unmodeled uncertainties, we put them in the lumped disturbance. In the later part, an adaptive robust control algorithm is proposed to address the modeling errors in human-machine interface dynamic model. Thus, the human-machine interface dynamic model is described as a spring with lumped disturbance and unknown stiffness: where the definition of all the terms can be seen in the notation list at the Appendix. Equation (28) is algebraic and devoid of dynamics. Using the integral of interaction force Ð t 0 F hm dτ in the controller design, then the following dynamic model can be obtained: 3. Human-Machine Interaction Force Control Schemes 3.1. Control Objective. Assuming that the wearer is able to achieve balance and locomotion, the control objective of a walking lower limb exoskeleton is to design a control law according to the multi-phase dynamic models that minimizes the human-machine interaction force so that accurate human motion tracking is achieved.
3.2. Overall Control Structure. Figure 3(a) shows the overall control structure. Foot switches are mounted at the exoskeleton foot to recognize whether the exoskeleton foot contact the ground or not, which can further help us recognize which walking phase the exoskeleton lies in. Six-axis force sensors are fixed at the back and the foot of the lower limb exoskeleton to measure the interaction force at these places. For different walking phase, different interaction force components are selected for force controller design. Many control techniques have been developed for hydraulic systems [19][20][21]. The ARC is verified to be effective in addressing various model uncertainties through lots of practical applications [22][23][24][25][26], especially for multijoint hydraulic manipulator [27,28]. Thus, it is adopted in the following controller design. Therefore, it will be adopted in our interaction force control algorithm design. For each walking phase, an adaptive robust cascade force controller (ARCFC) is designed so that a good robust force control performance can be achieved, as shown in Figure 3(b).

ARCFC Design for Double Leg Support.
For the walking phase of double leg support, there are three independent degree of freedoms, which can be seen from the dynamic model (17). Thus, in double leg support, we can only control three independent human-machine interaction force components. Since the exoskeleton feet are always in flat contact with the ground and will not hinder the movement of human feet, there is no need to reduce the interaction force at the foot contact points. Finally, three human-machine interaction force components at the back are minimized.
Since there are six control inputs, the system is overactuated.
In the later force controller design, three independent virtual control torques are figured out first, and then, a torque allocation method is proposed to specify the desired load force for six hydraulic cylinders. Finally, six control inputs of the hydraulic valves can be figured out.
The overall system dynamics for double leg support can be described by Applied Bionics and Biomechanics where K = diag fK ubx , K uby , K ubz g is the stiffness vector of the human-machine interface and B = diag fB 1 , B 2 , B 3 g is the damping ratio. The definitions of all the other terms can be seen in the notation list at the Appendix. The fifth equation of (30) has three properties: where _ q r and € q r are any vector and β is the parameter vector of the exoskeleton.
Define the following unknown parameters and lumped disturbances:Δ

Applied Bionics and Biomechanics
where Δ in and Δ i are the constant and the time-varying part ofΔ i . We assume that the bound of the uncertain parameters and disturbance is known. Define the following states: x 4 = P 1 = P 11 P 12 P 13 P 14 P 15 P 16 ½ T , x 5 = P 2 = P 21 P 22 P 23 P 24 P 25 P 26 ½ T ,

Applied Bionics and Biomechanics
The state space expression of dynamics (30) can be described as where 10 Applied Bionics and Biomechanics The control goal is to synthesize a control input u = u 1 u 2 u 3 u 4 u 5 u 6 ½ T based on (34) that minimizes the human-machine interaction force F hmub .

High-Level Human Motion Intent Inference. Treating
x eub as virtual control input, then a control law making the integral of force tracking error z 1 = x 1 − x 1d converge to zero or to be bounded is designed as where x ma is the model compensation term, x ms is the robust feedback term, K 1 = diag fK 11 , K 12 , K 13 g is the linear feedback gain, Γ 1 > 0 is a diagonal adaptation rate matrix, x msn is a nonlinear feedback item, We can treat x m as inferred human motion intent, and then, the desired joint positions can be obtained as Let z 2h = Kineðx 21 Þ − x m , where Kine means kinematics. Then, the first error dynamics is given as Similar to [18], in order to obtain the desired motion tra- c1m ), a output differential observer is adopted.

Low-Level MIMO Motion Tracking Controller.
In lowlevel controller design, a motion tracking control algorithm making the position tracking error z 2 = x 21 −q c1m converge to zero or to be bounded is proposed with the following design procedures.
Step 1. Specify the desired torque τ actd for B Dsp hP L that achieves accurate motion tracking (i.e., x 21 ⟶q c1m if B Dsp hP L = τ actd ).
Treating B Dsp hP L as the virtual control input in this part, the control law τ actd is given as τ actd = τ actda + τ actds ,

Applied Bionics and Biomechanics
where τ actda is the term for model compensation, τ actds is the robust feedback item, K 3s1 is the gain for linear feedback, K 3 > 0, g 3 > 0, Γ 2 > 0, and τ actdsn is a term for nonlinear feedback. Bx 31 = Y B ðx 31 ÞB θ .
Step 2. Torque allocation (i.e., specify the desired load force P Ld for each hydraulic cylinder such that B Dsp hP Ld = τ actd ).
With τ actd given in Step 4, the next task is to figure out the desired load force for each hydraulic cylinder such that the combined effort equals τ actd . Since the system is overactuated, there is an infinite number of solutions unless additional constraints can be added. Here, an intuitive scheme inspired from the CGA data during double stance is used to allocate the operational force between the two legs without relying on computationally expensive optimization methods. It is observed that the leg with foot lying closest to the torso center of mass takes a greater portion of the load [3]. Thus, the following constraints are added: where J ubL and J ubR represent the Jacobian matrix at the back point in the left leg and in the right leg, respectively.
x TR and x TL represent the horizontal distance from back to right ankle and the horizontal distance from back to left ankle. With (40), P Ld can be solved.
Step 3. Specify the desired flowQ Ld for Q L so that the actual load force tracks the desired load force synthesized in Step 5.
The same as [18], the joint velocity and acceleration used to compute _ P Ld for adaptive model compensation are estimated through an adaptive robust observer.
Define the observer errors as Then the nonlinear observer can be designed as where y and _ y r are the estimated joint positions and joint velocities, K o1 is any gain matrix, M Dsp , C Dsp , and G Dsp represent the estimated matrices using new parameter estimation θ q , K o2 is the gain for linear feedback, K o2s is the gain for nonlinear feedback, and T os is the robust feedback term.
Replacing _ q c1 with _ y r in P Ld , the estimated P Ld can be obtained whereP Ld = P Ld ðq c1 , _ y r , b θ q , tÞ. Letẑ 4 = P L −P Ld . Treating Q L as the control input in this part, the proposed 13 Applied Bionics and Biomechanics control law makingẑ 4 = P L −P Ld converge to zero or bounded is synthesized as follows: where Q Lda and Q Lds are the terms for model compensation and robust feedback, respectively, ϕ 4c = ϕ 4c1 ϕ 4c2 ϕ 4c3 ϕ 4c4 ϕ 4c5 ½ T is the vector with ϕ 4c1 = 0 6×16 , ϕ 4c2 = 0 6×6 , ϕ 4c3 = 0 6×6 , ϕ 4c4 = −q v x 3 , and ϕ 4c5 = I 6×6 , K 4s1 is the linear feedback gain, K 4 > 0, g 3 > 0, d 4 > 0, and Q Ldsn is a nonlinear feedback term. 3.5. ARCFC Design for Left Leg Support. For the walking phase of left leg support, there exist six independent degree of freedoms and six control inputs, which can be seen from the dynamic model (7). Thus, we can control six independent human-machine interaction force components. Here, we select to minimize six interaction force components at the back and right foot, as shown in Figure 4.
The dynamics of left leg support is a serial chain one, and there is no overactuated characteristic in the system. Thus, the torque allocation is not needed in the controller design for left leg support. Other steps are almost the same as the ARCFC for double leg support except that the independent degrees of freedoms becomes six compared to three in double leg support. Here, for simplicity, the detail controller design process for left leg support is omitted in this paper.

Simulation Setup.
Based on the dynamic model, a simulation model is constructed in MATLAB/Simulink. The simulation parameters are the same as those in [18]. At first, the estimates of system parameters are set to be real values. The lumped disturbances are set to zero. The sampling time is selected as t s = 0:0001s. The value for desired humanmachine interaction force is set to zero. Considering the similarity between human body and lower limb exoskeleton, it is better to use human clinical gait analysis (CGA) data as desired joint position for simulation. However, the CGA data is different for people with different height and walking speed. Since the sinusoid curves are often used as desired trajectory for simulation of control algorithms. For a preliminary validation of the proposed controller performance, we choose the desired joint motion trajectory as sinusoid curves. The amplitude and the frequency of the sinusoid curves can be selected on the same order of magnitude as that of CGA data. In the simulation, the following control algorithms are conducted: C1: The Low-Level PID Control with Velocity Feedforward The control law is given by In the simulation, a Z-N method with slight adjustments is used to obtain the control gains of PID controller. For left leg support, K p = diag f10, 10, 10, 10, 10, 10g, K I = diag f10 , 10, 10, 10, 10, 10g, K d = diag f30, 30, 20, 20, 15, 15g, and V f = diag f0:1, 0:1, 0:1, 0:1, 0:1, 0:1g. For double leg support, K p = diag f10, 10, 10, 10, 10, 10g, K I = diag f10, 10, 10 , 10, 10, 10g, K d = diag f30, 30, 20, 20, 30, 30g, and V f = diag f0:1, 0:1, 0:1, 0:1, 0:1, 0:1g.
To show the control performance, three sets are simulated: Set1: motion tracking control of two low-level controllers.
Set2: nominal interaction force control. Set3: interaction force control to load change. Set4: interaction force control to human-machine interface modeling errors.
In Set4, the human-machine interface dynamics is described as a spring-damper model which means in Equation (28), the modeling errors is described asD 1 x e Þ where B hm is the damping ratio at the humanmachine interface. From Figures 11 and 12, it can be seen that a consistent performance can be achieved for the proposed ARCFC (FARC+C2) to human-machine interface modeling errors both in left leg support phase and double leg support phase. For PID cascade force controller (FARC +C1), when adding the damping in the human-machine interface, the human interaction force becomes chattering. The reason is that the closed loop bandwidth and parameter adaptation rate of PID cascade force controller are limited leading to a poor disturbance rejection performance.
In this paper, only simulations are carried out, and also, the human motion trajectory is generated by sinusoid curves for a preliminary validation of the proposed controller performance. In the future, comparative experiments will be conducted on a real lower limb hydraulic exoskeleton platform to further validate the performance of the proposed interaction force controller in practical applications.

Conclusion
In this paper, a generalized multiphase dynamic modeling method is proposed for lower limb exoskeleton in which the dynamic model of each walking phase can all be obtained based on the dynamic model of a floating lower limb exoskeleton (with positions of the exoskeleton feet changed) and different holonomic constraints, which significantly simplify the dynamic modeling process of the multiphase lower limb exoskeleton. MIMO adaptive robust interaction force controllers with high level doing human motion intent inference while low level conducting human trajectory tracking are designed both for double leg support phase and single leg support phase. A torque allocation method is proposed to deal with the overactuated characteristic in double leg support. Comparative simulations show the effectiveness and better performance of the proposed multiphase human-machine interaction force controller. In our future research, we will do the modeling and controller design of underactuated lower limb exoskeleton systems. Stability analysis of uncontrolled internal dynamics and adaptive robust force control of underactuated exoskeleton systems will be conducted.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
The authors declare that they have no conflict of interest.