The motion of a disk rolling on a vibrating horizontal plane: Feasible control and path controllability

Mathl. Corn@.

Pergamon

Modelling Vol. 23, No. 4, pp. l-15, 1996

Copyright@1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/96 115.00 + 0.00

SO8957177(96)00001-S

The Motion of a Disk Rolling on a Vibrating Horizontal Plane: Feasible Control and Path Controllability Y. YAVIN AND C. FRANGOS Laboratory for Decision and Control, Department of Electrical and Electronic Engineering University of Pretoria, Pretoria 0002, South Africa (Received

June 1995)

and accepted

Abstract-Let

(X, Y, 2) be an inertial coordinate system and suppose that a horizontal plane is vibrating with an acceleration that is parallel to the (X, Y)-plane. A disk is rolling on the vibrating plane. Given two points A and B both fixed in the (X, Y)-plane. Open-loop strategies are computed for rolling the disk, on the vibrating plane, from A to a small neighbourhood of the point B, during a given time interval [0, tf] and subject to state and control constraints. In addition, the concept of path controllability is introduced and a condition is derived for the disk’s motion path controllability. The derivation of this condition enables one to design closed-loop control laws for the disk’s motion. Keywords-Feasible lability.

command strategies, Rolling disk, Nonholonomic constraints, Path control-

1. INTRODUCTION the control and guidance of the motion of a disk rolling on a vibrating horizontal plane. The dynamical model used here for the disk’s motion constitutes a simplified

This work deals with

model for the motion of a riderless moments, that is, “side inclination” of these moments moments and the the (X,Y)-plane (X, Y)-plane and considered “pedalling”

unicycle. The disk’s motion is controlled by two applied or “leaning” moment and a “pedalling” moment. For each

a servo dynamic model is assumed to interface between the commands to the moments themselves. Given an inertial coordinate system (X, Y, 2) in which is parallel to the vibrating plane. Let A and B be two points fixed in the let [0, tr] be a time interval where tf > 0 is a given number. The problem

in this work is as follows. Find command servo, respectively, such that:

(i) the disk will roll, during (ii) the motion

strategies

to the “inclination”

[0, tf], from A to a small neighbourhood

of the disk will be subjected

to state

and control

servo and the

of B;

constraints

during

[0, tf].

For more details on the constraints, see Section 3. Note that the motion of the disk is subjected to nonholonomic constraints [l]. Such a problem as posed above might arise, for example, in vehicle motion planning (see, for example, [2] for the vehicle motion planning problem) where the vehicle does not possess any sensors and the location of the obstacles on the ground is known and fixed. In general, such a dynamic nonlinear control problem would be treated by the maximum principle [3] or by some other optimal control technique [4,5]. The main feature of this work is that, instead of solving the above-mentioned problem by some optimal control technique, the problem is solved by computing feasible command strategies. The definition of feasible command Typeset by .A&?-?$$

Y. YAWN AND C. FRANGOS

2

strategies

and their computation,

(i) First

the command

strategy

parameters. (ii) A penalty function and the required

that

incorporates

c, a functional

is constructed

all the state

functions

all the state

constraints

function

of the system’s that

speaking,

by introducing

and control

The penalty

in such a manner,

and control

posed above, is, roughly

are parameterized

goals, is constructed.

but it is also, through function

for the problem

a vector

constraints

is a function

and all the required

c of

of the system, of the vector c,

state space trajectories.

it reaches

as follows:

This penalty

the value of zero if and only if goals are satisfied.

(iii) A gradient method is applied on the vector space in which c resides, to bring the value of the penalty function to zero. This yields a solution vector co which induces a feasible command The precise scribed

strategy. definition

in Sections

of a feasible

strategy

and its computation

are given and de-

3 and 4, respectively.

Let (20, YD, Zg) denote coordinate

command

system.

Denote

the coordinates y(t) = (zD(t),

of the center y,

of the disk in the inertial

(X, Y, Z)-

If for any two points

A’ and B’

yD(t), v).

in R4, and for any finite time interval [0, tr], an inclination moment and a pedalling moment can be found such that y will move from A’ to B’ during [O,tf], then we say that the disk’s motion is path controllable. The main goals of this work are (a) to introduce the concept of feasible command strategies and to demonstrate bility through the solution of the above-mentioned problem, and (b) to establish

a condition

The establishment cedure

under which the disk’s motion

of the condition

for the design of closed-loop

for path controllability control

