Fuzzy cerebellar model articulation controller-based adaptive tracking control for load-carrying exoskeleton

Abstract

Load-carrying exoskeletons need to cope with load variations, outside disturbances, and other uncertainties. This paper proposes an adaptive trajectory tracking control scheme for the load-carrying exoskeleton. The method is mainly composed of a computed torque controller and a fuzzy cerebellar model articulation controller. The fuzzy cerebellar model articulation controller is used to approximate model inaccuracies and load variations, and the computed torque controller deals with tracking errors. Simulations of an exoskeleton in squatting movements with model parameter changes and load variations are carried out, respectively. The results show a precise tracking response and high uncertainties toleration of the proposed method.

Introduction

The exoskeleton is a kind of wearable and anthropomorphic device that enhances the physical capabilities of a non-disabled wearer.^1,2 It looks like a humanoid robot or a part of that in appearance and is worn in parallel with human body in general. The significant feature of the mechanical device is that it combines the cognitive ability of human with powerful mechanical system, ³ which makes the human–machine system more intelligent and powerful. For example, lower limb exoskeleton can adapt complex terrains well without building map in advance. ⁴ Load-carrying exoskeleton (also called load-bearing or load transfer exoskeleton) is one kind of lower limb exoskeletons. It transfers substantial loads to the ground directly rather than augments human joints power. ¹ This type of leg exoskeleton is urgently needed in the military for long-distance locomotion with heavy loads ⁵ and is widely researched during the past two decades. However, one of the most difficult problems in designing an excellent load-carrying exoskeleton is the control system. On the one hand, the device should follow the pilot’s motion flexibly. On the other hand, it should overcome load variation and other disturbances and support the remainder loads steady at the same time. In most cases, the carried loads could be a few times heavier than the device itself, and the variation of loads is a tough challenge for control system design.

Much work has been done on the control of load-carrying exoskeleton by earlier studies. But most of them deal with fixed payload or small load variation, or only focus on the uncertainty of the dynamic model. Among them, more effective and practical methods include sensitivity amplification control (SAC), model-based control (MBC), sliding mode control (SMC), and neural network-based adaptive control (NNAC). The Berkeley Lower Extremity Exoskeleton (BLEEX) is the first functional load-carrying activated exoskeleton. It uses SAC method to acquire a high sensitivity of pilot’s forces or torques without any human side information, which is hard to be measured accurately.^4,6 But the control method depends on a better dynamic model and is unable to cope with parameter variations and model uncertainties. After BLEEX, SAC is seldom being used in other exoskeletons for its strick requirements of the model. MBC, such as direct force control^7,8 and impedance control,^9,10 is relatively easy to implement. HEXAR-CR50 developed by Hanyang University uses direct force control to calculate the desired joint torques in the stance phase and then applied a proportional–integral–derivative (PID) controller at each joint.^11,12 The exoskeleton developed by Massachusetts Institute of Technology (MIT) uses similar force control method under a state-machine control strategy with the support of plenty of sensor information. ¹³ Both of them showed wearer’s metabolism reduction during locomotion.

To enhance the robustness of the exoskeleton control system, SMC ¹⁴ and intelligent control NNAC ¹⁵ are mainly investigated. SMC has an advantage of parameter-insensitive, and NNAC is good at approximating uncertainties. Hence, the method of SMC with NNAC compensation is widely studied as well.^16–18 Though much progress has been achieved in adaptive control of load-carrying exoskeleton, few works touch on the situation of a broad range of load mass and inertia variation, which is practically unavoidable. Typical cases include loading and unloading, loads waggling during moving, impacting when weapon firing, and so on.

In this paper, a computed torque control (CTC) ¹⁹ associated with fuzzy cerebellar model articulation controller (FCMAC) ²⁰ method is proposed for load-carrying exoskeleton on tracking problems. The biomimetic reference trajectory is given in advance since we focus on control techniques rather than motion intention detection, which is another specific research point. ²¹ CTC, also called inverse dynamics control, is a unique application of feedback linearization of non-linear systems ¹⁹ and performs well in trajectory tracking of robot manipulators. ²² The CTC is made up of two feedback loops: an inner loop based on the dynamic model and an outer loop operating on the tracking error. The inner loop is used to obtain a linear and decoupled input/output relationship, whereas the outer loop is required to stabilize the overall system. To satisfy the precise model requirement of CTC, we introduced neural network compensation. FCMAC is used to approximate modeling deviation, overcome parameter variations as well as represent the dynamic payload model. This neural network is a local approximation network. For each input, only a few related weights should be adjusted, so the network can approximate non-linear terms exceptionally quickly.^23,24 A clear FCMAC structure which mainly aims at significant payload variation is provided, and the unified stability analysis of the whole control scheme is proved by the Lyapunov method. Simulation results show that the proposed control method performs well during tracking even there exist model uncertainties. Moreover, the method can tolerate a broad range of payload mass over one time bigger than the mass of the device itself. With the proposed adaptive controller, exoskeletons not only help soldiers or workers transfer heavy objects but also can maintain high stability when loading and unloading without bothering operators much.

Dynamic model

As Figure 1 shows, exoskeleton model can be simplified as a four-bar linkage which represents foot, shank, thigh, and torso that are parallel with wearer’s limbs. $O_a$ , $O_k$ , and $O_h$ are the rotation centers of the ankle, knee, and hip joint, respectively. Gravity center of each linkage is located at points $G_f$ , $G_s$ , $G_t$ , and $G_{\mathrm{ub}}$ separately. Other parameters and variables are defined in Table 1.

Figure 1.

Model of the load-carrying exoskeleton in the stance phase.

