Nonlinear Model Predictive Control with Evolutionary Data-Driven Prediction Model and Particle Swarm Optimization Optimizer for an Overhead Crane

Abstract

This paper presents a new approach to the nonlinear model predictive control (NMPC) of an underactuated overhead crane system developed using a data-driven prediction model obtained utilizing the regularized genetic programming-based symbolic regression method. Grammar-guided genetic programming combined with regularized least squares was applied to identify a nonlinear autoregressive model with an exogenous input (NARX) prediction model of the crane dynamics from input–output data. The resulting prediction model was implemented in the NMPC scheme, using a particle swarm optimization (PSO) algorithm as a solver to find an optimal sequence of the control actions satisfying multi-objective performance requirements and input constraints. The feasibility and performance of the controller were experimentally verified using a laboratory crane actuated by AC motors and compared with a discrete-time feedback controller developed using the pole placement technique. A series of experiments proved the effectiveness of the controller in terms of robustness against operating condition variation and external disturbances.

1. Introduction

Crane systems have been attractive for decades in the modeling and control field, mainly due to their underactuated properties that present challenging modeling and control problems that are theoretically interesting and significant from a practical viewpoint. Since cranes are widely used in the industrial, construction, and logistic sectors for transporting heavy and oversized loads by using wire ropes, cables, or chains, an excessive transient and residual cable-suspended payload oscillation caused by actuated horizontal motion degrades positioning performance and may pose safety hazards [1]. This problem is extensively studied in the literature, and various open-loop and closed-loop control strategies are comprehensively reviewed in [2], while the more current state of the art focusing on the modeling and control of overhead cranes is provided in [3]. Recent works include input shaping [4,5], time-optimal [6], command shaping [7], and feedforward [8] anti-sway crane control strategies. New concepts in the field of closed-loop control have been developed using PID-like coupling control [9,10], sliding mode control [11] combined with a state observer [12] or PD controller [13], passivity-based control [14,15], H-infinity [16], and back-stepping-based [17] control schemes.

An increase in computing power has resulted in growing applications of linear MPC and nonlinear MPC (NMPC) in a variety of processes. MPC has received attention in relation to crane control systems because it allows the implementation of a set of constraints in the optimization algorithm that keep the crane in a feasible and safe operating region. A combination of feedforward control and NMPC was used in numerical simulations for the payload sway suppression of a planar nonlinear model of a shipboard crane [18]. A real-time MPC method with a linearized model of a mobile harbor boom crane that is hydraulically actuated was developed in [19] and tested on a LIEBHERR harbor mobile crane with a load capacity of 104 tons and a rope length of up to 70 m. This work is further developed in [20], where the nonlinear behavior of the system is simplified by applying exact linearization and input–output linearization. Tysse et al. [21] proposed combining a Lyapunov-based controller to stabilize the pendulum dynamics and an NMPC scheme used to control the position of a small-scale knuckle boom crane. An NMPC-based scheme was designed in [22] to optimize the pendulum motion driven by a mechatronic cart–winch system. In [23], the dynamics of a 3D overhead crane system were decoupled into two second-order systems, and two separate NMPC schemes were developed for simultaneously positioning a payload along the two motion planes of a crane’s bridge and trolley, respectively. Experiments carried out on a small laboratory-scale overhead crane confirmed the effectiveness of this strategy. In [24], the MPC approach is applied to minimize an objective function that is formulated as the integration of energy consumption and swing angle. In [25], the authors developed a path-following controller for a tower crane using nonlinear model predictive control with a fast suboptimal solution and verified the algorithm in simulations and on a small experimental stand. In [26], the authors proposed a model predictive controller for an overhead crane. The model was linearized around the equilibrium point, and the sway constraints were converted to input constraints using the kinematic equation of the crane system. The MPC strategy is employed in [27] to optimize the crane velocity profile when the maximum position overshoot is given by the crane operator.

The aforementioned studies are mostly based on physics-based models of underactuated crane dynamics. However, the accurate analytical modeling of complex coupled crane dynamics can be difficult and time-consuming; thus, some recent works have proposed data-driven modeling approaches. An MPC strategy involving a discrete-time linear model, recursive least-squares parameter estimation, and a PSO solver of a multi-objective control optimization problem was proposed and tested on a laboratory overhead crane driven by DC motors in [28]. This linear MPC scheme is based on a straightforward approach that assumes a parametric model of crane dynamics, thereby reducing the identification task to parameter estimation using the RLS algorithm with parameter projection. The parametric identification procedure assumes there is an adequate model structure representation with uncertain parameters that are subject to optimization. Parametric identification requires fewer data to identify a model; however, some a priori knowledge and engineering experience with respect to the analyzed system are required to choose a suitable model structure [29]. Moreover, there is no unique model structure corresponding to different crane applications under different operating conditions. In this case, non-parametric methods and machine learning techniques are alternative approaches that can be useful for determining a model structure for a specific process and data set. Jakovlev et al. [30] analyzed several multi-layer perceptron neural network (NN) structures and learning approaches (e.g., variable learning rate backpropagation, scaled conjugate gradients, recursive prediction error, and the Levenberg–Marquardt algorithm) to solve a time-series prediction problem in an MPC scheme for a quay crane. In [31], an overhead crane’s dynamics were approximated by an adaptive neuro-fuzzy inference system (ANFIS) trained with operational data using a hairpin RNA genetic algorithm. Kim et al. [32] used the NARX-NN data-driven model trained offline and online using extreme learning machines (ELMs) in the predictive control of a small-scale overhead crane with a prediction and control horizon of $N_p=N_c=15$. The Koopman operator theory and deep learning are merged in [33] to identify a model of a gantry crane subsequently used to develop a linear quadratic controller tested on a laboratory stand. A data-driven active learning control scheme based on the Koopman operator was studied in simulations on a bridge crane model [34]. Different structures of artificial neural network (ANN) models of an overhead crane’s inverse dynamics were trained and validated using the Levenberg–Marquardt algorithm, with data computed in simulations carried out using a model derived analytically in [35]. A hybrid PID and MPC scheme with Bayesian optimization for the data-driven identification of linear-model and controller parameter tuning was tested on a laboratory-scale overhead crane in [36]. An interesting alternative to model-based control is data-driven model-free control, which has recently proven its effectiveness for complex dynamical problems [37,38]. Examples of such studies for underactuated crane systems are presented in [39,40].

