DESIGN AND IMPLEMENTATION OF MINIMUM TIME COMPUTER CONTROL SCHEMES
FOR STARTUP OF A DOUBLE EFFECT EVAPORATOR
By
SANTOSH NAYAK
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
1974
ACKNOWLEDGMENTS
The author wishes to express his sincere appreciation to
Dr. A. W. Westerberg, chairman of his supervisory committee and
principal research advisor, for his invaluable guidance and en
couragement.
The author also wishes to thank the following people who
have helped him at various times throughout the duration of his
research.
Dr. T. E. Bullock, a member of his supervisory committee,
for his valuable suggestions and for introducing the author to
optimal control theory.
Dr. U. H. Kurzweg who served as a member of his supervisory
committee.
The NorthEast Regional Data Center of the State University
system of Florida and the IBM Customer Engineers who helped out
with the working of the IBM 1070 interface and its communication
with the IBM 370/165.
The Department of Chemical Engineering for financial support,
and the faculty and technicians of the department for their suggestions
and help.
Mr. S. S. Sriram for his help in preparing the figures.
Miss Sara McElroy, the author's fiancee, for her encouragement
during the major part of the work and for her typing of the manuscript.
TABLE OF CONTENTS
ACKNOWLEDGMENTS ............................. .................
LIST OF TABLES .............................................
LIST OF FIGURES ................................................
NOMENCLATURE ...................................................
ABSTRACT .......................................................
CHAPTERS:
I. INTRODUCTION ..........................................
II. DESCRIPTION OF EVAPORATOR AND COMPUTER INTERFACING
EQUIPMENT .............................................
11.1 Evaporator Layout and Description ...............
11.2 Operating Notes .................................
11.3 Evaporator Instrumentation .....................
11.4 Transducing and Controlling Equipment ...........
11.5 IBM 1070 Interface ..............................
11.6 Software
III. DYNAMIC MODEL AND PARAMETER ESTIMATION ................
III.1 Dynamic Model ..................................
III.1.1 State Equations .........................
111.1.2 Connection Equations ....................
111.1.3 Heat Transfer Equations .................
111.1.4 Decision Variables ......................
111.1.5 Assumptions .............................
111.2 Parameter Estimation ...........................
111.2.1 Stochastic versus Deterministic
Estimation ..............................
Page
i
vi
viii
x
xiii
1
5
5
8
9
13
16
17
22
23
23
26
26
29
30
31
31
111.2.2 Experimental Work for Determining 1ea ... 34
111.2.3 Calculations and Results for 1a ........ 40
111.2.4 Experimental Work for Determining elb
and 02 ........................ .... .. 51
111.2.5 Calculations and Results for elb and
2 ............................... ...52
IV. MINIMUM TIME CONTROL POLICY ........................... 63
IV.1 Statement of the Problem for the Evaporator ..... 63
IV.1.1 State and Control variables .............. 63
IV.1.2 State and Control variable Constraints ... 66
IV.1.3 Control Scenarios ........................ 67
IV.1.4 Summary of the Problem Statement ......... 72
IV.2 A Minimum Time Algorithm ........................ 73
IV.2.1 General Problem .......................... 73
IV.2.2 Lagrange Formulation and Necessary
Conditions ............................... 74
IV.2.3 Comments on the Necessary Conditions ..... 79
IV.2.4 Minimum Time Algorithm ................... 81
IV.3 Solution to the Evaporator Problem .............. 85
IV.3.1 Problem 1.' Constraint on the Second
Effect Holdup ........................... 85
IV.3.2 Problem 2. Fixed Feed Rate .............. 102
IV.3.3 Problem 3. No Bound on the Second
Effect Holdup ........................... 110
IV.4 Experimental Runs ............................... 123
FOLDOUT NOMENCLATURE LIST ...................................... 148
V. COMMENTS AND RECOMMENDATIONS .......................... 149
V.1 Model ............................................ 149
V.2 Experimental Setup ............................... 150
V.3 Theory ......................................... 152
V.4 Conclusions .................................... 152
APPENDICES:
A. HEAT TRANSFER EQUATIONS AND OTHER RELATIONSHIPS ..... 154
A.1 Relation between Temperatures and Enthalpies ... 154
A.2 Heat Transfer EquationsFirst Effect .......... 154
A.2.1 Sensible Heating Zone ................... 155
A.2.2 Vaporizing Zone ......................... 157
A.3 Heat Transfer EquationsSecond Effect ......... 160
B. LISTING OF COMPUTER PROGRAMS ........................ 162
LITERATURE CITED ........................................... 206
BIOGRAPHICAL SKETCH ......................................... 208
LIST OF TABLES
Table
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
1 Data for Run A1..........................
1 Data for Run A2..........................
1 Data for Run A3..........................
1 Data for Run A4..........................
1 Data for Run A5..........................
versus Observed Values of T1 for Run Al....
versus Observed Values of T1 for Run A2....
versus Observed Values of T1 for Run A3....
versus Observed Values of T1 for Run A4....
versus Observed Values of T1 for Run A5....
1 Data for Run B1..........................
1 Data for Run B2..........................
1 Data for Run B3..........................
3.14 Calculated versus Observed Values of T, and T2 for
Run B1 ..............................................
3.15 Calculated versus Observed Values of T1 and T2 for
Run B2 ................................ ..............
3.16 Calculated versus Observed Values of T1 and T2 for
Run B3 ................................ ..............
4.1 State Variables for Problem 1, Iteration 1...........
4.2 Adjoint Variables for Problem 1, Iteration 1..........
4.3 State Variables for Problem 1, Iteration 2...........
4.4 Adjoint Variables for Problem 1, Iteration 2..........
Experimenta
Experimenta
Experimenta
Experimenta
Experiment
Calculated
Calculated
Calculated
Calculated
Calculated
Experimenta
Experimenta
Experimenta
Page
35
36
37
38
39
46
47
48
49
50
53
54
55
59
60
61
92
93
94
95
Table Page
4.5 State Variables for Problem 1, Iteration 3............ 96
4.6 Adjoint Variables for Problem 1, Iteration 3.......... 97
4.7 State Variables including Concentration Dynamics...... 100
4.8 State Variables for Problem 2, Iteration 1............ 107
4.9 State Variables for Problem 2, Iteration 4............ 108
4.10 State Variables for Problem 2, Iteration 5........... 109
4.11 Adjoint Variables for Problem 3, Iteration 1.......... 117
4.12 State Variables for Problem 3, Iteration 3............ 118
4.13 Adjoint Variables for Problem 3, Iteration 3.......... 119
4.14 State Variables for Problem 3, Iteration 5............ 120
4.15 Adjoint Variables for Problem 3, Iteration 5.......... 121
4.16 Experimental Data for Run C1 .......................... 130
4.17 Theoretical Minimum Time Simulation for Run Cl........ 131
4.18 Actual Minimum Time Simulation for Run C1............. 132
4.19 Experimental Data for Run C2........................... 134
4.20 Actual Minimum Time Simulation for Run C2............. 135
4.21 Experimental Data for Run C3.......................... 137
4.22 Actual Minimum Time Simulation for Run C3............. 138
4.23 Experimental Data for Run C4.......................... 140
4.24 Actual Minimum Time Simulation for Run C4............. 141
4.25 Experimental Data for Run C5.......................... 143
4.26 Actual Minimum Time Simulation for Run C5............. 144
4.27 Experimental Data for Run C6.......................... 146
4.28 Actual Minimum Time Simulation for Run C6............. 147
LIST OF FIGURES
Figure Page
2.1 Schematic Diagram of the Double Effect Evaporator..... 7
2.2 Evaporator Instrumentation.. ......................... 11
2.3 Evaporator Instrumentation.......................... 12
2.4 Layout of Transducer and Controller Cabinet and IBM
1070 Cabinets......................................... 15
2.5 Process Interface Computer Information Flow....... 18
2.6 Software Setup........................................ 20
3.1 Variables for Material and Energy Balances............ 24
3.2 Calculated versus Observed Values of T1 for Run Al.... 41
3.3 Calculated versus Observed Values of T1 for Run A2.... 42
3.4 Calculated versus Observed Values of T1 for Run A3.... 43
3.5 Calculated versus Observed Values of T1 for Run A4.... 44
3.6 Calculated versus Observed Values of T1 for Run A5.... 45
3.7 Calculated versus Observed Values of T1 and T2 for
Run B1 .............................................. 56
3.8 Calculated versus Observed Values of T1 and T2 for
Run B2 ............................................... 57
3.9 Calculated versus Observed Values of T1 and T2 for
Run B3 ............................................... 58
4.1 Control, State and Adjoint Variables for Problem 1,
Iteration 1 ......................................... 89
4.2 Control, State and Adjoint Variables for Problem 1,
Iteration 2 .......................................... 90
4.3 Control, State and Adjoint Variables for Problem 1,
Iteration 3 .......................................... 91
4.4 Optimal Simulation including Concentration Dynamics... 99
Figure Page
4.5 State Variables for Problem 2, Iterations 1, 4 and 5.. 106
4.6 Control, State and Adjoint Variables for Problem 3,
Iteration 1........ .................................. 114
4.7 Control, State and Adjoint Variables for Problem 3,
Iteration 3 .......................................... 115
4.8 Control, State and Adjoint Variables for Problem 3,
Iteration 5 ......................................... 116
4.9 Filtered versus Actual Flow Rate...................... 127
4.10 Experimental versus Actual Minimum Time for Run C1.... 128
4.11 Experimental versus Optimal Minimum Time for Run Cl... 129
4.12 Experimental versus Actual Minimum Time for Run C2.... 133
4.13 Experimental versus Actual Minimum Time for Run C3.... 136
4.14 Experimental versus Actual Minimum Time for Run C4.... 139
4.15 Experimental versus Actual Minimum Time for Run C5.... 142
4.16 Experimental versus Actual Minimum Time for Run C6.... 145
NOMENCLATURE
A = Heat transfer area, ft2
C = Solute concentration, Ibs/(lb solution)
C = Specific heat, Btu/lb
D = Diameter, ft
G = Mass velocity, lbs/(ft2)(hr)
Gr = Grashof number, dimensionless
H = Holdup, Ibs
I = Index set
K = Index set
L = Length, ft
N = Number of tubes
P = Pressure, lbs/(in2)
P(t) = Covariance of estimate, vector
Pr = Prandtl number, dimensionless
Q = Heat transfer rate, Btu/(hr)
Q(t) = Covariance of process noise, vector
R(t) = Covariance of measurement noise, vector
Re = Reynold's number, dimensionless
T = Temperature, F
AT = Temperature difference, F
U = Overall heat transfer coefficient, Btu/(hr)(ft2)(OF)
U(t) = Control vector
V(t) = Process noise, vector
V' = Vapor volume, ft3
V = Vapor flow rate, Ibs/min
W = Liquid flow rate, Ibs/min
X = State variables, vector
X = Estimate of state variables, vector
Xt = LockhartMartinelli factor
Y = Calculated observations, vector
Y = Actual observations, vector
f = Function of
g = Acceleration due to gravity, ft/(hr2)
h = Liquid enthalpy, Btu/lb
hv = Vapor enthalpy, Btu/lb
h = Film coefficient, (Btu)(ft)/(hr)(ft2)(OF)
k = Thermal conductivity, (Btu)(ft)/(hr)(ft2)(F)
p = Variance of estimate
t = Time, minutes
u = Control variable
v = Process noise
w = Measurement noise
x = State variable
x = Estimate of state variable
Subscripts
o = Outside world
1 = First effect
2 = Second effect
i j = From unit j to unit i
a = Before
c = Condensate
s = Sensible heating zone
t = Tube
B = Boiling zone
w = Wall conditions
F = Feed
i/j = At "i" given conditions at "j"
f = Film conditions or final condition
Superscripts
i = Inside
o = Outside
v = Vapor
* = Optimality
Greek Letters
a = Point constraint multiplier
B = Inequality constraint multiplier
o = Estimated parameters, vector
6 = Estimated parameter
p = Density, lbs/(ft3)
A = Lagrange multiplier or latent heat of vaporization,
Btu/lb
= Viscosity, Ibs/(ft)(hr)
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy in Chemical Engineering
DESIGN AND IMPLEMENTATION OF MINIMUM TIME COMPUTER CONTROL SCHEMES
FOR STARTUP OF A DOUBLE EFFECT EVAPORATOR
By
Santosh Nayak
March, 1974
Chairman: Dr. Arthur W. Westerberg
Major Department: Chemical Engineering
The application of optimal control theory to a chemical engi
neering problem is investigated by the development and implementation
of a control policy for the minimum time startup of a double effect
evaporator.
The particular evaporator on which experimental runs were made
was a laboratory scale double effect evaporator with reverse feed. It
was completely instrumented for control by the installation of orifices
for measuring flow rates, thermocouples for measuring temperatures,
pressure taps for measuring pressures and holdups, and pneumatic con
trol valves for manipulating flow rates. Transducing and controlling
instruments were installed. In order to do online computerized data
logging and control, interfacing of the process with the IBM 370/165
computer on campus was provided by an existing IBM 1070 interface. A
part of the existing user circuitry associated with the interface had
to be rewired and modified to function appropriately.
The dynamic model of the evaporator consisted of six differential
or state equations and about sixty algebraic equations. This latter
group consisted of connection equations between the effects, property
correlations and heat transfer equations. To overcome the uncertainty
in the empirical relationships for the inside film coefficient two
unknown parameters were introduced, one for each effect. These para
meters were estimated by correlating model predictions with data
collected on experimental runs. A nonlinear leastsquares technique
was utilized to get the best fit.
The algorithm used for the theoretical development of a minimum
time policy is one in which Hamiltonian minimizations result in control
policy updates on successive iterations. Control variable constraints
and point constraints are accounted for along the trajectory. The
utility of the algorithm is enhanced by assuming a control scenario and
determining whether it is optimal when compared to other likely scenarios
This approach keeps the number of active state and adjoint variables to
a minimum at any particular time resulting in simple Hamiltonians and
less computational expense for the integration of the state and adjoint
equations and Hamiltonian minimizations.
The minimum time algorithm was used on the evaporator model to
determine the optimal policy under three different sets of conditions.
In the first case it was assumed that there was a constraint on the
maximum value of the second effect holdup. The second case dealt with
a different set of control variables in that the feed rate to the second
effect was assumed to be fixed. In the third case the assumption of a
constrained maximum second effect holdup was done away with. The
simulation results indicated that the third case resulted in the smallest
startup time with the optimal policy calling for an overfilling of the
second effect followed by a gradual decrease in the second effect holdup
to the desired value which took place when boiling just started in
the second effect.
For all three cases it was found that the control policy is
bangbang in nature and that the control switches occur at times at
which the point constraints are met on the assumed scenario. Because
of this the switching times can be related to the state variables
and a feedback control policy is obtained.
Experiments were run to try out the optimal control policy and
to test the model. On an average, the simulations resulted in final
times which were between ten and fifteen percent within those obtained
experimentally. This accuracy was reasonable considering the experi
mental problems associated with holdup measurements and analog control
of the holdups and the theoretical problems associated with the as
sumption of the heat transfer mechanisms.
CHAPTER I
INTRODUCTION
Optimal control theory has been developed to a fairly sophis
ticated level in the field of electrical engineering. However, the
uses of the theory and possible applications in chemical engineering
have been virtually unexplored. The reasons for this rather limited
progress on both the theoretical and applied fronts are many (Foss,
1973).
The starting point in the applications of control theory is
a good dynamic model of the process. Most chemical processes have
been modeled poorly due to an incomplete understanding of the complex
interactions among numerous variables. High dimensionality and non
linearities in behavior require the use of sophisticated numerical
techniques for the simulation and design of control schemes. Many
chemical processes have inherently large time constants which make
them unsuitable for control. It is generally not possible to make
all the measurements that are required for a feedback control scheme,
and, in addition, measurements are subject to noise which is not
readily filtered in nonlinear systems. Since most of the states are
not measurable, there is a need for building reduced order observers
which again is a formidable problem for nonlinear systems.
In spite of the above drawbacks, quite a few articles on the
subject of optimal control have appeared in the chemical engineering
literature during the last five years. The early investigations dealt
with the control of simplified lumped parameter linear processes. The
linear quadratic loss problem resulted in feedback control which was
particularly useful in regulatory control; i.e. control in the face of
disturbances (Nieman and Fisher, 1970; Newell and Fisher, 1971; Newell
et al., 1972). The arbitrary nature of the correlation of the weighting
matrices to the actual physical problem often makes the linear quadratic
loss problem unrealistic. Later investigators extended these techniques
to the control of nonlinear lumped parameter systems by using various
forms of linearization on the system equations (Weber and Lapidus, 1971a,
1971b; Siebenthal and Aris, 1964; Tsang and Luus, 1973). Others worked
on nonlinear systems with one or two control and state variables (Joffe
and Sargent, 1971; Jackson, 1966). The control of distributed parameter
systems is still in its infancy. Some investigators have reported sub
optimal control of distributed systems in which some other criterion,
such as minimization of a Lyapunov functional, is used (Vermeychuk and
Lapidus, 1973a, 1973b; Chant and Luus, 1968). Simulated startup studies
have been made on plate distillation columns (Pollard and Sargent, 1966)
and autothermic reaction systems (Jackson, 1966).
This work was directed towards applying existing control theory
to a useful chemical process. The aim was to study the dynamics and
control of a simple, yet reasonably complex, piece of equipment commonly
found in the chemical industry. A double effect evaporator was chosen
as the subject of study for the above reason and also because a laboratory
scale double effect evaporator was available for experimental work. The
study involved,
a) Developing a nonlinear model for the evaporator.
b) Estimating model parameters to fit experimental data.
c) Developing a minimum startup time control policy taking into
account constraints on the state and control variables and putting
the optimal policy in feedback form in terms of switching times.
d) Experimentally determining the effect of the policy.
This approach differs from previous ones in a few respects.
The model is highly nonlinear and is treated as such. No linearization
is resorted to as startup involves large changes in the state variables
and linearized equations would be inaccurate. The mathematics involved
in obtaining the minimum time policy is simplified as the approach
adopted presupposes a startup scenario and then verifies that it is
optimal. The algorithm leading to the optimal policy handles con
straints on control and state variables in a logical fashion by directly
holding the state or control variable on the constraint and changing
the equation set and its solution procedure as a result. This avoids
the use of penalty function methods and the like. Finally, the control
policy is experimentally verified.
The minimum time objective was chosen primarily from the point
of view of economics. Control costs during startup are minimal com
pared to the startup time in batch processes. Reducing the startup
time results in reduced down time thus improving cycle efficiency and
increasing profits. The food industry is an example of an industry
which must shut down frequently to have the processing equipment cleaned.
Orange juice is concentrated in multiple effect evaporator systems, and
these systems are cleaned about three times a day. A second reason for
minimum time startup was more specific to the ultimate use of the parti
cular double effect evaporator investigated. It is to be used in an
undergraduate laboratory experiment in computer control, and past
experiences indicated that it took a very long time to bring it to
steady state under manual control. Thus, in order to reduce the
startup time and consequently to reduce the amount of online com
puter time for steady state observations, it was imperative to have
startup in a minimum time.
Chapter II contains a description of the experimental evapo
rator, the instrumentation and the interfacing equipment with the
IBM 370/165 (which is the main computer on campus). Chapter III deals
with the building of a dynamic model for the evaporator and also the
estimation of parameters from experimental data. A derivation of the
optimal control algorithm is given in Chapter IV. It also contains the
simulated and experimental results of the application of the control
algorithm. Some comments and proposals for further work are given in
Chapter V. Appendix A contains all the heat transfer equations which
supplement the main model equations in Chapter III. A listing and
description of the computer program implementing the optimal control
algorithm is the subject of Appendix B.
CHAPTER II
DESCRIPTION OF EVAPORATOR AND COMPUTER INTERFACING EQUIPMENT
II.1 Evaporator Layout and Description
The double effect evaporator is located in the unit operations
laboratory of the chemical engineering department. Figure 2.1 is a
schematic of the double effect, showing the arrangement of the two
effects, EV1 and EV2, and the basic process and vapor lines. Note
that backward feed is used; that is, the vapor flow and process fluid
flow are in opposite directions.
The first effect is a long tube vertical (LTV) evaporator.
It contains 3 tubes, each 9 feet, 6 inches (2.90 m) long and 1 inch
(0.0254 m) O.D. with heating steam at about 20 psig (2.39 bar) on the
outside of the tubes. The process fluid flows upward through the
tubes either by natural or forced circulation. The latter method is
almost always used because of the increased heat transfer coefficients
obtainable. The pressure on the process side is at or slightly above
atmospheric.
The mixture of process fluid and vapor formed in the first
effect enters a vaporliquid separator, SE1, which is at the same
pressure as the first effect. The liquid is drawn off the bottom of
the separator and is recirculated back into the first effect by pump
PU3 after some liquid product is withdrawn. Fresh feed to the first
4)
S
O
'
4
f0
LJ
S
S
a
Q.
41 r
uo
4
4
o ci
4*
4 C U Q
>l 4
LL(J
=:3
C)
E
ci
u
U Yu: >
LU
CZ C
V) UJ 
Z)
0 *
effect is pumped by pump PU2 from the second effect. The vapor from
separator SE1 is used as the steam input in the second effect. This
leads to steam economy as one pound of heating steam used in the first
effect should evaporate more than one pound of water from the first
and second effects combined.
The second effect, EV2, is a calandria type effect in which
there are 15 tubes, each 2 feet, 4 inches (0.711 m) long and 1 inch
(0.0254 m) O.D. The effect also has a 2 inch (0.0508 m) O.D. central
downtake. Heat transfer is by natural convection only, resulting in
much lower heat transfer coefficients compared to the first effect.
Fresh preheated feed is pumped into the bottom of the second effect
from the feed tank by pump PU1. The heating medium is the vapor from
the first effect on the outside of the tubes. Above the calandria is
a vapor body which separates the vapor from the liquid. The vapor is
drawn into a condenser, CDI, by means of a vacuum produced by a steam
jet ejector. The ejector maintains the pressure on the process side
in the effect at around 10 inches mercury vacuum (0.675 bar).
The vapor condensate from the first effect is collected in tank
T1 and that from the second effect is collected in tank T2, both of
which are maintained at a vacuum by the same steam jet ejector.
11.2 Operating Notes
There were a few precautions which had to be observed during
operation.
1) The feed rate was kept at around 23 gpm (0.12 to 0.18 kg/s).
2) The recirculation rate in the first effect was kept at a maximum
of 15 gpm (0.9 kg/s). A higher rate caused entrainment of liquid with
the vapor in the vaporliquid separator SE1. This separator has no
baffling of any kind and is very inefficient at high flow rates.
3) To avoid cavitation in the recirculation pump PU3, care was taken
to see that the vertical suction leg from the separator to the re
circulation pump was always filled with liquid. This was particularly
critical when the pump was first started. Incomplete filling of the
vertical leg led to pulsating flows resulting in large upsets in the
evaporator operation. A recirculation rate higher than 15 gpm (0.9
kg/s) also caused a high discharge head on pump PU3, much higher than
the maximum discharge head on pump PU2 (which is of a smaller capacity),
eliminating all flow of fresh feed to the first effect.
4) The liquid level in the second effect was maintained around the
top of the tubes for best utilization of the heat transfer area.
11.3 Evaporator Instrumentation
As part of this work the evaporator piping had to be modified
to accommodate the instrumentation required for control. The work con
sisted mainly of installing pneumatic control valves, orifices, pressure
taps, thermocouples and extra manual valves. Figures 2.2 and 2.3 show
the detailed instrumentation of the evaporator. The legend of Figure
2.2 also applies to Figure 2.3.
Three pneumatic control valves, CV1, CV2 and CV3, were installed
in the feed, intereffect and recirculation lines respectively. CV2
and CV3 were normally closed (airtoopen) valves and were installed in
bypasses on the lines, whereas CV1 was a normally open valve and was in
stalled in the feed line as such. The purpose of the bypasses was to
.4
3 00)
0 L
cu1 4. cr*
a LC
*r 0 41 C 1
*o 0 C 0 S* 0
CL 0 0 0 Q 0 0 0
3 0 a) 4 4) 4)
L, U L C ) V CO 0
) X > U F 1 D C 
= LU L) () L" LL F
0
3 4 0)
0) >
S) r
0)0)00)
U + S 1 
*o0 *O 0 C
> VI .
0L iIL U
00) W
4 41C V4 CM 41
I C C II II
S0 E E 0)
r 9C
VC E 0
Ii
E
C C:
Q: (A
C
>
C.
o
0
44
s n
o 0
_) S
.0
E 0)
*r 5
0) 4)
S f0
3 SL
LI 0)
0 H
4)
0 CC C C C
C  t 05
 o. H U. C H