Table 1.

Parameters and variables.

Notation	Meaning	Unit
$m_f$	Foot mass	kg
$m_s$	Shank mass	kg
$m_t$	Thigh mass	kg
$m_{\mathrm{ub}}$	Upper body mass	kg
$L_f$	Foot length	m
$L_s$	Shank length	m
$L_t$	Thigh length	m
$L_{\mathrm{ub}}$	Upper body length	m
$L_{G_f}$	Length from ankle to $G_f$	m
$L_{G_s}$	Length from ankle to $G_s$	m
$L_{G_t}$	Length from knee to $G_t$	m
$L_{G_{\mathrm{ub}}}$	Length from hip to $G_{\mathrm{ub}}$	m
$q_1$	Shank angle with respect to vertical	rad
$q_2$	Thigh angle with respect to shank extension	rad
$q_3$	Torso angle with respect to thigh extension	rad
${\tau }_1$	Ankle torque	Nm
${\tau }_2$	Knee torque	Nm
${\tau }_3$	Hip torque	Nm

As mentioned above, load-carrying exoskeleton saves wearer’s efforts by transferring heavy loads to the ground. The transfer progress only occurs in the stance phase. Thus, we mainly consider device activities in the stance phase of the sagittal plane. Figure 1 shows a state of two feet stand. A squat movement, which is a typical activity of the stance phase, is applied and researched in the following work.

The dynamic model of the stance phase exoskeleton is designed with Eural–Lagrange equation and is expressed as

$$\mathit{M}(\mathit{q})\stackrel{··}{\mathit{q}}+\mathit{C}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\mathit{G}(\mathit{q})=\mathit{u}+{\mathit{J}}^{\mathit{T}}(\mathit{q}){\mathit{F}}_{\mathit{int}}$$

(1)

where $\mathit{M}$ is the inertia matrix, $\mathit{C}$ is the Coriolis and centrifugal force matrix, $\mathit{G}$ is the gravity matrix, $\mathit{J}$ is the Jacobian matrix, ${\mathit{F}}_{\mathit{int}}$ is the interaction force between exoskeleton and its wearer, and $\mathit{u}=[{\tau }_1,{\tau }_2,{\tau }_3{]}^T$ is the control input torque. The interaction force ${\mathit{F}}_{\mathit{int}}$ is detected by a multidimensional force sensor, which is usually mounted on the back between human and machine. ²⁵

Note that the payload, which is usually mounted on the back of the exoskeleton torso, is not modeled considering its variation. That is, equation (1) is a non-loaded model and the existence of load will definitely influence the accuracy of $\mathit{M}$ , $\mathit{C}$ , and $\mathit{G}$ . While those inaccuracies can be pretty handled using the proposed control scheme. The friction term is not considered, as well. Assuming the friction is only related to the joint velocity, the existence of friction can be regarded as a deviation of $\mathit{C}$ . Besides that, we assume the human–robot interaction force is in linear relationship with displacement errors as expressed in equation (2)

$${\mathit{F}}_{\mathit{int}}={\mathit{K}}_{\mathit{M}}\mathrm{\Delta }\stackrel{··}{\mathit{x}}+{\mathit{K}}_{\mathit{D}}\mathrm{\Delta }\stackrel{·}{\mathit{x}}+{\mathit{K}}_{\mathit{P}}\mathrm{\Delta }\mathit{x}$$

(2)

If the exoskeleton tracks the reference trajectory decided by the wearer well, then the interaction force will be minimum, and that is what a load-carrying exoskeleton should achieve.

Method

The whole control structure is shown in Figure 2. The inputs of the exoskeleton include three parts: the physical human–robot interaction (pHRI) force ${\mathit{F}}_{\mathit{int}}$ , the control value of CTC ${\mathit{u}}_{\mathit{CTC}}$ , and the control quantity of FCMAC ${\mathit{u}}_{\mathit{CMAC}}$ . If the system model is precise enough, the CTC controller guarantees the tracking errors converging to zero by reasonable pole placement. The precise model is ensured by the FCMAC controller but cannot be ideal. Thus, a compensated controller ${\mathit{u}}_{\mathit{c}}$ is introduced.

Figure 2.

Block diagram of the proposed control scheme.

CTC controller

The control law is designed as the following form

$$\mathit{u}=\hat{\mathit{M}}(\mathit{q})\mathbf{\tau }+\hat{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\hat{\mathit{G}}(\mathit{q})-{\mathit{J}}^T(\mathit{q}){\mathit{F}}_{\mathit{int}}$$

(3)

$$\mathbf{\tau }={\stackrel{··}{\mathit{q}}}_{\mathit{d}}+{\mathit{K}}^T\mathit{E}$$

(4)

$$\mathit{E}=\left[\begin{array}{c}\mathrm{\Delta }\mathit{q}\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }\stackrel{·}{\mathit{q}}}\end{array}\right],\mathit{K}=\left[\begin{array}{c}{\mathit{K}}_1\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_2}\end{array}\right]$$

(5)

where $\mathrm{\Delta }\mathit{q}={\mathit{q}}_{\mathit{d}}-\mathit{q}$ , $\mathrm{\Delta }\stackrel{·}{\mathit{q}}={\stackrel{·}{\mathit{q}}}_{\mathit{d}}-\stackrel{·}{\mathit{q}}$ , ${\mathit{K}}_1,{\mathit{K}}_2\in {\mathbb{R}}^{3\times 3}$ .

Substituting equations (3) and (4) into equation (1), the system dynamic equation is updated as

