Control of Flexible Robot by Harmonic Functions

Abstract

This work deals with the control of flexible structures as underactuated systems. The invariant control method performs the control of a flexible robot as a representative of an underactuated system with zero dynamics. The control input is separated into two parts. The arbitrary part of the control input is designed to control the directly actuated part of the dynamic system. The invariant part of the control is selected to steer the system zero dynamics in the desired way. The harmonic functions create the base for the invariant part of the control function. The residual vibration cancellation is the target of the presented invariant control strategy. The harmonic function frequencies are overtaken from the so-called natural motion, amplitudes are the results of the optimization process. The main target of this paper is to show the invariant control approach and its application to the system with flexible elements.

1. Introduction

This paper deals with the control of flexible structures as underactuated systems. The definition of underactuation should be based on the number of blocks in the Brunovsky canonical form [1]. The invariant control—its idea is presented in [2]—is used for flexible robot trajectory tracking with minimal residual oscillations. The desired trajectory is parametrized, exact input–output linearization is applied to the system coordinates (modal coordinates are excluded), and outputs are substituted by function of parametrization variables. This transformation assures that the robot endpoint follows the desired trajectory. Zero dynamics (system flexibility) oscillations are realized only in the perturbation of parametrization variables. The invariant control, based on the harmonic functions [3], is applied to the model of a flexible robotic arm, and its parameters are obtained using predictive optimization.

Many authors have already investigated the control of underactuated mechanical systems. Neurodynamic programming, as a main part of the control in [4], is used for the control of the ball and beam system. In [5], a method of adding the actuators to the underactuated joints and running the swing-up motion with optimization minimizing the action of the added actuators is used. Sliding mode control extended with fuzzy logic is introduced in [6] to control the motion of overhead crane. An approach with optimal control is used in [7]. The adaptive approach used for tracking control is presented in [8]. Partly stable controllers derived using the dynamic model of the manipulator create basis for the optimal switching sequence control in [9]. An energy-based approach used in the design of stabilizing controllers is presented in [10]. Friction forces are used in the energy control strategy in [11]. Authors in [12] introduce an averaging theory for controlling the systems with and without drift terms. In [13], the trajectory tracking controller is used for following the upward trajectory created from the downward simulation (similar to the natural motion in [14] and the natural oscillations in [15] of a parallel robot) and stabilizing trajectory. An active disturbance rejection control is designed in [16] for a pendubot system. The walking-like trajectory tracking of an acrobot is continuously developed in [17]. A 3D bipedal gait walking robot is presented in [18], and the development of humanoid motion design is investigated in [19]. The robust sliding mode control is developed for stabilizing and trajectory tracking of a ballbot system in [20], and the wheel-legged biped robot motion is stabilized using active arm manipulator in [21]. The dynamic balancing of quadruped robot with optimal controller is presented in [22], and a different quadruped robot design is presented in [23]. Ref. [24] presents linear quadratic regulator applied to control vertical take-off and landing (VTOL) multi-rotor vehicle with passive rotor tilting mechanism. A control based on a function approximation technique using the variation of an approximated square system is presented in [25,26]. The cyclic control using sinusoids for nonholonomic systems is presented in [27] (stabilization), and motion planning is presented in [28]. The control of a wheeled, inverted pendulum is transferred into dynamic trajectory planning with a boundary value problem in [29]. State feedback fuzzy system controller, based on $H_{\infty }$, is investigated in [30]. The underactuated mechanisms in the gravity field can be controlled in their equilibrium positions (see [31]), but there are always singularity points between these equilibrium states. To move the mechanism through this region (see [32]) is a challenging task, and the work in [14] shows one possible way to overcome it.

The investigation of system flexibility as a system underactuation is usually investigated with linearized stiffness; the authors of [33] deal with dynamic coupling of such underactuated manipulators. Modelling flexibility in underactuated systems by a rotational spring is presented in [34], and a similar flexibility model is investigated in [35], where the spring stiffness is varied, and a point stabilization control is developed. A model with linear springs is investigated in [36], where input–output linearization is performed along with a PID controller. Sliding mode control is used to control the flexible robotic arm in [37]. The current work aims to introduce the invariant control strategy to the structures with compliancy (chain of flexible robotic arms). The nonlinear system with flexible modes is obtained using the procedure presented in [38].