41
0 )
0) 0)
U +4
X 0)
) >
= LU
30
C 
0 LI
0LI
0)
L 0
ES.
0
0)
> E

> +
L>C
0* )
>4)
U
10
0 )
C 0
10
0)
>
 c
> >
> s
V' 1
SI
Y
{29 <2^ 4
A ^ 8:s
~i 1
I___
~1
allow for complete or partial manual control of experimental runs
when desired.
Flow rates were measured with squareedged orifices, OR1,
OR2, OR3 and OR4, installed in the feed, intereffect, recirculation
and product lines respectively. The pressure drop across an orifice
is indicative of the flow through it.
Liquid levels (proportional to holdups) were measured by
taking the difference in total pressure between the bottom and top
of each of the two effects. Pressure taps were installed at both
ends of the sight glasses for this purpose. The upper taps were
also used to measure absolute pressure in the effects.
Temperatures were measured by jacketed copperconstantan
thermocouples, TC1, TC2, TC3 and TC4. These were installed in the
feed line, at the exit of the second effect, at the entrance to the
first effect, and in the steam chest of the first effect, severally.
All liquid lines from the pressure taps and air lines to the
valves were brought to a central panel in the front of the evaporator
with 1/4inch polyflo tubing. Quickconnect fittings were used at
the panel so that leads to the interfacing equipment could be con
nected quickly when required. The thermocouple wires also terminated
with special thermocouple outlets at the panel.
11.4 Transducing and Controlling Equipment
The transducing and controlling equipment was installed in a
19inch relay rack on casters. All air and liquid lines were of poly
flo tubing with quickconnect fittings. This made the rack very ver
satile as it can be moved to a number of different pieces of equipment
if desired. A layout of the cabinet is shown in Figure 2.4.
The pneumatic controllers have adjustable proportional and
reset action, motorized set point control, and indication facilities.
One Fischer and Porter model 51 and three Taylor model 662R control
lers were installed. The set point motor of the Fischer and Porter
controller operates on a pulse train input and Taylor controllers on
a 24volt DC signal. The controllers are equipped with feedback
potentiometers which indicate their set point positions. The pneu
matic input signal range to each controller (from the DP cells) is
315 psig and that of the output pressure to the associated valve
is also a 315 psig signal.
The EMF to pneumatic converters (not used in the current ex
periments) are Foxboro model 33A converters. They can transduce
either a millivoltage or voltage signal into a 315 psig pneumatic
signal.
The differential pressure (DP) cells used are Foxboro Model
13A DP cells. The adjustable range of the differential input signal
is 0500 inches water and the proportional pneumatic output is in the
315 psig range. These DP cells were used to transduce the pressure
drops across the orifices and the pressure differences corresponding
to the liquid levels. The outputs were thus proportional to the flow
rates or to the liquid levels (holdups). The DP cells were also
used for measuring absolute pressures by venting the high pressure
side when measuring pressures above atmospheric. The output pressure
in this case was proportional to the vacuum or aboveatmospheric
pressure.
 z w )I
