• TABLE OF CONTENTS
HIDE
 Front Cover
 Dedication
 Acknowledgement
 Table of Contents
 List of Tables
 List of Figures
 List of symbols
 Abstract
 Introduction
 Optimal control of a nuclear reactor...
 Optimal control of a nuclear rocket...
 Summary and conclusions
 Appendices
 References
 Biographical sketch
 Back Cover














Title: Computational techniques and performance criteria for the optimum control of nuclear system dynamics
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00097875/00001
 Material Information
Title: Computational techniques and performance criteria for the optimum control of nuclear system dynamics
Alternate Title: Optimum control of nuclear system dynamics
Physical Description: xvi, 148 leaves : ill. ; 28 cm.
Language: English
Creator: Saluja, Jagdish Kumar, 1934-
Publisher: University of Florida
Place of Publication: Gainesvill
Gainesvill
Publication Date: 1966
Copyright Date: 1966
 Subjects
Subject: Controlled fusion   ( lcsh )
Nuclear reactors   ( lcsh )
Nuclear rockets   ( lcsh )
Nuclear Engineering Sciences thesis Ph. D   ( lcsh )
Dissertations, Academic -- Nuclear Engineering Sciences -- UF   ( lcsh )
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Bibliography: Bibliography: leaves 145-147.
Additional Physical Form: Also available on World Wide Web
General Note: Manuscript copy.
General Note: Thesis - University of Florida.
 Record Information
Bibliographic ID: UF00097875
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000543376
oclc - 13142428
notis - ACW7084

Downloads

This item has the following downloads:

PDF ( 4 MBs ) ( PDF )


Table of Contents
    Front Cover
        Page i
        Page i-a
    Dedication
        Page ii
    Acknowledgement
        Page iii
        Page iv
    Table of Contents
        Page v
        Page vi
    List of Tables
        Page vii
    List of Figures
        Page viii
        Page ix
        Page x
        Page xi
    List of symbols
        Page xii
        Page xiii
    Abstract
        Page xiv
        Page xv
        Page xvi
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
    Optimal control of a nuclear reactor system
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
    Optimal control of a nuclear rocket engine
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
    Summary and conclusions
        Page 99
        Page 100
        Page 101
        Page 102
    Appendices
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
    References
        Page 145
        Page 146
        Page 147
    Biographical sketch
        Page 148
        Page 149
    Back Cover
        Page 150
        Page 151
Full Text












COMPUTATIONAL TECHNIQUES AND
PERFORMANCE CRITERIA

FOR THE OPTIMUM CONTROL OF
NUCLEAR SYSTEM DYNAMICS






By
JAGDISH K. SALUJA


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


UNIVERSITY OF FLORIDA
December, 1966


























DEDICATED TO


'om, Dad, My Brothers and

Sisters, and all my friends

for their constant ~ncourage-

ment.











ACKNOWLEDGMENTS


The author wishes to acknowledge the help of his

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

able for consultation and who has been a source of continued

encouragement and enlightening advice. The author also

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

"hrig for his encouragement during the course of his

pursuance of the doctoral degree and for being on his

supervisory committee. His help is greatly appreciated.

The concern and interest of the remaining members of the

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

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

upon the works progress and is sincerely appreciated.

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

for his help with the digital computer computation, using

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

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

Nair for many helpful discussions. The cooperation and

the assistance, especially of Mr. Kenneth Fawcett on the

Analog computer is greatly appreciated.

Furthermore, the author wishes to thank the depart-

ment of Nuclear Engineering Science for the financial

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

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











Mazelle Shadburn for her assistance in typing the disserta-

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

the drawings.













TABLE OF CONTENTS


Page


DEDICATICN . . . . . . . . . ii

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

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

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

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

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

CHAPTER


INTRODUCTION . .

OPTIMAL CONTROL CF
SYSTEM ....


A NUCLEAR REACT. .
A NUCLEAR REACTCR


A. Introduction and Model Description.

1 Statement of the Maximum Principle.

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

3 Gradient Technique . . .....

