https://ntrs.nasa.gov/search.jsp?

R=19710018599 2019-05-25T03:20:58+00:00Z

N A S A TECHNICAL NOTE
"
NASATN D-6218
e,/

A METHOD FOR REDUCING

THE SENSITIVITY OF
OPTIMAL NONLINEAR SYSTEMS
TO PARAMETER UNCERTAINTY

bY
JurrellR.Elliott
Lungley ReseurchCenter
U?Zd
Willidm F. T e u p e 1

University o f Kunsus . >

. .I

..I

N A T I O N AALE R O N A U T I CASNSDP A CAED M I N I S T R A T I O N WASHINGTON,

D. C. JUNE 1971

' l l l l l l I I Ill1
 TECH LIBRARY KAFB. NM

. . , ~ ""~
0233005
1. Repwt No.
NASA TN D-6218
4. Title and Subtitle
.
.. "~ 1 2. Government
Recipient's
Accession
No.

1 5. Report Date
Catalog No.

June 1971
A METHOD FOR REDUCING THE SENSITIVITY O F OPTIMAL I
NONLINEAR SYSTEMS TO
PARAMETER
. -
UNCERTAINTY
. "_ " "" .
I
1
6. Performing Organization Code

7. Author(s1

..
J a r r e l l R. ElliottandWilliam
~ - " .
9. Performing Organization Name and Address
.. .
"
F. Teague(University
-. --
I
of Kansas)L-7485. .
8. Performing Organization Report No.

rO:-Work Unit
. - .
No.
- ..

125-17-06-03
.. ". ~ ~

NASA Langley Research Center

Hampton,Va.23365 I 11. Contract or Grant No.

. " -
2. SponsoringAgencyNameandAddress
~ _" ~ " .

[-" 13. Type of Report and Period Covered

". TecgicalNote
-

"
National Aeronautics and Space Administration
Washington, D.C. 20546
. . -~
~
5. SupplementaryNotes
~" - ~-
- - L 14. SponsoringAgencyCode

"" " ..

___ ~ ~ -
~
I
~~
_

6. Abstract

Mathematical relationships are derived and used to establish a procedure for reshaping
the optimal solution so as to reduce the statistical uncertainty in the terminal conditionsof the
system due to system parameter uncertainties of known statistical properties. The procedure
introduces the use of an augmented performance index which contains a scalar measure of the
system sensitivity partial derivatives. A nonlinear multiparameter optimal-rocket-trajectory
problem was solved by using an algorithm based on the method of steepest descent to illustrate
the procedure.

" ~ ~
" .
17. Key Words(Suggested by Authoris) ) 18. Distribution Statement

Sensitivity Unclassified - Unlimited

Optimal control
Nonlinear systems

19. Security Classif.this

report)
(of20. Security Classif. (of this page) 21. NO. of Pages 22. Price'

Unclassified - - ..
Unclassified
~~~ . ~~
~. " ~~ ~ -

'For sale by the National Technical Information

Service, Springfield, Virginia 22151
 A METHOD FOR REDUCING THE SENSITIVITY
OF OPTIMAL NONLINEAR SYSTEMS
TO PARAMETER UNCERTAINTY

By Jarrell R. Elliott
Langley Research Center
and
William F. Teague
University of Kansas

SUMMARY

The parameters of a nonlinear dynamical system which is to be controlled optimally

a r e not always accurately known. It may therefore be desirable to accept a reduction in
the predicted nominal performance of the system in exchange for the ability to better pre-
dict the outcome of the system, or plant, operation.
In this paper relationships to predict mathematically the sensitivity of the system to
parameter errors are derived and used to establisha procedure for reshaping the optimal
solution to reduce the statistical uncertainty in the terminal conditions of the system due
to known statistical characteristics of the system parameters. The procedure requires
the introduction of an augmented performance index which is a linear combination of the
original performance index and a positive scalar measure of the system sensitivity com-
posed of the weighted s u m of the variances of the performance index and terminal con-
straints. The augmented performance index, the sensitivity partial derivatives, and the
original state variables are used inthe formulation of a new, higher dimensioned optimi-
zation problem of the same form as the original problem. The new problem introduces
certain weighting factors which permit different relative importance to be attached to
different types of sensitivity, such as position relative to velocity, and which allow for
adjustment of performance degradation and of sensitivity reduction.
The procedure developed was illustrated by solving a nonlinear multiparameter
rocket-trajectory problem. An algorithm based on the method of steepest descent was
used to solve the problem because of the widespread use and proven versatility of this
numerical technique. The example solutions serve to show the tradeoffs made possible
by changing the weighting factors and to illustrate the radically different solutions one
can obtain when sensitivity considerations are included in the problem formulation.
 INTRODUCTION

In calculating the open-loop control time history of a physical process or plant, it is

common to assume that the plant behavioro r outcome can be predicted with the aid of a
mathematical model of the plant with known or deterministic values of the plant param-
eters. More often than not, however, the plant parameters are stochastic, not determin-
istic, and serious discrepancies between predicted outcome and actual outcome may occur
as a result of plant-parameter variations of one of two types: (1)parameter variations
during plant operation, o r (2) errors in the estimates of fixed plant parameters. These
effects (that is, the sensitivity of a dynamical process to parameter variations) should be
a consideration in the design of any control system.
While many studies related to plant sensitivity have been reported in the literature
(see refs. 1 to 14 for a nonexhaustive list), most of them deal with linear problems. It
also appears that little effort has been devoted to computational aspects (a notable excep-
tion is ref. 14). The present paper deals with the stochastic nature of the parameters in
nonlinear problems and shows how a well-known computational algorithm in optimization
theory may be applied to obtain numerical results.
The problem considered is a multivariable, multiparameter optimal control problem
(with terminal constraints) whose mathematical model consists of a s e t of ordinary first-
order nonlinear differential state equations. The analysis is limited to parameter varia-
tions of type (2). The sensitivity of the performance index and the sensitivity of terminal
constraints to parameter variations, that is, the partial derivatives of these quantities
with respect to the parameters, areused to construct a s c a l a r m e a s u r e of plant sensi-
tivity. This sensitivity measure, which is the expected value of the weighted s u m of the
s q u a r e s of the sensitivity partials, is multiplied by a weighting factor and added to the
original performance index to form an augmented performance index. Minimization of
this augmented performance index, with appropriate weighting factors, results in solu-
tions which are less sensitive to parameter variations.
A simple nonlinear example problem, representative of rocket flight in a uniform
gravitational field, is worked out in detail to show the steps required in setting up and
solving plant sensitivity problems and to illustratehow reducing the sensitivity of a plant
to parameter variations can modify the open-loop control time history and state-variable
time history of the plant. The example also serves the purpose of demonstrating the use
of the cumputational algorithm in solving a typical problem.
This research was conducted a t NASA Langley Research Center, with William F.
Teague in residence under a grant arrangement with the University of Kansas.

2
 SYMBOLS

n X q sensitivity coefficient matrix, af/aa"

partition of matrix A (i=1,2, ...,n)

system parameter vector with components al,.. .,aq

mean (nominal) value of -

a

aerodynamic drag coefficient

aerodynamic drag

"length" of control step

control step-size measure

expectation of random variable ()

n X n state coefficient matrix

partition of matrix F (i=1,2, ...,n)

components of -
f

governing equations of state

governing equations of augmented state

gravitational acceleration

vector elements of augmented-state differential equations (i=1,2, ...,n)

I$+,IQ$,IQQ constant matrices used in steepest-descent algorithm and defined by

equations (A16)

JA augmented
performance index

3
KD modified drag coefficient

L
- vector set of terminal constraint values (in the example problem, L1 and
L2 are the values for altitude and vertical velocity component, respectively)

1 number of terminal constraints

m - vehicle mass
dimension of control vector u;

m0 initial vehicle mass

m time rate of change of m a s s

dimension of state vector -

x

number of augmented-state variables

covariance matrix

dimension of system parameter vector

S "

/
sensitivity matrix, ax aa; frontal surface area

Si partition of matrix S (i=1,2, ...,n)

T vehicle thrust

t time (independent variable)

tf final time

t0 initial time

U horizontal velocity component, 43

-
U m-dimensional control vector

V total velocity magnitude

V vertical velocity component, X4

W m X m (symmetric) control weighting matrix

relative weighting factors(i=1,2, ...,Z+l)

4
 sensitivity weighting factor

n-dimensional state vector

fi-dimensional augmented-state vector

components of -
ii.

range (with numerical subscript denoting a specific state variable)

altitude

flight-path angle

thrust-attitude control angle

n-dimensional adjoint vector associated with payoff

(2 + 1) X n variable matrix integrating factor

2 X n matrix of adjoint variables associated with terminal constraints

constant Lagrange multiplier

I-dimensional constant Lagrange multiplier vector

atmospheric density

sensitivity performance index

performance index (not including sensitivity)

sensitivity measure

2-dimensional terminal constraint vector

components of \I/ (i=1,2, ...,2 )

-

transpose of matrix ( )

A dot over a symbol indicates a derivative with respect to time. The symbol 6
denotes a variation in a quantity, and A denotes a first-order perturbation.

5
 PROBLEM STATEMENT

Consider the following fixed- time open-loop optimal control problem:

Minimize

for the system

with terminal constraints

”
+-L=O

where

@p is a scalar performance index,

-
x is an n-dimensional state vector,

-f is an n-dimensional state derivative vector,

-
u is an m-dimensional control vector,

-a is a q-dimensional parameter vector, whose mean and covariance matrix

(q.>
a r e known = 2 ; E(6a 6aT} = P), and

+
- is an 2-dimensional terminal constraint vector.

(All v e c t o r s a r e column vectors unless superscripted by T which indicates vector o r

matrix transpose, except -.
a vector)’ which is a row vector.)
It is desired to combine some measure of sensitivity G S to GP such that when the
combination is minimized, the performance index +
qP and *.he terminal constraints -
will be less sensitive to perturbations in the parameters - a. It should be realized that
this reduction in sensitivity will degrade the performance index qP. However, since
a system may have more performance capability than is needed, a system user may be
willing to sacrifice some performance in order to reduce terminal perturbations or