The selection of an appropriate model structure is a crucial point of identification that requires a tradeoff between model fitting and parsimony. Although the machine learning data-driven discovery of crane dynamics has been successfully explored in recent works, especially using an ANN framework, the problem of nonlinear model structure selection is rarely studied. To address this issue, a non-parametric data-driven identification of overhead crane dynamics using multi-gene genetic programming (MGGP) is proposed in [41]. However, symbolic regression can be prone to generating complex expressions and overfitting when aiming to decrease prediction error [42]. Thus, in [43], grammar-guided genetic programming (G3P) combined with l₀ sparse regression is applied to find a compromise between model complexity and prediction accuracy.

Following the work presented in [43] and in comparison to the linear parametric modeling-based MPC approach presented in [28], this paper develops an NMPC control strategy for an overhead crane based on a prediction model constructed using a regularized genetic programming (GP)-based approach. A polynomial NARX model structure is optimized using an evolutionary process and regularized least squares; the resulting prediction model is implemented in the NMPC scheme with a PSO algorithm used as the solver to find an optimal sequence of the control actions satisfying the multi-objective performance requirements and input constraints. The experimental verification of the proposed controller is carried out on a laboratory overhead crane by testing the NMPC performance at different operating points in the presence of external disturbances and comparing the results with those of a discrete-time feedback controller developed using the pole placement method. The main contributions of this paper are as follows:

(1)

A regularized genetic programming data-driven modeling approach is used to develop an NMPC scheme for an underactuated overhead crane system.

(2)

The proposed data-driven modeling approach and model-based control scheme are experimentally verified under different operating conditions and in terms of robustness with respect to external disturbances.

The rest of this paper is organized as follows: Section 2 introduces modeling assumptions, provides the identification procedure developed using GP with sparse regression, and describes the NMPC control strategy with the PSO solver. The hardware experimental setup and the results are presented in Section 3, while the final conclusions are made in Section 4.

2. Materials and Methods

2.1. Data-Driven Modeling of Crane Dynamics

The proposed data-driven modeling and NMPC approach are studied for a planar lumped-mass model of an overhead crane. A simplified planar single-pendulum system, often used as a benchmark case to develop and evaluate different anti-sway control strategies [1,2,3], is illustrated in Figure 1, where M, m, l, α, x, v, F are the masses of the trolley and pendulum, length of the pendulum, sway angle of the pendulum, trolley displacement and velocity, and force, respectively. The force consists of the actuating force (in real applications, usually produced by AC motors supplied by frequency inverters), disturbances, and friction forces. In this configuration, we assume that α, x, and v are measurable variables, the rope length l and the mass of the payload m are measurable or known parameters, and the actuating force is produced by an actuator controlled by a control signal u.

The crane is assumed to have high mechanical impedance, which allows for a decomposition of the crane dynamics into the actuated and unactuated parts and, therefore, can be modeled by two discrete-time nonlinear models, which can be presented in the NARX (nonlinear autoregressive with exogenous inputs) form with operating parameters $l$ and $m$ as

$$v\left(k\right)=f_1\left(u\left(k-1\right),\dots ,u\left(k-N_u\right),v\left(k-1\right),\dots ,v\left(k-N_v\right),l,m\right)+e_1\left(k\right)$$

(1)

$$\alpha \left(k\right)=f_2\left(v\left(k-1\right),\dots ,v\left(k-N_v\right),\alpha \left(k-1\right),\dots ,\alpha \left(k-N_{\alpha }\right),l,m\right)+e_2\left(k\right)$$

(2)

where $N_u,N_v,N_{\alpha }$ are the maximum lags for the system input and output, $e_1\left(k\right)$ and $e_2\left(k\right)$ are the errors at discrete time k, and the rope length l and the mass of the payload m are included in the nonlinear terms to take into account the variation in operating conditions. The nonlinear functions $f_1$ and $f_2$ can be reformulated in the linear-in-the-parameters regression functions as

$$v\left(k\right)={\varphi }_1^T\left(k\right){\theta }_1+e_1\left(k\right)$$

(3)

$$\alpha \left(k\right)={\varphi }_2^T\left(k\right){\theta }_2+e_2\left(k\right)$$

(4)

where ${\varphi }_1$ and ${\varphi }_2$ are the vectors of the nonlinear regressors composed of the lagged input and output and the operating parameters l and m, while ${\theta }_1$ and ${\theta }_2$ are vectors of the linear parameters that can be estimated using the least squares technique. The crane position at time step $k$ is approximated by an explicit Euler integration of the velocity to obtain

$$x\left(k\right)=x\left(k-1\right)+t_{s}v\left(k\right)$$

(5)

where $t_s$ is the sample time.

2.2. Genetic Programming-Based Identification

Under the modeling assumptions given in (1) and (2), the NARX model structure is selected with additional operating parameters in order to identify the input–output nonlinear crane dynamics, which can be presented in the general form

$$y\left(k\right)=f\left(\xi \left(k-1\right),\dots ,\xi \left(k-N_{\xi }\right),y\left(k-1\right),\dots ,y\left(k-N_y\right),p\right)+e\left(k\right)$$

(6)

where y is the output, ξ is the input excitation signal, e is Gaussian white noise, and $p$ is the operating parameter.

The coefficients of the linear-in-the-parameters model can be estimated using the least squares method, as follows:

$${\theta }^{\mathrm{*}}=\arg \underset{\theta }{\min }{‖\mathsf{\Phi }\theta -y‖}_2^2$$

(7)

where $\mathsf{\Phi }$ is a matrix, whose columns are the nonlinear regressors.

It is assumed that the dynamics of the system can be described by a subset of the regressors in $\mathsf{\Phi }$; therefore, the coefficients $\theta$ are sparse. In order to promote the sparsity of the solution, the $l_0$ penalty is added to (7) to obtain