$$\hat{\mathit{M}}(\mathit{q})\mathrm{\Delta }\stackrel{··}{\mathit{q}}+\hat{\mathit{M}}(\mathit{q}){\mathit{K}}^T\mathit{E}=\tilde{\mathit{M}}(\mathit{q})\stackrel{··}{\mathit{q}}+\tilde{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\tilde{\mathit{G}}(\mathit{q})$$

(6)

where $\tilde{\mathit{M}}(\mathit{q})=\mathit{M}(\mathit{q})-\hat{\mathit{M}}(\mathit{q})$ , $\tilde{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})=\mathit{C}(\mathit{q},\stackrel{·}{\mathit{q}})-\hat{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})$ and $\tilde{\mathit{G}}(\mathit{q})=\mathit{G}(\mathit{q})-\hat{\mathit{G}}(\mathit{q})$ . If

$$\tilde{\mathit{M}}(\mathit{q})\stackrel{··}{\mathit{q}}+\tilde{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\tilde{\mathit{G}}(\mathit{q})=0$$

(7)

Take equations (5) and (7) into equation (6), we get

$$\mathrm{\Delta }\stackrel{··}{\mathit{q}}+{\mathit{K}}_2\mathrm{\Delta }\stackrel{·}{\mathit{q}}+{\mathit{K}}_1\mathrm{\Delta }\mathit{q}=0$$

(8)

The system will be stable as long as the coefficient matrices ${\mathit{K}}_1$ and ${\mathit{K}}_2$ are positive definite, and the position error $\mathrm{\Delta }\mathit{q}$ will converge to zero. That is, the exoskeleton can follow wearer’s movement pretty well without hindering each other. However, it’s quite challenging to make equation (7) tenable because of the exists and the unpredictability of the dynamic uncertainties and external disturbances. Thus, the ideal control law equations (3) and (4) are practically limited used. To counteract those uncertainties and guarantee the tracking error $\mathrm{\Delta }\mathit{q}$ converging to zero, a compensation control quantity is added based on FCMAC network.

Architecture of FCMAC

FCMAC architecture is similar to the conventional CMAC. The difference is that FCMAC uses Gaussian basis function as receptive field basis function rather than local constant binary function. Figure 3 shows the architecture of an $n$ -dimensional input and $p$ -dimensional output CMAC. The network is constructed with five spaces: input space $\mathit{S}$ , association space $\mathit{A}$ , receptive field space $\mathit{T}$ , adjustable weight space $\mathit{W}$ , and output space $\mathit{Y}$ . And the network mainly has three mapping processes

$$\mathit{R}:\mathit{S}\Rightarrow \mathit{A},\mathit{U}:\mathit{A}\Rightarrow \mathit{T},\mathit{P}:\mathit{T}\Rightarrow \mathit{Y}$$

(9)

The function $\phi =\mathit{R}(\mathit{s})$ maps each input vector $\mathit{s}$ onto the association space $\mathit{A}$ . According to $\phi$ , function $\psi =\mathit{U}(\phi )$ can select specific hypercubes, which constitute the receptive field space $\mathit{T}$ . Each hypercube corresponds with a function that is obtained by multiplying related basis functions. At last, function $\mathit{y}=\mathit{P}(\psi )$ computes output vector $\mathit{y}$ by summing the weighted value of all the selected receptive field functions. The details are introduced as follows.

Figure 3.

Architecture of CMAC.

Quantization

Each data of the input vector s should be quantized into a given interval to map it onto the association space conveniently. The units and ranges of the inputs may be different, but the quantization converts them into uniform dimensionless values. The minimum scale of the quantization zone decides the side length of the hypercube, which represents the size of each receptive field. Usually, we define the minimum scale as 1, and the number of scales $n_s$ is called resolution. The quantization can be linear or non-linear, and the difference has not been systematically researched. In this study, we choose $\mathit{E}\in {\mathbb{R}}^{6\times 1}$ as input that can be measured relative easily in the physical system. Resolutions of all the $n=6$ inputs are defined the same as 9 with the minimum scale as 1. Figure 4 illustrates the quantized coordination of the inputs $s_1$ and $s_2$ .

Figure 4.

Quantization and mapping of CMAC with two inputs.

Receptive field function

After quantization, each input is mapped onto the association space and be represented by multilayer blocks. The association space of each input is constructed with $n_l$ layers. Each layer is divided into $n_b$ blocks. The maximum block length is defined as $l_b$ . The blocks in each layer are mismatched. In general, $n_s$ , $n_l$ , $n_b$ , and $l_b$ have a relation as shown in equation (10), so that the number of blocks in each layer is the same and the mismatch interval exactly is a minimum scale

$$n_l=l_b,n_b=\frac{(n_s-1)}{l_b+1}$$

(10)

Thus, any input data can be mapped to $n_l$ blocks and each minimum scale corresponding to $n_l$ specific blocks. For each input, there are totally $n_B=n_l{n}_b$ blocks.

Each block corresponds with a receptive field basis function. The basis function is effected only in its related block. Any continuously bounded function can be used as the basis function, and the Gaussian function is used in this work. The basis function of the kth block of the jth layer of the ith input $s_i$ is presented as

$${\varphi }_{\mathrm{ijk}}(s_i)=\exp \left[\frac{-{(s_i-m_{\mathrm{ijk}})}^2}{{\sigma }_{\mathrm{ijk}}^2}\right]$$

(11)

Note that $k$ is related with $i$ and $j$ .

Any input point s can be represented by a set of $n_l$ hypercubes in the receptive field. Each hypercube has a receptive field function that is composed of basic functions of the related $n$ blocks. Note that the related blocks are located in the same layer of different inputs. The jth layer hypercube function of the input s is defined as

$$b_{H_j}(\mathit{s},\mathit{m},\mathit{\sigma })=\underset{i=1}{\overset{n}{\Pi }}{\varphi }_{\mathrm{ijk}}=\exp \left[\sum _{i=1}^n\frac{-{(s_i-m_{\mathrm{ijk}})}^2}{{\sigma }_{\mathrm{ijk}}^2}\right]$$

(12)

where $\mathit{m},\mathit{\sigma }\in {\mathbb{R}}^{n{n}_B\times 1}$

$$\mathit{m}=[\underset{\mathrm{the}\mathrm{first}\mathrm{block}}{\underbrace{m_{11}\cdots m_{1n}}}\underset{\mathrm{the}\mathrm{second}\mathrm{block}}{\underbrace{m_{21}\cdots m_{2n}}}\cdots \underset{\mathrm{the}{n}_B\mathrm{th}\mathrm{block}}{\underbrace{m_{n_{B}1}\cdots m_{n_{B}n}}}{]}^T$$

(13)

$$\mathit{\sigma }=[\underset{\mathrm{the}\mathrm{first}\mathrm{block}}{\underbrace{{\sigma }_{11}\cdots {\sigma }_{1n}}}\underset{\mathrm{the}\mathrm{second}\mathrm{block}}{\underbrace{{\sigma }_{21}\cdots {\sigma }_{2n}}}\cdots \underset{\mathrm{the}{n}_B\mathrm{th}\mathrm{block}}{\underbrace{{\sigma }_{n_{B}1}\cdots {\sigma }_{n_{B}n}}}{]}^T$$

(14)

Define $\mathbf{\Gamma }\in {\mathbb{R}}^{n_b^n{n}_l\times 1}$ as the vector of hypercube functions

$$\mathbf{\Gamma }(\mathit{s},\mathit{m},\mathit{\sigma })=[\underset{\mathrm{the}\mathrm{first}\mathrm{layer}}{\underbrace{b_{11}\cdots b_{1n_b^{n_l}}}}\underset{\mathrm{the}\mathrm{second}\mathrm{layer}}{\underbrace{b_{21}\cdots b_{2n_b^{n_l}}}}\cdots \underset{\mathrm{the}{n}_l\mathrm{th}\mathrm{layer}}{\underbrace{b_{n_{l}1}\cdots b_{n_l{n}_b^{n_l}}}}{]}^T$$

(15)

There are $n_h=n_l\times n_b^n$ hypercubes in total. Different hypercubes can be overlapped with some same receptive field blocks. We called the overlap as local generalization, which can be understood as neighboring inputs producing similar outputs. If there is no overlap between the two inputs, there is no generalization between them and they will not influence each other when the network is training.

All networks have the problem of dimensional explosion, and CMAC is no exception. As mentioned above, a CMAC with $n$ inputs, if divided into $n_l$ layers and each layer into $n_b$ blocks, will have $n_l\times n_b^n$ hypercubes. The amount of computation could be horrible with the increase in $n$ and $n_b$ . However, using hash coding, the progress of mapping from association space to receptive field space can be compressed and does not influence the local generalization much at the same time.^26,27 Thus, the computation amount and the space memory of receptive field can be reduced to an acceptable range, $n_h=n_B$ , for example. ²⁸

Weight space and output

Each hypercube corresponds to a specific weight when mapping to a specific output. Suppose the output dimension is $p$ , the weight space W is $p$ times large as the receptive field space T

$$\mathit{w}=[{\mathit{w}}_1{\mathit{w}}_2\cdots {\mathit{w}}_p]=\left[\begin{array}{cccc}{w}_{11}\hfill & w_{21}\hfill & \cdots \hfill & w_{p1}\hfill \\ \multicolumn{1}{c}{w_{21}}& w_{22}\hfill & \cdots \hfill & w_{p2}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{w_{n_{h}1}}& w_{n_{h}2}\hfill & \cdots \hfill & w_{n_{h}p}\hfill \end{array}\right]$$

(16)

in which ${\mathit{w}}_{\mathbf{i}}(i\in 1,2,\cdots ,p)$ is the weight vector of the ith output $y_i$ . The output $\mathit{y}\in {\mathbb{R}}^{n_h\times p}$ is expressed as

$$\mathit{y}={[y_1,y_2,\cdots ,y_p]}^T={\mathit{w}}^T\mathbf{\Gamma }(\mathit{s},\mathit{m},\mathit{\sigma })$$

(17)

For example, if two input vectors are quantized to (2.1, 2.4) and (2.5, 2.6), respectively, they are both located at the square (3,3) and represented by the hypercubes ¹C²C, ¹D²D, ¹E²E, and ¹F²F.

FCMAC and compensated controller

Redesign the control law equation (3) as

$$\mathit{u}=^M\mathbf{\tau }+\hat{\mathit{C}}\stackrel{·}{\mathit{q}}+\hat{\mathit{G}}-{\mathit{J}}^T(\mathit{q}){\mathit{F}}_{\mathit{int}}+{\mathit{u}}_{\mathit{CMAC}}$$

(18)

Then equation (6) can be rewritten as

$$\mathrm{\Delta }\stackrel{··}{\mathit{q}}=-{\mathit{K}}^T\mathit{E}+\hat{\mathit{M}}(\mathit{q}{)}^{-1}\mathbf{\Psi }$$

(19)

in which $\mathbf{\Psi }$ represents the uncertainties and the neural network output term

$$\mathbf{\Psi }=\tilde{\mathit{M}}(\mathit{q})\stackrel{··}{\mathit{q}}+\tilde{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\tilde{\mathit{G}}(\mathit{q})-{\mathit{u}}_{\mathit{CMAC}}$$

(20)

The system state-space dynamic model can be written as

$$\stackrel{·}{\mathit{E}}=\mathit{AE}+\mathit{B}\mathbf{\Psi }$$

(21)

where

$$\mathit{E}=\left[\begin{array}{c}\mathrm{\Delta }\stackrel{·}{\mathit{q}}\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }\stackrel{··}{\mathit{q}}}\end{array}\right],\mathit{A}=\left[\begin{array}{cc}0\hfill & {\mathit{I}}^{3\times 3}\hfill \\ \multicolumn{1}{c}{-{\mathit{K}}_1}& -{\mathit{K}}_2\hfill \end{array}\right],\mathit{B}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{{\mathit{I}}^{3\times 3}}\end{array}\right].$$