I 1 1 U) 
Q F LLl F L C UV w CD
LU c 2 < C ) C)
11. C)I Z ) I Q C Q lC'
SCD o
UJ l < 0 oQ
.a:: C, V) F U) Im
= F 0
01 Z I I 0 QT I U
F o 0 C' >< rJ O
:K LU CD : L3 C
I __ a _ F 
F
0
LU .
J 2: I I I Z
mLU 0j 0 
) C > 0 C < O <
C 0 C F C) C
0 I l.  . Z u n
e a a 0
 (J
0. 0
U O U) . l. l L U0
i F F 0 a 0 0
U)l L0 l  V )
lu 1 2 I C F
0 F C)d L U
J LU
LU ) r J U)
a )
I CE __1 __U
LU LU I I
of 0 F 0 J (J Q)
F 0 t 0 
C) 0J Q 
L)
1+
0
41
=3
0
o
U
zl
= _C
LU )
C) C)
00
LU LU
F
t =
11.5 IBM 1070 Interface
An IBM 1070 interface was used which consisted of an IBM 1071
central terminal unit, two IBM 1072 multiplexing units, an IBM 1073
Model 3 digital to pulse converter and an IBM 1075 decimal display.
The interface with all the auxiliary equipment resides in 3 relay racks
as shown in Figure 2.4. These cabinets were originally assembled by
R. C. Eschenbacher and are described in his Ph.D. thesis (Eschenbacher,
1970). However, as part of this work, some of the equipment had to be
rewired to accommodate pulse duration outputs and most of the relays
had to be rewired and changed to double coil relays in order to isolate
the IBM equipment from the user circuitry. The entire set up is de
scribed in detail in the GIPSI (General Interface for Process Systems
Instrumentation) hardware manual (GIPSI, 1973) and is summarized very
briefly here.
1) Input
The pressurevoltage transducers converted the 315 psig air
signals from the DP cells to 05 volt DC signals which were fed as
analog inputs to the 1070. The millivoltage thermocouple signals
were also fed as analog inputs through the special thermocouple in
put feature in the 1070. The unit also has a facility for digital
input which was not used.
Another convenient input facility frequently used was a form
of digital input through predesignated demand functions which were
dialled into the system through rotary switches.
2) Output
Output from the 1070 was in digital and in pulse form. The
digital output was used mainly to ring a bell to alert the operator
to possible alarm conditions in the hardware and software. The
pulsed output was obtained from the 1073 and was used to move the set
points on the controllers. The digital to pulse converter (1073)
has outputs in the form of a pulse train as well as a duration pulse.
The 1075 decimal display is a feature which utilizes digital
outputs and was used to display particular variable values or error
codes.
3) Process Alerts
The 1070 interface is linked with the IBM 370/165 computer on
campus. The seven process alerts attached to a process alert (PA)
bus in the 1070 provide a hardware interrupt capability of the com
puter by the process. The software issues a conditional read of the
1070 terminal to the computer. It is then in a hardware wait stage.
When the PA bus is activated, by one of the process alerts on the
1070, the IBM 370/165 computer senses the closure and reactivates
the software which then determines which PA was set. The software
then resets the PA and executes the program associated with the PA.
Process alert 1 has the added facility of being set automatically by
a hardware poller on which the timing is adjustable. This enables
one to have PA1 periodically and automatically set after a predeter
mined time interval has elapsed.
Figure 2.5 is a schematic illustrating the flow of information
among the various hardware components of the experiment.
11.6 Software
The GIPSI software package was originally written by R. C.
Eschenbacher. Version 2, the version used, was written by L. A. Delgado
 LU
0:
PI
F
0
LO
c.
2C
O
CD
C>
LU
1 I
r3 ^
Q 0
41
1 0 3
c 0
0 0
CI *r C3
=3 a) a a
03 3 i: 0
43 1 CL. 0 4
4  ,u
S0 01 W ) 03
O 10 3 0 * E
Sro C C
0u a 0 10 r3 U) 0
U ELL'V ^
3 0 r
0 C > > 0 S
'
=3
C)
00
C
.) QW
4'
3 I
0
U)
w u)
0
n o
.' 3 S
:3 > *. 0o 2O 0
0W CC
C. S 0'V 0 0
I Q L) = m (L)
U C I :
E ) ) C) C.
Wooarr 0
 . > c: C 0