An invariant control is presented in [2]. The application of such control to the underactuated system (double pendulum, double pendulum on a cart) is shown in [14,39]. The current work does not include the flexible structures and does not deal with trajectory tracking. This paper aims to introduce the control strategy to the structures with compliancy—in our case, the flexible robotic arm. Another contribution of the paper is in its presentation of a control strategy which forces the undesired movements of the system into the desired system endpoint trajectory.

The paper is organized as follows. The exact input–output linearization is presented in Section 2 for underactuated systems. Next, the coordinates are separated into actuated and unactuated. Then, consecutive output parameterization is introduced. Section 2.1 introduces invariant control with harmonic functions and the natural motion used for the Eigen frequency evaluation. Section 3 introduces the presented method to flexible robotic arm: modelling and parameterization with input–output linearization is applied to the flexible robot. Next, the particular control application to the flexible robotic arm and simulation results are presented. Section 4 concludes the work.

2. Materials and Methods

The control of underactuated, flexible structure using invariant control is presented.

2.1. Exact Input–Output Linearization with Parametrization

Let us consider a nonlinear dynamic system (Eigen equations of motion):

$$\mathbf{M}\left(\mathbf{q}\right)\ddot{\mathbf{q}}=\mathbf{F}\left(\mathbf{q},\dot{\mathbf{q}}\right)+\mathbf{B}\left(\mathbf{q}\right)\mathbf{u}$$

(1)

where the matrix $\mathbf{M}[n\times n]$ is an inertia matrix; $\mathbf{q}[n\times 1]$ are coordinates describing the system; $\dot{\mathbf{q}}$ and $\ddot{\mathbf{q}}$ their time derivatives; $\mathbf{F}[n\times 1]$ is a matrix containing forces dependent on velocities and external forces besides control inputs; $\mathbf{B}[n\times k]$ is the control input distribution matrix; the vector $\mathbf{u}[k\times 1],k<n$ contains the control inputs.

The state vector $\mathbf{q}$ could be rewritten without loss of generality in the following way:

$$\mathbf{q}=\left[\begin{array}{c}{\mathbf{q}}_a\\ {\mathbf{q}}_z\end{array}\right]$$

(2)

where ${\mathbf{q}}_a$ has the same dimension as the rank of matrix $\mathbf{B}$. Remaining z coordinates in vector ${\mathbf{q}}_z$ represent not directly actuated states (zero dynamics). This coordinate selection is generally not unique. The system (1) could be rewritten as follows:

$$\begin{array}{c}\hfill {\mathbf{M}}_{aa}^I\left(\mathbf{q}\right){\ddot{\mathbf{q}}}_a+{\mathbf{M}}_{az}^I\left(\mathbf{q}\right){\ddot{\mathbf{q}}}_z={\mathbf{F}}_a^I\left(\mathbf{q},\dot{\mathbf{q}}\right)+{\mathbf{B}}_a\left(\mathbf{q}\right)\mathbf{u}\end{array}$$

(3)

$$\begin{array}{c}\hfill {\mathbf{M}}_{za}^I\left(\mathbf{q}\right){\ddot{\mathbf{q}}}_a+{\mathbf{M}}_{zz}^I\left(\mathbf{q}\right){\ddot{\mathbf{q}}}_z={\mathbf{F}}_z^I\left(\mathbf{q},\dot{\mathbf{q}}\right)+{\mathbf{B}}_z^I\left(\mathbf{q}\right)\mathbf{u}\end{array}$$

(4)

where the matrix ${\mathbf{B}}_a[a\times a]$ is the largest regular sub-matrix of the matrix $\mathbf{B}$ with the rank a (the actuator redundancy is not covered, so $a=k$). The number of ${\mathbf{q}}_z$ coordinates is $z=n-a$. The matrix dimensions in (3) and (4) are dependent on the lower indexes, and matrix ${\mathbf{B}}_z^I$ has dimensions $[z\times a]$. Let us consider the vector of system outputs $\mathbf{y}$ of the same dimension as the vector ${\mathbf{q}}_a$. There is a relation between system outputs and system coordinates:

$$\mathbf{y}=\mathbf{f}\left({\mathbf{q}}_a,{\mathbf{q}}_z\right)$$

(5)

with a vector $\mathbf{f}[a\times 1]$ representing the dependence of the output on the system coordinates ${\mathbf{q}}_a$ and ${\mathbf{q}}_z$. There exist the inverse relation ${\mathbf{f}}^{-1}$, such that the coordinates ${\mathbf{q}}_a$ and its second time derivative can be uniquely determined:

$$\begin{array}{cc}& {\mathbf{q}}_a={\mathbf{f}}^{-1}\left(\mathbf{y},{\mathbf{q}}_z\right)\end{array}$$

(6)

$$\begin{array}{cc}& {\ddot{\mathbf{q}}}_a={\mathbf{J}}_y\left(\mathbf{y},{\mathbf{q}}_z\right)\ddot{\mathbf{y}}+{\mathbf{J}}_z\left(\mathbf{y},{\mathbf{q}}_z\right){\ddot{\mathbf{q}}}_z+\mathbf{S}\left(\mathbf{y},\dot{\mathbf{y}},{\mathbf{q}}_z,{\dot{\mathbf{q}}}_z\right)\end{array}$$

(7)

where ${\mathbf{J}}_y[a\times a]$, ${\mathbf{J}}_z[a\times z]$ and $\mathbf{S}[a\times 1]$ are matrices resulting the time derivatives. The outputs $\mathbf{y}$ are often parametrized, especially in robotics applications. The parametrized trajectory is defined:

$$\mathbf{P}=\mathbf{y}=\mathbf{P}\left(\mathbf{p}\left(t\right)\right)$$

(8)

with the following output acceleration:

$$\ddot{\mathbf{y}}={\mathbf{J}}_p\left(\mathbf{p}\right)\ddot{\mathbf{p}}+\mathbf{T}\left(\mathbf{p},\dot{\mathbf{p}}\right)$$

(9)

where vector $\mathbf{P}[a\times 1]$ is output parametrization and ${\mathbf{J}}_p[a\times p]$, $\mathbf{T}[a\times 1]$ are time derivative resulting matrices with p equalling the number of parameters. The inputs $\mathbf{u}$ are expressed from (3) and substituted into (4), and then the second time derivative of ${\mathbf{q}}_a$, much like $\mathbf{y}$ from (7) and (9), are substituted, and the final set of equations is obtained:

$$\begin{array}{ccc}\hfill {\mathbf{M}}_{aa}\ddot{\mathbf{p}}+{\mathbf{M}}_{az}{\ddot{\mathbf{q}}}_z& =& {\mathbf{F}}_a+{\mathbf{B}}_a\mathbf{u}\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {\mathbf{M}}_{za}\ddot{\mathbf{p}}+{\mathbf{M}}_{zz}{\ddot{\mathbf{q}}}_z& =& {\mathbf{F}}_z\hfill \end{array}$$

(11)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{M}}_{aa}=& {\mathbf{M}}_{aa}^I{\mathbf{J}}_y{\mathbf{J}}_p\hfill \\ \hfill {\mathbf{M}}_{az}=& {\mathbf{M}}_{az}^I+{\mathbf{M}}_{aa}^I{\mathbf{J}}_z\hfill \\ \hfill {\mathbf{M}}_{za}=& \left({\mathbf{M}}_{za}^I-{\mathbf{B}}_z^I{\mathbf{B}}_a^{-1}{\mathbf{M}}_{aa}^I\right){\mathbf{J}}_y{\mathbf{J}}_p\hfill \\ \hfill {\mathbf{M}}_{zz}=& {\mathbf{M}}_{zz}^I+{\mathbf{M}}_{za}^I{\mathbf{J}}_z-{\mathbf{B}}_z^I{\mathbf{B}}_a^{-1}\left({\mathbf{M}}_{az}^I+{\mathbf{M}}_{aa}^I{\mathbf{J}}_z\right)\hfill \\ \hfill {\mathbf{F}}_a=& {\mathbf{F}}_a^I-{\mathbf{M}}_{aa}^I\left(\mathbf{S}+{\mathbf{J}}_y\mathbf{T}\right)\hfill \\ \hfill {\mathbf{F}}_z=& {\mathbf{F}}_z^I-{\mathbf{M}}_{za}^I\left(\mathbf{S}+{\mathbf{J}}_y\mathbf{T}\right)-\hfill \\ & -{\mathbf{B}}_z^I{\mathbf{B}}_a^{-1}\left({\mathbf{F}}_a^I-{\mathbf{M}}_{aa}^I\left(\mathbf{S}+{\mathbf{J}}_y\mathbf{T}\right)\right).\hfill \end{array}\end{array}$$