their applica-

is path controllable. is constructive

and it offers a pro-

laws for the disk’s motion.

The problem considered here is to a large extent an extension of the problems treated in [6,7]. In the problem dealt with here, we study of the influence of the vibrating plane on the constraints forces and (through this) on the disk’s maneuverability. This is done in Section 6 by solving numerically

a variety

of problems.

2. THE DYNAMICAL In this work, we consider and K be unit vectors

the control

along an inertial

MODEL

of a disk rolling on a vibrating (X, Y, Z)-coordinate

system.

horizontal Denote

plane.

Let I, J

by k

k = sin 8 cos ~$1+ sin 0 sin c,hJf cos OK,

(1)

a unit vector along the axis of the disk. Then, the vectors ie and i4, given by ie = $& and ig = (l/sin@@ are at all times in the plane of the disk. Define the following vectors i and j, i = cos $,ie + sin $$, j = -sin$&

(2)

+ cos+!$,,

which are always in the plane of the disk. Hence, we have now the following

(3) relation:

(4) where,

using the notation

E = E(8, 4, $),

The Motion of a Disk

cost9cos~cos~ E=

-cosf9cos+sin$

- sin$sin$ - sin4cos+

cosOsinf$cosII, + cosr#sin+ -cosOsin@sin$+cos#cos+

sin 8 cos 4 Note that are playing

sin 19sin #J

- sinBcos$ sin 8 sin T/J cos 8

*

(5)

E = E(0, @,$) is th e matrix of the Euler angles transformation [l]. Here, 8,$ and $J the following roles: 13is the leaning angle of the disk, that is, the angle between the

disk’s axis and the Z-axis, i4 represents

the direction

its axis. Also, the angular WD =

where for 0 = 4 the plane of the disk is vertical to the (X,Y)-plane; of the disk; and Q represents the angle of rotation of the disk around velocity

vector of the disk

wDli i- WDZ~-I- WDsk = W&

+ WD,J + tiozK,

(6)

is given by [l]

wD1

=

wD2 =

(s) sin+- ($) SineCOS$J,

(7)

($) COS+f ($) SineSin@,

(8)

and

In the sequel, the Lagrangian method [l] is applied to obtain the equations of motion of the disk. Using the notation given above, the Lagrangian function for the disk’s motion is given by

where mD is the disk’s mass, a is the disk’s radius, IDj, j = 1,2,3, denote the moments of inertia of the disk about the i, j and k axes, respectively; and (zD, MD, zg) are the coordinates of the center of the disk, relative to the (X, Y, Z)-coordinate system. Using the “thin wheel” apprOXiIIM+tiOn we have 101 = 102 = 0.25mDa2 Consequently, the discussion above and equations (7)-( 9) and (13) yield

where (x,3, r) =

(ZD,

YD, zD>-

and

103

=

0.5mDa2.

Y. YAVIN AND C. FRANGOS

4

In this work, it is assumed that the motion of the disk on the plane involves rolling without slipping. This leads to the condition vg+wDxrD=u,

(15)

at the point of contact between the disk and the plane, where vg = id+ &J + fK, WD = W& + WQ,J + LJD=K, and u = ~11 + uzJ denotes the velocity of the horizontal plane. From equations (l)-(3), it follows that at II, = 0, i = ie and j = i@,that is, at the point of contact, rD = aie. Hence, equation (15) yields a

(z)sinm]

fui-s=O,

(16)

a

($os4]

+,+o,

(17)

case-

2 =0,

(18)

Equations (16)-( 18) (or equation (15)) constitute nonholonomic constraints. Hence, the Lagrange equations, for the case of nonholonomic constraints [l], turn out here to be ---

dt

dt

de

=Ifj+XiasinBcosQ,+Xzasin8sinf$+XzacosB,

- = r@ + Xiacosesin4 lk$

7-w

dL = I? + Xiasin+

*

-

(19)

(20)

X2ac0st9c0sfp,

- Xzacos$

(21)

And, in a similar manner, the Lagrange equations for z, y and z are given by

(22) where Xi, i = 1,2,3 are the Lagrange multipliers, and I@, I’# and l?$ are the respective applied moments. That is, I’0 is the disk’s “inclination” applied moment, l?s is the disk’s direction applied moment, and I’tc,is the disk’s “pedalling” applied moment. Thus, by differentiating both sides of equations (16)-(M) with respect to t and using equations (22), we obtain r; =

Aiasinf?cos$

+ Xzasinesin$

+ X3acosB

=--mDa2[$+($)sin@($+($)cos@)] -mDasine r; = XiacosBsin$

[(%)

cos4+

(2)

sin’$] j

(23)

- X2ac06ec0S$

- mDa,Se[(~)Sin~-(~),OS,], r; = Ala

sin 4 - &?a cos C$= -mDa”

(24) coSe+$&-

-mDa[(%> where l?$, l?s and l?$ are the (generalized) constraint forces.

2 (z)

(2)

sine] (25)

The Motion of a Disk

Define (19)-(21),

q = (8,4,$jT the following

and

r

=

equation

(re,rd,r+)

T.

Hence,

5

by inserting

equations

(23)-(25)

into

is obtained:

M(q)g+h

(

q,$f

>

=r,

(26)

where

M(q) =

I1

0

0

101 sin2 B + 13~08~0

0

0

13 cos e

1scos0

,

(27)

I3 is given in the Appendix.

det M(q) equation

Since

= 111~113 sin2 0,

(28)

(26) yields

[-h(q.$)+r] =@(q,$$r),

d2s

-=M%) dt2 in any region that

excludes

0 = ix, i = 0, fl,

f2,.

(29)

...

Henceforward, it is assumed here that l?4 = 0. Define the following state variables: zi = 8, tidt,2s=$,2s= s, x7 = x, xs = y, 59 = z, 210 = ‘111and zll = ~2. Also, Zs=$.,Lrs=+,x4= that is, xl2 = l?e/Il and xl3 = I’+/13. we introduce here the following “normalized moments,” In this work, the following servo dynamic models, to interface moments and the moments themselves, are assumed: dx12 -=

(b(t)

dt

dx13 -=

(W)

dt where al and us are given positive

-

512)

a1

numbers

between

the commands

and '

to the

(30)

--x13)

a2

(31)

’

and C&= {&(t),

t > 0}, i = 1,2 are the command

signals.

(32)

Then, by using equations (15), (29) and (30)-(31), (assuming everywhere t > 0):

dxl

x

= 52,

the following set of state equations

is obtained

(33)

dxz = Fi(q),

(34)

dxs dt = x4,

(35)

dx4 = Fz(v),

(36)

dxs z =x6,

(37)

dt

dt

Y. YAVIN AND C. FRANGOS

6 dX6 -

dt dn dt= dxa -

dt da dt

=

(38)

F3(7I7),

a[x2 sin x1 cos x3 =

+ x4 co6 z1

sin x3

+ x6

sin 231

+ 210,

a[x.2sinx1 sin23 - x4cosx1 cosx3 - xfjcosxs] +x11,

= ax2cosx1,

(b(t)

-

(42) (43)

212)

dt

a1

dx13 -=

(62(t) - x13)

dt

(40) (41)

dxlo = u, sin(wot), dt dxll = 21acos(wot), dt dzl2 -=

(39)

’

a2

’

where ua and wo are given positive numbers. motion for the problem dealt with here.

Equations (33)-(45) constitute

3. FORMULATION

the equations of

OF THE PROBLEM

Henceforward, we will be interested in the motion of the disk only during a time interval [0, tr], where tf > 0 is a given number. Let 70 = 0 < 7-1 < 72 < . . . < TN-~ = tf be a partition of [0, tf] such that T~+~- TV= A,, i = 0,. . . , N - 2. In this work, the following class of command functions (61,s~) is dealt with. Consider the class of all functions 6, = (61,62) : [0, tr] + R2 such that

[q,~i+l],

b(t) = Ai(t)ci + &(+i+l,

tc

62(t)= Ai(+N+i + &(t)cN+i+l, ISI I Qmax,

t+,q+l],

=

cTi+l- t, A,

B+(t)

’

,..., N-2,

i=O

,...,

N-2,

i=0,...,2N-1,

where A&)

i=O

=

ct- ‘d A,

’

i=O,...,N-2,

(46) (47) (48)

(4%

and amax > 0 is a given number. Henceforward, we denote this class by A. Note that equation (48) implies I&(t)1 I am=, i = 1,2 for all t E [0, tr]. Hence, equation (48) constitutes the constraints on the command functions. Denote 11= (x1,22,... , x13) and let 6, be an element of A. Also, denote by {(t, 77;6,)) t 2 0, the solution to equations (33)-(45), whenever it exists, such that C(O,7; &) = 7, 7) E R13. Define the following sets:

where emin, emax, &, numbers.

, ~7 and

&s are given positive numbers, and XD and 90 are given real

The Motion of a Disk

7

Given a hxed time tf > 0, a partition {ri}Ei’, and a state ~0 E D,. The problem dealt with in this work is to find a command strategy SCE A such that

Equation (52) represents the state constraints imposed on the disk’s motion during [0, tf], whereas (53) represents the constraints imposed on the disk at the final time t = tf. A command strategy 6, E A for which equations (52), (53) are satisfied will be called here a

equation

feasible command strcrtegy.

4. SOLUTION

OF THE PROBLEM

Define the following functions: G(z, A) = [max(a - X,0) +

Gr2(z, &, Xl) =

[max(a

-

X2,0>

x > 0,

min(z + A,0)12,

+ min(.z

-

Xl, 0)12,

x2 >

Xl

(54 > 0,

(55)

for any z E R. Let c = (cc,cr ,..., cz~_i) E R2j”, and let SC be defined via equations (46), (47) and (49). Define the following function F(c) by 2N-1 F(c)

=

c

PoG(ci,

~max)

+ qG((i(tf,

rlo; 6,)

-

2D,E7)

+~8G(&dtf,rlo;~c)

-

YD,&8)

i=o tr +

EzPl2(Cdt,~o~~

J0

1

c ) 9 0 max,hn)

+

G(C6(t,

where Py, 7 = 0,7,8, a are given positive numbers. Note that here we use the following notation: <3(t, = 62

Ilo; z(t),

6,)

=

4(t),

<8(&m;

(t, rlo; &)

C4(t,f?o;&) &)

=

r&)/h

=

=

y(t),

9,

<5(t,r]0;

CQ(t,‘iD;&)

and

Cl&

= vo, 6,)

=

70;

cl(t,qo;6,)

6,)

=

ti(t),

610(&~0;&)

Z(t),

SC),

hoax)]

=

/3(t),

b(t,qO;&) =

(56)

dt,

&(t,vo; =

%(t),

6,) v,

= <7(t,

Cll(t,~O;&)

9, 70;

=

6,)

u2(t),

~&)P~.

In this section, the following problem is considered: find an element co E R2N such that F(c”) = 0. Note that an element co E R2N satisfies F(c”) = 0 if and only if the corresponding In this case, l&(t)1 I amax, i = 1,2, for all t E [O,tf] and equations (52), (53) are satisfied. S,+ defined via equations (46), (47), where c = co, is a feasible command strategy. The computation of 6,o was conducted by solving an unconstrained minimization problem on R2N. This was done by using a gradient method described in [8]. However, any other gradient method may be used; see, for example, [9,10]. At each stage, during the minimization process, the function F(c) was computed by solving equations (33)-(45) (using (46), (47) and (49)) on [0, tf]. Equations (33)-(45) were solved by using a fourth-order Runge-Kutta method.

5. PATH

CONTROLLABILITY

In this section, we derive a condition under which the disk’s motion is path controllable. the following functions: cos~cosr$-2(g)

(z)sinC?sind+

Define

(~)2cos*cos# (57)

Y. YAWN AND C. FRANGOS

8

cos*sin4+2 + 921 =

(g)

($)sinBcos$+

(g)lcmBsin$

($) (2) sinA

(42

(58)

+ h3 cos 59)

(101 sin2 0)

(5%

’

- 921 cos 0, sin~cos$+g2rcosgsin$+g22sin~+gr,

(60) and

sinf3sin$-g21cosec0s#-g22cos$+g2,

where hi q, $ , i = 1,2,3 are given in the Appendix. By inserting 9 ( > equations for $$ and 9, and using equations (57)-(62), we obtain

(61) (62)

(equation (29)) into the

Denote

(64)

(65) (66) and define the following transformation:

(67) Then, equation (63) yields

(68) Using the notations ur = z, ~2 = g,

vs = y, ~4 = $!, equation (68) leads to

(69) The system given by equation (69) is controllable [ll]. That is, for any two points A’ and B’ in R4 and any finite time interval [O,tf], one can find a control function v such that y(t) will move from $0) = A’ to y(tf) = B’, where ~(4 = @l(t), a(t), w(t), 4t>lT1 t 2 0.

The Motion of a Disk

9

However, the applied control function is I’ and not v. Thus, since det R = -sin8/(1i1s), equation (67) yields du I’ = o-Q-1 -ag-z+v ) (70) >

in any region that excludes 0 = in, i = 0, f 1, f2, . . . . Hence, the disk’s motion is path controllable in any region of (zi,zz,. . . ,x11) (that is, 212 and 21s are the control functions) that excludes sin21 = 0. Equations (69)) (70) enable one to design closed-loop control laws for the disk’s motion. However, this problem is out of the scope of this work.

6. EXAMPLES In all the examples solved here, the following set of parameters has been used: ai = 0.1 s, i = 1,2; N = 29, a = 0.4m, A = t f / in, where A is the time step used in the Runge-Kutta method and iT = 1400; ri+l - 7i = tf/(iV - l), i = 0,. . . , N - 2; where tf will be given later, cymax - 2g/arsd/s2, emin = r/4, emax = 3X/4, qrnax = 5/arad/s, E7 = Eg = 0.05 m, PO = 102, P7 = Ps = 102, and Pa = lo*, and no = (~,0,~,0,0,3.2/a,0,0,a,O,O,O,O). The following examples have been solved here. For the definition of u,, wc, see equations (42), (43). CASE a. tf = %s, zD = -lam,

90 = 10m.

Example a.1: Us = 0.005m/s2, ws = 2rad/s. Example a.2: u, = 0.005m/s2, ws = 4&/s. Example a.3: u, = 0.01 m/s2, wc = 4 rad/s. CASE

b. tf = 14s XD = -5X-q 30 = 5m. Example b.1: ua = 0.01 m/s2, wc = 8 rad/s. Example b.2: ua = 0.05m/s2, wo = 8&/s.

CASE c. tf =

Example Example Example Example Example

7s, XD c.1: u, = c.2: ua = c.3: Us = c.4: ue = c.5: ua =

= -2.5 m, 9~ = 2.5 m. 0.05 m/s2, we = 16 rad/s. 0.1 m/s2, we = 16 rad/s. 0.5 m/s2, we = 16 rad/s. 2m/s2, wc = 16rad/s. 4m/s2, wc = 16rad/s.

For all cases solved here, the procedure described in Section 4 was applied until the value of F(c) reached the value of zero in double precision. Thus, all the required goals have been met while satisfying the imposed constraints. The numerical study conducted via the solution of the above-mentioned examples suggests that the motion of the disk is more sensitive to an increase in the value of we than to an increase in the value of ua. Some of the results obtained are presented in Figures l-10. 7. CONCLUSIONS

AND

REMARKS

The control structure proposed here depends mainly on the “parameterization” procedure given by equations (46), (47) and (49). This “parameterization” procedure is not unique and one may choose other procedures as well. Also, any other gradient method for minimizing F(c), different from the one used here, can be applied. For example, see (121. Another important feature of the method proposed here is that only feasible command strategies, for which all the specifications of the system and the required goals are met, are computed. That is, here we look for an element co E R2N such that F(c”) = 0. By this, we circumvent the use of, optimal control, and avoid problems which arise in optimal control, such as the problem of existence of a local or a global minimum, or the solution of a two-point boundary-value problem. Also, once a problem has been solved for a prescribed maneuver, then only the values of co have to be stored.

10

Y. YAVINAND C. FFXANGOS 2.4

I

1

I

1

2

3

1

2

3

I

1

I

2

1.8

1.4

1.2

1

0

4 5 6 time (s) Figure 1. The values of e(t), t E [0, tr], for Example cl (solid line) and for Example c.5 (broken line).

7

5

4

3

2

1

G +

0

1 ::

-1

-2

-3

4

0

4

5

6

time (5) Figure 2. The WAITSof q, ple c.5 (broken line),

t E [0, tr], for Example c.1 (solid line) and for Exam-

7

The Motion of a Disk

11

6

5.5

5

4.5

4

2.5

0

1

2

4

3

5

7

6

time (s) Figure 3. The values of b(t), t E [0,tf],for Example ple c.5 (broken line). 10

I

I

c.1 (solid line) and for Exam-

1

I

I

I :: I: 1: II / \ : I

8

6

0

1

2

4

3

5

6

time (s) Figure 4. The values of 9, ple c.5 (broken line).

t E [0, tf], for Example

c.1 (solid line) and for Exam-

7

Y. YAWN AND C. FFLANGOS

12 f

1 1

0

f

I

I

I

I

I 2

! 3

1 4

I 5

t 6

I

time (s) Figure 5. The values of v, Example c.5 (broken line).

3.5

I

I

t E [0, tf],

!

I

for Example c.1 (solid line) and for

I

I

2

3

I

I

4

5

:___-----___“_-~------____

3

2.5

2

1.5

I

0.5

0

L

-3

-2

-1

0

1

X Figure 6. The values of y as function of z for Example c.1 (solid line) and for Example c.5 (broken line).

6

The Motion of a Disk Example c.1

1.2

1

0.6 i : :

0.6 / : 0.4

0.2 I

0 _...

-0.2 0

1

2

3

time (s)

4

5

6

7

Figure 7. The values of 61(t) (broken line) and zlz(t) (solid line), t E [O,tf], for Example c. 1. 1.5

Examplec.1 I

I

I

t 1

I 2

t 3

8

I

1

1

0.5

0

-0.5

-1

-1.5

-2

-2.5

-3 0

I 4

I 5

I 6

time (5) Figure 8. The values of 62(t) (broken line) and 213(6) (&id line), t E IO,+], for Example c. 1.

‘ 7

14

Y. YAWN

AND CFRANWS ExamPIe

CL5

0.6

0.4

0.2

0

,........,,...,...,....,....,..~."....,.......~.*.*....".,-.,.,,...........,...

-0.2

0

I

I

I

I

I

I

1

2

3

4

5

6

Figure 9. The values of 61(t) Example c. 5.

tima fs)