4 Quasilinearization Technique .....

B. Formulation of Specific Problems. .

1 Cost Function 1 . . . . .

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

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

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











TABLE OF CONTENTS (Continued)


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

C. Discussion of Results. . . .

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

Introduction and Model Description

Solution Procedure . .....

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

IV. SUMMARY AND CONCLUSIONS . . . .

Summary and results. . .

Suggestions for further study... .

APPENDICES . . . . . . . . .

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

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

LIST (F PEFFRFNCES . . . . . . .

BIOGRAPHICAL SKETCH . . . . . . . .


Page

58
59

67


77

77

85

85
86

94

99

99

102

103

104
104
110
126

138

138
140

145

148













LIST OF TABLES

Table Page

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

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

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

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

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

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

7. Scale Factors for the Four Control Problems. 113

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

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

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

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

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












LIST OF FIGURES


FIGURE PAGE

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

2. Cost function versus desired power level . 24

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

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

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

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

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

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

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

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









LIST OF FIGURES (Continued)

FIGURE PAGE

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

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

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

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

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

17. Cost function versus peak power attained 63

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

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

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

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

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

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











LIST CF FIGURES (Continued)

FIGURE PAGE

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

25. Cost function versus peak power desired. . 74

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

27. Schematic Diagram of a nuclear rocket engine 79

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

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

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

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

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

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

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

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

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











LIST OF FIGURES (Continued)

FIGTRE PAGE

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

38. Symbols used in analog computer diagrams. . 105

39. Logic symbols used in computer diagrams . 107

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

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

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

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

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

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

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













LIST OF SYMBOLS


n Neutron density.

nd Desired neutron density.

N Scaled neutron density.

Nd Scaled desired neutron density.
No Normalized neutron density.

c Delayed neutron precursor density.

C Scaled delayed neutron precursor density.

C' Normalized delayed neutron precursor density.

\ Nean effective neutron lifetime (sec).

P Core inlet pressure (psi).

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

T Maximum core surface temperature (OR).

Td Design maximum core surface temperature (oR).

T' Normalized maximum core surface temperature.

V Valve setting.
Vg Normalized valve setting.

a Error weighting term,

at Temperature coefficient of reactivity (OR-l).

al Normalized temperature coefficient of reactivity
(dollars).

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









LIST OF SYMBOLS (Continued)


a. Normalized propellant density coefficient of re-

activity (dollars).

B Delayed neutron fraction.

A Decay constant (sec ).

P Total reactivity.

P' Normalized total reactivity (dollars).

Pc.r Control poison reactivity.

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

-p Pressure time constant (sec).

A,B Ccnstants in the control law.

A.E Scaled constants in the control law.

H Hamiltonian.

J Cost function.

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

t Time.

to Initial tire.

ui Ccnponent of the control vector u.

xi Corpcnent of the state vector x.


xiii










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

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

by

Jagdish K. Saluja

December, 1966

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

The main objective of this study is to investigate

optimum control of the dynamics of nuclear systems and

associated computational problems. Two systems are con-

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

Engine System. For the first system, three problems are

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

ferred from one steady state condition to another with a

quadratic constraint on control rod movement. By varying

the operating time and placing an inequality constraint on

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

density is minimized with quadratic integral constraint

on control rod movement, c) a fixed configuration optimum

closed-loop control system is realized by letting control

rod rcvement be a proportional plus integral function

of neutron density error.









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

sulting by the application of optimization theory to point

reactor models are obtained. A modified gradient technique

is used to exploit the hybrid capabilities of modern analog

computers. The results are compared with the solutions

of this problem by the numerical technique of quasilinear-

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

no more than four iterations.

In the examples considered here, the cost function

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

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

of closed loop controller is used. The coefficients A

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

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

final state. Variation in system parameters may, however,

render the closed loop control better than the open loop

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

only a quadratic constraint on total reactivity is con-

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

core 25,000 megawatt, million pound thrust model developed

by Smith and Stenning. Here, the results obtained for

various start-up times and various desired power levels in

