Real-Time Weight Optimization of a Nonlinear Model Predictive Controller Using a Genetic Algorithm for Ship Trajectory Tracking

Abstract

This paper presents a weight optimization method for a nonlinear model predictive controller (NMPC) based on the genetic algorithm (GA) for ship trajectory tracking. The weight coefficients Q and R of the objective function in NMPC are obtained via the real-time optimization of the genetic algorithm instead of the trial and error method, which improves the efficiency and accuracy of the controller. In addition, targeted improvements are made to the internal crossover operator, mutation operator, crossover rate, and mutation rate of the genetic algorithm. The simulation comparison of trajectory tracking between NMPC with real-time-optimized weight coefficients and the one with constant coefficients is performed. Finally, the simulation result shows that the controller with real-time-optimized weight coefficients has a better trajectory tracking effect than that with constant weight coefficients.

1. Introduction

With the rapid development of deep-sea applications, ship automatic control has received unprecedented attention. Trajectory tracking is one of the key technologies used in ship intelligent control. Therefore, the research and design of trajectory tracking control to ensure the accurate tracking of the desired trajectory plays an important role in the safe and effective implementation of marine engineering applications [1,2,3].

Relevant scholars have made many efforts and achievements in the field of ship trajectory tracking control. Katayama and Aoki [4] used the characteristics of the Euler approximation model of ship dynamics to design a reduced order observer and output feedback controller, which can effectively track a preset straight-line trajectory. Qu et al. [5] proposed an unmanned ship exponential tracking controller based on an observer to estimate external disturbances and system uncertainties, with the controller able to ensure exponential convergence and strong robustness. Liu et al. [6] improved the nonlinear robust control algorithm based on sliding mode control theory to simulate the maneuverability of the surface ship.

However, the above control methods cannot take the system input constraints directly into account due to the physical limits of the thrusters. Considering that the model predictive control (MPC) has a superior ability to establish constraints on optimization problems, it has been widely used in various fields, such as semiconductor, energy, and aerospace, to achieve advanced control [7,8,9,10]. In the field of marine vessel control, Oh and Sun [11] combined the LOS guidance algorithm with MPC to implement the tracking control in the presence of input constraints. Li and Sun [12] proposed a new disturbance-compensated MPC algorithm for ship heading control. Kapetanović et al. [13] proposed applied linear model predictive control (LMPC) for the trajectory tracking of underactuated ships under disturbance and input constraint saturation. Furthermore, since ship dynamics is highly nonlinear, in order to improve the control performance, it stimulates the development of nonlinear model predictive control (NMPC) [14]. Abdelaal et al. [15] proposed a combined NMPC for the position and velocity tracking of underactuated surface vessels. Yang et al. [16] designed a robust NMPC scheme based on a nonlinear disturbance observer for the trajectory tracking control of a dynamic positioning ship.

It can be seen from the above research that NMPC can perform well for trajectory tracking. However, in traditional NMPC, the most important parameters of weight coefficients in objective function are generally obtained through trial or experience, which is time-consuming and easily falls into local optimization. Therefore, many studies have been conducted for the optimal weight tuning of model predictive control. Based on the real-time weight adjustment strategy, Zhao et al. [17] designed a spacing control law for the adaptive cruise control system, which can adjust the weight coefficient according to actual output demand. In addition, the development of intelligent algorithms also provides good ideas for this. For example, the genetic algorithm (GA) is combined with NMPC for parameter optimization because of its superior characteristics of global optimization [18]. In addition, Ramasamy et al. [19] used GA to determine the MPC weights to minimize overall energy utilization, reduce tracking error, and realize optimal control of the cement kiln. Essa et al. [20] adopted the combination of the Rhododendron search algorithm and GA to adjust the weight coefficients of MPC, and gained strong stability and robustness in the electro-hydraulic servo system. Zhang and Zhuan [21] proposed an improved genetic algorithm of the variable chromosome length coevolution method to optimize the weight matrices Q and R in the electromagnetic isolation system, and received better damping performance.