The new input variable $\mathbf{w}$ is equal to the second time derivative of parameters vector $\mathbf{p}$ in order to provide input–output linearization. Then, the system inputs $\mathbf{u}$ are derived from the (10) and (11):

$$\begin{array}{cc}\hfill \mathbf{u}=& {\mathbf{B}}_a^{-1}\left({\mathbf{M}}_{aa}-{\mathbf{M}}_{az}{\mathbf{M}}_{zz}^{-1}{\mathbf{M}}_{za}\right)\ddot{\mathbf{p}}+\hfill \\ & +{\mathbf{B}}_a^{-1}\left({\mathbf{M}}_{az}{\mathbf{M}}_{zz}^{-1}{\mathbf{F}}_z-{\mathbf{F}}_a\right)\hfill \end{array}$$

(12)

with an assumption of the invertibility of ${\mathbf{M}}_{zz}$, which is always fulfilled for flexible mechanical systems.

2.2. Invariant Control

In this section, the control variable $\mathbf{w}$ is assembled. There is a dynamic system to be controlled (derived from the previous section):

$$\begin{array}{ccc}\hfill \ddot{\mathbf{p}}& =& \mathbf{w}\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\ddot{\mathbf{q}}}_z& =& {\mathbf{M}}_{zz}^{-1}\left({\mathbf{F}}_z-{\mathbf{M}}_{za}\mathbf{w}\right).\hfill \end{array}$$

(14)

The input variable $\mathbf{w}$ directly influences parameters of $\mathbf{p}$. However, how can the zero dynamics coordinates ${\mathbf{q}}_z$ be controlled? The proposed way is to split the input variable into arbitrary part ${\mathbf{w}}_A$, which controls the parameters $\mathbf{p}$, and invariant part ${\mathbf{w}}_I$, which controls coordinates ${\mathbf{q}}_z$, with specific influence on parameters $\mathbf{p}$.

$$\mathbf{w}={\mathbf{w}}_A+{\mathbf{w}}_I$$

(15)

The arbitrary part is assembled to reach the final configuration at time, T. The invariant part must satisfy invariancy conditions (16). Hence it modifies parameters $\mathbf{p}$ and their time derivatives only inside the interval $\left(0,T\right)$. The invariancy conditions are the following:

$$\begin{array}{c}\hfill {\int }_0^T{\mathbf{w}}_I\mathrm{d}t=\mathbf{0}\\ \hfill {\int }_0^T{\int }_0^t{\mathbf{w}}_I\mathrm{d}t\mathrm{d}t=\mathbf{0}.\end{array}$$

(16)

If the input variable (15) is substituted into the (13), and the system of equations is integrated from $t=0$ to $t=T$, then the following is obtained:

$$\begin{array}{cc}\hfill \dot{\mathbf{p}}\left(T\right)-\dot{\mathbf{p}}\left(0\right)=& {\int }_0^T{\mathbf{w}}_A\mathrm{d}t+{\int }_0^T{\mathbf{w}}_I\mathrm{d}t={\int }_0^T{\mathbf{w}}_A\mathrm{d}t\hfill \\ \hfill \mathbf{p}\left(T\right)-\mathbf{p}\left(0\right)=& {\int }_0^T{\int }_0^t{\mathbf{w}}_A\mathrm{d}t\mathrm{d}t+{\int }_0^T{\int }_0^t{\mathbf{w}}_I\mathrm{d}t\mathrm{d}t=\hfill \\ \hfill =& {\int }_0^T{\int }_0^t{\mathbf{w}}_A\mathrm{d}t\mathrm{d}t\hfill \end{array}$$

