A High-Efficient Finite Difference Method for Flexible Manipulator with Boundary Feedback Control

The paper presents a high-efficient finite difference method for solving the PDE model of the single-link flexible manipulator system with boundary feedback control. Firstly, an abstract state-space model of the manipulator is derived from the original PDE model and the associated boundary conditions of the manipulator by using the velocity and bending curvature of the flexible link as the state variables. Then, the second-order implicit Crank-Nicolson scheme is adopted to discretize the statespace equation, and the second-order one-sided approximation is used to discretize the boundary conditions with excitations and feedback control. At last, the state-space equation combined with the boundary conditions of the flexible manipulator is transformed to a system of linear algebraic equations, from which the response of the flexible manipulator can be easily solved. Numerical simulations are carried out to simulate the manipulator under various excitations and boundary feedback control. The results are compared with ANSYS to demonstrate the accuracy and high efficiency of the presented method.


Introduction
Flexible manipulators are widely used in aerospace engineering for the applications of active debris removal, building large space structures, and servicing satellites, etc. [1][2][3]. In order to meet the requirements of high speed, high accuracy, and low energy consumption, the manipulators tend to be more and more lightweight and highly flexible. As a consequence, the vibration problem has become one of the most important considerations in the design of flexible manipulators [2][3][4]. As complex distributed parameter systems with coupling dynamics between the rigid and the flexible modes [5], the traditional controller design approach which uses the discretized ordinary differential equation (ODE) models is inaccurate to describe the spatiotemporally varying states and may cause spillover instability due to the neglect of higher-order modes [5,6]. Alternatively, the controller design for flexible manipulators using the partial differential equation (PDE) models without discretization has gained significant attention in recent years [5][6][7][8][9][10][11][12]. For example, Jiang et al. [5] designed an infinite dimensional disturbance observer for the boundary control of the flexible manipula-tor. Ge et al. [6] introduced nonlinear strain feedback to improve the performance of the boundary PD controller on the flexible manipulator. He et al. [12] designed a boundary controller based on the adaptive inverse algorithm for the flexible manipulator with input backlash.
In order to verify the controller by a numerical simulation or carry out a real-time state estimation in the control experiment, a high-efficient solving method for the PDE model of the flexible manipulator is in demand. The finite difference (FD) method is widely used for solving the PDE model of the control system [5][6][7][8][9][10][11][12]. As a classical numerical method for solving the PDE system, the FD method has been well developed for simulation of flexible structures such as strings [13,14] and beams [15][16][17]. However, the flexible manipulator incorporating hub inertia and payload is more complex than a single string or beam since it is a coupled PDE-ODE system. Specialized studies on developing the FD method for the single-link flexible manipulator (SLFM) were carried out by Tzes et al. [18] and Tokhi et al. [19][20][21]. They all used the central difference scheme to discretize the original fourth-order PDE model and the associated boundary conditions of the SLFM system. The methods presented by Refs. [18][19][20][21] have been used by many researchers for flexible manipulators in various applications such as response simulation, trajectory tracing, and parametric identification until now [22][23][24][25][26]. However, the above methods all involve iteration procedures at each time step, and the selection of time step Δt and space step Δx must satisfy the stability condition, i.e., ðΔt 2 EIÞ/ðΔx 4 ρAÞ ≤ 1/4 (EI and ρA are the flexural rigidity and the mass per unit length of the flexible link, respectively), which usually leads to a very small time step and a long computational time. Tokhi and Azad [19] investigated the real-time performance of their algorithm and suggested that the MATLAB package is not suitable for real-time simulation of their algorithm because of the need for large iteration loops. Alternatively, they implemented their algorithm on a parallel computation platform to achieve real-time performance [20]. Another method for solving the highorder PDE is to transfer it to a low-order PDE system by using the state variables; then, a wider variety of difference procedures are available, and high computational efficiency could be obtained. For example, Abhyankar et al. [15] used the bending moment as the state variable to reduce the order of the PDE of a nonlinear Euler-Bernoulli beam model and obtained a more efficient FD method. However, this approach has not been applied to the solution of the flexible manipulator till now. This paper presents a new FD method for simulation of the SLFM system by transferring the original PDE model to an abstract state-space model and then using the secondorder accurate and unconditionally stable Crank-Nicolson scheme to solve the state-space model. The presented method achieves a very high efficiency since it avoids the needs of large iteration loops and very small time steps as in the previous studies [18][19][20][21].
The paper is organized as follows: In Section 2, the PDE model and boundary conditions of SLFM with boundary feedback control are given, and an abstract state-space model is derived from the above PDE model by using the velocity and curvature of the flexible link as state variables. In Section 3, the Crank-Nicolson scheme and the second-order accurate one-sided approximation are adopted to discretize the statespace equation and the boundary conditions of the SLFM system. Finally, numerical simulations for the SLFM under excitations and boundary feedback control are presented in Section 4.  Figure 1, which consists of a rigid hub, a flexible link with length L, and a payload at the free end of the link, in which XOY is the inertial coordinate and X 1 OY 1 is the rotating coordinate attached to the link. Suppose that the manipulator only moves in the horizontal plane and a control torque τðtÞ and a control force FðtÞ are applied at the hub and the free end of the manipulator by the joint motor and a force actuator, respectively.
For an angular displacement θ of the hub and a small elastic deflection w of the flexible link, the total displacement y at an arbitrary point of the flexible link can be described as Modeling the flexible link as an Euler-Bernoulli beam with structural damping (using the Chen-Russell damping model) [27], the equation of motion for the flexible link can be derived as where EI and ρA are the flexural rigidity and the mass per unit length of the link, respectively, c s is the coefficient of structural damping, qðx, tÞ is a distributed spatiotemporally varying disturbance load along the flexible link, and the overdot and prime indicate the derivatives with respect to t and x, respectively. From the natural boundary conditions of the flexible link and the equilibrium equations for the rigid hub and the payload, the boundary conditions of the SLFM system can be derived as [6] w 0, t ð Þ= w′ 0, t ð Þ= w″ L, t ð Þ= 0, where I h is the inertia of the hub and m is the mass of the payload. Substituting Equation (1) into Equation (3) yields y 0, t ð Þ= 0, 2.2. Abstract State-Space Equation of the SLFM. Equation (2) is a fourth-order PDE; it can be discretized directly by using  Space: Science & Technology the central difference scheme as in Refs. [18][19][20][21]; however, this approach yields a recursive equation which needs an iterative algorithm. Note that the velocity and bending curvature of the flexible link can be written as Using the variable substitution method, the equation of motion of the manipulator can be written as where a 1 = EI/ρA, a 2 = c s /ρA, and qðx, tÞ = qðx, tÞ/ρA. The boundary conditions in Equation (4) can be expressed with the state variables as v 0, t ð Þ= 0, ð7aÞ The above transformation yields lower-order derivatives in the equations of motion as well as in the boundary conditions, which makes the equations much easier to solve by the finite difference method.
Denoting the state vector xðx, tÞ = fvðx, tÞ, κðx, tÞg T , then Equation (6) can be recast into the following abstract state equation [26]: where The abstract state-space model is widely used in the controller and observer design of distributed parameter systems [28]. Since the state variables in Equation (8) are the velocity and strain of the manipulator, this model is much appropriate for the controller design based on velocity feedback or strain feedback such as the studies in Refs. [29][30][31].