Aiming to improve the optimization performance of GA, the crossover operator and mutation operator have been improved upon in many studies. For example, Tang and Tseng [22] proposed an adaptive directed mutation operator, which can enhance the global search ability and accelerate the convergence of GA by integrating local directional and adaptive random search strategy. Chuang et al. [23] proposed the use of the direction-based crossover operator and the dynamic random mutation operator. The crossover operator uses the ranking information of chromosome fitness value in the population to generate offspring chromosomes. The mutation operator generates new chromosomes through an adaptively adjusted mutation step and random disturbance vector. However, the crossover operators cannot effectively produce offspring in some specific areas where the initial population is located, which reduces the performance of the algorithm. For the mutation operators, due to their random variability, on the one hand, the optimal solution of the current population may be mutated, thus destroying the optimal solution. On the other hand, it may not make any mutation to the worst solution, which cannot improve the population diversity. In addition, the crossover rate and the mutation rate are two important control parameters in GA, and their values affect the convergence of the algorithm. The traditional determination method is to take a constant value according to the results of many experiments, which cannot ensure the algorithm to finally converge to the global optimal solution.

On the basis of the above analysis, this paper proposes a novel real-time weight optimization method of a nonlinear model predictive controller using a genetic algorithm for ship trajectory tracking, and the main contributions are as follows:

(1)

NMPC is applied to the ship trajectory tracking control, and the genetic algorithm is used to optimize the weight coefficients of the objective function in NMPC in real time;

(2)

The crossover operator, mutation operator, crossover rate, and mutation rate in the genetic algorithm are improved to enhance the performance of the genetic algorithm.

The rest of the paper is organized as follows: Section 2 introduces the structure of NMPC for ship control. Section 3 puts forward the GA-based real-time weight optimization algorithm of NMPC. In Section 4, simulations are conducted and results are discussed. The last section gives the conclusions.

2. Nonlinear Model Predictive Control for Trajectory Tracking

2.1. Ship Mathematical Model

For surface ships, the heave, pitch, and roll motions are neglected, and only the horizontal motions of surge, sway, and yaw need to be considered. According to Fossen [24,25], the inertial earth-fixed frame and the body-fixed frame attached to the moving vessel are constructed to build the nonlinear mathematical model of a surface ship as:

$$\begin{array}{ccc}\hfill \dot{\mathbf{\eta }}& =& \mathit{J}\left(\psi \right)\mathbf{\nu },\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \mathit{M}\dot{\mathbf{\nu }}& =& -\mathit{D}\mathbf{\nu }+\mathbf{\tau }+{\mathit{J}}^{\mathit{T}}\left(\psi \right)\mathit{d},\hfill \end{array}$$

(2)

where $\mathbf{\eta }={[x,y,\psi ]}^T$ represents the ship’s position and heading vector in the inertial earth-fixed frame, $\mathbf{\nu }={[u,v,r]}^T$ is the velocity vector in the body-fixed frame, $\mathit{J}\left(\psi \right)$ is the transfer matrix, $\mathit{M}$ is the mass matrix, including the added mass, and $\mathit{D}$ is the damping matrix. These matrices are given as:

$$\begin{array}{ccc}\hfill \mathit{J}\left(\psi \right)& =& \left[\begin{array}{ccc}cos\left(\psi \right)& -sin\left(\psi \right)& 0\\ sin\left(\psi \right)& cos\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right], \mathit{M}=\left[\begin{array}{ccc}m-X_{\dot{u}}& 0& 0\\ 0& m-Y_{\dot{v}}& m{X}_g-Y_{\dot{r}}\\ 0& m{X}_g-Y_{\dot{r}}& I_z-N_{\dot{r}}\end{array}\right],\hfill \\ \hfill \mathit{D}& =& \left[\begin{array}{ccc}-X_u& 0& 0\\ 0& -Y_v& -Y_r\\ 0& -N_v& -N_r\end{array}\right],\hfill \end{array}$$

(3)

In addition, $\mathbf{\tau }={[{\tau }_u,{\tau }_v,{\tau }_r]}^T$ is the control force vector, composed of surge force ${\tau }_u$, sway force ${\tau }_v$, and yaw moment ${\tau }_r$, which are limited by the nonlinear saturation characteristic of the propeller:

$$\begin{array}{c}\hfill {\mathbf{\tau }}_{\mathit{min}}\le \mathbf{\tau }\le {\mathbf{\tau }}_{\mathit{max}},\end{array}$$

(4)

where ${\mathbf{\tau }}_{\mathit{min}}={[{\tau }_{umin},{\tau }_{vmin},{\tau }_{rmin}]}^T$, ${\mathbf{\tau }}_{\mathit{max}}={[{\tau }_{umax},{\tau }_{vmax},{\tau }_{rmax}]}^T$ are the minimal and maximal force vectors, respectively, that can be provided by the thrusters. Furthermore, $\mathit{d}={[d_1,d_2,d_3]}^T$ refers to the unknown time-varying disturbance vector composed of wind, wave, and current environmental interference and model uncertainties.