(17)

From (17), it is obvious that parameters $\mathbf{p}$ and their time derivatives at the end (final) time depend only on the arbitrary part of the control (15). The arbitrary part ${\mathbf{w}}_A$ could be arbitrarily chosen, just to fulfil the end position and velocities of the parameters $\mathbf{p}$. After the selection of the arbitrary part of the control, it is a time-dependent part, and (14) changes into (18).

$${\ddot{\mathbf{q}}}_z={\mathbf{M}}_{zz}^{-1}\left({\mathbf{F}}_z-{\mathbf{M}}_{za}\left({\mathbf{w}}_A\left(t\right)+{\mathbf{w}}_I\right)\right)$$

(18)

The invariant part of the control must fulfil the conditions of (16), and the control coordinates ${\mathbf{q}}_z$ from (18) must reach the desired states at the end time, T. The suitable functions are the following harmonic functions:

$${\mathbf{w}}_I=\sum _i{\mathbf{A}}_{i}sin\left(2\pi f_{i}t+{\varphi }_i\right)$$

(19)

The parameters ${\mathbf{A}}_i$, ${\varphi }_i$ of harmonic functions (19) are obtained from the optimization process, and the frequencies, $f_i$, are based on the so-called natural motion.

2.3. Natural Motion

According to [14], natural motion is derived from the vicinity of an unstable equilibrium to the vicinity of a stable equilibrium. The motion takes place in the gravity field with dissipative members instead of actuators. It is helpful to rearrange the formulation for the flexible systems. For obtaining the natural motion, it is used the simulation with zero invariant part ${\mathbf{w}}_I$, controlled only with the arbitrary part of the control ${\mathbf{w}}_A$. The system is moved from the initial position to the end position, where the system stays, and residual (zero dynamics) vibrations are exerted.

All system states are recorded during the motion, and the frequencies, $f_i$, are obtained from the frequency analysis of the recorded behaviour.

3. Results

The present section introduces the presented procedure with invariant control to the flexible robotic structure.

3.1. Model of Flexible Robot

A flexible robot is modelled as a double pendulum with flexible arms. It consists of two flexible beams, interconnected by a revolute joint, and coupled to the base frame with the second revolute joint. Every beam has flexibility modelled by three modes. There are actuators in each revolute joint. The model covers the mass of the actuator, and the endpoint is loaded by extra mass. The system is situated in the gravity field, see Figure 1.

Eigen equations of motion of the flexible double pendulum are obtained using the procedure derived in [38]. The implementation is presented in [40]. The Eigen equations of motion are in the following form:

$$\mathbf{M}\left(\mathbf{x}\right)\ddot{\mathbf{x}}=\mathbf{F}\left(\mathbf{x},\dot{\mathbf{x}}\right)+\mathbf{B}\left(\mathbf{x}\right)\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]$$

(20)

where

$$\mathbf{x}=\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\delta }_{11}\\ {\delta }_{12}\\ {\delta }_{13}\\ {\delta }_{21}\\ {\delta }_{22}\\ {\delta }_{23}\end{array}\right]\mathbf{B}=\left[\begin{array}{cc}1& -1\\ 0& 1\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]$$

Furthermore, matrix $\mathbf{M}$ is an inertia matrix, vector $\mathbf{F}$ covers velocity-dependent forces, gravity forces, and damping forces, and $u_1$, $u_2$ are actuator moments. Beam k endpoint position is determined by angular rotation ${\phi }_k$ and beam deflection $d_k$. Beam deflection is obtained according to the following equation:

$$d_k=\sum _{i=1}^3{d}_{ki}=\sum _{i=1}^3{\mathrm{\Phi }}_{ki}{\delta }_{ki}$$

(21)