If there exists ideal weight values ${\mathit{w}}^{*}$ and parameters ${\mathit{m}}^{*},{\mathit{\sigma }}^{*}$ that can make $\mathbf{\Psi }\to 0$ , then the uncertainties can be represented as

$$\begin{array}{c}{\mathit{u}}_{\mathit{CMAC}}^{*}=\tilde{\mathit{M}}(\mathit{q})\stackrel{··}{\mathit{q}}+\tilde{\mathit{C}}(\mathit{q},\stackrel{·}{\mathit{q}})\stackrel{·}{\mathit{q}}+\tilde{\mathit{G}}(\mathit{q})\hfill \\ ={\mathit{w}}^{*T}\Gamma (\mathit{s},{\mathit{m}}^{*},{\mathit{\sigma }}^{*})+\epsilon (\mathit{s})\hfill \end{array}$$

(22)

in which $\epsilon (\mathit{s})$ is the functional reconstructional error and $\parallel \epsilon (\mathit{s})\parallel \le {\epsilon }_M$ is bounded. ²⁹ Let $\hat{\mathit{w}},\hat{\mathit{m}}$ , and $\hat{\sigma }$ are the estimated values of ${\mathit{w}}^{*}$ , ${\mathit{m}}^{*}$ , and ${\mathit{\sigma }}^{*}$ , the output of FCMAC network ${\mathit{u}}_{\mathit{f}}$ is