Generally, the surface ship model can be expressed by the state-space equation as:

$$\begin{array}{ccc}\hfill \dot{\mathit{x}}\left(\mathit{t}\right)& =& \mathit{A}\left(t\right)\mathit{x}+\mathit{B}\left(t\right)\mathit{u}+\mathsf{\Gamma }\left(t\right)\mathit{d}\left(t\right),\hfill \\ \hfill \mathit{z}\left(t\right)& =& \mathit{C}\mathit{x}\left(t\right),t\ge 0,\hfill \end{array}$$

(5)

where $\mathit{x}={[x,y,\psi ,u,v,r]}^T$ is the state vector, $\mathit{u}=\mathbf{\tau }$ is the control input vector, and z is the output vector. In combination with Equation (1), $\mathit{A}\left(t\right)$, $\mathit{B}\left(t\right)$, $\mathit{C}\left(t\right)$, $\mathsf{\Gamma }\left(t\right)$ can be expressed as:

$$\begin{array}{cc}\hfill \mathit{A}\left(t\right)& =\left[\begin{array}{cc}{0}_{3\times 3}& \mathit{J}\left(\psi \right(t\left)\right)\\ 0_{3\times 3}& -{\mathit{M}}^{-1}\mathit{D})\end{array}\right],\mathit{B}\left(t\right)=\left[\begin{array}{c}{0}_{3\times 3}\\ {\mathit{M}}^{-1}\end{array}\right],\\ \hfill \mathsf{\Gamma }\left(t\right)& =\left[\begin{array}{c}{0}_{3\times 3}\\ {\mathit{M}}^{-1}{\mathit{J}}^T\left(\psi \left(t\right)\right)\end{array}\right],\mathit{C}\left(t\right)=\left[\begin{array}{c}{I}_{3\times 3},0_{3\times 3}\end{array}\right],\end{array}$$

(6)

2.2. Discrete-Time NMPC Scheme

In general, the NMPC scheme is designed based on the discrete-time model of the system. Therefore, referring to Zhu et al. [26] and Veerasamy et al. [27], the ship model (5) can be discretized utilizing a sampling time of $T_s$ as:

$$\begin{array}{ccc}\hfill \mathit{x}(\mathit{k}+\mathbf{1})& =& \mathit{A}\left(k\right)\mathit{x}\left(k\right)+\mathit{B}\left(k\right)\mathit{u}\left(k\right)+\mathsf{\Gamma }\left(k\right)\mathit{d}\left(k\right),\hfill \\ \hfill \mathit{z}\left(k\right)& =& \mathit{C}\mathit{x}\left(k\right),\hfill \end{array}$$

(7)

where k is the current time step with regard to the current time t, $\mathit{x}\left(k\right)={[x\left(k\right),y\left(k\right),\psi \left(k\right),u\left(k\right),v\left(k\right),r\left(k\right)]}^T$, $\mathit{u}\left(k\right)={[{\tau }_u\left(k\right),{\tau }_v\left(k\right),{\tau }_r\left(k\right)]}^T$, $\mathit{d}\left(k\right)={[d_1\left(k\right),d_2\left(k\right),d_3\left(k\right)]}^T$, $\mathit{z}\left(k\right)={[x\left(k\right),y\left(k\right),\psi \left(k\right)]}^T$ are corresponding vectors of the continuous-time model (5), and the discrete-time matrices $\mathit{A}\left(k\right),\mathit{B}\left(k\right)$, as well as $\mathsf{\Gamma }\left(k\right)$, can be determined at each time step as:

$$\begin{array}{ccc}\hfill \mathit{A}\left(k\right)& =& e^{A\left(t\right)T_s},\hfill \\ \hfill \mathit{B}\left(k\right)& =& {\mathit{A}}^{-1}\left(t\right)(\mathit{A}\left(k\right)-\mathit{I})\mathit{B}\left(t\right),\hfill \\ \hfill \mathsf{\Gamma }\left(k\right)& =& {\mathit{A}}^{-1}\left(t\right)(\mathit{A}\left(k\right)-\mathit{I})\mathsf{\Gamma }\left(t\right),\hfill \end{array}$$

(8)

For trajectory tracking, the control objective is to operate the ship by thrusters to track the predefined reference trajectory under external disturbances and input saturation constraints. That is, for a given reference trajectory ${\mathit{z}}_r={[x_r,y_r,{\psi }_r]}^T$, we need to determine the control input to make $\underset{x\to \infty }{lim}\Vert \mathit{z}\left(t\right)-{\mathit{z}}_r\left(t\right)\Vert =0$. Thus, the objective function of NMPC can be constructed as:

$$\begin{array}{ccc}\hfill \underset{U\left(k\right)}{min}\mathit{V}\left(k\right)& =& \sum _{n=1}^{N_p}(\Vert {\mathit{z}}_n\left(k\right)-{\mathit{z}}_r\left(k\right){\Vert }_Q^2+\Vert {\mathit{u}}_n\left(k\right){\Vert }_R^2),\hfill \\ \hfill s.t.{\mathit{x}}_{\mathit{n}+\mathbf{1}}\left(\mathit{k}\right)& =& \mathit{A}\left(k\right){\mathit{x}}_n\left(k\right)+\mathit{B}\left(k\right){\mathit{u}}_n\left(k\right)+\mathsf{\Gamma }\left(k\right){\mathit{d}}_n\left(k\right),n=0,\dots ,N_p-1,\hfill \\ \hfill {\mathit{z}}_n\left(k\right)& =& \mathit{C}{\mathit{x}}_n\left(k\right),n=1,\dots ,N_p,\hfill \\ \hfill {\mathbf{\tau }}_{\mathit{min}}& \le & {\mathit{u}}_{\mathit{n}}\left(\mathit{k}\right)\le {\mathbf{\tau }}_{\mathit{max}},\hfill \end{array}$$

(9)

where $\mathit{U}\left(k\right)={[u_0{\left(k\right)}^T,\dots ,u_{N_c-1}{\left(k\right)}^T]}^T$ is the control input sequence in the control horizon $N_c$, and $\mathit{Z}\left(k\right)={[z_0{\left(k\right)}^T,\dots ,z_{N_p}{\left(k\right)}^T]}^T$ is the predictive output sequence in the predictive horizon $N_p$. Since $N_c$ is less than $N_p$, it is assumed that, when $N_c\ge n\ge N_p$, $u_n\left(k\right)=u_{n-1}\left(k\right)$. In addition, $Q=diag(q_1,q_2,q_3)$ and $R=diag(r_1,r_2,r_3)$ are weight matrices that are positive and definite matrices, which proportionally allocate the state vector and the input vector in the objective function (9).

Note that the system is nonlinear; thus, the system matrices are time-varying during the predictive horizon. However, since the navigation state of the ship changes slowly, it can be assumed that $\mathit{A}\left(k\right)$, $\mathit{B}\left(k\right)$, and $\mathsf{\Gamma }\left(k\right)$ remain unchanged during the prediction horizon. Therefore, the predictive outputs can be calculated like a linear model, and the nonlinear optimization problem will be transformed into a relatively simple quadratic optimization problem. Thus, the optimal control sequence ${\mathit{U}}^{*}\left(k\right)={[u_0^{*}{\left(k\right)}^T,\dots ,u_{N_c-1}^{*}{\left(k\right)}^T]}^T$ is obtained by optimizing the objective function (9), and the first value $u_0^{*}\left(k\right)$ is substituted into the process.

The estimated value of disturbance is introduced into the online optimization calculation of NMPC, so as to eliminate the influence of disturbance and realize robust control. Based on Yang et al. [16], the calculation method of interference estimation value is as follows:

$$\begin{array}{ccc}\hfill \widehat{\mathit{d}}\left(t\right)& =& \mathit{p}\left(t\right)+{\mathit{L}}_0\mathit{J}\left(\psi \right)\mathit{M}\mathbf{\nu },\hfill \\ \hfill \dot{\mathit{p}}\left(t\right)& =& -{\mathit{L}}_0\mathit{p}\left(t\right)-{\mathit{L}}_0[\dot{\mathit{J}}\left(\psi \right)\mathit{M}\mathbf{\nu }+\mathit{J}\left(\psi \right)(-\mathit{D}\mathbf{\nu }+\mathbf{\tau })+{\mathit{L}}_0\mathit{J}\left(\psi \right)\mathit{M}\mathbf{\nu }],\hfill \end{array}$$

(10)

where $\widehat{\mathit{d}}\left(t\right)$ is the estimate of $\mathit{d}\left(t\right)$, $\mathit{p}\left(t\right)$ is an auxiliary state vector of the disturbance observer, and ${\mathit{L}}_0$ is a positive definite symmetric matrix. In addition, the predictive output sequence $\mathit{Z}\left(k\right)$ can be obtained from Equations (7) and (8) as:

$$\begin{array}{c}\hfill \mathit{Z}\left(k\right)=\mathsf{\Omega }\mathit{x}\left(k\right)+\mathsf{\Phi }\mathit{U}\left(k\right)+\mathsf{\Lambda }{\mathit{O}}_d\left(k\right),\end{array}$$

(11)

where ${\mathit{O}}_d\left(k\right)=\mathsf{\Gamma }\left(k\right)\widehat{\mathit{d}}\left(k\right)$, and:

$$\begin{array}{ccc}\hfill \mathsf{\Gamma }& =& \left[\begin{array}{cc}& \mathit{C}\mathit{A}\left(k\right)\\ & \mathit{C}{\mathit{A}}^2\left(k\right)\\ & ⋮\\ & \mathit{C}{\mathit{A}}^{N_p}\left(k\right)\end{array}\right],\mathsf{\Omega }=\left[\begin{array}{cc}& \mathbf{C}\\ & \mathit{C}(\mathit{A}\left(k\right)+{\mathit{I}}_{6\times 6})\\ & ⋮\\ & \mathit{C}\left(\sum _{i=1}^{N_p-1}{\mathit{A}}^{i}k+{\mathit{I}}_{6\times 6}\right)\end{array}\right],\hfill \\ \hfill \mathsf{\Phi }& =& \left[\begin{array}{cccc}\mathit{C}\mathit{A}\left(k\right)& 0& \cdots & 0\\ \mathit{C}\mathit{A}\left(k\right)\mathit{B}\left(k\right)& \mathit{C}\mathit{B}\left(k\right)& \cdots & 0\\ ⋮& \ddots & \cdots & ⋮\\ \mathit{C}{\mathit{A}}^{N_c-1}\left(k\right)\mathit{B}\left(k\right)& \mathit{C}{\mathit{A}}^{N_c-2}\left(k\right)\mathit{B}\left(k\right)& \cdots & \mathit{C}\mathit{B}\left(k\right)\\ \mathit{C}{\mathit{A}}^{N_c}\left(k\right)\mathit{B}\left(k\right)& \mathit{C}{\mathit{A}}^{N_c-1}\left(k\right)\mathit{B}\left(k\right)& \cdots & \mathit{C}(\mathit{A}\left(k\right)+{\mathit{I}}_{6\times 6})\mathit{B}\left(k\right)\\ ⋮& \ddots & \cdots & ⋮\\ \mathit{C}{\mathit{A}}^{N_p-1}\left(k\right)\mathit{B}\left(k\right)& \mathit{C}{\mathit{A}}^{N_p-2}\left(k\right)\mathit{B}\left(k\right)& \cdots & \mathit{C}\left(\sum _{i=1}^{N_p-N_c}{\mathit{A}}^{i}k+{\mathit{I}}_{6\times 6}\right)\mathit{B}\left(k\right)\end{array}\right],\hfill \end{array}$$

(12)

3. Genetic Algorithm-Based Weight Tuning of NMPC

The weight coefficients in objective function (9) are very important for the NMPC online optimization calculation. If not selected properly, it will lead to a non-optimal solution, or even non-convergence. However, the weight matrices Q and R are usually determined by trial or experience, which is usually time-consuming and cannot be adjusted in real time with the change of control state to maintain the optimal solution [28]. Therefore, it is necessary to adjust the weight matrices Q and R online in real time to improve the control accuracy of NMPC. In this way, the adjusted parameters will have the direct physical significance of making the optimization process faster and the results more accurate [29]. Considering that the genetic algorithm (GA) has the characteristics of strong adaptability and automatic search [30], it is selected to optimize the weight matrices Q and R.

3.1. General Principle of GA

GA is a method of searching an optimal solution by simulating the process of natural evolution. It transforms the problem-solving process into a process that is similar to the crossover and mutation of chromosome genes in biological evolution. The main components of GA are fitness function, selection operator, crossover operator, and mutation operator, etc. Fitness function is an index that is used to judge the quality of individuals in the group. It is evaluated according to the objective function of the problem. The function of the selection operator is to screen out individuals with excellent genes for follow-up. The crossover operator refers to the operation of replacing and reorganizing some structures of two parent individuals to generate new individuals. The content of the mutation operator is to change the gene value on some loci of the individual string in the population.

Generally, when using GA to solve an optimization problem, a population composed of multiple individuals should be constructed first. Each chromosome of individuals in the population represents a feasible solution to the problem. Additionally, each solution is substituted into the fitness function to evaluate the fitness of the individual. The solutions with low fitness are gradually eliminated, and those with high fitness are added. Then, the next generation is generated through replication, crossover, mutation, and other operations. In this way, individuals with high fitness will be obtained after n generations of evolution; thus, a relatively accurate solution to the problem can be obtained.

