1 Introduction

The multiple-flexible-link manipulator (MFLM) has become the key equipment of the space station and is widely used to dock spacecraft and transfer cargo. Meanwhile, due to the potential advantages of lower cost, higher operational speed, and safer operation, it has also received increasing attentions in the equipment manufacturing industry.

However, lightweight construction causes link flexibility, which results in bending and vibration of the flexible links during, at the end of, and after the maneuver. Then, the rigid motion and flexible vibration are coupled in the MFLM. Meanwhile, due to the interaction of these motions, the dynamics of the MFLM is highly coupled, nonlinear, and complex, which, in turn, leads the modeling and control of the MFLM to be a challenging task.

The modeling work is aimed at providing precise model to describe the actual MFLM's behavior. Several approaches, such as finite element technique, assumed mode method, and lumped parameter models, have been proposed^[1]. However, for the controller design application, the model is also required to be simple enough. The most common one is the adoption of the Euler-Lagrange formulation and assumed mode method, which yields closed-form expressions^[1-4].

However, even for the simplified model, the controller design is still a difficult task since the controller has to follow a designed joint trajectory and cancel arms' vibrations simultaneously. In the literature, two-loop structure control schemes were investigated. The inner loop tracks the designed trajectory and guarantees robustness of the system.While the outer loop damps links' the vibration and deal with the unmodeled dynamics.The control scheme with proportional derivative (PD) inner-loop controller and proportional integral derivative (PID) outer-loop controller was proposed in Ref.[5]. The author further proposed an adaptive algorithm to tune PD and PID gains^[6]. In Ref.[7], the outer-loop controller was designed based on the integral resonant control technique. Generalized proportional integral (GPI) controllers were proposed for both loops in Ref.[8], with disturbance compensated in the outer-loop based on an observer.The frequency identification and/or modal filter are also widely applied to the flexible-link manipulator. In Ref.[9], the author presented a control scheme, which combines the online natural frequency identification and the GPI output feedback controller, for a flexible robotic arm with unknown tip mass. A time-delay compensator combined with a modal filter was proposed in Ref.[10]. Other approaches, such as genetic learning approach^[11], distributed control strategy^[12], fuzzy logic control and H_∞ control^[13] were also investigated to control the flexible-link manipulator.

The MFLM is an under-actuated system since the single actuating torque has to control the joint's angular position and each link's vibration.Considering this feature and noting that the links' vibrations are relatively faster than the joints' rotation, singular perturbation technique is a candidate solution to the control problem of such under-actuated systems. Singular perturbation approach was employed to control the manipulator with flexible-joint^[14-16], flexible-link^[17-18], and both flexible-link and flexible-joint^[4]. It was also applied to control of the flexible manipulator actuated by pneumatic actuators^[19] and cooperate trajectory tracking of two robots^[20].

Active disturbance rejection control (ADRC), originally presented in Ref.[21], combines the advantages of both model-based and model-free methods and provides significant contribution to the model uncertain^[22-23]. The ADRC extends the system model to represent those internal and external disturbances (total disturbance). These extended states are estimated online and used to compensate the perturbation. However, the concept of ADRC was only used for rigid manipulators^[22] and flexible-joint manipulators^[23].

This paper proposes an integrated control scheme consisting of ADRC and linear quadratic regulator (LQR) technique is proposed to control MFLMs. The system of MFLM is separated into two subsystems via the singular perturbed technique. Then ADRC is applied to the slow subsystem. The proposed ADRC structure employs a tracking differentiator (TD), an extended state observer (ESO), and a feedback controller with feedforward term. The TD filter is used to preshape the desired trajectory, while the ESO is utilized to estimate the internal state as well as the total disturbance. Two feedback controllers were designed, namely PD and nonlinear PD (NPD) feedback controllers. For the fast subsystem, a fast stabilizing control is designed according to the standard LQR approach.

2 MFLM Model 2.1 Dynamic Model of MFLM

The structure of MFLM is shown in Fig. 1. It consists of n flexible links, n rigid revolute joints, and a payload. Assume the MFLM is operated horizontally. The i-th flexible link is with length l_i, linear density ρ_i, and rigidity EI_i. The driving torque τ_i is generated by the motors of the i-th joint. The rotation angle of the i-th joint is denoted as θ_i. The joints are considered as lumped-mass and their inertias are ignored. Thus, the i-th joint plays the role of a payload of the (i-1)-th link. The distal link is connected with a payload. The payload mass of the i-th link is m_i.

Utilizing Euler-Lagrange formulation and assumed mode method, the dynamics of the MFLM is modeled as^{[1, 3-4]}

$ \mathit{\boldsymbol{M}}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{q}}} \right)\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot \theta }}}\\ {\mathit{\boldsymbol{\ddot q}}} \end{array}} \right] + \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right) + \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&\mathit{\boldsymbol{K}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{\theta }}\\ \mathit{\boldsymbol{q}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{\tau }}\\ 0 \end{array}} \right] $

(1)

