309

1 INTRODUCTION

As far as the propulsion systems of marine vessels are

concerned, a number of thruster types have been

developed up to date for the requirements of the

maritime industry, but when it comes to their

application in dynamic positioning systems, the non-

retractable and retractable azimuth thrusters, in

eventual combination with tunnel thrusters, are

mostly used. This choice is quite understandable, and

is justified by the fact that due to the different

directions of external disturbances, thrusters should

be able to operate in all 360° at all times. Somewhat

less common application is the Voith Schneider

cycloid propulsor and combinations that include main

propeller(s) with rudder(s).

The orientation, i.e. the angle or azimuth of each

thruster, as well as the required thrust it generates, is

determined by the control logic of dynamic

positioning systems. This whole process is called the

thrust allocation and represents a very complex

mapping of the previously calculated or estimated

environmental forces and moment in the set of

referent states of the available thrusters.

Since dynamically positioned vessels usually have

fixed and azimuth thrusters, the vector of design

variables must have one variable for fixed (e.g.

tunnel) thrusters and two for azimuth thrusters. The

basic constraints on the objective function are (at

least) three equalities stemming from the fact that the

generated thrust forces in surge and sway and

moment in yaw should be equal with environmental

loads for all three horizontal degrees of freedom. This

simplified optimization approach can be reduced on

finding the conditional extremes of the Lagrange's

objective function. If additional constraints, usually

Derivative Free Optimal Thrust Allocation in Ship

Dynamic Positioning Based on Direct Search Algorithms

. Valčić

University of Rijeka, Rijeka, Croatia

. Prpić-Oršić

University of Rijeka, Rijeka, Croatia

Centre for Autonomous Marine Operations and Systems (AMOS), NTNU, Trondheim, Norw

ABSTRACT: In dynamic positioning systems, nonlinear cost functions, as well as nonlinear equality and

inequality constraints within optimal thrust allocation procedures cannot be handled directly by means of the

solvers like industry-standardized quadratic programing (QP), at least not without appropriate linearization

technique applied, which can be computationally very expensive. Thus, if optimization requirements are strict,

and problem should be solved for nonlinear objective function with nonlinear

equality and inequality

constraints, than one should use some appropriate nonlinear optimization technique. The current state-of-the-

art in nonlinear optimization for gradient-based algorithms is surely the sequential quadratic programing

(SQP), both for general applications and specific thrust allocation problems. On the other hand, in recent time,

one can also notice the increased applications of gradient-free optimization methods in various engineering

problems. In this context, the implementation of selected derivative free direct search algorithms in optimal

thrust allocation is proposed and discussed in this paper, and avenues for future research are provided.

http://www.transnav.eu

the

International Journal

on Marine Navigation

and Safety of Sea

Transportation

Volume 14

Number 2

June 2020

DOI:

10.12716/1001.14.02.05

310

expressed by the matrix inequalities (e.g. thruster

saturation, thruster efficiency, electrical power

limitation, etc.), are added to the basic constraints, the

optimization task becomes considerably more

complex and is usually solved by the quadratic

programming algorithms or so-called QP solvers

(Jenssen and Realfsen, 2006).

From a theoretical point of view, the problem of

the thrust allocation could be solved by linear

programming (LP solvers), but due to the

approximation of the relation between the electrical

power consumption (kW) and the generated thrust

(kN) by the quadratic function, some of the variants

of the quadratic programming are usually used

(Snijders, 2005; Wit, 2009).

If the problem of the quadratic programming of

the thrust allocation is set correctly, it can be explicitly

solved, i.e. it is possible to determine the global

minimum (Leavitt, 2008). In general, the QP consists

of the quadratic objective function and linear

equalities and inequalities representing the

conditions, i.e. the constraints. In addition to this,

Sørdalen (1997) has shown that the constraints on

azimuth thrusters can lead to singular configurations,

which he solved using the method of singular values

decomposition (SVD). This approach provided

significantly lower power consumption, effectively

eliminated the issue of the forbidden zones, reduced

tear and wear of thrusters. With the application of so-

called logical inequalities and Moore-Penrose pseudo-

inverse matrix (SVD method), it is possible to directly

determine the vector of demanded forces and moment

(Gierusz and Tomera, 2006; Yang et al., 2011b).

If the constraints in the quadratic optimization task

become nonlinear, it is no longer possible to use the

