Citation
Computational techniques and performance criteria for the optimum control of nuclear system dynamics

Material Information

Title:
Computational techniques and performance criteria for the optimum control of nuclear system dynamics
Added title page title:
Optimum control of nuclear system dynamics
Creator:
Saluja, Jagdish Kumar, 1934-
Place of Publication:
Gainesvill
Publisher:
University of Florida
Publication Date:
Copyright Date:
1966
Language:
English
Physical Description:
xvi, 148 leaves : ill. ; 28 cm.

Subjects

Subjects / Keywords:
Analog computers ( jstor )
Control loops ( jstor )
Control theory ( jstor )
Cost functions ( jstor )
Neutrons ( jstor )
Nuclear reactors ( jstor )
Nuclear rocket engines ( jstor )
Optimal control ( jstor )
Reactivity ( jstor )
Simulations ( jstor )
Controlled fusion ( lcsh )
Dissertations, Academic -- Nuclear Engineering Sciences -- UF ( lcsh )
Nuclear Engineering Sciences thesis Ph. D ( lcsh )
Nuclear reactors ( lcsh )
Nuclear rockets ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Bibliography: leaves 145-147.
Additional Physical Form:
Also available on World Wide Web
General Note:
Manuscript copy.
General Note:
Thesis - University of Florida.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
021507752 ( AlephBibNum )
13142428 ( OCLC )
ACW7084 ( NOTIS )
1143692801 ( OCLC )

Downloads

This item has the following downloads:


Full Text












COMPUTATIONAL TECHNIQUES AND
PERFORMANCE CRITERIA

FOR THE OPTIMUM CONTROL OF
NUCLEAR SYSTEM DYNAMICS






By
JAGDISH K. SALUJA


A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA
December, 1966


























DEDICATED TO


'om, Dad, My Brothers and

Sisters, and all my friends

for their constant ~ncourage-

ment.











ACKNOWLEDGMENTS


The author wishes to acknowledge the help of his

chairman, Professor A.P. Sage, who has always been avail-

able for consultation and who has been a source of continued

encouragement and enlightening advice. The author also

wishes to express his sincere appreciation to Dr. R.E.

"hrig for his encouragement during the course of his

pursuance of the doctoral degree and for being on his

supervisory committee. His help is greatly appreciated.

The concern and interest of the remaining members of the

committee, Professors M.J. Ohanian, R.B. Perez, A.E.S.

Green and P.G. Selfridge had a most beneficial effect

upon the works progress and is sincerely appreciated.

The author also wishes to thank Dr. P.P. Eisenberg

for his help with the digital computer computation, using

the quasilinearization technique and to Dr. T.W. Ellis,

Major J.T. Humphries, Mr. H. Atwater and Mr. Kesavan

Nair for many helpful discussions. The cooperation and

the assistance, especially of Mr. Kenneth Fawcett on the

Analog computer is greatly appreciated.

Furthermore, the author wishes to thank the depart-

ment of Nuclear Engineering Science for the financial

support, without which this work would not have been pos-

sible. It is also a pleasure to extend thanks to Mrs.











Mazelle Shadburn for her assistance in typing the disserta-

tion and to Mr. George Perry and Mr. Gerald Murdock for

the drawings.













TABLE OF CONTENTS


Page


DEDICATICN . . . . . . . . . ii

ACKNCWLEDGMENTS .................. iii

LIST OF TABLES . . . . . . . ... vii

LIST CF FIGURES . . . . . . . . viii

LIST OF SYMBOLS . . . . . . . . xii

ABSTRACT . . . . . . . . . . xiv

CHAPTER


INTRODUCTION . .

OPTIMAL CONTROL CF
SYSTEM ....


A NUCLEAR REACT. .
A NUCLEAR REACTCR


A. Introduction and Model Description.

1 Statement of the Maximum Principle.

2 Statement of the Specific Optimal
Control Problem (s.o.c) .. . ..

3 Gradient Technique . . .....

4 Quasilinearization Technique .....

B. Formulation of Specific Problems. .

1 Cost Function 1 . . . . .

Problem 1 via quasilinearization
technique . . . . . . .

2 Cost Function 2 ..........
Simulation of problem 2 .. ....

3 Cost Function 3 . . . ...
Simulation of problem 3 . . .











TABLE OF CONTENTS (Continued)


4 Cost Function 4. . . . .
Simulation of problem 4 . ...

C. Discussion of Results. . . .

III. CPTIVAL CCNTRCL CF A NUCLEAR ROCKET
ENGINE . . . . . . .

Introduction and Model Description

Solution Procedure . .....

Variable start-up time. . . .
Variable desired temperature. .
Influence of Parameter variation
on system response . . .

IV. SUMMARY AND CONCLUSIONS . . . .

Summary and results. . .

Suggestions for further study... .

APPENDICES . . . . . . . . .

Appendix A: Analog Simulation . . . .
Circuit Symbolism . ....
Problem Scaling. . .. ....
Simulation by Gradient Technique

Appendix B: Computer Codes . . . .
Solution of TPBVP via Quasi-
linearization . . . . .
Mimic Code . . . . . .

LIST (F PEFFRFNCES . . . . . . .

BIOGRAPHICAL SKETCH . . . . . . . .


Page

58
59

67


77

77

85

85
86

94

99

99

102

103

104
104
110
126

138

138
140

145

148













LIST OF TABLES

Table Page

1. Nuclear Reactor Constants for a Lumped System
and Initial Conditions . . . ... 20

2. Initial Conditions on the Adjoint Variables
And the Cost J, To Achieve Desired Power
Levels . . . . . . . . 25

3. Typical Results of Quasilinearization Tech-
nique Applied to Problem 1 (Nuclear Reactor
System). . . . . . . . . 40

4. Cost Functions as a Function of Weighting a
and Power level Nd for Problem 2 . .. 51

5. Cost Function as a Function of Power Level
F and Constants A and S for Problem 4
(Nuclear Reactor System) . . ... 61

6. Rocket Engine Constants. . . . . . 81

7. Scale Factors for the Four Control Problems. 113

8. Potentiometer Settings for Control Problems
1 and 2 (Nuclear Reactor System) . 115

9. Potentiometer Settings for Control Problem 2
(Nuclear Reactor System) . . . 117

10. Potentiometer Settings for Closed Loop Control
Problem 3 (Nuclear Reactor System) . .. 119

11. Potentiometer Settings for Control Problem 4
(Nuclear Reactor System) . . ... 123

12. P 's as a Function of Iteration Step for
Jo
Problem 1 (Nuclear Reactor System) . . 128












LIST OF FIGURES


FIGURE PAGE

1. Power level N and reactivity Pc versus time
for optimal control problem 1. (Nuclear
Reactor System). . . .. . .. 23

2. Cost function versus desired power level . 24

3. Power level and associated reactivity time
plots, leading to the desired trajectory
using gradient technique. (Nuclear Reactor
System). . . . . . . . 30

4. Number of iterations for convergence versus
P O. .. . . . . . . . . 31
5. Maximum value of N and P20 versus iteration
number in an iteration cycle . . . . 33

6. Power level versus time (A= 10-4) Comparison
of analog computer results with those ob-
tained by quasilinearization (Nuclear Reactor
System). .............. . 42

7. Power level versus time ( = 10-3) Comparison
of analog computer results with those ob-
tained by quasilinearization (Nuclear Re-
actor System) . . . . . . .. 43

8. Power level versus time (A= 10-3) Comparison
of analog computer results by gradient tech-
nique, with those obtained by quasilineari-
zation (Nuclear Reactor System) . . . 44

9. Total cost function versus peak power desired
(control problem 2) . . . . . 48

10. Power level N and reactivity Pc versus time for
problem 2 for various values of a and Nd
10 no (Nuclear Reactor System) . 49

11. Power level N and reactivity P versus time
for control problem 2 for various values of
d and Nd = 5no (Nuclear Reactor System) 50









LIST OF FIGURES (Continued)

FIGURE PAGE

12. Cost function associated with control rod
movement versus a Nd =10 no (control
problem 2) . . . . . .. 52

13. Total cost function versus a for control
problem 2 and Nd = 10no. ...... . 53

14. Variation of parameters in the control law
for s.o.c. problem (i.e., problem 3) and
for various values of Nd . . ... 56

15. Peak power level reached versus A (control
problem 3) . . . . . 57

16. Power level N and reactivity pc versus time
for specific optimal control problem 4
(Nuclear Reactor System) . . . . 62

17. Cost function versus peak power attained 63

18. Nd versus A and B. ........ . 64

19. Auxiliary functions, P2, P3, P4 and P5 ver-
sus tine for control problem 4; Nd = 10no
and 6no (Nuclear Reactor System) . 65

20. Auxiliary functions P2, P3, P4 and F5 ver-
sus time for control problem 4; Nd 5no
and 2no (Nuclear Reactor System) . 66

21. Cost function versus peak power attained
for open and closed loop optimal control
(Nuclear Reactor System) . . 69

22. Total cost function versus peak power attain-
ed for open loop control and cost varia-
tion for fixed and optimally varying
values of control parameters for closed
loop control (Nuclear Reactor System). 70

23. Open loop and closed loop optimal trajec-
tories N(t) versus time (Nuclear
Reactor System) . . ........ . 72











LIST CF FIGURES (Continued)

FIGURE PAGE

24. Open loop and closed loop optimal control
function to obtain ten times the initial
power level in a Nuclear Rocket System
73

25. Cost function versus peak power desired. . 74

26. Total cost function versus mean neutron life-
time . . . . . . .. . 75

27. Schematic Diagram of a nuclear rocket engine 79

28. Control poison reactivity and temperature
time plots, tf = 4 secs, 6 secs, and 11
secs (Nuclear Rocket Engine System). . 87

29. Pressure-time plot, tf = 4 secs, 6 secs and
11 sees (Nuclear Rocket Engine System) . 88

30. Optimal control poison reactivity and temper-
ature-time plots, Td = 0.4 (Nuclear Rocket
Engine System) .. . . ..... 89

31. Optimal pressure-time plot, 'd = 0.4 (Nuclear
Rocket Engine System) ... . 90

32. Optimal pressure and neutron density versus
tire. Td = 0.7 and 1.0 (Nuclear 'ocket
Engine System) . . . . . 91

33. Optimal pressure-time plot, Td = 0.7 and 1.0
(Nuclear Rocket Engine System) . . .. 92

34. Cost function, J versus desired average core
temperature (Nuclear Rocket Engine System) 93

35. Control poison reactivity versus time and
average core temperature versus time for
three cases of parameter variations (Nuclear
Rocket Engine System) . . . .. 96

36. Pressure versus tine response fcr three cases
of parameter variations, (Nuclear Rocket
Engine System) . . . . . . . 97











LIST OF FIGURES (Continued)

FIGTRE PAGE

37. Pressure versus time and average core temper-
ature versus time for three cases of para-
meter variations (Nuclear Rocket Engine
System) . . . . . .. . . ... 98

38. Symbols used in analog computer diagrams. . 105

39. Logic symbols used in computer diagrams . 107

40. Special circuits. . . . . . . ... 109

41. Analog computer diagram for control problems
1 and 2, for a Nuclear Reactor System. . 114

42. Analog computer diagram for control problems
3 and 4 (Nuclear Reactor Systen. . . .. 118

43. Computer control diagram. . . . ... 134

44. Updating circuit. . . . . . . ... 136

45. Box-oriented mimic diagram. ....... .. 141

46. Mimic instructions. . . . . . ... 141













LIST OF SYMBOLS


n Neutron density.

nd Desired neutron density.

N Scaled neutron density.

Nd Scaled desired neutron density.
No Normalized neutron density.

c Delayed neutron precursor density.

C Scaled delayed neutron precursor density.

C' Normalized delayed neutron precursor density.

\ Nean effective neutron lifetime (sec).

P Core inlet pressure (psi).

Pd Design core inlet pressure (psi).
P' Normalized core inlet pressure.

T Maximum core surface temperature (OR).

Td Design maximum core surface temperature (oR).

T' Normalized maximum core surface temperature.

V Valve setting.
Vg Normalized valve setting.

a Error weighting term,

at Temperature coefficient of reactivity (OR-l).

al Normalized temperature coefficient of reactivity
(dollars).

ah Propellant density coefficient of reactivity (R/psi).









LIST OF SYMBOLS (Continued)


a. Normalized propellant density coefficient of re-

activity (dollars).

B Delayed neutron fraction.

A Decay constant (sec ).

P Total reactivity.

P' Normalized total reactivity (dollars).

Pc.r Control poison reactivity.

Pc.r Normalized control poison reactivity (dollars)
T Time constant (sec).

-p Pressure time constant (sec).

A,B Ccnstants in the control law.

A.E Scaled constants in the control law.

H Hamiltonian.

J Cost function.

Pi i-th component of the adjoint vector P.
tf Time for solution, final tire.

t Time.

to Initial tire.

ui Ccnponent of the control vector u.

xi Corpcnent of the state vector x.


xiii










Abstract of Dissertation Presented to the Graduate Council
In Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy

COMPUTATIONAL TECHNIQUES AND PERFORMANCE CRITERIA
FOR THE OPTIMUM CONTROL OF NUCLEAR SYSTEM DYNAMICS

by

Jagdish K. Saluja

December, 1966

Chairman: Dr. Andrew P. Sage
Major Department: Nuclear Engineering Sciences

The main objective of this study is to investigate

optimum control of the dynamics of nuclear systems and

associated computational problems. Two systems are con-

sidered: 1) Nuclear Reactor System, 2) Nuclear Rocket

Engine System. For the first system, three problems are

analyzed: a) the neutron density or reactor power is trans-

ferred from one steady state condition to another with a

quadratic constraint on control rod movement. By varying

the operating time and placing an inequality constraint on

reactivity a form of minimum time control is obtained, b)
integral of quadratic form of error in desired neutron

density is minimized with quadratic integral constraint

on control rod movement, c) a fixed configuration optimum

closed-loop control system is realized by letting control

rod rcvement be a proportional plus integral function

of neutron density error.









The two-point boundary value problems (TPBVP), re-

sulting by the application of optimization theory to point

reactor models are obtained. A modified gradient technique

is used to exploit the hybrid capabilities of modern analog

computers. The results are compared with the solutions

of this problem by the numerical technique of quasilinear-

ization. In all of these, the convergence is obtained in

no more than four iterations.

In the examples considered here, the cost function

is higher for the closed loop case as opposed to the open

loop case. This is in general agreement with physical
insight and previous results, since an assumed form

of closed loop controller is used. The coefficients A

and B in the control law are highly dependent upon the final

time tf, the system parameters, and the initial and desired

final state. Variation in system parameters may, however,

render the closed loop control better than the open loop

control which does not change with changing system parameters.
For the case of the nuclear rocket engine system,

only a quadratic constraint on total reactivity is con-

sidered. The nuclear rocket engine is similar to the solid-

core 25,000 megawatt, million pound thrust model developed

by Smith and Stenning. Here, the results obtained for

various start-up times and various desired power levels in

terms of control needed and required cost, appear very










reasonable. Also it was observed that if the pressure and

the temperature time constants were increased from their

normal values, the temperature response became somewhat

sluggish. Possible improvements to this situation are

suggested.











CHAPTER I


INTRODUCTION


Just as all other fields of scientific thought have

undergone a major change in recent years, so also has the

field of rocketery. From the Chinese rockets of 1232

A.D. to those of the Germans during the second world war,

man has made a "case" for nuclear rocket propulsion as a

vital factor in the success of future space travel. This

has been pointed out in the recent literature (1,2,3,4)

and accepted in the field today, Its acceptance is fur-

ther strengthened by the successful testing of the KIWI

and NERVA reactors.

This sudden interest in space travel has led to
thoughts about improved system performances, as a result

of which a very instructive book on nuclear rocket system

optimization by Bussard (5) has appeared. An outgrowth

of this interest is the optimal control problem wherein

the control is such as to give the "best performance.

From the nuclear literature, however, it is clear that

the subject of optimal control of nuclear systems has

received little interest. Recently a number of papers

have appeared in this area, ranging from optimal reactor

shutdown programs for control of xenon poisoning (6,7),










and synthesis of optimal nuclear reactor systems (8,9) to
optimal control of nuclear reactor systems (10,11).

More recently Mohler and Weaver have reported some work

in the application of modern optimal control theory to

the dynamics of nuclear rocket systems (11,12) but much

still remains to be done.

The mathematical theory used for the study of

optimization problems is the calculus of variations. It

is of interest to note that Goddard (13) recognized this

technique as an important tool in the performance analysis

of rockets in his early paper published 40 years ago.

There are many possible variational methods fcr minimizing

or maximizing a functional over a function space. Among

the commonly used methods in control system design are:

(1) calculus of variations

(2) maximum principle
(3) dynamic programming

In all cases the object is to find the optimum con-

trol law such that the given functional of the performance
indices is minimized or maximized. These methods are re-

lated (5) in that the maximum principle can be derived from

the calculus of variations as well as from the method

of dynamic programming. Pontryagin's maximum










principle is a very elegant method, which is applicable

to a broader class of problems than the calculus of

variations. An outcome of this principle is the so called

two-point boundary value problem, abbreviated as "TPBVP,

with split boundary conditions.

This method will be applied in the study of dynamics

of nuclear systems. Since the dynamic response of nuclear

systems is represented by a set of nonlinear, time-varying

differential equations, conventional analysis is usually

too difficult to give a closed form solution. Develop-

ment of efficient computer techniques is therefore necessary.

A numerical approach which is currently quite

popular is the method of gradients (5). It proceeds by

solving a sequence of sub-optimal problems with the prop-

erty that each successive set of solution functions yields

an improved value for the functional being optimized.

Another numerical technique of considerable interest

is the approximating process of quasilinearization which

produces a sequence of approximations xn, n = 1,2,... that

converge quadratically (14) to the solution of the ori-

ginal two-point boundary value problem. This technique

was initially suggested by Bellman and Kalaba (15) as a

modification of the Newton-Raphson-Kantorovich method.

This converts the problem of solving a two-point boundary

non-linear differential equation into a problem of solving

a sequence of two-point boundary linear differential equation










problems. The solution of the latter equation is con-
siderably easier to find.
For the study of optimization problems, a precise
formulation of the performance index is necessary. The
performance index is assumed to be given by:

J(xotf,u) = t 0 (xU,t)dt
0

where 0 is a scalar, a function of the state vector x, the
control vector u and time.
Corresponding to this performance index, a control law
is found which is the best in the sense of maximizing
the performance index.
The performance indices considered in this study are:

i) J = 1/2 Pc2 dt

a) n(tf) = nd

b) T(tf) = Td
This performance criterion may be interpreted as an optimi-
zation of control rod movement, limiting erratic flux changes

(10). Clearly it tends to limit the maximum reactor period.

ii) 3 = 1/2 tf [ pc2 + a (nnd)2 ] dt
io +

Integral of the error squared of the reactor power and its
desired value is minimized along with the integral of the
square of the reactivity function. a is the error weight-
ing term











iii) J = 1/2 jf [ Pc2 + a(n-nd)2 ] dt


where PC = a + A(n-nd)

a and A are constants to be determined.

iv) J = 1/2 tf [ Pc2 + a(n-nd)2] dt
o

where Pc = A(n-nd) + Bf (n-nd) dt
o

A and B are constants to be determined.


The last two criteria give closed loop control and since

a specific form is chosen for the controller, it will

henceforth be referred to as the specific optimal control

problem.

The above four controls are analyzed and compared

with one another for two reactor models. The two models

analyzed are:

(1) Point model of a nuclear reactor system.

(2) Point model of a nuclear rocket engine.

Chapter II contains the description of the first model.
Chapter III concerns the second model and follows along





6





the same pattern as Chapter II. Chapter IV includes

a discussion of results and suggestions for further

study.












CHAPTER II


OPTIMAL CONTROL CF A NUCLEAR REACTOR SYSTEM


A. Introduction and Model Description


Need frequently arises in nuclear reactors as well
as in nuclear power plants for power level changes. These

changes are generally accomplished by semi-automatic

