1 Introduction

1.1 Underactuated system

The dynamic system described by a set of second-order ODE in form

$$\begin{aligned} \ddot{\varvec{q}}(t)=\varvec{f}_{1}+\varvec{f}_{2}\varvec{u}(t) \end{aligned}$$
(1)

with generalized coordinates \(\varvec{q}(t)=\left[ q_{1},q_{2},\ldots ,q_{n}\right] \) and inputs \(\varvec{u}(t)=\left[ u_{1},u_{2},\ldots ,u_{w}\right] \) is called underactuated if the unbounded control inputs cannot produce accelerations \(\ddot{\varvec{q}}\) in an arbitrary direction. This could be verified by the condition of \(\mathrm {rank}\varvec{f}_{2}<\mathrm {dim}\varvec{q}\). The situation when there is a fewer number of inputs than the number of degrees of freedom is called trivial underactuation. There are many systems with the underactuated property, e.g., acrobot, pendubot, cart-pole, the beam-and-ball, inertia-wheel pendulum, airplane and multicopter, hovercraft, surface vessel. A recent review of underactuated systems and their control is presented widely in [7].

1.2 Input coupling

One of the insufficiently studied problems related to underactuated systems is the input coupling. The system described by Eq. 1 has input coupling if at least one input acts on at least two accelerations [10]. This situation causes big problems in the trajectory tracking task—generalized coordinates cannot be controlled separately.

1.3 Trajectory tracking task

The most popular method of tracking control of underactuated systems with the input coupling effect is based on the change of variables that converts this problem into noncoupled [4]. Then, some accelerations are separately controlled by inputs and some stay without control (selective control method). This contribution focuses on a new method of full-state control, based on the pseudoinverse operation combined with a computed torque technique and PD (proportional-derivative) feedback presented in [5]. The inputs are then proposed as

$$\begin{aligned} \varvec{\varvec{\tau }}(t)=\varvec{f}_{2}^{+}\left( \ddot{\varvec{d}}(t)-\varvec{f}_{1}+\varvec{K}_{D}\dot{\varvec{e}}(t)+\varvec{K}_{P}\varvec{e}(t)\right) \end{aligned}$$
(2)

where \(\varvec{f}_{2}^{+}\) is a pseudoinverse of \(\varvec{f}_{2}\), \(\varvec{d}(t)\) is a vector of desired trajectory functions, \(\varvec{K}_{D}\) and \(\varvec{K}_{P}\) are diagonal matrices of positive coefficients, and \(\varvec{e}(t)=\varvec{d}(t)-\varvec{q}(t)\) is a vector of tracking errors. Moore–Penrose pseudoinverse defined by Tikhonov’s regularization can be calculated with formula

$$\begin{aligned} \varvec{f}_{2}^{+}=\underset{\delta \rightarrow 0}{\lim }(\varvec{f}_{2}^{T}\varvec{f}_{2}-\delta \varvec{I})^{-1}\varvec{f}_{2}^{T} \end{aligned}$$
(3)

Substitution of the control method proposed by Eq. (2) instead of \(\varvec{u}(t)\) in system of equations of motion (1) yields the errors dynamic equation

$$\begin{aligned} \varvec{\ddot{e}}+\varvec{f_{2}}\varvec{f}_{\varvec{2}}^{\varvec{+}}\varvec{K_{D}\dot{e}}+\varvec{f_{2}f_{2}^{+}K_{P}e}=\left( \varvec{f_{2}f_{2}^{+}}-\varvec{I}\right) \left( \varvec{f_{1}}-\varvec{\ddot{d}}\right) \end{aligned}$$
(4)

Analytic and numerical analysis of Eq. (4) can ensure stability of the system while tracking desired trajectory. This equation is also useful during tuning process of PD gains.

1.4 Chaotic behaviors of underactuated systems

Researchers show that some physical underactuated systems behave in a chaotic manner. One of these systems is a tethered satellite system [11]. In such a system, the mother satellite is treated as a rigid body, subsatellite as a point of mass and tether as an inelastic massless beam. The whole system moves in Kepler elliptical orbit. These ones can be controlled to convert chaotic motions into periodic motions using time-delay autosynchronization method. A time-delay feedback control strategy can also be used to stabilize unstable periodic orbits of a two-link planar manipulator [9]. Free-joint manipulator with one actuated and one unactuated joint was successfully controlled by using periodic inputs to obtain desired trajectories [12]. Complex chaotic behaviors were also observed in planar five-bar closed-chain mechanism [8]. In previous research, various types of underactuated hovercrafts were analyzed in terms of nonlinear control for trajectory tracking [1,2,3] but without input coupling effect and bifurcation analysis.

This research focuses on chaotic behaviors of an underactuated system controlled with novel full-state control method with pseudoinverse operation, where irregular behaviors stand from the input’s limitation, not from dynamics and control.

2 Nonlinear behaviors of a hovercraft model

2.1 Model formulation