QP solvers directly. One of the possible solutions to

this problem is the application of the so-called

sequential quadratic programming (SQP) technique

that is generally used to minimize an arbitrarily

selected objective function regarding the nonlinear

constraint set in the form of equalities and

inequalities. The possible applications of SQP

approach in optimum thrust allocation were

investigated by Liang and Cheng (2004) and Johansen

et al. (2004). Although tested only on simulation

models, the obtained results (Liang and Cheng, 2004)

indicate very good capabilities of the SQP solver

which in a computational sense can execute the

allocation very fast with a small thrusters' azimuth

change. Johansen et al. (2004) have further expanded

the application possibilities of the SQP approach with

the emphasis on avoiding possible singularities that

are unacceptable in control sense.

In addition to the SQP approach for solving the

problem of nonlinear constraints of the optimization

task, most recently the genetic algorithms (GA) have

been increasingly used as a robust solution that

ensures a good convergence of the global

optimization process (Yang et al., 2011a; Zhao et al.,

2010). The tests that have been carried out by Yang et

al. (2011a) indicate the promising results on using

these algorithms, although the authors point out the

problem of possible application of GA in thrust

allocation regarding the slow convergence.

In order to recap, one should notice that the

current state-of-the-art in nonlinear optimization for

gradient-based algorithms is surely the sequential

quadratic programing (SQP), both for general

applications and specific thrust allocation problems.

In comparison with e.g. Lagrangian multiplier

method (LMM) or pure QP algorithms, which are

both appropriate solutions for optimization problems

with linear equality and linear inequality constraints,

SQP approach is superior when dealing with

problems that have significant nonlinearities within

their constraints.

On the other hand, and in comparison with the

gradient-based optimization methods, derivative free

optimization methods usually does not need any

particular information about the gradient or Hessian

matrix of the objective function. Moreover, derivative

free methods can be applied even for objective

functions that are not continues nor differentiable,

which makes them particularly convenient in cases

when the objective function is not explicitly defined,

when evaluation of the objective function and/or its

derivatives is too much time consuming and thus not

acceptable, when objective function is noisy and

derivatives or finite difference approximations are not

reliable nor acceptable for further analysis, etc.

Although the field of derivative free optimization

is usually extended, or at least coupled with the so-

called black-box optimization methods, the focus in

this work is placed only to the family of the so-called

direct search (DS), or pattern search (PS) algorithms.

The main reason for this choice is that direct search

algorithms are much better supported and have very

detailed literature background on rigorous

convergence analysis (Audet and Hare, 2017; Conn et

al., 1997, 1991; Torczon, 1997). Therefore, the

applicability and implementation issues of selected

derivative free direct search algorithms in optimal

thrust allocation problems have been analysed and

discussed in this paper, and avenues for future

research are emphasized as well.

2 METHODOLOGY

2.1 General considerations on direct search algorithms

Direct or pattern search algorithms are based on a

common idea by which a sequence of points is

determined with the property that in each successive

point the value of objective function decreases. As

already mentioned, this sequence of points, which

defines directions from one point to another, is not

calculated by means of function gradient, but is rather

based on a set of points around the current point, in

which the objective function is evaluated. These

surrounding points are determined by polling and

thus they create a so-called poll set that presents a

mesh, i.e. all possible vector directions by which one

can shift from the current point to any other point

from the poll set. If some point within the poll set is

found that sufficiently decreases the objective

function at the current step, than that point becomes a

new current point for the next iteration. Otherwise,

the mesh should be redefined so the algorithm could

try to find a new direction on a smaller scale. In

general, one can differ three main direct search

algorithms as follows:

− generalized pattern search (GPS) algorithm,

311

− generating set search (GSS) algorithm,

− mesh adaptive direct search (MADS) algorithm.

Besides the all other significant properties,

generally speaking, the main difference between GPS,

GSS and MADS is the number of directions from the

current point to any other point from the poll set, as

well as direction geometrical characteristics. Other

important properties and differences are mostly

related to handling of linear and nonlinear

constraints.

2.2 Generating set search algorithm with augmented

Lagrangian

The generating set search (GSS) algorithm is very

similar to GPS, particularly for the problems of

unconstrained optimization in which their patterns

are identical. The main difference between GPS and

GSS is related to constrained optimization problems.

GSS, as an extension of GPS, is well suited for bounds