where ${\delta }_{ki}$ are dimensionless modal coordinates and ${\mathrm{\Phi }}_{ki}$ are coefficients of modal parameters beam deflection at particular point. The angular deflection of the first beam endpoint is needed for the determination of the second beam position. This endpoint angular deflection is indicated as $\theta$, and generally

$${\theta }_k=\sum _{i=1}^3{\mathrm{\Psi }}_{ki}{\delta }_{ki}$$

(22)

where ${\mathrm{\Psi }}_{ki}$ are coefficients of modal parameters, beam angular deflection at particular point. Angle ${\theta }_1$ is the relative angle between the un-deformed shape of the first beam and its endpoint angular deflection (or endpoint tangent to deformed beam shape).

3.2. Input–Output Linearization with Output Parametrization

It is necessary to split the state vector $\mathbf{x}$ into ${\mathbf{x}}_a$ and ${\mathbf{x}}_z$ according to (2). The claim on the regularity of submatrix ${\mathbf{B}}_a$ derived from matrix $\mathbf{B}$ in (20) directly implies the splitting of vector $\mathbf{x}$:

$${\mathbf{x}}_{\mathbf{a}}=\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right]{\mathbf{x}}_{\mathbf{z}}=\left[\begin{array}{c}{\delta }_{11}\\ {\delta }_{12}\\ {\delta }_{13}\\ {\delta }_{21}\\ {\delta }_{22}\\ {\delta }_{23}\end{array}\right]$$

(23)

The flexible double pendulum endpoint coordinates $x_2$ and $y_2$ are system outputs $\mathbf{y}$, and it corresponds with the length of vector ${\mathbf{x}}_a$. Input–output linearization can be accomplished. Outputs are obtained using the following relations:

$$\mathbf{y}=\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]={\mathbf{R}}_2\left(\mathbf{x},{\mathrm{\Phi }}_{11},\dots ,{\mathrm{\Phi }}_{23},{\mathrm{\Psi }}_{11},\dots ,{\mathrm{\Psi }}_{23}\right)$$

(24)

The coordinates in the vector ${\mathbf{x}}_a$ are now changed to the $\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]$ end effector (double pendulum endpoint) coordinates. Equation (24) is used to obtain $\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right]$ coordinates and its time derivatives according to (6) and (7) to be substituted into (20). Following one dimensional parametrization (typically endpoint trajectory, see (25)) ensures the endpoint position and leaves free movement along this parameter. Then, the zero dynamics control actions use this free parameter for realizing its unwanted movements.

$$\mathbf{y}=\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]={\mathbf{P}}_2\left(p\left(t\right)\right)$$

(25)

The system inputs $\mathbf{u}=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]$ are obtained from (12), but the following procedure could also be used. Outputs functions (24) and their parametrization (25) are the second-time derivate:

$$\begin{array}{ccc}\hfill {\ddot{\mathbf{R}}}_2& =& \left[\begin{array}{c}{\ddot{x}}_2\\ {\ddot{y}}_2\end{array}\right]=\frac{\partial {\mathbf{R}}_2}{\partial \mathbf{x}}\ddot{\mathbf{x}}+\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial {\mathbf{R}}_2}{\partial \mathbf{x}}\right)\dot{\mathbf{x}}\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\ddot{\mathbf{P}}}_2& =& \left[\begin{array}{c}{\ddot{x}}_2\\ {\ddot{y}}_2\end{array}\right]=\frac{\partial {\mathbf{P}}_2}{\partial p}\ddot{p}+\frac{{\partial }^2{\mathbf{P}}_2}{\partial p^2}{\dot{p}}^2\hfill \end{array}$$

(27)

In (26), the coordinates acceleration $\ddot{\mathbf{x}}$ are substituted from (20), and outputs acceleration ${\ddot{x}}_2$, ${\ddot{y}}_2$ are replaced from (27). The control inputs are obtained from the resulting alternative expression:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]=& {\left(\frac{\partial {\mathbf{R}}_2}{\partial \mathbf{x}}{\mathbf{M}}^{-1}\mathbf{B}\right)}^{-1}(\frac{\partial {\mathbf{P}}_2}{\partial p}\ddot{p}+\frac{{\partial }^2{\mathbf{P}}_2}{\partial p^2}{\dot{p}}^2-\hfill \\ & -\frac{\partial {\mathbf{R}}_2}{\partial \mathbf{x}}{\mathbf{M}}^{-1}\mathbf{F}-\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial {\mathbf{R}}_2}{\partial \mathbf{x}}\right)\dot{\mathbf{x}})\hfill \end{array}$$

