Issues of numerical accuracy and stability in inverse simulation

Simulation Modelling Practice and Theory 16 (2008) 1350–1364

Contents lists available at ScienceDirect

Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat

Issues of numerical accuracy and stability in inverse simulation Linghai Lu a,*, David J. Murray-Smith b, D.G. Thomson c a b c

Flight Science and Technology Research Group, Department of Engineering, University of Liverpool, C/O Room 107, 19 Abercromby Square, Liverpool L69 7ZG, UK Department of Electronics and Electrical Engineering, University of Glasgow, Glasgow G12 8LT, Scotland, UK Department of Aerospace Engineering, University of Glasgow, Glasgow G12 8QQ, Scotland, UK

a r t i c l e

i n f o

Article history: Received 3 June 2005 Received in revised form 26 May 2008 Accepted 3 July 2008 Available online 18 July 2008

Keywords: Inverse simulation Model inversion Newton–Raphson Nelder–Mead

a b s t r a c t Inverse simulation algorithms based on integration have been widely applied to predict the control input time histories required for aircraft to follow ideally defined manoeuvres. Several different inverse simulation algorithms are available but these different methods are all subject to a number of numerical and stability problems, such as high frequency oscillation effects and also lower frequency oscillatory phenomena termed ‘‘constraint oscillations”. Difficulties can also arise in applications involving discontinuous manoeuvres, discontinuities within the model or input constraints involving actuator saturation. This paper has shown that the dynamic response properties of the internal system are the cause of the so-called ‘‘constraint oscillation” phenomenon. In addition, a new inverse simulation approach based on the constrained derivative-free Nelder–Mead search-based optimisation method has been developed to eliminate problems of discontinuities and saturation. Simulation studies involving nonlinear ship models suggest that this new approach leads to improved properties in terms of convergence and numerical stability. Crown Copyright Ó 2008 Published by Elsevier B.V. All rights reserved.

1. Introduction Inverse simulation aims to determine the system inputs required to produce a given response, defined in terms of the system output variables. It is commonly carried out either by a direct approach based on differentiation or iteratively using integration methods. The first published accounts of the problems of inverse simulation were concerned with aircraft applications and included papers by Kato and Sugiura [1] and Thomson [2]. Their methods involved numerical differentiation of the vehicle state variables, with respect to time. Although this approach gives rapid convergence it is not generic and requires restructuring when the model is modified due to coupling between the mathematical description of the system and the inverse solution algorithm. Sentoh and Bryson [3] defined the inverse process as a Linear Quadratic (LQ) optimal problem that involves minimisation of the integral of a weighted sum of the squares of deviations from a straight flight path and the squared values of control surface deflections. However, this method suffered from significant practical limitations and was relatively cumbersome to apply. Moreover, both these classes of methods are not suitable for handling redundancy problems where the number of control inputs is greater than the number of path constraints [4]. In the early 1990s members of a research group at the University of California, Davis [5,6] proposed what is now the most commonly used approach to inverse simulation. In this methodology the inverse problem is solved by an iterative numerical scheme together with the Newton–Raphson (NR) algorithm. Unlike earlier methods based on the differentiation approach, the structure of this algorithm means that this integration-based method is less model-specific. It can therefore accommodate different models without any restructuring of the algorithm itself. One of the drawbacks of this approach is that it is an

* Corresponding author. Tel.: +44 151 794 8508; fax: +44 151 794 6841. E-mail address: [email protected] (L. Lu). 1569-190X/$ - see front matter Crown Copyright Ó 2008 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.simpat.2008.07.003

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

1351

order of magnitude slower than the method involving differentiation with respect to time. Other approaches to inverse simulation include the Broyden, Fletcher, Goldfarb and Shannon (BFGS) quasi-Newton method [7], the Levenberg–Marquardt (LM) least-squares optimisation approach [8], a global optimisation method proposed by Celi [9], the two timescale simulation method of Avanzini, de Matteis et al. [4,10], and a sensitivity-analysis (SA) based approach introduced recently by Lu et al. [11]. All of these approaches have shown advantages, to a greater or lesser extent, in terms of the solutions of specific inverse simulation problems. It is known that the most widely used integration-based approach to inverse simulation [6] may give rise to high-frequency oscillations superimposed on the desired variables. Lin [12] suggested that these oscillations could excite uncontrolled state variables and showed, through an illustrative example, that this could stimulate further high-frequency oscillatory effects. In addition, the redundancy problem which arises in some cases introduces a further difficulty for this approach due to the non-square nature of the Jacobian matrix. In addition to the high-frequency oscillations and the redundancy problem, a phenomenon termed ‘‘constraint oscillations” has been found which can lead to oscillatory phenomena involving relatively low frequency oscillations in the results. Furthermore, inverse simulation methods can also be used within output-tracking control systems and inversion-based feedforward controllers [3,4,13]. A more accurate model, which includes more physical information about the system such as input saturation, used to design a feedforward controller can lead to better tracking performance [14,15]. However, in comparison with the effort applied in the past to exploring these issues concerned with oscillatory phenomena in inverse simulation solutions, previous investigations have given less consideration to situations involving saturation constraints or to discontinuities within the model or in the manoeuvres [4]. In fact, these effects present a challenge to traditional approaches to the inverse simulation problem, involving not only the integration-based approaches [6] but also the differentiation-based methods [1], and many of the other approaches outlined above [4,7,8,11]. All of these techniques of inverse simulation require partial derivative calculations which depend on properties of continuity and smoothness of the model and of the manoeuvre. Therefore if a saturation situation is reached, the partial derivatives of all the tracked outputs with respect to the saturated variables become zero. The Jacobian matrix consequently will be singular and NR iterations will not converge. The same problem occurs in evaluation of the Hessian matrix used for the two-timescale method [4,7]. To avoid the above problems and to achieve better convergence properties, a new algorithm for inverse simulation based on the constrained Nelder–Mead (NM) method has been developed. It is well known that the NM algorithm, which is a search-based approach, can handle discontinuities satisfactorily, particularly if they do not occur near the optimum solution [16–18]. Furthermore, the derivative-free property can facilitate investigation of numerical issues that exist in the traditional inverse simulation methods. The paper begins, in Section 2, by reviewing the traditional integration method based on the NR algorithm. Section 3 discusses issues of numerical accuracy and stability for this algorithm. Section 4 investigates the constraint oscillation phenomenon on the basis of information gained from the application of analytical methods for the generation of inverse models. Section 5 presents the mathematical development of the new constrained NM method and describes applications involving two ship models. 2. The integration-based approach This section provides a summary of the integration-based method of Hess et al. [6] which forms a useful reference against which other methods of inverse simulation can be compared. A nonlinear system may be described by equations of the form:

x_ ¼ f ðx; uÞ y ¼ hðx; uÞ

ð1Þ ð2Þ

where f 2 Rm is the set of nonlinear ordinary differential equations describing the original system, h 2 Rp is the set of algebraic equations that construct the expected outputs, and u 2 Rq is the input vector. The vector x 2 Rm is the state-variable vector and y 2 Rp is the vector of output variables. Following the discretisation of the continuous model of Eqs. (1) and (2), the input-output relationship of this nonlinear system can be defined as follows:

xðt kþ1 Þ ¼

Z

t kþ1

_ xðtÞdt þ xðtk Þ

ð3Þ

tk

yðt kþ1 Þ ¼ h½xðt kþ1 Þ; uðtk Þ

ð4Þ

If the term on the left-hand side of Eq. (4) is now replaced by the ideal output values yd(tk+1), where the subscript ‘d’ is used to represent the desired value, Eq. (4) can be rewritten as:

h½xðt kþ1 Þ; uðt k Þ  yd ðtkþ1 Þ ¼ 0

ð5Þ

In the traditional algorithm, the NR-based method is used to find u(tk) by the following iterative relationship:

uðnþ1Þ ðtk Þ ¼ uðnÞ ðt k Þ  ðJfh½xðnÞ ðt kþ1 Þ; uðnÞ ðtk ÞgÞ1 f E ½xðnÞ ðtkþ1 Þ; uðnÞ ðtk Þ

ð6Þ

1352

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

where fE[x(n)(tk+1),u(n)(tk)] = h[x(n)(tk+1),u(n)(tk)]  yd(tk+1). In Eq. (6) the term J represents the Jacobian matrix of system outputs at the end of the time interval Dt (from tk to tk+1) with respect to input variables. If, in Eq. (2), there is a direct analytical relationship between input and output the Jacobian matrix may be obtained directly. Otherwise an approximation must be used as follows:

J ij ¼

oyi ðtkþ1 Þ yi ðuj þ Duj Þjtkþ1  yi ðuj Þjtkþ1  ouj ðtk Þ Duj

i ¼ 1; 2; 3; . . . ; p and j ¼ 1; 2; 3; . . . ; q

ð7Þ

where uj and yi are the jth and ith elements of the input and output vectors, respectively and Duj is the perturbation in uj at time tk. In this equation, the superscript n has been omitted for additional clarity. It should be noted that Rutherford and Thomson [19] have presented a slightly modified approach for calculation of the Jacobian matrix which involves perturbing Duj in both the negative and positive directions. When a redundant problem is dealt with, the Jacobian matrix is not square and it is not possible to use standard methods of matrix inversion in the NR iteration scheme, as shown in Eq. (6). Hess et al. [6] proposed the use of the pseudo-inverse matrix as a means of finding the roots of Eq. (5) when J is non-square. 3. Numerical issues and stability of inverse simulation with the integration-based NR iterative scheme It is not a trivial task to obtain the inverse response from the equations of motion of a vehicle or other system and many problems can be encountered. Because of the widespread adoption of the integration-based approach to inverse simulation [5,6] and the history of relatively successful applications with this type of method, this section focuses on numerical issues and the stability of this type of algorithm. 3.1. The problems of high-frequency oscillations and redundancy It has been commonly accepted that the integration-based iterative type of algorithm can provide good results with an appropriate choices of step tolerance limit and sampling rate and with appropriate initial values for the NR algorithm. However, as shown in Eqs. (6) and (7), the calculation of the Jacobian matrix J depends on dynamic properties of the state variables and its accuracy is determined by the discretisation interval. As a consequence, accuracy and stability problems may arise in the evaluation of J and previous investigations have shown that the practical application of this technique may involve failure of the inverse simulation process to converge, poor convergence or convergence with superimposed high-frequency oscillations. Also, redundant situations where there is a non-square J matrix are commonly encountered. The phenomenon of high-frequency oscillations in inverse simulation solutions has been extensively investigated in earlier work [5,6,11,12,19]. Hess et al. [6] described the phenomenon as involving multiple local-minima of the error equation fE in Eq. (6). They suggested that two modifications to the algorithm could alleviate this problem of numerical instability. The first of these involves adjusting the discretisation interval while the second relates to changes in the method for calculating the Jacobian matrix to achieve improved accuracy. In addition, they also demonstrated that the high-frequency oscillations could be smoothed by introducing a digital low-pass filter without degrading the values of the control inputs obtained. In contrast, Lin [12] considered this phenomenon in terms of unobservable motions that are excited during the inverse process. Furthermore, with respect to the first suggestion made by Hess et al. [6], Lin has emphasised that, instead of eliminating this problem, the use of a small time interval could excite uncontrolled state variables and this would further aggravate the high-frequency oscillation effect. Rutherford and Thomson [19] tried to solve this problem in a different way, specifically for helicopter flight mechanics modelling applications. Their approach was aimed at finding a better qualified Jacobian matrix from physical reasoning, by taking account of the difference of perturbations in the controls in terms of the blade aerodynamic force variables and velocity variables. In practice, the perturbations have an instantaneous effect on the blade aerodynamic forces and therefore on accelerations, but not on velocities and displacements. Changes in those latter variables cannot happen instantaneously and can be observed only after a finite period of time following a perturbation. This emphasis on physical reasoning and on the underlying dynamic properties of the model is also of direct importance in the two timescale method of Avanzini et al. which shows advantages in terms of numerical stability in that it eliminates the high-frequency oscillations by reducing the order of the model being investigated [4,10]. The redundancy problem introduces a further difficulty due to the non-square nature of the Jacobian matrix in such cases. Hess et al. [6] suggested that this issue could be avoided by means of the pseudo inverse. However, this could not always guarantee the convergence of the solutions. de Matteis et al. [7] demonstrated that their local optimisation algorithm can handle redundant problems by considering it as a particular case of a more general approach. In contrast to this, the pseudo-inverse solution minimises the norm of the control effort and the authors [7] suggest that other application-related cost functions could be defined within the framework of the local optimisation algorithm in order to properly tailor the resulting trajectory, while satisfying the considered constraints. Moreover, the use of computationally efficient minimisation routines such as the sequential quadratic programming (SQP) avoids the problems related to the singularities. As mentioned above, Avanzini et al. [4,10] overcame this problem in a helicopter application by reducing the model order through their two-timescale method with added constraint conditions. In addition, some other researchers [8,20] have provided feasible solutions

