ARCHEF
MATHEMATISCH INSTITUUT
UNI VERSITEIT
GRONINGEN
REPORT TW-56
The time optimal course changing of a ship
by
G. J. Olsder
Lab. y.
Schpbouw!uinc!e
Technische
HogschooI
Report TW-56
The time optimal course changing of a ship
by
G.J. Olsder.
I would like to express my gratitude to prof.dr. J.A. Sparenberg,under whose
guidance this investigation, which formed part of the "doctoral" examination, took place.
Also, I wish to acknowledge the discussions with mr. G. van Leeuwen and
prof.ir. J. Gerritsma from the Technical University of Delft.
-2--Contents. page
§i. Introduction 15
The equations of motion of a ship
The problem in the phase space 8
The maximum principle 9
Application of the maximum principle 13
Numerical results 22
Comparison with the first order simulation of Nomoto 27
Iterative procedure to determine i 31
The method of Powell
39
Appendix. Algol program for the method of Powell L1
§i. Introduction.
A ship follows a certain rectilinear course and wants to change to another
one. The initial and final course are prescribed. The problem we will discuss is
how to change the rudder angle with time to realize the prescribed change of course
as quickly as possible.
In order to define the mechanical behaviour of the ship we need its equations of motion [i], which contain the rudder angle 6 as a control variable. We will
consider the linearized equations. The coefficients of these have been determined in an experimental way [9,10].
For this steering problem a theory of trajectory optimization is needed. The
classical calculus of variations is not appropriate. Suitable theories are the
"Dynamic Programming" of Bellman and the "Maximum Principle" of Pontryagin and his associates. Because in general dynamic programming requires a large memory size of the computer [6], we use the maximum principle [13].
The problem is solved in two ways.
Firstly, a geometrical insight into possible optimal trajectories of the
ship is obtained. The desired optimal trajectory cannot be found in a straight-forward manner. Several trials must be made and the prescribed change of course has to be found by interpolation. This method is practically limited to linear
time invariant systems of second or third order. In our case we have a third
order problem.
Secondly, the optimal control is expressed explicitly in terms of solutions
of the related adjoint differential equations. The principal difficulty is that the initial values for these equations are not known beforehand. Iterative
procedures have been proposed for computing these values [5,11]. In these procedures the problem is reduced to the determination of the stationary value of a function of many variables. The gradient method and a modification of this method are used to find this value [11]. The second way of solving this problem is also suitable
for higher order systems of equations.
Our results are compared with those which are obtained by employing the
simplified equations of motion, developed by Nomoto [12]. It turns out that the results agree rather well for not too small changes of course.
§ 2. The equations of motion of a ship.
It is assumed that the ship moves through undisturbed water and air. Moreover we consider only its motion in the horizontal plane, hence for instance pitching and rolling are neglected. The number of degrees of freedom possessed by the ship is therefore limited to three: motion of the center of gravity in the two
dimensions of the horizontal plane and an angular orientation with respect to a fixed direction. We use a Cartesian coordinate system fixed to the ship. The origin of the system coincides with the center of gravity, the x-axis coincides with the axis of symmetry of the ship and the y-axis is perpendicular to this axis. This choice has the advantage that the force required to accelerate the
water surrounding the ship can be put into a particularly simple form [1i].
Fig. 2.1. Coordinates and velocities defining the state of the ship.
The following quantities, partly denoted in fig. 2.1, are introduced:
V : velocity of the ship.
components of V in x and y direction. driftangle.
9 : angular position of the x-axis with respect to the reference axis L. rudder angle.
L : length of the ship.
P : specific density of water.
t : time.
r : angular velocity, i.e. r =
-rn mass of the ship.
I : moment of inertia taken around the vertical axis through the center
of gravity.
Arrows indicate the senses in which the quantities are positive.
The equations of' motion are [l,]: de (2.1) dv (2.2) m( - r.q) X
(2.)
m( + r.v) = Y dt dr (2.1k)I--where X,Y are the forces and N the moment about the center of gravity acting on
- dv dr
the ship. In general X,Y and N depend on v,i,r,5, and
The ship is said to be in the equilibrium position if it moves forward with constant velocity and with driftangle zero, i.e. y = y, r = = 0. The disturbance of the forward speed V with respect to y will be denoted by E, , so
(2.5)
= V - V0
Because we consider only small disturbances of the equilibrium position we
dr
express X,Y,N in a Taylor series of
dt - dt and 6. In the following
relations only the linear terms have been written down.
(2.6) X = + X.TÌ + Xr.r + X5.6 + + X.fi + X.i' + ....
(2.7) Y = Y +
+ y.
+ Yr.r + Y6.6 + + Y.. + Y.. +(2.8) N N + + N.r + Nr.r + N6.6 + N.F, + N.. + N..r
A dot above a symbol denotes its derivative with respect to time. By symmetry
considerations it can be shown that X
= X = X = X. = X. = 0.
The constants6 r ii î r
Y and N enter the equations by the a-symmetric rotation of the screw; they
are very small and henceforth they will be neglected.
Now we can
write (2.2)-(2.li-) in the following form:
(m-X).
X+ (m-Y.). - y..j' =
+ Y.
+ Yr.r + Y5.6 - m.r.v ,-N,.E. -
+ (I-N).i =
+ + Nr.r +5
The quantities m,Y Y.,Y ,Y.,Y ,I,N ,N.,N ,N. and N are constants which have
r r 5 'q r r S
been measured for several types of ship models [9,10].
Often -Xe, is called the "added mass along the x-axisTt, -Y. is called the
-6-In this report it is assumed that E 0. Note that by constant power input
the speed of advance will decrease by manoeuvring, so that in case of a constant
forward speed the power input must be adjusted properly. Then equations (2.10) and (2.11) become
(2.12) (m-Y.). - y..' =
y.
+(Y-m).r +
(2.13) -N..j + (I-N.).
= N.
+ N..r + N.b
We now introduce dimensionless quantities. These new quantities, provided
with a bar, are obtained from the quantities with dimensions in the following
manner: (2.1L)
T=--=
1 ;.2
° Y-
rLm
I . r=
2 ,m
= 2 = ,=
2 V -p.L p.L p.v.L -Y. Y Y. Y-
-
r r b = p.L3 =p.v.L3
r = b = ::'Lp.L
= P.V.L'
r = p.v .L3 L oHaving introduced the dimensionless quantities, the bars are omitted and. formulae
(2.1), (2.12) and (2.13) can be written as
de
r
b
dt
3+br+b5b,
dt = 0311 + o r + c5b,4
where b., o., i 3,Li-,5 are constants which have been calculated from the
(dimensionless) coefficients in (2.12) and (2.13). If (2.18)
D=
then m-Y. -Y. 11 r -N. I-N. 11 r(2.19)
,-7
b3 =
i Y-Y.
r
b=
Y-m
r
-Y.
r
i N YI-N.
r
-Y.
r
Nr
m-Y.i
I-N.
r
Y nb5
CL1. = D N m-Y. n-N.
I-N.
r
Y -m
r
Nr
,
C3 ,C5 =
D-N.
n m-Y.i
-N.
n N n Y b Nr
b§3.
The problem in the phase space.We want to change the course of the ship. Therefore we are only interested
in the three coordinates e,1 and r. The initial situation will be (e,1,r) = (o,o,o)
and the prescribed final situation (y,O,O). The three quantities e,1 and r are described by the equations (2.15)-(2.17). We consider the motion of the ship in
a three dimensional phase space, spanned by the e-, and r-axes (fig. 3.1).
Fig. 3.1. The motion of the ship in the (O,1,r) space.
The question is how to reach the final position as quickly as possible from the initial position by steering the ship. So we have to determine the rudder
angle S as a function of time, 5 = 5(t). The rudder angle can be varied free
within a previously given angular region S
< 5. In our linearized theory we
assume that 5max is sufficiently small. From experiments
[9]
it can be concludedthat the linearized theory is rather accurate with respect to S for values of
S smaller than 20 degrees.
max C e y r
-8
§ !. The maximum rincipie.
We assume that the law of motion of the object to be governed - here the ship - is characterized by a system of first order autonomous differential
equations: = fi(xi,...,xn,U), i = or written in vectorform: (1..2) dx = f(x,u)
The quantities x1,. .. ,x are the phase coordinates of the object; they completely characterize the state of the object. For instance xi,.. Xr are the coordinates
und velocities. The vector u, which possesses \ components, is the control and determines the course of the process with time. For example, if \ 1, u may be the rudder angle of the ship. The functions f., i 1,...,n are continous with
respect to x1,...,x, u and are continuously differentiable with respect to
x1,...,xn.
ihe region U, in which u is allowed to vary, is called the control region.
In case of the rudder angle ô we have a closed interval 8(t)t < ô . We shall
= max
amiit piecewise continuous controls. At a point of discontinuity the rudder angle jumos from its old to its new position. It is clear that in practice this can never be realized. But if we restrict ourselves to ship manouvres which vary slowly
in relation to the time needed for turning the rudder angle from -8 to +8
max max
certainly true for a not too small change of course, the piecewise continuity
of 8(t) is a good. approximation.
The solution xt) of (L..2) is considered to be continuous and piecewise continuously differentiable. The points where the derivative of x(t) has a
discontinuity are the discontinuity points of u(t).
Several types of optimality exist. We consider here, as has been said
already, iime optimality. It is reauired. to find a control u(t) which realizes the transfer of the controlled. object from a previously given initial state x
to a previously givn final state
Xf
in a minimum of time, assuming the existence of such a control.In order to solve this problem we need. another system of equations, the adjoint system:
dî. n f (x,u)
= i =
ctl
iThis system is linear and homogeneous. Hence given x x(t) and u = u(t), for any initial conditions a unique solution exists on an interval t < t < tf [j].
As x(t), the vectorÍ'unction v(t) = (r1(t),...,c(t)) is also considered to be
continuous and. piecewise differentiable.
= E a.
x +b.0
in a. i
where A denotes the matrix (a. ) and b is the vector (b.).
in i
i1,...,n,
10
-Now we can construct a function H, called the Hamiltonian;
n
(L.L.) H(V,x,u) = (x,f) = E 4í .f
n c
In terms of H we can rewrite the equations -i.2) and (Li.)) in the form
(Li.5),(Li.6) = ;
-- = -
,i = 1,.,n.
Then, if u(t), to < t < tf and x0 are known, we can find the corresponding
trajectory x(t). After that the right hand side of (ii.3) is known. So if the initial conditions of (Li.3) are given, the solution r(t), to < t < tf is known.
If we consider H('4c,x,u) to be a function of only u (xand tr are fixed), we
can define the function M(r,x), i.e.
(1.7) M(4r,x)
sup H(,x,u).
u U
We now state the following theorem, called the maximum principle.
Let x(t) be the optimal trajectory which transfers the object from x to
Xf and let u(t) be the corresponding optimal control, then a nonzero continuous
vectorfunction i(t), corresponding to u(t) and x(t) by (Li.)), exists, in such a
way that for all t E [t
,tfIl
the function (Li.Li.), considered as a function of uonly, achieves its maximum for u(t), i.e.:
(Li.8) H(í,x,u)
= M(,x)
and at the terminal time the relation M(V(tf), x(tf)) > O holds. Furthermore
if r(t), x(t) and u(t) satisfy (Li.5), (!i.6) and (Li.8) then M(r(t), x(t)) is
independent of t and hence in M(V(tf), x(tf)) > O we may replace tf by any
t[t ,t
of
1.-The proof of the mximum principle can be found in [13]. -The proof is
rather tedious. In
[21
a proof is given based on the same ideas as in [13]; this proof gives more geometrical insight into the problem. A quite different proof isgiven in [7]; this proof is clear and etegant.
Next we restrict ourselves to linear time optimal processes:
dx. n r
(Li.9)
= E a. x + E b. u , i = 1,. ..,n.dt
cii
in n
31
In our problem u is the rudder angle, hence is one-dimensional. Then -i.9) can
be simplified to
dx. n
i
or written in vectorform
In the linear case the adjoint system has the form dí. n
(.12)
= -E a.
, i = 1,...,n, a=l or written in vectorform:(L .13)
.-A
dt=
iV,where AT is the transpose of
A.
Obviously ir(t) is independent of u(t). The function H(ic,x,u) has the form:(..i)4)
HEE
i4c ax +Eit
b u.a a3
a a
It is clear that H, considered as a function of the variable u [-u , +u ]
max max
achieves its maximum as E í b u achieves its maximum.
a a
The control u(t) is uniquely determined at all points of [t,tfl where
E iV b O, viz.:
a a
(L..15) n u(t) = u .sgn( E (t). b ). a a. max a= i nIt is easily proved.
[13]
that E has only isolated zeroes whenb,Ab,.. .,A b are linearly inependent. In that case (L.ii) is said to be normal.
We assume that our systems are normal.
If the derivative of E î .b with respect to t is O at the zeroes of
a a
E iV .b , then these zeroes are discontinuity-points of u(t). These points are
a.
a a
called the switching points of u(t).
So each control function u(t) is a piecewise constant function, whose valis
are only u and -u . At the points of discontinuity u(t) is not defined, but
max max
that has no consequences.
If we choose initial values for r, ì.e. ifr(t)
= '
the solution r(t) of (1...12) is kno. Note that in the
nonlinear
case we don't have this advantage because the right hand side of (Li.3) stillcontains
theunknown
x(t)and
(or)u(t). If we know iV(t) in the linear case, u(t) can be calculated from ().iL-)
and our problem has been solved. However, the trajectory x(t), which starts from
x at the time to will in general not pass through Xf So the difficulty is how
to choose iV so that the trajectory will pass through Xf An iterative procedure
to find has been developed by Neustadt, see
§8.
If we have found such air,
it need not be the optimal one! It can happen that more IV's exist, which transferx to x . All these trajectories are called extremals. The following theorem
o f
points out a case in which the extremal is unique and hence is the optimal
trajectory.
Suppose the origin u O is an interior point of the control region U and
let Xf be the origin in the phase space, then the extrernal is unique. For the
proof of this statement see [13]. In our case the condition that Xf coincides with the origin is superfluous, because a translation of e does not change the
-
12-Until now nothing has been said about the existence of optimal solutions.
Only in the case of the matrix A being stable, i.e. if all of its eigenvalues
have negative real parts, the following can be said about the existence.
Suppose A is stable and the origin u = O is an interior point of U, then for any point x an optimal control exists, which transfers x to the origin
in the phase space. A proof can be found in [131 and. in [8].
A very important theorem, whose proof can again be found in [13] is the
next one.
If, for a linear time optimal process, at least one control exists, which transfers x0 to Xf then there also exists an optimal control which transfers
X to x
(5.6) (5.7) -at
+Be
-13t + dr (5.3)= c3.+c4.r+c5.6
or written in vectorform(o
o
i\f'e\
df
(5.4)
q = O b3 b) TIr)
O e3rJ
C C3. C e4 b3- c3b4 a- b - at -b3 B e b3. C (5.8)-
Ae
+ C3 c4b3-c3b4, B and C are constants of integration and -a, -1 are the characteristic roots
of
/ -b3 C3
-b4 -C4
Te suppose a and are real and distinct. For all knqwn tes ç.
cand
are real. Now we can set up the Hanailtonian:
(5.9)
H = 1r + 2(b3TI + b4r + b58)+ r3(c3TI + c4r + c58)= H + 8(b52 + c53)
where H is that part of H that is independent of 8.
13
-Note that (5.2) and
(5.3)
don't containe;
they can be solved independently andafter this 8 can be calculated from (5.1).
From
(4.12)
we get for the adjoint system of equations:(o
oO\
(5.5) - O b3 C3
3)
i b4 c41The solution of this system is:
5. application of the maximum principle.
The three relevant differential equations are
(5.1) d8
(5.2)
The maximum principle states that if the solution of our problem exists, then
the function 1-I, considered as a function of 5, attains its maximum at the point
5= 3(t), for every t It,tf]. So, if the optimal trajectory exists, then
(5.10) S = 5 .sgn(b5r2 + c5t3).
max
Only when b5r2 + c5i3 = O, S is not defined, however this happens for two values
of t at most, as can easily be shown as follows. We have
-at -3t def.
(5.11) b5í2 + c5Ií3 =
ae
+ be + ewhere a,b and e are constants. The derivative of g(t) has the form
-3t
(5.12)
dt g(t)
-aae
-be
If we multiply (5.12) by
et
it is seen that - g(t) has at most one root. Bythe mean value theorem it is clear that g(t) has at most two roots.
The function g(t), defined in (5.11) contains the three constants of
integration A, B and C, which are introduced in (5.6)-(5.8). Note that sgn(g(t)) is independent of the length of the vector (A,B,C), because (5.5) is homogeneous. Choosing certain values for these constants, the vectorfunction 4r(t) is known.
We calculate 5(t) from (5.10) and then some optimal trajectory is known. However it is rather tedious to determine the values A,B and C which belong to the
trajectory that passes through Xf. This is discussed in section 8.
Now we try to get a geometrical insight into possible optimal trajectories of the ship. The solution of (5.1), (5.2) and (5.3) is
t (5.13) e
J
r dt,bc5-b5c
(5.1k)=
Ae +Be
+b3c- c3b
-b3Bet
+ b5c3- b3c5 (5.15) r=
bAe
+ bb3c- c3b
where A, B are constants of integration. We know that always SI
= 5max' so the solution of our problem is composed of two types of trajectories of the fora
t
e =
f
r dt atet±P
(5.16)Ae
+B
' r=MAetM2Bet±Q,
+ 1where the following abbreviations are used:
,
-bc5- cb5
b5c3- b3c5
(517)
max b3
C) - c3b1 'max b3 c)
c3b)ct-b3
ft-b3
M2=
b
We first discuss the case
the two dimensional q,r plane,
(n(t), r(t)).
See fig. 5.laN . The coordinates of N are
±
of the differential equations respectively.
M
b
in which ct and 13 are both negative. Considering
we can sketch the trajectories
((t), r(t)) and
and b. The trajectories have a node at the point(±
P,
±Q ).
The points N and N_ are singular points(5.2) and (5.3) with b = b and b =
-5
max max
n
(i,r)
plane-
15-(b)
Fig. 5.1. The trajectories (n, r±) in the (n,r) plane, in the case cx<
0, ¡3<0.
In the (n,r) plane the initial and final position have been settled in the
origin, so a closed curve must be constructed, running from the origin to the origin and consisting of parts of (+,r+) and (,r). Provisionally we do not take notice of e. Supposing the ship starts with b
= max at t = 0, the first
part of the trajectory is a part of the solution (1L,r+), in fig. 5.2 called 1,
passing through the origin. Along this solution the ship approaches N for t
and will never return to the origin. Therefore we need at least one switch of the rudder. Suppose this occurs at the time > 0, further arbitrarily chosen.
The second part of the trajectory is a part of the solution (,r), which passes through (n (r1), r('r1)). Along this solution the ship runs to N_ for t
_.
This part does not pass through the origin either, for then with the aid of the theorem on uniqueness of solutions it is seen that this part would coincide with trajectory n, which is symmetric with i, with respect to the origin. So we need yet another switching point to return to the origin. At most two switching points are admissible, so the third part has to pass through the origin and it iseasily seen that this part is the continuation of 1 to the right. At the point of intersection (which surely exists) of the second and third part at time T2 the
r
Fig. 5.2. The three parts of the trajectory.
The conclusion is that when n < O and. f3 < O, in which case the ship is
called stable, we can always go back to the origin in the (,r) plane. Let us now return to the three dimensional space (O,v,r). The final position will be
denoted by (y,O,O). Choosing T1 every-thing is determined and is a unique and
obviously continuous function of T1. If y is given beforehand, the right value
of T1 can be found. by interpolation.
As will be proved later on we have, in the case of an optimal trajectory,
the remarkable formula:
(5.18)
y = +
[T1-
(T2- T1) + (tf- T2)) , Jwhich is valid when the ship starts its change of course with + .
- max
Choosing the first switch of the rudder near N, T1 can be made arbitrarily
large, while (T2-'r1) and (tf--r2), the times necessary for the ship to run along
the second and third part of its trajectory, remain of the same order. So y may
be chosen arbitrarily; each change of course can be realized.
16
The shaded region in fig. 5.3 is the "attainable region" of the ship. It is
impossible for the helmsman to give his ship a transverse velocity and a
rotational velocity r belonging to points outside the shaded area, when he starts
from rj r = O. This can be proved as follows. We use the theorem [13] which
states that every point in the (rj,r) plane, which can be reached from the origin, can also be reached by an optimal control. Hence a point which cannot be reached
by an optimal control cannot be reached at all. However by optimal steering we cannot escape from the open region, which is bounded by the trajectories (n+,r), passing through N and (,r) passing through N.. Hence we remain in the shaded
region by any steering.
Some remarks for the case a > O, < O. In this case the solutions of (5.16) resemble hyperbolas [31. These solutions are drawn in fig. 5J--.
(+,r+) plane
-P
17
-(1,r) plane
Fig. 5.L. The trajectories (,r,) in the (,r) plane, in the case ct > O, < O.
The solutions (+,r) have a saddlepoint at N+, the solutions (,r) have a
saddlepoint at N. A difficulty arises into this case. Let us turn to fig. 5.5. If we choose T1 arbitrarily, we cannot always return to the origin. If we take the first switch e.g. at point B, we cannot find a second switch, in such a way that the third part of the trajectory passes through the origin. Being at point B, it is impossible to return to the origin, by any steering. Suppose the
contrary; but then also an optimal trajectory, which returns to the origin,
exists and that is not possible. We speak of an unstable ship. The linearized
eauations are no longer valid as the ship moves further away from the origin;
that is the case when is chosen too large. The shaded strip between the lines 1 and m is the attainable region.
is shown in fig.
5.6.
N-"s'
Fig.
5.5.
The possible optimal trajectories.Each change of course can again be realized. The time interval (T2T1) can be
made arbitrarily great in comparison with T1 and (tf-T2), choosing the first switch near, but to the right of point C. So by
(5.18)
each y can be realized.The case a < O, > O can be treated similarly. It is the retrograde motion
of fig.
5.5.
Some remarks for the case a > O, 3 > O. This is of pure theoretical interest. The quantities b3 and c24 are always negative for real ships. Then it is possible
to prove that one of the eigenvalues a, is positive, at most. Considering
(5.16)
it is clear that the case a > O, > O is obtained from the case a < O, < O by replacing t by -t. We have a retrograde motion, compared. with fig.5.2.
This': '"
A I t=o
lì t=tf t--T2 II rFig.
5.6.
The possible optimal trajectories of an unstable ship.-A difficulty steals into this case. If one chooses
t1
too large, e.g. at point Bof fig.
5.6,
an intersection with the third part of the trajectory does not exist,so we cannot find a suitable second switch; the ship will n.m away more and more
from the origin, we again speak of an unstable ship. Also by other steer manoeuvres it is impossible to return to the origin, for then it would be
possible by optimal steering,too. The attainable region is the whole plane but our linearized equations are no longer valid when the ship runs too far away
from the origin.
Each change of course can be realized. Choosing the first switching point near the "criticalt' point C, see fig.
5.6,
approaching C from the right side, thetime interval tf
t2
can be made arbitrarily great, while t1 and (t2-t1) remain of the same order and by (5.18) it is seen that every y can be reached.We now will derive some formulae,which give t2,tf and e as a ftnction of
T1.
These derivations are valid for all types of ships treated above.
The equations of the first part of the trajectory are
at 3t
(5.19)
(t)=Ee
+E2e
+p,
+ 1
(5.20) r (t) = M E
et
+ M2E2ePt ++
11
with M1 and M2 given by (5.17) and
f3-b3 b)4 a-b3
-b
E2={Q_
b)4}.
At time t T1 the rudder switches. The values of (t1) and r(T1) are known by (5.19) and (5.20). The equations of the second part are
at 3t (5.22) (t) = F1e + F2e - P (5.23) r_(t) M1 F1
et
± M2F2e - Q with -at (5.2)4) = E1(1-2e 1), 2 = E2(1-2e1)
The equations of the third part are the same as those from the first part within a translation with respect to t. Therefore the third part can be written as
t t
at 3t
(5.25) (t ) = E1e + E2e + P
(5.26)
r(t')
M1E1et
+M2E2et
+ Qwhere tt = ttf where the final time tf is still an unknown constant. The
trajectories (5.22),
(5.23) and
(5.25), (5.26) have two points of intersection: (t = t1, t' = t1) and (t = 2, ' = t14). Note that will be negative. For therelevant point of intersection the following formulae are valid:
19
OETIL 13T ar
(5.27) E1e + E2e - F1e 2 - F2e + 2P
=
0aT24 13T14
aT2
(5.28) M1E1e + M2E2e - M1F1e
-
M2F2e +2Q=
0.ctT2 3T
Solving e and. e from (5.27) and. (5.28) and syrnplifying these formulae as
far as possible with the aid of the definitions of E1
2 and F1 2 we get
, ,
(5.
29)(5.30)
Eliminating
T9,
by taking the logarithm of (5.29) and (5.30), we get(5.31)
From this formula we calculate 'r24 numerically by a Newton-Raphsontechnicjue, then
(5.32)
and tf =T2
-Finally we calculate y:(5.33)
7=
-
ina
aT
e-2
inT1
J 13 1-2e in 2a
rrdt
+J
T24 T=
f
(M1E1ea + M2E2e13t+Q)dt
+f
(M1F1eat+M2F2e13t_Q)dt T1M1E1 OET
T24\
M2E2 ( 13'r 13T24\\=
Q(2T1 -
- 'r24) +
(e 1 - e)+
e - e a M1F1 ¡aT2
T1 M2F2/
'r2 i3T__
-
)+
Çe-e
).
+ e e aT,aT2
e e = y-aT
1-2e 1 J3 e-2
e=
-1-2eI
r dt ± OET e -2-a'r1
J'
1-2e -e-2
1-2e=0.
tirdt
+ I"rlUdt
J
+ 20 -T1 T2where the indices 1,11
and
III denote the parts of the trajectory. Takingtogether the first and the final terms of the
right
hand side of(5.33)
we obtainusing (5.2L) and (5.32) it turns out that the sum of the last four terms in this
expression is zero, so
(5.))
7 =(2- T2T) =
( T1(T2T1)+
(tf-If the ship started its change of course with b -b . the right hand side of
max'
(5.3L1)
would have been provided with a minus sign.For the unstable ships it can easily be found that the upperbound for the
first switching time, in figures
5.5
and5.6
when the ship is at point C andassuminga>O, ct>O, is:
21
-(5.35)
in 2 T<
1 1.4
-c- rt-'-\rr-ç-
Çc
r..) (c. ¶LA1..4_#r
t
r
C
6,-
-
J
The eigenvalues a and. 3 are calculated from these coefficients and are
a
-0,3623
-0,0688
+O,O1)-7-3,2515
-2,7L-Li-2,715
We see that type (a) is stable, (b) is just stable and (e) is not stable. The
results are shown in fig. 6.1. In this figure T2 and tf have been plotted
against
y/a,
as , and therefore y, is linear in by(5.17)
and (5.3Ls-).In the linear theory ô 0,15 radials would be a good choice.
For small y and small time intervals [0,tfI it is not real to admit jumps
of the rudder as already has been remarked in section Li.. In the next section
this will be discussed. in greater detail. Interpreting fig.
6.1,
we willconsider only those y, of which y/6
>
It is difficult to get type (a), the stable one, out of its original course,
SO T1 lasts long in comparison with the "correction" times
(T2 - T1)
and (tf -T2).Type (a) does not reach large values of
Hiand
rl.Type (c), the unstable one, changes its angular position very quickly, but it also reaches large values of r and. )rl. So the "correction" time
(T2 -
T1)is rather long. By
(5.35)
the upperbound for the first switching time is 17(dimensionless) seconds. For O
< y < rt
radials and <0,3
radiais T does not= =
max=
1approach this limit.
In fig.
6.2
and6.3
some optimal trajectories of type (a) and. (b)respectively, in the physical horizontal plane, have been plotted. The coordinates
X and Y of the horizontal plane are fixed in space.
22
-§ 6.
Numerical results.The calculations have been made with regard to three types of ships, (a),
(b) and. (c). The measured quantities of (a) and (b) are mentioned in
[9]
and[io]
respectively. Type (e) has been derived from (b) by changing its coefficientssomewhat. coefficients (a) type of ship (b) (e) b3
-o,89i6
-0,576
-0,5
-0,2856
-0,283
-0,3
b5+O,iL.53
+0,1795
+0,1795
03-1i,392
-3,886
-3,8
c)4-2,719
-2,237
-2,2
C5-1,225
-1,307
-1,307
-
23-From section 2 it follows that
t
*
x(t)
f
sin(e(t) +o
t
y(t)
f
cos(e(t) + arctg(ï(t)4to *
assuming that V = i and that the process started at t = 0. Because x(t) and Y(t)
are not linear in , we have to make a
max choice for this quantity. In fig. 6.2
and 6.13 the maximum rudder deviation was chosen to be
0,15
radiais. The max(dimensionless) time. The radius of numbers along the trajectories denote the
circle, is 1/Q = L1,53 for type (a) and
the smallest turning circle, the limiting
type (a) has not been drawn in the
0,868
for type (b). The limiting circle oft
1'2tf
18 16 14 12 10 I I I I I I I I I/
/
/
-/
.-i-,. '2
,;..,
,
/
-- --.
.-, ._/
-
..
-- 1--
-
.
--.- #
type of sh!p Ca)
.
"
(b)
s, s, sa(c)
tftf
---.- t)
_-_---->
t
.
.
-I E I I I ¡ I i i I i 0 2 6 8 10 12 14 1618
20
ca44ats-
Y/6rrax
2L
-tf
ti
7-Y
n units of the ength
of the ship
-26-
o 3 4 5 6 7 t I 8 9Fig.
The oitirn8
tr:e.tori e
of tyoe () *n the hor5v.o'aL 1'L;.-.
to cì1n.ge; «
cour-e of
-?,
-,
4
lo
X in units of the tength
of the ship
Y in units of tP 'ength
f tt s-ip
6 t2,',tf/
11.timitrg crc(e
t2 tf-
26
-7 X tn units ot the Len,th
of the stup
Fig.
5..L'e optmnii trajectori-os of type (b) in the horzontnl
ir n
p1n-, oelonging to changes of course of -r-
y
nnd§7.
Comparison with the first order simulation of Nomoto.In [12] Nornoto gives a first order simulation of the steering problem. The equations of Nomoto do not contain the lateral velocity :
de
(7.1)
r (7.2) dr T + r = K.5 ,(7.)
T 27-T and K are constants.
Nomoto argues that K and T can be expressed in the constants b3,b)4,b5,c,c2 and c5
in the following way:
b5c3 - b3c5
(7.)
b3c1 - b1c3max-¡
b+c
05b3c
- bc3
b5c3 -b3c5)
Let us now consider again the problem of change of course by applying the maximum principle to system (7.1),(7.2). The same conditions as in the third
order problem of section 5. are assumed: the process starts at t
=
O at the point r=
O,=
O and ends at tf with r(tf)=
O,(t)
O, e(t)
= . Itfollows from the theory of the maximum principle that at most one switch of the
rudder exists. Choosing the switching time T arbitrarily, T > O the time
necessary for the change of course with angle y , tf can be calculated in the
same way as discussed in the foregoing sections:
(7.5)
tf= T
ln(2eT/T
-I
(7.6)
ye(tf)
= K {T- (tf-T))
.When the ship starts the change of course with 5 +5 , then the plus-sign is
max
valid in (7.6), otherwise the minus-sign is.
The results of the calculations are shown in fig. 7.1 . In fig. 7.2 the relative difference between the results of our method and Nomoto's have been plotted
against If the change of course y is small, the relative difference
becomes large. This can possibly be explained by considering the piecewise
continuity of 5(t). For small y and small time intervals [O,tf.] it is not real
28
-On an average the maximum anguLar velocity of the rudder is 0,07 radials/sec.,
expressed in the real time scale. On the dimensionless time scale, everywhere
used in this report, this quantity becomes i,)-i- radials/sec., when the length
of the ship L 161 meters and. the ship runs 15 knots. Then it can be stated
very roughly that the theory with the piecewise constancy of (t) is valid
when t > !0 b . Taking b = 0,15 radiais this coincides with ' >
0,6 radials.
± max max
For y < 0,6 radials,the inertia of the rudder must be considered. This is
ttf
16 14 12 10 I I I I I I I I Ii?i. (.1.
Hl
-29-
I 0 2 4 6 8 10 1214
16
1820
radiak Y/6max
0.3 0.2 0.1
tf
tf - tf
Nomototf
't
ttf
'
\
\
type of ship (a)
30
-I I I I I I
i
I I8 10 12
16
1618
20
rad tais Y/rDax
As has been discussed in §i4, the knowledge of the initial values for the adjoint differential equations determines the optimal trajectory directly. The
difficulty is that these initial conditions cannot be calculated in an easy way.
An iterative procedure to determine these has been developed by Neustadt [liii. This theory has been worked out for practical purposes by Fadden and Gilbert [5]. The theory will be reviewed briefly.
We start with the equation (!.11). The control variable u will be one
dimensional. The vector b is multiplied by u , so the admissible controls now
max
satisfy the condition tu(t) < 1 . As in §L4., the control function is assumed to
be piecewise continuous. The initial and final conditions x and Xf respectively are chosen in a way different from that in which they were chosen in
§5.
Weassume in this section that the point we have to reach, Xf is the origin of the
phase space and the starting point x is a point different from the origin. Because e does not appear in the right hand sides of (2.15), (2.16) and (2.17), the old initial and final conditions can easily be translated into the new ones. Further
we assume that the process starts at the time t = O, i.e. x(0) = x. Consider the system
(8.1) = + bu
A is an n X n matrix, b is a column vector with n elements. The solution of
(8.1) is given by t
(8.2) x(t) = (t). {x(0)
1
1(s)bu ds},O
where (t) is the fundamental matrix [115] satisfying the equations
(t) = A (t), (o) = I
where I denotes the identity-matrix.
Assuming Xf = x(t) o in (8.2) we can state the time optimal control problem
as: find the admissible control which satisfies t
(8.3) x(0) = x -
f
1()
bu(s)ds
O
for the smallest possible time t.
Define the subset c(t) of the Euclidian n-space by
t
(8.)
C(t) ={ f
(s) bu(s)ds, u(s) admissible )O
c(t) is the set of all initial conditions from which the origin can be reached
in t seconds. The following facts are known [8,13],
31
-
32
-0(t) is compact and convex.
0(t) ' 0(tt) whenever t > t'. 0(t) grows continuously with t.
If x is given, the least time in which the origin can be reached. is the
least t > O for which -x0 e 0(t). Denote this t by tf The point -x is a
boundary point of C(tf) and. since C(tf) is convex, at least one vector p exists, in such a way that
(8.5) -p.x > p.w, for all w C(tf).
p is a vector normal to the support hyperplane of 0(tf) at (-x); p is directed
away from C(tf). From (8.2) it follows that
t
(8.6)
p.w
=f
(p, 1(s)bu(s))ds.O
The interpretation of (8.5)
is,
that p.w takes on its maximum in C(tf) whenw = -x. It is clear that (8.6)
is
maximized by maximizing the integrand., so(8.7)
u(s)
= sgn( p,A difficulty is introduced in the case of
(p,1(s).b)
being zero on a subintervalof [O,tf] of positive length. This is the same difficulty as on page 11 . We
restrict ourselves to normal systems, i.e. (p,1 (s) b ) does not vanish on any interval of positive length (assuming p o).
From the theory of the maximum principle it follows that p is identical with the initial value î of the adjoint equations --.12).
Viz.
compare (8.7)with
()--.i5),
which are the same relations. Equation (8.7) can be written asu(s) = s
(b.T
(s),p)where T denotes the transpose. The vector T (s),p is the general solution of
AT
dt
with the initial condition 4r(0) - i4r = p.
Note that for a fixed t, each f = p determines by its optimal strategy
(8.7) and relation (8.L.) a unique situation w of the object on the boundary of 0(t). In general more than one 'V (even assuming = i) may determine
33
-Before we continue developing an iterative process to determine first something more has to be said. about normal systems. The geometrical interpretation of normal systems is that every support hyperplane contacts 0(t) at only onepoint. In that case c(t) is said to be strictly convex. In proving this fact, the contrary is supposed. A hyperplane contacts 0(T) at two points A and B. By
convexity it also contacts c(T) at all points of the interval AB. Let D be such a point. The points A and B can be reached in the optimal time T; in both cases
always u(t) = 1. From the uniqueness of solutions ( §
Li.) it follows that uA(t) uB(t) on an interval of positive measure on [0,T]. By the linearity of the differential equations D can also be reached in time T, but uD(t) = i
is not valid on the whole interval [0,T]. From §Li- and.
[8]
it is known thatif for a normal linear time optimal process at least one control transferring the origin to D, exists, there exists an optimal control which transfers O to D.
This optimal control satisfies u(tfl = 1 for almost every t e [0,T]. So 1 can
be reached in a time T-T with T > 0, because otherwise we are in conflict with the uniqueness theorems of control. Hence D also lies on the boundary of C(T-T)
and that leads to a contradiction with properties 2 and 3
of c(t).
An example of a system of differential equations which is not normal is:
=
u1, --
= u2, with u11 < i and u21 < 1.The attainable region in this case at time t is a square with sides of length 2t. Define the vectorfunction z(t,4r ) by
t
(8.8)
z(t,)
=
f
1b.sgn(,1(s).b)ds.
o
This function is on the boundary of c(t) and. it follows from (8.5) that
(8.9) V0.z(t»V0) > ,
for
all e c(t),The function t
.z(t,)
=0'
l,1(s).bds
is clearly nonnegative
and is strictly increasing with t.From the fact that the zeroes of (,(s).b) are continuous functions of
,
we conclude that z(t,f) is continuous in Consider the function
(8.10)
f(t,)
= (,(z(t,ï )
ox ))
ofor fixed it is
strictly increasing with
t.The point -x
lies on the boundary
of C(tf), hence it follows from (8.9)that
f(t,4r0) > 0. Let us restrict ourselvesto those
for which
f(0,4í)r.x
is negative. So by the strictly monotoniccharacter of f(t,V),
f(t,1110) =O for some t, O
< t < t. Denote this t by T(r),so that f(T(i' ),i' ) =
0. Since
f is continuous in its arguments, T is continuousin ijí . Hence T('V ) takes on its maximum for that r , which takes care of the
optimal trajectory through the point
-x0.
For an iterative procedure to maximize T(v) the gradient method is used:
(8.11)
y.
= y. + l.grad(T(v.)),j+1 3
with i > O. The initial value y must satisfy (v0,x0)
<
O. i may depend on v.grad(T(v)) = - T(v) =
In [11] it is proved that f(t,v) is continuously differentiable with respect to y
arid
-
34
-7'
t= (e.,x +
f
(s)bu ds),i O
where denotes the th component of the vector and e. denotes the vector
\VJ.
i
(0,0,...»,. . .,o) with the number i on the i place.
f
(
-i(t)b.sgn(v,_1.b))=
v,(t).b
Hence equation (8.11) becomes
(8.12)
T(v.) x+ f
(s)bu(s,v )dson
j y.=v.-1
3+1None of the preceding results depends essentially on the magnitude of y. Therefore after each iteration y. wilt be normalized to length 1.
The choice of i presents a problem. If it is conservatively small, convergence will be very slow. If i is large, the iteration process does not converge. In the
"trial and error" method i is chosen as follows:
i. 1. when T(v.) > T(v. ),
j j-1 j =
j1
li
i.i
when T(v.) < T(v.1).
One step of iteration consists of the following operations:
choose y with y -x < O,
o o o
suppose y. is known,
u(t,j) is calculated by (8.7): u(t,j) = sgn(v.,(t).b),
solve T. from
3 T.
(v.,x0 +
1
(s)bu(s,j)ds) o, 3. if T. > T. then
j =
j-1y.
v-1
j-1-1 3 if T. < T. then j j-1V.
j+1=V.
3-17
(z (T ( y.) T.x+ f
(s)bu(s,j)ds o -1 v., (T )bÌ i T.j1
i x+ f
(s)1u(s,j-i)ds li O 2 'O C(T(v.))\CTV.
3+1N
J
y.
,(T.
)bI J jThe geometrical interpretation is given in fig. 8.1.
V.
plane of support z(T(v
),v.1)
j+i
plane of support
Fig.
8.i.
Geometrical interpretation of the iterative procedure. The dotted lines give the case that T. < T.. Note thatj+i j
in the figure v.1 is not normalized.
The results of theUtrial and errortt-method, which will be discussed in greater detail in the next paragraph, are given in table i on page38. It seems that convergence is very slow. The cause of this bad convergence cannot be seen
from fig.
8.i
. In order to understand the bad convergence we consider the surfaceT(v) to be a mountain on the v-space. The mountain has its top in
í,
with heightT(T). In practice these kinds of mountains often seem te have sharp ridges. The path along which we come to the top will be along such a ridge, as can easily be understood from the theory of the gradient method. Many trials are necessary to
stay on the ridge andnotto slide down its sides. The path will have a zigzag
pattern. After each error the step length i is halved, so after a while we hardly
approach the top.
36
-Let us now return to the original problem as discussed in§3.
The point xwill be (0,0,7), with y = 0.15 rad.ials and
Xf
will be the origin. The calculations have been made for type (a). The vectory,
the initial condition of the functionr(t), defined by ()-i-.12) and being the vector p in (8.7), has three components.
These depend essentially on two parameters, because (Li-.12) is homogeneous. So the components of
y
can be written asV
-sinq
2V cos cp. sin y
V3 COS cp.cos y ,
-
< cp < + - , O < y < 2ii.The vector
y
satisfies lv 1. If we give cp and y certain values, T(v) can becalculated. In fig.
8.2
the contour-lines of equal height have been drawn.It must be noted that in the calculations, see the appendix, first the matrix A has been transformed on the mean diagonal, so that the lines in fig. 8.2 do not
represent the T(v) contours, but a transformation of them. However this does not
influence the ridge character of the contours.
We required that (v,x) was negative, therefore the region - - < p < O is excluded in the figure
8.
2. For p = - the vectory
points from -x in thedirection of the origin of the phase space. The iterative process started at
point A of fig.
8.2.
After 50 steps B was reached, after 100 steps C and after200
steps D. The top, with height
3,26,
has been situated at point P.In table 1 we have given the results of the gradient method, described in this section. The iterative procedure was stopped at the moment that
3
E lv - v
<
j+1 3 =
1=
1This occurred for j 14370. The computer time needed to realize these 14370 steps
was 20 minutes. The last row of the table gives the true quantities to which the
columns ought to converge. These have been determined by interpolation of the results of the geometrical method from section
5. In §9
a more useful iterative process is discussed and applied.Vm
rcJ ia Is
4 s s 2TL 3p n radiais
0.2 0.4 0.6 0.8 fl2.
10
L2 1.4-
38
-TABLE 1. first I second. componenti comp. third comp. O -0,914.7 +0,159 +0,278 +0,796 o +1 ,057 1 X 10 50 -0,961 -0,207 +0,101 +0, 908 +1 , 1156,3
x 10 loo -0,691 -0,719 +0,076+1 010
+1,2)-i-26,3
x 10_2 200 -O,14)40 -0,898 +0,0131,289
1,502
X 3oo -0,382 -0,92)-i-0,005
+1 14-57 +1,670 1,2 X 14.00 -0,359 -0,933 +0,003 +1, 571 +178)4
3,1
X 10 500 -0,314-6 -0,9138 +o,00ll4 +1,6711
8714.1,5
X 1O 600 -0,3)40 -0,9)-i-0 +0,0011 +1 71)-i. +1 9277,6
x 800 -0,332 -0,9)-i-3 +o,0006 +1,806 +2,019 1,9 X106
loco -0,328 -0,9)4-5 +0,000)4+1 867
+2,0809,5
X 10 1500 -0,322 -0,914.7 +0,00018 +1973
2,186
X 10 2000 -0,320 0,91.i-7 +0,00009 +2,0)-i-0 +2, 253 1,2 x 10 2500 -0,318 -0,914.8 +0,00005+2, 093
+2, 3076,o
x l0 3000 -0,317 -0,914.8 +0,00002 +2,114-6 +2, 13516,o
x3500
-0,316 -0,9)-i-9 +0,00000 +2, 1722,
385
3,0 X 14.000 -0,316 -0,914.9 -0,00001 2, 203 +3,Li-031,5
X 14.371 -0,31 56 -0,9)489 -0,00002 +2,229+3,
1357,5
x"exactt'
value s
-0,27)4 -0,962 -0,00003 +1,76)4 +3,056+3,26
numb er first second
the achieved
step-Of iterations vector y J switching time switching time
height
T(v.)
length
y- space
Fig. 9.1. The method of Powell.
We have to find the stationary value of T(v) as a function of y. Let the initial
estimate be point A in fig.9.1 . Determine the gradient line l through A and proceed
along 1 to point B, which is the extremal point on this line. Now proceed along the gradient line m in B to the extremal point on this line, point C. Join A and C by a line n and determine along n the stationary value D. From A to D is one cycle of iteration. The method has second order convergence, by which we mean that, if T(v) would be quadratic in its arguments, the stationary value of T(v) is
achieved in one step of iteration.
To find the stationary value on a line, we also use an iterative procedure. Let E = A+k.grad(A), see fig.9.1, k is called the steplength. When the innerproduct of grad(A) and grad(E) is positive, we have found a lowerbound for k. If the
innerproduct would be negative, e.g. at F in fig.9.1, we have an upperbound for k.
The new steplength will be the middle of the lower- and upperbound. The so found k
will be a new upper- or lowerbound. And so on.
We have three separated steplengths; k1 for 1, k2 for m and k3 for n. As initial value of the k. we take the k. of the previous cycle of the Powell
iteration.
One step of Powell lasts longer than one step of the Tttrial and error" method of section
8,
but still Powell seems to be more practicable. This method, applied to our problem, gives the results shown in table 2. The iterative process wasstopped when
39
-§9.
The method of Powell.Because of the bad results with the method of section
8
we were obliged to look for another iterative process. A useful method seems to be another gradient method, developed by Powell [iL]. We only give here the geometrical interpretation of the method of Powell, detailed information can be found in [i!-].The last row of the tables again gives the three quantities to which its columns ought to converge. The computer time to produce table 2 was 4,5 minutes.
TABLE 2.
T(v.) - T(v.1) T (y.)< 10.
-40
-number
ofiterations the three components of y
first switching
time
second switching time achieved heightT(v.)t
number of iterations to find ki,k2,k31o
-0,947 +0,159+0,278
0,796
+1,057
-0,719 -0,692+0,066
+1 070
+1,218
40
2-0,526
-0,850
+0,022
+1 237
+1 333
42
3 -0,407 -0,91)-i.+o,006
+1 453
+1,540
40
4
-0,355
-0,935
+0,0015
+1,679
1 766
46 5 -0,331
-0,944 +0,0004+1 ,903
+1 991
Li-6 6-0,319
-0,9).i.8
+0,00003
+2, 130
+2,222
Li-4 7-0,315
-0,649
-0,00004
+2,392
+2, 604
+2,502
Li-78
-0,276
-0,961
-0,00003
+1 785
+3,065
+3,265
48
9
-0,273
-0,962
-0,00003
+1 761
+3o6o
+3,266
"exact" value s-0,274
-0,962
-0,00003
1,764
3,056
3,26
Appendix. Some remarks about the Algol 60 program.
The program describing the method of Powell is written in Algol 60, with due
observance of the Alcor conventions, for the Telefunken TR Li computer of the ttRekencentrunfî of the Groningen University.
The program serves to print the steps of iteration as in table 2 on page
i-o.
The read instruction at line 7 expects the values of the real variables
b ,b,b ,c
C) ,c ,DM( = ), THETA (= the required change of course) andthe array EEl:)] ( = the initial values v0,v0,v0).
NORFIX is an output procedure: NORFIX ('(' ')') prints the words placed
between the string quotes; NORFIX (FIXPRT, 105,Z) prints the value of Z with one digit before and 5 digits after the decimal point; in general, if a value
should be printed with p digits before and a digits after the decimal point
(0< p < 11, 0< a <
lo),
then the integer following FIXPRT has to be 100 ±cr.The output-procedure NLCR means: new line and carriage return.
Description of the program.
In the first block, lines Lb-22, a transformation = Hx has been realized, so
that the matrix A has been brought on the mean diagonal in the new equation. The eigenvalues of A, different from zero, are A[l] and A[2]. Note that the condition
(v0,x0) < O has not automatically been satisfied after the transformation.
The procedure NOPM(VG), lines 27-15h-, normalizes the vector VG. After having used this procedure the vector VG has length 1
The procedure NULP (OG,BG,US,FPA,EPS,V), lines 35-Li-9, determines the root V
of the function PPA, so FPA(V) = 0, by means of a Newton-Raphson technique. The
root lies between the lowerbound. 0G and the upperbound BG. Only when the
Newton-Raphson technique does not converge, near a point where the second derivative
disappears, the interval between 0G and BG is halved. The new estimate will give us a new 0G or BG. The quantity US is the initial estimation of the root, which
has to be found with an accuracy of EPS.
The procedure HAP (ET,ETNtJ,RIcHT,BOE,I,EP), lines 62-20h-, determines the stationary value ETNLJ of the function T(v) on the line through point ET with
direction RIORT in the case of being true of the boolean BCE and otherwise with the direction of the gradient at point ET. The iterative procedure to determine EThTJ is stopped when the innerproduct of RICHT and the gradient vector at EThTJ is less than EP. The integer I takes on the values 1, 2 or
3,
which denotes thatthe stationary value on line 1(1 1), on line rn(I = 2) or on line n(I
= 3)
issought, see fig.9.1 . At lines 71-101 the switches of the control function (t),
L2
-exists, that vaLue is given to T2 and Ti becomes zero. The quantity T(v), the
height of the mountain, has been represented in the program by the real variable
TNN, which is calculated at lines i1i5-lLi7. The gradient vector at point TNN is
defined by the array
LEi :3],
which is calculated at Lines 171-172. At the Lines171i.-2O2 the steplength k, see
§9 ,
is determined. In the program k is called K...
i
'COMMENT'
G1J
OLSDER,OPTIíiIZATION.,
... 2
3
'ARRAY' B,XB,ET(/1..3/),A(/1..2/).'
4 'BEGIN'. 'EAL.'.
B3,B4,b5,Ç3,C,.C5,M,.g.,D,THTA.'
5
'INTEGER' IJ.,
6 . 'ARRAY' E(j.1...3/ ) . 7READB3,E4,B5,C3,C4,C5,DM,E,TNETA).,
. . a ..R.83*C4B4*C3..,
9KLAD.=53+C4.,
10
.D.=SQRTCKLA.DKLAD4*R1.,
... i i AAl =A C I i I ) . C KLAD+D ) *, 5., 12AA2.=A(/2j1...LKLD_L.*.5.,
. . 13H/1,1/) =R/63.
H(/1,2/).=C3/B3., H(/2,1/).H(/3,1/ ) =C.,
14.15
'FOR' I.1,2 'DO' H(/I+1,3/).(AC/I/)B3)/C3.,
16 .'FQR.. J.=1,2,3 ...!. DO'. .
17 '3EGIN' KLAD.0., t.FORt
I.=1 'STEP'
1 'UNTIL' 3 'DC'18
KLAD.=KLAD+bt(./J,I/)*EiI/).,
19ET(IJI).=KLAD., Xb(/J/).H(/J,1/)*THETA.,
20B(/J/).DM*(H/J,2/)*B5+H(/.J,3./)*C5).,
21 'END'.,
.22 lEND'.,23 'BEGIN'
,VW. , 24'INTEGER' M,D,W,UNUL,YPLON,YMICM,WW.,1M.,.
25'BOOLEAN' BB.,
26,ARRAYL,ETMIN,MIN,LTVO,LVOMIN,LM,ETNU,K5C/1..3/).,
27'PROCEDURE' N0RM(VG).
28'ARRAY' VG.,
29 'BEGIN' 'REAL' NOR,KLAD.,
30
'INTEGER'
.,31
NOR.=O., 'FOR' M.=1'5TEP'1'UNTIL'3.1
'DO'32 'BEGIN'
KLAD.=VG(/M/).,NOR.=NOR+KLAD*KLAD.,
33 'END'., NOR.SGRT(NOR).,
'FOR'M.=1,2,3'DO'VG(/M/) .VG(/M/)/N0R.,
34.'.END'.,
S'PROCEDURE' NULP(OG,BG,US,FPA,
EPS,V).,
36 'VALUE' OG,BG., .. .
37
'REAL' OG,BG,U5,EP.3,
V.38 'PROCEDURE.'....FPA.'.
39 'BEGIN' 'REAL' U,UN,F1,A1.
40 'BOOLEAN.'..BOOL.,
41
U.U5., BOOL.z'TRUE'.,
42. SlAP... FPA.(U,,.Fi,,A1,BOOL).,
43
'IF' F14A1 'LESS' o 'THEN' QG.=U 'ELSE'
BG.U.,
44
UN.=UF1/AI.,
S..
45
'IF' UN 'GREATER' BG 'OR' UN 'LESS' 0G 'THEN'
UN.=.5*(OG+bG).,
46. 'F' Ae5WUN) 'LE5S' EPS 'TN.'
V.UN 'ELSE'
47 'BEGIN'
'GOTO' 5TAP.,
...
49 'END'.
'PROCEDURE' FPA(U,F1,A1,BOOL).,
50...'REAL' .U.,F1..,.Aj.,... 51
'BOOLEAN' BOOL.,
5.2. BEG.1N' .!REAL!KLAD.1,KLAO2..., ...
53
KLAD1.P*EXP(_AA1*U)., KLAD2.=Q*EXP(_AA2*U).,
54 ...F1.=KLAD1+KL.AD2+R..,
...
55
'IF' BOOL 'THEN' A1.=(AA1*KLAD1+AA2*KLAD2).,
...,.
57
TNM.=O.,
58...KS(/ifl .=.KS.4/..2./).=KS(/3Í.).i..,RS..=RT.f4.,
79
NCRFIX( ' ( ' NO
VEKTORI
VEKTOR2
VEKTOR3
SWITCH1
' )
60
.(.1 . SW.ITCH2...TFI..NAL AA.TA.L' ).' ).
61
NORMEf).,
62 tBEGIN' 'PROCEDURE
,.ETNU,.R.J.CHT.,.BOE, L.EP) . ,63
'VALUE' ET.,
'REAL' EP.
.. .. 65'INTEGER' I.,
66
'BOOLEAN'....80E.., 67'ARRAY' ET,ETNU,RICHT.,
6 'BEGIN' 'REAL'INPR,KBG,KQG.,
59
'INTEGER'PP,1M.,
pp.=i.'
71 PAT..NOM.:T; .,
72
R.=ET(/i/)*B(/1/)., P.Eï(/2/)*B(/2/)s, Q.=:T(/3/*B(/3/.,
73
W.=W1.,.Y.M.OM.YPLON.=S1.GNf,P)..,
74
'IF' AAl
LES3' O ''-AND' AA2 'LESS' O 'THEN
YMION.=5IGN(R).,
.
7.5 ..
'IF' AAl 'C.REATER' Q !ANAA2.'GREATELR
O iTHP' YPLON.=Sj(R).,
'IF'AAl 'GREATER' O SAND' AA2 'LESS' O)OR'
.. 7._T (AAl 'GREAT-ER.'
2 '.AND.AA1 'LESS
QJ'TNE.'
78YPLON.=SIGN(Q).,
.79
'IF' AAÎ*P/(AA2*Q)
'NOTLESS.'.. O 'THEN.'60 'BEGIN'
'IF' YPLON*YMION 'GREATER' O 'THEM'
...&].
Ti.T2.0
82 'ELSE' .. 83 'BEGIN' KLAT.iO,,VW.=-1O.,
. . . . 84NULP(VW,KLAT,O,FPA,.î-6,T2} .,
."85 TEMD.,
.. .86 'END' 'EL5E'
.87 EGIN' TP4UL.=LNt._AA1*P/(AA2*Q)/CAA1-AA2.,
88US1.TNUL-2., U52.TNUL2.,
89
FpA(TNUL,.EX,KLAD,.FALSE').,
90
'IF' ABS(EX) 'LESS' f-15'THEN'
91
T1.T2.=TNUL 'ELSE'
92 'BEGIN' 'IF' Y1I0N*EX 'GREATER' O 'THEN' T1,=O 'ELSE'
93 iBEGIN' 'IF' TNUL 'GREAT.R' -10 'THEN' KLAT.=-1O 'ELSE'KLAT.=T,'UJ_-1O.,
94
NULP(KLAT,TNUL,US1,FPA,f-6,Ti).,
95 'END'.,
96
'iF' YPLON*EX 'GREATER' O 'THEN'
97
T2.=T1 'ELSE'
98 'BEGIN' 'IF' TNUL tLESS, 10 'THEN' KLAT.=1O 'ELSE'
LAT.=TNUL+1O.,
99
r4ULP(TNUL,KLAT,UB2,FPA,f-o,T2).,
100 'END'.,
101 'END'.,
102 'END'.,
103
'IF' Ti
'LESS' O 'TriENtT1.O.,.
104
'IF'T2 'LESS' o 'THEN'T2.=o.,
105
D.=1., 'IF' TI 'N0TGREATRz O
'AND'
T2 'GREATER.' Q 'THEN'. D.-1..,.
106 'BEGIN' 'REAL'
Yi1,Y22.,
107
'ARRAY' ETB(ui..3/).,
.108
'PROCEDURE' YIA(U,F1,A1,BCOL).,
109
'REAL' U,F1,A1..,
110
'BOOLEAN' BOOL.,
111 4SEGIN' 'ARRAY' KLAD(/1..2/1.
112
'FOR' M.=1,2 'DO' KLAD(/M/).=EXP(-A(//)*U).,
113
114
'FORv
M.=1,2
'DO'
Fi.FI+ET(/M+1/)*XB(/M+1/)+ETB(/'i+1/)*
115
( 1-KLAD( /4/ Y) /A( /4/ ) .116
'IF' BOOL 'THEN'
fl7....'BEGiN'AI.=ETBíJ/).,'FOR' M,i2 'DO' A1,A1+ETB(/M-i-i/)*KLAD(/M/).,
118 'END'.,
119 'END'., 'PROCEDURE' Y2A(U,F1,A1,BOOL).,....
120
'REAL' U,F1,A1.,
121
'BOOLEAN' BOOL.,
122 'BEGIN' 'ARRAY' KLAD(/1..2/).,
123
.'FOR'
M.1,2'DO' KLAD(/M/).EXP(-A(/M/).*U).,
124
F1.=Y11_D*ETB(/1/)*(U_T1).,
125
. . 'FOR' M.=1,2...'DO'Fi.=F1-D.*ETB.(./M+1/)*(Exp(-A(,M./)*T1)-126
KLAD(/M/))/A(/M/).,
127 ...'IF' BOOL 'THEN'
128 'BEGIN' A1.=_D*ETBU1/),, 'FOR' M.=1,2 'DO'
129
. A1.=A1_D*E.T5(/M+1/)*KLAD(I.M./)...*...
130 'END',,
131 'END.., 'PROCEDURE' y3A(U,Fi,A1,BOOL).,
132
'REAL' U,F1,A1.,
133
'BOOLEAN' BOOL.,
134
'BEGIN'
'ARRAY' KLAD(/1..2/.,
135
'FOR' M.=1,2'DO' KLAD(/M/).=E.XP.4-A4/M,')*tj.)..,
136
F1.=y22+D*ETB(/1/)*(U-12).,
137
'FOR' M.1,2 'DO'
. . ..
139
A1.=D*ETb(/1/).
'FOR' M.=1,2 'DO'
A1.=A1+D*ET(/+1/)*KLA(/)..,
140 t.Et4D.I..., . . . 141UNUL.SIGN(P+Q+R).,
. L5 -'FOR'.=i,2.,.3...
143
Y1A(T1,Y11,KLAD,'FALSE')., Y2A(T2,Y22,KLAD,'FALS').,
144
..VW,=Q.,. 'IF ...T.2 ...L.EGS'..lO...f.THEN.'KLA.T.=1O 'ELSE' KLAT.T2+1Q.,
14
'IF' Yll 'GREATER' O ,THENtNULP(V,T1,.5*T1,YiA,!-6,TN)
'EL.'
146
'IF' Y22 'GREATER'.. O.. 'T.H.EN.NUL.P....T1.,T2,.5*'T1+T2),Y2A,f-6,TNi).'ELSE'147
NULP(T2,KLAT,T2+1,Y3A,-6,TNN).,
148 'ENO'.,
149 'E3EGIN' 'INTEGER' M.,
150
. ,.REAL...,.P.ROCE.DURE' NE.T.1..ÇT.5...1..,... 151'REAL'TS.,
.1.52'BEGIN.!. 'REAL' K,
...
.... 153KL.=UNUL*B(/1/).,
...154 'IF' Ti 'GREATER'. .TNi
'THEN.'.NETl....= ...(/1/)±L*T).ELSEt
155
'IF' T2 'GREATER' INN 'THEN' NET1
(XB(/1/)+KL*(T1_-D*(15_T1)))
156
'ELSE'157
NET1.
(XB(/1/)+KL*(Ti_D*(2*T2_T1_T5H).,
158 'END'., 'REAL'
.!..P.RO.CED.URL' NET(I,T5 i..,159
'REAL' 15.,.
160
'INTEGER'
I...161 'BEGIN' 'IF'
T].'GREATER' INN 'THEN'
162
NET.=
(XBC/IJ)+UNUL*B(/I/)*(1_EXP(_A(/,Ii1)*T5))/A(/I_1/))
163
'ELSE' 'iF' 12 'GREATER' INN 'THEN'
164
NET.=
(XB(/I/)+UNUL*B(/I/)(1_EXP(_A(/i1/)4T1HD*
165