where M(θ, q) and K represent the matrixes of mass and stiffness; $ \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\theta }}{\rm{, }}\mathit{\boldsymbol{\dot \theta }}{\rm{, }}\mathit{\boldsymbol{q}}{\rm{, }}\mathit{\boldsymbol{\dot q}}} \right)$ represents the mass, the summation of the Coriolis, centrifugal, and viscous damping torques; θ^T=[θ₁, θ₂, …, θ_n] is the angular position vector; q^T=[q₁^T, q₂^T, …, q_n^T] with q_i^T=[q_i1, q_i2, …, q_im] being the modal displacement vector of the i-th link; q_ij is the j-th modal displacement of the i-th link; τ^T=[τ₁, τ₂, …, τ_n] is the input torque vector; 0 represents zero block matrix of appropriate dimensions.

The mass matrix M(θ, q) is symmetric and positive defined. Then M(θ, q) and its inverse N(θ, q) is partitioned as:

$ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_{\theta \theta }}}&{{\mathit{\boldsymbol{M}}_{\theta q}}}\\ {\mathit{\boldsymbol{M}}_{\theta q}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{\theta q}}} \end{array}} \right],\mathit{\boldsymbol{N}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}_{\theta \theta }}}&{{\mathit{\boldsymbol{N}}_{\theta q}}}\\ {\mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}}}&{{\mathit{\boldsymbol{N}}_{\theta q}}} \end{array}} \right] $

satisfying:

$ \mathit{\boldsymbol{M}}_{\theta \theta }^{ - 1} = {\mathit{\boldsymbol{N}}_{\theta \theta }} - {\mathit{\boldsymbol{N}}_{\theta q}}\mathit{\boldsymbol{N}}_{qq}^{ - 1}\mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}} $

By partitioning $ \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\theta }}{\rm{, }}\mathit{\boldsymbol{\dot \theta }}{\rm{, }}\mathit{\boldsymbol{q}}{\rm{, }}\mathit{\boldsymbol{\dot q}}} \right)$ into:

$ \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right) = \left[ \begin{array}{l} {\mathit{\boldsymbol{H}}_\theta }\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right)\\ {\mathit{\boldsymbol{H}}_q}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},\mathit{\boldsymbol{q}},\mathit{\boldsymbol{\dot q}}} \right) \end{array} \right] $

The dynamics in Eq.(1) is rewritten as:

$ \mathit{\boldsymbol{\ddot \theta }} = {\mathit{\boldsymbol{N}}_{\theta \theta }}\left( {\mathit{\boldsymbol{\tau }} - {\mathit{\boldsymbol{H}}_\theta }} \right) - {\mathit{\boldsymbol{N}}_{\theta q}}\left( {\mathit{\boldsymbol{Kq}} + {\mathit{\boldsymbol{H}}_q}} \right) $

$ \mathit{\boldsymbol{\ddot q}} = \mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}}\left( {\mathit{\boldsymbol{\tau }} - {\mathit{\boldsymbol{H}}_\theta }} \right) - {\mathit{\boldsymbol{N}}_{qq}}\left( {\mathit{\boldsymbol{Kq}} + {\mathit{\boldsymbol{H}}_q}} \right) $

2.2 Two-time Scale Singular Perturbation Model

A perturbation parameter μ is defined as^[17]:

$ {\mu ^{ - 1}} = \sqrt {{\mu _{\min }}\left( {{\mathit{\boldsymbol{N}}_{qq,0}}\mathit{\boldsymbol{K}}} \right)} $

(2)

where, N_{qq, 0}= N_qq|_{q =0} and μ_min(N_{qq, 0}K) are the smallest eigenvalues of the matrix (N_qqK) when θ varies in the specified operation range and q=0, i.e. there is no deformation in flexible links. In general, μ depends not only on the stiffness of links, but also on the configuration of the MFLM, namely θ. The MFLM tends to be a rigid manipulator if μ→0. If the bandwidth required in the angular position tracking system is much smaller than the first flexible mode, i.e. μ is actually small, the rigid body dynamics can be separated from the flexible dynamics. Scale K with the factor μ² as $ \mathit{\boldsymbol{\tilde K}}{\rm{ = }}{\mathit{\mu }^{\rm{2}}}\mathit{\boldsymbol{K}}$ and define a new state variables z=q/μ². Then the system in Eq.(1) is expressed in the two-time scale form:

$ \begin{array}{l} \mathit{\boldsymbol{\ddot \theta }} = {\mathit{\boldsymbol{N}}_{\theta \theta }}\left[ {\mathit{\boldsymbol{\tau }} - {\mathit{\boldsymbol{H}}_\theta }\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},{\mu ^2}\mathit{\boldsymbol{z}},{\mu ^2}\mathit{\boldsymbol{\dot z}}} \right)} \right] - \\ \;\;\;\;\;\;{\mathit{\boldsymbol{N}}_{\theta q}}\left[ {\mathit{\boldsymbol{\tilde Kz}} + {\mathit{\boldsymbol{H}}_q}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},{\mu ^2}\mathit{\boldsymbol{z}},{\mu ^2}\mathit{\boldsymbol{\dot z}}} \right)} \right] \end{array} $

(3)