$$F\left(\theta \right)=g\left(\theta \right)+h\left(\theta \right)=\frac{1}{2}{‖\mathsf{\Phi }\theta -y‖}_2^2+\lambda {‖\theta ‖}_0$$

(8)

where ${‖\theta ‖}_0$ is the $l_0$ pseudonorm, i.e., ${‖\theta ‖}_0={\sum }_{i=1}^N{𝟙}_{{\theta }_i\ne 0}$. Problem (8) is nonconvex and NP-hard [44], therefore several methods have been proposed to obtain sparse solutions, such as replacing the $l_0$ pseudonorm with the $l_1$ norm, smoothly clipped absolute deviation (SCAD), minimax concave penalty (MCP), sequentially thresholded least squares (STLS), and proximal gradient methods. In this paper, a suboptimal solution is found by using monotone accelerated proximal gradient descent (mAPG) [45], which converges to a critical point of $F\left(\theta \right)$. The proximal mapping of a function $h$ is given by the proximal operator

$$\mathrm{p}\mathrm{r}\mathrm{o}{\mathrm{x}}_h\left(q\right)=\underset{\mu }{\mathit{arg}\mathit{min}}h\left(\mu \right)+\frac{1}{2}{‖\mu -q‖}_2^2$$

(9)

which for $h\left(·\right)={\lambda ‖·‖}_0$ is separable and results in the hard thresholding operator with threshold τ, as follows:

$${\left[T_{\tau }\left(q\right)\right]}_s=\left\{\begin{array}{cc}0,& \left|q_s\right|\le \tau \\ q_s,& \mathrm{otherwise}\end{array}\right.,1\le s\le n$$

(10)

where $s$ is the element index and $n$ is the number of elements.

The monotone accelerated proximal gradient descent is initialized with $z^{\left(1\right)}={\theta }^{\left(1\right)}={\theta }^{\left(0\right)},t_0=0,t_1=1$, and updates ${\theta }^{\left(j+1\right)}$ according to

$${\theta }^{\left(j+1\right)}=\left\{\begin{array}{ll}{z}^{\left(j+1\right)},& F\left(z^{\left(j+1\right)}\right)\le F\left(v^{\left(j+1\right)}\right)\\ {\varrho }^{\left(j+1\right)},& \mathrm{otherwise}\end{array}\right.$$

(11)

where

• $z^{\left(j+1\right)}=\mathrm{p}\mathrm{r}\mathrm{o}{\mathrm{x}}_{{\gamma }_r\lambda {‖·‖}_0}\left(r^{\left(j\right)}-{\gamma }_r\nabla g\left(r^{\left(j\right)}\right)\right)=T_{\sqrt{2\lambda {\gamma }_r}}\left(r^{\left(j\right)}-{\gamma }_r\nabla g\left(r^{\left(j\right)}\right)\right)$,
• ${\varrho }^{\left(j+1\right)}=\mathrm{p}\mathrm{r}\mathrm{o}{\mathrm{x}}_{{\gamma }_{\theta }\lambda {‖·‖}_0}\left({\theta }^{\left(j\right)}-{\gamma }_{\theta }\nabla g\left({\theta }^{\left(j\right)}\right)\right)=T_{\sqrt{2\lambda {\gamma }_{\theta }}}\left({\theta }^{\left(j\right)}-{\gamma }_{\theta }\nabla g\left({\theta }^{\left(j\right)}\right)\right)$,
• $r^{\left(j\right)}={\theta }^{\left(j\right)}+\frac{t_{j-1}}{t_j}\left(z^{\left(j\right)}-{\theta }^{\left(j\right)}\right)+\frac{t_{j-1}-1}{t_j}\left({\theta }^{\left(j\right)}-{\theta }^{\left(j-1\right)}\right)$,
• $t_{j+1}=\frac{\sqrt{1+4t_j^2}+1}{2}$,
• where $j$ is the iteration number and ${\gamma }_{\theta }$ and ${\gamma }_r$ are the step sizes. The algorithm has been proven to converge to a critical point when ${\gamma }_{\theta },{\gamma }_r<\frac{1}{L}$ [45], where $L$ is the Lipschitz constant of $\nabla g$. However, in order to speed up convergence, the Barzilai–Borwein [46] method is used to initialize the step size, and a line search is used to find the step sizes ${\gamma }_{\theta }$ and ${\gamma }_r$ that result in sufficient descent.

$$F\left({\varrho }_{j+1}\right)\le F\left({\theta }_j\right)+\delta {‖{\varrho }_{j+1}-{\theta }_j‖}^2$$

(12)

$$F\left(z_{j+1}\right)\le F\left(r_j\right)+\delta {‖z_{j+1}-r_j‖}^2$$

(13)

where $\delta >0$ is a small positive constant.

The regressors are evolved using grammar-guided genetic programming, which overcomes the need for the closure property required in traditional genetic programming to produce valid offspring. Additionally, user knowledge of the system can be incorporated into the context-free grammar to restrict the search space. A context-free grammar $G$ is the 4-tuple $(S,N,\Sigma ,P)$, where $S\in N$ is the start symbol, $N$ is a finite set of nonterminal symbols, $\Sigma$ is an alphabet of terminal symbols, and $P$ is the set of production rules of the form $X\to \beta$, where $X\in N$ and $\beta \in {\left(N\cup \Sigma \right)}^{*}$.

For identifying the crane system dynamics, the grammar was constructed in such a manner that the model regressors are monomials, thereby resulting in a polynomial NARX model; however, the grammar allows for additional operations on the parameters $l$ and $m$, such as division and square root. The grammar, in Backus–Naur form (BNF), is given in

Table 1. Grammar used for polynomial NARX models.

〈S〉::=	〈exp〉 | 〈T〉
〈exp〉::=	×〈exp〉〈exp〉 | ×〈expp〉〈exp〉 | 〈T〉
〈expp〉::=	×〈expp〉〈expp〉 | ÷〈expp2〉〈expp2〉 | √〈expp〉 | 〈Tp〉
〈expp2〉::=	+〈expp2〉〈expp2〉 | ×〈expp2〉〈expp2〉 | √〈expp2〉 | 〈Tp〉
〈T〉::=	$\xi (k-1)|\cdots$$\left|\xi \right(k-N_{\xi }\left)\right|y(k-1)|\cdots$ | y(k − Nu)
〈Tp〉::=	m | l

.

In every generation, new individuals, called offspring, are produced by selecting individuals from the existing population, called parents, that pass through a stochastic variation operator. A total of four different variation operators were used: two binary variation operators, namely, subtree crossover and high-level crossover, and two unitary variation operators, namely, subtree mutation and point mutation.

Since every model term is represented by a separate tree, the selection of a parent occurs in two stages. In the first stage, the model is selected from the population by tournament selection, after which a model term is selected from the model to undergo variation. The model term is selected with the probability of selection given by

$$P=\frac{\exp {\overline{\theta }}_i}{{\sum }_{i=1}^M\exp {\overline{\theta }}_i}$$

(14)

where ${\overline{\theta }}_i=\frac{{\sum }_{i=1}^M\left|{\theta }_i\right|-\left|{\theta }_i\right|}{{\sum }_{i=1}^M\left|{\theta }_i\right|}$ is the normalized term coefficient.

In both crossover methods, two parents are selected; in subtree crossover, two root nodes of the subtrees are selected, and the subtrees are switched in a manner so that the resulting offspring is syntactically correct. In high-level crossover [47], the selected model terms are switched between the two parents. In both mutational methods, only one parent is selected from the population. Subtree mutation is similar to subtree crossover, but instead of a subtree from a second parent, a new subtree is randomly generated. In point mutation, a single terminal node is changed in a manner so that the resulting offspring is syntactically correct. The genetic programming parameters used for the identification of the overhead crane dynamics are given in

Table 2. GP parameters.

Parameters	Values
Population size	100
Number of generations	500
Initialization method	PTC2
Max tree depth during initialization	6
Max number of candidate model terms	50
Tournament size	2
Subtree crossover probability	0.75
High-level crossover probability	0.15
Mutation probability	0.1
Max number of proximal evaluation iterations	5000

.

The initialization method used was the probability tree creation 2 (PTC2) algorithm [48]. Due to the stochastic nature of genetic programming, each run could result in a possibly different solution; therefore, a total of 15 runs were carried out, and the solutions were compared in terms of accuracy, measured by using the mean square error, and complexity, measured by the number of model terms included in the model. One model from the Pareto frontier is then selected as the model used in the nonlinear model predictive control algorithm.

2.3. NMPC with PSO Solver

An NMPC controller with a PSO metaheuristic algorithm as the solver of the constrained optimal control problem was developed, implemented, and tested in real-time experiments on a laboratory stand. PSO is a stochastic optimization technique that was originally developed by Eberhart and Kennedy [49]. There have been several variations of the PSO algorithm, including utilizing different topologies, such as the star topology, which uses the swarm’s global best position, or the ring topology, in which a particle uses the local best position in its neighborhood. While the original PSO formulation used the star topology, it has been proposed to use the ring topology and limit the star topology to unimodal problems; however, according to [50], the star topology is still acceptable for multi-objective optimization problems. An empirical study by Engelbrecht [51] on 60 benchmark functions showed that the two topologies perform similarly, with the star topology performing slightly better in terms of success rate and efficiency, while the ring topology performed slightly better in terms of consistency. Additionally, PSO variants with sub-swarm populations have been proposed, such as dynamic multi-swarm PSO [52] and symbiotic particle swarm optimization [53]. A review of the different PSO variants is presented in [54].

In the proposed control scheme, at each time step, $k$, the NMPC-GP receives the current state of the plant and solves the constrained optimization problem to find a sequence of control actions $u\left(k+i-1\right)$ defined over the control horizon $N_u$ to minimize the cost function, which depends on the deviation of the process variables $x\left(k+i\right)$ and $\alpha \left(k+i\right)$ from the set point $x_r\left(k+i\right)$ as well as the control effort required to change the system states over the prediction horizon $N_p$. The PSO algorithm solves the optimization problem with the constraint $u(k+i+1)\in \mathbb{U}$ directly as one large optimization problem, and, therefore, there is no requirement for breaking down the problem into smaller subproblems as is the case with dynamic programming. The future trajectories of the states are predicted using (3–5) based on measurements at the current time instant $k$ and the candidate control sequence $u\left(k+i+1\right)$ obtained iteratively using the PSO algorithm (Figure 2).

The proposed cost function includes the sum of the squared prediction errors between the desired and predicted crane positions on the prediction horizon, which corresponds to shifting the equilibrium to the desired position; the sum of the weighted squares of the predicted payload angle on the prediction horizon; and the weighted squares of the control increments on the control horizon, as given in

$$J\left(k+i\right)=\sum _{i=1}^{N_p}{\lambda }_1{\left(e\left(k+i\right)\right)}^2+\sum _{i=1}^{N_p}{\lambda }_2{\left(\alpha \left(k+1\right)\right)}^2+\sum _{i=1}^{N_u}{\lambda }_3{\left(\mathsf{\Delta }u\left(k+i-1\right)\right)}^2$$

(15)

where

• $e\left(k+i\right)=x\left(k+i\right)-x_r\left(k+i\right)$,
• $\mathsf{\Delta }u\left(k+i-1\right)=u\left(k+i-1\right)-u\left(k+i-2\right)$,
• ${\lambda }_1=\left\{\begin{array}{c}\begin{array}{cc}0& if\left|e\left(k+i\right)\right|\le e_a\end{array}\\ \begin{array}{cc}1& if\left|e\left(k+i\right)\right|>e_a\end{array}\end{array}\right.$,
• ${\lambda }_2=\left\{\begin{array}{c}\begin{array}{cc}0& if\left|\alpha \left(k+i\right)\right|\le {\alpha }_a\end{array}\\ \begin{array}{cc}{\lambda }_{20}& if\left|\alpha \left(k+i\right)\right|>{\alpha }_a\end{array}\end{array}\right.$.

The position error, $e\left(k+i\right)$, is the difference between the predicted position $x\left(k+i\right)$ and the reference signal $x_r$ at time $k+i$; $\alpha \left(k+i\right)$ is the predicted sway angle of the payload; $\mathsf{\Delta }u\left(k+i-1\right)$ expresses the control effort, which is calculated as the predicted change in the control signal along the control horizon; and ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are the weighting coefficients. The control goal is to reach the desired set point within a specified tolerance, and therefore the cost function coefficients ${\lambda }_1$ and ${\lambda }_2$ are set to ${\lambda }_1=0$ and ${\lambda }_2=0$ when the expected crane positioning accuracy $e_a$ is satisfied and the payload sway oscillation is within the acceptable tolerance ${\alpha }_a$, respectively; otherwise, nonzero weights are applied (${\lambda }_1=1$, ${\lambda }_2={\lambda }_{20}$). The tolerance ${\alpha }_a$ of the pendulum oscillation corresponds to the range of the payload linear deviation $d_p$ from the stable equilibrium point, where the linear deviation of a payload depends on the rope length $d_p=l\sin \alpha$, and for small values of the sway angle, it can be approximated as $d_p\approx l\alpha$. In the experiments carried out in this work, the crane and payload positioning accuracy was assumed to be $e_a=d_p=0.02$ m. Thus, during the experiments, the desired range of magnitude of the residual sway was set as ${\alpha }_a=e_a/l=0.02/l$.

The PSO algorithm is adapted to determine the sequence of control signal values $u(k+i-1)$ over the control horizon $N_u$ to minimize the cost function (15). Consider a population of candidate solutions (a swarm of particles) in the constrained $N_u$-dimensional search space, where each particle’s position is represented by the coordinates corresponding to the sequence of control signals over the control horizon $N_u$, as follows:

$$u_j=\left[u\left(k+i-1\right),u\left(k+1\right),\dots ,u\left(k+N_u-1\right)\right]$$

(16)

satisfying the box constraints

$$u_{\mathrm{m}\mathrm{i}\mathrm{n}}\le u_j\le u_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

(17)

where $j=1,2,\dots ,n_s$; and $n_s$ is the size of the swarm. At each time step, $t$, called an iteration, the $j^{\mathrm{t}\mathrm{h}}$ particle’s position and velocity are updated in the hyperspace according to

$${\tilde{u}}_j^t=w{\tilde{u}}_j^{t-1}+c_1{r}_1\left(u_{j,p}-u_j^{t-1}\right)+c_2{r}_2\left(u_g-u_j^{t-1}\right)$$

(18)

$$u_j^t=u_j^{t-1}+{\tilde{u}}_j^t$$

(19)

where $u_j^{t-1}$ and ${\tilde{u}}_j^{t-1}$ are the particle’s position and velocity in the previous iteration, respectively; $u_{j,p}$ and $u_g$ are the personal and global best solutions, respectively; $w$ is called the inertia weight; positive constants $c_1$ and $c_2$ are called the acceleration coefficients and are responsible for the influence of the global and local attractors on the particle’s velocity; and constants $r_1$ and $r_2$ are drawn uniformly from the interval $\left[0,1\right]$. In this work, following [55], the inertia weight varies linearly between 0.9 and 0.4, together with the acceleration coefficients of 1.5, to ensure an efficient balance between exploration (the ability of the particles to search wide regions of the solution space) and exploitation (the ability of the particles to refine the solution in the neighborhood of an attractor).

Since the velocity of a particle refers to the change in a control signal, which is constrained by the lower and upper limits ($u_{min},u_{max}$), the velocity in formula (18) is updated by taking into account the control signal limitation:

$${\tilde{u}}_j^t=\max \left(\min \left({\tilde{u}}_j^t,{\tilde{u}}_{j,\mathrm{m}\mathrm{a}\mathrm{x}}^t\right),{\tilde{u}}_{j,\mathrm{m}\mathrm{i}\mathrm{n}}^t\right)$$

(20)

where ${\tilde{u}}_{j,\mathrm{m}\mathrm{i}\mathrm{n}}^t=u_{\mathrm{m}\mathrm{i}\mathrm{n}}-u_j^{t-1}$ and ${\tilde{u}}_{j,\mathrm{m}\mathrm{a}\mathrm{x}}^t=u_{\mathrm{m}\mathrm{a}\mathrm{x}}-u_j^{t-1}$.

The constraint described in (20) is the absorption constraint, i.e., if the particle’s velocity would cause the particle to leave the feasible region, then the boundary “absorbs” the excess velocity, causing the particle to stay on the boundary. The personal best position $u_{j,p}$ of the j-th particle and the global best position of swarm $u_g$ are updated at each iteration. The iterative process is repeated at each time step $k$, starting with a randomly initialized swarm of particles and their velocities satisfying condition (20). The initial population size is $n_s$ − 1 and is completed by the global best solution from the previous time step $k–1$ to enhance the convergence speed. The PSO searching process is terminated when the maximum number of iterations is reached and the first coordinate of the global best particle $u\left(k\right)$ is applied as the control signal at the current instant $k$.

2.4. Digital Feedback Controller

The NMPC-GP with the PSO solver is compared with a discrete-time feedback control system developed using the pole placement technique. The crane dynamics are decomposed in a similar manner as in Section 2.1 into the actuated and unactuated parts and identified using an output-error approach to obtain a discrete-time relation between the crane velocity v(k) and the control input u(k), approximated by a first-order discrete-time model (21), while the relation between the payload sway angle α(k) and the velocity v(k) is given by a second-order discrete-time model (22):

$$v\left(k\right)=\frac{D\left(z^{-1}\right)}{C\left(z^{-1}\right)}u(k)$$

(21)

$$\alpha \left(k\right)=\frac{B\left(z^{-1}\right)}{A\left(z^{-1}\right)}v(k)$$

(22)

where $z^{-1}$ is the backward shift operator and $A\left(z^{-1}\right)$, $B\left(z^{-1}\right)$, $C\left(z^{-1}\right)$, and $D\left(z^{-1}\right)$ are polynomials:

• $A\left(z^{-1}\right)=1+a_1{z}^{-1}+a_2{z}^{-2}$,
• $B\left(z^{-1}\right)=b_1{z}^{-1}+b_2{z}^{-2}$,
• $C\left(z^{-1}\right)=1+c_1{z}^{-1}$,
• $D\left(z^{-1}\right)=d_1{z}^{-1}$.

A digital control law at step k for parallel feedback loops of the crane position and payload sway angle is proposed as follows:

$$u\left(k\right)=\frac{\left(k_1+k_2{z}^{-1}\right)e\left(k\right)+\left(k_3+k_4{z}^{-1}\right)\alpha \left(k\right)}{1-k_5{z}^{-1}}$$

(23)

where $e\left(k\right)=x\left(k\right)-x_r\left(k\right)$, $x\left(k\right)$ is approximated by (5); $x_r\left(k\right)$ is the reference signal at time k; and $k_1$, $k_2$, $k_3$, $k_4$, and $k_5$ are the controller parameters calculated using the pole placement method. In this case, the 5th-order closed-loop characteristic polynomial can be presented as

$$\delta \left(z^{-1}\right)=qSK$$

(24)

where

• $q=\left[1,z^{-1},\dots ,z^{-5}\right]$,
• $S=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ a_1+c_1& 0& 0& 0& 0& -1\\ a_1{c}_1-a_1+a_2-c_1& t_s{d}_1& 0& -b_1{d}_1& 0& {-a}_1-c_1+1\\ {-a}_1{c}_1+a_2{c}_1-a_2& t_s{a}_1{d}_1& t_s{d}_1& \left(b_1-b_2\right)d_1& -b_1{d}_1& {-a}_1{c}_1+a_1-a_2+c_1\\ {-a}_2{c}_1& t_s{a}_2{d}_1& t_s{a}_1{d}_1& b_2{d}_1& \left(b_1-b_2\right)d_1& a_1{c}_1-a_2{c}_1+a_2\\ 0& 0& t_s{a}_2{d}_1& 0& b_2{d}_1& a_2{c}_1\end{array}\right]$,
• $K={\left[{1,k}_1,k_2,\dots ,k_5\right]}^T$.

The desired polynomial $\delta \left(z^{-1}\right)$ can be obtained by assuming the desired location of the closed-loop poles ${\lambda }_i$:

$$\delta \left(z^{-1}\right)=\prod _{i=1}^{n=5}\left(1-{\lambda }_i{z}^{-1}\right)$$

(25)

In [56], the desired closed-loop poles were assumed as ${\lambda }_i=exp\left(-{\omega }_n{t}_s\right)$, where ${\omega }_n=\sqrt{g/l}$ is the natural undamped frequency of the payload suspended on a rope of length l, while g = 9.81 ms⁻² is the acceleration gravity.

3. Results

The laboratory experimental setup hardware is illustrated in Figure 3. The double-girder crane bridge is driven by two 0.18 kW AC gear motors supplied by SV004iC5-1F 0.4kW frequency inverters (LS Industrial Systems Co., Ltd., Anyang-si, South Korea) that are controlled by the voltage signal u within the range ±10V. The incremental encoders with a resolution of 400, 2000, and 100 ppr (pulses per rotation) are used to measure position x, sway angle α, and rope length l, respectively. The encoders are installed on the wheels to measure position x relative to the crane bridge, under the trolley connected to a pair of fork arms embracing the hoisting cable to measure sway angle α, and on the hoisting drum to measure rope length l. The data from the sensors were sampled at 10 Hz using a PC (16 GB RAM, Quad Core 4 GHz Intel Core i7-6700 K CPU, Intel, Santa Clara, CA, USA) equipped with the IO card PCI1710HG, MS Windows 10 operating system, and Matlab software release R2020. During the experiments, the PLC was used to transmit the signals from the sensors to the PC and the control signal u to the actuators. A constant velocity Kalman filter was used to estimate the velocity of the crane v with the process noise and measurement noise that follow Gaussian distributions $N\left(0,{\sigma }_p^2\right)$ and $N\left(0,{\sigma }_m^2\right)$, where ${\sigma }_p=8$ and ${\sigma }_m=0.2$.

The identification model was obtained from input–output experimental data collected from a series of 10 experiments carried out at different operating points, with rope lengths l = {0.8, 1.1, 1.4, 1.7, 2.0} m and payload masses m = {10, 30, 50} kg. The input excitation signal u was a sequence of step functions with varying amplitude. The data from the experiment carried out with l = 1.1 m and m = 50 kg and the experiment carried out with l = 1.7 m and m = 10 kg were used as the testing data, while the rest were used as the training data to evolve the GP-NARX model. Similarly to [43], the training data are split randomly with a 70:30 ratio into two sets, where the first data set is used in the GP algorithm to determine the MSE used as the fitness function, while the second data set is used for model selection, i.e., the model with the lowest MSE on the second training data set is stored as the best model. The GP algorithm was run 15 times with the grammar and parameters given in Table 1 and Table 2, respectively, and the best model in each run was saved. The accuracy and complexity of the best models from each run are shown in Figure 4, with the green points corresponding to models on the Pareto frontier and the green pentagram representing the selected model. Nonlinear regressors ${\varphi }_1$ and ${\varphi }_2$ and coefficients ${\theta }_1$ and ${\theta }_2$ of the crane velocity (3) and pendulum sway (4) models, respectively, are presented in

Table 3. Crane velocity and sway model terms and coefficients.

Velocity Model (3)		Sway Model (4)
$\mathbf{Coefficients}{\mathit{\theta }}_1$	$\mathbf{Regressors}{\mathit{\varphi }}_1$	$\mathbf{Coefficients}{\mathit{\theta }}_1$	$\mathrm{Regressors}{\mathit{\varphi }}_2$
0.0464	$u\left(k-7\right)$	−29.5814	$\alpha {\left(k-2\right)}^3$
0.0508	$u\left(k-3\right)u\left(k-4\right)u\left(k-9\right)$	0.0208	$v\left(k-3\right)$
0.9514	$v\left(k-1\right)$	22.8790	$\alpha \left(k-1\right)\alpha \left(k-2\right)\alpha \left(k-3\right)$
0.1691	$v\left(k-4\right)u\left(k-1\right)u\left(k-4\right)$	0.0364	$l^2\alpha \left(k-2\right)$
0.1263	$v\left(k-4\right)u\left(k-1\right)u\left(k-2\right)u\left(k-3\right)u\left(k-6\right)$	9.9522	$l\alpha \left(k-1\right)\alpha \left(k-3\right)\alpha \left(k-4\right)$
−0.1294	$v\left(k-4\right)$	−7.8066	$\alpha \left(k-2\right)\alpha \left(k-3\right)\alpha \left(k-4\right)$
0.0608	$v\left(k-4\right)u{\left(k-5\right)}^2$	0.4542	$\alpha \left(k-2\right)$
−0.1743	$v\left(k-2\right)v{\left(k-4\right)}^{2}u\left(k-1\right)u\left(k-4\right)$	−0.0205	$v\left(k-1\right)$
−0.0904	$v\left(k-3\right)u\left(k-4\right)u\left(k-5\right)$	0.8585	$\alpha \left(k-1\right)$
−0.1309	$u\left(k-1\right)u\left(k-2\right)u\left(k-4\right)$	−288.7701	$l^{2}a{\left(k-1\right)}^2\alpha \left(k-2\right)\alpha \left(k-3\right)\alpha \left(k-4\right)$
0.1115	$u\left(k-4\right)$	33.9964	$l\alpha {\left(k-1\right)}^2\alpha \left(k-4\right)$
		−17.9259	$l^2\alpha {\left(k-1\right)}^2\alpha \left(k-4\right)$
		−0.5518	$\alpha \left(k-4\right)$

.

The GP-NARX model had 11 model terms selected for the velocity submodel (3), while the sway submodel (4) had 13 model terms selected. The mass of the payload did not appear in any of the terms, making the models independent of the payload mass. Figure 5 presents the input signal u applied to excite the underactuated part of the overhead crane in the two experiments carried out to obtain testing data at operating points {l = 1.1 m, m = 50 kg} and {l = 1.7 m, m = 10 kg}, and the 20-step-ahead prediction (crane velocity and payload sway angle) performance of the GP-NARX models. The models were then experimentally verified with the NMPC-GP control scheme in terms of robustness against rope length and payload mass variation, as well as anti-disturbance performance.

The NMPC with the GP-NARX prediction model and the PSO solver, denoted NMPC-GP controller, was tested in experiments carried out with a sample time $t_s=0.1$ s and under the constraints on the control signal in the range −8 V ≤ u ≤ 8 V. The prediction and control horizon were established as $N_P=20$ and $N_c=5$, respectively, while the swarm size and number of iterations of the PSO were set in the experiments as $n_s=40$ and $n_t=100$, respectively. The error tolerance for the crane and payload positioning was assumed to be $\pm$0.02 m; thus, the desired range of magnitude of the residual sway was approximated in the controller as ${\alpha }_a\approx 0.02/l$ rad. The weighting coefficients ${\lambda }_{20}=70$ and ${\lambda }_3=0.05$ of the fitness function (15) were selected empirically during the experiments. The NMPC algorithm with the PSO solver was implemented using a C-MEX S-function in Matlab R2020. The mean execution time and standard deviation for 100 runs of the NMPC algorithm with the above-mentioned settings were 2.72 and 0.0424 milliseconds, respectively.

Three cases of experiments were carried out to investigate the performance and prove the feasibility and effectiveness of the developed control approach in terms of the controller’s robustness against rope length and payload mass variation, in the case of an external disturbance affecting the pendulum, and by comparing the NMPC-GP with a linear discrete-time controller. The controller’s performance was tested with both the operating points considered during the identification experiments, with rope lengths and payload masses not included in the training data.

In Case 1, the experiments were performed for different operating points corresponding to varying rope lengths (l = {0.8, 1.4, 2.0}) m and payload masses (m = {10, 30, 50} kg), while the desired set point was $x_r=1\mathrm{m}$. The control performance of the NMPC-GP controller is presented in

Table 4. Performance of NMPC-GP in terms of rope length and payload mass variation.

Rope Length $\mathit{l}\left[\mathbf{m}\right]$	Payload Mass $\mathit{m}\left[\mathbf{k}\mathbf{g}\right]$	Rise Time ${\mathit{t}}_{\mathit{r}}\left[\mathbf{s}\right]$	Settling Time ${\mathit{t}}_{\mathit{s}}\left[\mathbf{s}\right]$	Maximum Sway Angle ${\mathit{\alpha }}_{\mathit{m}\mathit{a}\mathit{x}}\left[\mathbf{r}\mathbf{a}\mathbf{d}\right]$	Maximum Deviation ${\mathit{d}}_{\mathit{m}\mathit{a}\mathit{x}}\left[\mathbf{m}\right]$	Residual Sway ${\mathit{\alpha }}_{\mathit{r}}\left[\mathbf{r}\mathbf{a}\mathbf{d}\right]$	Residual Deviation ${\mathit{d}}_{\mathit{r}}\left[\mathbf{m}\right]$
0.8	10	5.7	7.2	0.080	0.064	0.004	0.0032
	30	5.5	6.5	0.079	0.063	0.012	0.0096
	50	5.4	6.4	0.076	0.061	0.006	0.0048
1.4	10	5.5	6.5	0.052	0.073	0.001	0.0014
	30	5.7	7.1	0.058	0.081	0.001	0.0014
	50	5.8	7.0	0.055	0.077	0.001	0.0014
2.0	10	5.7	6.6	0.044	0.088	0.003	0.006
	30	5.7	6.5	0.044	0.088	0.008	0.016
	50	5.8	6.6	0.046	0.092	0.008	0.016

in terms of rise time $t_r$, settling time $t_s$, maximum magnitude of the pendulum oscillation ${\alpha }_{max}$, and residual sway ${\alpha }_r$, which is also converted to linear deviation from an equilibrium point. In addition, Table 4 includes data on the maximum deviation and residual deviation of the rope-suspended payload ($d_{max}$ and $d_r$, respectively). The control input u, determined by the NMPC-GP controller, as well as the crane velocity, position, and payload sway angle are depicted in Figure 6 and Figure 7a. The crane position rise time is within 5.4–5.8 s, while the settling time range is within 6.4–7.2 s; after this time, the residual oscillation remains within the expected range (the linear payload deviation from equilibrium point is less than 0.02 m). The maximum value of the payload sway angle varies with rope length, while the influence of the payload mass is minor within the range tested during the experiments.

In Case 2, the performance of the NMPC-GP was verified in the presence of an external disturbance affecting the pendulum. It is presented in Figure 7b. The experiment was carried out for a rope length l = 2.0 m and a payload mass m = 10 kg. In the time interval between 1.8 and 2.8 s, the pendulum is displaced from −0.04 rad to 0.15 rad, the control signal varies between −4.8 V and 8 V (corresponding to a velocity change between −0.09 m/s and 0.24 m/s), and the crane and payload positions are settled at the desired value with a settling time of 9 s with the suppression of the suspended payload oscillation within the allowable tolerance range of below 0.01 rad.

In Case 3, the proposed NMPC-GP is compared with the DFC developed using the pole placement method in Section 2.3. The comparison was carried out for two operating points not included in the training of the GP-NARX model: rope length l = 1.1 m and l = 1.7 m, and payload mass m = 10 kg. Using input–output data obtained from the experiment carried out with a rope length l = 1.1 m, a discrete-time output-error model identification was performed to parametrize models (13) and (14), as follows:

$$v\left(k\right)=\frac{0.0031z^{-1}}{1-0.8942z^{-1}}u\left(k\right),\alpha \left(k\right)=\frac{-0.1409z^{-1}+0.1423z^{-2}}{1-1.9302z^{-1}+{0.9999z}^{-2}}u\left(k\right)$$

(26)

The parameters of control law (15) were determined using the pole placement method as $k_1=162.471$, $k_2=-141.906$, $k_3=394.426$, $k_4=-364.362$, and $k_5=0.1078$, by choosing the desired closed-loop poles as ${\lambda }_i=0.7865$ (i = 1, 2, …, 5) to obtain a similar settling time for the crane position as in the case of the NMPC-GP. The control performances are compared in

Table 5. The performances of NMPC-GP and DFC.

Rope Length $\mathit{l}\left[\mathbf{m}\right]$	Controller	Rise Time ${\mathit{t}}_{\mathit{r}}\left[\mathbf{s}\right]$	Settling Time ${\mathit{t}}_{\mathit{s}}\left[\mathbf{s}\right]$	Maximum Sway Angle ${\mathit{\alpha }}_{\mathit{m}\mathit{a}\mathit{x}}\left[\mathbf{r}\mathbf{a}\mathbf{d}\right]$	Residual Sway ${\mathit{\alpha }}_{\mathit{r}}\left[\mathbf{r}\mathbf{a}\mathbf{d}\right]$
1.1	NMPC-GP	5.7	7.1	0.063	0.005
	DFC	5.2	7.3	0.073	0.005
1.7	NMPC-GP	5.8	6.8	0.045	0.003
	DFC	5.7	7.7	0.045	0.012

and Figure 8. The DFC performs similarly to the NMPC-GP in the experiment with l = 1.1 m in terms of settling time of crane position and payload residual oscillation. However, the NMPC-GP outperforms the DFC in the reduction in payload oscillations in the transient state. The first two amplitudes of oscillation are reduced by the NMPC-GP to 0.063 rad and 0.041 rad, respectively, compared to 0.073 and 0.072 in the case of the DFC scheme. The NMPC-GP provides faster suppression of the payload sway below 0.02 rad in less than 3 s, while the DFC controller performs it in 5 s. In the experiment with rope length l = 1.7 m, there is performance deterioration of the DFC with respect to the settling time ($t_s=7.7\mathrm{s}$) when compared to the NMPC ($t_s=6.8\mathrm{s}$). The predictive controller performs better and is faster in reducing the pendulum oscillation with a residual sway of 0.003 rad compared to 0.012 rad in the case of the DFC.

4. Conclusions

This study conducted an experimental verification of a nonlinear model predictive control with particle swarm optimization as the solver for an overhead crane with a data-driven overhead crane dynamic model. Grammar-guided genetic programming with $l_0$ regularization was used to construct an input–output dynamic model of the overhead crane. The experimental data used for identification were obtained at 10 different operating points by varying the rope length and the payload mass. A total of 15 genetic programming runs were carried out and compared in terms of accuracy and complexity, and a model for the NMPC was selected from the Pareto frontier.

The NMPC cost function was selected to be a quadratic function of the payload positioning error, the sway angle, and the control input, which was solved with the PSO algorithm. Due to the physical limitations of the overhead crane, constraints on the control input were implemented as box constraints in the PSO algorithm. To prove the efficacy of the proposed method, three different experimental cases were carried out.

In the first case, the proposed nonlinear model predictive controller was tested at different operating points with the rope length ranging from 0.8 to 2.0 m and the payload mass ranging from 10 to 50 kg. The crane position rise time was within 5.4–5.8 s, while the settling time range was within 6.4–7.2 s. The residual oscillation was within the acceptable tolerance level of the linear payload deviation of $\pm$0.02 m for all operating points. In the second case, the payload was subjected to an external disturbance during positioning in the time interval 1.8–2.8 s. The proposed controller positioned the payload at the desired location with a settling time of 9 s and had residual oscillations within the acceptable tolerance level. In the third case, the nonlinear model predictive controller was compared with a digital feedback controller developed using the pole placement method with a data-driven linear output-error model. The NMPC had better oscillation suppression than the DFC in the transient state, as well as lower residual oscillations.

This study shows that the NMPC with PSO used to minimize the cost function of a data-driven dynamic model of an overhead crane is a viable alternative for the control of overhead cranes. The proposed method was able to position the payload with varying rope lengths and payload masses, as well as in the presence of external disturbances. A comparative study with a digital feedback controller was also performed, in which the proposed method performed better in the transient state and was able to achieve a residual oscillation within the desired tolerance.

The proposed identification method was only used to identify a NARX model; however, a NARMAX model would be able to model a larger class of systems by including lagged error terms. Additionally, the crane model was assumed to have high mechanical impedance, which in some industrial settings may no longer be valid. Therefore, in future work, we aim to expand the model structure from the NARX model to the more general NARMAX model structure and incorporate the coupling between the crane velocity and the payload sway into the model. We will also consider investigating the different prediction and control horizon lengths more comprehensibly.