and linear constraints, where directions in positive

spanning set



are determined using the nearby

active constraints from the working set (Kolda et al.,

2006, 2003). In other words, GSS is more efficient in

comparison with GPS for linearly constrained

optimization problems.

When it comes to nonlinear constraints, GPS is not

well suited, but with implementation of augmented

Lagrangian within the GSS algorithm, which was

introduced by Kolda et al. (2006), GSS can handle

optimization problems with both linear and nonlinear

constraints. However, this approach has been

analysed under the assumption that the objective

function and constraints should be twice continuously

differentiable, which is typical required property for

gradient-based methods.

The augmented Lagrangian pattern search (ALPS)

algorithm is primarily used for solving optimization

problems with nonlinear equality and inequality

constraints, which means that bounds and linear

constraints are handled differently, usually with

nearby active constraints strategy (Kolda et al., 2006,

2003). For some general optimization problem with

nonlinear equality and inequality constraints of the

following form

min ( )

s.t. ( ) 0,

( ) 0,

∈





= ∈





≥∈





(1)

that is considered to be solved by some pattern search

algorithm, associated sub-problem based on

augmented Lagrangian should be formed as follows

( , , , ) ( ) log( ( ))

() ()

ii i i

ii i

iE iE

∈

∈∈

=−+

∑

∑∑

x λ s ρ x λs s c x

λc x c x

(2)

where

0≥

are Lagrangian multipliers,

0≥

are slack variables, and

is positive penalty

parameter. One should notice that each sub-problem

(2) presents one iterative step, which makes this

approach with nonlinear constraints highly

computationally expensive. During each iteration,

values of

,λ

and

are kept constant, until the

sub-problem (2) is minimized, whereupon are all

updated. Otherwise, the penalty parameter

increased and a new sub-problem is formed. These

steps are repeated until the termination, which is

based on some predefined stopping criteria.

2.3 Mesh adaptive direct search algorithm

The MADS algorithm, which was introduced by

Audet and Dennis (2006), primarily as direct search

algorithm for solving constrained optimization

problems of the general form

min ( )

∈Ω

(3)

does not require any assumptions related to the

smoothness of objective function nor to constraints

that could be either linear, or nonlinear, or both. If

RΩ=

in (3), than the previous optimization

problem becomes unconstrained.

General constraints with MADS algorithm are

usually handled by the so-called extreme barrier

strategy (Audet and Hare, 2017), which is based on

extreme barrier function

: {}

fR R

Ω

→ ∪∞

defined as

( ), if

, if .

Ω

∈Ω





∞ ∉Ω



(4)

The associated principle is very simple and is

based on the fact that optimization is performed using

the barrier function

()

Ω

f x

as the objective, rather

than

( ).f x

More advanced approaches can take

into account two-phase extreme barrier strategy, filter

methods (Audet and Dennis, 2004; Dennis Jr. et al.,

2004), progressive barrier strategy (Audet and Dennis,

2009; Le Digabel, 2011) or mixture between extreme

barrier and progressive barrier called progressive-to-

extreme barrier strategy (Audet et al., 2010).

However, the MADS algorithm is primarily

orientated to inequality constraints with bounds,

which means that equality constraints could be

challenging. For this purpose, one can substitute one

equality constraint with two equivalent inequality

constraints, although this approach could be

cumbersome in some optimization problems and

algorithm efficiency could be questionable or even not

acceptable. This issue is also related to the complexity

of equality constraints and to the selection of initial

point

Possible alternatives for handling this issues could

be closely related to approaches introduced with GPS

and GSS algorithms, i.e. to equality constraint

handling by means of the nearby active constraints or

augmented Lagrangian method (Kolda et al., 2006,

2003). Recent research directions are also aimed

towards the combining of gradient-based methods,

like sequential quadratic programing, and derivative

312

free optimization with equality constraints (Tröltzsch,

2016).

3 OPTIMAL THRUST ALLOCATION

3.1 Problem definition

From the optimization point of view, thrust allocation

problem usually comes down to the minimization of

total power consumption or some other appropriate

objective function

:,→ fX

in terms of thrust

Hence, this optimization task can be defined as

follows

* arg min( )

s.t. ( ) 0

() 0

∈







≤



≤≤



l xu

(5)

where

∈ Xx

and

{ | ( ) 0,

X =∈=x hx

1,2,..., ; ( ) 0, 1,2,..., }