terms of control needed and required cost, appear very










reasonable. Also it was observed that if the pressure and

the temperature time constants were increased from their

normal values, the temperature response became somewhat

sluggish. Possible improvements to this situation are

suggested.











CHAPTER I


INTRODUCTION


Just as all other fields of scientific thought have

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

field of rocketery. From the Chinese rockets of 1232

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

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

vital factor in the success of future space travel. This

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

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

ther strengthened by the successful testing of the KIWI

and NERVA reactors.

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

of which a very instructive book on nuclear rocket system

optimization by Bussard (5) has appeared. An outgrowth

of this interest is the optimal control problem wherein

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

From the nuclear literature, however, it is clear that

the subject of optimal control of nuclear systems has

received little interest. Recently a number of papers

have appeared in this area, ranging from optimal reactor

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










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

More recently Mohler and Weaver have reported some work

in the application of modern optimal control theory to

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

still remains to be done.

The mathematical theory used for the study of

optimization problems is the calculus of variations. It

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

technique as an important tool in the performance analysis

of rockets in his early paper published 40 years ago.

There are many possible variational methods fcr minimizing

or maximizing a functional over a function space. Among

the commonly used methods in control system design are:

(1) calculus of variations

(2) maximum principle
(3) dynamic programming

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

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

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

the calculus of variations as well as from the method

of dynamic programming. Pontryagin's maximum










principle is a very elegant method, which is applicable

to a broader class of problems than the calculus of

variations. An outcome of this principle is the so called

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

with split boundary conditions.

This method will be applied in the study of dynamics

of nuclear systems. Since the dynamic response of nuclear

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

differential equations, conventional analysis is usually

too difficult to give a closed form solution. Develop-

ment of efficient computer techniques is therefore necessary.

A numerical approach which is currently quite

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

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

erty that each successive set of solution functions yields

an improved value for the functional being optimized.

Another numerical technique of considerable interest

is the approximating process of quasilinearization which

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

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

ginal two-point boundary value problem. This technique

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

modification of the Newton-Raphson-Kantorovich method.

This converts the problem of solving a two-point boundary

non-linear differential equation into a problem of solving

a sequence of two-point boundary linear differential equation










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

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

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

i) J = 1/2 Pc2 dt

a) n(tf) = nd

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

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

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

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











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


where PC = a + A(n-nd)

a and A are constants to be determined.

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

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

A and B are constants to be determined.


The last two criteria give closed loop control and since

a specific form is chosen for the controller, it will

henceforth be referred to as the specific optimal control

problem.

The above four controls are analyzed and compared

with one another for two reactor models. The two models

analyzed are:

(1) Point model of a nuclear reactor system.

(2) Point model of a nuclear rocket engine.

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





6





the same pattern as Chapter II. Chapter IV includes

a discussion of results and suggestions for further

study.












CHAPTER II


OPTIMAL CONTROL CF A NUCLEAR REACTOR SYSTEM


A. Introduction and Model Description


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

changes are generally accomplished by semi-automatic

means in small research reactors. In large power re-

actors and in nuclear rocket engines, however, automatic

control is desired. If this control is optimal in some

specified sense, economic or other advantages could be

gained.

The reactor model considered here is that of a one

delayed-neutron group without temperature feedback and

without spatial and angular dependence, described by the

following system of differential equations:


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


S= n Xc (2.2)

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

tion are the state variables and the reactivity p is








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

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

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

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

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


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

in fixed time tf


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

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


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


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

in fixed time tf.









1. Statement of the P'aximum Principle

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

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

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

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

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



Assume a performance criterion of the type:

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

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

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


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










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


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

Auxiliary or adjoint, or costate functionsPi are defined by


p (t) = (2.8)
axi

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


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

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

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

t(to) = Xo (2.10)

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

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










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

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

It is desired to determine the unknown parameters in a

control law of the form

= h(,) = 0 (2.12)

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


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

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

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











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

x = f(x.F.t)

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

and b = o

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