$ \begin{array}{l} {\mathit{\boldsymbol{\mu }}^2}\mathit{\boldsymbol{\ddot z}} = \mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}}\left[ {\mathit{\boldsymbol{\tau }} - {\mathit{\boldsymbol{H}}_\theta }\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},{\mu ^2}\mathit{\boldsymbol{z}},{\mu ^2}\mathit{\boldsymbol{\dot z}}} \right)} \right] - \\ \;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{N}}_{qq}}\left[ {\mathit{\boldsymbol{\tilde Kz}} + {\mathit{\boldsymbol{H}}_q}\left( {\mathit{\boldsymbol{\theta }},\mathit{\boldsymbol{\dot \theta }},{\mu ^2}\mathit{\boldsymbol{z}},{\mu ^2}\mathit{\boldsymbol{\dot z}}} \right)} \right] \end{array} $

(4)

Let μ=0, then the boundary layer correction and the slow time scale dynamics are derived as:

$ \begin{array}{l} {{\mathit{\boldsymbol{\ddot \theta }}}_{\rm{s}}} = {\mathit{\boldsymbol{N}}_{\theta \theta ,{\rm{s}}}}\left[ {{\mathit{\boldsymbol{\tau }}_{\rm{s}}} - {\mathit{\boldsymbol{H}}_{\theta ,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right] - \\ \;\;\;\;\;\;{\mathit{\boldsymbol{N}}_{\theta q,{\rm{s}}}}\left[ {\mathit{\boldsymbol{\tilde K}}{\mathit{\boldsymbol{z}}_{\rm{s}}} + {\mathit{\boldsymbol{H}}_{q,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right] \end{array} $

(5)

$ \begin{array}{l} \mathit{\boldsymbol{0}} = \mathit{\boldsymbol{N}}_{\theta q,{\rm{s}}}^{\rm{T}}\left[ {{\mathit{\boldsymbol{\tau }}_{\rm{s}}} - {\mathit{\boldsymbol{H}}_{\theta ,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right] - \\ \;\;\;\;\;{\mathit{\boldsymbol{N}}_{qq,{\rm{s}}}}\left[ {\mathit{\boldsymbol{\tilde K}}{\mathit{\boldsymbol{z}}_{\rm{s}}} + {\mathit{\boldsymbol{H}}_{q,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right] \end{array} $

(6)

where the subscript 's' represents the slow time scale quantities, i.e. for μ=0. The slow manifold is solved from Eq.(6)as:

$ \begin{array}{l} {\mathit{\boldsymbol{z}}_s} = {{\mathit{\boldsymbol{\tilde K}}}^{ - 1}}\left[ {\mathit{\boldsymbol{N}}_{qq,{\rm{s}}}^{ - 1}\mathit{\boldsymbol{H}}_{\theta q,{\rm{s}}}^{\rm{T}}\left( {{\mathit{\boldsymbol{\tau }}_{\rm{s}}} - {\mathit{\boldsymbol{H}}_{\theta ,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right) - } \right.\\ \;\;\;\;\;\;\;\left. {{\mathit{\boldsymbol{H}}_{q,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \right] \end{array} $

(7)

Substituting Eq.(7) into Eq.(5) yields the slow subsystem as:

$ {{\mathit{\boldsymbol{\ddot \theta }}}_{\rm{s}}} = - \mathit{\boldsymbol{M}}_{\theta \theta ,{\rm{s}}}^{ - 1}{\mathit{\boldsymbol{H}}_{\theta ,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\rm{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right) + \mathit{\boldsymbol{M}}_{\theta \theta ,{\rm{s}}}^{ - 1}{\mathit{\boldsymbol{\tau }}_{\rm{s}}} $

(8)

Letting the slow variables be constants and defining the fast variable and the fast control to be z_f:= z-z_s and τ_f:= τ-τ_s, the fast subsystem in Eq.(4) becomes:

$ \frac{{{{\rm{d}}^2}{z_{\rm{f}}}}}{{{\rm{d}}{\sigma ^2}}} = \mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}}{\mathit{\boldsymbol{\tau }}_f} - {\mathit{\boldsymbol{N}}_{qq}}\mathit{\boldsymbol{\tilde K}}{\mathit{\boldsymbol{z}}_f} $

where σ=t/μ is the fast time scale.

3 Design of the Composite Control Scheme

The composite controller is illustrate in Fig. 2. Then, the control torque τ is given by:

$ \mathit{\boldsymbol{\tau }} = {\mathit{\boldsymbol{\tau }}_{\rm{s}}} + {\mathit{\boldsymbol{\tau }}_f} $

where, the slow control τ_s is derived via ADRC method and the fast control τ_f is designed based on LQR approach.

3.1 Active Disturbance Rejection Controller

The slow subsystem in Eq.(8)can be reformulated as:

$ {{\mathit{\boldsymbol{\ddot \theta }}}_{\mathit{s}}} = \mathit{\boldsymbol{f}}\left( \cdot \right) + \mathit{\boldsymbol{u}} $

(9)

where

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{f}}\left( \cdot \right) = {{\left[ {\begin{array}{*{20}{c}} {{f_1}\left( \cdot \right)}&{{f_2}\left( \cdot \right)}& \cdots &{{f_n}\left( \cdot \right)} \end{array}} \right]}^{\rm{T}}} = }\\ { - \mathit{\boldsymbol{M}}_{\theta \theta ,{\rm{s}}}^{ - 1}{\mathit{\boldsymbol{H}}_{\theta ,{\rm{s}}}}\left( {{\mathit{\boldsymbol{\theta }}_{\mathit{s}}},{{\mathit{\boldsymbol{\dot \theta }}}_{\rm{s}}},0,0} \right)} \end{array} $

represents the total disturbance, namely both the internal and external uncertainties,

$ \mathit{\boldsymbol{u}} = {\left[ {\begin{array}{*{20}{c}} {{u_1}}&{{u_2}}& \cdots &{{u_n}} \end{array}} \right]^{\rm{T}}} = \mathit{\boldsymbol{M}}_{\theta \theta ,{\rm{s}}}^{ - 1}{\mathit{\boldsymbol{\tau }}_{\rm{s}}} $

(10)

is the virtual control.

The model is extended with an additional state to represent the total disturbance, which is estimated via ESO and further compensated.

Let x_i1=θ_{i, s} and x_i2= $ {{\dot \theta }_{\mathit{i}, {\rm{s}}}}$, the i-th subsystem of the slow system in Eq.(9) can be rewritten in state space form as:

$ \left\{ \begin{array}{l} {{\dot x}_{i1}} = {x_{i2}}\\ {{\dot x}_{i2}} = {f_i}\left( \cdot \right) + {u_i} \end{array} \right. $

(11)

Assuming that f_i(·) is m times differentiable and introducing augmented states x_i3=f_i(·), x_i4=$ {{\mathit{\dot f}}_\mathit{i}}\left( \cdot \right)$, …, x_i(2+m)=f_i^(m-1)(·), the plant is extended as:

$ \left\{ \begin{array}{l} {{\dot x}_{i1}} = {x_{i2}}\\ {{\dot x}_{i2}} = {x_{i3}} + {u_i}\\ {{\dot x}_{i3}} = {{\dot f}_i}\left( \cdot \right)\\ {{\dot x}_{i4}} = {{\ddot f}_i}\left( \cdot \right)\\ \; \vdots \\ {{\dot x}_{i\left( {2 + m} \right)}} = f_i^{\left( m \right)}\left( \cdot \right) \end{array} \right. $

Assuming that $ \mathit{\dot f}_\mathit{i}^{\left( \mathit{m} \right)}\left( \cdot \right)=0$, an m order linear ESO for the above system can be designed to estimate states x_i1 and x_i2, disturbance and its derivatives f_i(·), $ {{\mathit{\dot f}}_\mathit{i}}\left( \cdot \right)$, …, f_i^(m-1)(·), as:

$ \left\{ \begin{array}{l} {{\dot {\hat x}}_{i1}} = {{\hat x}_{i2}} - {\beta _{i1}}{{\hat e}_i}\\ {{\dot {\hat x}}_{i2}} = {{\hat x}_{i3}} - {\beta _{i2}}{{\hat e}_i} + {u_i}\\ {{\dot {\hat x}}_{i3}} = {{\hat x}_{i4}} - {\beta _{i3}}{{\hat e}_i}\\ {{\dot {\hat x}}_{i4}} = {{\hat x}_{i5}} - {\beta _{i4}}{{\hat e}_i}\\ \; \vdots \\ {{\dot {\hat x}}_{i\left( {2 + m} \right)}} = - {\beta _{i\left( {2 + m} \right)}}{{\hat e}_i} \end{array} \right. $

where $ {{\mathit{\hat x}}_{\mathit{i}{\rm{1}}}}$ and $ {{\mathit{\hat x}}_{\mathit{i}{\rm{2}}}}$ are estimates of the i-th link's angular position and its velocity, respectively; $ {{\mathit{\hat x}}_{\mathit{i}{\rm{3}}}}$, …, $ {{\mathit{\hat x}}_{\mathit{i}\left( {2 + \mathit{m}} \right)}}$ are the estimates of the total disturbance and its derivatives, respectively; β_i1, β_i2, …, β_i(2+m) are the observer gains, $ {{\mathit{\hat e}}_\mathit{i}} = {\mathit{x}_{\mathit{i}1}} - {{\mathit{\hat x}}_{\mathit{i}{\rm{1}}}}$ stands for the estimation error of the state x_i1.

The ESO gains can be obtained by comparing the roots of the characteristic polynomial,

$ \lambda \left( s \right) = {s^{2 + m}} + {\beta _{i1}}{s^{2 + m - 1}} + \cdots + {\beta _{i\left( {2 + m - 1} \right)}}s + {\beta _{i\left( {2 + m} \right)}} $

to the following polynomial

$ G\left( s \right) = {\left( {s + {\omega _0}} \right)^{2 + m}} $

Which placed all the observer poles in the left half plane at -ω₀, making the characteristic polynomial a Hurwitz-type^[24].

Assuming that $ {{\mathit{\hat x}}_{\mathit{i}{\rm{3}}}}$ estimates f_i(·) closely, i.e. $ {{\mathit{\hat x}}_{\mathit{i}{\rm{3}}}}$=f_i(·), and substituting the selected virtual control input

$ {u_i} = - {{\hat x}_{i3}} + {{\bar u}_i} $

(12)

into Eq.(11), the i-th subsystem becomes a double integrator system

$ \left\{ \begin{array}{l} {{\dot x}_{i1}} = {x_{i2}}\\ {{\dot x}_{i2}} = {{\bar u}_i} \end{array} \right. $

Then, a simple controller combining PD feedback and a feedforward term can be designed as:

$ {{\bar u}_i} = {{\ddot \theta }_{{\rm{d}}i}} + {k_{{\rm{p}}i}}{e_{{\rm{p}}i}} + {k_{{\rm{d}}i}}{e_{{\rm{d}}i}} $

(13)

where

$ {e_{{\rm{p}}i}} = {\theta _{{\rm{d}}i}} - {{\hat x}_{i1}},{e_{{\rm{d}}i}} = {{\dot e}_{{\rm{p}}i}} = {{\dot \theta }_{{\rm{d}}i}} - {{\hat x}_{i2}} $

with θ_di being the desired angular position of the i-th link, k_pi and k_di being the proportional and derivative gains, respectively. Under the assumption of perfect estimation, i.e. $ {{\mathit{\hat x}}_{\mathit{i}{\rm{1}}}}$=θ_i and $ {{\mathit{\dot {\hat x}}}_{\mathit{i}{\rm{1}}}} = {{\dot \theta }_\mathit{i}}$, one has the following error dynamic equation:

$ {{\ddot e}_{{\rm{p}}i}} + {k_{{\rm{d}}i}}{{\dot e}_{{\rm{p}}i}} + {k_{{\rm{p}}i}}{e_{{\rm{p}}i}} = 0 $

(14)

Thus, k_pi and k_di can be designed to make the above system stable and shape its response.

The desired angular position, θ_di, can be measured directly from the encoder placed on the motor shaft. However, the rate and acceleration of the angular position, $ {{\dot \theta }_{{\rm{d}}\mathit{i}}}$ and $ {{\ddot \theta }_{{\rm{d}}\mathit{i}}}$, need to be estimated from a differentiator. Herein, the TD filter is employed. The discrete form TD filter is given by^[21]:

$ \left\{ \begin{array}{l} {x_1}\left( {k + 1} \right) = {x_1}\left( k \right) + h{x_2}\left( k \right)\\ {x_2}\left( {k + 1} \right) = {x_2}\left( k \right) + h{f_{{\rm{fh}}}}\left( {{x_1} - u,{x_2},r,{h_0}} \right) \end{array} \right. $

where x₁, x₂, and u represent state variables and the control input, respectively; h is the step size; h₀ and r are designed filter and speed factor, respectively; and the function f_fh(x₁, x₂, r, h₀) is defined as:

$ \left\{ \begin{array}{l} d = r{h_0}\\ {d_0} = {h_0}d\\ y = {x_1} + {h_0}{x_2}\\ {a_0} = \sqrt {{d^2} + 8r\left| y \right|} \\ a = \left\{ \begin{array}{l} {x_2} + 0.5\left( {{a_0} - d} \right){\rm{sign}}\left( y \right),\;\;\;\;\left| y \right| > {d_0}\\ {x_2} + y/{h_0},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| y \right| \le {d_0} \end{array} \right.\\ {f_{{\rm{fh}}}} = - \left\{ \begin{array}{l} r{\rm{sign}}\left( a \right),\;\;\;\;\left| a \right| > d\\ ra/d,\;\;\;\;\;\;\;\;\;\left| a \right| \le d \end{array} \right. \end{array} \right. $

Remark 3.1  Thanks to the separation principle, the parameters tuning for the TD filer, the ESO, and the PD controller can be taken place independently. To make the separation principle holds, the bandwidths of the TD filter and the ESO need to be relatively higher than that of the PD controller.

Replace the real rate and acceleration of the angular position by the estimated ones, which are denoted as $ {{{\hat{\dot{\theta }}}}_{\rm{d}\mathit{i}}}$ and $ {{{\hat{\ddot{\theta }}}}_{\rm{d}\mathit{i}}}$. Then, substituting Eq.(13) into Eq.(12) yields the virtual control of the i-th subsystem as:

$ {u_i} = \left( {{{\hat {\ddot \theta} }_{{\rm{d}}i}} - {{\hat x}_{i3}}} \right) + {k_{{\rm{d}}i}}\left( {{{\hat {\dot \theta }}_{{\rm{d}}i}} - {{\hat x}_{i2}}} \right) + {k_{{\rm{p}}i}}\left( {{{\hat \theta }_{{\rm{d}}i}} - {{\hat x}_{i1}}} \right) $

Then the control input of the slow subsystem can be calculated from Eq.(10) as:

$ {\mathit{\boldsymbol{\tau }}_{\rm{s}}} = {\mathit{\boldsymbol{M}}_{\theta \theta ,{\rm{s}}}}\mathit{\boldsymbol{u}} $

The terms e_pi and e_di in Eq.(13) can also be replaced by a nonlinear function f_fal(x, α, δ)^[21]:

$ {f_{{\rm{fal}}}}\left( {x,\alpha ,\delta } \right) = \left\{ \begin{array}{l} {\left| x \right|^\alpha }{\rm{sign}}\left( x \right),\;\;\;\;\left| x \right| > \delta \\ \frac{x}{{{\delta ^{1 - \alpha }}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\left| x \right| \le \delta \end{array} \right. $

Where δ > 0, 0 < α < 1. Then a NPD controller is obtained as:

$ {{\bar u}_i} = {{\ddot \theta }_{{\rm{d}}i}} + {k_{{\rm{p}}i}}{f_{{\rm{fal}}}}\left( {{e_{{\rm{p}}i}},{\alpha _{{\rm{p}}i}},{\delta _{{\rm{p}}i}}} \right) + {k_{{\rm{d}}i}}{f_{{\rm{fal}}}}\left( {{e_{{\rm{d}}i}},{\alpha _{{\rm{d}}i}},{\delta _{{\rm{d}}i}}} \right) $

Remark 3.2  As shown in Fig. 3, for small error, the function f_fal(x, α, δ) provides high gain. In the neighborhood of origin the gain is limited by choosing a small δ. By tuning the parameters α and δ, different gains are obtained as tracking error varies.

3.2 LQR Controller for the Fast Subsystem

The state space representation of the fast subsystem is:

$ {{\mathit{\boldsymbol{\dot x}}}_f} = {\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{x}}_f} + {\mathit{\boldsymbol{B}}_f}{\mathit{\boldsymbol{u}}_f} $

where

$ {\mathit{\boldsymbol{A}}_f} = \left[ {\begin{array}{*{20}{c}} 0&\mathit{\boldsymbol{I}}\\ { - {\mathit{\boldsymbol{N}}_{qq}}\mathit{\boldsymbol{\tilde K}}}&0 \end{array}} \right],{\mathit{\boldsymbol{B}}_f} = \left[ {\begin{array}{*{20}{c}} 0\\ {\mathit{\boldsymbol{N}}_{\theta q}^{\rm{T}}} \end{array}} \right],{\mathit{\boldsymbol{x}}_f} = \left[ {\begin{array}{*{20}{c}} {{z_f}}\\ {\frac{{{\rm{d}}{z_f}}}{{{\rm{d}}\sigma }}} \end{array}} \right] $

Since the pair (A_f, B_f) is completely controllable, optimal control technique can be applied to achieve a fast vibration depression. The object function is chosen as:

$ J = \frac{1}{2}\int_0^\infty {\left[ {\mathit{\boldsymbol{x}}_f^{\rm{T}}{\mathit{\boldsymbol{Q}}_f}{\mathit{\boldsymbol{x}}_f} + \mathit{\boldsymbol{\tau }}_f^{\rm{T}}{\mathit{\boldsymbol{R}}_f}{\mathit{\boldsymbol{\tau }}_f}} \right]{\rm{d}}t} $

where Q_f and R_f are weighting matrices. Following the standard LQR design approach, the fast control is derived as:

$ {\mathit{\boldsymbol{\tau }}_f} = - {\mathit{\boldsymbol{K}}_f}{\mathit{\boldsymbol{x}}_f} + {\mathit{\boldsymbol{K}}_f}\left[ {\begin{array}{*{20}{c}} {{z_{\rm{s}}}}\\ 0 \end{array}} \right] $

(15)

where

$ \mathit{\boldsymbol{K}} = - \mathit{\boldsymbol{R}}_f^{ - 1}\mathit{\boldsymbol{B}}_f^{\rm{T}}{\mathit{\boldsymbol{P}}_f}\left( t \right) $

and P_f is solved from the Riccati equation:

$ {\mathit{\boldsymbol{P}}_f}\left( t \right) = - {\mathit{\boldsymbol{P}}_f}\left( t \right){\mathit{\boldsymbol{A}}_f} - \mathit{\boldsymbol{A}}_f^{\rm{T}}{\mathit{\boldsymbol{P}}_f}\left( t \right) + {\mathit{\boldsymbol{P}}_f}\left( t \right){\mathit{\boldsymbol{B}}_f}\mathit{\boldsymbol{R}}_f^{ - 1}\mathit{\boldsymbol{B}}_f^{\rm{T}}{\mathit{\boldsymbol{P}}_f}\left( t \right) - {\mathit{\boldsymbol{Q}}_f} $

Remark 3.3  The first term in the fast controller is derived from LQR approach, while the second term is employed to compensate for the slow control and referred to as the boundary layer correction (BLC).

Implementing the full state feedback control in Eq.(15) requires z_f to be available. Thus the links' deflections and their time derivatives, i.e. q and $ {\mathit{\boldsymbol{\dot{q}}}}$, have to be measured or estimated. The links' deflections can be constructed via measurements of strain gage.Their rate $ {\mathit{\boldsymbol{\dot{q}}}}$ is estimated via the TD filter.

4 A Case Study 4.1 Simulation Parameters

The proposed control scheme is verified via numerical study. The manipulator considered has two flexible links and each link has two assumed modes. The dynamic model is derived according to Section 2 and using the physical parameters given in Table 1, where m_{1, 0} and m_{2, 0} refer to initial payloads. The fourth-order Runge-Kutta solver is utilized to solve the simulation. The sampling time is 0.000 5 s.

Let the operation range be {(θ₁, θ₂)|θ₁∈[0, 2π], θ₂∈[0, 2π]}, then the perturbation parameter in Eq.(2) is calculated as μ=0.016 7 s. The smallest natural frequency is 59.81 rad/s. The poles is placed at -2.5. Then, the feedback gains are calculated as k_di=5 and k_pi=6.25. With the above parameters, the bandwidth of the closed slow subsystem (without the feedforward term $ {{{\ddot{\theta }}}_{\rm{d}\mathit{i}}}$ is 6.20 rad/s. The observer poles are placed at ω₀=20 rad/s to make the observer relatively fast enough. Then the observer's gains can be calculated as β₁=3ω₀, β₂=3ω₀², β₃=ω₀³. The parameters for the TD filter is set to be r₁₁=r₁₂=3 000, h₁₁₀=h₁₂₀=0.04, r₂₁=1 500, h₂₁₀=0.03, r₂₂=2 000, and h₂₂₀=0.04. The parameters for the NPD controller are chosen as α_p1=0.3, α_d1=0.5, δ_p1=0.2, δ_d1=0.7, α_p2=0.2, α_d2=0.4, δ_p2=0.2, and δ_d2=0.9. The weighting matrices are selected as:

$ {\mathit{\boldsymbol{Q}}_f} = {\rm{diag}}\left( {100,10,10,10,10,10,10,10} \right) $

$ {\mathit{\boldsymbol{R}}_f} = {\rm{diag}}\left( {100,10} \right) $

Then the feedback gain of the fast controller is calculated as:

$ \mathit{\boldsymbol{K}}_f^{\rm{T}} = \left[ {\begin{array}{*{20}{c}} { - 0.034\;2}&{0.046\;5}\\ {0.055\;8}&{ - 0.082\;7}\\ { - 0.040\;1}&{ - 0.033\;1}\\ { - 0.006\;5}&{0.077\;4}\\ { - 0.593\;4}&{0.630\;2}\\ {0.208\;6}&{ - 0.205\;4}\\ { - 0.364\;2}&{ - 0.349\;2}\\ { - 0.111\;4}&{0.293\;3} \end{array}} \right] $

The maneuver trajectory is designed as:

$ {\theta _i} = \left\{ \begin{array}{l} \frac{{{a_i}}}{2}\left[ {\sin \left( {{\omega _i}t - \frac{{\rm{ \mathsf{ π} }}}{2}} \right) + 1} \right],\;\;\;\;{\omega _i}t \le {\rm{ \mathsf{ π} }}\\ {a_i},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\omega _i}t > {\rm{ \mathsf{ π} }} \end{array} \right. $

where i=1, 2, a_i and ω_i are the desired final angular position and maneuver speed, respectively.

4.2 Relatively Slow Speed Maneuver

In the first set of simulations, relatively slow speed maneuver is tested by setting ω₁=ω₂=2 rad/s. The desired angular positions of both joints are a₁=π/3 rad and a₂=π/4 rad. With the initial payload only, the tracking performance of different controllers is shown in Fig. 4. Herein and after, ADRC-PD refers to the controller utilizing ADRC structure and PD feedback controller and ADRC-NPD refers to that with ADRC structure and NPD feedback controller. For comparison, the PD controller proposed in Ref.[14] is designed and results are also presented in Fig. 4. The superiority of the ADRC controllers is significant in reducing tracking error and vibration.Using ADRC-PD, the overshot of θ₁ and θ₂ are 4.75% and 4.82%, respectively. And the settling times are respectively 2.49 s and 2.29 s. The performance of PD controller is worse. The overshots are 8.76% and 13.20%, and the settling times are 2.99 s and 3.17 s, respectively. The ADRC-NPD controller improves the performance slightly. The overshot decrease to 4.74% and 3.54%, and the settling times reduce to 2.05 s and 1.79 s, respectively. Seen from Figs. 4(c)-(d), all the three controllers have similar performance in reducing the oscillations of both links.

If the payload increases, namely m₁=m_{1, 0}+ΔM_p, m₂=m_{2, 0}+ΔM_p, ΔM_p being the increase of the payload, the performance is comprised in Fig. 5. The responses of joints' positions are similar. However, there is a slight vibration in the links, especially for the PD controller.

4.3 Relatively Fast Speed Maneuver

For the second set of simulations, a relatively fast speed maneuver is tested. The parameters for the desired trajectory are ω₁=ω₂=50 rad/s. As shown in Figs. 6(a)-(b), due to bandwidth limitation, all three controllers cannot follow the desired trajectory at the initial stage since the desired angular positions increase rapidly. Similar as the slow speed maneuver case, the ADRC-NPD has the smallest overshot and settling time, i.e. 0.94 s, 3.60% for θ₁ and 0.74 s, 7.24% for θ₂. The PD controller has the worst performance, the overshot for θ₁ and θ₂ are 11.25% and 22.08% and the settling time are 2.14 s and 2.35 s, respectively. Meanwhile, seen from Figs. 6(c)-(d), the ADRC-PD and ADRC-NPD smooth down the oscillations of both links more quickly than the PD controller.

As shown in Fig. 7, with an increase in the payload, the overshot of θ₁ increases slightly for ADRC-PD controller and ADRC-NPD controller, whereas there is no overshot by utilizing the PD controller. However, the response of the PD controller becomes slowly. The reason is that, with these parameters, the system with PD controller becomes an over damping system. Seen from the response of θ₂, the overshot decreases for all the three controllers. However, as seen from the response of link 2, the oscillation is much larger when the PD controller is employed.

4.4 Performance at Different Maneuver Speeds

The performance of the three controllers are also tested at different maneuver speeds and summarized in Fig. 8, where the maneuver speed is drawn by ω=ω₁=ω₂. Generally, the settling time of all the three controllers decreases with ω. Especially, when ω < 10 rad/s, the settling time decreases rapidly. This indicates that the controllers have sufficiently wide bandwidth. However, the maneuver speed is limited by the desired maneuver trajectory. When ω > 10 rad/s, the maneuver speed is limited by the bandwidth of the controllers. Then, a fast desired maneuver trajectory can only increase the maneuver speed slightly. The ADRC-NPD can achieve the smallest settling time while that of the PD controller is the largest. Seen from Figs. 8(c)-(d), for the ADRC-PD and PD controller, a faster maneuver speed generally causes a larger overshot. However, for the ADRC-NPD, the overshot increases with ω firstly and then decreases with it. This is caused by the nonlinear function f_fal(x, α, δ), which yields a large gain when x is small. When it is added to the D-term of the controller, a larger damping is obtained even when the rate of the tracking error is small, see Fig. 3. Compared with PD controller, both of the ADRC-PD and ADRC-NPD controllers can realize a smaller overshot.

To further compare the tracking performance of three controllers, the index mean square error (MSE) is introduced as:

$ M\left( x \right) = \sqrt {\frac{1}{{{t_{\rm{f}}}}}\int_0^{{t_{\rm{f}}}} {{{\left( {{x_{\rm{d}}} - x} \right)}^2}{\rm{d}}t} } = \sqrt {\frac{1}{N}\sum\limits_N {{{\left( {{x_{\rm{d}}} - x} \right)}^2}} } $

(16)

where t_f is the total time of the simulation and N is the number of samples. The MSEs for three controllers at different maneuver speed are shown in Fig. 9. Generally, the MSEs of three controllers increase with ω and ADRC-PD and ADRC-NPD have similar performance. For θ₂, the performance of ADRC-PD and ADRC-NPD is always better than that of PD controller. However, for θ₁, the ADRC-PD and ADRC-NPD controller has smaller MSE only when the maneuver speed is small. As the desired maneuver speed increases, the bandwidth limit of ESO and TD filter will cause an obvious delay at the initial stage when the desired angular position increases rapidly, see Fig. 6. As shown in Fig. 10 and noting the squire operator in Eq.(16), the error caused by the delay at the initial stage results in a larger MSE index.

The performances of these controllers with different payloads are shown in Fig. 11. For the relative slow maneuver, i.e. ω≤10 rad/s, the MSEs of θ₁ and θ₂ is robust to the payload mass for all the three controller. However, the MSEs of θ₁ and θ₂ increase with maneuver speed. For the response of θ₁, the ADRC-NPD controller has a better performance if the maneuver speed is small but worse than PD controller if the speed increases. The reason has been explained above. However, for the response of θ₁, the ADRC-NPD controller has better performance, no matter what the payload mass and the maneuver speed are.

5 Conclusions and Discussions

This paper focuses on application of the active disturbance rejection control (ADRC) on the multiple-flexible-link manipulator (MFLM). The dynamics of the MFLM is modeled and separated into two-time scale based on singular perturbation technique. The ADRC structure consisting of a tracking differentiator (TD) filter, an extended state observer (ESO), and a feedback controller is applied to the slow subsystem. The TD filter is employed to generate the derivation of the desired trajectory and the ESO is utilized to estimate the states as well as the total disturbance. The feedback controller can be a linear one, proportional derivation (PD) for example, or a nonlinear one, i.e. nonlinear PD controller. The controller for the fast subsystem is designed by following standard linear quadratic regulator approach. The simulation demonstrates the superiority of controllers based on ADRC structure. Generally, the ADRC structure based controllers have a superior performance on settling time, overshot, vibration cancelling, and etc. It is robust to the maneuver speed. However, the stability proof is a difficult task. And more parameters need to be adjusted. In the future, more detailed work will be carried out on the experiment research on the application of ADRC structure based controllers on the MFLM. The method how to adjust the parameters and how to handle the uncertain payload will be focused on in the next step work.