ip jq= ≤=gx

presents a set

of feasible thrust region that depends on equality and

inequality constraints

() [(), ()... ()],

hh h=

hx x x x

and

() [ (), (),... ()],

gg g

gx x x x

respectively,

and

are lower and upper boundaries of

and

power function

is commonly assumed to be

twice-continuously differentiable, i.e. sufficiently

smooth.

Therefore, nonlinear optimization problem (5) for

thrust allocation can be redefined in the following

form

,max

eq eq

FZ FZ

22 2

,max

* arg min ( )

()

s.t. 0

ik ik

ix iy

iii

ix iy i

uuT

∈





= +











−=





≤



+≤



∑

b Au

Au0

(6)

where objective function is total delivered power for r

thrusters (Leavitt, 2008; Wit, 2009), defined in terms

of individual thrust components

and

resultant individual thrust

of each thruster in

body reference frame {b},

,maxi

and

,maxi

indicate

maximum power and maximum thrust for any i-th

thruster, respectively,

1 2,<≤

and u presents the

space of thruster states, which is for r thrusters

defined as

11 2 2

[ , , , ,..., , ,..., , ] .

x y x y ix iy rx ry

uuu u uu uu= ∈u

(7)

Matrices

and

in (6) are defined as

eq eq

[, ]=A BC

and

eq eq

[, ],

= −

b τ0

which take

into account thrust allocation problem

= −

Bu τ

and

additional equality constraints for tunnel thrusters in

form of

eq eq

,=Cu 0

if there are any. Matrix

∈B

is a well-known configuration matrix and

although

is usually a control vector calculated by

a DP controller, in this quasi- static analysis it

presents a vector of environmental loads in the

horizontal plane that are calculated on the basis of

model tests according to the usual design

recommendations. Matrices

and

are

defined as

FS FS

FZ FZ

22 2

start start

FS FS

end end

cos sin

, ,

cos sin

ik ik

ϕϕ



−



= ∈=∈



−



(8)

and are basis for modelling forbidden zones in terms

of circle sectors bounded with two radii at angles

start

and

end

where k indicates what feasible set

is selected according to some predefined

criteria. Final inequality equations

22 2

,max

+≤

ix iy i

uuT

in (6) are related to saturation of thrusters and in this

form they also present nonlinear thrust regions for

each thruster.

Alternative approach, for

1.5,

in which the

relationship between delivered power and generated

thrust is based on thrust and torque coefficients

and

is very similar to (6). The only difference is

in the form of the objective function, while all other

constraints are the same. In this case, objective

function of nonlinear optimization problem can be

defined as

0, 0,

* arg min ( )

ix iy

i Ti Ti

DK K

∈





= +





∑

(9)

where

0,Ti

and

0,Qi

are thrust and torque

coefficients at bollard pull conditions, respectively,

is propeller diameter,

is (sea) water density,

and

1,2,..., .=

3.2 Numerical example and analysis of results

Straight applications of direct search algorithms in

optimal thrust allocation do not require any

additional transformation of optimization tasks. In

other words, optimization tasks like (6) or (9) are

already fully prepared in order to be solved by means

of pattern search algorithms, which is very

convenient, particularly if one wants to perform

appropriate comparisons between these algorithms

and any other algorithm that could be of interest, like

QP, SQP, etc.

313

Table 1. Thruster configuration with basic data

__________________________________________________________________________________________________

# Thruster

(m)

/PD

-1

(min )

max

(kN)T

min

(kN)T

max

(kW)P

__________________________________________________________________________________________________

1 Tunnel 82.0 0.0 2.0 1.2 330.0 165.0 -165.0 ±1200.0

2 Azimuth 57.0 4.5 2.5 1.2 900.0 390.0 0.0 2400.0

3 Azimuth 52.0 -4.5 2.5 1.2 900.0 390.0 0.0 2400.0

4 Azimuth 28.0 -15.0 2.5 1.2 900.0 390.0 0.0 2400.0

5 Azimuth -22.0 15.0 2.5 1.2 900.0 390.0 0.0 2400.0

6 Azimuth -60.0 15.0 3.6 1.2 630.0 760.0 0.0 4500.0

7 Azimuth -60.0 -15.0 3.6 1.2 630.0 760.0 0.0 4500.0

__________________________________________________________________________________________________

For the purpose of this paper, a simple example is