3. Gradient Technique

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










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


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

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

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

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


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

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


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










4. Quasilinearization Technique

This technique,developed by Kalaba, Bellman and others

is particularly suited for the type of optimization problems

being considered here. It provides an effective computa-

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

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

technique consists of producing a sequence of approxi-

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

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

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

asymptotically the number of correct digits in each approx-

imation is doubled.

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

x = (;) (2.19)

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

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

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

ential of (2.19) and obtain












ax


or iik+l 1 ik =


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

i = 1.2....


From (2.19) above


i k =
Therefore
ik+l 1
Xj


(2.23)


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


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

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

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

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


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


(2.25)


(2.21)


%t( k)










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

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

Subject to the initial conditions

p(o) =o

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

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


ti = hi (2.27)
m"1 "axm

The initial vector has the ith row unity and the others

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

leads to the following linear algebraic equations:

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

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


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










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

B. Formulation of Specific Problems

1. Cost Function 1

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

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

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

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

Then by (2.8)

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

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

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


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


(2.33)










and

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

Setting = 0, yields


Pin (2.35)

The following system of equations are therefore obtained.

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

= n .c (2.37)

l = ( c -B P P2 (2.38)

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

with boundary conditions:

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

P2(tf) =0.0









Simulation of problem 1

The simulation of the problem as described above, for

nuclear reactor system dynamics, was performed with an

Applied Dynamics analog computer. A detailed computer

wiring diagram used for programming the computer is

given in Figure 41 Table 1 gives a list of system

constants and the initial conditions used. The system

equations described in Section 1, including the constants

of Table 1 are given below.

Problem 1:

n = 1000 pcn 6.4n + Ic (2.42)

S= 6.4n .lc (2.43)

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

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

With boundary conditions:

n(o) = 5.0

c(o) = 64n,

n(tf) = 5.0a

P2(tf) =0.0 (2.47)













TABLE I


NUCLEAR REACTOR CONSTANTS FOR A LUMPED SYSTEM AND INITIAL
CONDITIONS




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

Decay constant, 0.10 sec-1

Average neutron life time,A 103 sees

Delayed neutron fraction,B 0.0064



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

Meutron density or power level, no 5K.W

Precursor population, co 64no











Simulation of the system of equations on an analog

computer requires amplitude scaling. This scaling is

necessary because all voltages must be kept within + 100

volts in order to prevent overloading of the machine

components. This overloading will lead to inaccurate

results. The basic set of scaled equations are rewritten

here from Appendix A for convenience. The bars indicate

scaled quantities.

Problem 1:

N =-100 N1P 6.4 N + 10 C

C = 0.064N 0.1 C

P1= 100 N P + 6.4 P1 64 P2


P2 = 0.1 P + 0.1 P2

Pc =-0.1 N P1

With boundary conditions:


S= 0.05

Co = 0.032


V(tf) = aNo


P2(tf) = 0.0










Typical solution procedure: With the equations scaled, the

modified gradient method will be used for optimization on

the analog computer and the results verified with the num-

erical technique of quasilinearization.

The gradient technique is preferred as opposed to

the trial and error procedure which can be quite cumber-

some and extremely time consuming, particularly for high

order systems. A brief description of this technique is

given in Section A.4.

The equations used for simulation are on page 21,

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

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

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

in one second. Figure 1 shows the neutron density plots

and the reactivity function plots. Table 2 shows the

initial values of the adjoint variables and the cost

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

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

the reactivity will have to be programmed differentially

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

kinetic equations ( 10,11).

The cost function is plotted versus the peak power

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

increase in cost with increase in power level and this in

general is to be expected.







Nd=0no Nd= no \














o



4
0 Nd= 5no










\ o
40















4 0 0

0 Ndd 2 no





So o
50
1 ^ Nd [O no
40-





















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

















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

































0








F-
Z







0







0.5


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

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








TABLE 2

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


Nf J P10 P20


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

5no 1.263x10-5 0.631x105 0.268x105

6no 1.51 x105 0.552x105 0.235x105

1Ono 1.945x105 0.457x105 0.176x10-5




The gradient technique is studied in greater detail

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

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

produces the desired optimal control u* (reactivity in

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

all Pi values are scanned, much more general information

concerning the control of the system may be obtained.

In particular the reachable set of the system may be

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

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

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

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

with admissible controls. By taking a suitable set of










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

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

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

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

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

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

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











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

above conditions are satisfied are held on the n analog

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

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

the increments 8Pi are selected by the potentiometers.

The convergence of the iterations is very sensitive to

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

overshoots in the desired trajectory (see comment column

in Tablel2). Smaller values result in an increased

number of iterations to converge Figure 4 The

tolerance on the end conditions during these tests are

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

obtain high accuracy in the terminal conditions as this

would require use of more analog and logic elements and

this is not warranted at this stage. Figure 3 shows

the sequence of trajectories converging to the optimal

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

convergence was obtained in twenty-nine trials. Total

time for convergence was 60 seconds. A graph of number

of iterations for convergence versus Pl10 is shown in

Figure 4 for various step sizes. Note that when P10

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

tions for convergence is infinite; Figure 5










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

and the value of P20 versus the iteration number in a

particular iteration cycle. With successive iterations

the value of P20i decreases until an optimal value is

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

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

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

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

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

lated to maximum n.

The steps in the iterative procedure are sufficiently

complex that it is wise to break the programming down into

a series of sub-problems.

An iterative procedure consists of the following steps:

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

take on arbitrary values.

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

will help to determine sign of B P2. During

the iteration process, the addition is per-

formed by the analog memory so that the new

value


P2 201 + 8P21










where


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

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

Pj+l = P hP1


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




30

8 --2-




'0 6 1 /
/ ITERATION

5 a


-

[ 3
0





2 2




40 DESIRED TRAJECTORY


2727h
-J 0


o th
-J 20t








10 0 o
20 -





2 4 6 LO
TIME (Sec.)

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


















































Figure 4. Number of iterations for convergence versus
P10












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





4o o
0


I-








S30









"20





0 .








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









28

26

24

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

18 C-15A

16

14

12

.0
0 5 10 15 20 25 30 35

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

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

40




c 30




20








0
5 0 5 20 25 30 35

ITERATION NO

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








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

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

c = 6.4n 0.1c (2.51)

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

P2 =-0.1P1+0.1P2 (2.53)

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

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

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

But

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

Therefore

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

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

+ 0.1 (ck+l_-k)









Hence

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

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

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

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

6.4x11+0.1x21 (2.55)


Linearizing Equation (2.51) there results

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

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

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

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

x21 = 6.4x11 O.1x21

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

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

+ 6.4x31 6.4x41 (2.57)









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

The above equations may be written in the following form:

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

21 = a21x + a22x21 (2.60)

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

x2 b2
x3 b3
x4 b4

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

b = 106(x0)2 x 30

b2 = b4 = 0

b3 : 106x0(x3o)2

all = 106 x 10x3-6.4

a12 = 0.1

a13 = 5x10 5(xo)2











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

a21 = 6. 4
a22 = 0.1

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

a34 = -6.4

a43 = -0.1

a44 = 0.1



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

the coefficients are time varying. Equation (2.63) may

also be rewritten as

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

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


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

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


P2 = a21i1 + a22P2


(2.67)










p3 = b3 + a31 + a33 (2.68)

P4 = a43P3 + a44P4 (2.69)

Appendix B shows incorporation of equation (2.65) in the

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

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

or X = A(t)5 (2.71)

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

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

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

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











The complete solution is therefore

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



The following initial and final conditions were used



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

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


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

p20 = 0.18 x 10-5

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

and p20 were taken from the analog computer results

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

may be summarized in the graphs. The graphs also indicate

comparison with the analog computer results, which in

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

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

second case only. Comparison with analog computer results

is the subject of Figure (6).



















c
0

M n














or
C x









- E
Cp




C
--4
+














OI
+'
C







*C







aC
0


if


o r
oaoo0


0
X
o r"

0 C' 0


0
x
00
In '
(0
0


Io
0

So X
S0 r-


o C o




o c
4f) o "0
ocuo


0 0





41




The initial and final conditions used in Figure 7

are:


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

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

Figure 8 shows power level versus time plot of run

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

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

second trial.
















- Quasilineorizotion Data Points

-- Data Points From Analog
Computer Results


0 .02 .04 06 OB 0.1

TIME SECSS)

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






















50 -o- Quasilinearization Data Points

-r- Data Points From Analog
Computer Results


40





S30
-J

nii
.0
o-20





10 -
I0





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

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









50




40 st Trial By Quasilinearization-

Lost Trial By Quasilinearization /




-J


o /
S. 20 0



-29th Iteration By
SGradient Technique On

I0
^'' Anolog Computer





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





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









2. Cost Function 2


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


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

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

c = 6.4n 0.1 c (2.77)

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

Pc =-1000 nP1 (2.80)

with boundary conditions:

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


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


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











where = 1/2 pc2 dt
0


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


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


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

system from the desired value.


Simulation of Problem 2


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

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

Figure 9 Neutron density and reactivity function plots

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

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

that for increasing values of a the rate of control rod

withdrawal is also increased. However, in order to limit

the power level to a certain desired value the

control rod is inserted back into the reactor at a faster

rate, thus shifting the reactivity function peak to the

left.










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

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


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


J2 = 5 x 10-4 J21

where

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

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

that the cost function increases with increasing a This

appears reasonable as increased b would decrease the error

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

The total cost function increases rapidly with an increase

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

from graph of J versus d This result is again quite

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

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































5
4 5


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












0 5 -" 2

4
=5




















50
40
3 a=0--\Q













































time for problem 2 for various values
20




I0





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




of and Nd = 10 no (Nuclear Reactor
System)



















4-







SI ,I I I I

0.2 0.4 06 08 10
TIME SECSS)


50




40






OJ


2








II 0 -i
LlL





0.2 0.4 0.6 0.8 1.0
TIME SECSS)

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











TABLE 4

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



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

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

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

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

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

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

5no 0.696x10-5

6no 1.323x10 -

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

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

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

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

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

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






















1.2









1.0




0.9 -




0.8 -




0.7




0.6



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


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




















































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









3. Cost Function 3


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


pc = a + A(n-nd)

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

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

S= 6.4n 0.lc (2.83)

A = 0. a = 0 (2.84)

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

P2 = 0.1P1 + 0.1 P2 (2.86)

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

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

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

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











Simulation of Problem 3


The equations used for simulation are given in

Appendix A. and the procedure used is similar to that

followed in Problem 1. The circuit diagram is shown in

Appendix A. The trial and error procedure is adopted in

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

obtain the desired power level response. Figure 14 shows

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

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

power levels, the desired control is obtained. The initial

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

varied to satisfy the end conditions on the adjoint vari-

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


It is observed in this investigation and justified by

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

desired power level is reached only for an infinitely

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

also suggests a strong dependence on the desired level.

These observations therefore suggest that the form chosen

for this control will not yield a desirable closed loop

control. The form that accomplishes this is given in

Problem 4.







21

20

19

18

17

16

15

14

13

12




S10

9A




7 o

6

5




3

2





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

















9n,


LUSn

U
ir



-j
LLJ

LUJ
-J
Ir




L'
6n0

4


4no,


0 2 4 6 8 10
VALUE OF A


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









4. Cost Function 4


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

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

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

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

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

c = 6.4n O.lc (2.92)

x = n-nd (2.93)

A = 0 B = 0 (2.94)

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

P2 = -0.1P1 0.1P2 (2.96)

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

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

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

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










The boundary conditions are:

n(o) = no

c(o) = co

x(o) = 0

nl(tf) = nld =. ano


P2(tf) = 0

P3(tf) = 0

P4(tf) = 0

P5(tf) = 0

P4(o) = 0

P (o) = 0

Simulation of Problem 4


(2.101)


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











that such a system would be stable and the response would

be dependent on the coefficients A and B.

The equations used for simulation are given on page

63 and the procedure used is similar to that followed in

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

16 shows the power level and the reactivity function

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

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

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

down to the desired value in fixed time tf. Corresponding

rod withdrawal is also much faster in the initial stages

and then gradual before it reaches a constant value, to

maintain the above level. Figure 17 shows variation

of cost function versus peak power attained. A graph is

also plotted showing behavior of constants A and B with

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

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

B which is independent of the desired power level. This

would imply that the control parameters A and B should

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

sired final state. Hence one would conclude that some sort

of adaptive control may be required. A control system is

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

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

an optimal or nearly optimal fashion.











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

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

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




TABLE 5

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

(NUCLEAR REACTOR SYSTEM)


Nf A B a J

10no 2.123 42.12 3.178 1.201 x 105

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

2no 3.836 128.4 0.105 x 10-5






Nd = 1Ono


Nd = 6 no



Nd= 5no

c-- -----

0 Nd = 2no


Nd = IOno










Nd = 6no o- -o



Nd 5no-

Nd= 5no


50






40




-j
> 30
LIJ
-J

w
_J



0
20






10


.2 .4


.8


1.0


TIME (Sec.)

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


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


Figure 16.


I I I I I ~














13

12




10

0.9

0 0.8

0- .7

0.6

0.5

0.4

0.3

02

01 -


0 2 3 4 5 6 7 8 9 10 11

Nd
Figure 17. Cost function versus peak power
attained












A
4.0




3.5
B
130

120 3.0




100 2.5 -

90 -

80 2.0 -









40 -1.0 -

30

20 .5

10

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

Figure 18. Nd versus A and B









.10
















NdP lOno
S.05 P T.89













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

















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






P 3 P4 T I.Osec

T \.94
S 0

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







.6
ST .90 sec.


.2

0
0

-r T=.9sec.

-2
-.4

-.6 P3

-.8


T= .84sec


.05
P5







.04 T= .O sec

.02 P4 --

-0

Nd = 5n,


T =.95 sec.


T= .8 sec


P5=0, 0 < t sec.


-- Nd = 2no

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











Discussion of Results


Results obtained in the four problems attempted in

this chapter are analyzed and a comparison is given of

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

open loop control system and 2) closed loop control or