(28)

The right-hand side of (28) also depends on the state vector $\mathbf{x}$ and its time derivative $\dot{\mathbf{x}}$. Therefore, inputs could also be evaluated according to (12) in the form with only p and ${\mathbf{x}}_z$ coordinates.

The advantage of such a strategy is that flexible robot endpoint never leaves the prescribed trajectory, and zero dynamics (${\mathbf{x}}_z$) control actions are realized only along the parameter, p.

3.3. Control of Flexible Robot

The task is to move the flexible robot endpoint along the circle arc, see Figure 2. The circle centre is in the joint between the base frame and the first arm. The radius is r in (29), and the circle parameter p is angular position measured from the robot’s initial position.

$$r=\sqrt{L_1^2+L_2^2+2L_1{L}_{2}cos\frac{\pi }{4}}$$

(29)

At first, the arbitrary part of the control (${\mathbf{w}}_A$) is prepared. According to (13), the input is the second time derivative of parameter p. Design of arbitrary part of the control is based on the parameter p behaviour. The initial position and end position have velocities and accelerations selected to be zero. Such conditions satisfy a fifth-order polynomic in function p, see Figure 3. Movement duration is one of the free parameters, but it will be determined by the desired speed of the robot endpoint.

3.3.1. Natural Motion Frequencies

Natural motion is performed only with arbitrary part of the control ${\mathbf{w}}_A$ with the duration of 1 s. When the robot reaches the end position, residual vibrations are measured, see Figure 4, left. It is evident that all system coordinates vibrate at just one natural frequency ($f_N=4.2915$ Hz), and this is verified by the frequency analysis in Figure 4, right. For the different simulation times, the resulting natural motion frequency is the same.

3.3.2. Invariant Control Parameters

The invariant part of the control is constructed according to (19). The optimization process is performed to obtain the number of intervals ${f_N}^{-1}$, natural frequency multiplication factors $k_i$ for obtaining particular $f_i$, phase ${\varphi }_i$, and corresponding amplitudes ${\mathbf{A}}_i$. The phase ${\varphi }_i=-0.5\pi$ is same for all frequencies $f_i$.

The objective function $OF$ aims to fulfil the prescribed trajectory endpoint position and velocity. Physical coordinates ${\phi }_i$ and modal coordinates ${\delta }_{jk}$ at end time should be ideally equal to the trajectory endpoint static position and their time derivatives (velocities) should be zero, see (30).

$$\begin{array}{cc}\hfill OF=& a_1\sum _i\left|{\phi }_i\left(t_{END}\right)-{\phi }_{iEND}\right|+\hfill \\ & +a_2\sum _j\sum _k\left|{\delta }_{jk}\left(t_{END}\right)-{\delta }_{jkEND}\right|+\hfill \\ & +a_3\sum _i\left|{\dot{\phi }}_i\left(t_{END}\right)\right|+a_4\sum _j\sum _k\left|{\dot{\delta }}_{jk}\left(t_{END}\right)\right|\hfill \end{array}$$

(30)

The optimization process minimizes the objective function in (30), where $a_1,\dots ,a_4$ are variable weights, $t_{END}$ is prescribed movement end time, ${\phi }_{iEND}$ and ${\delta }_{jkEND}$ are prescribed movement endpoint coordinates.

The extensive optimization brings three frequency multiplication factors $k_1=1$, $k_2=0.5$ and $k_3=2.5$, with resulting frequency $f_i=k_i{f}_N$, the 16 intervals of ${f_N}^{-1}$, and amplitudes ${\mathbf{A}}_i$ shown in

Table 1. Amplitudes.

	${\mathbf{A}}_1$	${\mathbf{A}}_2$	${\mathbf{A}}_3$
