Ind. Eng. Chem. Process Des. Dev. 1982, 21, 297-301
297
On-Line Process Identification Randall 0. Brantley,' Robert A. Schaefer,2and Pradeep B. Deshpande" Department of Chemical and Environmental Engineering, University of Louisville, Louisviile, Kentucky 40292
A computer-oriented procedure for process identification of overdamped second-order plus dead-time systems is presented. The procedure is readily adaptable in on-line process identification work and has excellent potential in adaptive tuning of digital control loops. The method is based on an experimental test on the open-loop process. The experimental output is compared to the predicted output which results from the computer solution of an assumed process model. The model is subjected to the same disturbance as the experimental process. Any difference between the experimental and predicted outputs represents the error in the assumed process parameters. The
parameters are adjusted using a unidimensional nonlinear regression procedure until the error is reduced to an acceptable level. This identification procedure has been applied to some systems reported in the literature. The resutts show the excellent capability of the technique in arriving at the optimum values of the process parameters.
Introduction Dynamic information plays an important role in the design of process control systems. This information is generally available in one of three forms, viz. process reaction curves resulting from open-loop step tests, frequency response curves based on pulse testing, or as transfer functions derived from mathematical models. The frequency response data are used in finding the parameters of conventional controllers and in the design of advanced control algorithms. A survey of literature shows that there are many techniques for finding the dynamic parameters of processes. Ziegler and Nichols (19421, Miller (1967), and Sundaresan and Krishnaswamy (1978) have presented methods for finding first-order plus dead-time models from process reaction curves, while Oldenbourg and Sartorius (1948), Smith (1959),Meyer, et al. (1967),Moore et al. (19691,Sten (1970), and Sundaresan et al. (19781, have dealt with second-order with dead-time models. The theory and methods of application of the pulse testing technique for dynamic identification have been presented by Dreifke and Hougen (1963),Hougen (1964),Murrill et al. (1969),Patke and Deshpande (1977),Shork and Deshpande (1978),and others. The above literature deals with open-loop stable processes. Deshpande (1980) recently proposed a technique for process identification of open-loop unstable processes. The importance of on-line identification and tuning methods has increased with the increasing number of computer control applications in industry. It is well known that dynamic parameters change with changing operating conditions, equipment fouling, etc. One objective of on-line techniques is to identify the parameters of the process from time to time so that controller parameters may be adjusted, if necessary. In this paper we present an efficient on-line procedure for dynamic identification of overdamped second-order processes with dead-time. The Method (Schaefer, 1978; Brantley, 1981) The dynamic behavior of many chemical engineering processes may be approximated by second-order process models with time delays. The transfer function of this class of processes is of the form
The objective of the process identification procedure is to determine the parameters of the model K , Od, T ~ and , 72. The procedure consists of the following steps. (1)A suitable regulated disturbance (e.g. a pulse) is applied to the input (i.e., the manipulated variable) of the process while it is operating at steady state in manual control. (2) The input and the output data are recorded (the output is the controlled variable). (3) Trial values of the parameters K , dd, T ~ and , T~ are assumed. (4) Using these trial values and the input data of step 1, the output data are predicted using the time-domain equivalent to eq 1. (5) The difference between the experimental output and the predicted output is an indication of error in the trial parameters. (6) The parameters are updated and steps 4 and 5 are repeated until the error is minimized to an acceptable level. The detailed procedure is presented in the following paragraphs. The time-domain differential equation corresponding to eq 1 can be solved in one of two ways (Smith, 1972). One way is to begin with the differential equation corresponding to eq 1 which is
and express the derivatives in the form of finite differences, giving an algebraic equation relating the output Xi to the input Vi at the ith time instant. The other way is to use a zero-order hold and process combination to describe the.relationship between the input U and the output X. That is (3)
Equation 3 upon inversion gives an equation relating X i to Vi. The equation relating Xi to Vicontains the model parameters which are to be determined. The latter approach has been used in this investigation. The expressions for Gb(s) and G,(s) are substituted in eq 3 to give
'E. I. du Pont de Nemours and Company, Louisville, KY 40201. Exxon Chemical Americas, Baytown, TX 77520.
where T = sampling period.
0196-4305/82/1121-0297$01.25/0
0 1982 American Chemical Society
298
Ind. Eng. Chem. Process Des.
Dev., Vol. 21, No. 2,
1982
Using the translation theorem, eq 4 can be written as
Now the dead-time 8d is expressed as 8d
=8
+ NT
(6)
where N is the integer number of sampling periods in 8d (thus if 8d = 0.76 and T = 0.5, then N = 1 and 19 = 0.26). Equation 5 can now be written as
where 2, represents the process of taking the modified 2-transform and m is the modified 2-transform operator which is given by m = 1-(B/T) (8) The modified 2-transform of the expression in eq 7 may be found in Tables (Deshpande and Ash, 1981) which upon substitution and some simplification gives
X ( Z ) / W )= [blZ+“+l)+ b22-(N+2) b32-(N+3)] / [ 1 - a12-l
+
+ a22-2] (9)
where
al = e-T/rl + e-T/rz a2 = e-T/rl x
e-T/r:,
b, = K ( 1 + [72e-mT/r2- 7 1e-mT/rl 1/[71 - ~ ~ 1 ) + e-T/rz)- ~ ~ e - ~ ~+ / ~ 2 ( 1 b2 = Kf-al e - T / r l ) l / [ ~-l r211 (10)
+
and b3 = K(a2- [71e-mT/r1 X e-T/r2- 72e-mT/7z x e-T/r1]/[71- ~
~ 1
Cross-multiplying the terms in eq 9 and inverting gives a difference equation relating X’s to Us. Thus xi = a1Xi-1 - azXi-2 + blUi-N-1 + bzUi-N-2 + b ~ U i - ~ - 3 (11) where the subscript i represents the current sampling instant, i - 1 represents the last sampling instant, etc. Equation 11 shows that the problem of modeling the process has been reduced to finding the values of K , ed, T ~ and , r2 that best fit the data. With a given (i.e., assumed) set of these parameters, the values of the constants al, u2, bl, b2,and b3 are calculated using eq 10. Then each output data point Xiis predicted using eq 11. The difference between the experimental output and the predicted output provides a measure of fit which may be expressed as ERROR = :(Xi, - X i J 2 i=O
(12)
The lower the value of ERROR, the more accurately the parameters fit the data. Thus to find the optimum set, , 7 2 must be adjusted the regression variables K , @d,T ~and until ERROR is minimized. This involves numerous iterations. Finding optimum values of the regression variables can be accomplished by adjusting all the variables simultaneously or by adjusting one variable at a time. The first method usually involves first and second partial derivatives of ERROR with respect to the regression variables. The calculation of these derivatives involves very bulky equations which are prone to significant round-off error. Consequently the single-variable approach is used in this work. The method is easy to use and has worked well provided K , Od, T ~ and , 7 2 are positive real numbers. Two different search techniques are applied to each of the four variables. The first of these was developed by Davies, Swann, and Campey (Himmelblau, 1972) which will be referred to as the DSC search. This method is excellent in that is quickly finds the approximate location of the optimum value of a variable. In other words, the DSC search “brackets” the optimum value. The second search developed by Powell (Himmelblau, 1972) rapidly zeroes in on the optimum once the variable has been bracketed. Since the cyclic regression operation is unidimensional in nature, only one variable is operated on at a time while the remaining three variables are held constant. First K is optimized. Next, 8d is optimized using the new value of K and the old values of 7 1 and 72. Then r1 is optimized using the new values of K and 8d and the old value of 7 2 . Similarly, 7 2 is optimized. At this point K is reoptimized since there are new values of 8d, 7lSand 72. Then ed, r1 and 7 2 are reoptimized. The cyclic process is repeated until a preselected “goodnessof fit” criterion is satisfied as will be explained. A step by step procedure for optimization of K is presented here. Exactly the same operations are required for optimizing 8d, 71, and r2except that a different variable is searched while the other three variables are held constant. Let F ( K ) = ERROR evaluated at K with Od, r l , and r2 held constant. Let DK be an increment size. (1)calculate F(K) with the initial values of K , od, rl, and 72.
(2) Calculate F(K + DK). If F(K + D K ) 5 F ( K ) ,go to step 3. If F(K DK) > F ( K ) ,let DK = -DK and recalculate F(K + DK). (3) Let the index j = 1. ) (4) Let K.+l = K . + DK. (5) Calcuiate p ( k j + ~ ) . ( 6 ) If F(Kj+,)IF(Kj)double DK, increase index j by 1, and return to step 4. If F(Kj+,)> F(Kj),denote Kj+l as K,, Ki as Km+ and Kj-, as Km-2. Let DK = -0.5DK and return to step 4 and 5 for one more calculation. (7) Denote the last value of K as K,+l. A graphical representation of the procedure is shown in Figure 1. (8) From the four equally spaced values of K in the set (Km-2,Km-l,K,, K,+,) discard either K , or Km-2,whichever is the farthest from the K corresponding to the smallest value of ERROR in the set. Let the three remaining values of K be denoted as K,, Kb, and K,, where Kb is the center point and K , = Kb - DK and K , = Kb + DK. (9) Carry out a quadratic interpolation to estimate a pseudo-optimum value of K , given by the equation
+
.
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 2, 1982
299
1
READ
+
DATA
K OPTlMlZATlON
!
osc
SEARCH - Part
2
P O W E L L SEARCH
VALUES K, TI ,T> €&
PRINT,,
"CVC LE,
b
Q 'm-2
'm-I
'm+l
K
Km
r
I
osc
1
SEARCH - P a r t 2
Figure 1. Illustration of DSC search.
(10) Calculate F(XS0PT) and denote it and FSOPT. These ten steps comprise the DSC search. The Powell search is now applied. (11)Let XSTAR = XSOPT, FSTAR = FSOPT; and FA = F(K,), FB = F(Kb),FC = F(Kc). (12) If XSTAR and whichever of the set (K,, Kb, K,) corresponding to the smallest ERROR differ by less than the prescribed accuracy in K , XACCUR, terminate the K search. Otherwise, discard from the set (K,, Kb,K,) the one that corresponds to the greatest value of ERROR (unless the bracket on the minimum of ERROR will be lost by doing so, in which case discard instead the K so as to maintain the bracket). (13) Let XSTAR replace that K which was discarded from the set (K,, Kb, K,). Let FSTAR be the corresponding value of ERROR. Thus there will again be three elements in the set. (14) Calculate a new XSTAR according to the equation
+
& POWELL SEARCH PRINT'WTIMUM VALUES:' K,T,,TzOd
1
(15) Calculate F(XSTAR) and denote it a FSTAR. (16) Go to step 12. The last five steps complete the Powell search. Thus the variable K has been optimized. The next step is to optimize 6d using the new K value (i.e., XSTAR) and the original 71 and 72. Then 71 and T~ are optimized as explained earlier and the cycle is repeated as described earlier. At the end of each cycle, the values of K , dd, 71, and 7 2 are subtracted from those of the previous cycle. If the differences between subsequent cycles are within prescribed accuracy, YACCUR, the procedure is terminated. When YACCUR accuracy level is achieved the optimum values of K , 6d, T ~ and , 72 are obtained and the modeling procedure is complete. Program Testing and Results A flow chart of the computer program which implements the DSC and Powell searches is shown in Figure 2. A digital computer program for determining the optimum process parameters has been developed. The inputs to the program are the values of disturbance input vs. time, values of the output vs. time, initial trial parameters, and the sampling period (time between successive samples) which equals 0.5 time units in all the runs.
@d
DX
OPTIMIZATION SEARCH- P a r t 1
I
1 DSC SEARCH - P a r t 2
L_T___
Q POWELL SEARCH
Notation
VACCUR = T e r m i n a t i o n c r i t e r i a , m i n i m u m a l l o w a b l e d i f f e r e n c e o f r e g r e s s i o n v a r i a b l e s o f s u b s e q u e n t CyCleS
IC
= Cycle
number
Figure 2. Computer flowchart for implementation of DSC and Powell searches. Table I. Comparative Assessment of Modeling least-squares test 1, test 2, fit from step anal, pulse anal, Sundaresan Sundaresan this work this work et al. (1978) et al. (1978)
+
XSTAR = 1/2[(Kb2- K:)FA (K: - K2)FB ( K 2 Kb2)Fc]/ [ (Kb - Kc)FA + ( K , - K,)FB + (K, - Kb)FC] (14)
1
K 7,
r2
1.0008 0.3407 1.921 1.249
0.99708 0.3309 1.731 1.441
1.0=
1.oa
0.33 1.67 1.50
0.315 1.69 1.49
By definition,
To check the program capability, several tests were conducted. The first test is based on a known third-order transfer function G,(s) = 1/[(0.5s
+ l)(s + 1)(2s+ 111
(15)
The step response of the process whose transfer function is given by eq 15 is determined by Laplace transform techniques. The step input and the resulting output data are read into the computer program as inputs. The execution of the program resulted in the optimum process parameters which are shown in Table I. Also shown in Table I are the process parameters reported by Sundaresan et al. (1978) for this transfer function and those resulting from a least-squares fit. A comparison of the step response has shown that there is no significant difference between the three methods tested. The second test is also based on the same third-order transfer function but uses a pulse as the disturbance input. The resulting process parameters are also shown in Table I.
300
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 2, 1982
Table 11.
Effect of Accuracy Levels, Starting Values on Model Parameters test 3
starting values 1.0 1.0
K Od
final values in 212 cycles
7 2
test 4
test 5 final values in 11 cycles
starting values 1.0 1.0
1.008 1.814 3.947 3.423
6.0 2.0
71
__.-
starting values
1.031 1.828
6.0
5.460
2.0
2.311
1.0 1.0 2.0 0.5
0
,I -
/ ,Model
Model A Model B
C
Accuracy 4ccuracy
Single
b b s o l v t e (212 c y c l e s ) I*/. i l l c y c l e s '
Time
Model
COn5Ianl
I
fin a1 values in 37 cycles
1.045 3.442
0.5
6.545
= 0.9 w , = 0.5
0.02
0.78
3L
-C'O
I
parameters from Smith (1972)
1c TIME
UNITS, t
Figure 4. Error between fourth-order process and models. selecting a suitable data acquisition system and the reader is referred to his text for details. The type of input used in experimental testing is also equally important. All things considered, a pulse input is probably the most suitable input in experimental dynamic testing based on the proposed identification method. Finally, the sampling period selected must be small enough to permit a good representation of the continuous response of the process by a sampled-data response of the model. Conclusions This paper has presented a useful technique for process identification of overdamped second-order plus dead-time processes. It can be used off line but, more importantly, it can be adapted in on-line digital control applications. In such cases the approach would be to bypass the control algorithm and perturb the process with a suitable input and record the input and output data. The computer program residing in background memory would analyze the data and determine the optimum model parameters. These parameters would then be used to update the controller parameters. This procedure would provide the capability of adjust controller parameters as and when necessary to maintain the controlled variable at set point. With some modifications the program presented in this paper can be used to determine optimum tuning constants of digital controllers. When this is done, it should be possible to automate the entire identification and controller adjustment procedure in digital control systems. Nomenclature a,, a2 = constants in eq 10 b,, b2, b3 = constants in eq 10 C, = predicted value of process variable C = process variable db= zero-order hold G = process transfer function K"= process gain m = variable in eq 8, 1 - (813") N = integral number of sampling periods in Bd
Ind. Eng. Chem. Process Des. Dev. $982, 27, 301-308
n = number of samples Laplace transform operator T = sampling period t = time U = input variable X = output variable Xi,= actual or experimental value of output Xi,= predicted value of output 2 = 2 transform operator 2 , = modified 2-transform operator
s =
Greek Letters 8 = fractional dead-time Bd = dead-time T = time constant w, = natural frequency { = damping factor
Literature Cited Brantley, R. 0. M.S. Thesis, Unhrersity of Louisville, Louisville, KY, 1962. Deshpande, P. B.; Ash, R. H. “Elements of Computer Process Control with Advanced Control Applications”; Instrument Society of America: Research Traingle Park, 1981; Appendix A. Deshpande. P. B. AIChE J. 1980. 26. Dreifke, G. E.; Hougen, J. 0. Fourth Joint Automatic Control Conference, June 19, 1963.
301
Himmelblau, D. M. “Applied Nonlinear Programming”, McGaw-Hill: New York, 1972; p 44. Hougen, J. 0. “Experiences and Experiments with Process Dynamics”; CEP Monograph Series, American Institute of Chemical Engineers: New York, 1964; Vol. 60. Meyer, J. R.; Whitehouse, G. D.; Smith, C. L.; Murriil, P. W. Instrum. Control Syst. 1987,40, 76. Miller. J. A.; Lopez, A. M.; Smith. C. L.; Murriii, P. W. Contro/Eng. 1987. 14, 72. Moore, C. F.; Smith, C. L.; Murrlll, P. W. 64th Natlonal AIChE Meeting, New Orleans, LA, 1969. Murrili, P. W.; Pike, R. W.; Smith, C. L. Chem. Eng. 1989,76, 105. Oklenbourg. R. C.; Sartourius. H. “The Dynamics of Automatic Controls”; American Society of Mechanical Engineers: New York, 1948 p 276. Patke, N. G.; Deshpande, P. B. ISA Annual Conference, Nlagara Falls, NY, 1977. Schaefer, R. A. M.Eng. Thesis, University of Louisville, Louisville, KY, 1978. Schork, F. J.; Deshpande, P. B. Hydrocarbon Process. 1978,57, 113. Smith, C. L. “Digital Computer Process Control”; Intext Publishers: Scranton, 1972; Chapter 7. Smith, 0. J. M. ISA J. 195B,6,28. Sten, J. W. Instrum. Techno/. 1970, 17, 39. Sundaresan, K: R.; Krishnaswamy. P. R. Can. J. Chem. Eng. 1978, 56, 257. Sundaresan. K. R.; Prasad, C. C.; Krishnaswamy. P. R. I d . Eng. Chem. Process Des. D e v . 1978, 17, 237.
Received for review February 25, 1981 Accepted November 30,1981
Effects of Organic and Aqueous Phase Composition on Uranium Extraction from Phosphoric Acid with Octylphenyl Acid Phosphate Wesley D. Arnold,’ Donna R. McKamey, and Charles F. Baes, Jr. Oak Ridge National Laboratory, Oak Ridge, Tennessee 37830
Octylphenyl acid phosphate (OPAP) is a commercially available mixture of monooctylphenylphosphoric acid (MOPPA) and dioctylphenyl phosphoric acid (DOPPA) that is used to extract U’” from wet-process phosphoric acid. The uranium extraction coefficient and the organic phase iron loading were not affected by changes in the organic phase composition due to hydrolysis of the extractant or by different MOPPA/DOPPA ratios in two OPAP batches. The effect of aqueous phase composition was studied in a factorial design test and also by varying the concentrations of the components individually over a range typical of their concentrations in wet-process phosphoric acid. Uranium extraction was decreased by increasing the concentrations of total phosphate, ferric iron, hydrofluoric acid, or total sulfate, and increased by increasing the concentrations of aluminum, magnesium, or silicon. The extraction of iron was decreased by increasing the concentration of total phosphate but was not affected by changing the concentrations of the other aqueous phase components.
Introduction The wet-process phosphoric acid that is produced and converted to fertilizer each year in the United States contains about 2.7 Gg of UBOsthat is available for recovery as a byproduct, and this is projected to increase to more that 7 Gg by the end of the century (McGinley, 1972). About half the domestic phosphoric acid production capacity is now committed to uranium recovery; i.e., plants to recover the uranium are in operation, under construction, or being designed. The half that is not committed represents smaller plants (P205capacity less than 180 Gg/y) or plants that process phosphate rock with low uranium concentration, where uranium recovery is less attractive economically. A development program at Oak Ridge National Laboratory has included work to increase the efficiency of the uranium recovery processes so that they are more attractive for use in the smaller plants. The program has also included studies of uranium recovery from the more-concentrated phosphoric acids produced by
the hemihydrate processes. These processes produce calcium sulfate hemihydrate solids and phosphoric acid containing 40 to 50% Pz05,compared to calcium sulfate dihydrate and 30% P205acid prepared by the conventional processes. A t least four hemihydrate process versions are in commercial use or under development ( C h e m . Eng. News,1978; Goers, 1978; OrB, 1978; Schwer and Crozier, 1978) and they are attracting considerable interest from the phosphate industry because of the potential savings in capital and energy costs they offer. Uranium recovery from the hemihydrate process acids will be more difficult because of the higher concentrations of P205and impurities. Octylphenyl acid phosphate (OPAP), a mixture of mone and di(p-1,1,3,3-tetramethylbutyl)phenylphosphoric acids, is the extractant in one of the recovery processes developed at ORNL (Hurst and Crouse, 1974; Arnold et al., 1980). I t is a powerful extractant for uranium(1V) from phosphoric acid that could provide the desired efficiency for
0196-43051a211 121-0301t01.2510 0 1982 American Chemical Society