(GIPSI, 1973) and was marginally modified as part of this work to
extend its capabilities. It is written entirely in Fortran (with
the exception of certain input/output routines in BTAM provided by
IBM), has an extensive debugging facility, and an extensive error
handling facility to flag user software and hardware errors. A
simplified flow chart is reproduced from Westerberg and Eschenbacher
(1971) and is shown in Figure 2.6. It is described in greater detail
in Westerberg and Eschenbacher (1971) and GIPSI (1973).
The heart of the software is the concept of the execute and
delay stacks. When a PA is set, it is identified and the program
(or programs) associated with it are stacked by the program stacker
in the order of priority on the execute stack. Control then passes
to the Execute subprogram which then examines the execute stack and
passes control one at a time to the programs that are due for execu
tion. If the sequence began with PAl, its response program CLOCK
removes programs from the delay stack if their delay time has expired
and puts them on the execute stack. The Execute subprogram then finds
additional programs on the execute stack which it continues to remove
and cause to be executed. This is done until the execute stack is
empty whence control returns to the PA handler which issues a condi
tional read and the IBM 370/165 again waits for a PA to be set to
start the cycle again. All delay times are compared to the computer
clock. Data on program priorities, delay time, etc. are specified
in the Program Descriptive Data.
All the user has to do is to provide his specific programs
associated with the various process alerts, the program descriptive
data for all his programs, and a subroutine GOTO which the Execute
II
II
L~II
><  )
bj I
auJ
I , \
w ~
(I)I
21
subprogram uses to pass control to CLOCK and the user and system
subprograms.
The computer costs are extremely low when based mainly on
central processing unit (CPU) time. The software utilizes very
little CPU time. Typical costs are in the range of $35 an hour
provided no elaborate computations are called for in the user pro
grams. However, costs for core residency charges dominate as the
basic software package requires around 20,000 words or 80,000
bytes of core.
CHAPTER III
DYNAMIC MODEL AND PARAMETER ESTIMATION
The dynamic modeling of multipleeffect evaporators has been
extensively investigated in recent years at the University of Alberta
(Andre and Ritter, 1968), (Newell, 1970). In simulation and experi
mental work high order, linear models have been found to be satisfactory.
However, linear models are not realistic when the operating conditions
change drastically as in startup. In the first part of this chapter
a nonlinear, first order, lumped parameter model is proposed. The
first order and lumped parameter nature of the model was resorted to
for two main reasons:
1) The model was simple and adequately described the data.
2) The model was used to devise an optimal control policy for minimum
time startup. Optimal control theory has been rigorously developed
for lumped parameter systems and its extension to distributed systems
has not yet been extensively investigated.
In addition, the model presented here takes into account heat
transfer dynamics from the viewpoint of film coefficients. Although
this leads to complicated algebraic equations, it has the advantage of
leading to a better understanding of the heat transfer dynamics. It
also gives rise to two constant correction parameters. The necessity
for these parameters is due to the uncertain coefficients that are used
in the film coefficient equations. The second part of this chapter
deals with the estimation of these parameters to fit the experimental data.
III.1 Dynamic Model
The dynamic model is a collection of the material and energy
balances for each effect. For a double effect evaporator concen
trating a solution with one major solute, there are two material
balances (one for the solution and one for the solute) and one energy
balance for each effect, giving rise to a total of six dynamic or
state equations for the two effects. In addition, there are dynamic
equations for the vapor phases and metal but the time constants of
these are negligible compared to the six mentioned earlier (Andre,
1968) so that these dynamic equations could be reduced to be alge
braic equations. This procedure of setting the derivatives of the
equations with small time constants to zero reduces the order of
the system. The full model will be presented here. In later chapters,
appropriate simplifications will be applied as some of the model
states are held fixed (for example, as boiling does or does not take
place). A summary of all the assumptions made is presented at the
end of the model. Refer to Figure 3.1 for the symbols used for the
flows, holdups and temperatures. There is also a foldout nomenclature
list on page 148.
III.1.1 State Equations
1) First and second effect holdups, H1 and H2.
dH1
dt 12 V21 1 l 301
dHl
dt =IF V2 2 (3.2)
2) First and second effect enthalpies.
L)
u
ra
cr
a
U.
ct
L
0
r
r
C
VI
Since the evaporator was to be used ultimately for the con
centration of dilute solutions, it was assumed that there would be
no boiling point elevations in the effects. Also, perfect mixing
is assumed which would be close to the case for small holdups and
dilute solutions.
d(Hlh)
dt = L'12h2 V21h (1 11+ l01)h1 + 1
where Q1 is the heat transferred from the steam. Simplifying this
using equation (3.1), we get
dh1 1
dt I 12(h2 h1) + (h] (3.3)
Similarly, an energy balance on the second effect under the
same assumptions gives rise to
dh2
[W2 F (hF h 2(h h 2 (3.4)
at h2) + V02h2 2+Q2
3) First and second effect solute material balances.
Again, assuming perfect mixing we have for the first effect
solute,
d(H1C1)
dt 12C12 (11 + "01)C1
Simplifying this with equation (3.1) we have
dC I
dC1 1
dt 1 1[2(C12 C1) + V21C1] (3.5)
Similarly, a balance on the second effect solute yields,
dC2 1
t = T [ F(CF C2) + V02C2] (3.6)
111.1.2 Connection Equations
In addition to the state equations listed above, there are
algebraic equations which arise due to the mixing of the two streams
between the second and first effects. One energy and two material
balances describe the mixing as follows:
W12 = '2 + W11 (3.7)
W12C12 = 12C2 + "UIIC1 (3.8)
1 12h2 = 2h2 + lllhl (3.9)
111.1.3 Heat Transfer Equations
These equations arise in computing the terms Q1 and Q2 that
arise in the enthalpy equations (3.3) and (3.4).
The heat transfer rate Q1 is a function of the steam temper
ature, the first effect temperature, the inside and outside film
coefficients, the wall resistance including fouling and the heat
transfer area in the first effect. The inside film coefficient is
a function of the flow rate through the tubes, the vaporization,
the inner wall and bulk temperatures and the entrance temperature.
The heat transfer mechanism is initially simplea combination of
the DittusBoelter equation for the inside and the Nusselt equation
for the condensing steam. However, when boiling takes place two
phase heat transfer occurs because of the vapor formed. The complete
boiling mechanism is a topic for further investigation. Approximate
correlations were obtained from (Fair, 1960, 1963a, 1963b), (Hughmark,
1969) and more recently from (Tong, 1965). The complete list of
equations leading to the determination of 01 from the state variables
and flow rates is given in Appendix A. Due to the uncertainty in
the empirical equations which predict the inside film coefficients,
a parameter 01 was introduced in the overall heat transfer equations
(A.17) and (A.40). It is assumed that the outside film coefficient
is predicted by the Nusselt equation, (A.14) and (A.28) to a rea
sonably high degree of accuracy as is borne out later by experiment.
It is also assumed that the parameter 91 has two different values
depending upon the heat transfer mechanism in the first effect. This
depends upon the stage of startup as follows:
1) e1 = Ola, when the first effect liquid is being heated. The
dominant equation for the inside film coefficient is solely the Dittus
Boelter equation (A.16).
2) 61 = alb, when the liquid in the first effect is boiling. The
inside film coefficient is a combination of many factors including a
coefficient due to nucleate boiling (A.38) and a twophase convective
coefficient (A.33).
The equation for the first effect heat transfer rate is given
in functional form as:
Q1= Q (Ts'T1T 12,'H W12',V21' 1) (3.10)
where the temperatures, T, are functions of the enthalpies, h, in the
form
T = f(h)
Equation (3.10) has implicitly used the fact that the overall
coefficient, film coefficients and heat transfer area are functions of
temperatures (enthalpies) and holdups.
Strictly speaking, the equations describing the steam (vapor)
temperatures or enthalpies, hs, hI and h2 are differential equations,
dp
1 21 c2
v v
d(p hl) v
Vi E = V21h1 Wc2hc2 2
where Vi is the volume of the vapor space in the first effect, the
vaporliquid separator and the tubes of the second effect. Note that
two assumptions have been made herethe steam (vapor) is saturated
and there is no subcooling of the condensate.
However, it has been shown by Andre and Ritter (1968) that
the response rate of the steam enthalpy is negligible compared to
that of the holdup, concentration and enthalpy equations (3.1) to
(3.6) in the two effects. The differential equations describing
the steam density and temperature are so replaced with the steady
state equations
V21 =c2
and = V21(h hc2
or Q2 = V21 (3.11)
where = f(T)
The heat transfer rate in the second effect, Q2, is also a
function of the film coefficients, wall resistance including fouling,
area and temperatures in the second effect. The heat transfer
mechanism is purely natural convection. Here again, it is assumed
that the Nusselt equation (A.47) is reasonably accurate in predicting
the outside film coefficient. The inside film coefficient is pre
dicted by the natural convection equation (A.51) and the overall
coefficient (A.52) has an undetermined parameter 82 which again can
29
have two values. One value (e2a) is for the heating of the liquid
and the other (e2b) is for the boiling of the liquid in the second
effect. The functional form for Q2 is:
Q2 = Q2(TIT2TFH2'82) (3.12)
where the temperatures have been determined from the corresponding
enthalpies.
Note that there are no liquid flow and vapor flow terms here
as the natural convection overall coefficient is not a function of
these variables.
A subroutine called HEAT has been written to calculate the
heat transfer rates Q1 and Q2 from the temperatures, holdups and
flow rates. It is included in Appendix B. The rates 01 and 02 are
found by using all the heat transfer equations in Appendix A. The
inner and outer wall temperatures that figure in the film coefficient
calculations are unknown. These temperatures are initially guessed
and the film coefficients are calculated. The wall temperatures are
adjusted until the equations predict the same heat transfer rate
per unit area across the inside and outside films and the wall. With
the final wall temperatures, the film coefficients and heat transfer
rates are estimated using the appropriate equations depending upon
the nature of boiling and the mechanism.
111.1.4 Decision Variables
Equations (3.1) to (3.6) are the differential equations and
(3.7) to (3.12) are the algebraic equations describing the dynamic
response of the evaporator. After enumerating the number of variables
and the number of equations it is found that there are eight more
variables than there are equations. These eight decision variables
are chosen in a natural way making them manipulative or control
variables. These are:
1) feed flow rate, WF
2) feed temperature, TF
3) feed concentration, CF
4) intereffect flow rate, Wi2
5) recirculation in first effect, W11
6) product flow rate, W01
7) steam temperature in first effect, T
8) total pressure in second effect, P2
111.1.5 Assumptions
A summary of the assumptions made in writing the model is pre
sented here.
1) Time responses of the vapor phase and the tube walls are negligible
compared to that of the liquid phase. This results in simpler algebraic
equations for the vapor phase and tube walls and also decreases the
dimensionality of the model.
2) The vapor is saturated and is in equilibrium with the liquid at the
same temperature.
3) Condensate on the vapor side of the first and second effects is not
subcooled and condensate holdup is negligible.
4) Boiling point elevations due to the presence of solute in the two
effects is negligible. This is justified in the case of dilute solutions.
5) There is perfect mixing in the two effects resulting in lumped
parameter concentration and heat transfer equations.
6) The heat transfer mechanism in the first effect is single phase con
vection followed by twophase convection and nucleate boiling and that
in the second effect is natural convection.
7) Slug flow is the predominant flow pattern in the first effect when
boiling takes place.
8) Heat losses are negligible.
9) The inaccuracy of the heat transfer rates is due mainly to the
uncertainty of the inside film coefficient leading to undetermined
parameters to correct for the inside coefficients alone.
111.2 Parameter Estimation
111.2.1 Stochastic versus Deterministic Estimation
An extensive review of parameter estimation techniques in
differential equations is available in Nieman et al., (1971). In
the deterministic case the simplest and most effective method is a
least squares fit. The problem is stated as follows:
Given the state equations
X = f[X(t), U(t), B(t)]
where U(t) are the control or manipulative variables and o(t) are the
parameters to be estimated from experimental data Y(t). The obser
vations are related to the states and controls by
Y(t) = h[X(t), U(t)]
The problem is to determine the parameters o(t) such that the model
"best fits" the given experimental data Y(t).
Assuming that the parameters are constant, as in the case of
the evaporator, e(t) = 0, the problem can be reduced to a leastsquares
estimation
S2
Min i (Yi Yi)2
0 1=1
where Yi is a calculated observation at time t. (i=1 ,...,N) and Y. is
the actual observation.
The alternative to the above deterministic estimation is the
problem of stochastic estimation. It seemed that the model would fit
the data far better if the parameters 0 were updated as each measure
ment was made. Further, if the states and measurements were subject
to process and measurement noise it would be necessary to estimate the
states using a nonlinear or a linearized Kalman filter. It was apparent
all along that this would require an appreciable amount of computer
time; however, the estimation technique was investigated offline from
an academic viewpoint.
The approach follows that of Padmanabhan (1970). The model
is assumed to have process and measurement noise v(t) and w(t)
X(t) = f(X(t), t) + V(t)
Y(t) = h(X(t), t) + W(t)
where V(t) and W(t) are white Gaussian noise sequences with zero mean
and covariances Q(t) and R(t) respectively. Note that the control
vector U(t) is expressed in terms of the state vector X(t) in the state
and observation equations above. To the state equations could be
augmented the parameter equations
e(t) = 0 + V'(t)
thus treating the parameters 0 as states.
The problem reduces to estimating the states X(t) and the
parameters 0(t) from the experimental data Y(t). A recursive estimate
X(tk/tk) representing an estimate of X at time tk based on all the
data collected until time tk is given by
d t) f(X(tk/tk), t) P(tk)(tk/tk)
with X(O/O) = p = initial estimate of X(t )
where Z(tk/tk) = hR(tk)1 ((tk) h(X(tk/tk), tk))
and the covariance of the estimate P(tk) satisfies a matrix Ricatti
equation
P (tk) = P(tk) + P(tk)fx + P(tk)ZxP(tk) + Q(tk)
with P(0) = n = initial covariance of estimate X(t0). In the linear
case the covariance equation can be integrated offline and stored for
use by the state estimation equations.
To study the effectiveness of the scheme an example was run
offline in which it was assumed that the first and second effect hold
ups and the first effect temperature were constant. It was desired to
estimate the second effect temperature and the parameter 6' in the
equation
= hdh2) Q2
S = H DF(hF h2 +e2
while measurements were obtained on the second effect temperature
T2 = f(h2)
The state equation and the Ricatti equations were integral
for eight minutes of real time and this took 100 seconds of IBM 3i
CPU time. The results for the first four minutes were as follows
Time
30.7
31.7
32.7
33.7
34.7
X4,measured
118.36
122.63
124.32
125.82
126.62
ted
70/165
x4 x4,predicted e'
1.0
121.4 118.73 0.08
124.1 122.37 0.18
126.76 125.05 0.19
127.34 126.35 0.17
This example showed that it was not practical to try online
stochastic estimation for a problem of this nature especially when
there are more state variables. Because of this the deterministic
leastsquares estimate was resorted to and the results obtained were
acceptably good.
To account for noisy flow variables a linear filter was used
on the flow measurements in the form
ui(j) = ci(j 1) + (1 a)ui(j)
where O
ui(j) was the filtered value of the flow ui at t.. ui(j) was the
measurement of the flow u. at t.. a=0.3 was a value commonly used.
It was assumed that the temperatures were measured with a high degree
of accuracy.
111.2.2 Experimental Work for Determining 8la
The experimental runs conducted for determining 81 and 02 were
made in two sets, A and B. The runs in set A were made when the first
effect was being heated and hla was estimated.
Effects 1 and 2 were initially filled to their steady state
holdups. The product flow, intereffect flow and feed to the second
effect were cut off. The recirculation flow, W ,ll was held between
10 and 15 gpm (0.6 to 0.9 kg/s). Heating steam was then started to
the first effect. Recordings of the intereffect temperature T12,
steam temperature, Ts, and recirculation flow, W,11, were made along
with the other variables but it was only these variables that figured
in the ensuing calculations. Tables 3.1 to 3.5 contain the data from
Runs Al to A5.
0 000 0 ~ 4' 0 Q0 o 0
0ON
Loo00o Ao LA LA LA0 oLoLo LoAo4 t4o4ooo
I r P. r ?I C' F P r 1 C P t" P r. el P t P C
LD
W NO'O'4LALAC~.L NOC
Lu C! 4 4o Co  0' N 4 '0 C Co 0' 0t o N L 0
S 00 a I n0 0 I0 00 0 000000 4m
0 WM444 4,000400 0000
S   iN N f N N N rN N
: O0000000000000000000
Oli f
J 0
CO
4C o4 IOII It aItOLaO 4OO''LAIOtm
LU
I LA r00Nr4 IN IN cN M(N cl ooiN
. . . . . .. NNNN NNN
< c 000000000)00000000000
S0. )
C 0 F 0.
on: 3000000ocoooo000oocoo
u LJ LU
uUJ C0o00000LU0)00oo0mnr CT)0Q
S. . ............ .....
UQ Z c0
W 0) L r 0 a, 0
In mA n 4 A A Co C .m  N 4 L 0 0 0 N 0 '0
m mm L G mo 4 4It4 4 4 44o4 4 IA InO LA L0
~ o r W t ,>t t ,.iiiiii
N N N co rsj> ^ .>. ri p I r N r rm 
L 10 1; 0' 0' 0' 0; 0' 0' 0 C 0' 00' 00000' a 0
LL .  
LU 4N 4404e00mmemo0
N0 N c 0 0cfn U N W .N t 4 Or fC. MC oi00
w U% 000000 w0 ma00000000000
LU
a. 'r .4 .4 It It It. n .
SNN NNN c N "N N N N c NN
co coc r' o ooo ooag
0
0 0 V) oD/ o 3
0
J a S Ln*
T 00 0000o000000000000
2 < .  . . . . . . 
m _C', n0I0tC0000000000000
Ij Q n ri
3 I IU
I
4Nc0000NC0N'4000000000
i o 6 u . .... **
>< 'L''33S} sjm \Jt ^ L C s CT U ^fN^ C OO
h^ ^ ^ f^co ^O ^ ~j 'ir ^coF o ^ru ^
^^ ^ ^' ^ l ()f nj+ T fi(jr^ rlfi n f
r r r I r r Ir r r r r r r 10 z 0 1Fi 
*u oO 'Tocccco d  s~t
: 1 0 0 10 P r  r F c m m m w a, m 0 l '0.0
. . . . . .. . .. .. . .
0)0)0r0)0'0'00000000@
Lij
0. (j'0)00rinin 0mm) NMNNNNNN
SNNNNN NNNNNNNNNNNNNNN
00C0ON000000C0CO 00000
F . L I j..... .
4 4)
uL
0 L OC000000 00000000000
0
unw CD
N N N \jN "N"NNMmmmm
CO I I
eC I S ~ s s~Pm~~~
l c o ^ ^J4 o 'dlo^^o ^tM co ^ ^'**^ oo OLr
z~ oo^ cmNNNoN ooorMO' iriorN ooNjiro
t 11
w UU
LUU 00000000000000000000~~
U. U u. . . . . . . . .
UL XOOO 0000000000000000
LUIIIc
I
r/ 0~v
0^0 r s ^' '^' ^i^ wrorio ~o f
c< aj/ 4o<^ i o o i ^(oO ^ n
U J4 l 'C "
^uOtoooj^Of~~~^^OT r c ~'
O I / ~ oo r r000 0n m' rio'o o ir' l
^ I ^ ^r^r jOjjrlt sr MCjcjs n m r ^f
(c n mw co a m a Mrm mN "m m O z 0 "N
LL. L in iCs 4 4 I Ln Un It IT I Id LC afr tr m mr U) un un
Scr' m' c cr CC M' Cr a' Cr' a, a,' O' CO a' M OC C' a' cra a)
0
w U)4M4 M 44't ~4 ,4n
S0000 '0' @0' 0' 0'' 0' a 000000000a
z H Co CL Ma t f I r r r r r I H N t t t H r t r
J  4 4 l  
10 N a 4 O 0 d4 N IU rH co 0' 40 4 a' It L '0 t ao 0r' 0
w r 'n'U')rOOOOOOOOOOODl
0 000 0i0 0 0 N Poso 0M o
10
I 00 N0040'' 0000 000 0Na''00
Lio o o
 ED
