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