provided for heavy lift DP vessel Saipem 3000, which

was selected as a reference vessel with the length

overall of

162.0 m,

L =

beam

38.0 mB =

displacement

24000 t,∆=

and is equipped with

seven thrusters, one of which is a bow tunnel thruster

and the rest are azimuth thrusters. Their approximate

positions on the hull regarding the body reference

frame are given in Table 1, together with basic

thruster data.

In order to illustrate the problem, only the

allocation results with environmental loads at wind

speed of

10 m/s,

wind

significant wave height

3.21 m,=

wave peak period

8.41 s=T

and

sea current velocity

0.5 m/s=

are presented. The

angle of attack was the same for all disturbances at

any time and varied from

= °

360°

with the

rate of

10 .°

Environmental loads were calculated on

the basis of model tests according to usual design

recommendations. Finally, the vector

,, ,

[ ],

x loads y loads z loads

FFN

which represents the

total environment load for some angle of attack

was calculated quasi-statically in order to be used in

(6).

The optimal thrust allocation is performed by

MADS algorithm with augmented Lagrangian and

obtained results with MADS are compared with

results obtained by SQP. For this analysis, results of

SQP algorithm served as reference values, and

comparisons were made in terms of average RMSE of

optimal solutions between MADS and SQP, as well as

in terms of average convergence time for these two

approaches. The target PC configuration was based

on Intel(R) Core (TM) i5-7500 CPU @ 3.40 GHz, 16 GB

RAM, x64-based processor, 64-bit operating systems

(MS Windows 10). Optimization was performed using

optimization task (6) for

1.5=

and with

MathWorks MATLAB R2017b as a support software.

In order to additionally simplify this analysis,

forbidden zones were not included, so the problem of

non-convex thrust regions could be omitted.

Hence, after all 36 optimization tasks had been

solved, the average convergence time for MADS was

0.587459 s, while for SQP average convergence time

was 0.020125 s. RMSE between optimal solutions

for MADS and SQP was equal to 4.0748·10

-4

. These

results clearly indicate that there are no significant

differences between optimal solutions obtained with

MADS and SQP, but one can notice that SQP is

relatively much faster. However, the reason for this is

also in relative simplicity of optimization task (6),

particularly in this case, i.e. without forbidden zones

included and without some additional nonlinear

constraints. Moreover, objective function in this case

is also relatively simple and convex, which presents

favourable conditions for gradient-based algorithm

like SQP. Nonetheless, MADS showed overall very

promising results, especially in comparison with

other derivative free algorithms like genetic

algorithms (GA) for which convergence time is

usually the biggest issue (Yang et al., 2011a; Zhao et

al., 2010).

4 CONCLUSION

Although the results of optimal thrust allocation

problem obtained with direct search algorithm are

more than satisfactory, particularly in comparison

with state-of-the-art algorithm like SQP, it should be

pointed out that the gradient-based algorithms could

and probably should be a better choice for

optimization problems where the gradient and/or

Hessian of the objective function is known or at least

it can be obtained in sufficient amount of time.

However, this will probably be true for most

unconstrained optimization problems, but on the

other hand, when a large number of nonlinear

constraints is added into the optimization task, and

even more when the objective function is constantly

changing due to some external disturbances, then the

calculation of associated Lagrangian and its gradient,

Hessian or Hessian approximation could be very

demanding within gradient-based methods. Thus, in

these cases direct search algorithms could present a

good or even better alternative. In order to better

evaluate the possibilities of direct search algorithms in

optimal thrust allocation, future analyses and

comparisons should take into account forbidden

zones, i.e. non-convex thrust regions, thrust loss

effects and complete environmental envelope.

Moreover, additional procedures in order to enable

faster convergence of direct search algorithms should

be identified as well.

ACKNOWLEDGEMENTS

This work has been fully supported by the Croatian Science

Foundation under the project IP-2018-01-3739 and by the

University of Rijeka under the project numbers uniri-tehnic-

18-18 and uniri-tehnic-18-266.

REFERENCES

Audet, C., Dennis, J.E., 2009. A Progressive Barrier for

Derivative-Free Nonlinear Programming. SIAM J.

Optim. 20, 445–472.

Audet, C., Dennis, J.E., 2006. Mesh Adaptive Direct Search

Algorithms for Constrained Optimization. SIAM J.

Optim. 17, 188–217.