means in small research reactors. In large power re-

actors and in nuclear rocket engines, however, automatic

control is desired. If this control is optimal in some

specified sense, economic or other advantages could be

gained.

The reactor model considered here is that of a one

delayed-neutron group without temperature feedback and

without spatial and angular dependence, described by the

following system of differential equations:


S = (p B )n + ).c (2.1)


S= n Xc (2.2)

where n and c, the neutron density and precursor popula-

tion are the state variables and the reactivity p is








the control variable. The dot represents derivative
with respect to time.

A is the average neutron lifetime
is the decay constant, and
j is the delayed neutron function.

The system is at steady state for t < o and the
initial conditions are: n(o) = no and c(o) = co.
S It is desired to optimize system performance such as to
minimize the following integrals.

tf
1) J = 1/2 f Pc2(t) dt

with n(tf) = nd and I Pc < Pmax


2) J = 1/2 f [Pc2 + a (n-nd)2] dt

in fixed time tf


3) J = 1/2 jf [ Pc2 + a(n-nd)2 ]dt

where pc = a + A(n-nd)
in fixed time tf


4) J = 1/2 [ PC2 + a (n-nd)2 ] dt


tf
where pc = A(n-nd) + B /f (n-nd) dt

in fixed time tf.









1. Statement of the P'aximum Principle

Consider the control processes to be described by
a vector differential equation:

x = T(x,ut) x(o) = xo (2.3)

where T(x,u,t) is a vector with coordinates fl(x,u,t),

f2(x,u,t) ....... fn(x,u,t).
The control functions have to be piecewise continuous
and can be restricted by the inequalities:

gj(u .... ur) 5 0 (j = 1,...m) (2.4)



Assume a performance criterion of the type:

tf
J = 0[x(t)(),u(tf). +f+ f (x,u,t) dt (2.5)
to

where "J" is generally referred to as the performance
criterion or performance index or cost function.
The Hamiltonian of the system is now defined by.

H(x,u,p.t) = (x,5,t) + pT(t) f(x,u.t) (2.6)


This Hamiltonian becomes the Hamiltonian of classical
mechanics only if the control u becomes optimal control
u* defined by equation (2.9).










The variable appearing in Equation (2.6) is called
the auxiliary vector or the adjoint vector or momentum
function. The system equations can now be written in the
following manner.


xi = H( ,; ,p,t) (2.7)
aPi

Auxiliary or adjoint, or costate functionsPi are defined by


p (t) = (2.8)
axi

The Hamiltonian is now minimized with respect to the control
variable u, such that


HO(x,p,t) = Min H(x,p,u,t) (2.9)

where uO(x,p,t) minimizes H and the symboln- refers to

the admissible control space. If ~-Lis infinite (2.9) be-
comes -, = 0. The boundary conditions are given by:

t(to) = Xo (2.10)

and P(tf) = 7x (x,u,t) ,where ( Vx 8 )= 6x
t = tf
2. Specific Optimal Control

The desirability of closed loop control laws has led
to the development of "specific optimal control."










The specific optimal control (S.O.C.) problem is
defined in the following manner (17): given a plant with
a state equation of the form

x = T(;,G,t) (to) = xo (2.11)

It is desired to determine the unknown parameters in a

control law of the form

= h(,) = 0 (2.12)

In the above equation y is a p-dimensional vector which is
a known function of the state and E is a q-dimensional
constant vector to be determined such that an index of
performance of the form


J(u) = g(x,u) dt (2.13)
Jo

is minimized.
The specific optimal control problem can be solved
by a simple re-formulation of the above problem. The re-
formulation is carried out by substituting equation (2.12)
into equations (2.11) and (2.13), leading to

x = f(x,h(y,b),t)
tf
J(;) = f g(7,h(y,E)) dt
Jo











Since y is a known function of x and b is a constant vector

x = f(x.F.t)

J(b) = f gR(x.b)dt (2.14)
fo

and b = o

Thus the S.O.C. problem has been embedded in the pro-
blem of finding the vector E., which minimizes the cost
function and differential constraints imposed by equation
(2.14). The above is now in the form of the typical opti-
mal control problem. Using Pontryagin's maximum principle
then allows one to solve the S.O.C. problem by solving a
two-point boundary-value problem.

3. Gradient Technique

This technique is a successive approximation pro-
cedure in which an iteration is made on the control vari-
able u(t) so as to improve the function to be minimized
at each new iterate, subject to the constraining differ-
ential equations.
Consider, for example, those problems in which it
is required to extremalize performance indices of the
form










tf
J(x,ut)= o 0 (x,u,t) dt (2.15)


subject to the constraint in the form of the first order
differential equations

ei = xi Ji(C,u,t) = o (2.16)

It is required to find u*(t) which at the same time ex-
tremalizes J and satisfies equation (2.16). The overall
algorithm which simultaneously satisfies both conditions
gives (i+l)st iteration of u in terms of the ith iteration:

ui+l ui + K [J m ]+ A (2.17)
= I + j u


where m denotes the number of equations 0 and K and A are
chosen arbitrarily. It was seen in the section
that the maximum principle applied to the above set of
equations, resulted in a coupling equation, coupling the
state vector and the adjoint vector via the control vari-
able u. This result therefore can modify equation (2.17)
to yield the following algorithm which approximates some of
the characteristics of (2.17).

pi+l = i + K [ P (tf) Pd(tf) (2.18)


This algorithm was used in the problem solved in section
B.1. This is an excellent approximation to (2.17) which
is most useful on a repetitive analog computer.










4. Quasilinearization Technique

This technique,developed by Kalaba, Bellman and others

is particularly suited for the type of optimization problems

being considered here. It provides an effective computa-

tional tool for a wide class of non-linear two-point and

multi-point boundary value problems (15,18.19). The

technique consists of producing a sequence of approxi-

mating functions xO(t), x(1)(t).... x(n)(t) which converge

rapidly to the solution of the original problem. In (14.20)

it is shown that this convergence is quadratic, i.e.,

asymptotically the number of correct digits in each approx-

imation is doubled.

To illustrate the technique, consider a 2N-dimensional
dynamical vector equation:

x = (;) (2.19)

subject to the boundary conditions
xi(o) = bi i = 1,2,...,N
(2.20)
xi(tf) = bi i = N+1,...,2N

Assuming that the nth approximation is known,find
the (n+l)st approximation as the solution of the linear

two-point boundary value problem. To do this take a differ-

ential of (2.19) and obtain












ax


or iik+l 1 ik =


a F1-(k I (xm(k+l)-xm(k)),
m=1 m (2.22)

i = 1.2....


From (2.19) above


i k =
Therefore
ik+l 1
Xj


(2.23)


k(k) a3(xk) (x (k+l)_ (b) (2.24)
mi=l axm xm


Note that the terms on the right-hand side of equation
(2.24) are the first two terms in the power series expan-

sion of the function i( (k+l))about the point xk.
The computational scheme consists in the numerical

production of a particular solution and two independent
homogeneous solutions of equation (2.24) on the interval

(otf) and the determination of the constant multiplier
of the homogeneous solutions. The solution may be repre-
sented in the following form:


xk+l(t) = p (t) + ci hi(t)
i=1


(2.25)


(2.21)


%t( k)










where the vector function p(t) is the particular solution

of equation (2.24), i.e.:
-k
P= i ax k (p xm(k)) (2.26)
m=l 3xm

Subject to the initial conditions

p(o) =o

and the vector function Fi(t) is the solution of the homo-

geneous form of equation (2.24), i.e.:


ti = hi (2.27)
m"1 "axm

The initial vector has the ith row unity and the others

zero. These vectors, p(t) and Fi(t), are all considered to
be determined on the interval o < t < tf.
The solution of equation (2.24)i.e. equation (2.25)
together with the boundary conditions in equation (2.20)

leads to the following linear algebraic equations:

x+1(o) = b4= p(o) + 1 cjhj(o) i-=. 2 (2.28)
3=1

and x+l(tf) = b,= p(tf) + cjh;(tf), .-.-i-, (2.29)


In equation (2.28) since p(o) and Hi(o) are known
the scalars ciappearing in equation (2.28) and (2.29)
are determined in terms of the boundary conditions.










Similarly the function k+2 is determined from the
-k+I
known values of the function xk+. In this manner a
large system of non-linear differential equations subject
to the boundary values may be resolved.

B. Formulation of Specific Problems

1. Cost Function 1

S= 1/2 j Pc2(t)dt subject to the constraints that
0 n(tf) = ano

TPBVP: By using the statement of the maximum principle
as in A.2, write the Hamiltonian as:

H(n,c,pl,P2,p ,t) = 1/2 c2 + Pl(t)n+P2(t)c
= 1/2 pc2 + Pl(t) L P B )n+ xc]

+ p2(t)[ n c] (2.30)

Then by (2.8)

= H(n,c,pl,P2,p ,t)
pl (t)=--- (2.31)
1 3n

and
bH(n,c,Pl,P2, P ,t)
P2(t) =- (2.32)

From equations (2.31) and (2.32) get the auxiliary equa-
tions


1 -
Pl(t) = (P-B)PI(t) T P2(t)


(2.33)










and

P2(t) = .Pl(t) + .P2(t) (2.34)

Setting = 0, yields


Pin (2.35)

The following system of equations are therefore obtained.

n = ( Pc ) n + .c (2.36)

= n .c (2.37)

l = ( c -B P P2 (2.38)

P2 = .P + )P2 (2.39)
Pn
PC --- (240)

with boundary conditions:

n(o) = no
C(o) = c
n(tf) = an, (2.41)

P2(tf) =0.0









Simulation of problem 1

The simulation of the problem as described above, for

nuclear reactor system dynamics, was performed with an

Applied Dynamics analog computer. A detailed computer

wiring diagram used for programming the computer is

given in Figure 41 Table 1 gives a list of system

constants and the initial conditions used. The system

equations described in Section 1, including the constants

of Table 1 are given below.

Problem 1:

n = 1000 pcn 6.4n + Ic (2.42)

S= 6.4n .lc (2.43)

Pi =-1000 PcP1 + 6.4 P1 6.4P2 (2.44)

P2 = .1 P1 + .1 P2 (2.45)
Pc =-1000 nP1 (2.46)

With boundary conditions:

n(o) = 5.0

c(o) = 64n,

n(tf) = 5.0a

P2(tf) =0.0 (2.47)













TABLE I


NUCLEAR REACTOR CONSTANTS FOR A LUMPED SYSTEM AND INITIAL
CONDITIONS




Nuclear Reactor Constants
----------------------------------------------------------

Decay constant, 0.10 sec-1

Average neutron life time,A 103 sees

Delayed neutron fraction,B 0.0064



Initial Conditions
..........................................................

Meutron density or power level, no 5K.W

Precursor population, co 64no











Simulation of the system of equations on an analog

computer requires amplitude scaling. This scaling is

necessary because all voltages must be kept within + 100

volts in order to prevent overloading of the machine

components. This overloading will lead to inaccurate

results. The basic set of scaled equations are rewritten

here from Appendix A for convenience. The bars indicate

scaled quantities.

Problem 1:

N =-100 N1P 6.4 N + 10 C

C = 0.064N 0.1 C

P1= 100 N P + 6.4 P1 64 P2


P2 = 0.1 P + 0.1 P2

Pc =-0.1 N P1

With boundary conditions:


S= 0.05

Co = 0.032


V(tf) = aNo


P2(tf) = 0.0










Typical solution procedure: With the equations scaled, the

modified gradient method will be used for optimization on

the analog computer and the results verified with the num-

erical technique of quasilinearization.

The gradient technique is preferred as opposed to

the trial and error procedure which can be quite cumber-

some and extremely time consuming, particularly for high

order systems. A brief description of this technique is

given in Section A.4.

The equations used for simulation are on page 21,

and the circuit diagram is shown in Figure 41,in Appendix

A. The neutron density or the power level is forced to

different desired values ranging from nd = 2n to nd = 10no

in one second. Figure 1 shows the neutron density plots

and the reactivity function plots. Table 2 shows the

initial values of the adjoint variables and the cost

function to attain various peak power levels. It is ob-

served that n(tf) =O for all t > tf. To accomplish this,

the reactivity will have to be programmed differentially

for all t > tf. This can be easily determined from the

kinetic equations ( 10,11).

The cost function is plotted versus the peak power

attained. This is shown in Figure 2. The graph indicates

increase in cost with increase in power level and this in

general is to be expected.







Nd=0no Nd= no \














o



4
0 Nd= 5no










\ o
40















4 0 0

0 Ndd 2 no





So o
50
1 ^ Nd [O no
40-





















0N1------- _______ _..
-J30

















.2 .4 .6 .8 10
TIME (Sec.)
Figure 1. Power level N and reactivity Q versus time
for ptirral control problem 1. Nuclear
Reactor System)

































0








F-
Z







0







0.5


2n% 4no Sno 8no IOr
DESIRED POWER LEVEL, Nd

Figure 2. Cost function versus desired po'.er level








TABLE 2

INITIAL CONDITICr"' ON THE ADJCINT VARIABLES AND THE COST
J. TO ACHIEVE DESIRED POWER LF"PLS


Nf J P10 P20


2no 0.454x10-5 0.629x10l5 0.408xl0-5

5no 1.263x10-5 0.631x105 0.268x105

6no 1.51 x105 0.552x105 0.235x105

1Ono 1.945x105 0.457x105 0.176x10-5




The gradient technique is studied in greater detail

for only one case, namely that of nd = 10no. A procedure

is developed whereby changing Pi, a search for Pi* which

produces the desired optimal control u* (reactivity in

this case) and trajectory n* can be made. Actually, if

all Pi values are scanned, much more general information

concerning the control of the system may be obtained.

In particular the reachable set of the system may be

determined. The reachable set R(t,no) is the set of all

states which may be obtained at time t (starting from no

at time t = 0) by means of admissible controls. Thus

R(t.no) determines the extent of possible trajectory

with admissible controls. By taking a suitable set of










Pi's enough boundary points may be obtained to define
R(tf,no), (tf is a particular value of t, one second in
this case). This can easily be seen in Table 3. Run
c-12A for instance is plotted in Figure (3). Converg-
ence of the trajectory to the desired one is shown.
Only the last four iterations are plotted. Correspond-
ing control functions are also shown.

Updating Circuit: The equation for the iterative pro-
cedure used are obtained from Equation (2.18), and the
computer circuit is given in Figure 44 Appendix A.
The equations are:

Pljl = P (o) + KI1[ Njt ) Nd(tf)]

+ K12 p (t )-pd ) (2.48)

P2j+l p2j() L K21[ N(t) Nd(tf)]
SK22[ ptf) P2d(tf)-] (2.49)

The above search algorithm seems to be the logical
choice since it is desired to satisfy the following end
conditions:

N(tf) = 10no
P2(tf)= 0











The last updated values of the Pij's for which the

above conditions are satisfied are held on the n analog

memories (number two integrators of the track-and-hold

system. Figure 44 ). The size K11, K12, K21, K22 of

the increments 8Pi are selected by the potentiometers.

The convergence of the iterations is very sensitive to

the values of K's. Larger values of K result in large

overshoots in the desired trajectory (see comment column

in Tablel2). Smaller values result in an increased

number of iterations to converge Figure 4 The

tolerance on the end conditions during these tests are

of the order of a few volts. No attempts are made to

obtain high accuracy in the terminal conditions as this

would require use of more analog and logic elements and

this is not warranted at this stage. Figure 3 shows

the sequence of trajectories converging to the optimal

one. In a typical case with Pi(o) far from optimal,

convergence was obtained in twenty-nine trials. Total

time for convergence was 60 seconds. A graph of number

of iterations for convergence versus Pl10 is shown in

Figure 4 for various step sizes. Note that when P10

is very close to the optimal value, the number of itera-

tions for convergence is infinite; Figure 5










shows the plot of the maximum value of n(neutron density)

and the value of P20 versus the iteration number in a

particular iteration cycle. With successive iterations

the value of P20i decreases until an optimal value is

obtained, i.e., value for which n(tf) = C0no and P2(tf)

= 0. Note that in C-12A since one is closer to the

optimal case, the change in P20 is gradual as also with

maximum n. Moving away from optimal setting of P20, the

change in P20 is much more and this is similarly trans-

lated to maximum n.

The steps in the iterative procedure are sufficiently

complex that it is wise to break the programming down into

a series of sub-problems.

An iterative procedure consists of the following steps:

(a) Choose an initial PlJ(o) to be the same as
that obtained by trial and error. Let P2j(0)

take on arbitrary values.

(b) Notice response of the power level (n). This

will help to determine sign of B P2. During

the iteration process, the addition is per-

formed by the analog memory so that the new

value


P2 201 + 8P21










where


SP2 = 1 K21 n(t) nd(tf)Ji K22 [P(tf) P2d(tf)]

can be transferred to the P20 memory, when up-
dating of this parameter is required. The
values of K21 and K22 are chosen to satisfy the
tolerance and oscillation conditions mentioned
earlier.
(c) Now reverse the above procedure. Chcose an
initial P23(O) as obtained by trial and error.
Let P1J(0) take on arbitrary values.
(d) Observe response of power level to determine
sign of hP1. During the iteration process
addition is performed according to the follow-
ing equation

Pj+l = P hP1


Procedure in (b) is repeated for suitable
selection of K11 and K12.
(e) Having determined the signs of hPl and hP2,
P13(o) and P23(o) are now assigned arbitrary
values and convergence to the true values is
obtained by observing the neutron density re-
sponse during each operation cycle to the end
of the iterations. The computer circuit dia-
gram for updating P3 is shown in Figure




30

8 --2-




'0 6 1 /
/ ITERATION

5 a


-

[ 3
0





2 2




40 DESIRED TRAJECTORY


2727h
-J 0


o th
-J 20t








10 0 o
20 -





2 4 6 LO
TIME (Sec.)

Figure 3. Power level and associated reactivity
time plots, leading to the desired
trajectory using gradient technique.
(Nuclear Reactor System)


















































Figure 4. Number of iterations for convergence versus
P10












STEP SIZE
-# I 11 K 12 K 13 K 14
S .01 .01 .01 .01
2 .02 .02 .02 .02
50 3 .03 .03 .03 .03
4 .04 .04 .04 .04 0





4o o
0


I-








S30









"20





0 .








0
5 15 25 35
010
C-)
< / i
0 l









28

26

24

22
"-o0 C-13A
20 C1 C-12A

18 C-15A

16

14

12

.0
0 5 10 15 20 25 30 35

50 ------- -- -- -- 0%

C-12A
C-14A C-15A /C-13A

40




c 30




20








0
5 0 5 20 25 30 35

ITERATION NO

Figure 5. Maximum value of N and P20 versus iteration
number in an iteration cycle








Problem 1 "ia Quasilinearization Technique: The quasi-
linearization technique is described on page 14 and will
now be applied to this problem.
For convenience rewrite equations (2.42 -2.45) and
modified by (2.46) as:

n = 5x105pln2-6.4n+0.lc (2.50)

c = 6.4n 0.1c (2.51)

P1 = 5xl05Pl2n+6.4P,-6.4P2 (2.52)

P2 =-0.1P1+0.1P2 (2.53)

Linearizing the above equations, as on page 15, there results

hk+l = hk + 5x105 [plk2nk(nk+lnk)+(nk)2(plk+l-plk)

6.4(nk+l-nk) + 0.1 (ck+l-ck)

But

nk = 5xl05pk(nk)2 6.4nk + 0.1 ck

Therefore

hk+l = 5x105plk(nk)2-6.4nk + O.1ck + 5x105 [p1k2nk

(nk+l nk) + (n)2lk+l k)] 6.4(nk+l-nk)

+ 0.1 (ck+l_-k)









Hence

nk+l = 106 k nnk+l 106plk(nk)2 + 5xlO(nk)2plk+1

6.4nk+l + O.lck+1 (2.54)

Now let n = xl, c = x2, Pl = x3 and P2 = x4, and let
k = o and this leads to

x11 = 106xox11-106 3lO o(xO)+5x105x3 (x10)2

6.4x11+0.1x21 (2.55)


