Stability domain determination for nonlinear dynamical syste

Abbate, J. V.

Geiss, G. R.

Mc Gill, R.

English

NASA Technical Reports Server

https://ntrs.nasa.gov/search.jsp?R=19670004280 2020-03-16T17:13:35+00:00Z
Grumman Research Department Memorandum RM-348 
STUDY ON DETERMINING STABILITY DOMAINS 
FOR NONLINEAR DYNAMICAL SYSTEMS - I1 
J. V. Abbate 
G. R. Geiss 
R. M c G i l l  
Systems Research Section 
Grumman Aircraft Engineering Corporation 
Bethpage, New York 11714 
Quarterly Progress Report No. 2 
Covering t h e  Period 1 August 1966 t o  31 October 1966 
Prepared under Contract No. NAS 8-20306 
f o r  
National Aeronautics and Space Administration 
Astrodynamics and Guidance Theory Division 
Aero-Astrodynamics Laboratory 
Huntsville, Alabama 
George C. Marshall Space F l igh t  Center 
15 November 1966 
Approved by: 
Charles E. Mack, Jr. 
Director of Research 
I .  
--. ABSTRACT 
T h i s  second quarter ly  progress repor t  discusses r e s u l t s  ob- 
ta ined during t h e  period August 1, 1966 t o  October 31, 1966 under 
Contract NAS 8-20306, "Study on Determining S t a b i l i t y  Domains f o r  
Nonlinear Dynamical Systems ." 
experimental r e s u l t s  obtained via  the  algorithm developed during 
the  previous quarter ;  reformulation of t h e  algorithm t o  avoid 
computational d i f f i c u l t i e s  encountered during t h e  experiments; 
development of a parameterization of t h e  set of pos i t i ve  d e f i n i t e  
n x n 
closed loop guidance. 
of t h e  contract  a r e  outlined. 
Pa r t i cu la r  i t e m s  discussed are: 
matrices; and analysis  of a system w i t h  time dependent 
Plans f o r  research during t h e  las t  phase 
TABLE OF CONTENTS 
I t e m  -
Introduct ion .............................................. 
Experimental Results ...................................... 
Reformulation of t h e  Numerical Algorithm ............. 
Parameterization of t h e  Set of 
Posi t ive Defini te  Matrices ........................... 
S t a b i l i t y  Analysis of a System w i t h  
T i m e  Dependent Closed Loop Guidance .................. 
Plans f o r  Future Work ..................................... 
References ................................................ 
iii 
Page 
1 
2 
2 
4 
8 
17 
18 
- INTRODUCTION 
T h i s  second quar te r ly  report  discusses t h e  work ca r r i ed  out 
during t h e  period August 1, 1966 t o  October 31, 1966 under Contract 
NAS 8-20306, "Study on Determining S t a b i l i t y  Domains f o r  Nonlinear 
Dynamical Systems." 
voted t o  continuing t h e  experimental inves t iga t ion  of the  algorithm 
described i n  Ref .  1, reformulating the  algorithm t o  circumvent 
observed computational d i f f i c u l t i e s ,  developing a parameterization 
of t h e  set of pos i t i ve  de f in i t e  
a system w i t h  t i m e  dependent closed loop guidance. 
Our a c t i v i t i e s  during t h i s  period were de- 
n x n matrices, and analyzing 
1 
EXPERIMENTAL RESULTS 
I 
I 
I 
1 
1 
I 
A s  i n  t h e  previous quarter ly  repor t  (Ref. l), w e  w i l l  only 
r epor t  on t h e  qua l i t a t ive  aspects of t h e  r e s u l t s  of our numerical 
experiments, During t h e  past quarter  w e  devised a subroutine f o r  
obtaining an i n i t i a l  search point on t h e  cons t r a in t  locus 9 = 0. 
This a l l ev ia t ed  some of t h e  s e n s i t i v i t y  of t he  success of t h e  
computation t o  t h e  i n i t i a l  search point. 
of t h e  Van der  Pol equation, t h e  locus has two p a i r s  of 
branches t h a t  are not symmetric, and t h e  success of the  computa- 
t i o n  depends on which p a i r  of branches t h e  i n i t i a l  point i s  on. 
In an attempt t o  circumvent t h i s  s i t ua t ion ,  w e  have reformulated 
t h e  function t o  be minimized. (This reformulation i s  described 
i n  the  following section.) 
However, i n  ( t h e  case 
$ = 0 
W e  a l s o  invest igated whether  using the  quadratic p a r t  of the 
der iva t ive  of t h e  optimal quadratic V as the  a r b i t r a r y  function 
i n  Zubov's equation would yield an improved estirnate of t h e  domain 
of a t t r a c t i o n  via t h e  fourth order so lu t ion  t o  Zubov's equation. 
The r e s u l t  was negative f o r  t he  Duffing equation and pos i t i ve  f o r  
t h e  Van der  P o l  equation. 
Reformulation of t h e  Numerical Algorithm 
I n  R e f s .  1 and 2 an algorithm f o r  computing an estimate of 
t h e  domain of a t t r a c t i o n  of a quasi- l inear  dynamical system was 
given, 
funct ion parameters a (the elements of Q) ,  compute x such 
t h a t  V(x,cx) i s  minimum on $(x,cx) = 0 (i.e., 
The algorithm consisted of two s teps:  1) given a set of 
2 
I .  
i s  used t o  solve t h i s  s tage  numerically); and 2) f ind  a such 
t h a t  t h e  f i gu re  described by 
volume (i.e., f i nd  a" such t h a t  
V(x,a") = , t (a") has  t h e  l a rges t  
/ 
( d e t  P(ao))' (de t  P(a))* 
f o r  a l l  allowable a). 
c a l l y  by describing t h e  Q-matrix i n  a parametric form, i n  
terms of a ,  
n i t e  (given appropriate r e s t r i c t i o n s  on a )  and by computing 
,!(a) from Eq. (1) and P(a)  from 
This second s t ep  i s  ca r r i ed  out numeri- 
which guarantees t h a t  t h e  matrix i s  pos i t i ve  def i -  
T 
A P(a )  + P ( a ) A  = - Q(a) , (3) 
where 
equations . 
A i s  t h e  matrix of t h e  l i n e a r  p a r t  of t h e  dynamical 
T h i s  two s tage optimization, i n  wbich one computes suc- 
cess ive ly  b e t t e r  values of x and a ,  does not work w e l l  f o r  
t h e  Van der P o l  problem. Apparently, a relative i n s e n s i t i v i t y  
of t h e  area of t h e  estimate t o  a 
branches (two symmetric pa i r s )  of t h e  locus o(x,a) = 0 
cause t h e  observed d i f f i c u l t i e s  i n  t h e  computation of an 
estimate of t h e  domain of a t t rac t ion .  I n  order t o  alleviate 
these d i f f i c u l t i e s  and t o  make t h e  algorithm more e f f i c i e n t ,  
w e  have reformulated it as a s ing le  minimization over 
a ,  simultaneously. 
and t h e  existence of four  
x and 
Consider t h a t  t h e  problem i s  t o  compute a constrained mini- 
mum of V relative t o  x and a constrained minimum of the  area 
r e l a t ivz  t o  a. This can be accomplished by recognizing t h a t  
I 3 
b 
t h e  content of t h e  f igu re  described by V(x,a) < j ( a )  i s  pro- 
1 por t iona l  t o  ( , ! ? ( c x ) ) ~ ’ ~  (de t  P(a) -2 . Then we  can compute e i t h e r  
1 ! 
o r  
where w e  have used t h e  penalty formulation f o r  handling t h e  con- 
s t r a i n t s  +(,,a) = 0 and Q pos i t i ve  d e f i n i t e  (accounted f o r  by 
g(a)).  Cursory examination ind ica tes  t h a t  x = 0, t h e  t r i v i a l  
solut ion,  i s  t h e  global solut ion of Eq. (5) ;  t h i s  i s  not true f o r  
Eq. ( 4 ) ,  and thus only Eq. ( 4 )  i s  considered fur ther .  
A computer program h a s  been wr i t t en  t o  compute Eq. (4 )  for 
the  Van der Pol equation by using t h e  Min-All algorithm developed 
by M c G i l l  and Taylor. 
da t a  t o  draw conclusions about t h i s  formulation a t  t h i s  t i m e .  
Unfortunately, w e  do not have s u f f i c i e n t  
- Parameterization of t h e  Set  of Pos i t ive  Defini te  Matrices 
The s t a b i l i t y  analysis  algorithm described i n  Refs. 1 and 2 
requi res  t h e  generation of pos i t ive  d e f i n i t e  
t h e  ana lys i s  of 
The b ru te  force approach i s  t o  form an n x n symmetric matrix 
and then apply t h e  determinantal t es t  (Ref .  3) t o  determine i f  i t  
n x n matrices f o r  
n-dimensional quasi-l inear dynamical systems, 
i s  pos i t i ve  def in i te .  
matrix elements and then t h e  evaluation of t h e  de- of 
terrninants of t h e  n-principal minors of t h e  matrix. However, it 
This procedure requires  the  a r b i t r a r y  choice 
n n + l  
2 
4 
does not  provide information on how t o  co r rec t  a candidate matrix 
t h a t  f a i l s  t h e  t es t  f o r  posi t ive def ini teness .  Therefore i t  would 
be des i rab le  t o  generate t h e  matrix by a procedure t h a t  guarantees 
t h e  matrix i s  pos i t ive  de f in i t e  and spans t h e  e n t i r e  set of pos i t i ve  
d e f i n i t e  matrices. 
based upon t h e  work of Murnaghan (Ref. 4 )  on t h e  parameterization 
of t h e  group of un i t a ry  matrices. 
It i s  w e l l  known (Ref. 3) t h a t  a l l  real  symmetric'matrices are 
I n  t h i s  section we  develop such a procedure 
orthogonally s i m i l a r  t o  a diagonal matrix, and t h a t  a l l  pos i t i ve  
d e f i n i t e  (pd) matrices are then orthogonally s imilar  t o  a diagonal 
matrix w i t h  pos i t ive  diagonal elements; i.e., l e t  P be pd, then 
T P = S A S  , 
where 
hi > 0 , i = 1, 2, ...I n 
T S S = I  . 
Thus t h e  parameterization of a l l  pd matrices P i s  
reduced t o  t h e  parameterization of t h e  group of orthogonal 
matrices S. 
I n  Ref. 4 Murnaghan proves t h a t  t h e  parameterization 
of t h e  group of n x n uni tary matrices U i s  accomplished 
by t h e  f ac to r i za t ion  
5 
(7) 
where 
n-1 
.t=k+l 
+ (n - 1) , 2n - k - 2)(k - 1) I-L = (  2 
[uii = 1 , i f k , j  
u = + e+i0 s i n  e , i Jk 
The fac tor iza t ion  of t h e  group of orthogonal matrices i s  i m -  
mediately obtained by requiring U t o  be r e a l ;  i.e., 6 = 0 = 0 , 
'n 
K 
= t K , -K 5 cpk < 7~ k # n , and -3 < 8 < 9 In  pa r t i cu la r ,  
6 
n-1 
k=l  
n-1 
- 2, thetas and 2 T h i s  f ac to r i za t ion  contains ( + 1 parameters. The n-phis, o r  a t o t a l  of - 
n-lambdas i n  Eqo (7) r a i s e  the number of parameters t o  
nu + 1 - one more than required. 
t h e  number of parameters w i l l  be 
quired t o  represent an a rb i t r a ry  symmetric matrix. 
'n 
wish t o  r o t a t e  and scale the  e l l i p s o i d  associated w i t h  the  
mat r ix  and do not want t o  reflect  coordinates o r  change t h e  
handedness of the coordinate system. 
n l n  - 1) 
2 
Thus i f  w e  re- 2 
str ict  S t o  be a ro t a t ion  matrix (i.e,,, choose 'pn = 01, 
n n + l  u;the  number re= 2 
The choice 
= 0 i s  i n t u i t i v e l y  motivated by t h e  consideration t h a t  w e  
pd 
The f ac to r i za t ion  of a pd matrix of dimension t h r e e  i s  
thus given by 
T P = S A S  , 
where  
7 
(15) 
1 
I * -  
1 
I 
I 
I 
1 
I 
I 
I 
I 
I 
I 
I 
8 
I 
I 
I 
I 
and 
il 
s12 = I"' 1 
0 
'$ 
0 o \  
', 
I 
-SO 1 0 '\ 
0 l /  
O C T l ,  i 
The fac tor iza t ion  given here  is similar t o  t h a t  conjectured i n  
Ref, 1, 
S t a b i l i t y  Analysis of a System w i t h  Time Dependent Closed Loop 
Guidance 
I n  t h e i r  simplest  fora, booster guidance control  systems are 
described f o r  f i n i t e  in te rva ls  of t i m e  by nonlinear, nonautonomous 
d i f f e r e n t i a l  equations. In  t h e  more general  case, t h e  cont ro l  laws 
8 
I 
I * *  
1 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
a r e  determined by i t e r a t i v e  procedures and cannot be expressed 
as simple functions of the s ta te ,  Because t h e  standard w e l l -  
known techniques of s t a b i l i t y  analysis  a r e  not readi ly  appl icable  
t o  such systems, new techniques, o r  var ia t ions  of present Liapunov 
methods a re  being sought. 
I n i t i a l  research i n  t h i s  area has consis ted of determining 
t h e  s t a b i l i t y  of a s imple  control  system t h a t  contains  character-  
i s t i c s  representat ive of t he  more complex booster guid'ance systems. 
Specif ical ly ,  w e  a r e  concerned w i t h  guiding t h e  motion of a p a r t i c l e  
moving a t  constant speed i n  a plane i n  t h e  presence of a constant 
disturbance. The dynamics of t h e  system a re  given by 
Gl = v cos u 
f2 = v0 + v s i n  u , 
where 
r e l a t i v e  t o  t h e  disturbance, V i s  t h e  disturbance and u 
i s  t h e  d i rec t ion  of t h e  velocity.  The control  l a w  u ( t )  i s  
such t h a t  t h e  p a r t i c l e  i s  guided from t h e  i n i t i a l  point  
face  of disturbances i n  t h e  i n i t i a l  conditions. 
v i s  the magnitude of t h e  ve loc i ty  of t h e  p a r t i c l e  
0 
(xl,x2) 0 0  t o  t h e  f i n a l  point (xl,x2) f f  i n  minimum t i m e  i n  t h e  
The problem 
h a s  been made more spec i f i c  by l e t t i n g  
x1 = x2 = 0 
v = 1, Vo = 1/2 ,  
0 0 f f and x1 = 2, x2 = 1. 
I n  t h e  absence of any disturbance i n  the i n i t i a l  conditions,  
w i t h  the  
* t h e  optimum control  law f o r  minimum t i m e  i s  u (t) = 0, * * corresponding t r a j ec to ry  x1 = t, x2 = t / 2 ,  and nominal t i m e  
Tf = 2,  
buted w i t h  a b iva r i a t e  normal d i s t r ibu t ion  of errors i n  i n i t i a l  
condi t ions w i t h  mean value 0 and standard deviat ion 0.1, The 
con t ro l  law i s  assumed t o  be l inear ,  time varying, and of t h e  form 
The i n i t i a l  disturbance i s  assumed t o  be randomly d i s t r i -  
9 
For t h e  spec i f i c  problem considered, with 
t h e  mean square m i s s  of t h e  ta rge t  point  (2 , l )  
t h e  r e su l t i ng  expected f i n a l  t i m e  i s  2.00287. 
i s  0.00139; 
The s t a b i l i t y  problem for t h e  system described above may 
be s t a t e d  as follows: From w h a t  set of i n i t i a l  states (xy,xy) 
w i l l  the  p a r t i c l e  reach a point i n  some €-neighborhood of t h e  
f i n a l  s t a t e  (2,1), L e . ,  
subjec t  t o  t h e  cons t r a in t  Tf 5 M, where M i s  some constant. 
I n  order t o  determine any characteristic of the  problem 
t h a t  might be usefu l  i n  t h e  development of a general  s t a b i l i t y  
technique, an ana lys i s  of the system response w a s  undertaken, 
It w a s  i n i t i a l l y  assumed tha t  t h e  con t ro l  l a w  u ( t )  w a s  time 
invar ian t ,  i.e., 
Thus t h e  system i s  described by I * '  
Under t h e  t r ans l a t ion  
Eqs. (25) b me 
1 + cos $1 = - 
= s i n  'p y 92 L 10 1 
with the  equilibrium solut ion y1 = y2 = 0. 
The loci of equilibrium points  f o r  Eqs. (27) are given by 
l ies on t h e  Thus i f  t h e  i n i t i a l  point  (yl 0 = x1 0 , y2 0 = x2> 0 
locus,  t h e  solut ion f o r  y,(t) and y2 ( t )  w i l l  be constant 
i n  t i m e ,  and t h e  t r a j e c t o r i e s  i n  the  xl,x2-plane w i l l  be given by 
X l ( t )  = t + X1(O) 
x 2 ( t )  = t / 2  + x2(0) 
(29) 
0 
11 
These t r a j e c t o r i e s  are p a r a l l e l  t o  t h e  optimum t r a j ec to ry  and 
do not  converge t o  t h e  target  point. 
To f ind  t h e  complete t r a j e c t o r i e s  i n  t h e  yl,y2-plane, l e t  
then from Eqs.  (27) 
2 ( t )  = ayy(t)  + byy(t) = - a + a cos z + b s i n  z . (31) 
Solving f o r  z ( t )  , 
-1 z ( t )  = 2 t a n  
where 
p20 1 -- 
1 - ke -P20t 
By subs t i t u t ing  Eq. (32) i n t o  E q s ,  (27), and in tegra t ing ,  
1 2  
Y,(t) 
It i s  
i s  on 
1 - 1) - tan" 2 (k - 1) - 2p20 tan" (ke -P20t p20 2 - [ p20 - 2  p10 + p20 
I n  p20 
p10 + p20 
- 2  2 
(34 Cont.) 
0 
evident from Eqs. (33) and (34) t h a t  i f  t h e  i n i t i a l  point 
t h e  equilibrium locus, Eqs. (34) reduce t o  
Y p  = Y p )  , Y p  = Y2(0) 
It can a l so  be seen from Eqs. (34) o r  Eq. (32) t h a t  i f  t h e  
i n i t i a l  point (yl(0), ~ ~ ( 0 ) )  does not l i e  on the  equilibrium 
i s  r e a l  and negative, then the  f i n a l  point locus,  and i f  p 20 
(yl(m), y2(w), \ w i l l  l i e  on t h e  equilibrium locus. Thus i f  
i s  real  and negative, every solut ion y ( t )  tends t o  t h e  p20 
set  of points  where 
An approximate analysis  of t h e  response of t h e  system when 
t h e  cont ro l  l a w  i s  t i m e  varying, as given by Eqs. (20) and (21) ,  
can be ca r r i ed  out i f  p,(t) and p2( t )  are approximated by 
staircase" functions such t h a t  I 1  
13 
P 2 W )  = P20 + P 2 p T  
t h  I n  t h e  n in te rva l ,  i .e.,  (n - l )T < t < nT 
for 0 5 z 5 T o  The constant kn i s  given by 
L 
By using Eqs .  (36) and (37), t he  following l i nea r  difference 
equation f o r  kn r e s u l t s :  
where 
Solving Eq. (38) yie lds  
14 
. 
It . *  
I * *  
I 
i R k n = if Fi]-ko + - 
1 f i F j  
1 i 
I 
o r  
n c 
i=l 
The approximate locat ion of t h e  p a r t i c l e  a t  any i n s t a n t  of 
t i m e  i s  then specif ied by 
yin(T) = Pin('t) + y. (0) + ), P . . ( l )  , i = 1, 2 , (41)  10 1.l 
j =I 
where 
- 2P2 
e 
15 
The work car r ied  out t o  date has e s sen t i a l ly  consisted of 
analyzing t h e  given system and has  resul ted i n  Eqs. (41) and (42). 
It w i l l  now be necessary t o  study these equations along w i t h  t h e  
t r a j e c t o r i e s  t h a t  r e s u l t  from a computer simulation i n  order t o  de- 
termine cha rac t e r i s t i c s  t h a t  w i l l  enable development of a s t a b i l i t y  
analysis  technique t h a t  does not require  a solut ion of t h e  system 
equations. To achieve t h i s ,  our e f f o r t s  i n  the next quarter  w i l l  
be directed toward examining the p o s s i b i l i t y  of applying standard 
Liapunov s t a b i l i t y  techniques over each in t e rva l  (n - 1)T< t 2 nT, 
and using t h e  information obtained from t h e  sub-intervals  t o  con- 
clude s t a b i l i t y  f o r  t h e  e n t i r e  i n t e rva l  0 < t 2 Tf. 
7.6 
1. 
2. 
3. 
4. 
P M S  FOR FUTURE WORK 
The  emphasis during t h e  las t  phase of t h i s  cont rac t  w i l l  be 
upon formulating other examples t h a t  are representat ive of t h e  
booster guidance s t a b i l i t y  problem. There are, however, other 
p a r t i c u l a r  areas t h a t  w i l l  receive some a t t en t ion ,  v iz , :  
Formulation of a function whose global  
minimum i s  achieved f o r  t h e  values of 
s ta te  x and parameters a t h a t  produce 
t h e  optimal quadratic estimate of t h e  do- 
main of a t t rac t ion .  
Development of  a procedure f o r  successively 
improving a given quadratic estimate by 
modification of t h e  Liapunov function via 
geometric considerations. 
Development of an e f f i c i e n t  procedure f o r  
solving t h e  Liapunov matrix equation 
(Eq. (4 )  of Ref. 2) .  
Numerical invest igat ion of t h e  guidance ex- 
ample discussed earlier,  and an ana ly t i ca l  
invest igat ion of whe the r  s t a b i l i t y  r e s u l t s  
over segments of t h e  i n t e r v a l  of operation 
of t h e  nonautonomous system can be combined 
i n t o  a r e s u l t  for  t h e  e n t i r e  interval of 
concern. 
1 7  
REFERENCES 
1. Geiss, G. R., Study on Determining S t a b i l i t y  Domains f o r  
- Nonlinear Dyn amical Systems, Research Department Memorandum 
RM-334, Grumman Aircraf t  Engineering Corporation, August 1966. 
2. Geiss, G. R., Estimation of t h e  D,omain of Attract ion,  Research 
Department Memorandum RM-3165, Grumman Aircraf t  Engineering 
Corporation, March 1966. 
3 .  Gantmacher, F. R., The Theory of Matrices, Vol. I, Chelsea 
Publishing Co., New York, 1959, pg. 337, Theorem 19. 
4 ,  Murnaghan, F. D., Lectures on Applied Mathematics, Vol. 111: 
The Unitary and Rotation Groups, Spartan Books, Washington, 
D.C., 1962. 
18 


Study on determining stability domains for nonlinear dynamical systems, II Quarterly progress report, 1 Aug. - 31 Oct. 1966

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server