Finite Difference Method for the State-Space
Model of SLFM

Discretization of the Abstract State Equation.
The abstract state Equation (8) is a parabolic PDE system when the damping coefficient c s > 0 or a Petrovskii parabolic sys-tem when c s = 0 [32]. For solving this equation with the boundary condition (7a), (7b), (7c), (7d), the implicit and unconditionally stable Crank-Nicolson difference scheme is suitable [32]. The Crank-Nicolson scheme can be used to solve PDEs that are first order in time and arbitrary order in space with mixed partial derivatives; it applies central differencing to the PDEs about the point (x, t + Δt/2) and evaluates the time derivate and the space derivate of an arbitrary function uðx, tÞ as [32] _ u x, t It has second-order accuracy in both time and space. This stability property can be utilized to take large time steps in the simulation.
For the abstract state Equation (8), discretizing the spatial variable as x m = mΔx (m = 0, 1, 2, ⋯, N + 1), and the time as t n = nΔt (n = 0, 1, 2, ⋯), then by utilizing the Crank-Nicolson scheme, it can be discretized as where x n m ≡ xðx m , t n Þ, q n m ≡ qðx m , t n Þ. Letting μ = Δt/Δx 2 , Equation (11) can be rewritten as where I is a 2 × 2 identity matrix. At the time step n + 1, the state vector at the previous time step and the load vector at the current time step are assumed known, so Equation (12) can be simplified as where is known at the current time step.
3 Space: Science & Technology By discretizing the state equation at each interior point x m (m = 1, 2, ⋯, N), the following system of equations can be obtained: Equation (15) is a system of linear algebraic equations containing 2N equations and 2N + 4 unknown variables, which cannot be solved directly. However, there are four unknown variables that can be determined directly or calculated in terms of their neighboring values from the four boundary conditions in Equation (7a), (7b), (7c), (7d). (7c) give

