
Citation 
 Permanent Link:
 https://ufdc.ufl.edu/UF00097875/00001
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:
 1966
 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 )
nonfiction ( marcgt )
Notes
 Bibliography:
 Bibliography: leaves 145147.
 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 nonprofit 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 startup 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= 104) Comparison
of analog computer results with those ob
tained by quasilinearization (Nuclear Reactor
System). .............. . 42
7. Power level versus time ( = 103) Comparison
of analog computer results with those ob
tained by quasilinearization (Nuclear Re
actor System) . . . . . . .. 43
8. Power level versus time (A= 103) 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. Pressuretime plot, tf = 4 secs, 6 secs and
11 sees (Nuclear Rocket Engine System) . 88
30. Optimal control poison reactivity and temper
aturetime plots, Td = 0.4 (Nuclear Rocket
Engine System) .. . . ..... 89
31. Optimal pressuretime 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 pressuretime 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. Boxoriented 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 (ORl).
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 ith 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
closedloop control system is realized by letting control
rod rcvement be a proportional plus integral function
of neutron density error.
The twopoint 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 startup 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
twopoint 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, timevarying
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 suboptimal 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 twopoint boundary value problem. This technique
was initially suggested by Bellman and Kalaba (15) as a
modification of the NewtonRaphsonKantorovich method.
This converts the problem of solving a twopoint boundary
nonlinear differential equation into a problem of solving
a sequence of twopoint 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(nnd)2 ] dt
where PC = a + A(nnd)
a and A are constants to be determined.
iv) J = 1/2 tf [ Pc2 + a(nnd)2] dt
o
where Pc = A(nnd) + Bf (nnd) 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 semiautomatic
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
delayedneutron 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 (nnd)2] dt
in fixed time tf
3) J = 1/2 jf [ Pc2 + a(nnd)2 ]dt
where pc = a + A(nnd)
in fixed time tf
4) J = 1/2 [ PC2 + a (nnd)2 ] dt
tf
where pc = A(nnd) + B /f (nnd) 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 pdimensional vector which is
a known function of the state and E is a qdimensional
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 reformulation 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
twopoint boundaryvalue 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 nonlinear twopoint and
multipoint 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 2Ndimensional
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
twopoint 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 righthand 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 nonlinear 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) = (PB)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 sec1
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.454x105 0.629x10l5 0.408xl05
5no 1.263x105 0.631x105 0.268x105
6no 1.51 x105 0.552x105 0.235x105
1Ono 1.945x105 0.457x105 0.176x105
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
c12A 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 trackandhold
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 twentynine 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 C12A 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 subproblems.
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 C13A
20 C1 C12A
18 C15A
16
14
12
.0
0 5 10 15 20 25 30 35
50     0%
C12A
C14A C15A /C13A
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 = 5x105pln26.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+lplk)
6.4(nk+lnk) + 0.1 (ck+lck)
But
nk = 5xl05pk(nk)2 6.4nk + 0.1 ck
Therefore
hk+l = 5x105plk(nk)26.4nk + O.1ck + 5x105 [p1k2nk
(nk+l nk) + (n)2lk+l k)] 6.4(nk+lnk)
+ 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 = 106xox11106 3lO o(xO)+5x105x3 (x10)2
6.4x11+0.1x21 (2.55)
Linearizing Equation (2.51) there results
ak+l = k + 6.4(nk+lnk) 0.1(ck+lck)
= 6.4nkC0.ck+6.4(nknk+lk) 0.(ck+lck)
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 ox315xlO5(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 10x36.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(ttf) = 0
Pl0 = 0.5 x 105 Pl(tf) = 0.30xl05 (2.75)
p20 = 0.18 x 105
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, = 103secs and A = 104secs. 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 = 103
Figure 8 shows power level versus time plot of run
C12A 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 = 104)
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
o20
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  103)
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 = 103)
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 (nnd)2 dt
The twopoint 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(nnd) 1000 pcP1 + 6.4 P16.4P2 (2.78)
P2 =0.1 P1 + 0.1 P2 (2.79)
Pc =1000 nP1 (2.80)
with boundary conditions:
n(o) = 50
c(o) = 64no (2.81)
n(tf) = 5Oa
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 (nnd)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 (nnd)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 (Nd)2dt
2 =0
J2 = 5 x 104 J21
where
tf 2
2' = a J O.l(NNd) 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.0149x105 0.243x105
5no 0.696x105 0.527x105 0.875x105
6no 1.323x105 0.823x105 1.484x105
1Ono 1.914x105 2.88x1075 3.84x105
.
a = 10
..........................................................
2no 0.456xl05 0.0297x105 0.258x105
5no 0.696x105
6no 1.323x10 
10no 2.058x105 4.95x105 6.0x105
a= 20
..........................................................
2no 0.456x1i05 0.0595x105 0.288x105
5no 0.745x105 1.92x105 2.30x105
6no 1.445x105 2.81x05 3.53x105
1Ono 2.126x105 9.2x105 10.26x105
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 (nnd)2 dt
pc = a + A(nnd)
The twopoint boundary value problem resulting by
the application of the maximum principle is:
1 = 1000 an + 1000 A(nnd) 6.4n + O.lc (2.82)
S= 6.4n 0.lc (2.83)
A = 0. a = 0 (2.84)
p, = (A2 + a )(nnd) aA 1000 Pl(a+2AnAnd)
+ 6.4P1 6.4P2 (2.85)
P2 = 0.1P1 + 0.1 P2 (2.86)
P3 = A(nnd)2 1000 npl(nnd) a(nnd) (2.87)
P4 = a A(nnd) lO00npl (2.88)
Pc = a + A(nnd) (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(NNd) 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 (nnd)2 dt
tf
PC = A(nnd) + B 1 (nnd)dt
tf
Let x = J (nnd)dt Hence x = (nnd). and p = A(nnd)+Bx
0
The twopoint boundary value problem resulting by the appli
cation of the maximum principle is:
n = 1000 A n(nnd) + 1000 Bxn 6.4n +0.1c (2.91)
c = 6.4n O.lc (2.92)
x = nnd (2.93)
A = 0 B = 0 (2.94)
PI = (A2+a )(nnd) ABx 1000 P1 [2AnAnd+BxJ
+ 6.4P1 6.4P2 P3 (2.95)
P2 = 0.1P1 0.1P2 (2.96)
P3 = B2x 1000 BnPl AB(nnd) (2.97)
P4 = A(nnd)2 1000nPl(nnd) Bx(nnd) (2.98)
P5 =Bx2 1000 xnP1 A x (nnd) (2.99)
PC = A(nnd) + 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 105
5no 3.695 79.78 11.34 0.294 x 105
2no 3.836 128.4 0.105 x 105
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(nnd)
tf
and Pc = A(nnd) + B J (nnd)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 103
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 closedloop 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 nonlinear
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 solidcore, 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 closedloop control is desired. The maximum
principle however, leads to an openloop control. The
closedloop 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 SmithStenning
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 sec1
106 oR1
0.3 OR/psi
104 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 MohlerPerry ( 2 )
formulation. In view of this modification, equation (3.5)
is rewritten 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 UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EUY1EAJF9_DGS0KF INGEST_TIME 20170719T20: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 startup 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. Pressuretimie plot, tf = 4 sees, 6 sees and 11 sees (Nuclear Rocket Engine System) . . 88 30. Optimal control poison reactivity and temperaturetime plots, T(j = 0.4 (Nuclear Rocket Engine System) . . o . o o . Â» . c . . = . 89 31. Optimal pressuretime 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 pressuretime 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. Boxoriented 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^ ith 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 closedloop control system is realized by letting control rod nrcvement be a proportional plus integral function of neutron density error. xiv
PAGE 16
The twopoint 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 startup 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, timevarying 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 suboptimal 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 NewtonRaphsonKantorovich method. This converts the problem of solving a twopoint boundary nonlinear differential equation into a problem of solving a sequence of twopoint 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(nnd)^ ] dt where ^c = a + ACnn^) a and a are constants to be determined, iv) J = 1/2 ) ^ [ ^c^ + a(nn.)^J dt where P^ = Atnn^j) + b) (nn^) 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 semiautomatic 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 delayedneutron 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{nn^f] dt in fixed time tj 3) J = 1/2 J^ \^ Pc2 + a(nnd)2 1 dt where p^ = ^ + A(n=n^) in fixed time tÂ£ 4) J = 1/2 J^ [ Pc2 + a (nn^)^ ] dt where p^ A{nnci) + B f (nn^) 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 pdimensional vector which is a known function of the state and F is a qdimensional 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 reformulation 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 twopoint boundaryvalue 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 nonlinear twopoint and multipoint 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 2Ndimensional 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 twopoint 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 (PXm^^)) (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.: m1 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 nonlinear 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^n6.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 C12A 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 trackandhold 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 twentynine 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 C12A 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 subproblems. 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% CI3A 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.422.45) and modified by (2.46) as : (2.50) n = 5xl0^pin26.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^^lnk).(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^*ln' + 0.1 (c'^^lc'')
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^l106x3O(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 Â•N1 =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^X36.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 C12A 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(nn^)^ 1 dt The twopoint boundary value problem resulting by the application of the maximum principle is: (2.76) c = 6.4n 0.1 c (2.77) P^ = a(nnd) 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) = 50 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 (nnj)^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 (nnd)^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 (NNd)2dt J2 = 5 X lO""^ J2 where J2" = a J O.l(NNjj) 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 (nn^) o dt p Q = a + A(nnci) The twopoint boundary value problem resulting by the application of the maximurr. principle is: n = 1000 an + 1000 A(nn^) 6.4n + 0.1c (2.82) c = 6.4n 0.1c (2.83) (2.84) A = 0. a = P^ = (a2 + a )(nn^) aA 1000 P^ (a+2AnAnj) + 6.4?^ 6.4P2 P2 = O.lPi + 0.1 P2 (2.85) (2.86) P^ = A(nnj)^ 1000 np^(nn^) a(nnd) (2.87) P4 = a A(nnd) lOOOnpj^ (2.88) Pc a + A(nnd) (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(NN(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(nn^) + B J (nnd)dt o Let X = J (nn^^)dt Hence x = (nn^). and p^ = A(nnc) + Bx o The twopoint boundary value problem resulting by the application of the maximum principle is: n = 1000 A n(nnd) + 1000 Bxn 6.4n +0.1c (2.91) c = 6.4n 0.1c (2.92) X = nn^j (2.93) A = . B = (2,94) P2 = (A^+a )(nn(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(nnj) (2.97) P4 = A(nn^)^ lOOOnP^(nn^) Bx(nn^) (2.98) P5 =Bx2 1000 xnP^ A x (nn^) (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 Systeni)
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(nnd) + B j (nnj)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 nonlinear 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 solidcore, 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 SmithStenning 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 startup 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 Startup 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 reinserted to maintain the final temperature at the desired level. This reinsertion 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. Pressuretime plot, tf = 4 sees, 6 sees and 11 sees (Nuclear Rocket Engine System)
PAGE 106
TIME (sec) Figure 30. Optimal control poison reactivity and temperaturetime 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 temperaturetime 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 = lOYpn .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 openloop control and the specific optimal control of nuclear reactor processes in a nuclear reactor system and openloop 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 closedloop control suffers in performance as compared to the openloop 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 closedloop controller parameters as a function of neutron density and allowed startup time, as well as system parameters, if performance near optimum is to be achieved for varying startup 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 openloop control problem for the nuclear rocket engine system was solved in three different parts: a) variation in the startup 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 startup 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 reinserted and this reinsertion is very gradual for longer startup times as opposed to shorter startup 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 startup 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 startup 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 squareroot and division are formed by suitable combinations of quartersquare 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 Â— hAAA, 1 SUhiWU'JG Ar.lPLiFlSR AAA.^Â— Rg _ ^AA.\. r .t n c / E liec't IMTEGRATUJG AMPLIFIER 0, te 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 AB AB 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 (NNj) P2 = 0.1 P^ + 0.1 P2 Pc =0.1 N Pi J = 1/2 J p^2 dt + 5x10Â°^ a J O.l(NNd) 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^)(NNj) 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 (NN^)^Wi (NNj) 0.001 a (NNj) P4 = 10""^ a 0.1 N P^ Pc = 10"^ a + 0.01 A (NNj), A is negative Problem 4 ; Fi = AN (NFid) +0.1 BxN 0.64 N + C C = 0.0064 N 0.01 C X = 0.1 (NNd) P^ = 0.01 (A^+a )(NN^) = 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'^ (NN^) P4 = 0.01 A(NN^)2 NPx(N=N^) lO"^ B x(NNd) Pp^ = =0.01 A(NNd)^ 0.01 (O.lB)x^ xNPj^ P^ = 0.01 A(NNd) + 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(NN ) dt o o J = J, + 5x10Â°^ a JoÂ° h = 1/2 1 pM where ^^ 'o and J2' = J 0.1(NN.) 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 rH 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 Repop 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 Repop 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 trackandhold 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 readout (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 online 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 = 002 ^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 n0 CIC 5100 CO MP G22 "[> fr Â•MOO PULSHR .T *l GIG iC\To OP OPJCriC lir;'. u 4J r~ p.?"'^ ;^Â• E OP TiiuG T = p/Q= ISec 65 t CDE n{T) ^Â•n(T) n(t) _.. J : Qh) A!0 cor. IP Gl DFFI i ; 100 Nop'" Nop ri(T) = !Ono BIO 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 nonlinear 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 boxoriented 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, Boxoriented mimic diagram 10
PAGE 159
142 The boxoriented 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 RungeKutta 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. NS9. 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. 318323 (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. 244245 (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 WuclearReactor 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 CR193, '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 RM4038=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 TREE65=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 RM4284PR, 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 TwoPoint 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. 169172. 1960. 25 H P. Smith, "Closed Loop Dynamics of Nuclear Rocket Engines with Bleed Turbine Drive," Nuclear Science and Engineering, 14. 371379, 1962. 26. H.P. Smith, "Closed Loop Dynamics of Nuclear Rocket Engines with Topping Turbine Drive." Nuclear Science and Engineering, 18, 508513. 1964. 27. A.H. Stenning, and H.P. Smith, "The Stability and Control of Nuclear Rocket Engines," Presented at the ARS SemiAnnual 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 6512, WrightPatterson 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/^^