Linearizing Equation (2.51) there results

ak+l = k + 6.4(nk+l-nk) 0.1(ck+l-ck)

= 6.4nkC-0.ck+6.4(nknk+l-k) 0.(ck+l-ck)

or ck+1 = 6.4nk+l O.ck+ (2.56)

let k = o,n = xl, etc, and get

x21 = 6.4x11 O.1x21

Repeating above procedure for equations (2.52) and (2.53)
obtain:

31 = 0l6x1(x30)2 106xlX3 ox31-5xlO5(x30 2xl

+ 6.4x31 6.4x41 (2.57)









and x 41 = -0.1x31 + 0.1x41 (2.58)

The above equations may be written in the following form:

x = bl + all1 + a2x21 + a13x3 (2.59)

21 = a21x + a22x21 (2.60)

X31 = b3 + a31x1l + a33x31 + a34x41 (2.61)
X41 = a43x3 + a44x41 (2.62)
AN+l N N+ N+1
In vector notation x = A(t) N + b(t) (2.63)
where
x = x b = b

x2 b2
x3 b3
x4 b4

and the components of the A matrix and b vector are:

b = 106(x0)2 x 30

b2 = b4 = 0

b3 : 106x0(x3o)2

all = 106 x 10x3-6.4

a12 = 0.1

a13 = 5x10 5(xo)2











a14 = a23 = a24 = a32 = a41 = a44 = 0

a21 = 6. 4
a22 = 0.1

a31 = -5x10(x30)2
a33 = -106xlx3o + 6.4

a34 = -6.4

a43 = -0.1

a44 = 0.1



It should be noted here that in vector equation (2.63),

the coefficients are time varying. Equation (2.63) may

also be rewritten as

n
i = aji(t)xj + b'(t), i = 1,2....,n. (2.641
j=l

If p(t) is the particular solution of equation (2.63). then
n
p = aji(t)pJ + bi(t) (2.65)


subject to the initial conditions p(o) = 0
or

p1 : h 4 aJlP1 + al2P2 + al3P3 (2.66)


P2 = a21i1 + a22P2


(2.67)










p3 = b3 + a31 + a33 (2.68)

P4 = a43P3 + a44P4 (2.69)

Appendix B shows incorporation of equation (2.65) in the

main programme.
The homogeneous solution of (2.64) is the solution of
equation

n
xi = ~ a(t) x i = 1,2,.....,n (2.70)
j=1

or X = A(t)5 (2.71)

Pontryagin (20) points out that the solution of the above
equation with time varying coefficients is

x(t) = 0 (t to)x(to) (2.72)
where

Ot) = A(t) o(t), 0(o) = I (2.73)

where o(t) is an nxn matrix
A(t) is an nxn matrix
I is a unit nxn matrix











The complete solution is therefore

x(t) = D (t)x(o) + p (2.74)



The following initial and final conditions were used



no = 5 k.w or 0.5 n (t=tf) = 1Ono = 5

co = 64no = 32 p2(t--tf) = 0


Pl0 = 0.5 x 10-5 Pl(tf) = 0.30xl0-5 (2.75)

p20 = 0.18 x 10-5

The initial guesses for the adjoint variables, i.e. pl0

and p20 were taken from the analog computer results

described earlier. The above case is for: neutron life
-3
time \ = 10cs
secs'
The numerical results of the first control problem

may be summarized in the graphs. The graphs also indicate

comparison with the analog computer results, which in

general is very good. Two cases were run average
neutron life time, = 10-3secs and A = 10-4secs. A

typical table, i.e., Table 3, shows the results of the

second case only. Comparison with analog computer results

is the subject of Figure (6).



















c
0

M n














or
C x









- E
Cp




C
--4
+














OI
+'
C







*C







aC
0


if


o r
oaoo0


0
X
o r"

0 C' 0


0
x
00
In '
(0
0


Io
0

So X
S0 r-


o C o




o c
4f) o "0
ocuo


0 0





41




The initial and final conditions used in Figure 7

are:


pi(tf) = 0.0000020, p2(tf) = 0.0

n(o) = 0.5, c(o) = 32 and A = 10-3

Figure 8 shows power level versus time plot of run

C-12A of Table 12. The initial conditions on the adjoint

variables in this run are: pl0 = 0.54x105 and P20 =
0.18x105. Convergence in this case is obtained in the

second trial.
















- Quasilineorizotion Data Points

-- Data Points From Analog
Computer Results


0 .02 .04 06 OB 0.1

TIME SECSS)

Figure 6. Power level versus time ( A = 10-4)
Comparison of analog computer re-
sults with those obtained by quasi-
linearization (Nuclear Reactor
System)






















50 -o- Quasilinearization Data Points

-r- Data Points From Analog
Computer Results


40





S30
-J

nii
.0
o-20





10 -
I0





0 1 I I I t I I I
0 .2 .4 .6 .8 1.0
TIME SECSS)

Figure 7. Power level versus time ( A -- 10-3)
Comparison of analog computer results
with those obtained by quasilinearization
(Nuclear Peactor System)









50




40 st Trial By Quasilinearization-

Lost Trial By Quasilinearization /




-J


o /
S. 20 0



-29th Iteration By
SGradient Technique On

I0
^'' Anolog Computer





0 .2 .4 .6 .8 1.0
TIME SECSS)





Figure 8. Power level versus time ( A = 10-3)
Comparison of analog computer results by
gradient technique, with those obtained
by quasilinearization (N:uclear Reactor
System)









2. Cost Function 2


J = 1/2 [pc2 + a (n-nd)2 dt


The two-point boundary value problem resulting by
the application of the maximum principle is:

S= 1000 pcn 6.4n+0.1 c (2.76)

c = 6.4n 0.1 c (2.77)

P1 = a(n-nd) 1000 pcP1 + 6.4 P1-6.4P2 (2.78)
P2 =-0.1 P1 + 0.1 P2 (2.79)

Pc =-1000 nP1 (2.80)

with boundary conditions:

n(o) = 5-0
c(o) = 64no (2.81)
n(tf) = 5-Oa
P2(tf) = 0.0


The cost function of this problem may be written as
the sum of the two cost functions as shown.


J = 1/2 [c2 + a (n-nd)2] dt = 31 + 32
0 / ~op2+











where = 1/2 pc2 dt
0


and may be thought of as the cost of the control action and


tf 2
J2 = 1/2 J a (n-nd)2dt
0


is the cost of the deviation of the output of the dynamical

system from the desired value.


Simulation of Problem 2


This problem reduces to problem 1 for a = 0. A graph

of the cost function for various values of a is shown in

Figure 9 Neutron density and reactivity function plots

are shown in Figures 10,11 for 5 (0,20) with the final

time tf fixed at 1 sec. It is clear from these figures

that for increasing values of a the rate of control rod

withdrawal is also increased. However, in order to limit

the power level to a certain desired value the

control rod is inserted back into the reactor at a faster

rate, thus shifting the reactivity function peak to the

left.










Table 5 shows the peak power attained, a and the cost

functions J1, J2 and J. With appropriate scaling J2
is given by:


J2 = -- a 0.1 (N-d)2dt
2 =0


J2 = 5 x 10-4 J21

where

tf 2
2' = a J O.l(N-Nd) dt
o

From the graph of J1 versus a (Figure 12) it appears

that the cost function increases with increasing a This

appears reasonable as increased b would decrease the error

between n and nd and hence more control would be needed.

The total cost function increases rapidly with an increase

in b as is apparent from the graph of the total cost J
versus Nd the peak power attained (Figure 13). and also

from graph of J versus d This result is again quite

logical as one would expect to pay a higher cost in order

to drive the error between n and nd to zero at a faster rate.































5
4 5


PEAK POWER DESIRED, Nd
Figure 9. Total cost function versus peak power
desired (control problem 2)












0 5 -" 2

4
=5




















50
40
3 a=0--\Q













































time for problem 2 for various values
20




I0





0.2 0.4 0.6 0,8 1,O
TIME (SEC)




of and Nd = 10 no (Nuclear Reactor
System)



















4-







SI ,I I I I

0.2 0.4 06 08 10
TIME SECSS)


50




40






OJ


2








II 0 -i
LlL





0.2 0.4 0.6 0.8 1.0
TIME SECSS)

Figure 11. Power level N and reactivity e ver-
sus time for control problem 2 for
various values of and Nd = 5no
(Nuclear Reactor System)











TABLE 4

COST FUNCTIONS AS A FUNCTION OF WEIGHTING a AND POWER
LEVEL Nd FOR PROBLEM 2



Nf a = 5.0
2J1 J2 J=Jl+J2

2no 0.456x10"5 0.0149x10-5 0.243x10-5

5no 0.696x10-5 0.527x10-5 0.875x10-5

6no 1.323x10-5 0.823x10-5 1.484x10-5

1Ono 1.914x105 2.88x1075 3.84x105
---------------------------------------------------------.
a = 10
..........................................................

2no 0.456xl0-5 0.0297x10-5 0.258x10-5

5no 0.696x10-5

6no 1.323x10 -

10no 2.058x10-5 4.95x10-5 6.0x10-5

a= 20
..........................................................

2no 0.456x1i05 0.0595x10-5 0.288x10-5

5no 0.745x10-5 1.92x10-5 2.30x10-5

6no 1.445x10-5 2.81x0-5 3.53x10-5

1Ono 2.126x10-5 9.2x10-5 10.26x10-5






















1.2









1.0




0.9 -




0.8 -




0.7




0.6



0.5 -- I- i 1- -!i l il _l_--1 --
0 5 10 15


Figure 12. Cost function associated with control rod
movement versus Nd = 10 no (control
problem 2)




















































Figure 13. Total cost function versus for control
problem 2 and M!d = lOn









3. Cost Function 3


J = 1/2 c2 + a (n-nd)2 dt


pc = a + A(n-nd)

The two-point boundary value problem resulting by
the application of the maximum principle is:

1 = 1000 an + 1000 A(n-nd) 6.4n + O.lc (2.82)

S= 6.4n 0.lc (2.83)

A = 0. a = 0 (2.84)

p, = (A2 + a )(n-nd) aA 1000 Pl(a+2An-And)
+ 6.4P1 6.4P2 (2.85)

P2 = 0.1P1 + 0.1 P2 (2.86)

P3 = A(n-nd)2 1000 npl(n-nd) a(n-nd) (2.87)

P4 = -a A(n-nd) lO00npl (2.88)

Pc = a + A(n-nd) (2.89)

The boundary conditions are:
n(o) = 5.0 P2(tf) = 0
c(o) = 64no n(tf) = 5.0 a (2.90)
P3(o) = P4(o) = 0 P3(tf) = P4(tf) = 0











Simulation of Problem 3


The equations used for simulation are given in

Appendix A. and the procedure used is similar to that

followed in Problem 1. The circuit diagram is shown in

Appendix A. The trial and error procedure is adopted in

the sense that the constants "a" and "A" are varied to

obtain the desired power level response. Figure 14 shows

the variation of "A" with "a" to achieve different power

levels ranging from 3no to 8no, Corresponding to these

power levels, the desired control is obtained. The initial

conditions of PI and P2 and the weighting factor a are now

varied to satisfy the end conditions on the adjoint vari-

ables, namely P2(tf) = P3(tf) = P4(tf) = 0.


It is observed in this investigation and justified by

solution of equation: N = AN(N-Nd) 0.64N + C that the

desired power level is reached only for an infinitely

large value of A. Figure (15). The power level response

also suggests a strong dependence on the desired level.

These observations therefore suggest that the form chosen

for this control will not yield a desirable closed loop

control. The form that accomplishes this is given in

Problem 4.







21

20

19

18

17

16

15

14

13

12




S10

9A




7 o

6

5




3

2





I 2 3 4 5 6 7 8 9 10 11
A
Figure 14. Variation of parameters in the con-
trol law for s.o.c. problem (i.e.,
problem 3) and for various values
of ild

















9n,


LUSn

U
ir



-j
LLJ

LUJ
-J
Ir




L'
6n0

4


4no,


0 2 4 6 8 10
VALUE OF A


Figure 15. Peak power level reached versus A (control
problem 3)









4. Cost Function 4


J =1/2 C [ 2 + a (n-nd)2 dt

tf
PC = A(n-nd) + B 1 (n-nd)dt

tf
Let x = J (n-nd)dt Hence x = (n-nd). and p = A(n-nd)+Bx
0

The two-point boundary value problem resulting by the appli-
cation of the maximum principle is:

n = 1000 A n(n-nd) + 1000 Bxn 6.4n +0.1c (2.91)

c = 6.4n O.lc (2.92)

x = n-nd (2.93)

A = 0 B = 0 (2.94)

PI = (A2+a )(n-nd) ABx 1000 P1 [2An-And+BxJ
+ 6.4P1 6.4P2 P3 (2.95)

P2 = -0.1P1 0.1P2 (2.96)

P3 = -B2x 1000 BnPl AB(n-nd) (2.97)

P4 = -A(n-nd)2 1000nPl(n-nd) Bx(n-nd) (2.98)

P5 =-Bx2 1000 xnP1 A x (n-nd) (2.99)

PC = A(n-nd) + Bx (2.100)










The boundary conditions are:

n(o) = no

c(o) = co

x(o) = 0

nl(tf) = nld =. ano


P2(tf) = 0

P3(tf) = 0

P4(tf) = 0

P5(tf) = 0

P4(o) = 0

P (o) = 0

Simulation of Problem 4


(2.101)


It was seen in Problem 3 that the power level re-
sponse was highly dependent on the setting of the desired
value of the neutron density. In order to remove this
restriction a control function of the type shown above.
namely a combination of proportional and integral control,
was tried. A simple transfer function analysis indicates











that such a system would be stable and the response would

be dependent on the coefficients A and B.

The equations used for simulation are given on page

63 and the procedure used is similar to that followed in

Problem 3. Circuit diagram is shown in Appendix A. Figure

16 shows the power level and the reactivity function

response for values of nd ranging fron 2no to 1Ono with

the final time fixed at one sec. Note in this case that

the power level rise is much faster, even though it settles

down to the desired value in fixed time tf. Corresponding

rod withdrawal is also much faster in the initial stages

and then gradual before it reaches a constant value, to

maintain the above level. Figure 17 shows variation

of cost function versus peak power attained. A graph is

also plotted showing behavior of constants A and B with

respect to peak power attained (Figure 18). It is clear

from this graph that there is no value of optimum A and

B which is independent of the desired power level. This

would imply that the control parameters A and B should

not be fixed but should be varied as a function of the de-

sired final state. Hence one would conclude that some sort

of adaptive control may be required. A control system is

said to be adaptive if it is capable of changing its con-

trol object in such a way as to operate at all times in

an optimal or nearly optimal fashion.











Figures 19 and 20 indicate the time response of
the auxiliary functions. Since these functions are highly

unstable, they are very sensitive to amplifier drift.
Because of this reason it is very difficult to satisfy and
maintain the end conditions, namely, P2(tf) = P3(tf) =

P4(tf) = p5(tf) = 0, and n(tf) = ano.




TABLE 5

COST FUNCTION AS A FUNCTION OF POWER LEVEL
Nt AND CONSTANTS A AND B FOR PROBLEM 4

(NUCLEAR REACTOR SYSTEM)


Nf A B a J

10no 2.123 42.12 3.178 1.201 x 105

6n, 2.709 55.50 5.667 0.495 x 10-5
5no 3.695 79.78 11.34 0.294 x 105

2no 3.836 128.4 0.105 x 10-5






Nd = 1Ono


Nd = 6 no



Nd= 5no

c-- -----

0 Nd = 2no


Nd = IOno










Nd = 6no o- -o



Nd 5no-

Nd= 5no


50






40




-j
> 30
LIJ
-J

w
_J



0
20






10


.2 .4


.8


1.0


TIME (Sec.)

Power level N and reactivity .c versus
time for specific optimal control pro-
blem 4 (Nuclear Reactor System)


Nd = 2n,
oo -----
0-~T*~
0-,~"~C~a


Figure 16.


I I I I I ~














13

12




10

0.9

0 0.8

0- .7

0.6

0.5

0.4

0.3

02

01 -


0 2 3 4 5 6 7 8 9 10 11

Nd
Figure 17. Cost function versus peak power
attained












A
4.0




3.5
B
130

120 3.0




100 2.5 -

90 -

80 2.0 -









40 -1.0 -

30

20 .5

10

0 0 1 L L
I 2 3 4 5 6 7 8 9 IOno
DESIRED Nd

Figure 18. Nd versus A and B









.10
















NdP lOno
S.05 P T.89













T =.90
.10
P3 P 4 .
T-.92 T=.9

















Nd= Ino
T=.90 T=.88
)F- 0






P 3 P4 T I.Osec

T \.94
S 0

Nd= 6no
Figure 19. Auxiliary functions, P2, P3, P4 and
P5 versus time for control problem
4; Nd = 1Ono and 6no (Nuclear Reactor
System)







.6
ST .90 sec.


.2

0
0

-r T=.9sec.

-2
-.4

-.6 P3

-.8


T= .84sec


.05
P5







.04 T= .O sec

.02 P4 --

-0

Nd = 5n,


T =.95 sec.


T= .8 sec


P5=0, 0 < t sec.


-- Nd = 2no

Figure 20. Auxiliary functions P2, P3, P4 and P5
versus time for control problem 4;
Nd = 5n, and 2no (Nuclear Reactor System)











Discussion of Results


Results obtained in the four problems attempted in

this chapter are analyzed and a comparison is given of

the two basically different control systems, i.e., 1) an

open loop control system and 2) closed loop control or

specific optimal control system.

The first two problems analysed represent open loop

control. Problem 1 is simulated in Section B.1 and the

operating time is set equal to one second. The perform-

ance criterion used here could be interpreted as an opti-

mization of control and movement, limiting erratic flux

changes (10).

Another useful performance index is the subject of

problem 2 wherein the integral of the error squared of

the reactor power and its desired value is minimized along

with the integral of the squared of the reactivity function.

In this problem the control function required to achieve

the desired state for example n(tf) = ano is dependent

on the value of a. The higher the value of a, the faster

the control action and the faster the power rise. This

variation is exhibited in Figure 10. Also it is observed

that the value of the cost function increases with in-

creasing value of a (Figure 12). This is to be expected

since the higher value of a would imply desirability of











driving the error to zero at a faster rate. Naturally more

control is needed to accomplish this.

Problems 3 and 4 represent closed loop or specific

optimal control. The specific forms chosen for the

controller are:


pc = a + A(n-nd)
tf
and Pc = A(n-nd) + B J (n-nd)dt



The performance index used is the same as in Problem 2.

The results of problem 3, point out a strong
dependence of the peak power attained to the desired

value of the power level. This result led to the invest-

igation of Problem 4. In this problem. one observes that

the cost function increases with increase in the desired

peak power (Figure 17), a result which is intuitively

obvious.

In problem 4. the control laws are determined to
achieve four different power levels and it is found in

each case that the closed loop or specific optimal control

suffers in performance in the sense that the cost function

is higher than in the open loop case (Figure 21). If the

values of the control parameters A and B are fixed as those

obtained for the case Nd = 10no and the reactor is brought

to various power levels, then the dotted curve in Figure 22

















































2no 4no 6no

DESIRED Nd


8no IOno


Figure 21. Cost function versus poak power attained
for open and closed loop optimal control
(Nuclear Reactor System)


1.8 I


1.4

1.2


0.6


0.4


0.2














4.0

[C.L.C.]

3.0 A = 2.13 x 0
B 4.25 x 10-3


2.0



1.0 -




0 2no 4no 6no Bno Ion0

DESIRED Nd






Figure 22. Total cost function versus peak power
attained for open loop control and cost
variation for fixed and optimally varying
values of control pararreters for closed
loop control (Nuclear Reactor System)










is obtained for the total cost function. Note here that the

cost function is higher in every case. These results com-

pare favorably with those obtained by Eisenberg and Sage

(21) power level trajectories and the corresponding time

variation of the control functions are shown in Figures 23

and 24. Figure 25 depicts the variation of the cost func-

tion versus the peak power desired for varying values of