314

Audet, C., Dennis, J.E., 2004. A Pattern Search Filter Method

for Nonlinear Programming without Derivatives. SIAM

J. Optim. 14, 980–1010.

Audet, C., Hare, W., 2017. Derivative-Free and Blackbox

Optimization, Springer Series in Operations Research

and Financial Engineering. Springer International

Publishing, Cham.

Audet, C., Dennis, J.E., Le Digabel, S., 2010. Globalization

strategies for Mesh Adaptive Direct Search. Comput.

Optim. Appl. 46, 193–215.

Conn, A.R., Gould, N.I.M., Toint, P.L., 1997. A Globally

Convergent Augmented Lagrangian Barrier Algorithm

for Optimization with General Inequality Constraints

and Simple Bounds. Math. Comput. 66, 261–288.

Conn, A.R., Gould, N.I.M., Toint, P., 1991. A Globally

Convergent Augmented Lagrangian Algorithm for

Optimization with General Constraints and Simple

Bounds. SIAM J. Numer. Anal. 28, 545–572.

Dennis Jr., J.E., Price, C.J., Coope, I.D., 2004. Direct Search

Methods for Nonlinearly Constrained Optimization

Using Filters and Frames. Optim. Eng. 5, 123–144.

Gierusz, W., Tomera, M., 2006. Logic thrust allocation

applied to multivariable control of the training ship.

Control Eng. Pract. 14, 511–524.

Jenssen, N.A., Realfsen, B., 2006. Power Optimal Thruster

Allocation, in: Proceedings of the Dynamic Positioning

Conference, Houston, TX, USA.

Johansen, T.A., Fossen, T.I., Berge, S.P., 2004. Constrained

Nonlinear Control Allocation With Singularity

Avoidance Using Sequential Quadratic Programming.

IEEE Trans. Control Syst. Technol. 12, 211–216.

Kolda, T.G., Lewis, R.M., Torczon, V., 2006. A generating set

direct search augmented Lagrangian algorithm for

optimization with a combination of general and linear

constraints. SANDIA REPORT SAND2006-5315, Sandia

National Laboratories, Albuquerque, New Mexico,

California, USA.

Kolda, T.G., Lewis, R.M., Torczon, V., 2003. Optimization by

Direct Search: New Perspectives on Some Classical and

Modern Methods. SIAM Rev. 45, 385–482.

Le Digabel, S., 2011. Algorithm 909. ACM Trans. Math.

Softw. 37, 1–15.

Leavitt, J.A., 2008. Optimal Thrust Allocation in a Dynamic

Positioning System. SNAME Trans. 116, 153–165.

Liang, C.C., Cheng, W.H., 2004. The optimum control of

thruster system for dynamically positioned vessels.

Ocean Eng. 31, 97–110.

Snijders, J.G., 2005. Wave Filtering and Thruster Allocation

for Dynamic Positioned Ships. M.Sc. Thesis. Delft

University of Technology, Delft, the Netherlands.

Sørdalen, O.J., 1997. Optimal thrust allocation for marine

vessels. Control Eng. Pract. 5, 1223–1231.

Torczon, V., 1997. On the Convergence of Pattern Search

Algorithms. SIAM J. Optim. 7, 1–25.

Tröltzsch, A., 2016. A sequential quadratic programming

algorithm for equality-constrained optimization without

derivatives. Optim. Lett. 10, 383–399.

Wit, C. de, 2009. Optimal Thrust Allocation Methods for

Dynamic Positioning of Ships. M.Sc. Thesis. Delft

University of Technology, Delft, the Netherlands.

Yang, S.-Z., Wang, L., S., Z., 2011a. Optimal Thrust

Allocation Based on Fuel-efficiency for Dynamic

Positioning System. J. Sh. Mech. 15, 217–226.

Yang, S.-Z., Wang, L., Sun, P., 2011b. Optimal thrust

allocation logic design of dynamic positioning with

pseudo-inverse method. J. Shanghai Jiaotong Univ. 16,

118–123.

Zhao, D.-W., Ding, F.-G., Tan, J.-F., Liu, Y.-Q., Bian, X.-Q.,

2010. Optimal Thrust Allocation based GA for Dynamic

Positioning Ship, in: Proceedings of the 2010 IEEE

International Conference on Mechatronics and

Automation, Xi'an, China. pp. 1254–1258.