Discretization of the Boundary Conditions. The boundary conditions (7a) and
According to the Crank-Nicolson scheme, the boundary condition (7b) is differenced about the point ð0, t n + Δt/2Þ, which yields It follows from Equation (18) that Furthermore, calculating the spatial derivative at the boundary x = 0 by using the second-order accurate onesided approximation method yields [32] v′ Substituting Equation (20) into (19) gives where Taking a similar difference procedure for the boundary condition (7d) yields Noticing that the spatial derivative at the right boundary is approximated by [30], At last, the difference formula for the boundary condition (7d) can be obtained as where b 2 = EIΔt/mΔx.

Solution of the System of Equations.
For making use of the boundary condition at the left end of the flexible link, the first equation in Equation (15) is rewritten as where t 1 = f1, 0g T , t 2 = f0, 1g T . Substituting the boundary conditions (16) and (20) into Equation (25) yields where For making use of the boundary condition at the right end of the flexible link, the last equation in Equation (15) is rewritten as Substituting the boundary conditions (17) and (24) into Equation (28) gives With the above transformation, the systems of equations in Equation (15) can be replaced by new systems of equations which contains 2N equations and 2N unknown variables as where Equation (32) is a tridiagonal system of equations, which can be solved directly by evaluating the inverse of matrix A d or be solved more efficiently by using the Thomas algorithm [32].
From the velocity of an arbitrary point at the flexible link, the displacement y with second-order accuracy can be obtained by using the average velocity as [15] y n+1 m = y n m + Furthermore, the hub angle can be obtained by using Equation (1) and the second-order accurate one-sided approximation as

Solution of the SLFM with Boundary Feedback Control.
Considering the position and vibration control of the SLFM using linear boundary feedback control, the proportional derivative and strain (PDS) control [29] is selected. The deriving torque is obtained as where θ d is the designed angular position of the manipulator, k p and k v are the feedback gains of the PD controller, and k e is the strain feedback gain. Substituting Equation (37) into the boundary condition (7b) yields Differencing Equation (38) about the point (0, t n + Δt/2) yields From Equation (39) and utilizing Equation (35), it can be obtained that Calculating the spatial derivative at the boundary x = 0 in Equation (40) by the second-order accurate one-sided approximation method as in Equation (20) yields where Substituting Equation (41) into Equation (25) yields where in Equations (32) and (34) by using Equations (44), (45), and (46), the system of equations for the SLFM with boundary feedback control is obtained.