1353

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

for cases involving the redundancy problem but these methods are recognised as having other disadvantages, such as difficulties in determining the necessary performance indexes. Finally, the recently introduced SA method [11] can deal with the redundancy problem as a natural consequence of the structure of this special algorithm. In addition, it also can solve the problem of high-frequency oscillations by increasing the number of integration time steps within the discretisation interval for the inverse simulation process but at the cost of increased computational complexity. 3.2. The constraint-oscillation phenomenon Constraint oscillations involve relatively low-frequency oscillatory behaviour compared with the higher frequency numerical phenomena discussed above. Thomson and Bradley [21] showed that this phenomenon exists not only in data from inverse simulations but also in the control input behaviour of pilots measured during flight. In addition, the occurrence and the properties of this oscillatory mode can be predicted by analysing the system matrix of the corresponding inverse linearised system. Anderson [22] has suggested that the root cause of this phenomenon is associated with the severity of the constraints and the model order and that it may result from the excitation of internal zero dynamics. Anderson [22] and also Thomson and Bradley [21] used dynamic analysis of the inverse linearised model to reach conclusions about whether an oscillation observed in an inverse simulation solution was a constraint oscillation or not. However, the system model linearised around an equilibrium point cannot represent the whole envelope of a complex system such as a manoeuvring rotorcraft, especially in cases involving severe manoeuvres. 3.3. Problems associated with input saturation and discontinuities In the case involving the inclusion of saturation limits, the convergence characteristics of the inverse simulation algorithm based on the NR approach may not be as simple as the situation without saturation limits. The inclusion of saturation limits allows the inverse simulation to be more closely related to the characteristics of the actuators and other related mechanical or electrical subsystems, thus providing more physical insight but introducing additional difficulties. The problem may be explained by considering Fig. 1 which is a block diagram illustrating what happens in the kth discretised interval for the mth iteration. In this diagram the input u0,k represents the initial value of u at the first iteration for the kth discretised interval. When the chosen manoeuvres are demanding, particular subsystems such as actuators and control surfaces may reach their limits. This means that the amplitudes of the inputs required to perform the manoeuvre would be larger than the saturation level umax or juminj. Therefore, if the mathematical model represents the real physical system accurately enough, the amplitudes of the calculated input values would be larger than the saturation levels. However, as shown in Fig. 1, before being fed into the model block, the absolute values of input variables have to be limited to umax or juminj. This will make the elements of the corresponding column of the Jacobian matrix zero, as can be seen from Eq. (7). Hence, in this case, the NR algorithm fails to converge because the Jacobian matrix is singular. This non-convergence, or no-solution problem, is also consistent with failure of the system to perform such a demanding manoeuvre. As a consequence, it is more physically meaningful to consider the effects of inclusion of input saturation. However, even when the saturation limits are not reached, the inverse simulation process may not be as simple as when saturation limits are not accounted for. Satisfactory convergence of the NR method is known to depend on a proper choice of starting condition [23]. From the mathematical viewpoint, if u0,1 > umax or u0,1 < umin, the actual values u are equal to umax or umin. Therefore, initial values, which the designer considers as good values to assure satisfactory convergence of the inverse simulation process, may be inappropriate and undesirable in this case with saturation limits present. After a one-step forward simulation, the saturation effect is likely to cause the NR algorithm to fail to converge. Hence, initial values must be selected more carefully in such situations, compared with the case without saturation. Trial and error methods may be necessary in dealing with practical applications.

One- Step Forward Simulation

umax Initial Values

u 0 ,k

45º

u

Model

umin y u Updated Inputs

NR Algorithm

yd (kth)

m th iteration Fig. 1. The kth discretised interval of inverse simulation with input saturation.

1354

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

Discontinuous point

y y (t k)

y (t k+1)

tk

t k+1

t, s

Fig. 2. Illustration of a discontinuous point.

In addition to failure to converge at the first-step, unexpectedly large inputs occurring during the iterative process may also cause non-convergence. Although the chosen manoeuvre may be assumed to be smooth and not severe, the NR algorithm may give a result u which is out of the interval between umin and umax at some stage during the iterative process. This computed value may not be abnormal but arises simply as a result of the numerical process. Such a value would not affect the iterative consistency if the model being considered did not include input saturation. However, the inverse simulation process may break down when saturation is included since the updated input has to pass through this nonlinearity. As a consequence, the output from the one-step forward simulation will diverge from what the NR algorithm is ’expecting’. Finally, a discontinuous manoeuvre can also lead to non-convergence of the inverse simulation. This is because a discretisation process is involved in the traditional inverse simulation approach, as illustrated in Fig. 2, where a discontinuous point is located within the time interval tk to tk+1. When the inverse simulation process meets this special point, the initial input values u0,k and the calculated Jacobian matrix may not guarantee that a solution can be found for the large transient value Dy = y(tk+1)  y(tk) and, as a result, the inverse simulation process may not converge. 4. The investigation of constraint-oscillation phenomena Lu et al. [13] have shown that inverse simulation and model inversion share the same objective in terms of calculation of the input from a defined output trajectory. This introduces the possibility of analysing inverse simulation techniques using the methods of model inversion and, thus, of explaining the constraint-oscillation phenomenon which has been found to exist in the inverse simulation process. The illustration presented here can be considered as an extension of the work of Thomson and Bradley [21] and Anderson [22], showing that this phenomenon is related to the internal zero system, or in other words, it depends on dynamic properties of the model rather than on numerical issues. In this section a SISO system is presented to illustrate the potential advantage of using the model inversion approach to explain constraint oscillations. For example, consider the following SISO system:

2

3 2 32 3 2 3 0 1 0 1 x1 x_ 1 6_ 7 6 76 7 6 7 0 1 54 x2 5 þ 4 5 5u 4 x2 5 ¼ 4 0 x3 6 11 6 69 x_ 3 y ¼ x1

ð8Þ

The zeros of its internal system are as follows:

z1;2 ¼ 0:5  7:0534i

ð9Þ

After transformation, the following inverse model can be obtained [24,25]:

n_ ¼ y_ d g_ 1 ¼ 5g1 þ g2  5n_