$${\mathit{u}}_{\mathit{f}}={\hat{\mathit{w}}}^T\hat{\mathbf{\Gamma }}(\mathit{s},\hat{\mathit{m}},\hat{\sigma })$$

(23)

To minimize the approximation error of the network, a compensation quantity ${\mathit{u}}_{\mathit{c}}$ is added. The approximation error $\tilde{\mathit{u}}$ is expressed as

$$\tilde{\mathit{u}}={\mathit{u}}^{*}-{\mathit{u}}_{\mathit{CMAC}}={\tilde{\mathit{w}}}^T{\mathbf{\Gamma }}^{*}+{\hat{\mathit{w}}}^T\tilde{\Gamma }+\epsilon -{\mathit{u}}_{\mathit{c}}$$

(24)

where $\tilde{\mathit{w}}={\mathit{w}}^{*}-\hat{\mathit{w}}$ and $\tilde{\mathbf{\Gamma }}={\mathbf{\Gamma }}^{*}-\hat{\mathbf{\Gamma }}$ . Linearize $\tilde{\mathbf{\Gamma }}$ using Taylor expansion

$$\begin{array}{c}\tilde{\mathbf{\Gamma }}={{[\frac{\partial b_1}{\partial \mathit{m}},\frac{\partial b_2}{\partial \mathit{m}},\cdots ,\frac{\partial b_{n_h}}{\partial \mathit{m}}]}^T|}_{\mathit{m}=\hat{\mathit{m}}}({\mathit{m}}^{*}-\hat{\mathit{m}})\hfill \\ +{{[\frac{\partial b_1}{\partial \mathit{\sigma }},\frac{\partial b_2}{\partial \mathit{\sigma }},\cdots ,\frac{\partial b_{n_h}}{\partial \mathit{\sigma }}]}^T|}_{\mathit{\sigma }=\hat{\mathit{\sigma }}}({\mathit{\sigma }}^{*}-\hat{\mathit{\sigma }})+{\mathit{O}}_{\mathit{t}}\hfill \\ \equiv {\mathit{G}}^{\mathit{T}}\tilde{\mathit{m}}+{\mathit{H}}^{\mathit{T}}\tilde{\mathit{\sigma }}+{\mathit{O}}_{\mathit{t}}\hfill \end{array}$$

(25)

where ${\mathit{O}}_t\in {\mathbb{R}}^{n_h\times 1}$ is a higher-order term and $\mathit{G},\mathit{H}\in {\mathbb{R}}^{n{n}_B\times n_h}$ are the coefficient matrices with elements expressed as (symbol “ $\stackrel{\mathrm{c}}{\Rightarrow }$ ” means the block left of the symbol is one of the elements that construct the hypercube right of the symbol)