(broken line) and x12(t)

{solid line), t E 10, tf],

-

7 for

Example c.5 1.5

I

1

I

I

1

0.5 0

,...,..

I

*..**”

,,............,..

. ...,.....

. . .

...”

. . .

. . . .

. . . .

. . . . . . . . .

. .

-0.5

-1

+1.5

-2

-2.5

-3

0

1

t

I

I

I

t

1

2

3

4

5

6

Figure 10. The values of 62(t) Example c.5.

time (s)

(broken line) and zyg(t) (solid line), t E [O,tp], for

7

15

The Motion of a Disk

APPENDIX

hl (q,$)

= -IDI

(~)2sinOcosO+I~

+mDasine hz

q, 2 (

=

2(1Dl

-

[(%)

13)

)

+mDacOse h3 (q,z)

= --I3 ($)

l

(z cos4+

(%)

sin+]

sinecose-1D3

(g(s)

[(%)

+ (g)cosQ)

sin&

($sine-mDa2

(%)

cos$] ,

kmDa[(%)sin$-(%)cos$].

(71)

, (2)

(g)

($)sinB+mDgacOsQ

(z)

sine (72)

(s)sinO (73)

REFERENCES 1. E.T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, Cambridge University Press, Cambridge, UK, (1917). 2. J.C. Latombe, Robot Motion Planning, Kluwer Academic Publishers, Boston, MA, (1991). 3. L.S. Pontryagin, V.G. Boltyansky, R.V. Gamkrelidze and E.F. Mishchenko, The Mathematical Theory of Optimal Processes, Interscience Publishers, John Wiley and Sons, New York, (1962). 4. A.E. Bryson and Y.C. Ho, Applied Optimal Control, Blaisdell Publishing Company, Waltham, MA, (1969). 5. D.E. Kirk, Optimal Control Theory, Prentice-Hall, Englewood Cliffs, NJ, (1970). 6. Y. Yavin and C. Frangos, Open loop strategies for the control of a disk rolling on a horizontal plane, Computer Methods in Applied Mechanics and Engineering (to appear). 7. Y. Yavin and C. Frangos, Feasible strategies for the control of a disk rolling on a moving horizontal plane, Mathl. Comput. Modelling 20 (12), 81-95 (1994). 8. Y. Yavin and C. Frangos, Computation of feasible control trajectories for the navigation of a big ship around an obstacle in the presence of a sea current, Mathl. Comput. ModelZing 21 (3), 99-117 (1995). 9. J.A. Snyman, A new and dynamic method for unconstrained minimization, Applied Mathematical Modelling 6, 449-462 (1982). 10. J.A. Snyman, An improved version of the original leapfrog dynamic method for unconstrained minimization: LFOPl(b), Applied Mathematical Modelling 7, 216-218 (1983). 11. B. Friedland, Control System Design, McGraw-Hill, New York, (1987). 12. M.S. Bazaraa, H.D. Sherali and C.M. Shetty, Nonlinear Programming Theory and Algorithms, John Wiley, New York, (1993).