Results and Discussions
To verify the correctness and efficiency of the presented method, a numerical study is performed. The system  The damping coefficient c s is assumed to be 0.05 N·s. Since the presented method is based on the unconditional stable Crank-Nicolson scheme, it has no stability requirement on the time step and the space step. The selection of these two simulation parameters is determined by the excitation frequency and the requirement of accuracy.
The presented method is carried out in the MATLAB platform. Its accuracy and performance are verified by comparing with the finite element method. The finite element model of the flexible manipulator is established in ANSYS, in which the flexible link is discretized by planar beam (Beam3) elements and the rigid hub and the payload are modelled as MASS21 elements. The mesh size of the flexible link in the ANSYS model is identical to the space step used in the finite difference method, and the dynamic analysis in ANSYS uses the same time step as in the finite difference method. Since the Chen-Russell damping model is difficult to model in ANSYS, during comparison of the presented method and ANSYS, the damping is neglected.
To demonstrate the correctness of the presented method for SLFM under various loading conditions, both of the low-frequency excitation and the high-frequency excitation are applied simultaneously on the manipulator. For the low-frequency excitation, a driving torque τðtÞ = 5sinð2π · 5tÞ N·m and a distributed disturbance load qðx, tÞ = 0:5 sin ð2π · 10tÞ N/m are applied. For the high-frequency excitation, a broadband impulse load FðtÞ is applied on the payload with [33] where HðtÞ is the Heaviside function and f is the bandwidth. In this simulation, F 0 = 10 N and f = 1 kHz are selected; the corresponding impulse load is shown in Figure 2. The time step and space step for this simulation are set as Δt = 1 × 10 −4 s and Δx = L/100 = 0:01 m, respectively. The responses of the flexible manipulator under the above excitations are shown in Figure 3. Comparing the results of the presented method and ANSYS for the undamped flexible manipulator, it can be found that they agree very well, which verifies the correctness of the presented method for simulating the SLFM under both lowand high-frequency excitations. Furthermore, it can be found that the damping of the flexible link has a negligible influence on the hub angel but has an obvious effect on the vibration of the flexible link.
The comparison of computational efficiencies of the presented method and ANSYS using different simulation parameters is listed in Table 1. It can be found that the computational times of the presented method are much less than those of ANSYS in all simulation cases, which shows the extremely high efficiency of the presented method. Moreover, it can be found that if a large time step is used when there is no high-frequency impulse load applied on the flexible manipulator, the computational time of the presented method is smaller than the simulation time, which indicates that the presented method can satisfy the requirement of real-time simulation of the flexible manipulator if no highfrequency loads are applied.
Next, the position and vibration control of the SLFM using the boundary controller in Equation (39) are simulated. The desired angular position of manipulator is θ d = 0:5 rad. The gains of the PD controller are set as k p = 14:58 N·m/rad and k v = 5:83 N·m·s/rad [6]; the strain feedback gain is selected as k e = −EI = −2 N·m 2 . The time step for this simulation is selected as Δt = 1 × 10 −3 s. In order to verify the results of the presented method by ANSYS, the damping of the flexible link is neglected here. The simulation results of the hub angle, the tip displacement, the deflection of the flexible link,   Figure 4. It can be found that the manipulator is driven to the desired position and the vibration of the manipulator is suppressed effectively by both of the two controllers, and the PDS controller can obtain better performance than the PD controller due to the incorporation of strain feedback.
To verify the results of the SLFM with boundary control evaluated by the presented method, the driving torque τðtÞ of the PDS control obtained by the presented method is applied to the ANSYS model, and the response of the manipulator under control was evaluated. The hub angle and tip deflection evaluated by ANSYS are shown in Figure 5. It can be found that the results of the ANSYS model agree well with those of the presented method, which demonstrates the correctness of the presented method for simulating the SLFM with boundary control.

Conclusions
In this research, a high-efficient finite difference method for solution of the PDE model of single-link flexible manipulators with boundary feedback control was proposed based on the unconditional stable Crank-Nicolson scheme. Numerical simulations verify the correctness and extremely high computational efficiency of the present method for simulating the single-link flexible manipulators under both low-and high-frequency excitations and the PDS (proportional derivative and strain) boundary feedback control. The presented method can satisfy the requirement of real-time simulation for the single-link flexible manipulators if no high-frequency excitations are applied. In the future, the extension of the presented method to the solution of more sophisticated control systems such as the nonlinear rigidflexible coupled models with advanced controllers can be studied [34,35]. Meanwhile, the application of the presented method in real-time state estimation of the flexible manipulator can be investigated.

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

Conflicts of Interest
The authors declared that they have no conflicts of interest to this work.