$$\frac{\partial b_k}{\partial m_{\mathrm{ij}}}=\{\begin{array}{c}\frac{2(s_j-m_{\mathrm{ij}})}{{\sigma }_{\mathrm{ij}}^2}{b}_k(\mathit{s},\mathit{m},\mathit{\sigma }),\mathrm{if}{\varphi }_{\mathrm{ij}}\stackrel{\mathrm{c}}{\Rightarrow }{b}_k\hfill \\ \multicolumn{1}{c}{0,\mathrm{else}}\end{array}$$

(26)

$$\frac{\partial b_k}{\partial {\sigma }_{\mathrm{ij}}}=\{\begin{array}{c}\frac{2(s_j-m_{\mathrm{ij}}{)}^2}{{\sigma }_{\mathrm{ij}}^3}{b}_k(\mathit{s},\mathit{m},\mathit{\sigma }),\mathrm{if}{\varphi }_{\mathrm{ij}}\stackrel{\mathrm{c}}{\Rightarrow }{b}_k\hfill \\ \multicolumn{1}{c}{0,\mathrm{else}}\end{array}$$

(27)

The optimal hypercube function set is represented as

$${\mathbf{\Gamma }}^{*}=\hat{\mathbf{\Gamma }}+\tilde{\mathbf{\Gamma }}=\hat{\mathbf{\Gamma }}+{\mathit{G}}^{\mathit{T}}\tilde{\mathit{m}}+{\mathit{H}}^{\mathit{T}}\tilde{\mathit{\sigma }}+{\mathit{O}}_{\mathit{t}}$$

(28)

Substitute equations (20), (22), (24), and (28) to equation (21)

$$\begin{array}{c}\stackrel{·}{\mathit{E}}=\mathit{AE}+\mathit{B}\tilde{u}\hfill \\ =\mathit{AE}+\mathit{B}({\tilde{\mathit{w}}}^T{\mathbf{\Gamma }}^{*}+{\hat{\mathit{w}}}^T\tilde{\mathbf{\Gamma }}+\epsilon -{\mathit{u}}_{\mathit{c}})\hfill \\ =\mathit{AE}+\mathit{B}({\tilde{\mathit{w}}}^T\hat{\mathbf{\Gamma }}+{\tilde{\mathit{w}}}^T\tilde{\mathbf{\Gamma }}+{\hat{\mathit{w}}}^T\tilde{\mathbf{\Gamma }}+\epsilon -{\mathit{u}}_{\mathit{c}})\hfill \\ =\mathit{AE}+\mathit{B}[{\tilde{\mathit{w}}}^T\hat{\mathbf{\Gamma }}+{\hat{\mathit{w}}}^T({\mathit{G}}^T\tilde{\mathit{m}}+{\mathit{H}}^T\tilde{\mathit{\sigma }})+\mathit{\eta }-{\mathit{u}}_{\mathit{c}}]\hfill \end{array}$$

(29)

in which $\mathit{\eta }$ is the error of the network and is bounded ${\mathit{\eta }}_0$ . The FCMAC and the compensated control laws are directly given using the Lyapunov method

$$\stackrel{·}{\hat{\mathit{w}}}={\gamma }_1\hat{\mathbf{\Gamma }}{\mathit{E}}^T\mathit{PB}$$

(30)

$$\stackrel{·}{\hat{\mathit{m}}}={\gamma }_2\mathit{G}\hat{w}{\mathit{B}}^T\mathit{PE}$$

(31)

$$\stackrel{·}{\hat{\mathit{\sigma }}}={\gamma }_3\mathit{H}\hat{w}{\mathit{B}}^T\mathit{PE}$$

(32)

$${\mathit{u}}_{\mathit{c}}=\hat{\mathit{\eta }},\stackrel{·}{\hat{\mathit{\eta }}}={\gamma }_4{\mathit{E}}^T\mathit{PB}$$

(33)

where ${\gamma }_1,{\gamma }_2,{\gamma }_3$ , and ${\gamma }_4$ are some constants greater than zero. $\mathit{P}$ is a symmetric positive definite matrix and fulfills the following Lyapunov equation

$$\mathit{PA}+{\mathit{A}}^T\mathit{P}=-\mathit{Q}$$

(34)

and $\mathit{Q}$ is a positive definite matrix.

As for the load-carrying exoskeleton described in equation (1), the control laws, which correspond with the control scheme in Figure 2, are entirely designed as

$$\mathit{u}={\mathit{u}}_{\mathit{CTC}}+{\mathit{u}}_{\mathit{CMAC}}={\mathit{u}}_{\mathit{CTC}}+{\mathit{u}}_{\mathit{f}}+{\mathit{u}}_{\mathit{c}}$$

(35)

in which ${\mathit{u}}_{\mathit{CTC}}$ is designed as equations (3)–(5), ${\mathit{u}}_{\mathit{f}}$ is expressed in equation (23) with parameters presented in equations (30)–(32), and ${\mathit{u}}_{\mathit{c}}$ is defined in equation (33).

Stability analysis

Define the Lyapunov function as

$$\begin{array}{c}\multicolumn{1}{c}{V(\mathit{E},\tilde{\mathit{w}},\tilde{\mathit{m}},\tilde{\mathit{\sigma }},t)=\frac{1}{2}{\mathit{E}}^T\mathit{PE}+\frac{1}{2{\gamma }_1}\mathrm{tr}\parallel \tilde{\mathit{w}}{\parallel }_F^2+\frac{1}{2{\gamma }_2}{\tilde{\mathit{m}}}^T\tilde{\mathit{m}}}\\ +\frac{1}{2{\gamma }_3}{\tilde{\mathit{\sigma }}}^T\tilde{\mathit{\sigma }}+\frac{1}{2{\gamma }_4}{\tilde{\mathit{\eta }}}^T\tilde{\mathit{\eta }}\hfill \end{array}$$

(36)

Taking the derivative of the Lyapunov function and substituting equations (29) and (34) successively, we get

$$\begin{array}{c}\stackrel{·}{V}=\frac{1}{2}{\stackrel{·}{\mathit{E}}}^T\mathit{PE}+\frac{1}{2}{\mathit{E}}^T\mathit{P}\stackrel{·}{\mathit{E}}+\frac{1}{{\gamma }_1}\mathrm{tr}({\stackrel{·}{\stackrel{~}{\mathit{w}}}}^T\tilde{\mathit{w}})-\frac{1}{{\gamma }_2}{\stackrel{·}{\hat{\mathit{m}}}}^T\tilde{\mathit{m}}\hfill \\ -\frac{1}{{\gamma }_3}{\stackrel{·}{\hat{\mathit{\sigma }}}}^T\tilde{\mathit{\sigma }}-\frac{1}{{\gamma }_4}{\stackrel{·}{\hat{\mathit{\eta }}}}^T\tilde{\mathit{\eta }}\hfill \\ =\frac{1}{2}{\mathit{E}}^T({\mathit{A}}^T\mathit{P}+\mathit{PA})\mathit{E}+{\mathit{E}}^T\mathit{PB}\mathbf{\Psi }+\frac{1}{{\gamma }_1}\mathrm{tr}({\stackrel{·}{\stackrel{~}{\mathit{w}}}}^T\tilde{\mathit{w}})\hfill \\ -\frac{1}{{\gamma }_2}{\stackrel{·}{\hat{\mathit{m}}}}^T\tilde{\mathit{m}}-\frac{1}{{\gamma }_3}{\stackrel{·}{\hat{\mathit{\sigma }}}}^T\tilde{\mathit{\sigma }}-\frac{1}{{\gamma }_4}{\stackrel{·}{\hat{\mathit{\eta }}}}^T\tilde{\mathit{\eta }}\hfill \\ =-\frac{1}{2}{\mathit{E}}^T\mathit{QE}+\frac{1}{{\gamma }_1}[\mathrm{tr}({\gamma }_1{\mathit{E}}^T\mathit{PB}{\tilde{\mathit{w}}}^T\hat{\mathbf{\Gamma }})+\mathrm{tr}({\stackrel{·}{\stackrel{~}{\mathit{w}}}}^T\tilde{\mathit{w}})]\hfill \\ +\frac{1}{{\gamma }_2}[{\gamma }_2{\mathit{E}}^T\mathit{PB}{\hat{\mathit{w}}}^T{\mathit{G}}^T\tilde{\mathit{m}}-{\stackrel{·}{\hat{\mathit{m}}}}^T\tilde{\mathit{m}}]\hfill \\ +\frac{1}{{\gamma }_3}[{\gamma }_3{\mathit{E}}^T\mathit{PB}{\hat{\mathit{w}}}^T{\mathit{H}}^T\tilde{\mathit{\sigma }}-{\stackrel{·}{\hat{\mathit{\sigma }}}}^T\tilde{\mathit{\sigma }}]\hfill \\ +\frac{1}{{\gamma }_4}[{\gamma }_4{\mathit{E}}^T\mathit{PB}(\mathit{\eta }-{\mathit{u}}_{\mathit{c}})-{\stackrel{·}{\hat{\mathit{\eta }}}}^T\tilde{\mathit{\eta }}]\hfill \end{array}$$

(37)

According to equations (30)–(33), equation (37) is simplified as

$$\stackrel{·}{V}=-\frac{1}{2}{\mathit{E}}^T\mathit{QE}+0+0+0+0=-\frac{1}{2}{\mathit{E}}^T\mathit{QE}\le 0$$

(38)

where $\stackrel{·}{V}$ is zero only occurs at the situation of $\mathit{E}=0$ . According to Barbalat’s Lemma, $\mathit{E}\to 0$ as $t\to \infty$ .

Simulation

The proposed method is verified through simulations. The leg on the stance phase rather than the swing phase is selected as the application object because the load variation has a slight influence on the swing leg. A squat movement is used in the simulation. The desired trajectories of the ankle, knee, and hip joints are obtained through the fitting results of human squatting data, which are collected by the motion analysis equipment. Figure 5 is the squatting progress during the simulation.

Figure 5.

Exoskeleton squat progress simulation. The payload is fixed on the back.

The average mass $m$  = 25 kg of the existing exoskeleton products ³⁰ is chosen as the total weight of the load-carrying exoskeleton simulation model. The parameter values in Table 1 are specified as $m_f$  = 3.0 kg, $m_s$  = 6.0 kg, $m_t$  = 6.0 kg, $m_{u}b$  = 10.0 kg, $L_f$  = 270 mm, $L_s$  = 520 mm, $L_t$  = 450 mm, and $L_{\mathrm{ub}}$  = 530 mm. Assume all the links are homogeneous, and the mass point is located at the geometrical center.

As for FCMAC network, joint angle displacements and their differential $\mathit{E}$ are chosen as the inputs and the outputs ${\mathit{u}}_{\mathit{CMAC}}$ are the joint torque control quantity. Assume the three joint angle displacements $\mathrm{\Delta }\mathit{q}$ have a range of $[-\pi /90,\pi /90]$ , and the differential of angular deviations $\mathrm{\Delta }\stackrel{·}{\mathit{q}}$ changes within the range of $[-\pi /15,\pi /15]$ . All of the six inputs are quantized to the section $[0,9]$ with a minimum scale of 1. The block length $b_l$ is chosen as 4. Hence, there are four layers for each input in the association space, and each layer has three blocks. The mean values of the basis function $\mathit{m}$ are initialized as the center value of each block and the standard deviations $\mathit{\sigma }$ as zero. The weight $\mathit{w}$ is initialized as zero as well.

Three groups of simulation in different conditions are designed in this paper. The first is without model uncertainties. The second is with model parameter changes, and the last is with load changes.

Simulation without uncertainties

In this condition, no load is added on the exoskeleton. The dynamic system model is definite and is used in the CTC controller design. The parameters used in the controller are designed as $K_1=600{\mathit{I}}^{3\times 3},K_2=100{\mathit{I}}^{3\times 3},Q=100{\mathit{I}}^{6\times 6},{\gamma }_1=2000,{\gamma }_2=100,{\gamma }_3=100$ , and ${\gamma }_4=10$ . Figure 6 is the trajectory tracking results during a period of squatting. The system tracks the desired trajectory well whether using the neural network or not, which indicates that the CTC controller is adequate. Moreover, from the tracking errors, we can see the CMAC controller can restrain errors and plays a positive role during tracking.

Figure 6.

Squatting control without uncertainties. (a) Trajectory tracking. (b) Tracking errors.

Simulation with model parameter changes

The model parameter changes come from the modeling simplification and the measuring errors. In this condition, the load is not considered for the moment. For simplicity, we set the model parameter changing by adjusting the value of $\mathit{M}$ , $\mathit{C}$ , and $\mathit{G}$ in equation (1). The maximum variable quantity is $\pm 40\%$ . The controllers are designed identical as in the above subsection. The tracking results are displayed in Figures 7 and 8.

Figure 7.

Squatting control with model parameters increase 40%. (a) trajectory tracking and (b) trajectory errors.

Figure 8.

Squatting control with model parameters decrease 40%. (a) Trajectory tracking. (b) Trajectory errors.

In the condition that the model parameters increase 40%, the maximum tracking errors of the control without FCMAC and with FCMAC are $2.12\times {10}^{-2}$ and $1.8\times {10}^{-3}$  rad, respectively. Both of the maximum errors occur at the knee joint. Favorable tracking response is obtained when using the FCMAC controller. When the model parameters decrease 40%, a similar control effect appears that the maximum tracking errors of control without FCMAC and with FCMAC are $-3.84\times {10}^{-2}$ and $-1.6\times {10}^{-3}$  rad, respectively. From the above two simulations, it can be seen that without the FCMAC controller, though the joints can track the desired trajectory, the tracking errors are relatively high, which increases the interaction force between the pilot and exoskeleton. However, with the FCMAC controller, more precise tracking responses can be guaranteed.

Simulation with load changes

The load adheres to the exoskeleton torso. During the working, the load supported by the exoskeleton can be changed frequently; for example, loading and unloading equipment, weapon launching, fuel consumption, and so on. Thus the control system should have the ability to cope with load variations. In this condition, only load mass variation is considered, and the load position and orientation variations are not taken into account, which can be regarded as the model parameter changes in the above subsection.

The initial state of the system used in the simulation is non-loaded. Three groups of simulation are designed with 10, 20, and 30 kg load, respectively. The controllers stay the same as subsection “Simulation without Uncertainties.” The tracking responses are illustrated in Figure 9. For comparison, the tracking responses without FCMAC controller are displayed in Figure 10.

Figure 9.

Squatting control of load changes with FCMAC. (a) Trajectory tracking. (b) Trajectory errors.

Figure 10.

Squatting control of load changes without FCMAC. (a) Trajectory tracking. (b) Trajectory errors.

Along with the adding of load, the tracking error of both the controls with FCMAC and without FCMAC is increasing. However, the control without FCMAC is hard to track the desired trajectory, especially when the load is heavy. The maximum tracking errors of the three joints under the condition of 30 kg load in Figure 10 are 0.05, 0.21 and −0.33 rad, respectively. While using the control with FCMAC, the maximum tracking errors are reduced to $-1.3\times {10}^{-3}$ , $2.6\times {10}^{-3}$ and $1.7\times {10}^{-3}$  rad, respectively. The proposed control scheme can well deal with the load variation.

Conclusion

A trajectory tracking controller associated with FCMAC compensation for load-carrying exoskeleton in the stance phase is proposed in this paper. The control scheme is mainly composed of CTC and FCMAC. The CTC is the primary controller for the trajectory tracking and has a benefit of smooth control quantity compared with the sliding mode controller. The FCMAC is used to approximate the dynamic model uncertainties. The system stability is proved by the Lyapunov function. The simulation results show that the proposed control strategy has a perfect tracking response in squatting movements and can better deal with situations, such as model uncertainties and load variations.

The proposed controller has a robust ability of system uncertainties. As to the flexibility of the load-carrying exoskeleton, it can be affected by many factors besides the robust control. First, the flexibility is related with the load weight and the motion state. In practical applications, the power of the actuators of the exoskeleton is limited. Second, the load-carrying exoskeleton can be regarded as a biped robot. During walking, the flexibility of the exoskeleton is mainly affected by the swing phase and the swing–stance–switch phase. Then, if the system uncertainties are too large, the proposed control method cannot guarantee the tracking effect. Because the higher-order term $O_t$ in equation (25) will have large impact on the controller. Also, in this paper, the desired trajectory is given directly. Actually, the acquisition of the desired trajectory is a challenge of intention perception of human motion. The accuracy of motion intention perception decides the control effect of the exoskeleton. So, if the power of the actuators is large enough, the system uncertainties are limited (the load is less than 30 kg and the errors of the model parameters are below 40%) and the desired trajectory is accuracy enough, the proposed FCMAC-based adaptive tracking control impacts the flexibility of the exoskeleton less.

Future work includes designing proper association space of FCMAC, researching non-linear quantization of FCMAC, and simulating the control strategy on various movements. Considering that the proposed control strategy has a high tolerance of uncertainties, the swing leg of the exoskeleton can be regarded as disturbances of the stance phase model, thus the control strategy of the complex activities can be simplified. The method will be evaluated on a physical system in the future.