Herein, to obtain online the optimal weight factors of the weight matrices Q and R in the NMPC objective function, the weight factors $q_1,q_2,q_3,r_1,r_2,r_3$ can be continuously optimized in GA. To improve the solution efficiency, a novel GA algorithm is proposed, considering the actual needs of ship track control, and utilizing the improved crossover operator, mutation operator, crossover rate, and mutation rate.

3.2. Application and Improvement of GA for Real-Time Weight Tuning

3.2.1. Fitness Function

Fitness function determines the direction of GA evolution, and is set according to practical problems. In NMPC, the best choice of Q and R is required to minimize the objective function (9) for the best trajectory tracking effect. Therefore, the Q and R optimized by the fitness function in this novel GA needs to meet the requirements of minimizing the objective function, and the fitness function is defined as follows:

$$\begin{array}{c}\hfill Fitness=\frac{1}{\mathit{v}\left(k\right)}=\frac{1}{{\displaystyle \sum _{n=1}^{N_p}}(\Vert {\mathit{z}}_n\left(k\right)-{\mathit{z}}_r\left(k\right){\Vert }_Q^2+\Vert {\mathit{u}}_n\left(k\right){\Vert }_R^2)},\end{array}$$

(13)

3.2.2. Improvement of Crossover Operator

In order to ensure that the genes of better individuals can be inherited to the offspring, the tournament selection strategy is adopted for the genetic operator. The selected individuals will generate a new individual through the crossover operator, and the standard crossover algorithm can be expressed as follows:

$$\begin{array}{c}\hfill {\mathit{Q}}_{\mathbf{1}}=({\mathit{T}}_{\mathbf{1}}+{\mathit{T}}_{\mathbf{2}})/2+\beta \times ({\mathit{T}}_{\mathbf{1}}-{\mathit{T}}_{\mathbf{2}})/2,\end{array}$$

(14)

where ${\mathit{Q}}_{\mathbf{1}}$ is the offspring chromosome; ${\mathit{T}}_{\mathbf{1}}$ and ${\mathit{T}}_{\mathbf{2}}$ are parent chromosomes. $\beta$ is a random number between [0,1]. However, there is a certain chance that the tournament strategy will repeatedly choose parent chromosomes, which indirectly leads to insufficient diversity. In order to solve this problem, by introducing ${\mathit{T}}_{\mathbf{3}}$, the improved crossover strategy is as follows:

$$\begin{array}{c}\hfill {\mathit{Q}}_{\mathbf{2}}=({\mathit{T}}_{\mathbf{1}}+{\mathit{T}}_{\mathbf{2}}+{\mathit{T}}_{\mathbf{3}})/3+\beta \times (2{\mathit{T}}_{\mathbf{1}}-{\mathit{T}}_{\mathbf{2}}-{\mathit{T}}_{\mathbf{3}})/3,\end{array}$$

(15)

Comparing Equations (14) and (15), the former is conducive to individual stability and less individual change range, due to less utilization of parent information. In Equation (15), each individual uses more information from the surrounding individuals, so that it is easy to quickly search a certain area in the early stage of the algorithm. Therefore, if the current number of iterations is less than $T\times Ge{n}_{max}$, the crossover algorithm of Equation (15) is adopted. Otherwise, Formula (14) is adopted. The specific T and $Ge{n}_{max}$ will be determined in subsequent experiments.

3.2.3. Adoption of Mutation Operator

Aiming to solve the problems of mutation operator, an improvement is made on the basis of the method proposed by Tsoulos et al. [31]. Firstly, chromosomes in the population are sorted according to their fitness values from high to low, and then the mutation scale is determined according to the mutation rate $P_m$. In order to purposefully mutate the chromosome with low fitness rather than the optimal one, the mutation operator equation calculated using the PSO optimization method is as follows:

$$\begin{array}{c}\hfill c_i^{{}^{\prime }}=l_1{m}_1(c_i^b-c_i)+l_2{m}_2(c_i^b-c_i),\end{array}$$

(16)

where $l_1$ and $l_2$ are two positive constants, $m_1$ and $m_2$ are uniform random numbers between [0,1], and $c^b$ represents the chromosome with largest fitness currently found. The improvement can ensure that only the worst solutions will be mutated. This will improve the gene, increase the diversity of the population, and prevent the algorithm from falling into a local optimal solution.