1	−0.0256	0.0504	0.0050
2	−0.0251	−0.0502	0.0051
3	0.0253	0.0500	0.0050
4	0.0255	−0.0506	−0.0050
5	−0.0253	0.0523	−0.0050
6	−0.0264	−0.0494	−0.0050
7	0.0251	0.0514	0.0051
8	0.0263	−0.0407	0.0051
9	−0.0251		0.0050
10	−0.0252		−0.0051
11	0.0253		−0.0051
12	0.0251		−0.0050
13	−0.0253		0.0050
14	−0.0252		0.0050
15	0.0244		0.0050
16	0.0225		−0.0050
17			−0.0050
18			−0.0051
19			0.0051
20			0.0051
21			0.0052
22			−0.0050
23			−0.0050
24			−0.0051

.

3.4. Simulation

The applied control forces the robot endpoint to move along the desired trajectory. The difference between the control methods (computed torques method and invariant control) is in the behaviour of the zero dynamics movements from the end time of motion. The computed torque method cannot effectively cancel out zero dynamics oscillations after reaching the final position. The invariant control implies slightly different claim on the system torques $u_1$ and $u_2$ than the computed torque method, but after reaching the final time, the system is still. The influence of invariant control to the desired parameter p behaviour is on Figure 5. The amplitude of the parameter, p, invariant control modification is significantly low in comparison with the desired motion change.

Figure 6 and Figure 7 show the behaviour of $u_1$ and $u_2$ input torques, respectively. The $u_1$ plot of Figure 6 shows the first input torque needed for satisfying the desired motion. There are two curves, one is dedicated for the computed torques method, and the second one corresponds to the torque claim of the invariant control. In the cut is the detailed behaviour after the prescribed motion final time. Invariant control modifies input torque so that there is no residual vibration of zero dynamics, contrary to the computed torque method. The second part of Figure 6 graph $\Delta u_1$ shows the difference between computed torques and invariant control $u_1$ behaviour. The red curve represents modification of computed torques $u_1$, claimed due to the invariant part of the control, and the blue line shows the residual (zero dynamics) vibrations of the system with control using the computed torques method. The corresponding behaviour of $u_2$ input torque is shown on Figure 7.

The advantage of invariant control is that the system moves from one equilibrium to another equilibrium position. The input torque is modified within $5\%$ of the reference torque (computed torque method serves as the reference). The movement within parameter, p, is the only disadvantage of the invariant control, but it is of small amplitude, as shown on Figure 5 (bottom).

4. Conclusions

The control of a flexible robot as a representative of underactuated system with zero dynamics is performed using invariant control method. The basis of invariant control are harmonic functions, which are used for feedforward control of particular flexible system. The harmonic functions provide smooth motion as is needed for vibration cancellation when reaching the desired position. Presented here is an invariant control of flexible underactuated structure in connection with observing the desired trajectory.

The control strategy uses motion frequencies inspired by so-called natural motion. The residual oscillations provide the natural motion frequency, which is used to cancel out unwanted vibrations. The control input is separated into two parts. The arbitrary part of the control input (${\mathbf{w}}_A$) is designed to control the directly actuated part of the dynamic system. The invariant part of the control (${\mathbf{w}}_I$) is selected to steer the system zero dynamics in the desired way. The invariant part of the control fulfils the invariant conditions with respect to overall control input behaviour. The amplitudes of chosen periodic functions are obtained by optimization process with the final position in objective function. The achieved accuracy depends on the optimization process inputs, the trajectory length, and the number of used frequencies. The control parameters are predicted using the optimization process and resulting tracking accuracy is the same for the same trajectory.

The main aim of this work was to present the functionality of an invariant control approach to a system with flexible elements. Further research will focus on the extension of the invariant control by feedback to avoid the optimization process. The proposed concept will be enhanced by the model predictive, wave-based or sliding mode control to secure missing feedback. To fulfil the invariance conditions (16) the control strategy will be formulated in discrete form.. The presented control can be practically applied in serial robotics, i.e., robotic manipulators or machine tools.