In this paper, one of the simplest underactuated models is described [6]. It can be used for a basic representation of a hovercraft, rocket or sliding vehicle. Consider a planar rigid body moving on a plane (Fig. 1). The object has mass m and inertia \(I_{ C}\) in the center of mass (point C). Coordinates x(t) and y(t) describe its position; \(\varphi (t)\) denotes the angle between the object symmetry line and the X axis of the global coordinate system \(O_{XY}\). The vector of force \(\overrightarrow{\varvec{F}}\) acts on the object in a point away from point C by distance a. Force \(\overrightarrow{\varvec{F}}\) as a system’s input is described by its magnitude f(t) and direction angle \(\beta (t)\). Constant drag coefficients c and \(c_{\varphi }\) are used. The equations of motion for the system are as follows

$$\begin{aligned} m\ddot{x}(t)+c\dot{x}(t)= & {} f(t)\cos \big (\varphi (t)+\beta (t)\big ) \end{aligned}$$
(5)
$$\begin{aligned} m\ddot{y}(t)+c\dot{y}(t)= & {} f(t)\sin \big (\varphi (t)+\beta (t)\big ) \end{aligned}$$
(6)
$$\begin{aligned} I_{c}\ddot{\varphi }(t)+c_{\varphi }\dot{\varphi }(t)= & {} af(t)\sin \beta (t) \end{aligned}$$
(7)
Fig. 1
figure 1

Planar rigid body moving on a plan subjected with eccentric force

Right-hand sides of Eqs. (57) make them coupled with respect to inputs. Description of the system in a form of Eq. (1) leads to

$$\begin{aligned} \varvec{q}=\left[ \begin{array}{c} x(t)\\ y(t)\\ \varphi (t) \end{array}\right] \text{, } \varvec{f_{1}}=\left[ \begin{array}{c} -\frac{c}{m}\dot{x}(t)\\ -\frac{c}{m}\dot{y}(t)\\ -\frac{c_{\varphi }}{I_{ C}}\dot{\varphi }(t) \end{array}\right] \text{, } \varvec{f_{2}}=\left[ \begin{array}{cc} -\frac{\sin \varphi (t)}{m} &{} \quad \frac{\cos \varphi (t)}{m}\\ \frac{\cos \varphi (t)}{m} &{}\quad \frac{\sin \varphi (t)}{m}\\ \frac{a}{I_{ C}} &{} \quad 0 \end{array}\right] \text{, } \varvec{u}=\left[ \begin{array}{c} f(t)\sin \beta (t)\\ f(t)\cos \beta (t) \end{array}\right] \end{aligned}$$
(8)

A task of tracking control presented here focuses on an eight curve in a parametric form

$$\begin{aligned} \varvec{d}(t)=\left[ \begin{array}{c} R\sin (2\theta t)\\ 2R\sin (\theta t)\\ \mathrm {Arg}(\cos {2\theta t}+\sqrt{-1}\cos {\theta t}) \end{array}\right] \end{aligned}$$
(9)

with size parameter R and velocity parameter \(\theta \). Third component of Eq. (9) refers to desired tangential orientation of the object with reference to the trajectory. Substitution of desired trajectory functions (9) and the system’s description matrices into algorithm proposed by Eq. (2) results with control functions. The most interesting behavior of tracking errors occurs in the situation of bounded inputs. Limitation of the force magnitude f(t) results in maximum possible velocity of the system. Limitation of the force direction \(\beta (t)\) may cause irregular behaviors.

2.2 Bifurcation diagrams

For the presented system parameters were arbitrarily chosen: \(m=2\,\mathrm {kg}\), \(I_{ C}=0.1\,\mathrm {kg\,m^{2}}\), \(a=0.4\,\mathrm {m}\), \(c=0.6\,\mathrm {N\,s\,m^{-1}}\) and \(c_{\varphi }=0.1\,\mathrm {N\, m\, s}\). For desired eight-shaped trajectory and proposed control functions, numerical simulations of the system’s tracking errors were precisely analyzed. Figure 2 presents bifurcation diagrams for x, y and \(\varphi \) tracking errors. Limitation parameter \(\beta _\mathrm{max}\) of the force direction \(\beta (t)\) was chosen as a bifurcation parameter. Only the most interesting region of parameters is presented in these pictures. For \(\beta _\mathrm{max}>19.07^{\circ }\) region of periodic motion is observed. Example steady-state trajectory, phase portraits and Poincare maps are presented in Fig. 3. After a bifurcation at \(\beta _\mathrm{max}=19.07^{\circ }\), period doubling effect is visible only on the x tracking error phase portrait and Poincare map. Trajectory period stays unchanged, but trajectory loses symmetry in relation to x axis. Second bifurcation occurs at \(\beta _\mathrm{max}=17.85^{\circ }\). For its lower values period doubling for all tracking errors occurs (Fig. 4b, c, d), which involves trajectory qualitative change (Fig. 4a). Another bifurcation happens at \(17.59^{\circ }\) and \(17.50^{\circ }\). Then, a region of irregular motion starts. Figure 5a presents the system’s behavior for \(\beta _\mathrm{max}=17.4^{\circ }\). One can observe some one-dimensional sets on the Poincare maps. The system’s behavior for \(\beta _\mathrm{max}=17.2^{\circ }\) (Fig. 5b) presents interesting shapes on the Poincare maps—probably fractals. The last two presented examples give us probability of chaotic motion, which needs proof by proper values of Lyapunov exponents and shapes of Fourier spectrum. Region of \(\beta _\mathrm{max}\in (16.85^{\circ },16.96^{\circ })\) presents a situation of periodic motion with period five times longer than for \(\beta _\mathrm{max}>19.07^{\circ }\) (Fig. 6a). Starting from \(\beta _\mathrm{max}=16.85^{\circ }\) bifurcations occur up to the region of irregular motion that ends with qualitative change of solution at \(\beta _\mathrm{max}=16.65^{\circ }\) (Fig. 6b). Further decrease of \(\beta _\mathrm{max}\) is not reasonable because of a lack of similarities between the desired (tracked) trajectory and real trajectory.