specific optimal control system.

The first two problems analysed represent open loop

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

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

ance criterion used here could be interpreted as an opti-

mization of control and movement, limiting erratic flux

changes (10).

Another useful performance index is the subject of

problem 2 wherein the integral of the error squared of

the reactor power and its desired value is minimized along

with the integral of the squared of the reactivity function.

In this problem the control function required to achieve

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

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

the control action and the faster the power rise. This

variation is exhibited in Figure 10. Also it is observed

that the value of the cost function increases with in-

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

since the higher value of a would imply desirability of











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

control is needed to accomplish this.

Problems 3 and 4 represent closed loop or specific

optimal control. The specific forms chosen for the

controller are:


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



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

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

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

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

the cost function increases with increase in the desired

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

obvious.

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

each case that the closed loop or specific optimal control

suffers in performance in the sense that the cost function

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

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

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

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

















































2no 4no 6no

DESIRED Nd


8no IOno


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


1.8 I


1.4

1.2


0.6


0.4


0.2














4.0

[C.L.C.]

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


2.0



1.0 -




0 2no 4no 6no Bno Ion0

DESIRED Nd






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










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

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

pare favorably with those obtained by Eisenberg and Sage

(21) power level trajectories and the corresponding time

variation of the control functions are shown in Figures 23

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