Figure 26 shows the variation of the cost function with

variation in system parameters (neutron lifetime). It is

noted here that the closed-loop control is better than the

open loop control for large variation in lifetime about its

normal value.

It may be pointed out at this stage that the coeffi-

cients in the control law are highly dependent upon the

final time tf. the system parameters, and the initial state

of the system. The performance of the system, thus is a

function of the form chosen for the control law. The form

of the control law is very instrumental in determining how

well the performance of the s.o.c. compares to that of the

true optimal.

Figure 18 shows the variation of the parameters A

and B in the control law, as a function of the peak

power. Since this variation is considerable, adaptive

control is visualized as a means of further improving the

performance. If adaptive control is used, the example at

















SIOn,





8no





6no





4no





2no






0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

TIME SECSS)

Figure 23. Open loop and closed loop optimal
trajectories N(t) versus time
(Nuclear Reactor System)


* 0.3
2













S- lOn, [O.L.C]



10no [C.L.C]




























0 I I --- i I l-- I J
0 .2 .4 .6 .8 1.0
TIME SECSS)
Figure 24. Open loop and closed loop optimal control
function to obtain ten times the initial
power level in a Nuclear Reactor System








II



4 b




0.9




0.8




07




06




0.5




04




03




0.2




0.1
0 2 4 6 8 10 12
DESIRED POWER LEVEL Nd

Figure 25. Cost function versus peak power desired
























C.L.C]


[O.L.C]


.5 1.0 1.5 2.0
103 x NEUTRON LIFE TIME

Figure 26. Total cost function versus mean neutron
lifetime











hand could be that of the so called performance criterion

sensing or extremum adaptive system. An overall performance

index is measured and an attempt is made to find a control

law (by trial and error) which optimizes this performance

index. The learning states here are the parameters des-

cribing the control law, and the learning process consists

of the trial and error adjustment of these parameters. If

the learning states can be changed rapidly and accurately

the adaptive controller may be able to handle nonlinear

control systems the learning states will represent the

best instantaneous linear approximation to the non-linear

system. This approach is currently under study (22).

The gradient technique employed yields good re-

sults. The accuracy of the results can be considerably

improved at the expense of more analog and logic elements

and adjusting values of KOs. The convergence of the intera-

tions is very sensitive to the values of the K's. Larger

values of K's result in excessive overshoot about the de-

sired trajectory. Smaller values result in an increased

number of cycles to converge.

The Quasilinearization technique is next employed

to compute the results already obtained in Problem 1.

Figures 6 and 7 indicate close agreement between the

results obtained by the two techniques.











CHAPTER III


OPTIMAL CONTROL OF A NUCLEAR ROCKET ENGINE


Introduction and Model Description

The nuclear rocket engine to be considered is

similar to the solid-core, 25,000 megawatt, million

pound thrust model developed by Smith and Stenning (23,

24,25,26). The system of equations describing such a

model are based on a lumped parameter analysis. Such an

analysis does not describe the system completely and

furthermore the analysis is affected by a lack of agree-

ment among investigators concerning the temperature re-

activity effects. Stenning and Smith (27 ) contend that

temperature reactivity is directly proportional to the

square root of the core exit stagnation temperature, while

Mohler and Perry ( 2 ) hold that temperature reactivity

is directly proportional to this temperature. Weaver (12)

observes, that both contentions could give fairly accurate

results simply by choosing appropriate reactivity coeffi-

cients. Recent rocket engine tests (28 ) however indicate

that temperature reactivity may be proportional to square

of this temperature. All this adds up to the fact that the

nuclear rocket engines possess characteristics which need

exhaustive analysis. Thus, since the system is not well










understood, a closed-loop control is desired. The maximum

principle however, leads to an open-loop control. The

closed-loop control is a complex function of the state

variables and due to the analog computer limitations, the

treatment is restricted to open loop control only.

Basic components of the engine are depicted in

Figure 27 Liquid hydrogen, stored in a low pressure

tank, is pumped to high pressure, and is used to regenera-

tively cool the nozzle and reflector. The propellant gas

is further heated as it passes through the reactor core

and is exhausted through the nozzle. A portion of the hot

gas is bled off from the main stream and is used to power

the turbo pump, a feature giving rise to the name, "Bleed

turbine system." The bleed flow, and consequently the sys-

tem pressure are controlled by the valve in front of the

turbine.

The following equations describe the Smith-Stenning

nuclear rocket engine in terms of the time dependent nor-

malized quantities. An account of the normalization pro-

cedure is given in references (23) and (29) and therefore

will not be described here.

1 N' = P'~' N' + C' (3.1)
B


(3.2)





79

























o ~4,
C


















(v 0

mcy
-,






L4
o 4-

-'4
-4- 0






U.
-4. -'-4
4.4 441 o





C 04 4141"
.440 444
.0 O.,
4, 4, 0
f-.a 0.041
0.










TtTq = N pO(Tf)1/2 (3.3)

,p = p'V (T.)1/2 (P1 2
P (T )1/2 (3,4)


+ P (T')1/2 + h (3.5)
t ah TV


where N = C' -
Nd Cd

T P


VV V
0' =_ v d

and atTd ahPd
and at = a- BTd


Equations (3.1) and (3.2) are the reactor kinetics
equations. Equation (3.3) gives the time rate of change
of maximum core surface temperature, amd equation (3.4)
gives the time rate of change of the core inlet stagna-
tion pressure. Reactivity is given by equation (3.5)
Table 6 contains important constants for the nuclear
rocket engine. Note that the delayed neutron fraction is
smaller than the value typically associated with a uranium-
235 reactor. This reduced value was used by Smith and
Stenning to account for diffusion of some precursors into
coolant channels and their subsequent expulsion from the
system before emission of delayed neutrons.











TABLE 6


ROCKET ENGINE CONSTANTS


Design reactor power

Design core surface temperature

Design core inlet stagnation pressure

Thermal time constant

Pressure time constant

Control poison time constant

'alve time constant

Decay constant (one group)

Temperature reactivity coefficient

Hydrogen reactivity coefficient

Neutron lifetime

Delayed neutron fraction

Thrust, F

Specific impusle Isp


25,000 MW

5,000R

1.500 psi

0.68 sec

2.50 sec

2.00 sec

0.20 sec

0.10 sec-1

-10-6 oR-1

0.3 OR/psi

10-4 sec

0.005

1,500,000 lb

640 sec


~ _~I___~~


~_










The specifications in Table 6 significantly exceed
those of nuclear engines presently under development and
will probably be first achieved through clustering of small
engines. However, since the above system has received con-
siderable attention in the open literature, it makes it an
attractive system for an optimal control study.

It may be pointed out here that the reactivity
equation (3.5) is modified, particularly with respect to

the temperature coefficient of reactivity. The (T')1/2

term is replaced by T', as in the Mohler-Perry ( 2 )
formulation. In view of this modification, equation (3.5)
is re-written as :

p'
Pq = p + a( TV + a (3,6)
c t h T

The above system of equations with the constants
of Table 6 reduce to:

N' = 50 ( P'N'-NI + C') (3.7)
C' = 0.1 (NI C') (3.8)
T' = 1.471 (N' P'(T')1/2) (3.9)

P =0.4 (P' V' (T')1/2 ) (30)

P' = POC T' + 18 p2 (3.11)
Tu

The system is in an idle condition for T 0, the
idle or the initial condition being :To = 10000R and Po
275 psi. The value setting is fixed at 0.915. The full








power condition is T = 5000R and P = 1500 psi. It is
desired to optimize system performance such as to mini-
mize the following integral.


J = 1/2 tf P2(t)dt


with T(tf) = Td and pl < Pmax

This performance index can be interpreted as tending to
limit the maximum reactor period as is easily seen by
considering equation (3.7).
Since the effect of delayed neutrons is small in the
interval of operation, N'/N' = 50 P'. This indeed is the
inverse reactor period.
By using the statement of the maximum principle
as in Chapter II, the following two point boundary value
problem is obtained.
N' = 50( P N' NO + CO)
C' = 0.1(N' Cc)
iT = 1.47 (NO P'(T')1/2)
P' = 0.4(P'V'(T')1/2 () 2
(Tu)1/2
P = P'C TO + 18 L

Pi =-50 P'P1 + 50 P1 0.1 P2 1.46 P3
P2 =-50 Pl + 0.1 P2
P3 = 50 N'P1 + 900 ~LE- + 0.735 P'(T')-1/2 aP4




Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EUY1EAJF9_DGS0KF INGEST_TIME 2017-07-19T20:25:39Z PACKAGE UF00097875_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

COMPUTATIONAL TECHNIQUES AND PERFORMANCE CRITERIA FOR THE OPTIMUM CONTROL OF NUCLEAR SYSTEM DYNAMICS By JAGDISH K. SALUJA A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTUL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA December, 1966

PAGE 2

UNIVERSITY OF FLORIDA 3 1262 08646 415 2

PAGE 3

DEDICATED TO f'orr,, Dad, My Brothers and Sisters, and all my friends for their constant encourage^ ment .

PAGE 4

ACKNOWLEDGMENTS The author wishes to acknowledge the help of his chairman, Professor A. P. Sage, who has always been avail= able for consultation and who has been a source of continued encouragement and enlightening advice. The author also wishes to express his sincere appreciation to Dr. R.E. nhrig for his encouragement during the course of his pursuance of the doctoral degree and for being on his supervisory corrmittee. His help is greatly appreciated. The concern and interest of the remaining m^embers of the committee, Professors M.J. Ohanian, R.B. Perez, A.E.S. Green and R.G. Selfridge had a most beneficial effect upon the work's progress and is sincerely appreciated. The author also wishes to thank Dr. E.P. Eisenberg for his help with the digital computer computation, using the quasilinearization technique and to Dr. T.W. Ellis, Major J.T. Humphries, Mr. H. Atwater and Mr. Kesavan Nair for many helpful discussions. The cooperation and the assistance, especially cf Mr. Kenneth Fawcett on the Analog computer is greatly appreciated. Furtherm.ore , the author wishes to thank the department of Nuclear Engineering Science for the financial support, without which this work would not have been pos= sible. It is also a pleasure to extend thanks to Mrs. iii

PAGE 5

Mazelle Shadburn for her assistance in typing the disserta^ tion and to Mr. George Perry and Mr^ Gerald Murdock for the drawings. iv

PAGE 6

TABLE CF CCMTENTS Page DEDICATICN „ , , ii ACKNOWLEDGMENTS iii LIST OF TABLES , , . vii LIST OF FIQJRES . . . viii LIST OF SYMBCLS .................. xii ABSTRACT. , , . „ „ xiv CHAPTER I. INTRODHCTION , . „ . „ 1 II. OPTI^^AL CCNTRCL OF A NUCLEAR REACTOR SYSTEM. ....... 7 A. Introduction and Model Description^ . 7 1 Statement of the Maxirrum Principle. . 9 2 Statement of the Specific Optimal Control Problem (s.o.c) 10 3 Gradient Technique. . . „ 12 4 Quasilinearization Technique 14 B. Formulation of Specific Problems. . . 17 1 Cost Function 1 .... o ..... . 17 Problem 1 via quasilinearization technique 34 2 Cost Function 2 45 Simulation of problem 2 46 3 Cost Function 3. . „ 54 Simulation of problem 3 55

PAGE 7

TABLE CF CONTENTS (Continued) Pag( 4 Cost Function 4. ,..,.„.. . 58 Simulation of problem 4 59 C. Discussion of Results. ...... 67 III. CPTIN'AL CCNTRCL CF A NUCLEAR R(XKET ENGINE . . . . „ 77 Introduction and Model Description 77 Solution Procedure ....<,... 85 v^ariable start-up time, ..... 85 Variable desired temperature. <, . 86 Influence of Parameter variation on system response ....<,.. 94 IV. SUMMARY AND CONCLUSIONS ........ 99 Sumnary and results. . . „ o . . = 99 Suggestions for further study. -. . 102 APPENDICES o . . 103 Appendix A: Analog Simulation . . <, 104 Circuit Symbolism. . . . c . . . 104 Problem Scaling. ......... 110 Simulation by Gradient Technique 126 Appendix B: Computer Codes 138 Solution of TPBVP via Quasilinearization 138 Mimic Code 140 LIST CF REFER FNCES 145 BIOGRAPHICAL SKETCH ... 148 vi

PAGE 8

LIST CF TABLES Table Page 1. Nuclear Reactor Constants for a Lumped Systerr and Initial Conditions 20 2. Initial Conditions on the Adjoint Variables And the Cost J, To Achieve Desired Power Levels . . , , 25 3. Typipal Results of Quasilinearization Technique Applied to Problem 1 (Nuclear Reactor System). ........... 40 4. Cost Functions as a Function of Weighting a and Power Level N^ for Problem 2 . . . o » 51 5. Cost Function as a Function of Power Level f'f and Constants A and R for Problem 4 (Nuclear Reactor System) » . . . 61 6. Rocket Engine Constants « . . ^ . 81 7. Scale Factors for the Four Control Problems, 113 8. Potentiometer Settings for Control Problems 1 and 2 (Nuclear Reactor System) . . . „ c 115 9. Potentiometer Settings for Control Problem 2 (Nuclear Reactor System) .<,.,„,... 117 10. Potentiometer Settings for Closed Loop Control Problem 3 (Nuclear Reactor System) . . , . 119 11. Potentiometer Settings for Control Problem 4 (Nuclear Reactor System) ..„.....= 123 12. P^ «s as a Function of Iteration Step for Problem 1 (Nuclear Reactor System) .... 128 vii

PAGE 9

LIST OF FIGURES FIGURE PAGE 1. Power level N and reactivity ^c versus time for optimal control problem 1. (Nuclear Reactor System) „ » .« . . . . » c 23 2. Cost function versus desired power level ... 24 3. Power level and associated reactivity time plots, leading to the desired trajectory using gradient technique. (Nuclear Reactor System) „ . . . . ^ o » . . . o o , 4. Number of iterations for convergence versus PlO... o ........ e . ...» . 9. Total cost function versus peak power desired (control problem 2) .... = <,..<,... 11. Power level N and reactivity Pr versus time for control problem 2 for various values of a and Nj = 5n (Nuclear Reactor System) . 30 31 5. Maximum value of N and P20 versus iteration number in an iteration cycle .......<> 33 6. Power level versus time (A= lO""^) Comparison of analog computer results with those obtained by quasilinearization (Nuclear Reactor System). . ^ c . o , . . o c . . 42 7. Power level versus time ( \= 10"^) Comparison of analog computer results with those obtained by quasilinearization (Nuclear Reactor System). ... c ... <:.-<>».' •= ^-^ 8. Power level versus time (A= 10"^) Comparison of analog computer results by gradient technique, with those obtained by quasilinearization (Nuclear Reactor System). . . . . c 44 48 problem 2 for various values of a and Nj = 10 n^ (Nuclear Reactor System). . . . . » ^9 50 vlii

PAGE 10

LIST OF FIGirRES (Continued) FIGURE • PAGE 12. Cost function associated with control rod moverr.ent versus a N^j =10 n (control problenn 2) 52 13. Total cost function versus a , for control problem 2 and N^ lOn^ 53 14. Variation of parameters in the control law for s.o.c. problem (i.e., problem 3) and for various values of N(j >. . 56 15. Peak power level reached versus A (control problem 2) ....... o .... o . . 57 16. Power level N and reactivity p^ versus time for specific optimal control problem 4 (Nuclear Reactor System) 62 17. Cost function versus peak power attained , 63 18. N(j versus A and B o . « . . « 64 19. Auxiliary functions, P2. P3 , P4 and P5 versus tine for control problem 4; M^ = l.OnQ and 6nQ (Nuclear Reactor System) .... 65 20. Auxiliary functions P2 , P3, P4 and P5 versus tirtie for control problem 4; N(j 5no and 2no (Nuclear Reactor System) .... 66 21. Cost function versus peak power attained for open and closed loop optimal control (Nuclear Reactor System) . . . c » . . c 69 22. Total cost function versus peak power attained for open loop control and cost variation for fixed and optimally varying values of control parameters for closed loop control (Nuclear Reactor System) o . 70 23. Open loop and closed loop optimal trajectories N(t) versus time (Nuclear Reactor System) . . ^ . » o . . o . . . 72 ix

PAGE 11

LIST CF FIGIJRES (Continued) FIGURE PAGE 24. Open loop and closed loop optimal control function to obtain ten times the initial power level in a Nuclear Rocket System . „ 73 25. Cost function versus peak power desired. „ ^ 74 26. Total cost function versus m.ean neutron lifetime .... c ... , . o . 75 27. Schematic Diagram of a nuclear rocket engine 79 23. Control poison reactivity and tem.perature time plots, tf = 4 sees, 6 sees, and 11 sees (Nuclear Rocket Engine System). . . » 87 29. Pressure-timie plot, tf = 4 sees, 6 sees and 11 sees (Nuclear Rocket Engine System) . . 88 30. Optimal control poison reactivity and temperature-time plots, T(j = 0.4 (Nuclear Rocket Engine System) . . o . o o . » . c . . = . 89 31. Optimal pressure-time plot, "^ ^ — 0.4 (Nuclear Rocket Engine System) .o, c ,..,.. . 90 32. Optimal pressure and neutron density versus tire. T^ = 0.7 and 1.0 (Nuclear ^ccket Engine System) » . . » 91 33. Optimal pressure-time plot, T^^ = 0.7 and 1.0 (Nuclear Rocket Engine System) . 92 34. Cost function, J versus desired average core temperature (Nuclear Rocket Engine System) 93 35. Control poison reactivity versus time and average core temperature versus time for three cases of parameter variations (Nuclear Rocket Engine System). ... o ..<>.. . 96 36. Pressure versus tine response for three cases of parameter variations, (Nuclear Rocket Engine System) o ...... » 97

PAGE 12

IIST OF FIGirpES (Continued) FIG'TRP PAGE 37. Pressure versus time and average cere tenperature versus time for three cases of parameter variations (Nuclear Rocket Engine System). ... c ...<,......... , 98 38. Symbols used in analog computer diagrams. . . 105 39. Logic symbols used in computer diagrams . . . 107 40. Special circuits. ... ..... 109 41. Analog computer diagram for control problems 1 and 2, for a Nuclear Reactor System. . . . 114 42. Analog computer diagram for control problems 3 and 4 (Nuclear Reactor Systen , . 118 43. Computer control diagram. „ , . . , 134 44. Updating circuit. ..„..., 136 45. Box-oriented mimic diagram. ..,.„.... 141 46. Mimic instructions 141 xi

PAGE 13

LIST OF SYMBOLS n Neutron density. nj Desired neutron density. N Scaled neutron density. N^ Scaled desired neutron density. N' Normalized neutron density. c Delayed neutron precursor density. C Scaled delayed neutron precursor density. C» Noriralized delayed neutron precursor density. A Kean effective neutron lifetime (sec). P Core inlet pressure (psi). P(j Design core inlet pressure (psi). P» Normalized core inlet pressure. T Maximum core surface temperature (°R). 1^ Design maximum core surface temperature ( R). T" Normalized maximum core surface temperature. V Valve setting. V° Normalized valve setting. O Error weighting termo at Temperature coefficient of reactivity (°R ). a| Normalized temperature coefficient of reactivity (dollars ) . a^ Propellant density coefficient of reactivity (°R/psi) xii

PAGE 14

LIST CF SYMBOLS (Continued) a^ Morrralized propellant density coefficient of re< activity (dollars), /? Delayed neutron fraction. ^ Decay constant (sec ). P Total reactivity. P Normalized total reactivity (dollars). Pc • r Control poison reactivity. P'cT Normalized control poison reactivity (dollars) T Time constant (sec). Tp Pressure time constant (sec)o A,B Constants in the control law. A.E Scaled constants in the control law. H Hamiltonian. J Cost function. Pj^ i-th component of the adjoint vector F, tf Time for solution, final time. t Time. t^ Initial time. Uj_ Component of the control vector u. x^ Component of the state vector x. xiii

PAGE 15

Abstract of Dissertation Presented to the Graduate Council In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPITATIONAL TECHNIQUES AND PERFORMANCE CRITERIA FOR THE OPTIMUM CONTROL OF NUCLEAR SYSTEM DYNAMICS by Jagdish K. Saluja December, 1966 Chairman: Dr. Andrew P. Sage Major Department: Nuclear Engineering Sciences The main objective of this study is to investigate optimum control of the dynamics of nuclear systems and associated computational problems. Two systems are considered: 1) Nuclear Reactor System, 2) Nuclear Rocket Engine System. For the first system, three problems are analyzed: a) the neutron density or reactor power is trans= ferred from one steady state condition to another with a quadratic constraint on control rod movement. By varying the operating time and placing an inequality constraint on reactivity a form cf minimum time control is obtained, b) integral of quadratic form of error in desired neutron density is minimized with quadratic integral constraint on control rod movement, c) a fixed configuration optimum closed-loop control system is realized by letting control rod nrcvement be a proportional plus integral function of neutron density error. xiv

PAGE 16

The two-point boundary value problems (TPBVP) , resulting by the application of optimization theory to point reactor models are obtained. A modified gradient technique is used to exploit the hybrid capabilities of modern analog computers. The results are compared with the solutions of this problem by the numerical technique of quasilinear= ization. In all of these, the convergence is obtained in no more than four iterations. In the examples considered here, the cost function is higher for the closed loop case as opposed to the open loop case. This is in general agreement with physical insight and previous results, since an assumed form of closed loop controller is used. The coefficients A and B in the control law are highly dependent upon the final time tf, the system parameters, and the initial and desired final state. Variation in system parameters may, however, render the closed loop control better than the open loop control which does not change with changing system parameter; For the case of the nuclear rocket engine system, only a quadratic constraint on total reactivity is con= sidered. The nuclear rocket engine is similar to the solidcore 25,000 megawatt, million pound thrust model developed by Smith and Stenning. Here, the results obtained for various start-up times and various desired power levels in terms of control needed and required cost, appear very XV

PAGE 17

reasonable. Also it was observed that if the pressure and the temperature tirre constants were increased frotr their normal values , the temperature response became somewhat sluggish. Possible improvements to this situation are suggested. xvi

PAGE 18

CHAPTER I INTRODUCTION Just as all other fields of scientific thought have undergone a major change in recent years, so also has the field of rocketery. From the Chinese rockets of 1232 A.D. to those of the Germans during the second world war, man has made a "case" for nuclear rocket propulsion as a vital factor in the success of future space travels This has been pointed out in the recent literature (1,2,3,4) and accepted in the field today.Its acceptance is fur= ther strengthened by the successful testing of the KIWI and NERVA reactors. This sudden interest in space travel has led to thoughts about improved system performances, as a result of which a very instructive book on nuclear rocket system optimization by Bussard (5) has appeared. An outgrowth of this interest is the optimal control problem wherein the control is such as to give the "best performance". From the nuclear literature, however, it is clear that the subject of optimal control of nuclear systems has received little interests Recently a num.ber of papers have appeared in this area, ranging from optimal reactor shutdown programs for control of xenon poisoning (6,7) ,

PAGE 19

and synthesis of optimal nuclear reactor systems (8,9) to optimal control of nuclear reactor systems (10. 11 ) . More recently Mohler and Weaver have reported some work in the application of modern optimal control theory to the dynamics of nuclear rocket systems (ll.,12) but much still remains to be done. The mathematical theory used for the study of optimization problems is the calculus of variations. It is of interest to note that Goddard (l_3) recognized this technique as an important tool in the performance analysis of rockets in his early paper published 40 years ago. There are many possible variational methods fcr minimizing or maximizing a functional over a function space. Among the commonly used methods in control system design are: (1) calculus of variations (2) maximum principle (3) dynamic programming In all cases the object is to find the optimum control law such that the given functional of the performance indices is minimized or maximized. These methods are re= lated (5) in that the rr.aximum principle can be derived from the calculus of variations as well as from the method of dynamic programming^ Pontryagin's maximum

PAGE 20

principle is a very elegant method, which is applicable to a broader class of problems than the calculus of variations. An outcome of this principle is the so called two»point boundary value problem, abbreviated as "TPBVP',' with split boundary conditions. This method will be applied in the study of dynamics of nuclear systems. Since the dynamic response of nuclear systems is represented by a set of nonlinear, time-varying differential equations, conventional analysis is usually too difficult to give a closed form solution. Development of efficient computer techniques is therefore necessary. A numerical approach which is currently quite popular is the method of gradients (5). It proceeds by solving a sequence of sub-optimal problems with the property that each successive set of solution functions yields an improved value for the functional being optimized. Another numerical technique of considerable interest is the approximating process of quasilinearization which produces a sequence of approximations x^, n = l,2,..c that converge quadratically (14) to the solution of the ori= ginal two=point boundary value problem. This technique was initially suggested by Bellman and Kalaba (15) as a modification of the Newton-Raphson-Kantorovich method. This converts the problem of solving a two-point boundary non-linear differential equation into a problem of solving a sequence of two-point boundary linear differential equation

PAGE 21

problems. The solution of the latter equation is considerably easier to find. For the study of optimization problems, a precise formulation of the performance index is necessary. The performance index is assumed to be given by: ftf J(Xo,tf,u) (x,u,t)dt where is a scalar, a function of the state vector x, the control vector u and time. Corresponding to this performance index, a control law is found which is the best in the sense of maximiizing the performance index. The performance indices considered in this study are: f ( i) J = 1/2 t: dt a) n(tf) = n^i b) T(tf) = Td This performance criterion may be interpreted as an optimization of control rod movement, limiting erratic flux changes (10). Clearly it tends to limit the maximum reactor period. J = 1/2 ( ^ ^c2 + a (n=nd)^ dt L Integral of the error squared of the reactor power and its desired value is minimized along with the integral of the square of the reactivity function, a is the error weighting term ii

PAGE 22

iii) J = 1/2 j ^ [ ^c^ + a(n-nd)^ ] dt where ^c = a + ACn-n^) a and a are constants to be determined, iv) J = 1/2 ) ^ [ ^c^ + a(n-n.)^J dt where P^ = Atn-n^j) + b) (n-n^) dt o A and B are constants to be determined. The last two criteria give closed loop control and since a specific form is chosen for the controller, it will henceforth be referred to as the specific optimal control problem. The above four controls are analyzed and compared with one another for two reactor models. The two models analyzed are: (1) Point model of a nuclear reactor system. (2) Point model of a nuclear rocket engine. Chapter II contains the description of the first model. Chapter III concerns the second model and follows along

