## 2 Materials and Methods
### 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(k)=f\big{(}\xi(k-1),\ldots,\xi\big{(}k-N_{\xi}\big{)},y(k-1),\ldots,y\big{(} k-N_{y}\big{)},p\big{)}+e(k) \tag{6}\]  
where \(y\) is the output, \(\xi\) 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^{*}=\arg\!\min_{\theta}\|\Phi\theta-y\|_{2}^{2} \tag{7}\]  
where \(\Phi\) is a matrix, whose columns are the nonlinear regressors. Figure 1: Schematic diagram of an overhead crane in 2D operating space. It is assumed that the dynamics of the system can be described by a subset of the regressors in \(\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(\theta)=g(\theta)+h(\theta)=\frac{1}{2}\|\Phi\theta-y\|_{2}^{2}+\lambda\| \theta\|_{0} \tag{8}\]  
where \(\|\theta\|_{0}\) is the \(l_{0}\) pseudonorm, i.e., \(\|\theta\|_{0}=\sum_{i=1}^{N}\mathbb{I}_{\theta;\
eq 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(\theta)\). The proximal mapping of a function \(h\) is given by the proximal operator  
\[\text{prox}_{h}(q)=\underset{\mu}{argminh}(\mu)+\frac{1}{2}\|\mu-q\|_{2}^{2} \tag{9}\]  
which for \(h(\cdot)=\lambda\|\cdot\|_{0}\) is separable and results in the hard thresholding operator with threshold \(\tau\), as follows:  
\[\left[T_{\tau}(q)\right]_{s}=\left\{ \begin{matrix}0,&|q_{s}|\leq\tau\\ q_{s},&\text{otherwise}\end{matrix} ,1\leq s\leq n\right. \tag{10}\]  
where \(s\) is the element index and \(n\) is the number of elements. The monotone accelerated proximal gradient descent is initialized with \(z^{(1)}=\theta^{(1)}=\theta^{(0)},t_{0}=0,t_{1}=1\), and updates \(\theta^{(j+1)}\) according to  
\[\theta^{(j+1)}=\left\{ \begin{matrix}z^{(j+1)},&F\Big{(}z^{(j+1)}\Big{)}\leq F \Big{(}v^{(j+1)}\Big{)}\\ e^{(j+1)},&\text{otherwise}\end{matrix} \right. \tag{11}\]  
where \(j\) is the iteration number and \(\gamma_{\theta}\) and \(\gamma_{\tau}\) are the step sizes. The algorithm has been proven to converge to a critical point when \(\gamma_{\theta},\gamma_{\tau}<\frac{1}{L}\)[45], where \(L\) is the Lipschitz constant of \(\
abla 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_{\tau}\) that result in sufficient descent. \[F\big{(}q_{j+1}\big{)} \leq F\big{(}\theta_{j}\big{)}+\delta\big{\|}\phi_{j+1}-\theta_{ j}\big{\|}^{2} \tag{12}\] \[F\big{(}z_{j+1}\big{)} \leq F\big{(}r_{j}\big{)}+\delta\big{\|}z_{j+1}-r_{j}\big{\|}^{2} \tag{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.