r) 3,: C\ C 'CQ r ( ri rj N ()j 4 N C4 f\j
HO
j 0 Q0000 000000 000000
u L LL
LU
UW Z 0)00aM L 00 0 0r M0
a'U, o'o0' Itc H '0 '0 1(M( a'C' r m4 z1(m A
o C'0 m 'T 4 'D N 000 N C%4 m n '0 o 1(0
_N N N N: N N N N N m fl a' a' a' a' r' n en 4.
01 3000COOOOOCOOOOOOOOOO
I 3~n
I 0 9N~~00 ~ ~ 9
> <> )O ~ m~~~~
U 01 U~mm ~mAm
OONNvOOOO^OOOOOlman^O^
LL N .4 4 4 .4 4 N 4 4 4 .4 < ~  44 4
N N . N N N N N N N N N . .4 .
W m00N0'em w omNNN
D N 4 r m 0 t 0 0 m 0 0
SIrO M0100a, c0' 00000 0000 o
<  4 44mm 4 t M j
cNNNNNNNNNNNNNN NNN
I3
o C000 cOOOOOOOOOOOOOCOO
I L/
cO
I 00C
I u m ( n 0oM M 0 r f co Md( lr) I 0 0 0 fn' C n
C C00 N000 C
SNC S004000000000000100
4444 4 4 4 4  4 44 4 4
00
ES it o
ClI C
S(L 30C000 00000C0000000
I 3
U N OL0 0 0 0 0 000 '0 0
0LL LL ***..............
0L IU U
U NOO O O OOOOON 0 N t 0' N cC 0
.4 N N N Ni N^* NI N4* (I C^j N Nq N~ N
il jLUt/iri O ^' ^J fvO c09 t T'c n o ^uc
U.L^ .l) * *
tlUJ M co i ~N fi'Ofj ^ r ^^^ocoo ltftr M 6 O T
t/i ^ r'o 't^t^ t^ ^'^'^o ^ "10^*^ '' '
^ri f ojrN () Nj\JfM(\ rio cj McnrQr^rOf^m
111.2.3 Calculations and Results for 01a
A computer program was written to estimate ,la from the data
collected in Runs Al to A5. The leastsquares program used was sub
routine RMINSQ (Westerberg, 1969). This program was based on a program
coded by M. J. D. Powell and described in Powell (1964). This routine
has the capability of performing a leastsquares search over several
functions in several variables. The search routine does not require
evaluation of derivatives.
The equation describing the enthalpy rise in the first effect
is equation (3.3).
dh
dh1 [12(h2 hl) + h +) + 1(3.3)
dt [12(h12 + V21 (h1 Q1]
Runs Al to A5 were conducted when the first effect was not
boiling and with constant holdups in the first and second effects.
Thus, the liquid entering the first effect was only the recirculated
liquid, W12 = 11 and the vaporization was zero, V21 = 0. The obser
vation T12 = T1 since all the liquid entering the first effect was
recirculated. For every value of 61a which subroutine RMINSQ searched
over, equation (3.3) was integrated from the initial to the final
time for each run. Ten functions of the form (T,calc T,observed)2
over the time span were minimized by RMINSQ. The functions toward the
end of the time interval were weighted one hundred times more than those
at the start. This was because the final time at which the first effect
liquid started to boil was more important from the point of view of
possible switching times in the control variables at this point in time.
The results of the minimization are tabulated in Tables 3.6 to 3.10 and
these values are plotted in Figures 3.2 to 3.6. It can be seen that the
L)
LU
1
_i
O
0 0
W t I
II
I
0
0
0
O
u)
LIJ
1
r
uj
LD
I
(Jo)3aniva3dW3I
Lu
j
c c
II
0 2
o
4
0 o
oo
OO 0
H 
3
0 O
0 a)
\ Uo
\ i"s
0\
0 \
\ ,
CC)
0 0,
00
00 0 00
\ Q)
\ >~n~~d3
1
U)
=3
om
ou o
0 l
0 7
Ln
(do)Hanivy'jdr i
On in
r
m
44
0
u
S00
LU
,0 50
II o
con
\O
\LU 
0a,
0
 O