PAGE 23

the same pattern as Chapter II. Chapter IV includes a discussion of results and suggestions for further study.

PAGE 24

CHAPTER II OPTIMAL COKTRCL CF A NUCLEAR REACTOR SYSTEM A. Introduction and Model Description Need frequently arises in nuclear reactors as well as in nuclear power plants for power level changes. These changes are generally accomplished by semi-automatic means in small research reactors. In large power reactors and in nuclear rocket engines, however, automatic control is desired. If this control is optimal in some specified sense, economic or other advantages could be gained. The reactor model considered here is that of a one delayed-neutron group without temperature feedback and without spatial and angular dependence, described by the following system of differential equations: n = -^ (P /3 )n + /.c (2.1) c = -fn AC (2.2) where n and c, the neutron density and precursor population are the state variables and the reactivity ^j is

PAGE 25

the control variable. The dot represents derivative with respect to time. A , is the average neutron lifetime /. , is the decay constant, and /3 , is the delayed neutron function. The system is at steady state for t < o and the initial conditions are: n(o) = n^ and c(o) = c^. ' It is desired to optimize system performance such as to minimize the following integrals. 1) J = 1/2 f Pc^it) dt •^0 with n(tf) = n^ and I P^ I < Pmax 2) J = 1/2 [^ [Pc^ + o{n-n^f] dt in fixed time tj 3) J = 1/2 J^ \^ Pc2 + a(n-nd)2 1 dt where p^ = ^ + A(n=n^) in fixed time t£ 4) J = 1/2 J^ [ Pc2 + a (n-n^)^ ] dt where p^ A{n-nci) + B f (n-n^) dt in fixed time tf.

PAGE 26

1 . Staterrent of the K^axitnum Principle Consider the control processes to be described by a vector differential equation: X = f (x,u,t) , x(c) = Xo (2.3) where f(x,u,t) is a vector with coordinates f^(x,u,t), f2(x.u,t) fn(x.U,t). The control functions have to be piecewise continuous and can be restricted by the inequalities: gj(u^,...u^) < (j = l....m) (2.4) Assume a perforrrance criterion of the type: tf J= 0[x(tf).u(tf), ^^]+ /
PAGE 27

10 The variable appearing in Equation (2.6) is called the auxiliary vector or the adjoint vector or momentum function. The system equations can now be written in the following manner. aH(x.u.p.t) a Pi (2.7) Auxiliary or adjoint, or costate functionsp^^ are defined by (2.8) p.(t) = '^H(x.u,p,t) d ^x The Hamiltonian is now minimized with respect to the control variable u, such that HO(x.p,t) = Min H(x,p,u,t) (2.9) where u°(x,p,t) minimizes H and the symbol^x refers to the admissible control space. If rx is infinite (2.9) becomes -^ = 0. The boundary conditions are given by : x(to) = Xq (2.10) and P(t^) = V a (x.u.t) 2. Specific Optimal Control where ( Vx ^ )^ = t = tr be The desirability of closed loop control laws has led to the development of "specific optimal control."

PAGE 28

11 The specific optimal control (S.O.C.) problem is defined in the following manner (17) : given a plant with a state equation of the form X = T{x,u,t) , x(to) = Xq (2.11) It is desired to determine the unknown parameters in a control law of the form u = h{y,b) . b = (2.12) In the above equation y is a p-dimensional vector which is a known function of the state and F is a q-dimensional constant vector to be determined such that an index of performance of the form ^f J(u) = /• g(x,u) dt (2.13) Jo is minimized. The specific optimal control problem can be solved by a simple re-formulation of the above problem,, The reformulation is carried out by substituting equation (2.12) into equations (2.11) and (2.13), leading to x = f (x,h(y,b),t) tf _ _ _ J(b) = r g(x.h(y.b)) dt •'o

PAGE 29

12 Since y is a known function of x and b is a constant vector X = f (x.b.t) ^f J(b) = f g(x.b)dt (2.14) •'o and b = Thus the S.O.C. problem has been embedded in the problem of finding the vector b, which minimizes the cost function and differential constraints imposed by equation (2.14). The above is now in the form of the typical optimal control problem. Using Pontryagin* s m.aximum principle then allows one to solve the S.O.C. problem by solving a two-point boundary-value problem. 3. Gradient Technique This technique is a successive approximation procedure in which an iteration is made on the control variable u(t) so as to improve the function to be minimized at each new iterate, subject to the constraining differential equations. Consider, for exam.ple, those problems in which it is required to extremalize performance indices of the form

PAGE 30

13 J(x,u,t) = f (t> (x.u.t) dt (2.15) •'o subject to the constraint in the form of the first order differential equations 0i = Xi "^^(x.H.t) = o • (2.16) It is required to find u*(t) which at the same time extremalizes J and satisfies equation (2.16). The overall algorithm which simultaneously satisfies both conditions gives (i+l)st iteration of u in terms of the ith iteration: m '-^ = "^tK[^. ., l-^] (2.17) where m denotes the number of equations 9 and K and A are chosen arbitrarily. It was seen in the section that the maximum principle applied to the above set of equations, resulted in a coupling equation, coupling the state vector and the adjoint vector via the control variable u. This result therefore can modify equation (2.17) to yield the following algorithm which approximates some of the characteristics of (2.17), pi+1 _ pi + K [ pi (tf) P^(tf)] (2.18) This algorithm was used in the problem solved in section B.l. This is an excellent approximation to (2.17) which is most useful on a repetitive analog computer.

PAGE 31

14 4 . Quasllinearization Technique This technique, developed by Kalaba, Bellman and others is particularly suited for the type of optimization problems being considered here. It provides an effective computational tool for a wide class of non-linear two-point and multi-point boundary value problems (15.18.19). The technique consists of producing a sequence of approximating functions xO(t). x^^ ^ (t ) . . . , x^") (t ) which converge rapidly to the solution of the original problem. In (14.20) it is shown that this convergence is quadratic, i.e., asymptotically the number of correct digits in each approximation is doubled. To illustrate the technique, consider a 2N-dimensional dynamical vector equation: ^= |(>^) (2.19) subject to the boundary conditions Xi(o) = bi i = 1.2, ..,,N (^ \ u (2.20) Xi(tf ) = bj_ i = N+1,... ,2N Assuming that the nth approximation is known, find the (n+1)^ approximation as the solution of the linear two-point boundary value problem. To do this take a differential of (2,19) and obtain

PAGE 32

15 A7^'(^)

PAGE 33

16 where the vector function p(t) is the particular solution of equation (2.24), i.e.: P= ^i^^*") ^ ^ ^HJ (P-Xm^^)) (2.26) Subject to the initial conditions p{o) = and the vector function hj^(t) is the solution of the homogeneous form of equation (2.24), i.e.: m-1 o^m The initial vector has the ith row unity and the others zero. These vectors, p(t) and hi(t), are all considered to be determined on the interval o £ t < t£ . The solution of equation (2.24) i.e. equation (2.25) together with the boundary conditions in < quation (2.20) leads to the following linear algebraic equations: x't^-^(o) = b^= p(o) + E cjhj(o) , -L=i,^, (2.28) * j=l and xl^"^l(tf) = b,= pjtf) + E c^h^^itf), ^'-^^\y-,^^ (2.29) In equation (2.28) since p(o) and hi(o) are known the scalars c^ appearing in equation (2.28) and (2.29) are determined in terms of the boundary conditions.

PAGE 34

17 Similarly the function x is determined from the — k+1 known values of the function x . In this manner a large system of non-linear differential equations subject to the boundary values may be resolved. B. Formulation of Specific Problems 1 . Cost Function 1 J = 1/2 J Pc^(t)dt subject to the constraints that o n(tf) = ano TPBVP: By using the statement of the maximum principle as in A. 2, write the Hamiltonian as: H(n,c,pi.P2. p ,t) 1/2 p^^ + P2(t)n+P2(t)c = 1/2 Pc^ + Pi(t)[ ;\( P^B)n^ /.c] + P2(t)[ ;f n Ac] (2.30) Then by (2.8) aH(n,c,p, ,pp, p ,t) P, (t) = -^— ^ (2.31) -Idn and dH(n.c,Pi,P2, P ,t) P2(t) = ^-^ (2.32) From equations (2.31) and (2.32) get the auxiliary equations Pl(t) = ^ (p-/3)Pi(t) 4 P2(t) (2.33)

PAGE 35

18 and P2(t) = /.P^d) + ;.P2(t) (2.34) Setting -^ = 0, yields -^c = !l! (2.35) A The following system of equations are therefore obtained. r^ = -i Pc B ) r\ -^ /.c (2.36) 6 = -^ n AC (2.37) Pi = " A ^ ^c ^ ) Pi T ^2 (2.38) P2 = aP^ + /.P2 (2.39) P,n Pc =A (2.40) with boundary conditions: n(o) = no C(o) = c^ n(t^) = an^ (2.41) P2(tf) =c.C

PAGE 36

ly Simulation of problem 1 The simulation of the problem as described above, for nuclear reactor system dynamics, was performed with an Applied Dynamics analog computer. A detailed computer wiring diagram used for programming the computer is given in Figure 41 . Table 1 gives a list of system constants and the initial conditions used. The system equations described in Section 1, including the constants of Table 1 are given below. Problem 1 ; n = 1000 p^n-6.4n + , Ic (2.42) c = 6.4n .Ic (2.43) F^ =-1000 ^cPi + 6.4 P^ 6.4P2 (2.44) P2 = .1 Pi + .1 P2 (2.45) Pc =-1000 nP^ (2.46) With boundary conditions: n(o) = 5.0 c(o) = 64n^ n(tf ) = 5.0a P2(tf) =0,0 (2.47)

PAGE 37

20 TABLE I N'JCLEAR REACTOR CONSTANTS FOR A LUMPED SYSTEM AND INITIAL CONDITIONS Nuclear Reactor Constants Decay constant, /. OclO sec~^ Average neutron life time, A 10~'^ sees Delayed neutron fraction,^ 0.0064 Initial Conditions •^'eutron density or power level, n 5K.W Precursor population, c 64n_,

PAGE 38

21 Simulation of the system of equations on an analog computer requires amplitude scaling. This scaling is necessary because all voltages must be kept within + 100 volts in order to prevent overloading of the machine components. This overloading will lead to inaccurate results. The basic set of scaled equations are rewritten here from, Appendix A for convenience. The bars indicate scaled quantities. Problem 1 : N =-100 N^PjL 6.4 N + 10 C C r 0,064N 0.1 C P = 100 Np2 + 6.4Pj^-6.4P2 P2 = 0.1 P^ + 0.1 P2 Pc =-0.1 N P^ With boundary conditions: K = ^'

PAGE 39

22 Typical solu t ion procedure : With the equations scaled, the modified gradient rricthod will be used for optiaization on the analog computer and the results verified with the numerical technique of quasilinearization. The gradient technique is preferred as opposed to the tria] and error procedure which can be quite cumbersome and extremely tim.e consuming, particularly for high order systems, A brief description of this technique is given in Section A. 4. The equations used for simulation are on page 21, and the circuit diagram is shown in Figure ^1 , in Appendix A. The neutron density or the power level is forced to different desired values ranging from n^ = 2n to n^ = lOno in one second. Figure 1 shows the neutron density plots and the reactivity function plots. Table 2 shows the initial values of the adjoint variables and the cost function to attain various peak power levels. It is observed that n(tf) ^0 for all t > tf. To accomplish this, the reactivity will have to be programmed differentially for all t > tr. This can be easily determined from the kinetic equations ( 10,11) . The cost function is plotted versus the peak power attained. This is shown in Figure 2. The graph indicates increase in cost with increase in power level and this in general is to be expected.

PAGE 40

I. ,4 .6 TIME (Sec.) Figure 1. Power level N and reactivity Q^ versus tirr.e for cptirral control problori', 1. (Nuclear Reactor System)

PAGE 41

24 in o DESIRED POWER LEVEL, N^j Figure 2. Cost function versus desired pov;e: level

PAGE 42

25 TABLE 2 INITIAL CCMDITir^'S PN THE ADJCI^^T VARIABI ES AND THE COST J. TO ACHIE\'E DESIRED PCWER LF\'FI,S N^

PAGE 43