3.2.4. Improvement of Crossover and Mutation Rate

According to the Darwinian evolution principle, in the early stage of genetic algorithm evolution, the crossover operator can enhance the global search ability of the algorithm, which can avoid falling into the local optimal solution. In the later stage of evolution, when getting closer to optimal solution, the crossover operator will affect the convergence speed. At this time, the mutation operator plays a more important role. Herein, the crossover rate ${\mathit{P}}_{\mathit{c}}$ and the mutation rate ${\mathit{P}}_{\mathit{m}}$ can be defined as follows:

$$\begin{array}{c}\hfill {\mathit{P}}_c={\mathit{P}}_1-0.35\times \left(\frac{gen}{Ge{n}_{max}}\right),\\ \hfill {\mathit{P}}_m={\mathit{P}}_2+0.05\times \left(\frac{gen}{Ge{n}_{max}}\right),\end{array}$$

(17)

where ${\mathit{P}}_1$ is the initial crossover rate and ${\mathit{P}}_2$ is the initial mutation rate. $Ge{n}_{max}$ is the maximum number of iterations and $gen$ is the current number of iterations. With the progress of the search process, the crossover rate is getting smaller and smaller, while the mutation rate is getting larger and larger.

3.3. Optimization Procedure of GA-NMPC

The genetic algorithm-based nonlinear model predictive control (GA-NMPC) for vessel trajectory tracking is shown in Figure 1. The weight matrices Q and R of NMPC are optimized in real time using the genetic algorithm, and the implementation steps are concluded as follows.

(1)

Initialize the weight matrices Q and R, and obtain $\mathit{U}\left(k\right)$ and $\mathit{Z}\left(k\right)$ through Equations (9) and (10), respectively.

(2)

Define parent population size as $N_{pop}$, and $q_1,q_2,q_3,r_1,r_2,r_3$ are the genes carried by each individual’s chromosome.

(3)

Encode chromosomes and calculate the fitness values through Equation (13), and sort the chromosomes from high to low according to fitness values.

(4)

Determine the number of chromosomes that need to be crossed according to crossover rate $P_c$. Select a corresponding number of chromosomes with high fitness values through the tournament strategy, and then cross these chromosomes according to Equations (14) and (15) to obtain new chromosomes. Finally, the newly generated chromosomes are used to replace the chromosomes with the lowest fitness values in the parent population.

(5)

After crossover, the chromosome will be mutated according to the mutation rate $P_m$. The individuals with low fitness values in the parent population are mutated based on Equation (16). Thus far, a new parent population has emerged.

(6)

If the number of iterations does not reach the maximum number of iterations $Ge{n}_{max}$, continue the cycle. Otherwise, the individual with largest fitness value is introduced into objective function (9) as weight matrices Q and R. The optimal solution sequence ${\mathit{U}}^{*}\left(k\right)={[u_0^{*}{\left(k\right)}^T,\dots ,u_{N_c-1}^{*}{\left(k\right)}^T]}^T$ is obtained by solving the objective function (9), and then the first element $u_0^{*}\left(k\right)$ is input into the closed-loop system to complete the ship trajectory tracking.

4. Simulation Evaluation and Discussion

In this section, the simulation evaluation of the trajectory tracking performance using the proposed GA-NMPC is implemented. The CSII ship, which is a replica of the supply ship with a ratio of 1:70, is chosen as the simulation object. Its a priori identified dynamic parameters are shown in

Table 1. CSII ship parameters [32].

Parameter	Value	Unit	Parameter	Value	Unit
m	23.8	kg	$Y_r$	–0.1	kg· m/s
$I_z$	1.76	kg· m${}^2$	$Y_{\dot{v}}$	–10.0	kg
$X_g$	0.046	m	$Y_{\dot{r}}$	–0.0	kg· m
$X_u$	–2	kg/s	$N_v$	–1.0	kg· m/s
$X_{\dot{u}}$	–2.0	kg	$N_r$	–0.5	kg· m${}^2$/s
$Y_v$	–0.7	kg/s	$N_{\dot{r}}$	–1.0	kg· m${}^2$

.

4.1. Main Parameters Setting

The maximal and minimal thrust limit in (4) is set as ${\mathbf{\tau }}_{max}=$ [2.0N, 2.0N, 1.5N·m]${}^T$, ${\mathbf{\tau }}_{min}=$ [–2.0N, –2.0N, –1.5N·m]${}^T$, respectively. The sample time step is set as $T_s=0.1$. In GA-NMPC, the population size $N_{pop}=30$, and the initial crossover rate and mutation rate are $P_c=0.9$, $P_m=0.2$, respectively. The parameters defined in Section 3.2.2 are set as $T=0.6$, $Ge{n}_{max}=100$.