\ ) LU ,
O CO *
0 a,
0
\ I S
0 u,
0 n O O
CD CD 00
\ 1 '0
0 \S
\U US3
(do)3aniva3dw3I
O ~
\o o C
o 00
0
O
o
0
0 L >
\ / aj
O 21
O C
0o>
0
O1!2
0 =
0 \o
0
C
00
\ <?
3 ^~n~~d3
46
o
wLOONmO ommo OMNt O
c . . . .
w It w a, It ooorco aooomrc
Ujuln rW COCOOO' 'O O C OOOO4
Q; a
0 Z
U
C L*
uij crr
4 < L
o4nVn Jm C) M10sM L
11 Q ui m Cmo oo eCo:om o
o f a< L r u r 00 co u LrI (P C, LL u 0 u Q
LU II J
C)
Sw
HU <
U W O'
) U
0
o f . . . .. . . . . . .
=: uj(A 0 ODo0oo Ch CT'Tcr Qoooooooooo
: 0
C
u u
> i
LU
_1 LLJ II 0 i . . . . . . . . . .
z < <
Cr)
L L
L W
LI)
CLJ
j I 'CCC Ln M L) r o C), M rM~i a %0 M r tA 3
I x rs" f m m c re) M to UNLC 'o
QW m~~OVN~
We ,Ce ee~~NNN N
0
0 COO 1 14 
a 0N:mS 0 0 4 0
C Z
4
t O
I
C3 I N m 00 m 00U
L U
LUj M Lu < co fl) C1 1 O OM co 0 e <1 1 JL CC) rc i N
_J Li It
0 LU '
u
< 4 0 n 0 0,o 0 o0 o Co U
LU
> WuJ
n I
I
0
< Umo m m m e a O 4
0 lcomm N N ONOOOOOooo
.4UJ
LU
a a0
< UJ
L Nrma .,a f OOm 4
c 0 _j) N N N ..s 0..o
V) <
LU 1
LO
LItl
zJ mm rr O o4 0^ c*O 4o It Q
50
0
SU It m m 4 ^ fN 4 mmNUN wO
LU MCCI0 a', a,0'0 0'0 00000000 00o 
of 0
o~ O
uO
LV
LU
w < m^ r 0 0 f C, m, r^ 0 'T ( f Mo w n
co Of x i u, r r ro M wo c^ o~M 00 V, 0' 0' 0 (P 0 0000 W 0 00
i A < <
L m
Lo
I
< z
S
3) m m mmmn1 ar T ItIt '. r n LnmnriLn n
calculated values are indeed leastsquares values (with the weighting
indicated) which converge on the observed values toward the end times.
On the basis of these results a value for ela was taken as
0.1395.
111.2.4 Experimental Work for Determining 61b and e2
The runs in Set B were made with the first effect boiling and
the second effect initially in the heating and later in the boiling stage.
Three parameters were actually estimated from the data for each run.
81b was a parameter for the boiling of the liquid in the first effect,
62a was a similar parameter for the heating of the liquid in the second
effect and 02b was a parameter for the boiling of the liquid in the
second effect. Note that 62b is, in effect, the revised value of 02a
when there is boiling in the second effect.
The experimental procedure consisted of bringing the holdups
in the two effects to their steady state values. The controllers in
the Transducing and Controlling cabinet (Chapter II) were used to main
tain the holdups constant. One controller was used on the first effect
and another on the second. The pneumatic input signal to the first
controller was the output from the DP cell which measured the height
(holdup) in the first effect. The pneumatic output from this controller
was directed to valve CV2 (Figure 2.2). Thus analog control of the first
effect holdup was achieved by manipulating flow W12 which is the liquid
input stream to the first effect from the second effect. Likewise, analog
control of the second effect holdup was achieved by manipulating flow
WF which is the input stream to the second effect.
When the holdups were about constant, steam was let in to the
first effect for heating and a vacuum of around ten inches'mercury
(0.675 bar) was maintained in the second effect. Data were recorded when
the first effect temperature was near boiling. The results for three
runs B1, B2 and B3 are shown in Tables 3.11, 3.12 and 3.13 respectively.
111.2.5 Calculations and Results for 81b and 02
A similar leastsquares search using subroutine RMINSQ was used
to estimate Olb' 82a, and 82b from the data obtained in Runs B1 to B3.
The function evaluations for RMINSQ entailed integration of all four of
the differential equations (3.1) to (3.4) as neither the holdups nor
the temperatures were held constant. Whenever the second effect tem
perature corresponded to the temperature of boiling in the second effect
(which was found from the pressure observed) parameter 82b was used
instead of 02a. The calculated temperature of the first effect solution
T1 and that of the second effect solution T2 were obtained from the
integration of the state equations (3.1) and (3.3) respectively. The
criteria for minimization were the functions
10
f(l) = (Tlcalc Tlobserved) (3.3)
10
f(2) = (T T )2 (3.4)
i= 2,calc 2, observed
Ten functions of each type were evaluated in the time span of each run
resulting in a total of 20 functions for the evaluation of 1lb' 02a and
82b.
The correspondence between the observed and calculated values is
shown in Tables 3.14 to 3.16, while Figures 3.7 to 3.9 are plots of these
values. The minimization took an average of three minutes CPU time on
the IBM 370 for each run. This was mainly due to the threedimensional
53
(N r 0%' 4 0 'r co (0 0o IN 14 Ln LIN Ln m wo Cl) m N P 4 l 4
In C0 cn r .4 '0 LA 0 m 0' 0' i^ f 'n (N (N ( n M a M IN %N
L4 ****
. 4*' 1 It CIO 110 '0 N N O C0 m O m' 0 1 c0 I 1 00 1 a, m 0 1 0
I rN N N N N N N N N N N N N N N N Nt ND M N ND O Ln N NN
a'a' Cy' r  4 74 
< LL 1 1 N N ( 1< 4 ( a'CO C 40'
n r o N N r
W N N N N . N N N . . "N N N N N N N.
S 0oN(NN 00NN 0NNNNNmmmn m 0 A
a Z 0< m 44oo4 444 4444 oooo4o rr
0NoO000 000000000000000000
SL 4 M T N 0000 L000000 00000 000 0  N
m CC rj 0 ( (NC (li C m rN A It? 0.t a, a, 0' 0' a, 0' a, a t 0 0 m0' CO 'I C)
S V)
=3 I l ) T
N u0 0 Ca NO a, a, a 0O Om r 4( Ln( 10 1 4 0 0a,0 c( 0( N0 co Ln l
CO CC 'Nt N0 ()10
a >
u 0 t o nC 1 A L N 'z z a ,O c (3 "1' N 4NIOCN O
Do ee OQW^fwO mmm memm cwmm cor
4
SI > 0C O ;CN0 mN(fLC0(
u C CONNNWACOw m m a ("I w 1rCOCOCOCOCOCOF
2 Li< M N mi r O ms 4 o r 1 1 frl 10 0 It N ^^ I* mro
I *
0 @'Z 0' 0
W O Z ( O000 ? 00 O C O0 0 0 0 0
flu.0 0OD0QCOi LC'OOOOOOOOOOOCOOOc
 UJ 0'( LA N 0 lr '00' 0 t' N0 U4 (LA 0 (N
54
S0000000 0 0 0 0 0 0 '4
LM f M 0 CM 4 ' r 00 o 0 0 4 010 Z 04 04 t4 *
SLn 0 0 W r m r mmr mf mM N
a W OO0, )^ O C CO CT C!N
O N
W. )o4NNN4N' .J.L.NN 4'Q
SNNNNNNNNNNNNNNNNNNNNNNNNN
I m 1: * * * * * o** cr o
II I) fO U r' U ) r'r '4 M C M' a' 0 00 U) U) U) C) a (Ma CO N N
NNORCN N % 0NNNNNNNN:N
Uo ) U) U) U)'7 I 4 44 U) M U) U) U) U) U) U) U) U) U) U) n V) U U )
Z)Z mmUNNNNNNU m mmmmmm)U ) )mU))
Wg a\ m C ol w n 0 zm mOm m N
IuE 0'N0 N00000a0'0 )0'000'0U)0
CL LO 4 '7UU ) ()U % U)U ) UU )U )U ) U CU )
Sm (' 3 0' 03 0 ) 0 0 0 00 7 co r N I N
<
0 N'7' U)0 a, P o0 co U)4 0 s 74 U) ( oa)'0' aN
N3 ftO4n 0 N')4U00U 0nf00000'000n
U z N 30 0 00 C: 0 000 0
I 0
LL c7 ro 0
Sco u co
WZNOOOC4C)N4 4000000OOOOOOOOCO0
9.ULJ 0004u0'CT'70Ur t'OOOOOC)000000
0
L /) U) rU U 'o '0 Co o r^ C C C C^ C r a) a) a) c) C) c co C) 0'
(71O4m w mmmz 0(7%mzmmme
LL . . . . . . . .
cH fl m HH H n H It rn 10 r0 H r H rr r r o a: o OD
O
LA 0OC404A0H40'0' OO O >44
W 0 O ''0rA04'00' L0 0 t444 0L
u 00 00000m0'0'0'0'0' 00' 0*'0 '0'0'
m 0 W Nn .4t N O ' LA N 0 M' M 4 4 T 00 4 0'r 0) 0fl
SNN moco CO N 0OOO C 0 CM' '0 0
NNNNNNNNNNNNNNNNNNNNNNNNN
1 0 0 0D r It N N 0 0' M a, l Ln m 0 4 % u A Lco LA LA o co
m44 4H 000m0 LAAt'
S*0CL u r 00*co *0co0000i0cnci I l M It 0 00 It 00oco0
O 99 C;O 0 4 1 o 9 O OO O Z C; C; C; Z CO Q CO C;DCO CO aO a3 C
Lo LfA In W% LAr V) Ln LA LA L L LnA LA Ln LA LA LA LA L LAr LA LA LA Lr LU
" N N (C N^ C N N N NjJ rJ cI N (j cNj NIj cN Ns N\I N\I r\J
V)
0 **s 4 ( t0 0r H L A t 00 (1 00* Lr% LA0 0 0 r M0
0 :r . . . ..J X..
L <. . . . . . . . . . ..* .* *
S 0 cr' o0fo o o r M NNNN r ' r L00
OZ mJZ C'LtmQfmlNNNN ('LCOLN
M m mfLwre 4H0N'(N 00C''m i
J *. .. *******.. * *
 *R n OLA co m00'HO L0 '1 N COO4 lHO' 0
SJ L)CL m (3m1 To m r r mm01 0, U CC)mm m Nr lo r 0m
j j 00LAN Nm0'NJ m m r000'0'ON N'
H L
1 10 0 r U r 00 0 D 000 M M M L
o 
CU LU LA M In 0 M LA 0 LA 0' 0 N LA 01 0 00' 0 1 LM 0 M
.... ... ... ... ...
0 A40H0 0 4'0 4L A''Hrl0w
U: LA'
O 0 0'O04HO0)r O0'00'0LAH0N NN
ILli HN04m '0'040'4LA'0'LA~t 0H0000'('
w J. rl.. . . . . . . .
UNL ~~9r4) 0 L 0 0 H0'0 1) H0'0 4('
I LA000HH 0 0 0 0' ' 0 0
0
0
0
0
LA_
C
LU
._J
LA
to
DN
0
w~j 1 1
> >
aJ LUI L
V A Ln
CD CO
0 0 r
00
C(M
IC
ro
CL
CD
CD
0
CI
0
F
0)
0 ..
CD c >
rS
3 .1

u
U
C0
a,
U.
LA
U>
O
0 0
r r
0 0
0 Ci
(^o 3ci d i
0
0
0
0
0
0
O
c'J
c' N
(j)n dW3
(jo)]WnVW3JN3
LJ
H_
< co C)
CD
o a
Oo
1t 00
SLn
O o
\ o
II
0
0
Io
S o
[]
[]
[]
O o
0 I I o
oo 3
w3
U0 0
I] o o
c 
0 O
O no
0 0 0
0 a o
0
o 0
0 0 D 0
O 0
0
o 00
O
0
0 0
0
0 cD
O
0
SI 0
0
0a C0d
UlI
LUJ
H
I
59
U. 0
uw
` j
N4 ( mNNNmm r
nZu 0' 0 0' 0'
O* cc
r" <
c .1
Z 0
> It O
ii r ooo
0  
u
0
*A < W
LL
LLJ 1o L o' o iv rv iv N N
ol C
LU y ((U L0t .( rsCMI
i mc < c Ln ff a r% 10 P rO *4
Q C c >
0 
r F 0
Q: Uo>4L
w
LUJ a
O( 0 U 0 ,.4 (l) 0 CC) 0 MA
0 ^
I
I n< ^
60
1. 0
0 Omom NO4 mOONCO
O < 0, 00 N m0 0' N 0 o U s LA Wn 00 aN'
.4 _4 4 .4 .4 . 4 4 4
4 MJ ONUSN~mm~b
CN <' NO
SLU
: 0
o 4 NCJlAmWNfl '4h'OOO040 '0
<4a' >Ofl4 OLrrrmNNN flr NLA0
LA
F m
o 
L0
eo n N N4
C) er
0 LU
as ww aiQo i(Op m C) 0% L%
01 44 ;_:44 4 4 N N
UJ
w V) w
IT) L0 N
i oi > C^CT (looo m
I 0I
0' @0 00 4o Nfl en o 4 LA LA
4444444444444
61
0U
F11
C i 0 0 M M 0' N LA .. 4 .
S4 0 .40 m f Nm N N mN N n
en z Q cT' 0' a, 0, 0a 0' a, a, a, a 0' 0' a, 0' 0Y C c 0 a'
0
r <
C m oIr
F 0 ;
N
o C.j > . . . C4
 w <w n LAO C70' m 1 "0 m 0 .D
Nx > 10 w  U rr ra,m
II0
< 3
0 .
o m
m () m 0 m
1
...... ..... .....
N .0 '0 N N C C C****
4 444444  
03 UJ
o o;~m ~ fu~~~~~
W nWN^?0N~N NN N~
t/^ m 3LJ^\t^^\* *o T ~' ^ t ^^ 4 ~ ^
62
search and the integration of the four differential equations for
every function evaluation. Because of this only three runs were
analyzed. The results for runs B1 and B2 were much better than those
for run B3. On this basis the mean values for elb, 62a and e2b were
taken to be 0.0928, 0.1359 and 0.1763 respectively. Physically,
this meant that the model predicted an inside film coefficient which
was from 80 to 90 percent higher than that obtained experimentally.
This could be either due to the assumptions made in the model or to
the empirical nature of the film coefficient correlations (A.16), (A.39)
and (A.51).
CHAPTER IV
MINIMUM TIME CONTROL POLICY
This chapter deals with the development of a minimum start
up time control policy for the double effect evaporator using the model
equations of Chapter III. The problem is stated in Section IV.1 and
this involves identifying state and control variables, equality
constraints in the form of algebraic equations, state and control
variable inequality constraints, and possible startup scenarios.
Section IV.2 contains the derivation of a general algorithm useful
for solving minimum time problems similar to that for the evaporator.
The actual use of the algorithm for solving the present problem is
described in Section IV.3. It contains the results of model simulations
in arriving at the optimal policy for three problems, all of which are
minimum time startup problems under various conditions. Experimental
verification of one of the minimum time policies and the effectiveness
of the model is presented in Section IV.4.
Refer to the foldout nomenclature list of the more important
symbols at the end of this chapter (page 148) to aid in interpreting
statements made using these symbols in the succeeding sections.
IV.1 Statement of the Problem for the Evaporator
IV.1.1 State and Control Variables
The state variables, which are necessary to describe completely
the state of the process at any particular time, are the differential
variables in the model differential equations (3.1) to (3.6) of
Section III.1.1. For uniformity, x. will be used for a state variable
and X will be the state vector with components xi. These are assigned
as follows:
x1 First effect holdup (H1)
x2 First effect liquid enthalpy (hl)
x3 Second effect holdup (H2)
x4 Second effect liquid enthalpy (h2)
X5 First effect solute concentration (C1)
x6 Second effect solute concentration (C2).
The decision variables listed in Section 111.1.4 have to be
defined so that the model is complete. The control or manipulative
variables are chosen from this set depending upon the controllability
of the process and the physical realizability of the control. For
example, the second effect vacuum pressure is not capable of being
manipulated physically on this system and so it is not chosen to be a
control variable. The feed temperature and concentration are not used
as control variables in this problem either. The remaining decision
variables, comprising four flow rates and the steam pressure to the
first effect, can be easily manipulated physically and can be used to
force the process in any desired direction. For example, the feed to
the second effect, WF, and intereffect flow rate, Wi2, can control
the inventories in the first and second effects. The steam temperature,
Ts, and recirculation rate, W11, have an effect on the first effect
temperature and the rate of increase of the second effect temperature.
The product flow rate, W01, is used to control the product concentration.
These control variables are capable of keeping the process at steady
state and the steady state values for these variables are governed by
the steady state solution of the differential equations (3.1) to
(3.6).
For uniformity, let ui denote a control variable and let U be
the control vector with components ui. The assignment is as follows:
uI Feed to the second effect (1F)
u2 Intereffect flow rate (WI2)
u3 Recirculation flow rate (U11)
u4 Temperature of steam to first effect (TS)
U5 Product flow rate out of first effect (H01).
Rewriting the state equations (3.1) to (3.6) and the algebraic
equations (3.7) to (3.12) in terms of the state and control variable
nomenclature defined above, we have
1 = W12 V21 u3 u5 (4.1)
1
x = W2 (h2 x2) + V21 (x2 h) + Q1] (4.2)
3 = Ul u2 V02 (4.3)
1 2 (4.4)
4 = [l(hF x4) + V2(x4 h) + Q] (4.4)
x5 12(C12 x5) + V21x5]
X1
x6= [u(CF x6) + V02x6] (4.6)
3
Connection equations,
W12 = u2 + u3 (4.7)
+ux(48
W12h12 = u3x2 + u2x4
(4.8)
V12C12 = u3x5 + u2x6 (4.9)
Heat transfer equations,
Q1= Ql(x1' x2 u3' u4, h12' '12' V21) (4.10)
Q2 2(x2, x3' x4, hF) (4.11)
V21 = Q2/A (4.12)
IV.1.2 State and Control Variable Constraints
It is evident that in a real system the control variables
cannot take on all values as there are physical limitations on the
maximum and minimum flow rates and temperatures. The lower limit for
all the flow variables is zero. The lower limit for the steam
temperature is 2120F as steam cannot be supplied at a lower pressure
than atmospheric in the first effect. The upper limit depends upon
the pipe size and the valve size for the flow rates and on the steam
supply pressure for the steam temperature. Thus, all the control
variables are subject to lower and upper bounds of the form
u. < u. < u. (4.13)
Ui,min i i,max (4.13)
In a like manner some of the state variables are constrained.
At steady state all the state variables should be greater than or equal
to their steady state (desired) values (xi)the steady state holdups
are the desired operating holdups, the steady state temperature in
the first effect should be at least the boiling temperature of water at
1 atmosphere, the steady state temperature in the second effect should
be the boiling temperature of water at the pressure in the second
effect, and the steady state concentration in the first effect should
be equal to the desired concentration. The upper bounds are less
clearly defined; for example, the liquid level for the second effect
(holdup) should not exceed the overflow limit. The upper bounds
on the temperatures are dictated by the design specifications and by
characteristics of the solution being concentrated. In general, the
state constraints are given by
x X i,max (4.14)
IV.1.3 Control Scenarios
In the startup of the evaporator, it is useful to visualize
the change with time of the state variables for certain values of the
control variables. It is evident that the control variables determine
the order in which the state variables reach their steady state or
desired values. Intuitively, the optimal policy will endeavor to
force each state variable directly to its final steady state value
and then maintain it while bringing the others to steady state. The
order in which the state variables reach their steady state values
completes the scenario. Theoretically, there should be a total of n
factorial scenarios for a system with n state variables. But most of
these are not possible as certain state variables can reach steady
state only if certain others have. For example, in the case of the
evaporator, the first effect liquid has to boil before the second
effect liquid does. The equations describing a stage in a scenario
are different from those describing another stage. For each stage they
are simplifications of the general equations. A typical scenario
for the startup of the evaporator (and the resulting simplifications
in the general equations (4.1) to (4.12)) is shown below.
Stage A: tO s time t < t1. Filling and heating first effect.
Control variables:
Feed to second effect at maximum flow rate, ul = ul,max
All feed delivered directly to first effect, u2 = u1
Temperature of steam for first effect at maximum value, u4 = u4,max
No recirculation or product flow possible, u3 = u5 = 0.
Resulting state equations:
First effect holdup increasing; x1 = u2, x1(t0) = 0
1
Heating of first effect; x2 = ) + h
No increase in second effect holdup; x3 = 0, x3 = x3(to) = 0
No heating of second effect; x4 = 0, x4 = x4(t0) = hF
No concentration occurring in either effect; x5 = x6 = 0;
x5 = x6 = 5(t0) = x6(t0) = CF'
Time tl, signifying the end of Stage A, is determined when the first
effect is filled, i.e., when xl(tl) = xl"
Stage B: t,1 time t < t2. First effect solution is being heated
and second effect is being filled.
Control variables:
Feed flow to second effect at maximum, ul = ul,max
Flow to first effect stopped, u2 = 0
Recirculation flow set at maximum, u3 = u3,max
Temperature of steam to first effect set at maximum value,
u4 = "4,max
No product withdrawn, u5 = 0.
Resulting state equations:
No change in first effect holdup; xl = 0, x1 = x1(tl) = xl
Heating of first effect; x2 = Q0/xl
Second effect holdup increasing; x3 = ul, x3(tl) = 0
No heating of second effect; x4 = 0, x4 = x4(t1) = hF
No concentration occurring in either effect; x5 = x6 = 0,
x5 = x6 = x5(tl) = x6(tl) = CF'
Time t2, signifying the end of Stage B, is determined when the second
effect is filled, i.e., when x3(t2) = x3.
Stage C: t2 < time t < t3. First effect solution is heated to boiling.
Control variables:
Feed flow stopped, u1 = 0
Flow to first effect stopped, u2 = 0
Recirculation flow set at maximum, u3 = u3,max
Temperature of steam to first effect set at maximum value,
u4 = U4,max
No product withdrawn, u5 = 0.
Resulting state equations:
No change in first effect holdup; x1 = 0, x1 = xl(t2) = x
First effect is being heated; x2 = Ql/xl
No change in second effect holdup; x3 = 0, x3 = x3(t2) = x3
No heating of second effect; x4 = 0, x4 = x4(t2) = hF
No concentration occurring in either effect; x5 = x6 = 0;
x5 = x6 = 5(t2) = x6t2) = CF'
Time t3, signifying the end of Stage C, is determined when the first
effect starts to boil, i.e., when x2(t3) = x2'
Stage D: t3 < time t < t4. First effect solution is boiling (and
becoming concentrated). Second effect solution is being heated.
Control variables:
Feed to first effect is equal to vaporization from first
effect to maintain constant holdup, u2 = V21
Feed to second effect set to maintain constant holdup in
second effect, ul = u2
Recirculation at maximum, u3 = U3,max
Temperature of steam to first effect at maximum, u4 = u4,max
No product withdrawn, u5 = 0.
Resulting state equations:
No change in first and second effect holdups; x = x3 = 0;
x = x (t3) = 1 X; x2 = x2(t3) = 2
First effect is being heated; 2 = [u2x4 + V21h ,
Xl
x2(t3) = x2
Second effect is being heated; x4 = [ul(hF x) + ];
x4t3) = hF 3
First effect solution is being concentrated; x (u2CF);
x5(t3) = CF
No concentration of second effect solution; x6 = 0; x6 = x6(l
Time t4 signifying the end of Stage D is determined when the second
effect starts to boil, i.e., when x4(t4) = x4'
Stage E: t4 < time t < tf. Solution in both effects is boiling and
being concentrated.
Control variables:
Feed to first effect set equal to vaporization in first
effect to maintain constant holdup, u2 = V21
3) = CF
Feed to second effect set to maintain constant holdup in
second effect, ul = u2 + V02
Recirculation at maximum, u3 = u3,max
Temperature of steam to first effect at maximum, u4 = u4,max
No product withdrawn, u5 = 0.
Resulting state equations:
No change in first or second effect holdups; xl = x3 = 0;
Xl = xl(t4) = X', x3 = x3(t4) = x3
First effect is being heated; = (u2x + Vp1h)
xl.
No change in second effect enthalpy; x4 = 0, x4 = x4(t4) = x4
First effect solution is being concentrated; x5 = (lx6)
x1
Second effect solution is being concentrated;
6 = (ulCF u2x6), x(t4) = CF.
x3
Time tf, signifying the end of Stage E, and consequently the final
time, is determined when the first effect concentration reaches a
desired value, i.e., when x5(tf) = x5'
After the final time, tf found above, the control variables
Ul, U2 and u5 are obtained from the steady state solutions of the
differential equations (4.1), (4.3) and (4.5). The feed to the first
effect is such that a constant holdup is maintained in the first
effect, i.e., u2 = u5 + V21. The feed to the second effect is such that
a constant holdup is maintained in the second effect, i.e., ul = u2 + V02.
The product rate is such that a constant product concentration is
x6
obtained, i.e., u5 = u2 x5 Simplifying the above three relationships
we can obtain a sequence for determining the control variables ul,
u2 and u5 at steady state as follows:
x6 / x6
u2 V21 1 = u2 + 02; u5 = u2 5
Note that certain point constraints of the form xi(tj) = xi
separate the various stages of the scenario. The control variables,
state and algebraic equations change when these points in time are
encountered.
It is assumed that certain state variables such as the first
and second effect holdups will be maintained at their steady state
values once these are attained. That is, the point constraint
xi(t.) = xi is followed by the equality constraint xi = 0. This
latter equality is maintained by calculating the required value of an
appropriate control variable. This is similar to the concept of a
first order state variable inequality constraint (Bryson et al., 1963),
(Bryson and Ho, 1969).
It is clear that another scenario will give rise to a different
ordering of the point constraints, xi(t.) = xi. depending upon the
order in which the state variables reach their steady states. The
optimal scenario is the one which will result in a minimum final time.
This scenario approach simplifies the mathematical problem
to a great extent as a minimum number of state equations have to be
integrated. The equations simplify considerably leading to simpler
adjoint equations and Hamiltonian minimizations.
IV.1.4 Summary of the Problem Statement
The problem is to minimize the final time, that is
Min J = tf
U
subject to the state equations (4.1) to (4.6), the connection equations
(4.7) to (4.9), and the heat transfer equations (A.1) to (A.52) where
the control variables are constrained,
U < U < U
min  max
and point constraints pertaining to the particular scenario are to be
satisfied
xi(tj) = xi ; i = 1,n
j = l,n
IV.2 A Minimum Time Algorithm
IV.2.1 General Problem
The objective is to minimize the final time by selection of
the controls U
Min J = tf (4.15)
U
subject to:
the state equations
X = f(X, U, Z) ; X(O) given (4.16)
the algebraic equations
g(X, U, Z) = 0 (4.17)
the control constraints
Umin < U < Umax (4.18)
or
hli(u.) = ui i,max
1 i i = 1 ,mXi i
(u Uiin h(U) < 0 (4.19)
h(u,) = Ui.min Ui < 0
and the point constraints
Xo(to) = X
Xkl(tkl) = Xkl
xk2(tk2) = k2 (4.20)
Xf(tf) = xf
t = fixed
IV.2.2 Lagrange Formulation and Necessary Conditions
Assume that the point constraints are met at times tO, tkl ...'
tf. Also assume the times ti correspond to all those times when one
(or more) of the controls or states reaches or leaves a constraint.
Next, form the following index sets,
10 = {0, 1, 2, ...., kl, k2, ......... f)
I1 = {1, 2, ..... .. kl, k2, ......... f}
12 = {0, ......., kl, k2, ........, fl}
Within each set, order the indices so the times ti, i e I0 or I1 or 12
are in increasing order. Also form the index set
K = {kl, k2, ........ }
Introducing multipliers a, and A we can write the Lagrangian
for the problem
L(X, U, ti, i E I, a, B, X) tf + af(xf Xf(tf)) + 'k(xk xk(tk))
t keK
+ (AT(f x) + BTh)dt (4.21)
isl1 il
The algebraic equations g(X, U, Z) = 0 are not included.
These will always be satisfied by solving for the dependent variables
Z and substituting the resulting values into the state equations.
Also the equations at tO, X(t0) = XO, will always be satisfied and are
not included. Note that the fourth term accounts for the changes in
the state equations, algebraic equations, and control variable
constraints along the trajectory. This term has been written to allow
for possible discontinuities in the Lagrange multipliers, A, at the
entry and exit corners of constraints.
We define the Hamiltonian H = ATf and rewrite equation (4.21)
L(X, U, ti, i 10c a, a, x) = tf + af(xf Xf(tf))
+ ak(xk Xk(tk))+ i I (H TX + gTh)dt (4.22)
kcK iel
Sil
On an extremal solution, the Lagrangian L must be stationary
with respect to small, arbitrary perturbations of the times ti(i e Il);
the multipliers a, , and X; and the states X and controls U. To
derive the necessary conditions for a stationary point of the Lagrangian,
we take first order variations of L with respect to ti(i E Ii), a, B,
A, X and U,
6L = (1 + H H T + Th fxf)t tf af xf(tf) + (x^ xf(tf))6af
+ k [ ck6xk(tk) + (xk xk(tk))6k (akXk)tk6tk
+ [(H AT + JTh) t6ti (H ATX + OTh)+ 6t l
ieI2 1 il 1
ieI + DX I X T
iI
2T ah Ta 6X
aX
+ BT Uh JU + hT6b dt
3U1 1
Integrating the term
ATsxdt = (AT6x) +
ii
+ HuT
(4.23)
 AT6xdt by parts,
 (T6x)
ti
+
t+
iI
AT6xdt
(4.24)
Substituting (4.24) into (4.23) and collecting coefficients of
6c, 6B, 6x, 6X, 6U, 6xk and 6ti separately, we have
6L = (1 + H T + Th a fxf)tftf
af xf(tf) AT(tf)6X(tf)
+ AT(t0)X(t0)
+ (xf Xf(tf))6af
i l2
ilK
(AT) 6X(ti) (AT) 6X(t.)
it ti
+ I [(H ATx + Th) 6ti (H 
id82 ti
i/K
+ I [(T ) +6X(tk) (AT) 6X(tk) +
kcK tk tk
ATX + Th) +6ti
t.
(f)
ck6xk(tk)] (g)
+ I (xk xk(tk))6k
kEK
+ I [(H XT + BTh) _6tk (H XT + B Th) +tk
kK tk tk
+ (kk)t tk]
+ i H 6
ti
+ 
ti
ti
+
i1
[iT +2 T + T ah ]6xdt
T T 2XT
ax ax
[ H + BT 3h ]Udt
aU aU
t
+ 1 hTSadt (s)
1+
i1 (4.25)
The KuhnTucker conditions arise from the KuhnTucker multipliers B and
the inequality constraints h,
B.h. = 0 (t)
> 0 (u)
For stationarity of the Lagrangian all the coefficients of 6X,
6U, 6xk, 6a, 62, 6X and 6ti occurring in equation (4.25) must be zero.
By equating these coefficients to zero we arrive at the following
necessary conditions (NC).
From term (q);
From term (r);
From term (p);
From term (t);
From term (u);
From term (h);
From term (d);
From term (b);
From term (c);
From term (e);
From term (g);
From
From
From
T aH T ah +
A T t < t < t
1 T 1 i 1
aH T h +
  = 0; ti1 < t
SaU
X = f(X, U); t < t
.jh = 0; t < t( t
Sj i 1
> 0; t < t
i1l t
^k xk(tk) = 0; k K
xf Xf(tf) = 0
f(tf) = af(tf)
Xj(tf) = 0; j=1,...,n; j/f
xj(tO) = unknown since 6x(to) = 0
Xj(tt) = xj(ti); j=1,.. ,n
icI2, ijK
x.(t+) = Xj(tk); j=1,...,n; j/k; kcK
xk(tk) = k(tk) ak; ksK
(NC4), (NC5) and (NC8); H(tf) = 1
(NC4) and (NC11); H(t) = H(tI); iCI2, i/K
(NC4) and (NC12); H(tk) = H(tk); keK
(NC1)
(NC2)
(NC3)
(NC4)
(NC5)
(NC6)
(NC7)
(NC8)
(NC9)
(NC10)
(NC11)
(NC12)
(NC13)
(NC14)
IV.2.3 Comments on the Necessary Conditions
The state equations and point constraints are the necessary
conditions NC3, NC6 and NC7. The KuhnTucker conditions NC4 and
NC5 indicate that, when a control constraint h is encountered, that
constraint is held by choosing U such that h(X, U) = 0 provided that
the KuhnTucker multiplier B is nonnegative.
The problem is set up such that the point constraints arise
when there is a change in the form of the state equations and,
consequently, a change in the adjoint equations and Hamiltonian.
For example, before time t1 (where x,(t) = X) let f(1) represent
the state equations, or
X = f(1)(X, U) t < t1
After time tI let f(2) represent the state equations, that is
= f(2)(X, U) t > t1
The Hamiltonians are
H() = Tf(1) and H(2) = XTf(2)
The functional forms f(1) and f(2) reflect the change in the state
equations. Condition NC14 shows that the Hamiltonian is continuous
across the times tk where point constraints are encountered. For
example, at t = tl, xl(tl) = xl and I() = H(2)
Condition NC11 states that the Lagrange multiplier (or adjoint
variable), corresponding to the particular state on which there is a
point constraint at ti, has a discontinuity at ti, whereas the
multipliers corresponding to the other states are continuous at ti.
At time tI, for example,
x (t+) = (t) 
Xj(t+) = Xj(t) ; j=2,...,n
The value of the "jump," a1, in the multiplier Xl,is readily determined
from the condition H() = 1(2)
At times ti, i E I, when control constraints have to be
satisfied, all the multipliers are continuous as the trajectories
enter and leave the constraintsIC10. The Hamiltonian is also
continuous at such points in time as shown by condition NC13.
Necessary conditions NC8 and NC12 provide values for the adjoint
variables and Hamiltonian at the final time. The adjoint variables
corresponding to the states that are unconstrained at the final time
have a value zero. If only one state is constrained at the final
time, the adjoint variable corresponding to this state can be
determined from the final condition on the Hamiltonian
H(tf) = 1
or ff
However, if more than one state variable is constrained at the final
time, the values for all but one adjoint variable have to be guessed
and the remaining one adjoint variable can then be determined from the
final condition on the Hamiltonian. Let us suppose that the adjoint
variable corresponding to the stopping condition Xf(tf) = xf is
determined from the final value of the Hamiltonian and that another
variable A. has been guessed. On subsequent iterations this value of
A. is updated by noting that the gradient of L with respect to Ai(tf)
is (xi xi(tf)). The Saddle Point Theorem requires us to maximize
L with respect to Xi(tf) which is equivalent to driving the gradient to
zero. So the value of Ai(tf) on following iterations should be in such
a direction that at the final time the deviation of xi(tf) from the
desired value xi, i.e., (xi xi(tf)) is driven to zero.
This same iterative technique is used if more than one point
constraint is encountered at an intermediate point in time tk. One
of the adjoint variables at time tk can be determined by using the
continuity of the Hamiltonian at time tk. The other variables will
have to be guessed initially and then updated using the gradients
available on subsequent iterations.
Necessary condition NC2, along with the KuhnTucker conditions
NC4 and NC5 enables us to replace condition NC2 by another condition
made possible by invoking the strong minimum principle. Denn (1969)
has lucidly shown how the Hamiltonian takes on its minimum value for
the optimal decision function U(t), both on and off the inequality
constraints. Utilizing this result, we can replace the stationary
condition NC2 by
Min H(X, U, X)
U
subject to Umn U Umax
min  max
An algorithm for arriving at the minimum time policy was
IV.2.4 Minimum Time Algorithm
developed based on the necessary conditions arrived at in Section
IV.2.2. It can be classified as a "Minimum H" algorithm as it deals
with direct Hamiltonian minimizations along the trajectory. Proposing
a control scenario to start with is essential to simplify the problem
and to use the algorithm effectively.
Let the various elements of the state vector be divided into
3 groups as follows:
Group A: States which remain unconstrained through startup.
Group B: States which meet their steady state values during
startup and which define the point constraints.
Group C: States which meet their steady state values at the
final time.
The state equations are modified along the trajectory when
the point constraint conditions are met by the variables in Group B.
Also, all algebraic equations are satisfied throughout the trajectory,
and this is implied in the state equations. When a variable belonging
to Group B arrives at its steady state value, it is assumed that for
subsequent time the state equation is replaced by the corresponding
algebraic equation x = 0. Thus the Hamiltonian H = T f is different
for points on the trajectory depending upon which of the state
equations are active; so also are the adjoint equations
3H T ah +
^1 B ; t t < t.
aXT 3XT il t I 1
The algorithm proceeds as follows:
1) Guess a nominal control policy U which will cause a stopping
condition, xf = xf, to be satisfied at the final time.
2) Integrate the state equations forward until the stopping
condition is satisfied. This implicitly determines the final time.
The proper set of state equations should be integrated depending upon
the point constraints, the control constraints, and the algebraic
equations which have to be satisfied along the trajectory.
3) At the final time determined in step (2), let all the
multipliers of the variables in Group A be zero. Guess the multipliers
of all the variables in Group C with the exception of that corresponding
to the stopping condition. Determine this latter multiplier from the
final condition on the Hamiltonian,
fxf = 1 ixi
icC
iff
4) With the values of the adjoint variables at the final time
determined in step (3), integrate the adjoint equations in the reverse
direction. At times tk, when point constraints were met on the forward
integration, determine the values X(tk) by utilizing the continuity
and jump conditions on the adjoint variables and Hamiltonian. For
example, if only one point constraint of the form xk(tk) = k is met
at time tk, then Ak(tk) can be determined from the continuity of the
Hamiltonian
Xk(tk)fk(tk) + Xj (tk)fj(tk) = kxj(t )fj k)
jEn jen
jik jUk
and Xj(tk) = Aj(tk); jen, jUk. If more than one point constraint is
met at time tk, then the values of all but one of these multipliers
should be guessed at time tk and the last one determined from the
continuity of the Hamiltonian.
5) Simultaneously on the reverse iteration, minimize the
T *
Hamiltonian H = Tf at each point to determine the optimal U
Min H(X, U, X)
U
6) On reaching time tO, update the control policy used on
the forward integration with that found on the reverse integration in
step (5)
Ui+ = Ui + C(t)(U* Ui)
c(t) is chosen to limit the change in U if too large a change is
indicated.
7) Integrate the state equations forward as in step (2).
At the final time determine the difference in the states for the
variables in Group C from their desired values and update the guess on
the multipliers .i such that the gradient of L with respect to
Xi(tf) is driven to zero, i.e., to drive (xi xi(tf)) to zero. As
before, determine the X corresponding to the stopping condition from
the final value of the Hamiltonian.
8) Integrate the adjoint equations in the reverse direction
as in step (4). Update the guess on Xj(tk); jen, j/k in a similar
manner as in step (7). Determine Xk as before from the continuity
of the Hamiltonian.
9) Repeat steps (5)(8) until
(1) 6J < eI no significant improvement in the final time
(2) IUi+l Ui 1< 2 no significant change in the control policy
and (3) 116Al I < n at tf
and (4) 116XXII < n at tk.
.k'
The optimal policy is chosen as the one which satisfies the
above conditions.
IV.3 Solution to the Evaporator Problem
IV.3.1 Problem 1. Constraint on the Second Effect Holdup
The problem was solved for the scenario described in
Section IV.1.3. The concentration equations (4.5) and (4.6) were not
used in this simulation. The Hamiltonians and adjoint equations for the
various stages of the scenario are given below.
The general Hamiltonian for the problem is:
H = (W12 V21 u3) + [W12(h12 x2) + V21(x2 h ) + 1
+ 13(u1 u2 V02) + 4 [u(hF x4) + V2(x4 h) + 2
This simplifies for the various stages as follows:
Stage A: t < t < t
H = Au +2
Ha = 2 + [ 2(4 x2) + Q1] + 3(ul "2)
Stage B: t1 < t < t2. x (tl) = x
A2
Hb = X + 3(u1 u2)
x1
W12 = u3 ; u2 = 0 or fI = 0
Stage C: t2 < t < t3 x3(t2) = x3
x1
ul = u2 or f3 = 0 and f = 0