26 Pi's enough boundary points may be obtained to define R(t^,nQ), (tf is a particular value of t, one second in this case). This can easily be seen in Table 3. Run C-12A for instance is plotted in Figure (3). Convergence of the trajectory to the desired one is shown. Only the last four iterations are plotted. Corresponding control functions are also shown. Updating Circuit : The equation for the iterative procedure used are obtained from Equation (2.18), and the computer circuit is given in Figure 44 ^ Appendix A. The equations are: P^J^^ = PJ(0) 1 Kii[N^(t ) Nd(tf)] ^^4^^^f^-P2d(^f)J ''-''' p^J-^l = P2^(0) t K2i[ ^'^(t) Nd(tf)] 1 K22[ P^Ctf) -P2d(tf)-] (2.49) The above search algorithm seem.s to be the logical choice since it is desired to satisfy the following end conditions : N(tf) = lOn^ P2(tx) =

PAGE 44

27 The last updated values of the P^^'s for which the above conditions are satisfied are held on the n analog memories (number two integrators of the track-and-hold system. Figure 44 ). The size K^xKx2. K2X , K22 of the increments 8 P^ are selected by the potentiometers. The convergence of the iterations is very sensitive to the values of K's. Larger values of K result in large overshoots in the desired trajectory (see comment column in Tablel2). Smaller values result in an increased number of iterations to converge Figure 4 . The tolerance on the end conditions during these tests are of the order of a few volts. No attempts are made to obtain high accuracy in the terminal conditions as this would require use of more analog and logic elements and this is not warranted at this stage. Figure 3 shows the sequence of trajectories converging to the optimal one. In a typical case with P^(o) far from optimal, convergence was obtained in twenty-nine trials. Total time for convergence was 60 seconds. A graph of number of iterations for convergence versus Piq^ is shown in Figure 4 for various step sizes. Note that when P.„^ is very close to the optimal value, the number of iterations for convergence is infinite; Figure 5

PAGE 45

28 shows the plot of the maximum value of n(neutron density) and the value of Poq versus the iteration number in a particular iteration cycle. With successive iterations the value of P20 decreases until an optimal value is obtained, i.e., value for which n(tf) ICn^ and P2(t£) = 0. Note that in C-12A since one is closer to the optimal case, the change in P20 is gradual as also with maximum n. Moving away from optimal setting of Poq* ^^® change in P20 is much more and this is similarly translated to maximum n. The steps in the iterative procedure are sufficiently complex that it is wise to break the programming down into a series of sub-problems. An iterative procedure consists of the following steps (a) Choose an initial Pj^J(o) to be the same as that obtained by trial and error. Let P2'^(0) take on arbitrary values. (b) Motice response of the power level (n). This will help to determine sign of 8P^. During the iteration process, the addition is performed by the analog m.emoiry so that the new value ^2 " 20 2

PAGE 46

29 where SP2^ = t K21 rnJ(t) n^(tf) 1 K22rP^(tf) " P2d(tf can be transferred to the P20 memory, when updating of this parameter is required. The values of K22 and K22 are chosen to satisfy the tolerance and oscillation conditions mentioned earlier. (c) Now reverse the above procedure. Choose an initial P2-^(0) as obtained by trial and error. Let P,-^(0) take on arbitrary values. (d) Observe response of power level to determine sign of ^ Pj^ . During the iteration process addition is performed according to the following equation J + 1 10 + ^P, Procedure in (b) is repeated for suitable selection of K^^^ ^^^ ^ 12' (e) Having determined the signs of hP^ and i^P^, P2^-^{o) and P2"'(o) are now assigned arbitrary values and convergence to the true values is obtained by observing the neutron density response during each operation cycle to the end of the iterations. The computer circut diagram for updating P-" is shown in Figure

PAGE 47

30 L^^ I 50 40 DESIRED TRAJECTORY 20 10 ,^^.-...^,,:=,£L Figure 3. Power level and associated reactivity time plots, leading to the desired trajectory using gradient technique. (Muclear Reactor System)

PAGE 48

31 Figure 4. Number of iterations for convergence versus PlO

PAGE 49

60 50

PAGE 50

33 — ^/_ 10% C-I3A 15 20 ITERATION NO. Figure 5. Maximum value of M and p2Q versus iteration number in an iteration cycle

PAGE 51

34 Problem 1 ^la Quasilinearization Technique : The quasilinearization technique is described on page 14 and will now be applied to this problerr. For convenience rewrite equations (2.42-2.45) and modified by (2.46) as : (2.50) n = 5xl0^pin2-6.4n+0.1c c = 6.4n 0.1c P^ = 5xl0^Pi2n+6.4P^-6.4P2 ?2 =-0aPi+0.1P2 (2.51) (2.52) (2.53) Linearizing the above equations, as on page 15, there results pk+l = nk + 5x10^ 'pik2nk(n^^l-nk).(n^)2(p^k.l.p^k)j 6.4(n'^^i = nk) + 0.1 {c^^^-c^) But n*' = 5xl0^pi''(n'')^ 6.4n'' + 0.1 c^ Therefore k+1 = 5xl0^pi'^(n'^)^-6.4n'^ + O.lc"^ + 5x10' k^ k Pi 2n I k+l kv , / (n -n ) + ( n'^)'(Pl'^^-p/)] 6.4(n^*l-n' + 0.1 (c'^^l-c'')

PAGE 52

35 Hence 6.4n''^^ + O.lc'''^^ (2.54) Mow let n = x^, c = X2, p^ = X3 and p2 = x^, and let k = o and this leads to xi^ = 10^xOx^l-106x3O(xiO)2+5xl05x3^(x^°)^ 6.4x^^+0. 1X2^ (2.55) Linearizing Equation (2.51) there results ek+1 ^ -k + 6.4(n'^^l.n'^) OA{c^^^-c^) = 6.4n'<-C.lc»^+6.4(n''"'^-n'^) 0. Kc'^^'^-c'^) or c'<^i = 6.4nk^l Clc^+l (2.56) let k = o,n = XI, etc, and get X2^ = 6.4x^1 0.1x2^ Repeating above procedure for equations (2.52) and (2.53) obtain : ^3! = 10^xi°(x30)2 . 106x^°X3Ox3l»5xl0^(x3C)2x/ + 6.4x3! 6.4x4! (2.57)

PAGE 53

36 and X • .1 4-^ = -O.lxol + 0.1: (2.58) The above equations tray be written in the following form: ^1^ = b^ + aiix/ + ai2X2^ + ai3^3^ i + <3^ = b3 + 331X^1 + 333x3^ + 334x4! <4'^ = 343X3^ + 344X4! In vector notation x where •N-1 =A(t) x^^^l. b(t) ^^1 (2.59) (2.60) (2.61) (2.62) (2.63) X = b = and the components of the A inatrix and b vector are: 106(xi°)2 X 3° b2 = b4 = b3 r 10^i°(x3°)2 a^i = 10^ X 1^X3-6.4 ^12 = 0.1 ai3 = 5xlO^(xi°)^

PAGE 54

37 ^14 ^23 = ^24 = 332 = ^41 = 844 = 321 = 6. 4 322 =0.1 331 = -5xl0^(x3°)^ 333 = -10^x1^x3° + 6.4 334 = -6.4 343 = -0.1 344 =0.1 It should be noted here thst in vector equ3tion (2.63), the coefficients 3re time vsrying. Equ3tion (2.63) m3y also be rewritten 3S x^ = L a/(t)xJ + b^(t). i = l,2.._,n. (2.64} If p(t) is the psrticulsr solution of equation (2.63). then P^ = E 3/(t)pJ + bi(t) (2.65) subject to the initial conditions p(o) = or p^ = h^ + a^^n + ai2P2 ^ ^13^3 (2.66) P2 = ^21^1 ^ a22P2 (2.67)

PAGE 55

38 P3 = ^"3 -^ 331 "" ^33 ^^-^^^ P4 = a43P3 ^ a44P4 (2.69) Appendix B shows incorporation of equation (2,65) in the main programme. The homogeneous solution of (2.64) is the solution of equation i. = £ a.^(t) xJ . i = 1.2... .,n (2.70) ^ j=l ^ or 3c = A(t)x (2.71) Pontryagin (20) points out that the solution of the above equation with time varying coefficients is x(t) = (t tjx(to) (2o72) where i(t) = A(t) :'(t), -(o) = I (2,73) where c(t) is an nxn matrix A(t) is an nxn matrix I is a unit nxn matrix

PAGE 56

39 The complete solution is therefore x(t) = 4) (t)x(o) + p (2.74) The following initial and final conditions were used n^ = 5 k.w or 0.5 n (t=tf) = lOn^ = 5 32 P2(t=tf) = 10 0.5 X 10-^ Pi(tf ) = 0,30x10=^ ^^-"^^^ P20 = 0.18 X 10= The initial guesses for the ad joint variables , i.e. p^Q and P20 were taken from the analog computer results described earlier. The above case is for: neutron life time A = 10;^^3 . The numerical results of the first control problem may be summarized in the graphs. The graphs also indicate comparison with the analog computer results, which in general is very good. Two cases were run average neutron life time, y\ = 10~^secs and A = 10~^secs. A typical table, i.e., Tables, shows the results of the second case only. Comparison with analog computer results is the subject of Figure (6).

PAGE 57

40 OQ o CC a, O H a M a, S " »-.:tt nj CO a; 'o 'o •—I I— i X X c^ o 'o 'o •H a.+j c are

PAGE 58

41 The initial and final conditions used in Figure 7 are: Pj^(t£) = 0„0000020, P2(tf) = 0.0 n(o) = 0.5,c(o) = 32 and a = 10"^ Figure 8 shows power level versus time plot of run C-12A of Table 12. The initial conditions on the adjoint variables in this run are: p^Q = 0^54x10 and P20 = —5 0<,18xlO . Convergence in this case is obtained in the second trial.

PAGE 59

42 Quasilinearizotion Data Points 55 45 35 2 5 I 5 Data Points From Analog Computer Results .04 .06 TIME (SECS) Figure 6. Power level versus time ( A = 10 ) Corrparison of analog computer results with those obtained by quasilinearization (Nuclear Reactor System)

PAGE 60

43 50 40 20 -

PAGE 61

44 50 40 > 30 bJ o 020 10 r 1st Trio! By Quasilineorization »/ Last Trial By Quasilineorization^ / — 29}h Iteration By Grcdient Tec^lnlque On Analog Computer .4 .6 TIME (SECS) .8 1.0 Figure 8. Power level versus time( A = IQ-^) Comparison of analog computor results by gradient technique, with those obtained by quasilinearization (Nuclear Reactor System)

PAGE 62

45 2. Cost Function 2 = 1/2 f^[Pc2 + a(n-n^)^ 1 dt The two-point boundary value problem resulting by the application of the maximum principle is: (2.76) c = 6.4n 0.1 c (2.77) P^ = a(n-nd) 1000 p^P^ + 6.4 p^-6.4P2 (2.78) P2 =-0.1 P^ + 0.1 ?2 (2.79) Pc =-1000 nPi (2.80) with boundary conditions: n(o) = 5-0 c(o) = 64n^ (2.81) n(tf ) = 5.0a P2(tf) 0.0 The cost function of this problem may be written as the sum of the two cost functions as shown. J = 1/2 j [ ^c2 + a (n-nj)^j dt = Ji + ^2

PAGE 63

46 where J^ = 1/2 J P c^ (it o and may be thought of as the cost of the control action and ^f 9 J2 =1/2 j a (n-nd)^dt o is the cost of the deviation of the output of the dynamical system from the desired value. Simulation of Problem 2 This problem reduces to problem 1 for a = 0, A graph of the cost function for various values of a is shown in Figure 9 , Neutron density and reactivity function plots are shewn in Figures 10,11 for ae (0,20) with the final time tf fixed at 1 sec. It is clear from these figures that for increasing values of a , the rate of control rod withdrawal is also increased. However, in order to limit the power level to a certain desired value the control rod is inserted back into the reactor at a faster rate, thus shifting the reactivity function peak to the left.

PAGE 64

47 Table 5 shows the peak power attained, a and the cost functions Jj^, J2 and J. With appropriate scaling J2 is given by: 10"3 _ tf o -T— a j 0.1 (N-Nd)2dt J2 = 5 X lO""^ J2 where J2" = a J O.l(N-Njj) dt From the graph of Ji versus O (Figure 12) it appears that the cost function increases with increasing a . This appears reasonable as increased a would decrease the error between n and n^ and hence more control would be needed. The total cost function increases rapidly with an increase in a . as is apparent from the graph of the total cost J versus N, the peak power attained (Figure 13). and also from graph of J versus a , This result is again quite logical as one would expect to pay a higher cost in crder

PAGE 65

o Figure 9, PEAK POWER DESI RED , N^ i"otal cost function versus peak power desired (control problem 2)

PAGE 66

49 0.4 0,6 TIME (SEC) 50

PAGE 67

50 0.4 0,6 TIME (SECS) Figure H. Power level N and reactivity ^g versus time for control problem 2 for various values of ^ and n , — 5n (Nuclear Reactor System)

PAGE 68

51 TABLE 4 COST FUNCTIONS AS A FUNCTION OF WEIGHTING AND POWER LE>/EL N^ FOR PROBLEM 2 2Ji a =z 5.0 J=Jl+J2 2no 0,456x10"^ 5no 0.696x10"^ 6no 1.323x10"^ lOrio 1.914x10"^ -5 0.0149x10 0.527x10"^ 0.823x10"^ 2.88x10"^ a = 10 0.243x10"^ 0.875x10"^ 1.484x10""^ 3.84x10" ,-5 2nQ 0.456x10 ^ 5no 0.696x10"^ .-5 0.0297x10 ,-5 0.258x10"^ IOpq 2.058x10"^ 4.95x10 ^ a = 20 6,0x10"^ 2no

PAGE 69

52 Figure 12. Cost function ? movement versur> piobi em 2) iSociated with contio]. rod M.10 n., (control

PAGE 70

53 Figure 13. Total cost function versus ^. , for control problerr. 2 end n^ =: lOn^^

PAGE 71

54 3. Cost Function 3 f _ 9 . ^ / ^2 J = 1/2 J Pc^ + a (n-n^) o dt p Q = a + A(n-nci) The two-point boundary value problem resulting by the application of the maximurr. principle is: n = 1000 an + 1000 A(n-n^) 6.4n + 0.1c (2.82) c = 6.4n 0.1c (2.83) (2.84) A = 0. a = P^ = (a2 + a )(n-n^) aA 1000 P^ (a+2An-Anj) + 6.4?^ 6.4P2 P2 = O.lPi + 0.1 P2 (2.85) (2.86) P^ = A(n-nj)^ 1000 np^(n-n^) a(n-nd) (2.87) P4 = -a A(n-nd) lOOOnpj^ (2.88) Pc a + A(n-nd) (2.89) The boundary conditions are: n(o) = 5.0 P2(tf) = c(o) = 64nQ n(tf) = 5.0 a (2.90) P3(o) = P4(o) z: P3(tf) = P4(tf) =

PAGE 72

55 Simulation of Problem 3 The equations used for simulation are given in Appendix A, and the procedure used is similar to that followed in Problem 1. The circuit diagram is shown in Appendix A. The trial and error procedure is adopted in the sense that the constants "a" and "A" are varied to obtain the desired power level response. Figure 14 shows the variation of "A" with "a" to achieve different power levels ranging from Sn^ to Sn^. Corresponding to these power levels, the desired control is obtained. The initial conditions of P^ and P2 and the weighting factor a are now varied to satisfy the end conditions on the adjoint variables, namely P2(tf) = Paltf) = P4(tf) = 0. It is observed in this investigation and justified by solution of equation: N = AN(N-N(j) 0.64N + C that the desired power level is reached only for an infinitely large value of A, Figure (15). The power level response also suggests a strong dependence on the desired level. These observations therefore suggest that the form chosen for this control will not yield a desirable closed loop control. The form that accomplishes this is given in Problem 4.

PAGE 73

Figure 14. V/ariation of parameters in the control law for s.o.c. problem (i.e., problem 3) and for various values of V.^

PAGE 74

57 9n, 8n^ 26n, 4n, VALUE OF Figure 15. Peak power level reached versus A (control problem 3)

PAGE 75

58 4„ Cost Function ^ J = 1/2 l^ [Pc^ '' ^ ("-"d)^J dt t£ p^ = A(n-n^) + B J (n-nd)dt o Let X = J (n-n^^)dt Hence x = (n-n^). and p^ = A(n-nc|) + Bx o The two-point boundary value problem resulting by the application of the maximum principle is: n = 1000 A n(n-nd) + 1000 Bxn 6.4n +0.1c (2.91) c = 6.4n 0.1c (2.92) X = n-n^j (2.93) A = . B = (2,94) P2 = (A^+a )(n-n(j) ABx 1000 P^ F 2An=Anj+BxJ + 6.4?^ 6.4P2 P3 (2.95) P2 = -O.IP^ _ O.IP2 (2.96) P3 = -B^x 1000 BnPi AB(n-nj) (2.97) P4 = -A(n-n^)^ lOOOnP^(n-n^) Bx(n-n^) (2.98) P5 =-Bx2 1000 xnP^ A x (n-n^) (2,99) p^ = A(n=nd) + Bx (2.100)

PAGE 76

59 (2.101) The boundary conditions are: n(o) = no c(o) = Co x(o) = n<^(tf ) = n^d = ano P2(t£) = P3(tf) = P4(tf) = Psdf) = P4(o) = Pp^(o) = Simulation of Probletr 4 It was seen in Problem 3 that the power level response was highly dependent on the setting of the desired value of the neutron density. In order to remove this restriction a control function of the type shown above, namely a combination of proportional and integral control, was tried. A simple transfer function analysis indicates

PAGE 77

60 that such a system would be stable and the response would be dependent on the coefficients A and B. The equations used for simulation are given on page 63 and the procedure used is similar to that followed in Problem 3. Circuit diagram is shown in Appendix A. Figure 16 shows the power level and the reactivity function response for values of n^j ranging fron 2nQ to lOn^ with the final time fixed at one sec. Note in this case that the power level rise is much faster, even though it settles down to the desired value in fixed time tf . Corresponding rod withdrawal is also much faster in the initial stages and then gradual before it reaches a constant value, to maintain the above level. Figure 17 shows variation of cost function versus peak power attained. A graph is also plotted showing behavior of constants A and B with respect to peak power attained (Figure 18). It is clear from this graph that there is no value of optirrum A and B which is independent of the desired power level. This would imply that the control parameters A and B should not be fixed but should be varied as a function of the desired final state. Hence one would conclude that some sort of adaptive control may be required. A control system is said to be adaptive if it is capable of changing its control object in such a way as to operate at all times in an optimal or nearly optim.al fashion.

PAGE 78

61 Figures 19 and 20 indicate the time response of the auxiliary functions. Since these functions are highly unstable, they are very sensitive to amplifier drift. Because of this reason it is very difficult to satisfy and maintain the end conditions, namely, P2(t£) = Po(tf ) = P4(tf ) = P5(t£) =: 0, and n(tf) = an^. TABLE 5 COST FUNCTION AS A FUNCTION OF POWER LEVEL N^ AND CONSTAhTTS A AND B FOR PROBLEM 4 (NUCLEAR REACTOR SYSTEM) ^1 lOno

PAGE 79

62 Nd=IOno A .6 TIME (Sec.) Figure 16. Pov/er level ^J and reactivity c^ versus time for specific optimal control problem 4 (Nuclear Reactor System)

PAGE 80

63 1.3 L2 10 0.9 I 0.8 0.7 0.6 0.5 0.4 0.3 0.2 1 I 23456789 10 I Figure 17. Cost function versus peak power attained

PAGE 81

64 B 130 120 no 100 90 80 7060 50 40 30 20 10 A 4.0 3.5 3.0 2.5 2.0 I.O oL J I I 1 — J_l I I ^L 23456 7 89 10 n^ DESIRED Nd Figure 18. Nj versus ^ and B

PAGE 82

65 .10 .OS ^ .10 r Nd= lOno Nd= 6no Figure 19. Auxiliary functions, P2 , P3, F4 and PC) versus time for control problem 4; N(j lOn^j and tn^ (Muclear Reactor Systen-i)

PAGE 83

65 Nh = 5n, "

PAGE 84

67 Discussion of Results Results obtained in the four problems attempted in this chapter are analyzed and a comparison is given of the two basically different control systems, i.e., 1) an open loop control system and 2) closed loop control or specific optimal control system. The first two problems analysed represent open loop control. Problem 1 is simulated in Section B.l and the operating time is set equal to one second. The performance criterion used here could be interpreted as an optimization of control and movement, limiting erratic flux changes (10) . Another useful performance index is the subject of problem 2 wherein the integral of the error squared of the reactor power and its desired value is minimized along with the integral of the squared of the reactivity function. In this problem the control function required to achieve the desired state for example n(t£) = an^ is dependent on the value of a. The higher the value of a, the faster the control action and the faster the power rise. This variation is exhibited in Figure 10. Also it is observed that the value of the cost function increases with increasing value of a (Figure 12). This is to be expected since the higher value of a would imply desirability of

PAGE 85

68 driving the error to zero at a faster rate. Naturally more control is needed to accomplish this. Problems 3 and 4 represent closed loop or specific optimal control. The specific forms chosen for the controller are : p^ = a + A(n=nj) and p^ = A(n-nd) + B j (n-nj)dt Th^ performance index used is the same as in Problem 2. The results of problem 3, point out a strong dependence of the peak power attained to the desired value of the power level. This result led to the investigation of Problem 4. In this problem, one observes that the cost function increases with increase in the desired peak power (Figure 17), a result which is intuitively obvious. In problem 4, the control laws are determined to achieve four different power levels and it is found in each case that the closed loop or specific optimal control suffers in performance in the sense that the cost function is higher than in the open loop case (Figure 21). If the values of the control parameters A and B are fixed as those obtained for the case N^ = lOnQ and the reactor is brought to various power levels, then the dotted curve in Figure 22

PAGE 86

69 2.0

PAGE 87

70 4.0 3.0 2.0 1.0 A = 2.13 X 10 [C.L.C.] 2n 4n, 6n, o -»"o ""0 DESIRED N^ [O.L.C] 8n, ion. Figure 22. Total cost function versus peak power attained for open loop control and cost variation for fixed and optimally varying values of control paran eters for closed loop control (Nuclear Reactor System)

PAGE 88

71 is obtained for the total cost function. Note here that th< cost function is higher in every case. These results compare favorably with those obtained by Eisenberg and Sage (21) power level trajectories and the corresponding time variation of the control functions are shown in Figures 23 and 24. Figure 25 depicts the variation of the cost function versus the peak power desired for varying values of Figure 26 shows the variation of the cost function with variation in system parameters (neutron lifetime). It is noted here that the closed^loop control is better than the open loop control for large variation in lifetime about its normal value. It may be pointed out at this stage that the coefficients in the control law are highly dependent upon the final time t^, the system parameters, and the initial state of the system. The performance of the system, thus is a function of the form chosen for the control law. The form of the control law is very instrumental in determining how well the performance of the s.o.c. compares to that of the true optimal. Figure 18 shows the variation of the parameters A and B in the control law, as a function of the peak power. Since this variation is considerable, adaptive control is visualized as a means of further improving the performance. If adaptive control is used, the example at

PAGE 89

72 0.4 0.5 0.6 TIME (SECS) Figure 23, Open loop and closed loop cptirnal trajectories N(t) versus time (Nuclear Reactor System)

PAGE 90

73 Figure 24. Open loop and closed loop optimal control function to obtain ten tirres the initial power level in a Nuclear Reactor System

PAGE 91

74 1.0 0.9 0.8 in 0.6 o 0.5 0.4 0.3 0.2 -J ^ \ ^ \ I \ I I 2 4 6 8 10 DESIRED POWER LEVEL N, Figure 25. Cost function versus peak power desired

PAGE 92

75 6.0

PAGE 93

76 hand could be that of the so called performance criterion sensing or extremum adaptive system. An overall performance index is measured and an attempt is made to find a control law (by trial and error) which optimizes this performance index. The learning states here are the parameters describing the control law, and the learning process consists of the trial and error adjustment of these parameters. If the learning states can be changed rapidly and accurately the adaptive controller may be able to handle nonlinear control systems the learning states will represent the best instantaneous linear approximation to the non-linear system. This approach is currently under study (22) . The gradient technique employed yields good results. The accuracy of the results can be considerably improved at the expense of more analog and logic elements and adjusting values of K"s. The convergence of the interations is very sensitive to the values of the K's. Larger values of K's result in excessive overshoot about the desired trajectory. Smaller values result in an increased number of cycles to converge. The Quasilinearization technique is next employed to compute the results already obtained in Problem 1. Figures 6 and 7 indicate close agreement between the results obtained by the two techniques.

PAGE 94

CHAPTER III OPTIMAL CONTROL OF A NUCLEAR ROCKET ENGINE Introduction and Model Description The nuclear rocket engine to be considered is similar to the solid-core, 25,000 megawatt, million pound thrust model developed by Smith and Stenning (23, 24,25,26). The system of equations describing such a model are based on a lumped parameter analysis. Such an analysis does not describe the system corrpletely and furthermore the analysis is affected by a lack of agree= ment among investigators concerning the temperature reactivity effects. Stenning and Smith (27 ) contend that temperature reactivity is directly proportional to the square root of the core exit stagnation temperature, while Mohler and Perry ( 2 ) hold that temperature reactivity is directly proportional to this temperature. Weaver ( 12) observes, that both contentions could give fairly accurate results simply by choosing appropriate reactivity coeffi^ cients. Recent rocket engine tests (28 ) however indicate that temperature reactivity may be proportional to square of this temperature. All this adds up to the fact that the nuclear rocket engines possess characteristics which need exhaustive analysis. Thus, since the system is not well 77

PAGE 95

78 understood, a closed=loop control is desired. The maximum principle however, leads to an open«=>loop control. The closed=loop control is a complex function of the state variables and due to the analog computer limitations, the treatment is restricted to open loop control only. Basic components of the engine are depicted in Figure 27 . Liquid hydrogen, stored in a low pressure tank, is pum.ped to high pressure, and is used to regenera = tively cool the nozzle and reflectoro The propellant gas is further heated as it passes through the reactor core and is exhusted through the nozzle„ A portion of the hot gas is bled off from the main stream, and is used to power the turbo pum.p , a feature giving rise to the name, "Bleed turbine system." The bleed flow, and consequently the sys= tern pressure are controlled by the valve in front of the turbine. The following equations describe the Smith-Stenning nuclear rocket engine in terms of the time dependent nor= malized quantities. An account of the normalization P^o ^ cedure is given in references (23) and (29) and therefore will not be described here. i M5 B N' = P'N» N« + C» (3a) 1 r A C» = N« C« (3„2)

PAGE 96

79 en

PAGE 97

80 N» PHT)^^'^ (3.3) p« = P'V (T« ) 1/2 .lElli P ' . '' ' (T )l/2 (3^4) >= Pc " a; (T«)^^2 'a' ~° (3^^) where N» = 77. C« = ~ Nd Cd ^ --fe and Equations (3,1) and (3,2) are the reactor kinetics equations. Equation (3..3) gives the time rate of change of maximum, core surface temperature, and equation (3,4) gives the time rate of change of the core inlet stagnation pressure. Reactivity is given by equation (3^5) Table 6 contains important constants for the nuclear rocket engine. Note that the delayed neutron fraction is smaller than the value typically associated with a uranium235 reactor. This reduced value was used by Smith and Stenning to account for diffusion of some precursors into coolant channels and their subsequent expulsion from the system before emission of delayed neutrons.

PAGE 98

81 TABLE 6 rcx:ket engine constants Design reactor power 25,000 MW Design core surface temperature 5,000°R Design core inlet stagnation pressure 1,500 psi Thermal time constant 0.68 sec Pressure time constant 2.50 sec Control poison time constant 2.00 sec ^^alve time constant 0.20 sec Decay constant (one group) 0.10 sec"^ Temperature reactivity coefficient -10~"" 0R~^ Hydrogen reactivity coefficient 0,3 °R/psi Neutron lifetime 10~ sec Delayed neutron fraction 0.005 Thrust, F 1,500.000 lb Specific impusle Isp 640 sec

PAGE 99

82 The specifications in Table 6 significantly exceed those of nuclear engines presently under development and will probably be first achieved through clustering of small engines. However, since the above system has received considerable attention in the open literature, it makes it an attractive system, for an optimal control study^ It may be pointed out here that the reactivity equation (3.5) is modified, particularly with respect to the temperature coefficient of reactivity. The (T")^^^ term is replaced by T\ as in the Mohler=Perry ( 2 ) formulationo In view of this modification, equation (3.5) is re=written as : The above system of equations with the constants of Table 6 reduce to: N' 50 ( P«N«-N« + C») (3.7) C» = 0.1 (N« CO (3.8) f' = 1.471 (N» P«(T«)^'^^) ^ (3.9) (T»)V2^ (3^10) 2 P' =: 0.4 (P« V» (T«)^/2 IP;)^ ,) p, =p,^ T« + 18 ^] (3.11) The system is in an idle condition for T^O, the idle or the initial condition being -.T^ = 1000°R and P^ = 275 psi. The value setting is fixed at 0.915. The full

PAGE 100

83 power condition is T = 5000°R and P = 1500 psi. It is desired to optimize system performance such as to mini= mize the following integral. tr J = 1/2 i p2(t)dt o with T(tf) = T^ and ip, < p^^^ This performance index can be interpreted as tending to limit the maximum reactor period as is easily seen by considering equation (3.7). Since the effect of delayed neutrons is small in the interval of operation, N«/N« = 50 P « . This indeed is the inverse reactor period By using the statement of the maximum principle as in Chapter II, the following two point boundary value problem is obtained. N» = 50( P «M« N» + C") C» = 0.1(N» CM T« = 1.47 (N» P»(T«)^'^^) (P^^2 P« = 0.4(P«V«(T")^/2 T^^flT^) P» = p»^ T« + 18 !^ Pi =-50 P'P^ + 50 P^ Oa P2 1.46 P3 P2 =-50 P^ + 0.1 P2 P3 = 50 N'Pi + 900 ^'^^ ' 2' ^ ^•'^^^ PMT")" aP4

PAGE 101

84 P4

PAGE 102

85 Solution Procedure : Since the useful life of a nuclear rocket per shot is approximately 2 to 15 minutes, the nuclear second stage must begin to supply full thrust almost immediately, following a short start-up period. Thus the rocket reactor must be brought up to full power, of hundreds of megawatts, in a very short time. In light of this fact the optirrization problem was solved for various ranges of desired final time, 4 $ t^ ^ 11 sec, desired final state 0.4 ^T. ^1.0 and variations in the a system parameters, 2 t^ < t ^10 t and 5 T^p ^ T^f; 10 t^^ ^^ariable Start-up Tjme : 4 $ tf $: 11 sees. The final time t^, in which the reactor J.s brought to full power is varied between 4 sec and 11 sec. Three different trajectories, as shown in Figure 28 are used in the optimization procedure. The corresponding control poison reactivity and the pressure responses are shown in Figures 28 and 29. These are obtained by forcing the system to follow the desired trajectory in an optimum manner. It is noted from the figure, that shorter the temperature response time, the faster is the control rod withdrawal and the faster is the pressure rise in the system. The rod has to be re-inserted to maintain the final temperature at the desired level. This re-insertion is very gradual for longer temperature response times.

PAGE 103

86 ''ariable desired temperature: 0.4 ^ j. .; 1.0 The final time tf, is arbitrarily fixed to 10 sec and the system is driven to three different desired average core temperatures from the idle condition of 0.2 (1000°P). In each case it is desired to drive the system to some final state in 10 seconds such that 1/2 ;' 2 P dt is minimum o The optimal temperature response and the corresponding control poison reactivity response for T^ = 0.4 (2000°F) is given in Figure 30. Figure 31 shows the optimum pressure response for this case. Figures 32 and 33 show the optimum temperature, control poison reactivity, neutron density, and pressure responses for T^ = 0,7 (3500OF) and Td = 1.0 (SOOOOF). With the optimal trajectories determined and the cost function J computed for these trajectories, it is now possible to make a plot of the cost function J versus desired average core temperature (Figure 34). The cost to drive the system to desired core temperatures inr creases slowly at first and rapidly after Td = 0.4.

PAGE 104

'to

PAGE 105

88 6 8 10 12 14 TIME (sec) Figure 29. Pressure-time plot, tf = 4 sees, 6 sees and 11 sees (Nuclear Rocket Engine System)

PAGE 106

TIME (sec) Figure 30. Optimal control poison reactivity and temperature-time plots, Tjj = 0.4 (Nuclear Rocket Engine System) 10

PAGE 107

90 N < 2 4 6 TIME (sec) Figure 31. Optirral pressure = time plot Rocket Engine Systerr) Tj = C,4 (Nuclear

PAGE 108

Td=i.o -9 > o < UJ or CO o Q. -15 1.0 ir < 0.8 en Id =1.0 •-0.6 Q UJ M 0.4 2 1 = 0.7 4 6 TIME (sec) 10 Figure 32. Optimal control poison reactivity and temperature-time plots Rocket Engine System) 0.7 and 1.0 (Nuclear

PAGE 109

92 Tw -1.0 TIME (sec) Figure 33. Cptirral pressure and neutron density versus time, Jd = 0.7 and 1.0 (Nuclear Rocket Engine Systerr)

PAGE 110

93 .2 .4 .6 .8 NORMALIZED TEMPERATURE , 1.0 Figure 34. Cost function, J versus desired average core terrperature (Nuclear Rocket Engine System)

PAGE 111

94 3 ) Influence of Parameter v^ariation on System Response The following four cases were examined with respect to the normal temperature and pressure time constants, the normal values being, y. = 0.68 sec and T ^n ~ ^'^ sec. The reactivity input used was optimal one obtained for the case of 7^ = 1.0. a) Tp = 2 Jpr^ , :^ = t„ b) yp = lO-Ypn .Tt = tn c) J p = Tpn . Tt = 5 Jtn d) Tp = Tpn . U = 10 tn The results of the variation of the pressure time constant from the normal value is illustrated in Figure 35. As the pressure time constant is increased, the temperature response becomes quite sluggish. The corresponding pressure response is shown in Figure 36. If the pressure time constant is held constant and the temperature time constant is increased, the system response is not influenced as greatly as in the previous case (Figure 37). Humphries (30) has reported similar results on model reference adaptive control of a nuclear rocket engine. In light of the results obtained here and

PAGE 112

95 Humphries' results, it appears that the next appropriate step would be to drive the system with a control whose specific form is known. The controller gains will have to be varied, for the system to optimally adapt to the desired response. If this is not satisfactory an optimal adaptive controller will have to be used. Such a scheme is currently under development ( 22) .

PAGE 113

96 Figure 35, Control poison reactivity versus time and average core terrperature versus time for three cases of parameter variations (Nuclear Rocket Engine System)

PAGE 114

97 0.6 0.5 < nn en Z) CO 000.4 LlI (T QP0.3 aTpn < §0,2 '0 Td 4 6 TIME (sec) 10 Figure 36. Pressure versus time response for three cases of parameter variations, (Nuclear Rocket Engine Systerr)

PAGE 115

98 Figure 37, Pressure versus time and average core temperature versus time for three cases of parameter variations (Nuclear Rocket Engine System)

PAGE 116

CHAPTER IV SUMMARY AND CONCLUSIONS Summary and Results : Optimal open-loop control and the specific optimal control of nuclear reactor processes in a nuclear reactor system and open-loop control in a nuclear rocket engine system, has been analysed in this study. Optimal control theory is applied to the problems and the analog computer is used for simulation and in some cases optimization by the gradient technique is used. In certain examples, quasilinearization technique is used for the solution of the TPBVP. In all of these, the convergence is obtained in no more than four iterations. In the examples considered here, for nuclear reactor system, the closed-loop control suffers in performance as compared to the open-loop control in that the cost function is higher. This is in general agreement with physical insight and previous results (21), The coefficients A and B in the control law are highly dependent upon the final time t£. the system parameters, and the initial and desired final state. Thus, it is desirable to vary the closed-loop controller parameters as a function of neutron density and allowed start-up time, as well as system parameters, if performance near optimum is to be achieved for varying start-up time and final neutron density. This 99

PAGE 117

100 suggests the desirability of an adaptive type control, as a means of further improving the performance for such applications. Such a scheme is currently under investigation by Sage and Ellis (22, 3i). The open-loop control problem for the nuclear rocket engine system was solved in three different parts: a) variation in the start-up time, b) variation in the desired average core temperature, c) influence of parameter variation on system response. The results obtained appear quite reasonable. For example, in the first part it is observed that shorter the start-up time, the faster the control rod withdrawal and the faster the pressure rise. Moreover, in order to main= tain the system at the desired final state, the control rod has to be re-inserted and this re-insertion is very gradual for longer start-up times as opposed to shorter start-up times. In the second part it is seen that the cost function, J, increases, if it is desired to drive the system to higher and higher power levels. In the third part, it is observed that if the pressure and the temperature time constants were increased from their normal values , the temperature response became somewhat sluggish. In the light of this observation and Humphries (29) results, it appears that the next appropriate step would be to drive the system with a control whose specific

PAGE 118

101 form is known. The controller gains will have to be varied, for the system to optimally adapt to the desired response. Then try optimum adaptive or adaptive control. Since it is possible to reduce the start-up time considerably using optimal control theory and still maintain smooth transition to full power, this can -result in large amounts of propellant savings and hence reduction in payload. Does computation of optimal control and its inclusion in driving the system involve additional equipment to render it impractical? This is not considered to be the case as it is anticipated that computational equipment will be available in the total payload. Also in this study, the control valve opening is fixed at 0.915 as opposed to the maximum of 1.33 to 3.0 considered by Smith and Stenning (23) and 1.10 considered by Humphries (29). This would in turn lower the turbine size because of lower flows and thus allow additional reduction in weight. The results obtained here are therefore quite encouraging. However, this control concept must be evaluated within the framework of the total body of knowledge pertaining to the development of nuclear rocketry. How these short start-up times affect the attitude and guidance control problems and various other questions pertaining to the equipment and their costs cannot be answered here. The results are sufficiently encouraging to prompt further investigation.

PAGE 119

102 Suggestions for Further Study A logical sequel to this work would involve the study of optimal adaptive control of the two systems in order to improve the system performance. This is currently undertaken by Sage and Ellis. Moreover, the control emphasis should be shifted from the control poison to the hydrogen propellant. An optimal control of more complicated reactor nodels and other constraints should be studied, since the optimal process is only as good as the constraints and the model used. In particular the nuclear reactor system with spatj.al dependence should be given consideration and the reactor system should be modified to include effect of reflector.

PAGE 120

APPENDICES

PAGE 121

APPENDIX A ANALOG SIMULATION Circuit Symbolism The analog computer wiring diagrams for the simulation of the four control problems is shown in Figures 41 and 42. . The explanation of the circuit diagrams is given in Figures 38 and 39. The special circuits used are shown in Figure 40 . The values of the potentiometer settings (K) are contained in tables rather than in the computer diagram. Only an identification number is used in the computer diagram. Diodes are used in the square root circuits and in limiting circuits. Input and feedback resistances for summing amplifiers and input resistances for integrating amplifiers are expressed in megohms. Values of either 0.1 or 1.0 megohms are used. Feedback capacitances of integrating amplifiers are expressed in microfarads, with values of 0.1 land 1.0 microfarads being used here. An identification number is assigned to each amplifier in the computer circuit diagram. Special circuits for square-root and division are formed by suitable combinations of quarter-square multipliers and summing amplifiers (Figure 40 ). With the symbolism thus 104

PAGE 122

105 Figure 38. Symbols used in analog computer diagrams

PAGE 123

106 -POTENTIOMETER ^A^.©2 — J\AA. R„ r Rf — hAA-A, 1 SUhiWU'JG Ar.lPLiFlSR AAA.^— Rg _ ^AA.\. r .t n c / E liec't IMTEGRATUJG AMPLIFIER -0, t-e O.i M QUARTER -CJQUARE' f.iULTlFl.iER iir r\L DIODE

PAGE 124

107 Figure 39, Logic symbols used in corrputer diagrams

PAGE 125

lOB A-B A-B AND GATE A D ^-. — A 48 — A+B OR GATE B A C

PAGE 126

109 — vAAA/ .1 M -oiI — -w' \j \/ \.« — — I CIRCUIT FOR SQUARE ROOT s = 10 yio .IM

PAGE 127

no defined, it is now possible to describe the various siniu= lations , performed in the course of this investigation. Problem Scaling (Nuclear Reactor Systefm) Amplitude and/or time scaling is usually a necessary part of analog computer work to contain amplitudes of all variables to within + 100 volts and time to within measurable quantities on the equipment used. The scaling used for the four problems is given in Table "^ . The bars above the letters denote scaled quantities. With this scaling the TPBVP equations for the four control problems become: Problem 1 N =-100 N^P^ 6.4 N + IOC t = 0,064N O.IC Pi = IOC NP^^ + 6,4P^-6.4P2 P2 = -0.1 P^ + 0.1 P2 Pc = 0.1 N P^ with boundary conditions: N(o) = 0.05 ; C(o) = 0.032 Pj^(tf) * and 'P2(tf ) Problem 2 N = lOON^P^ 6.4 N + 10 C C = 0.064N 0,1 C

PAGE 128

Ill P^ = lOON P^^ + 6,4 P^ 6.4 P2 + a (N-Nj) P2 = -0.1 P^ + 0.1 P2 Pc =-0.1 N Pi J = 1/2 J p^2 dt + 5x10°^ a J O.l(N-Nd) dt o o J = J^ + 5x10^"^ a J2' with boundary conditions same as in problem 1. Problem 3 : N=0.1 aN+ AN (N=Nj) 0.64 N + C t = 0.0064 N 0.01 C Pi = -0.01( a + A^)(N-Nj) 0,1 a P^ 2A (NP^) + A N^ Pi + 0.64 P^ 0.64 P2 0.001 a A P2 = -0.01 P^ + 0.01 P2 P3 = -0.01 A (N-N^)^Wi (N-Nj) 0.001 a (N-Nj) P4 = -10""^ a 0.1 N P^ Pc = 10"^ a + 0.01 A (N-Nj), A is negative Problem 4 ; Fi = AN (N-Fid) +0.1 BxN 0.64 N + C C = 0.0064 N 0.01 C X = 0.1 (N-Nd) P^ = -0.01 (A^+a )(N-N^) = 0.2 A (lOMPj^) + ANdPi 0.1 BxP^ + 0.64pj 0.64P2 O.lP^ 10'^ ABx

PAGE 129

112 ?2 = -0.01 P^ + 0.01 P2 P3 =-0.0l(0.01B^)x 0,1 BNP-^ AB lO'^ (N-N^) P4 = -0.01 A(N-N^)2 NPx(N=N^) -lO"^ B x(N-Nd) Pp^ = =0.01 A(N-Nd)^ -0.01 (O.lB)x^ xNPj^ P^ = 0.01 A(N-Nd) + 0.01 (OaB)x A and B are both negative The boundary conditions on the state variables are: N(o) = 0.05 , C(o) = 0.032 x(o) = 0; P3 (t^) = P4(tf) = P5(tf) = P4(o) = P4(o) = r^ ^ f^f 2 J = 1/2 J p^2 dt + 5x10°^ O J 0.1(N-N ) dt o o J = J, + 5x10°^ a Jo° h = 1/2 1 pM where ^^ 'o and J2' = J 0.1(N-N.) dt With the system equations scaled, and the initial conditions determined, it is now possible to proceed with the analog computer simulation. The analog computer circuits for simulating the above problems are represented

PAGE 130

113 in Figures 38 and 39. Potentiometer settings are listed in Tables 8, 9, 10 and 11. It should be mentioned here that great care and ingenuity is needed to arrange the analog circuit, so that the noise is not amplified at any particular stage. This is particularly brought about by the fact that in the square root circuit, the output is multiplied by a factor of ten and in the division circuit the output is multiplied by 100 . TABLE 7 SCALE FACTORS FOR THE FOUR COrJTRCL PROBLEMS Problem

PAGE 131

114 E

PAGE 132

115 TABLE 8 POT E NTT IOMETER SETTINGS FOR CONTROL PROBLEMS 1 AND 2 (NUCLEAR REACTOR SYSTEM) PotentioParameter Description +100v Problem 1 Problem 2 meter 1 2 3 Initial Condition on N +100 4 Initial Condition on C +100 0.032 5 6 7 8 9 Weighting factor for square of the neutron density error term i.e, 10 0,64

PAGE 133

116 TABLE 8 (CONTINUED) PotentioParameter Description +100v Problem 1 Problem 2 meter 19 Kj2 varies varies from 0.01 from 0.01 to 0.04 to 0.04 20 Koi varies varies ^^ from 0.01 from 0.01 to 0.04 to Oo04 varies varies from 0.01 from 0.01 to 0.04 to 0.04 21 K22

PAGE 134

117 TABLE 9 POTENTIOMETER SETTINGS FOR CONTROL PROBLEM 2 (NUCLEAR REACTOR SYSTEM)

PAGE 135

lis

PAGE 136

119 o c o c iT) O -4 o o o o c c c if) n n CNI

PAGE 137

120

PAGE 138

121 in

PAGE 139

122 lO

PAGE 140

123 1 in

PAGE 141

124 +J c 0) CM CO -I O 00 o ^ m o in ^ CO o 00 n o +j o -fJ c a> iT) (A II o lO >o CM CO O O CO 00 ao o (^ rH [^ -t-> o 0) o o o rro CM in in in o in in -H in c •H O P C PO 0) ^ CM O ^ in CM r-H O CM -H I I p a e >H

PAGE 142

125 p o +J c (A -p -o "^ 2 +J o +J c -o" C XI c •rl O t-> c -PO (/) CUD CX) ^ 00 CM O CN n n \0 o o ^ "^ ^ CM CM CM O

PAGE 143

126 S imulation by Gradient Technique The algorithm used in this technique is given by equation (2.18), The problem is described on page 12 and the analog circuit is shown in Figure 41. it must be borne in mind that the integrators must all be started together; thus the operate (op) inputs are all joined on the logic patchboard. Each time the op input receives a logic one, a solution n(t) is generated. By systematically changing P^, a search for P^^*, which produces the desired optimal control u* (reactivity in this case) and trajectory n* can be made. The computer control diagram for iterative computation is shown in Figure 43 whereas the updating circuit is in Figure ^'^ . Before the start of the computation the Rep-op is off and the integrators in the analog circuit are in the I.C. or reset mode. The initial conditions P^"'" and P2"'~ a^e now set. The first run begins when Rep-op pushbutton is depressed. At t = t^ , the comparator output triggers the memory units and transfers the first boundary point to their output (i.e to the output of the second amplifier in the track-and-hold system of the updating circuit (Figure -^4 ). since the outputs are held until t = tf (or the interval set by the pulser), there is ample time for read-out (on the brush recorder or oscilloscope), before

PAGE 144

127 the next tun begins. This is a great advantage of the analog=digital computer because it can show the human operator as to actually what is happening via an on-line picture on the oscilloscope. The iterations can be watched step by step and the source of trouble can be predicted with ease. Suppose for example, it is desired to satisfy the following end conditions: Mlt^) =: 10 n^ and N(tf ) = 0. Then the p^"*"^ for which the above conditions are satisfied joS for various iteration step sizes K^^, K12. ^21, and K22 is given in "^able 12, It is clear from the table that large values of K result In large overshoots in the desired trajectory„ Sm.aller values however result in increased number of iterations to converge^

PAGE 145

128 TABLE 12 P^°S AS A FUNCTION OF ITERATION STEP FOR PROBLEM 1 JO (NUCLEAR REACTOR SYSTEM) Run Iteration step, K^^ = K22 = ^21 = K22 = 0.01 ^10 ^20 Pfo ^00 ^ °^ Comments I Iterations 1 0.2169 0.1762 0.42 0.1762 19 No overshoot 2 0.19 0.1762 0.50 0.1762 37 No overshoot 3 0.25 0„1762 0.50 0.1762 35 No overshoot 4 0.30 0.1762 0.53 0,1762 48 No overshoot 5 0.35 0,1762, 0.1762 131 Problem terminated 6 0.40 0.1762 0.1762 very large Problem term= inated 7 0.45 0.1762 0.1762 very large Problem term„ ^ ors inated 8. 0.52 0.1762 0.1762 very large Problem term= ^ ^ inated 9. 0.55 0.1762 0.1762 very large Very close to desired trajectory Iteration step, K^^ = K^2 = ^2^ = K22 = 0.02 11 0.19 0.1762 0.46 0.1762 17 No overshoot 12 0.21 0.1762 0,505 0.1762 15 No overshoot 13 0.25 0.1762 0.51 0.1762 13 No overfehoot 14 0.30 0.1762 0.54 0.1762 19 No overshoot 15 0.35 0.1762 0.48 0,1762 : 75 Problem terminated

PAGE 146

129 TABLE 12 (ContinuecJ) Run Iteration step, K^^j^ = K22 = K21 = K22 = 0.02 # — P{n ^on P(o Pqo ^ °f Comments t_ ^^ ^" ^" Iterations 16 C.40 C.1762 0.393 0.1762 very large Problem terminated 17 0.45 0.1762 0.444 0.1762 very large Problem terminateci 18 0.50 0.1762 0.50 0.1762 very large Problem terminateci 19 0.55 0.1762 0.55 0.1762 very large Problem termi= nated

PAGE 147

130 TABLE 12 (Continued) Run # Iteration step, K^^ = Ki2 = K21 = K22 = 0.01 ^10 '^20 ^10 ^20 ^ °^ Comments Al

PAGE 148

131 TABLE 12 (Continued) Run Iteration step, K^ = K2^2 ^^21 = ^^22 ^.03 # pf^ pin p!^ p!_ # of Comments 10 20 '^10 '^20 Iterations A21

PAGE 149

132 TABLE 12 (Continued) Run

PAGE 150

133 TABLE 12 (Continued) Run # Iteration step, K^-^ = Kx2 = K21 = K22 = 0-02 ^lo "20 ^^0 4o iJ.°^,o„3 """"-^"^^ Does not iterate Does not iterate 12% overshoot No overshoot excessive overshoot no overshoot A% overshoot No overshoot No overshoot Does not iterate Does not iterate Does not iterate Does not iterate Does not iterate Does not iterate Does not iterate Does not iterate No overshoot (not plotted) Compare with B24. C7

PAGE 151

134 Figure 43. Computer control diagram

PAGE 152

13t n-0 C-IC -5-100 CO MP G22 "[> fr •MOO PULSHR .T *l G-IG iC\To OP OPJCriC lir;'. -u -4-J r~ p.?"'^ ;^• E OP TiiuG T = p/Q= ISec 65 t C-D-E n{T) ^•n(T) n(t) _.. J : -Qh) A-!0 cor. IP Gl -D-FFI i ; -100 Nop'" Nop ri(T) = !Ono B-IO I' cor^p G2 0FF2 IC — * to OP Sin3 Op-^ — ^, to iC line

PAGE 153

136 Figure 44. Mpdating circuit

PAGE 154

137

PAGE 155

APPENDIX B COMPUTER CODES 1. Solution of TPBVP via Quasilinearization The FORTRAN program and the various subroutines which produced some of the results in Section B.l are given in the pages that follow. A brief summary of the code is given below. After the necessary data have been read in, the initial approximation is generated by integrating the non-linear differential equations. The (n+l)st approximation is obtained by integrating the particular and homogeneous equations, determining the unknown constants, and forming the linear combination of the particular and homogeneous solutions. This (n+l)st approximation is printed and stored as the nth approximation in preparation for the next iteration. If the initial values of the (n+l)st approximation agree with those of the nth approximation to within five or more decimal places, the problem may be considered solved. It has been mentioned before that the convergence is quadratic ahd hence the number of correct digits is roughly doubled with each , iteration. One fact that stands out 138

PAGE 156

139 about this method is that if the initial approximation is too far removed from the solution, the iterations may not converge. The initial approximation used in this work is determined from the analog computer results, ^his is found to be very reasonable and involves large savings in computing costs.

PAGE 157

140 2. Mimic Code It is desirable in any analog computer simulation to have a means for checking results. A digital computer program called Mimic ( 32 ) is used for this purpose. Mimic, a digital simulator program was developed for use on an IBM 7090/7094 computer with a FORTRAN IV IBSYS monitor. The University of Florida Computing Center personnel modified the program for use with the IBM 709 computer. Mimic's input language endows the digital computer, with the parallel nature of the analog computer and at the same time eliminates amplitude and tim.e scaling. If the various analog operations are represented by black boxes, then the equivalent box-oriented circuit diagram may be used for writing the mimic program. The following example will illustrate this fact. Suppose that the following equation is to be solved: X + X + X = where x(o) = 2 and x(o) = To solve the equation on the analog computer, it is rewritten in the following form X = -X X

PAGE 158

141 NEG 2. INT 0. INT h^^^ NEG ' "" Figure 45, Box-oriented mimic diagram 10

PAGE 159

142 The box-oriented circuit diagram is given in Figure 45 with the initial conditions specified on the top of the boxes . Figure 46 illustrates the format in which Mimic instructions are punched on an IBM card. "Results" are punched beginning in column 10, with "expressions" beginning in column 19. The variable names 2DX and IDX represent the second and first derivatives respectively, while X represents the desired solution. The IMT function integrates the first argument with respect to the independent variable, and the initial condition is the second argument. Mimic employs a high order Runge-Kutta subroutine for integration. Note that the Mimic instructions proceed very much like the corresponding analog computer program. With appropriate control instructions, the solution X may be obtained as a function of the independent variable. Mimic contains a large number of function subroutines, including logical functions. The ease with which its instructions can be presented and its great similarity to analog programming make it an ideal instrument for checking analog computer scaling and setup. Some of the analog computer results checked by this code are given in the MIMICsource language program in the pages that fellow.

PAGE 160

143 MIWIC PROGRAM NEI riRCN FUTX OPTIMISATION (N R S) CON ( A, B. ALPHA) PAR(XND) SOLVE EQHATIONS

PAGE 161

144 UIUIC PROGRAM NEirPRON FL'JX CPTIMISATION (N R E)

PAGE 162

LIST CF REFERENCES 1, J. Grey, "Nuclear Rockets," Nucleonics, Vol. 16, No. 7, July, 1958. 2 R R Mohler and J.E. Perry, Jr., "Nuclear Rocket Engine Control," Nucleonics, Vol. 19, No. 4, April, 1961. 3 Notes for the 1962 National Science Foundation Advanced Subject Matter Institute on "Nuclear Rocket Propulsion," Edited by C.W. Watson, University of Florida, 1964. 4. G. Seaborg, "Modern Control Theory," IRE Trans. Nucl. Sc. , Vol. NS-9. No. 1, January, 1962. 5. G. Leitman, Editor; Optimization Technique Wi th Applications to Aerospace Systems , Academic Press, 1962. 6. Z.R. Rosztoczy. and I.E. Weaver, "Optimum Reactor Shutdown Program for Minimum Xenon Buildup, "Nucl. Sci. Eng,, 20. 318-323 (1964). 7. G.R. Woodcock, and A.L. Babb, "Optimal Reactor Shutdown Programs for Xenon Poisoning," Trans. American Nucl. Soc. , 8, 1. p. 235 (1965). 8. P. A. Seeker, and L.E. Weaver, "Synthesis of Optimal Reactor Feedback Control Using the First Variation," Trans. Amer. Nucl. Soc, 8. 1, p. 244-245 (1965). 9. I. Kliger, "Synthesis of an Optimal Nuclear Reactor System," University of Arizona, Symposium on Neutron Dynamics and Control, April 5=7, 1965. 10. Z.R. Rosztoczy, A. P. Sage, and L.E. Weaver, "Application of Pontryagin"s Maximum Principle to Fluxstate Changes in Nuclear Reactors," TID=7662 AEC Symposium Series, 2, Reactor Kinetics and Control, 1964. 11 R.R. Mohler, "Optimal Control of Wuclear-Reactor Processes," Ph.D. Thesis (Michigan, 1964). 145

PAGE 163

146 LIST OF REFERENCES (Continued) 12. L.E. Weaver. D.G. Schultz . et al. , "Research in and Application of Modern Automatic Control Theory to Nuclear Rocket Dynamics and Control," NASA CR-193, 'fniversity of Arizona. Tucson, Arizona, March, 1965, 13. R.H. Goddard, "A Method for Reaching Extreme Altitudes," Smithsonian Misc. Collections, Vol. 71, No. 2. 1921. 14. R.E. Kalaba. "On Nonlinear Differential Equations, the Maximum Operation and Montotone Convergence," J. Math, and Mech., Vol. 8 (1959). 519=574. 15. R.E. Bellman and R.E. Kalaba, "Dynamics Programming, Ivariant Imbedding and Quasilinearization : Comparison and Interconnections," Rand Memo RM-4038=PR , March, 1964, AD 431 866. 16. L.S. Pontryagin, et al. , The Mathematical Theory of Optimal Processes . Interscience Publishers, John Wiley and Sons, New York, 1964. 17. R. Sridhar. et al . , "Investigation of Optimization of Altitude Control Systems," JPL Report TR-EE65=3, Purdue University, School of Electrical, Lafayette. Indiana, 1965. 18. R.E. Bellman and H.H. Kagiwada and R.E. Kalaba, Quasilinearization, Boundary Value Problems and Linear Programming, Memo RM-4284-PR, September, 1964. 19. R.E. Bellman and R.E. Kalaba, Quasilinearization and Boundary ^^alue Problems . American Elseview Publishing Co. . New York. 20. R. McGill and P. Kenneth, "A Convergence Theorem on the Iterative Solution of Nonlinear Two-Point Boundary ^'alue Systems, "XlVth International Astronautical Federation Congress, Paris, France, September. 1963. 21. A. P. Sage and B.R. Eisenberg, "Closed Loop Optimization c Pixed Configuration Systems," Int -^rnatlonal Journal of control, 'cl. 3. ^Io. 2, 1966. 22. A. P. Sage and T.W. Ellis, "Sequential Suboptirral Adaptive Control of Nonlinear Systems," Proceedings National Electronics Conference, Cctober, 1966.

PAGE 164

147 LIST OF REFERENCES (Continued) 23. H.P. Smith, Jr., "Dynamics and Control of Nuclear Rocket Engines," Ph.D. Dissertation, Massachussets Institute of Technology, August, 1960. 24. A.H. Stenning, "A Rapid Approximate Method for Analyzing Nuclear Rocket Performance," J. Am. Rocket Soc, 30. 169-172. 1960. 25 H P. Smith, "Closed Loop Dynamics of Nuclear Rocket Engines with Bleed Turbine Drive," Nuclear Science and Engineering, 14. 371-379, 1962. 26. H.P. Smith, "Closed Loop Dynamics of Nuclear Rocket Engines with Topping Turbine Drive." Nuclear Science and Engineering, 18, 508-513. 1964. 27. A.H. Stenning, and H.P. Smith, "The Stability and Control of Nuclear Rocket Engines," Presented at the ARS Semi-Annual Meeting, Los Angel.es, May, 1960. 28. Personal Communications between R.E. l.fhrig and J. Saluja. 29 J."^. Humphries, "Model Reference Adaptive Control of a Nuclear Rocket Engine," Ph.D. Dissertation, University of Florida, April, 1966. J.T. Humphries, R.E. IJhrig, and A. P. Sage, "Adaptive, Con30. trol Concepts for Nuclear Rocket Propulsion Systems," AIAA Second Propulsion Joint Specialist Conference. C Colorado Springs, Colorado, June, 1966. 31. A. P. Sage and T.W. Ellis, "Optimal and Suboptimal Guidance and Control for Low Thrust Orbital Transfer," Proceedings, Symposium in Missiles and Aerospace Vehicle Sciences, American Astronautical Society, 1966, 32. H.E. Peterson and L.M. Warshawsky, "MIMIC A Digital Simulator Program," SESCA Internal Memorandum 65-12, Wright-Patterson Air Force Base, 1965.

PAGE 165

BIOGRAPHICAL SKETCH Jagdish Saluja was born January 14, 1934 in district Jhelum, now in Pakistan. He graduated from St. John the Baptist High School in June, 1951 and received the Bachelor of Science degree in Physics froiTi the "niversity of Bombay in June, 1955» He went to the I'niversity of Michigan in 1956, receiving his Bachelor of Science degree in Electrical Engineering In August, 1957 and a Master of '^^cience degree in Nuclear Engineering in August, 1959„ He then joined the staff of the Reactor Engineering Division at Argonne National Laboratory. He began graduate study at the University of Florida in September, 1963, and until the present time has pursued work toward the degree of Doctor of Philosophy, He is the son of Mr. and Mrs. Kirparam Saluja and the third of their eight children, six boys and two girls. 148

PAGE 166

This dissertation was prepared under the direction of the chairman of the candidate's supervisory committee and has been approved by all members of that committee. It was submitted to the Dean of the College of Engineering and to the Graduate Council, and was approved as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December, 1966 Dean, Graduate School

PAGE 168

^y/^-^