6
 P
$b.
e r r o r s i n @ and - In order to provide different degrees of performancelossand
sensitivity reduction, a weighting constant is put on the sensitivity measure @s, and the
sensitivity performance index to be minimized is then defined as

The sensitivity measure should reflect the fact that the system is stochastic. The
sensitivity measure will therefore be taken to be an expected value of the weighted sum
of the squares of the perturbations in @D and each Gi due to parameter variations.
(It is necessary to square the perturbation in GP and in each +i (i=1,2, ...,Z) to insure
that the measure is positive)., Now since the parameter variations are of type (2)
described in the introduction, the perturbations to first o r d e r are

and

(i=1,2, ...,Z)

The squares of these perturbations may be expressed as

and

(i=l,2, ...,Z)

The sensitivity measure contains certain weighting factors, to be called relative weighting
factors, which are required because one A2 quantity may be considered more, o r less,
important than another (for example, one may be more concerned with velocity errors
than position e r r o r s in some application although both are used as terminal constraints).
Finally, the sensitivity measure is
On a trajectory defined by some open-loop control time history, the partial derivatives
inside the expected value are constant. By well-known identities involving the expected
value (see ref. 15), equation (1) may be rewritten as

Since E (6. haT) P, thisequationbecomes

where P is the covariance matrix of the random parameters which is assumed to be

known along with the mean value, that is