To demonstrate the effectiveness of the proposed real-time weight optimization algorithm for NMPC, the control effects of NMPC with different Q and R will be compared. Herein, three cases are simulated, which are:

(1)

The proposed GA-NMPC algorithm with real-time online weight optimization with initial $Q=diag(50,50,50)$, $R=diag(1,1,1)$;

(2)

NMPC with constant weight matrices $Q=diag(100,100,100),R=diag(1,1,1)$;

(3)

NMPC with constant weight matrices $Q=diag(5,5,5),R=diag(1,1,1)$.

In addition, two different reference trajectories are used to verify the performance, and they are set as follows:

(1)

Circle trajectory, with the expression described as:

$$\begin{array}{ccc}\hfill x_r& =& 3-3cos\left(0.025t\right),\hfill \\ \hfill y_r& =& 4sin\left(0.025t\right),t\ge 0,\hfill \end{array}$$

(18)

(2)

Sine curve trajectory, with the expression described as:

$$\begin{array}{ccc}\hfill x_r& =& 0.05t,\hfill \\ \hfill y_r& =& 3sin\left(0.025t\right),t\ge 0.\hfill \end{array}$$

(19)

The desired heading ${\psi }_r$ of the two trajectories are both taken as the direction of the tangent vector along the path.

4.2. Results Analysis

The simulation results of the three cases based on the two trajectories are shown in Figure 2, Figure 3, Figure 4 and Figure 5, where the results of cases (1)–(3) are represented by red, blue, and pink solid lines, respectively, while the reference trajectory is represented by green dotted lines. Figure 2a and Figure 4a are ship position and heading. Figure 2b and Figure 4b represent ship velocities in surge, sway, and yaw. Figure 3 and Figure 5 show the trajectory tracking results under different cases. For the circle trajectory, it can be seen from Figure 2a and Figure 3 that the trajectory tracking results of the three cases are ideal. However, by observing the local enlarged picture, we can see that the error of the result obtained by the controller using constant weight matrices is slightly larger than that of GA-NMPC. In addition, the heading angle of NMPC with constant matrices Q and R cannot return to zero in time after completing a circle, so there will be obvious fluctuations when starting a new voyage. This problem can also be reflected in Figure 2b that, at the end of a tracking cycle, the ship velocities will fluctuate significantly. However, the adoption of real-time-optimized Q and R solves this problem well.

Figure 4 and Figure 5 show the tracking results of the sine curve trajectory. In Figure 4a and Figure 5, the three cases have a good tracking effect on the trajectory. However, in contrast, the control result of NMPC with constant Q and R is slightly unstable compared with that of real-time-optimized Q and R. Additionally, when the element values in matrix Q are small, as can be seen from Figure 4b, the velocity $\nu$ fluctuates greatly. Therefore, by receiving ship feedback to optimize the weight matrices in real time, the tracking error can be reduced and obtain a better control effect.

Figure 2c,d and Figure 4c,d respectively show the matrices Q and R optimized in real time by GA-NMPC during ship tracking of the two trajectories. It can be seen that, at every step, the genetic algorithm adjusts elements in the Q and R according to the tracking error to improve population fitness, which can achieve the best control effect. Although using constant weight matrices can also obtain a relatively ideal result, it needs many attempts and adjustments, which is time-consuming and laborious. In addition, the final control effect may not be the best, due to the invariance of weights.

5. Conclusions

In this paper, a real-time weight tuning method based on the genetic algorithm (GA) is proposed and applied in a ship nonlinear model predictive controller for ship trajectory tracking. The weight matrices Q and R to be optimized in the NMPC objective function are taken as individuals in GA. The individuals with the highest fitness are obtained as weight factors in the matrices of Q and R; thus, the NMPC will solve the objective function online and obtain the optimal control action. In addition, the crossover operator, mutation operator, crossover rate, and mutation rate of the genetic algorithm are improved to enhance the diversity of the population and the efficiency of searching for the optimal solution. This optimization method adjusts the weight proportion of the input value in the objective function in real time according to the actual navigation and makes up for the defects of the traditional trial and error method, as well as the empirical method. The controller optimized by this method is simulated and compared with the traditional NMPC. The results show that the proposed GA-NMPC method makes it easy to implement NMPC and improves the control efficiency of ship trajectory control.