g_ 2 ¼ 80g1  6g2  6n þ 69n_ u ¼ g1 þ n_

ð10Þ

where yd is the desired trajectory, and g1 and g2 are the state variables for the internal system. The simulated results using inverse simulation obtained for the case of a sinusoidal input applied to the given model equation (8) are shown in Fig. 3. Fig. 4 shows results obtained through simulation using the inverse model Eq. (10). Figs. 3 and 4 both show the constraint-oscillation phenomenon within the early part of the records. In the case when Dt = 0.001 s (dashed line), it can be seen that both methods provide essentially the same results. In addition, compared with model inversion, where the variation of the value Dt only slightly affects the accuracy while keeping nearly the same shape, the value of Dt has a bigger influence on the inverse simulation process. With Dt values such as 0.2 s, the inverse simulation will not provide good results due to the fact that a lot of potential information is lost due to the larger sampling interval. This

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

1355

1.5

Δ t=0.2s

1

Δ t=0.05s Δ t=0.001s

u

0.5

D

0

-0.5

-1 0

2

4

6

8

10

12

Time, s Fig. 3. Constraint oscillations from inverse simulation.

1.5

Δ t=0.2s

1

Δ t=0.05s Δ t=0. 001s

u

0.5

0

-0.5

-1 0

2

4

6

8

10

12

Time, s Fig. 4. Constraint oscillations from model inversion.

is illustrated in Fig. 3 (dash dot and solid lines), in which the inverse simulation results for the larger sampling intervals are inconsistent with results from the inverse model. All these results show that, in practice, it is necessary to select relatively small values of Dt in inverse simulation to accurately capture the system dynamic properties and this may conflict with requirements in terms of the choice of Dt to eliminate the high-frequency oscillation [5,6]. Thus, selection of the Dt value for inverse simulation of a specific system should involve consideration of the internal dynamic properties of that system as well as the convergence characteristics of the NR algorithm in the inverse simulation process. Finally, it can be concluded that it is not the inverse simulation algorithm which leads to the oscillation or involvement of the pilot [21], but an inherent property of the dynamic system. This can be further clarified, using the example presented above, again for the case of a sinusoidal trajectory. In this case it may be shown that the response of the inverse model may be approximated by the equation:

uðtÞ  0:2 cosðtÞ þ 1:13 e0:5t sinð7t þ p=4Þ

ð11Þ

The second term forms the constraint oscillation component and this proves that, for the linear case, the constraint oscillation results from the superimposition of sinusoidal functions and functions with exponential coefficients. The frequency of the oscillations is determined by the internal-system properties since it relates closely to the known locations of the pair of zeros of the internal system Eq. (9), being seven times larger than the frequency of the input. Within the first four seconds, the large amplitude and rapidly changing shape of this component of the response significantly influences the inverse simulation process. In addition, if the value of Dt is greater than 0.4 s, for the example system considered, the inverse simulation will not converge, whereas the model inversion technique still can provide good results. The methodologies on which these two methods are based are responsible for this difference. The inverse simulation process involves the NR algorithm and larger interval values will lead to a failure to converge. Conversely, the results from the inverse model are obtained using a conventional forward simulation method and there are no problems of convergence other than those normally associated with integration algorithms and the fact that, in general, less accurate results are obtained when larger intervals are used.

1356

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

A more detailed analysis of the differences between inverse simulation and model inversion have been previously carried out by Lu et al. [13]. In the case of a nonlinear system the situation becomes more complex but the properties of the system in the vicinity of equilibrium points can be analysed as above.

5. Development and applications of the constrained NM method 5.1. Development of the constrained NM method The determination of the elements of the Jacobian matrix or Hessian matrix has already been discussed in Sections 1 and 2. Optimisation methods based on direct search algorithms, being derivative-free, eliminate the need to determine elements of these matrices and thus alleviate the difficulties associated with discontinuities and input saturation. This introduces the possibility of an alternative approach to inverse simulation that does not require gradient information and might therefore have more general applicability. Lewis et al. [26] have reviewed the history and development of direct search methods and point out that they remain popular because of their simplicity, flexibility, and reliability. Among direct search methods, the most widely used is the downhill simplex method of Nelder and Mead [18]. It is a popular method for minimising a scalar-valued nonlinear function of q real variables using only function values, without any derivative information (explicit or implicit). The latest developments of this method [17,27,28] have expanded its functions so that it can be used to tackle multimodal, discontinuous, and constrained optimisation problems but these developments inevitably make the algorithm more complex. As with the NR method, the NM approach is applied over the interval [tk, tk+1]. One of the main differences between the NR and NM methods is that the former updates the input values by means of Eq. (6) while the latter method relies exclusively on values of a cost function to find the optimal solution [26]. Hence, it is important for the NM method to define a good form of cost function. One suitable and widely used form may be described by the following:

8 min L½uðt k Þ; where > < u2Rq p P > : L½uðtk Þ ¼ fhi ½uðtk Þ; xðt kþ1 Þ  ydi ðt kþ1 Þg2

ð12Þ

i¼1

subject to



umin;j 6 uj ðtk Þ 6 umax;j ;

j ¼ 1; 2; . . . ; q

xðtkþ1 Þ ¼ F½xðt k Þ; uðtk Þ

ð13Þ

where the symbol L[  ] is the cost function and F represents the integration process of the vector field f from the initial conditions x(tk+1). If the NM algorithm fails for the quadratic cost-function form of Eq. (12), the following equation based on the absolute value provides a potentially useful alternative:

8 min L½uðt k Þ; where > < u2Rq  p  P   > : L½uðtk Þ ¼ hi ½uðt k Þ; xðtkþ1 Þ  ydi ðtkþ1 Þ

ð14Þ

i¼1

Many other possible penalty functions exist and some of these may be equally appropriate, or perhaps even better, for this application. The reason for selecting the two particular functions in Eqs. (12) and (14), which are two different norms, is the fact that they are relatively simple to apply. The second constrained condition in Eq. (13), which involves the first-derivative term, can be handled by using the structure of the integration-based approach so that the process of finding solutions is divided into two sub-processes. The first of these involves forward simulation to obtain x(tk+1) and the second stage of the process involves calculation of the solution u(tk) from Eqs. (12) or (14) with the available x(tk+1) values. The inequalities in Eq. (13) are solved by means of two transformations before the solution process of Eqs. (12) or (14). These transformations are based on concepts outlined by Errico [29] and Stevens and Lewis [30] as follows: Step 1: Transformation of input constraints The purpose of this step is to transform the original domain of the input variables into a new space before searching for the solution using the NM algorithm. The unconstrained input variables will be left alone. If an input variable is constrained by only a lower or an upper bound, a quadratic transformation will be performed by means of the following equations:

(

If uj ðt k Þ 6 umin;j or uj ðt k Þ P umax;j ; ua;j ¼ 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi otherwise; ua;j ¼ uj ðt k Þ  umin;j or ua;j ¼ umax;j  uj ðt k Þ

ð15Þ

where ua is the transformed input vector. If both the lower and upper bounds are required, a sin transformation can be defined as follows:

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

    uj ðtk Þ  umin;j ua;j ¼ arcsin 1; 1; 2  1 umax;j  umin;j min max

1357

ð16Þ

The details of both transformations have been fully presented by Errico [29]. Step 2: Transformation back to the original domain with the constraints This step is used to transform the new input domain back into the original domain with the constraints for each evaluation of the cost function. Thus, the search domain for input values is based on the values transformed from Eqs. (15) and (16) and the actual function values evaluated by the NM algorithm have to be obtained by application of a second transformation. In the approach adopted the unconstrained input variables remain unchanged. For the input variables that are constrained in terms only of a lower or an upper bound, the transformation is applied as follows [29]:

(

For the lower bound : ub;j ¼ umin;j þ u2a;j ðt k Þ For the upper bound : ub;j ¼ umax;j  u2a;j ðt k Þ

ð17Þ

where ub is the vector of the transformed input values during NM iterations or at convergence. If inputs are constrained in terms of both lower and upper bounds, a sin transformation can be applied. This transformation is defined as follows:

ub;j ¼

1  fsin½ua;j ðt k Þ þ 1g  ðumax;j  umin;j Þ þ umin;j 2

ð18Þ

The constraint conditions in the cost functions of Eqs. (12) or (14) may be handled by the above two steps so that the resulting values of the control variables ub,j transformed back in the original control domain are bounded for the NM algorithm. Step 3: Finding of a solution by means of the NM algorithm A particular form of NM algorithm, which is the modified version described by Lagarias et al. [16], is used to find the solutions in the constrained domain but with the function values evaluated in the original domain. The algorithm first characterises a simplex in q-dimensional space by q + 1 distinct vertices. Then, based on four rules that involve processes of reflection (q), expansion (v), contraction (c) and shrinkage (r), a new point in or near the current simplex is generated at each step of the search. Then a new simplex can be constructed by replacing a vertex in the old simplex, after the function value from Eqs. (12) or (14) at the new point is compared with the function’s values at the vertices of the old simplex. This process is repeated until the diameter of the simplex is less than the specified tolerance. Optimum solutions are thus found for the step under consideration. If each step converges successfully, the complete input time histories u(t) can be formed by combining together the solutions obtained over each interval. The values of the four coefficients q, v, c, and r used are those recommended by Lagarias et al. [16]. These are also almost universal choices for the standard NM algorithm and are q = 1, v = 2, c = 0.5, and r = 0.5. Step 4: Transformation back to the original domain with the constraints The final solutions from the NM algorithm have to be transformed back to the original domain. This process is similar to Step 2.

5.2. Comments on the NM method The complete computational process is illustrated by the flow chart shown in Fig. 5. Luersen et al. [17] describe a more sophisticated NM method that can globalise a local search by probabilistic restarts. The essence of this approach is to summarise the topology of the basins of attraction in which a fixed total cost can be achieved using a Gaussian Parzen-Windows algorithm [31]. For the method proposed in this paper, the initial guess values for uk+1,0 are the calculated values uk from the previous step. Thus, if the manoeuvre is smooth and continuous, this could be a good starting point. Even in presence of discontinuous points, the probability of the NM algorithm finding a global solution is still high because a piece-wise constant input is assumed within the time interval between tk and tk+1, rather than a set of time-varying signals u(t), which would span a higher dimensional space. Hence, problems associated with discontinuous points in the manoeuvre, which may harm convergence of standard NR schemes as discussed in Section 3, are successfully tackled by the NM method. 5.3. Applications of the constrained NM method To compare the new algorithm with the NR method, two case studies have been selected. They relate to two nonlinear mathematical models used in studies of ship steering control – a nonlinear Norrbin model [32] and a nonlinear container ship model [33]. The Norrbin model is a SISO model and the container ship model has two inputs and two outputs. These models are different from those considered in most previous inverse simulation investigations in that they include input saturation effects for each input control channel and also rudder rate constraints, as shown in Fig. 6. These limitations tend, in practice, to degrade the performance of the vessel and result in limited controllability of the system [34]. In addition to the

1358

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

Original Input Domain uj (t k ) Step 1

New Input Domainua, j (t k)

Step 2

Evaluation of Cost Function by the NM Algorithm in the Original Domain

ub,j(t k)

text

Step 3

NM algorithm

Find Solution ?

No

Yes

Step 4

Move to Next Interval

Fig. 5. Flow chart for the kth interval for the constrained NM algorithm.

δc

.

δ max

δ max

1 s

1 τ Rudder Limiter

δ

Rudder Rate Limiter

Fig. 6. Diagram illustrating rudder amplitude and rate limit [36].

obvious performance degradation that occurs whenever an input ceases to have any further influence on the output response, a further difficulty may be encountered due to the fact that the integral term often included in controllers used within ship steering systems may grow unboundedly when the rudder saturates [35]. In addition, some coefficient values of the models being investigated are not fixed and change according to the sea conditions. This may produce a discontinuity that leads to convergence problems for an inverse simulation based on the conventional NR algorithm. In Fig. 6 the value of the time constant s is selected to be 1 s in the applications considered in this paper. Table 1 gives the numerical values of the saturation limits for each input channel as well as the rudder-rate.

Table 1 Input saturation values for the different ship models

dmax (°) d_ max (°/s) nmax (rpm)

Norrbin

Container

35 7 –

10 5 160

1359

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

5.3.1. Application to a nonlinear Norrbin type model of a ship Inverse simulation has been carried out for a Norrbin type model of the ROV Zeefakkel (RZ) ship [32] for a number of different values of forward speed in the range between U = 3 m/s and U = 10 m/s and a demanded heading angle of 20°. A thirdorder reference model has been used to generate the desired heading response:

Wd cm ¼ Wr s3 þ am s2 þ bm s þ cm

ð19Þ

where am, bm, and cm are constants. In the current application, these values are selected as follows [32]: am = 0.9341, bm = 0.2040 and cm = 0.0182. This choice of reference model can guarantee sufficient smoothness of the heading acceleration. The quality of the results from inverse simulation has been validated by carrying out a conventional forward simulation (FFS) using the calculated inputs. As shown in Fig. 7, the results from the inverse simulation and the forward simulations agree well with each other for manoeuvres for three forward speeds (U = 3, 5, 10 m/s). These results demonstrate the successful application of inverse simulation to the nonlinear RZ ship model without saturation limits. From these results in Fig. 7 it may be seen that the rudder angles are actually beyond the limits defined in Table 1 for part of the time history for the case where U = 3 m/s. Moreover, further tests with the NM algorithm have also been successful with the forward speeds varying from 1 m/s to more than 15 m/s. Inverse simulation has also been investigated with the saturation limits included. As shown in Table 1, for the RZ vessel the rudder limit is 35° and the rudder-rate limit is 7°/s. Results for this situation are quite different from those found without the limiters. The inverse simulation based on the NR algorithm fails to converge for the manoeuvre involving a change of heading angle of 20° if the forward speed U is less than 9 m/s. This problem of convergence failure is a feature of the NR algorithm when saturation limits are reached, as discussed in Section 3.3. In this application it arises, to a considerable extent, from the fact that the inputs or the rate of change of inputs required to track the ideal manoeuvres become larger as the forward speed decreases, as is shown very clearly in the results of Fig. 7. However, the inverse simulation methodology based the NM algorithm, which has been developed to overcome the problems of input saturation, has been found to deal very effectively with this situation. Two typical cases (the most serious situation, purposely selected) with forward speeds 3 and 8 m/s have been investigated. As mentioned above, the NR algorithm failed for both cases. The results for inputs determined using the NM algorithm are shown in Fig. 8. Fig. 8 shows that the NM-based approach achieves good convergence. For a heading angle of 20° and the forward speed of 8 m/s, the required input is within the saturation limits, as shown in Fig. 8a and the trajectory from the FFS also complies well with the ideal manoeuvre in Fig. 8b. However, the rate limit value for the rudder on this vessel is 7°/s and the rates encountered without rate limiting reached 34°/s. Therefore, it is not the rudder amplitude limiter but the rudder-rate limiter that leads to the NR-based approach failing to converge. For a heading angle change of 20° and the forward speed of 3 m/s, the amplitude of the required input reaches the saturation level during two time intervals, for 4 < t < 9 s and 18 < t < 22 s. These saturation effects are responsible for the fact, shown in the results of the corresponding forward simulation in Fig. 8b, that the vessel is unable to follow the desired manoeuvre. In addition, the input result of Fig. 8a for a forward speed of 3 m/s has a damped oscillatory shape, which is also shown in the results of the forward simulation in Fig. 8b. As shown in Fig. 8a, this oscillation is damped and disappears after some time. This is a numerical issue that arises from the NM searching process when the saturation situation is triggered and is not related to coupling with the internal dynamics.

a

b

60

25

20 40

IS -U = 3 m/s 15

Ψ ,deg

δ ,deg

IS -U = 5 m/s IS -U = 10m/s 20

FFS -U =3 m/ s

10

FFS -U =5 m/ s FFS -U =10 m/s

0 5

Idea l

-20

0 0

10

20

30

Time,s

40

50

0

10

20

30

40

Ti me,s

Fig. 7. Inverse simulation of the RZ ship without saturation limits for three forward speeds (Dt = 0.2 s, Wd = 20°, NM method).

50

1360

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

a

60

b 25 IS-U= 3 m/s IS-U= 8 m/s

20

20

15

Ψ ,deg

δ ,deg

40

0

10

-20

5

FFS-U= 3 m/s Ideal

-40

FFS-U= 8 m/s

0 0

10

20

30

40

50

0

10

20

Time,s

30

40

50

Time,s

Fig. 8. Inverse simulation of the RZ ship with saturation limits for two forward speeds (Dt = 0.2 s, Wd = 20°, NM method).

5.3.2. Application to a nonlinear container ship model The nonlinear container ship model [33] involves two inputs – the rudder angle and propeller speed. In this case two kinds of manoeuvre are investigated and these are a turning circle and a pullout manoeuvre [37]. The parameters configured to generate the manoeuvres by a forward simulation process are: (a) the time point for rudder execution (which, for the turning circle, is 10 s after the start); and (b) the set values for the rudder angle and propeller speed, which are 35° (dmin = 10°) and 80 rpm (nmax = 160 rpm), respectively. The cost function for the turning circle and pullout manoeuvres are based on Eq. (14), and take the form shown in the following equation:

L2 ½uðtk Þ ¼ fg_ 1 ½uðtk Þ; xðtkþ1 Þ  y_ d1 ðt kþ1 Þg2 þ fg_ 2 ½uðtk Þ; xðt kþ1 Þ  y_ d2 ðt kþ1 Þg2

ð20Þ

The first-derivatives of the desired values in Eq. (20) are obtained through simulation. The two output tracking quantities for the turning circle are the surge velocity (u) and the sway velocity (v) and the two output quantities of interest for the pullout manoeuvre are r (yaw velocity) and W (yaw angle). Therefore, the inverse simulation problems being investigated involve a number of outputs equal to the number of inputs. The reasons for selecting these outputs are based on consideration of the output variables that are most influenced by the chosen manoeuvre. Apart from the three constraint conditions shown in Table 1, the expression for shaft acceleration in the container ship model also involves a discontinuity as follows:

(

na na > 0:3; n_ a ¼ 5:65 ðn  na Þ60 1 ðn  na Þ60 otherwise; n_ a ¼ 18:83

ð21Þ

where na is the actual shaft velocity. Two sets of tests, without and with input saturation and the constrained function, have been carried out. In Tables 2 and 3 p these results are compared with those obtained from application of the NR method. In these tables the symbol stands for Table 2 Convergence of the NM and NR methods without input saturation for the turning circle and pull-out manoeuvres (container ship)

Dt (s) NM

Turning circle Pullout

NR

Turning circle Pullout

10 p

8 p

6 p  p

4 p p

3 p p

2 p p

 p

 p





1 p

0.5 p p

0.2 p

p

p

p

 p

p

 p















Table 3 Convergence of the NM and NR methods with input saturation for the turning circle and pull-out manoeuvres (container ship)

Dt (s) NM

Turning circle Pullout

NR

Turning circle Pullout

20 p p ? p

10 p p ? p

8 p p p







?

6 p p

4 p p

3 p p

2 p p

1 p p

0.5 p p

p

p

p

p

p

p













?

1361

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

p convergence and ? represents convergence but poor consistency, while  is used to indicate no convergence. The term ‘‘consistency”, as used here, relates to the difference between the results from the forward simulation using the calculated inputs and the corresponding values that define the desired manoeuvre. These symbols are also used in the sections that follow. Results for some selected cases shown in Tables 2 and 3 are plotted in Figs. 9–13. Table 2 shows the convergence qualities of the inverse simulation methods without input saturation. The NM and NR methods achieve similar convergence qualities for the turning circle manoeuvre, as shown in Figs. 9a and 10. The calculated rudder input (d) in Fig. 10 for this case is equal to the value of 35° used in the forward simulation that defined the manoeuvre. However, for the pullout manoeuvre, the NR method fails completely to converge. This is due to the fact that the algorithm cannot deal with the effects of the transient around t = 350 s, as shown in Figs. 12 and 13, and as discussed more generally in Section 3.3 in the context of the problems associated with manoeuvre discontinuities. In addition, the different convergence property for the NR and NM methods without input constraints is associated with the existence of a further quantity which has a limit – the rudder-rate. This limit (5°/s) is present in all the cases considered and it has been found that during the inverse simulation process the rudder rate may sometimes be above this limit for these manoeuvres, both for the NR and NM methods. When the model being considered includes input saturation (Table 3), good convergence is still obtained for both the NR and NM methods for the turning-circle manoeuvre shown in Figs. 9b and 11. The calculated input (d) in this case is limited to the saturation level of 10°. The main reason that the NR method achieves good convergence in this case is probably due to the smooth properties of the turning circle compared with the pullout manoeuvre. However, for the pullout manoeuvre, the

a

b

100

100

-300 y, m

y, m

-200

FFS Ideal

- 900

-500

-1500

-800 0

200

400

600

-500

100

700

x, m

1200

x, m

Fig. 9. Results obtained from the FFS of the container ship using inputs found from inverse simulation without (a) and with (b) saturation limits compared with ideal manoeuvre (Dt = 1 s, turning circle, NM method).

a

0

b 85 80

-10

n, rpm

δ,deg

NM -20

75

70

-30 65

-40

60 0

100

200

300

400

Time,s

500

600

700

0

100

200

300

400

500

600

700

Time,s

Fig. 10. Inputs obtained from inverse simulation of the container ship without saturation limits (Dt = 1 s, turning circle, NM method).

1362

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

a

b

20

85

80

δ, deg

n,rpm

0

-20

75

70

NM -40 65

-60 0

100

200

300

400

500

600

60

700

0

100

200

300

Time,s

400

500

600

700

Time,s

Fig. 11. Inputs obtained from inverse simulation of the container ship with saturation limits (Dt = 1 s, turning circle, NM method).

a

b 140

30

120 20 100

n,rpm

δ ,deg

NM

10

80

60 0

40 -10 0

100

200

300

400

500

600

700

0

100

200

300

Time,s

400

500

600

700

Time,s

Fig. 12. Inputs from inverse simulation of the container ship without saturation limits (Dt = 2 s, pullout manoeuvre, NM method).

a

30

b 120 100

20

n,rpm

δ,deg

NM

10

80

60 0 40 -10 0

100

200

300

400

Time,s

500

600

700

0

100

200

300

400

500

600

700

Time,s

Fig. 13. Inputs from inverse simulation of the container ship with saturation limits (Dt = 1 s, manoeuvre, NM method).

NR method fails to converge for all the values of Dt considered. The convergence of the NM method for this manoeuvre also becomes worse and its good convergence can only be achieved for Dt values smaller than 4 s. This again emphasises the

1363

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

a

b

0.7

250

0.6 200

150

0.4

Ψ ,deg

r, deg/s

0.5

0.3

FFS

100

Ideal 0.2 50 0.1

0

0 0

100

200

300

400

Time,s

500

600

700

0

100

200

300

400

500

600

700

Time,s

Fig. 14. Forward simulation using calculated inputs from inverse simulation for the container ship with saturation limits compared with the ideal manoeuvre (Dt = 1 s, pullout manoeuvre, NM method).

difficulties that can arise due to input saturation. Furthermore, the results for the shaft speed channel (n) in Figs. 12 and 13 show oscillations which begin around t = 350 s. This phenomenon is clearly associated with the discontinuity in the manoeuvre, as discussed already in Section 3.3. In addition, the results in Fig. 14, which relate to forward simulation using the inputs calculated from the inverse simulation (and shown in Fig. 13), demonstrate that good consistency with the ideal manoeuvre has been achieved regardless of the oscillations found in the calculated inputs. Although oscillations are shown in Figs. 12 and 13, the NM method still shows significant benefits taking into account the total failure of the NR method for these cases. The results show that the rudder input channel, which is usually the feedforward control input generally achieves good results and, in practice, this is the important channel in terms of ship steering control. Thus, it can be concluded that the NM method provides a useful alternative to the NR type of approach to investigate sources of difficulty in inverse simulation by excluding some factors such as derivative information. This can help investigators to focus their investigations on issues associated with physical properties of the underlying model and the required manoeuvre and helps to separate these problems from issues of numerical accuracy and stability. 6. Conclusions The paper has described systematically the issues of numerical accuracy and stability in inverse simulation. Firstly, the constraint oscillation phenomenon has been explained with the aid of model inversion methods. Analysis has shown that constraint oscillations relate to internal system properties. Secondly, instabilities and failure to converge for methods of inverse simulation based on the NR algorithm has been discussed in terms of discontinuous manoeuvres, discontinuities within the model and input constraints. Methods for avoiding these issues have been developed using a new derivative-free procedure, based on use of the constrained NM algorithm. This approach avoids use of the augmented Lagrangian method to deal with the constrained conditions and instead applies one-step forward simulation input transformations before applying a pattern-search form of optimisation method. Simulation studies with two nonlinear ship models show that the new approach gives improved convergence and numerical stability properties compared with the NR algorithm for cases that include input saturation in the model or involve a discontinuous manoeuvre. Acknowledgements The authors wish to acknowledge helpful communications from Professor Lin Kuo-Chi of Department of Mechanical Materials and Aerospace Engineering (MMAE) at the University of Central Florida (UCF) and to the anonymous reviewers for all their helpful suggestions. Linghai Lu gratefully acknowledges the award of a University of Glasgow Postgraduate Scholarship and an Overseas Research Studentship from the British Government. D.G. Thomson and D.J. Murray-Smith acknowledge support from the UK Engineering & Physical Sciences Research Council through Grant GR/S91024/01. References [1] O. Kato, I. Sugiura, An interpretation of airplane general motion and control as inverse problem, J. Guid. Contr. Dynam. 9 (1986) 198–204. [2] D.G. Thomson, Evaluation of helicopter agility through inverse solution of the equation of motion, Ph.D. thesis, Faculty of Engineering, University of Glasgow, Scotland, UK, 1987.

1364 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

L. Lu et al. / Simulation Modelling Practice and Theory 16 (2008) 1350–1364

E. Sentoh, A. Bryson, Inverse and optimal control for desired outputs, J. Guid. Contr. Dynam. 15 (1992) 687–691. G. Avanzini, G. de Matteis, L.M. de Socio, Two-timescale-integration method for inverse simulation, J. Guid. Contr. Dynam. 22 (1999) 395–401. C. Gao, R.A. Hess, Inverse simulation of large-amplitude aircraft manoeuvres, J. Guid. Contr. Dynam. 16 (1993) 733–737. R.A. Hess, C. Gao, S.H. Wang, Generalized technique for inverse simulation applied to aircraft maneuvers, J. Guid. Contr. Dynam. 14 (1991) 920–926. G. de Matteis, L.M. de Socio, A. Leonessa, Solution of aircraft inverse problems by local optimization, J. Guid. Contr. Dynam. 18 (1995) 567–571. S. Lee, Y. Kim, Time-domain finite element method for inverse problem of aircraft maneuvers, J. Guid. Contr. Dynam. 20 (1997) 97–103. R. Celi, Optimization-based inverse simulation of a helicopter Slalom manoeuvre, J. Guid. Contr. Dynam. 23 (2000) 289–297. G. Avanzini, G. de Matteis, Two-timescale inverse simulation of a helicopter model, J. Guid. Contr. Dynam. 24 (2001) 330–339. L. Lu, D.J. Murray-Smith, D.G. Thomson, Sensitivity-analysis method for inverse simulation, J. Guid. Contr. Dynam. 30 (2007) 114–121. K. Lin, Comment on generalized technique for inverse simulation applied to aircraft maneuvers, J. Guid. Contr. Dynam. 16 (1993) 1196–1197. L. Lu, D.J. Murray-Smith, E.W. McGookin, Investigation of inverse simulation for design of feedforward controllers, Math. Comput. Modell. Dyn. Syst. 13 (2007) 437–454. M.B. Tischler, Advances in Aircraft Flight Control, Taylor & Francis Ltd., London, 1996. L. Lu, Advances inverse modelling and inverse simulation for system engineering and control applications, Ph.D. thesis, Department of Electronics and Electrical Engineering, University of Glasgow, September 2007. J.C. Lagarias, J.A. Reeds, M.H. Wright, P.E. Wright, Convergence properties of the Nelder–Mead simplex method in low dimensions, Siam. J. Optim. 9 (1998) 112–147. M.A. Luersen, R. Le Richem, F. Guyon, A constrained, globalised, and bounded Nelder–Mead method for engineering optimization, Struct. Multidiscip. Optim. 27 (2004) 43–54. J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313. S. Rutherford, D.G. Thomson, Improved methodology for inverse simulation, Aeronaut. J. 100 (1996) 79–86. K.M. Yip, G. Leng, Stability analysis for inverse simulation of aircraft, Aeronaut. J. 102 (1998) 345–351. D.G. Thomson, R. Bradley, Prediction of the dynamic characteristics of helicopters in constrained flight, Aeronaut. J. 94 (1990) 344–354. D. Anderson, Modification of a generalised inverse simulation technique for rotorcraft flight, Proc. Inst. Mech. Eng., (G) J. Aero. Eng. 217 (2003) 61–73. W. Cheney, D. Kincaid, Numerical Mathematics and Computing, fifth ed., Brooks/Cole Publishing Company, America, 2004. A. Isidori, Nonlinear Control Systems: An Introduction, second ed., Springer, Berlin, London, 1989. S. Sastry, Nonlinear Systems: Analysis Stability and Control, Springer-Verlag Inc., New York, 1999. R.M. Lewis, V. Torczon, M.W. Trosset, Direct search methods: then and now, J. Comput. Appl. Math. 124 (2000) 191–207. R. Chelouah, P. Siarry, Genetic and Nelder–Mead algorithms hybridized for a more accurate global optimization of continuous multiminima functions, Eur. J. Oper. Res. 148 (2003) 335–348. S. Wolff, A Local and Globalized, Constrained and Simple Bounded Nelder–Mead Method [Ver. 2.0], Computer Program, Bauhaus University Weimar, Germany, 2004. J.D. Errico, Bound Constrained Optimization, Computer Program, MATLAB and Simulink Centre, the MathWorks, Inc., 2005. B.L. Stevens, F.L. Lewis, Aircraft Control and Simulation, second ed., John Wiley & Sons, Inc., New York, 2003. R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classification, second ed., John Wiley & Sons, Inc., New York, 2000. M.A. Unar, Ship steering control using feedforward neural networks, Ph.D. thesis, Department of Electronics and Electrical Engineering, University of Glasgow, UK, 1999. T.I. Fossen, T. Perez, Q.N. Smogeli, A.J. Sorensen, Guidance, Navigation and Control Toolbox, Norwegian University of Science and Technology, Trondheim, 2005. E.W. McGookin, D.J. Murray-Smith, Y. Li, T.I. Fossen, Ship steering control system optimisation using genetic algorithms, Contr. Eng. Pract. 8 (2000) 429–443. D.C. Donha, D.S. Desanj, M.R. Katebi, M.J. Grimble, H infinity adaptive controllers for auto-pilot applications, Int. J. Adapt Contr. Signal Process 12 (1998) 623–648. T.I. Fossen, Guidance and Control of Ocean Vehicles, John Wiley & Sons Ltd., UK, 1994. E. Lóez, F.J. Velasco, E. Moyano, T.M. Rueda y, Full-scale manoeuvering trials simulation, J. Maritime Res. 1 (2004) 37–50.