This expression for & (eq. (2)) may be recognized as the weighted sum of the vari-
ances of the perturbations in performance and constraints. Finally the sensitivity per-
formance index for minimization is

Now returntotheproblem of computing A+ and Aqi. Sincetheseperturbations are

P
explicitly functions of the terminal state, they may be written as

]
(i=1,2, ...,2)
(4)

The remaining problem is that of determining

ax
aa
-
-
(tf). In reference 16 it is shown that
for parameter variations of type (2) and f o r the system differential equations
ax
x = f(x,s(t),a,t) = f(x_,a,t),
- $t) is the solution of

8
By defining

and

aa

equation (5) may be written as

d S=
- FS +A ; S(t0) = 0
dt

Then @s becomes

A simple, easy-to-work-with form for @s may be determined by reordering and pos-

sibly introducing some new variables such that

and

This procedure may, i f new variables are required, introduce a new system of differential
equations. The notation used will remain unchanged, however, and the system of differ-
ential equations will continue to be called ir = f. "

Now partition S, where SiT (i=1,2, ...,n) denotes the rows of S, as follows:

9
Then
(
$I= x1 + WISITPSl + w2s2 T PS2 + . . . + wl+lsz+lTPSkl) I
tc1
T
Notethat SiT = Fi S + AiT; F and A have been partitioned in the same manner as S,
that is,

FIT
""
AIT
""

FaT A2T
F= ""
; A= ""

"" ""

FnT AnT

Observe that this optimization problem involving the sensitivity index has the same form
as the original problem involving only the performance index (but with more dimensions)
and may be solved by any of several algorithms. Restated, the problem is to minimize

subject to the differential constraints

2 = f = fG&2,t) ; x(to)
- = go

SI = STF1 + A1 = gl(Z,Si7g,B,t) ; Sl(to) = 0

10
and the terminal conditions

k(Z(tf)) = 4

The augmented-state terms are defined as follows:

4

dx
Then the differential constraints are = = f(ii,u,Z,t). Thus it may be seen that the only
dt - - - -
difference between this problem and the original problem is that it is now possible to
control, at the expense of performance index Cp the sensitivity of the trajectory, and
P'
the problem is now a n n ( l + q)-dimensional state-variable problem rather than an
n-dimensional one. Also several weighting factors (ws and ~ 1 . ~..., 2 W, Z + ~ ) have been
introduced into the problem. These weighting factors are to be used in applications to
control the loss in the performance index C#J and to properly emphasize the importance
P
of one type of terminal error relative to another type. The contribution and effect of
these weighting factors will be illustrated in an example problem.
It should be pointed out that solving the differential equations for the sensitivity
partials gives all the information needed to perform a f i r s t - o r d e r e r r o r a n a l y s i s on any
of the state variables of the original system. In the example problem, a comparison will
be made between first-order error analysis results obtained by using the sensitivity
partials and first-order error analysis results for which one-sigma errorswere intro-
duced one a t a time into the equations of motion.
The discussion herein deals with the fixed-time problem. However, fixed time w a s
used for convenience and clarity only and is not a limitation of the formulation.
While almost any of the algorithms for solving optimization problems would be
applicable here, it was decided to use the steepest-descent algorithm as outlined in ref-
erence 16 because of its widespread use and proven versatility. The procedure used in
this algorithm is to linearize about a solution - ?*(t) provided by some reasonably chosen
controltimehistory -u*(t)(which ingeneralneitherminimizes Cp nor satisfies
@ - 4 = 0) andtosolvefor6u(t)(which
- - improves Cp(tf) and Ic/(tf)). A new solution
provided by - u"(t) + 6u(t)
- is obtained and the procedure is repeated successively until
@
- - 4 is sufficiently close to zero and @(tf) can no longer be decreased. This final
solution is said to be optimal although no necessary conditions for optimality have been

11
 satisfied. A brief derivation of the necessary relationships with notation common to
many steepest-descent programs is provided in appendix A. The algorithm was applied
to the following example problem.

A NUMERICAL EXAMPLE

Problem Statement
The example problem is a fixed-time problem in which it is required to determine
the thrust-attitude program of a single-stage rocket vehicle starting from rest and going
to specified terminal conditions of altitude and vertical velocity which will maximize the
final horizontal velocity. The idealizing assumptions made a r e the following:
(1)A point-mass vehicle
(2) A flat, nonrotating earth
(3) A constant-gravity field, g = 9.8 m/sec2 (32.2 ft/sec2)
(4) Constant thrust and mass- loss rate
(5) A nonlifting body in a nonvarying atmosphere with a constant drag parameter
1
KD = -PC$, where S is the frontal surface area.
2
The coordinate system and pertinent geometric relations and t e r m s are shown in
figure 1. The differential equations of motion needed in the algorithm setup a r e

du = - (T c o s 8 - K p V ) = 21 = f l
”

mdt
dY -
”
-v=x2=f2
7
dt

*
dt
1
= -(T sin e - KDUV) - g = k3
J
= f3

where m is the vehicle mass and where

and

m(t) = mo + mt

An equationfordx/dt is notincludedintheequations of motionbecuasexdoes

not enter into the problem. However the equation k = u was integrated separately to
obtainrange.Theequationfordm/dt is analyticallyintegrable,sincedm/dt is a
constant, and is therefore not included in the set of differential equations.

12
 The parameters of the problem, which will be considered fixed constants whose
precise values are unknown, are thrust level, mass-loss rate, initial vehicle mass, and
modified drag coefficient. The parameter vector is defined as

These parameters are assumed to be statistically independent and to have a normal dis-
tribution function with mean and one-sigma values as shown in table I.

TABLE I.- SYSTEM PARAMETER VALUES

__
Parameter Mean value One-sigma value,
fraction of mean valu
Thrust,
T, kN (Ibf) .................... 17.8(4000) 0.0067