tion versus the peak power desired for varying values of

Figure 26 shows the variation of the cost function with

variation in system parameters (neutron lifetime). It is

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

open loop control for large variation in lifetime about its

normal value.

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

cients in the control law are highly dependent upon the

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

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

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

of the control law is very instrumental in determining how

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

true optimal.

Figure 18 shows the variation of the parameters A

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

power. Since this variation is considerable, adaptive

control is visualized as a means of further improving the

performance. If adaptive control is used, the example at

















SIOn,





8no





6no





4no





2no






0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

TIME SECSS)

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


* 0.3
2













S- lOn, [O.L.C]



10no [C.L.C]




























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








II



4 b




0.9




0.8




07




06




0.5




04




03




0.2




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

Figure 25. Cost function versus peak power desired
























C.L.C]


[O.L.C]


.5 1.0 1.5 2.0
103 x NEUTRON LIFE TIME

Figure 26. Total cost function versus mean neutron
lifetime











hand could be that of the so called performance criterion

sensing or extremum adaptive system. An overall performance

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

law (by trial and error) which optimizes this performance

index. The learning states here are the parameters des-

cribing the control law, and the learning process consists

of the trial and error adjustment of these parameters. If

the learning states can be changed rapidly and accurately

the adaptive controller may be able to handle nonlinear