Fig. 2
figure 2

Bifurcation diagrams for tracking errors: a x error, b y error, c \(\varphi \) error

Fig. 3
figure 3

Object motion for \(\beta _\mathrm{max}=19.15^{\circ }\): a trajectory on the XY plane, b phase portrait and Poincare map for x error, c phase portrait and Poincare map for y error, d phase portrait and Poincare map for \(\varphi \) error

Fig. 4
figure 4

Object motion for \(\beta _\mathrm{max}=17.7^{\circ }\): a trajectory on the XY plane, b phase portrait and Poincare map for x error, c phase portrait and Poincare map for y error, d phase portrait and Poincare map for \(\varphi \) error

Fig. 5
figure 5

Poincare maps for x, y and \(\varphi \) errors: a for \(\beta _\mathrm{max}=17.4^{\circ }\) (500 periods of desired trajectory), b for \(\beta _\mathrm{max}=17.2\) (1500 periods of desired trajectory)

2.3 Lyapunov exponents

Values of the system’s Lyapunov exponents were obtained with numerical simulation using method presented in [13]. For \(\beta _\mathrm{max}>19.07^{\circ }\) all exponents are less than zero that proves stability of the periodic solution. For \(\beta _\mathrm{max}<19.07^{\circ }\) there are regions where the highest exponent is greater than zero and other which are less than zero—that means chaotic behavior in this region occurs. Collocation of the highest exponent with one bifurcation diagram presented in Fig. 7 shows exactly the regions of chaos, periodic orbits and bifurcation points. Around the bifurcation points, numerically calculated exponents are very close to zero, but an exact zero value cannot be obtained without more precision and longer evaluation time.

Fig. 6
figure 6

Trajectory on the XY plane: a for \(\beta _\mathrm{max}=16.9^{\circ }\), b for \(\beta _\mathrm{max}=16^{\circ }\)

Fig. 7
figure 7

Collocation of the highest Lyapunov exponent values with the bifurcation diagram of x error, with respect to the force’s angle limitation \(\beta _\mathrm{max}\)

2.4 Fourier spectrum

To satisfy all formal proofs of chaos occurrence, Fourier spectrum from errors’ time series was calculated (Fig. 8). Basic motion harmonic (0.125 Hz) related to the tracked path period and its poliharmonics are visible. Additional sub- and poliharmonics become visible when \(\beta _\mathrm{max}\) parameter passes bifurcation points (Fig. 8a). Exemplary Fourier spectra are presented in Fig. 8b—for \(\beta _\mathrm{max}=19.4^{\circ }\) many dominant harmonics are visible, for \(\beta _\mathrm{max}=19.1^{\circ }\) spectrum becomes more smooth, for \(\beta _\mathrm{max}=17.0^{\circ }\) harmonics are not visible. Chaotic behaviors for \(\beta _\mathrm{max}<17.5^{\circ }\) are related to smooth spectra without distinguishable dominant harmonics.

Fig. 8
figure 8

Fourier spectrum for the x error, calculated with FFT algorithm and presented in dB: a colormap at low frequencies w.r.t. \(\beta _\mathrm{max}\), b exemplary graphs for given \(\beta _\mathrm{max}\) values

3 Conclusions

In this contribution nonlinear behaviors of an underactuated hovercraft model were presented. Trajectory tracking method, based on the Moore–Penrose pseudoinverse operation combined with a computed torque technique and PD feedback, gives the possibility to control full-system configuration (hovercraft position and rotation). It was presented on an eight-shaped trajectory that the limitation of force direction can change system behavior—from periodic solution close to the desired trajectory, through a sequence of bifurcation, up to chaotic situations. These chaotic behaviors were verified by the Lyapunov exponent values and Fourier spectra of time series. There is a possibility to use the Lyapunov exponents for system’s behavior prognosis, but online operation could be complicated due to a long evaluation time.

Example presented in this paper shows that the limitation of driving force direction may cause not only huge tracking errors but also irregular behaviors of underactuated systems. Systems controlled to track periodic trajectories with relatively high speed, like tethered satellite system, are particularly vulnerable to these behaviors. Irregular behaviors of the presented hovercraft model stand from the limitations of input, not from dynamics and control.

Future research could be focused on general analysis of error dynamics Eq. (4) and 6DoF underactuated models of airplanes and rockets in control task. Influence of system damping properties onto chaotic behavior and its existence should be analyzed.