Mass-lossrate, m, kg/sec(slugs/sec) .......... -9.1

(-0.62) 0.0167

Initialmass, m,, kg (slugs) . . . . . . . . . . . . . . . . 1433.6

(98.259) 0.010

Modified dragcoefficient, KD, kg/m(slugs/ft) ...... 0.048(0.001), 0.479 (0.01) 0.0167

These values are considered representative of the current state of knowledge i n

(
solid rocketry. The mean values used were chosen to conform to a flight of 100 seconds
for a rocket with specific impulse of 200 seconds, an average acceleration

l J l o o
100 t=O m(t>
x dt of 2g, and a ratio of initial to final mass m(to)/m(tf) equal to e, the
\

base of the natural logarithm.

The boundary conditions a t t = 0 andt = 100 seconds are shown in table II.

TABLE 11.- BOUNDARY CONDITIONS

Variable Initial conditions Terminal conditions

(t = 0) (t = 100 seconds)

i
u,m/sec(ft/sec) .......... Maximum
y, m (ft) . . . . . . . . . . . . . . . 15 240 (50 000)
v, m/sec(ft/sec) ..........

In the notation previously introduced

13

I
 @ - y(tf)
1-
- L1 = 0 ; L1 = 15240 m (50000 ft)

and the vector control variable -u in the general formulation is the thrust-attitude
angle 8 . This completes the problem statement (without trajectory-sensitivity
considerations).

Sensitivity Relations
The sensitivity-matrix differential equation is

@=FS+A
dt

where

i
"
au
am

S=

-
aT
av

and

fl

if3 + 9)
14
 T dSiT
-= Fi S + A iT (i=1,2,3)
dt

or in column vector form

- = S TF i + +
dSi
(i=1,2,3) (16)
dt

The sensitivity measure is

qs = E {wl AU2(Q) + w2 AY2(tf) + w3 Av2 (Qj)

which may be written as

where
-
0 0

R2
1
0 0 0

0 0 0

Problem Statement Including Sensitivity

The augmented-state vector is

15
with the initial conditions

-
';T(to) = Eta,. . .,q
The problem restatement with sensitivity considerationsis to minimize

subject to

and

Ql = fr2 (4)- L1 = 0 ; L1 = 15 240 m (50000 ft)

q2 = ii.3(tf) - L~ = o ; L~ = o m/sec (0 ft/sec)

Numerical Results
As previously mentioned, the steepest-descent algorithm of reference 17 was imple-
mented to obtain numerical results. These results were computed with the use of a
fourth-order integration subroutine with a fixed step size of 0.5 second. For this example
problem each iteration in the solution required integrationof 15 augmented state differen-
tial equations plus the forward integration of the range equation and backwards integration
of 45 adjoint variable differential equations. Each iteration required about 22 seconds on
the Control Data 6600 computer system used with no special effort having been made to
keep computer run time down. Several check solutions were computed with the use of a
0.125-second fixed step size. These check solutions always agreed to at least five, and
usually to at least seven, significant figures with those computed with the use of the
0.5-second step size.
Numerical results for a variety of cases were obtainedby changing the values of
KD (themodified dragcoefficient), ws (thesensitivityweightingfactor),and w1,w2,w3
(the s e t of relative weighting factors). The different combinations for which results were
obtained a r e shown in table III.

16
 TABLE ID.- CATALOG OF CASES

I Relative-weighting-factor sets
I KD I WS

I Set (1): w1 = w3 = I; w2 = 10-2 I 0 1 0, 0.1, 0.396, 1.0, 10.0

Set (2): w 1 = w3 = 1; w2 = 10-4 I{ :.ool I 0, 0.1, 0.3, 1.0, 10.0

0, 0.1, 1.0, 10.0

The relative-weighting-factor sets were chosen on the basis of the relative impor-
tance attached to velocity (wl applies to horizontal velocity and w2 applies to vertical
velocity) and position (w3 applies to altitude errors) . Set (1) gives equal weight to a
0.305-m/sec (1-ft/sec) velocity error and a 3.05-meter (10-foot) altitude error, while
set (2) gives equal weight to a 0.305-m/sec (l-ft/sec) velocity error and a 30.5-meter
(100-foot) altitude e r r o r .
The modified drag coefficient KD was set at zero, simulating vacuum flight, and
a t 0.001, a representative value for small rockets in the earth's atmosphere. Setting
KD at zero reduced the number of parameters to three. It also made possible, for
ws = 0, an analytical calculus-of -variations solution for the optimal thrust-attitude time
history; namely, the well known linear tangency law (ref. 18).

tan O(t) = (Y + pt

where CY and p a r e constants determined by theboundaryconditions of the problem.

This calculus-of -variations solution was used to validate the steepest -descent-algorithm
programing and solution. The comparison of results showed negligible differences.
Also, for KD = 0, it can be shown that

and
av = -
- 1 V(tf) +tfg
aT T

Therefore, these partials will not change with ws f o r KD = 0.

As an independent check on the sensitivity partials, perturbations in u, v, and y
a t tf due to plus and minus one-sigma errors in each parameter, taken one at a time,
were computed. For example, the one-sigma (lo) perturbation in u due to parameter
errors was computed by using the relation

17
 whereu(9)meansuat tf on a trajectory
with T increased
byits
one-sigma
+loT
value,whereut meansuat tfon a trajectorywith & decreased by itsone-
( f > -I*&
111

sigmavalue,and s o forth.Similarcomputationsweremadefor Av!? (2)

andAylo.One-
sigma perturbation values computed by this method were compared with one-sigma values
computed by using the sensitivity partial derivatives. For example, the one-sigma (lo)
perturbation in u was computed as follows by using the sensitivity partials:

where6Tlo,6mloand s o forth are one-sigma perturbations in the parameters. Simi-

lar computations were made for Avf? andAylo.(1 1
Superscript (1) refers to 10.pertur-
bations obtained by using sensitivity partials, and superscript (2) indicates that the method
of computation was to obtain an average perturbation value by assuming first plus and then
minus la variations in each parameter.

Discussion of Numerical Results

Numerical results for the converged trajectories are summarized in figures 2 to 11,
which give control time histories and trajectory plots, and in tables that give sensitivity
partials,one-sigmaperturbationvalues, and @ s and @p.
Figure 2 shows the control time history used initially and the converged control
time history after 10 iterations for the conditions indicated (KD = 0; ws = 0) to illustrate
how a reasonably chosen control time history is modified by the algorithm to obtain an
"optimal" solution. Optimal here means that a gradient function has been reduced several
o r d e r s of magnitude o r to some reasonable value. A plot of tan O(t) as a function of time
was made for this converged trajectory, and the intercept and slope of a straight line

1%
fairing to that plot are compared in table IV with the values of the constants CY and p
required in the calculus-of -variations (C.O.V.) solution.

TABLE 1V.- LINEAR-TANGENCY-LAW CONSTANTS

..
Constants Faired values I
I I
C.O.V. values
--
CY..... 3.49

"
p ..... -0.0401 ._

The agreement is considered to be good, in view of the crude fairing method and
other factors relating to how the C.0.V. values of a, and p were obtained. Also the
trajectories showed good agreement.
Table V(a) and figures 3, 4, and 5 summarize converged numerical results for
KD = 0 and s e t (1) relative weighting factors for several values of the sensitivity
weighting factors. The sensitivity partials of table V(a) show consistent trends, with
altitude and horizontal-velocity sensitivity partials decrea.sing in magnitude while
vertical-velocity sensitivity partials are increasing in magnitude with increases in ws.
These trends are also reflected in the root-sum-square perturbations; one-sigma (la)
perturbation values of altitude and horizontal velocity decrease while those of vertical
velocity increase. The percentage of the original value (value with ws = 0) of Aula,
Aylo,Avla, C$s, and C$p is plottedinfigure 3 againstsensitivityweightingfactor.
For values of ws larger than about 0 . 3 , little change takes place in Avla, AylU, and
qS, but theperformance C$ continuestodegradealongwith a decrease in Aula.The
( p)
control time histories for these cases, shown in figure 4, exhibit an interesting charac-
teristic. As ws becomes larger, the control tends toward a bang-bangtype of control
where the thrust is either directed straight up (vertical) o r straight down. While i t
appears that the steepest-descent program indicates the existence of a bang-bang optimal
control law, attempts to predict this behavior analytically have been unsuccessful. Alti-
tude is plotted against range in figure5 where it may be observed that the trajectory
becomessteeper as ws increases.
Table V(b) and figures 6 to 9 summarize the results for KD = 0 and s e t (2) rela-
tive weighting factors. Set (2) puts less emphasis on the altitude sensitivity partials than
s e t (1). Thus altitude sensitivity partials increase for set (2) rather than decrease as
they did f o r s e t ( l ) , and both horizontal and vertical-velocity partials decrease. Figure 6
clearly illustrates these results in the plots of Ayla,Avla, and Aula. By comparing
this figure with figure 3, it may be seen that Ayla and Avlo have switched positions
on the plots. Also it appears that $s begins to level off and remain essentially con-
stant at ws = 10.0 in figure 6 whereas it leveled off and became essentially constant at

19
 TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION

(a) KD = 0; w1 = w3 = 1, w2 = 10-2(Set (1))

.

t
~"
~" - - _ .

Sensitivity partials Values for sensitivity weighting factor ws, of -

- " .-
and related terms
(*I 0 0.1 0.396 1.0 10.0
. . -

au/aT 1.106 0.973 0.372 0.163 0.02c

au/alin -6 322 -5 970 -2 044 -838 -101
au/amo -84.98 -77.34 -28.07 -11.93 -1.44
(1) 110.1 101.3 36.1 15.2
Au la 1.8
(2) 110.2 101.4 36.2 15.3 1.8

ay /aT 52.75 52.75 52.75 52.75 52.75

ay/alin -96 140 -82 950 -7 3 640 -72 790 -72 610
ay/amo -2 755 -2 672 -2 613 -2 608 -2 606

AY g 3 209 3 100 3 025 3 019 3 017

AY 3 210 3 100 3 025 3 019 3 017

av/aT 0.805 0.805 0.805 0.805 0.805

%/a& -1 419 -1 646 -1 805 -1 820 -1 823
av/amo -41.74 -43.17 -44.18 -44.27 -44.29
(1) 48.57 50.50 51.91 52.04 52.07
Av la
(2) 48.57 50.51 51.93 52.12 52.18
Av 10.

4 423 3 892 1489 651 79

@P
@S 117 500 108 900 95 500 94 100 93 700
" "
- -
* Superscripts
~

(1)and (2) refer to one-sigma (lo) perturbations obtained by using the

sensitivity-partial-derivative method and the parameter method, respectively.

20
 TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION - Continued

.- ~

Values for sensitivity weighting factor, ws, of -

~

0 0.1 0.3 1.o 10.o

1.106 1.084 0.723 0.448 0.104
-6 322 - 5 869 -2267 -103 174
-84.98 -81.22 -43.76 -18.88 -3.40

110.1 104.4 52.6 22.1 4.5

110.2 104.4 52.6 22.1 4.5

52.75 52.75 52.75 52.75 52.75

-96 100 -102 100 -108 500 -112100 - 126 400
-2755 -2793 -2833 -2856 -2946

3 209 3 260 3 315 3 347 3 474

3 210 3 260 3 316 3 347 3 475

0.805 0.805 0.805 0.805 0.805

-1419 -1316 - 1 206 -1 146 -899
-41.74 -41.09 -40.40 -40.02 -38.46

48.57 47.71 46.82 46.35 44.45

48.57 47.72 46.82 46.33 44.40

4423 4 336 2 891 1791 417

15 510 14 240 6060 3 760 3 200
* Superscripts (1) and (2) refer to one-sigma (la)perturbations obtained by using
:he sensitivity-partial-derivative method and the parameter method, respectively.

21
ws = 1.0 in figure 3. Part of the reason is the change in the magnitude of @s due t o
changing only the relative valuesof the relative weighting factors without regard to their
magnitude. This is clearly illustrated by the near-order-of-magnitude difference in GS
f o r set (1) and set (2) at ws = 0. (Compare tables V(a) and V(b).) The trajectories for
the two sets are the same; @s changes with the change from set (1) t o set (2). Reducing
the values of @s in table V(a)by dividing by the constant '17 '0° = 7.57 so that both
15 510
cases have the same GS at ws = 0, and incorporating this constant into the weighting
factor by multiplying each wsof table V(a) by 7.57, allows a more reasonable compari-
son. This comparison is shown in figure 7 where the ordinate is called ws (adjusted).
This figure shows the differences which come about due to different relative-weighting-
factor sets.
The control time histories for relative-weighting-factor s e t (2) are shown in fig-
u r e 8. Again there is a tendency toward a bang-bang type of control, as ws increases,
which may be noted by observing that the angle difference between the nearly constant
attitude portion of the control time history at the beginning and near the end of flight is
about 180° f o r both ws = 1.0 and ws = 10.0. Figure 9 shows that in comparison with
the set (2) trajectory for ws = 0, the other set (2) trajectories are generally less steep.
The opposite result is shown for set (1) in figure 5; in comparison with the set (1) trajec-
tory for ws = 0, the other set (1)trajectories are more steep. These results indicate
the importance of the relative weighting factors in shaping the trajectories.
Data for a nonzero value of KD,KD = 0.001, and relative-weighting-factor s e t (2)
a r e shown in table V(c) and figures 10 and 11. Data f o r ws = 0.1 a r e shown in table V(c)
but not in figures 10 and 11 because the control time histories and trajectories essen-
tially coincide with those for ws = 0. In table V(c) it may be observed that while consis-
tent data trends are shown in the first three ws columns, the data in the column for
ws = 10.0 are not consistent. The consistent data trends in the first three columns are
very similar to those in table V(b). The inconsistency of the last column may be explained
by examination of figures 10 and 11 where the radically different character of the control
time history and trajectory for ws = 10.0 may be seen. This result for ws = 10.0was
so unusual that it was believed necessary to verify the answer. Accordingly, verification
was obtained by iterating to the same result (essentially) from an additionaltwo different
nominal trajectories. These results therefore appear correct. It is believed that the
data inconsistency results from the different character of the trajectory - that is, the
bending back of the trajectory as seen in figure 11. With the exception of this ws = 10.0
case, the trajectories are l e s s s t e e p with increasing ws just as they were in figure 9
for a similar case with zero drag.
For all of the cases discussed herein, good agreement was observed (see
tables V(a), V(b), and V(c)) between the one-sigma perturbations computed by using
sensitivity partials and those computed by using the parameter-perturbation method.
22
 TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION - Concluded

(c) KD = 0.001; w1 = w3 = 1, w2 = (Set (2))

Sensitivity partials
and related terms
Values for sensitivity weighting factor, ws, of
~
- 1
(*I 0 0.1 1.0 10.0
~

&/aT 25.35 25.24 24.35 12.07

au/a& -756.9 -747.1 -668.1 -002.8
au/amo -8.11 -7.99 -7.04 -1.27
au/aK,, -687000 -6 88 900 -697 500 -359 900
(1) 17.37 17.28 16.52 6.61
(21 17.37 17.28 16.52 6.61
Au la

ay /a T 31.94 31.96 32.11 31.81

ay/al;l -44 640 -44 840 -46 037 -44 800
ay/am, -1 438 - 1 440 -1456 -1 428
ay/aKD 14 190
000 -14 200
000 -14 220
000 -14 350
000

1730 1732 1751 1712

1729 1732 1751 1712

16.18 16.13 15.81 17.01

-7.8 4.6 101.8 34.5
-3.80 -3.68 -2.81 -3.52
-278700 -280 400 -293 500 -313 700

7.358 7.309 7.102 7.766

7.357 7.311 7.100 7.766

1 824 1824 1811 911

@P
@S 655 652 630 403
* Superscripts (1) and (2) refer to one-sigma (lo)perturbations obtained by using
the sensitivity-partial-derivativemethod and the parameter method, respectively.

23

I
 CONCLUDING REMARKS

It has been shown how a sensitivity measure, composed of the weighted sum of the
variances of the performance index and terminal constraints, may be added to the perfor-
mance index of a stochastic optimal control problem to achieve a reduction in performance
and constraint sensitivity due to parameter perturbations.
It was necessary to assume that the parameters of the system remained fixed during
system operation and that the stochastic nature of the problem came about because of
inexact knowledge of these fixed values.
It was shown how this technique increased the dimensionability of the optimization
problem and introduced weighting factors o r constants for use as design parameters.
These weighting factors permit different relative importance to be attached to different
types of sensitivity, such as position relative to velocity, and allow for adjustmentof
performance degradation and of sensitivity reduction.
The feasibility of solving nonlinear multiparameter problems by using this technique
was illustrated by solving an example rocket-trajectory problem through application of a
steepest-descent algorithm. The example solutions also served to illustrate the tradeoffs
made possible by changes in the weighting factors.

Langley Research Center,

National Aeronautics and Space Administration,
Hampton, Va., April 14, 1971.

24
 APPENDIX A

STEEPEST-DESCENT DERIVATION

The problem is to find the control time history u(t), to 5 t 6 tf, which minimizes

subject to the differential constraints

-
2 = f(z,z,t) ; x(to) = zo (given)

and the terminal constraints

Thesolution is obtainediteratively.Choose a reasonable time history u*(t)

- and
obtain x* (t), a solution to equation (A2). This solution, in general, is such that neither
@ is minimumnor "
Q = L.
Linearizing about this solution gives

or

Now, for convenience, let A(t)bethesolution of

25

Il I I 1 1l1l 1 I
 APPENDIX A - Continued
with

where
N

Q
- = ];[
Then equation (A5) becomes

82
- =-
or, letting b+ - ,
ax 6x(tf>

Partition A (t) as-follows :

Then, since -
x (to) - (to) = 0,
is given and 6x

and

-
Now for some measure of allowable control perturbation

(dP)2 = ltf
t0
- dt ; W = W T > 0
6gTW 6u (A13

minimize 64 and choose - such that

6+ +( x(tf)) = L
" will be satisfied (or more nearly
so). Form an augmented function to be minimized

26
 APPENDIX A - Concluded
where p and - v a r e constantLagrangemultipliers. For JA tobeanextremum,its
variation with respect to - must be zero - that is,
6u

This requirement implies that

for all t. Solving for 6u,

-

In o r d e r t o solvefor v thisexpressionfor
p and - - is substitutedintoequa-
6u
tions (A10) and (All) to obtain

where the sign on the radical is minus because @ is being minimized and where
\

Now let

- = u*(t)
u(t) - + 6u(t)
-

be the new control time history used to obtain a solution to equation (A2) and repeat the
same procedure until - ) - L and the gradient of the payoff-constraint surface
I
-

a r e sufficiently close to zero.

27
 REFERENCES

1. Dorato, P.: On Sensitivity in Optimal Control Systems. IEEE Trans. Automat. Contr.,
vol. AC-8, no. 3, July 1963, pp. 256-257.
2. Pagurek, Bernard: Sensitivity of the Performance of Optimal Control Systems to Plant
Parameter Variations. IEEE Trans. Automat. Contr., vol. AC-10, no. 2, Apr. 1965,
pp. 178-180.
3. Witsenhausen, H. S.: On the Sensitivity of Optimal Control Systems. IEEE Trans.
Automat. Contr., vol. AC-10, no. 4, Oct. 1965, pp. 495-496.
4. Sinha, N. K.; and Atluri, Satya Ratnam: Sensitivity of Optimal Control Systems.
Proceedings Fourth Annual Allerton Conference on Circuit and System Theory,
William R. Perkins and J. P. Cruz, Jr., eds., Univ. of Elinois and IEEE, 1966,
pp. 508 -5 16.
5. Breakwell, John V.; Speyer, Jason L.; and Bryson, Arthur E.: Optimization and
Control of Nonlinear Systems Using the Second Variation. J. SOC. Ind. Appl. Math.
Contr., ser. A, vol. 1, no. 2, 1963, pp. 193-223.
6. Kreindler, 'Eliezer: Synthesis of Flight Control Systems Subject to Vehicle Parameter
Variations.AFFDL-TR-66-209, U.S. Air Force, Apr. 1967.(Availablefrom DDC
as AD 653 600.)
7. Higginbotham, Ronald E.: Design of Sensitivity-Constrained Optimal Linear State-
Regulator Control Systems. Conference Record of Second Asilomar Conference on
Circuits & Systems, 68 C 64-ASIL, IEEE, 1968, pp. 330-335.
8. H a a s , V. B.; and Steinberg, A.M.: Minimum Sensitivity Optimal Control f o r Nonlinear
Systems. TR-EE67-10 (NGR 15-005-021), School Elec. Eng., Purdue Univ., Aug.
1967.(Available as NASA CR-88863.)
9. Rohrer, R. A.; and Sobral, M., Jr.: Sensitivity Considerations in Optimal System
Design. IEEE Trans. Automat. Contr. vol. AC-10, no. 1, Jan. 1965, pp. 43-48.
10. Chan, S. Y.; and Chuang, K.: A Study of Effects on Optimal Control Systems Due to
Variations in Plant Parameters. Int. J. Contr., First ser., vol. 7, no. 3, Mar. 1968,
pp. 281-298.
11. Patrick, L. Benjamin; and D'Angelo, Henry: Optimal Control of Plants With Random
Slowly-Varying Parameters. Proceedings Fifth Annual Allerton Conference on
Circuit and System Theory, J . B. Cruz, Jr., and T. N. Trick, eds., Univ. of Illinois
and IEEE, 1967, pp. 711-720.

28
12. Holtzman, J. M; and Horing, S.: The Sensitivity of Terminal Conditions of Optimal
Control Systems to Parameter Variations. IEEE Trans. Automat. Contr.,
vol. AC-10, no. 4, Oct. 1965, pp. 420-426.
13. Kokotovic, P. V.; and Rutman, R. S.: Sensitivity of Automatic Control Systems
(Survey). Automat. Remote Contr., vol. 26, no. 4, Apr. 1965, pp. 727-749.
14. Rosenbaum, Richard C.; and Willwerth, Robert E.: Launch Vehicle Error Sensitivity
Study. NASA CR-1512,1970.
15. Lindgren, B. W.; and McElrath, G.W.: Introduction to Probability and Statistics.
Second ed., Macmillan Co., 1966.
16. Coddington, E a r l A.; and Levinson, Norman: Theory of Ordinary Differential Equa-
tions.McGraw-Hill BookCO., Inc.,1955.
17. Bryson, A. E.; and Denham, W. F.: A Steepest-Ascent Method for Solving Optimum
Programming Problems. Trans. ASME, s e r . E: J. Appl.Mech., vol. 29, no. 2,
June 1962, pp. 247-257.

18- BrYson, Arthur E., Jr.; and Ho, Yu-Chi:Applied Optimal Control. Blaisdell Pub. CO.,
c.1969.

29
 T

Figure 1.- Coordinate system and geometric relations.

30
I

\
II ._L- .. . -1- I I
0 20 40 60 80 100
Time from lift-off, t, sec

Figure 2.- A comparison of the initial guessed thrust-attitude

control time history with the steepest-descent generated
optimal thrust-attitude control time history. KD = 0; ws = 0.

31
w
N
I20

IO0

80

60

40

20I c
I

0-
.o I .I I IO
WS
Figure 3.- One-sigma (lo) e r r o r components and sensitivity and performance indices
plotted against sensitivity weighting factor. KD = 0; w1 = w3 = 1, w2 = 10-2.
 \

-c
F -40

-7 0
-80
-9oL
0
1 __I"L_ L
1- 1~ .I"1 - 1
I O 20 30 40 50 60 70 80 90 100
Time from lift-off, t, sec
Figure 4.- Comparison of optimal control time histories for different sensitivity
weighting factors KD = 0; w1 = w3 = 1, w2 =

33
 I< lo3
ws = 10.0

1
50 75 25
i
I50x103
Range, x, ft
J I 1 I I 1
0 IO 20 30 40 50
Range, x , km
Figure 5.- Comparison of optimal trajectories for different sensitivity
weighting factors. KD = 0; w1 = w3 = 1, w2 = 10- 2 .

34
 120-

0
“m 80
3
4-
0
0)
-
3

> 60
Y-
O
t
c
a
40
a

20

Figure 6. - One-sigma (lo) e r r o r components and sensitivity and performance indices

plotted against sensitivity weighting factor. KD = 0; w1 = w3 = 1, w2 =
 \ \ \

\\ \
tP
- - - - - Set( I): WI = w 3 = I, w2=
0
> \ \ *
\
\
Set (2):WI = w 3 = 1 , w = IO"'

I I , I I I I I 1 I I I I 1 1 1 1 I I I 1 I I l l I
-1"-~ II
01
.01 .I I 10 100
ws (adjusted)

Figure 7.- Set (1)and set (2) sensitivity and performance indices plotted against
a sensitivity weighting factor adjusted for equal sensitivity at w2 =O. KD=O.
 \
-

O
+L
-40-
2 -50-
v,

e -60-

I
/

-I 10
-120
- I 30
-140 t
I l l I l l l l l I
102030405060708090 100
Time from Iift-off,t, sec
Figure 8.- Comparison of optimal control time histories for different
sensitivity weighting factors. KD = 0; w1 = w3 = 1, w2 =

37
 40 -

i 1
IO0 I25 l50x IO3
Range, x, f t
I I I"
0 '0 20
R a n g e , x, k m
30
~

40
.J
50

38
 I 40r

- 1801
-225I \ws = 10.0
" l I I I 1 l I I I ~
0 20 40 60 80 100
Time from lift-off, t, sec

Figure 10.- Comparison of optimal control time histories for different

sensitivity weighting factors. KD = 0.001; w1 = w3 = 1, w2 = 10-4 .
i
0

16
50X lo3
rws = ‘Oa0

40 - / I L ws = o

c
w- 30-
r:
c
Q)
U
3
E 20-
-
a
4

0
-20 -10 0 IO 20 30 40 50 60x IO3
Range, x, ft
I I I I 1 I I
-4 0 4 8 12 16 20
Range, x, km
NATIONAL
AERONAUTICS
AND SPACE ADMINISTRAI
ION
WASHINGTON,D. C. 20546
OFFICIAL BUSINESS FIRST CLASS MAIL
PENALTY FOR PRIVATE USE $300
POSTAGE AND FEES PAID
NATIONAL AERONAUTICS ANI:
SPACE ADMINISTRATION

1 1 U 001 44 51 3DS 71166 00903

A I R FORCE WEAPONS LABORATORY /WLOL/
K I R T L A N O A F B ? NEW M E X I C O 87117

AT1 E. LOU BOWMAN? CHIEFtTECH- L.IBRARY

If Undeliverable (Section 158

Postal Manual) Do Not Return

“ T h e aeronaz~tical and space activities of the United States shall be

conducted so as to contribute , . . t o the expansion of h u m a n knowl-
edge of phenomena in the at?tlosphere and space. T h e Adnzinistration
shall provide for the widest prncticable and appropriate dissemination
of information concerning its actizdies and the resdtz thereef.”
-NATIONALAERONAUTICS
AND SPACE ACT OF 1958

NASA SCIENTIFIC AND TECHNICAL’ PUBLICATIONS

TECHNICAL REPORTS: Scientific and TECHNICAL TRANSLATIONS: Information
technical information considered important, published in a foreign language considered
complete, and a lasting contribution to existing to merit NASA distribution in English.
knawledge.
SPECIAL PUBLICATIONS: Irdormation
TECHNICAL NOTES: Information less broad derived from or of value to NASA activities.
in scope but nevertheless of importance as a Publications include conference proceedings, .*,
contribution to existing knowledge. monographs, data compilations, handbooks,
sourcebooks, and special bibliographies.
TECHNICAL MEMORAnjDUMS:
Information receiving limited distribution TECHNOLOGY UTILIZATION
because of preliminary data, security classifica- PUBLICATIONS: Information on technology
tion, or other reasons. . used by NASA that may be of particular
interest in commercial and &her non:aerospace
CONTRACTOR REPORTS: Scientific and applications. Publications include Tech Briefs,
technical information generated under a NASA Technology Utilization Reports and
contract or grant and considered an important
Technology Surveys. ’

contribution to existing knowledge.

Details on the availability of these publications may beobtained from:

SCIENTIFIC AND TECHNICALINFORMATIONOFFICE

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION

Washington, D.C. PO546