control systems the learning states will represent the

best instantaneous linear approximation to the non-linear

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

The gradient technique employed yields good re-

sults. The accuracy of the results can be considerably

improved at the expense of more analog and logic elements

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

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

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

sired trajectory. Smaller values result in an increased

number of cycles to converge.

The Quasilinearization technique is next employed

to compute the results already obtained in Problem 1.

Figures 6 and 7 indicate close agreement between the

results obtained by the two techniques.











CHAPTER III


OPTIMAL CONTROL OF A NUCLEAR ROCKET ENGINE


Introduction and Model Description

The nuclear rocket engine to be considered is

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

pound thrust model developed by Smith and Stenning (23,

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

model are based on a lumped parameter analysis. Such an

analysis does not describe the system completely and

furthermore the analysis is affected by a lack of agree-

ment among investigators concerning the temperature re-

activity effects. Stenning and Smith (27 ) contend that

temperature reactivity is directly proportional to the

square root of the core exit stagnation temperature, while

Mohler and Perry ( 2 ) hold that temperature reactivity

is directly proportional to this temperature. Weaver (12)

observes, that both contentions could give fairly accurate

results simply by choosing appropriate reactivity coeffi-

cients. Recent rocket engine tests (28 ) however indicate

that temperature reactivity may be proportional to square

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

nuclear rocket engines possess characteristics which need

exhaustive analysis. Thus, since the system is not well










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

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

closed-loop control is a complex function of the state

variables and due to the analog computer limitations, the

treatment is restricted to open loop control only.

Basic components of the engine are depicted in

Figure 27 Liquid hydrogen, stored in a low pressure

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

tively cool the nozzle and reflector. The propellant gas

is further heated as it passes through the reactor core

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

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

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

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

tem pressure are controlled by the valve in front of the

turbine.

The following equations describe the Smith-Stenning

nuclear rocket engine in terms of the time dependent nor-

malized quantities. An account of the normalization pro-

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

will not be described here.

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


(3.2)





79

























o ~4,
C


















(v 0

mcy
-,






L4
o 4-

-'4
-4- 0






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





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










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

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


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


where N = C' -
Nd Cd

T P


VV V
0' =_ v d

and atTd ahPd
and at = a- BTd


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











TABLE 6


ROCKET ENGINE CONSTANTS


Design reactor power

Design core surface temperature

Design core inlet stagnation pressure

Thermal time constant

Pressure time constant

Control poison time constant

'alve time constant

Decay constant (one group)

Temperature reactivity coefficient

Hydrogen reactivity coefficient

Neutron lifetime

Delayed neutron fraction

Thrust, F

Specific impusle Isp


25,000 MW

5,000R

1.500 psi

0.68 sec

2.50 sec

2.00 sec

0.20 sec

0.10 sec-1

-10-6 oR-1

0.3 OR/psi

10-4 sec

0.005

1,500,000 lb

640 sec


~ _~I___~~


~_










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

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

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

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

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

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

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

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

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

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








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


J = 1/2 tf P2(t)dt


with T(tf) = Td and pl < Pmax

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

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




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs