= reactor or exchanger temperature, = condensate temperature, O K
O K
feed temperature, O K inlet temperature of heating fluid, O K time, sec over-all heat t p n s f e r coefficient, cal/ om2 sec . K = system inputs = reactor volume, cc = velocity, cm/sec = dimensionless feed composition = state variable = dimensionless reactor composition (Equation 13) = dimensionless reactor temperature (Equation 13) = functions in approximate solution = total exchanger length = exchanger length, cm = functions in approximate solution
= = = =
z Zi
GREEKLETTERS =
Y
damping coefficient (Equation 25)
I.1
= functions in approximate solution = small parameter
P
=
PH
= density of heating fluid, g/;c = dimensionless time (Equation 13)
80,
81, 82
r
density, g/cc
70 +l,
constant defined by Equation 25
= l/w,
=
42
= phase angles (Equation 25) = frequency of oscillation and natural
w , w*
frequency SUBSCRIPTS av S
= average = steady state
literature Cited
Douglas, J. M., Ind. Ens. Chem. Proc. Design Develop. 6 , 34 (1967). Douglas, J. M., Rippin, D. W. T., Chem. Eng. Sci. 21, 305 (1966). Gaitonde, N. Y., Douglas, J. M.,A.Z.Ch.E. J. 15, 902 (1969). Koppel, L. B., IND.ENG.CHEM.FUNDAMENT.4LS 1, 131 (1962). Minorsky, N., “Nonlinear Oscillations,” Van Nostrand, New York, 1962. Penrod, D. r).,Crandall, E. D., IND.ENG.CHEM.FUNDAMENTALS 5, 581 (1966). Peters, M.S., Timmerhaus!,K. D., “Plant Design and Economics for Chemical Engineers, 2nd ed., p. 308, McGraw-Hill, New York, 1968. Ray, W. H., IND.ENG.CHEM.FUNDAMENTALS 5, 138 (1966). Wilhelm, R. H., Rice, A. W., Bendelius, A. R., IND.ENG.CHEM. FUNDAMENTALS 5, 141 (1966). RECEIVED for review May 9, 1968 ACCEPTEDJuly 22, 1969
Determination of Dynamic Model Parameters Using Quasilinearization Brian 1. Ramaker,l Cecil 1, Smith, and Paul W. Murrill Department of Chemical Engineering, Louisiana State University, Baton Rouge, La. 70803
Quasilinearization was applied to determining numerical values of the parameters in a dynamic model. It can b e used to identify a dead time (time delay or transport lag) along with other parameters such as time constants and gains. The region in parameter space from which convergence will occur can b e increased by incorporating Marquardt’s method into this portion of quasilinearization. The model equation instead of the fundamental solutions can b e used to calculate the output and its derivatives for each iteration. The effect of data density and the number of significant figures in the data on the convergence of quasilinearization is discussed, and the performance of quasilinearization is compared to multidimensional search techniques used on similar problems.
THE
DEVELOPMENT of a model for a given system has become important in all fields of engineering. The possession of an adequate model facilitates implementation of control and optimization of the system and gives insight into its operation. T h e forms of the differential equations describing various systems often are fairly well documented; however, values of the parameters are usually specific to each system. These parameters generally cannot be measured directly but must be determined from experimental input-output records. The problem is to determine in some “best” sense the values of the parameters in the differential equations. The method of system identification called quasilinearization (Bellman and 1
28
Present address, Shell Oil Co., Deer Park, Tex. I&EC FUNDAMENTALS
VOL. 9 NO. 1 FEBRUARY 1970
Kalaba, 1965) treats this problem as a multipoint boundary value problem. Starting with an initial trial solution, the problem is formulated as a sequence of linear problems whose solutions converge to the solution, whether linear or nonlinear. This paper discusses the application of quasilinearization to the specific problem of determining the parameters in a dynamic model, which includes parameters related to the system time ConStantS, the gain, and the dead time (transportation lag or time delay). The objectives of this study are to show that quasilinearization can be used t o determine all model parameters including a dead time (transportation lag); to minimize the number of calculations per iteration by calculating the solution in terms of parameter
changes rather than the parameters themselves, which removes the necessity of calculating the particular solution, and by using the model equation instead of the fundamental solutions to update the output and its derivatives; and to investigate various methods to increase the region in parameter space from which convergence will occur (hereafter called the convergence region), including incorporation of Marquardt’s method (Marquardt, 1963) into quasilinearization, and taking only a portion of the calculated step. T h e specific identification problem is to determine the values of the parameters in the nth-order linear differential equation generally used for control purposes to describe system dynamics (including dead time) : an
dnx(t) dtn ~
+- + *
d 4 ) dt
+ z(t) = 9 m(t - e)
K
Z(t)
+ 6(t)
f(.t+l)
=
f(.i)
+
f’(.i)
(Q+l
-4
to Equation 4 we obtain:
+ Gtm(t - 0,) + + ( G ~ +-I Gt) m(t +
x‘t+i = -cizi
(Q+l
zt)(-ct)
et)
- Ct) (-4 + (Oi+i - e t ) Gt[-m’(t - e d l (Ci+l
(5)
T h e remaining differential equations which describe this problem are:
dci+i - dGi+i - doi+i dt dt dt
(1)
T h e a’s are related to the time constants of the process, g is the process gain, and is the dead time. This equation can be written in state variable form by choosing the n state variables to be z ( t ) and its derivatives: Z(t) =
value of the variable on the i t h iteration)
=
and are already linear. B y defining the 12 state variables to be the terms in the i 1s t solution
+
(2)
It is desired to determine parameters ai through a n , g, and 6 from discrete observations of only the first state variable, z ( t ) , and the manipulated variable, m(t). The argument t is omitted from the vectors from this point on. Since in the original problem parameters ai, a2, . . . a n , g, and e are unknowns, they can be treated as additional state variables, and since these state variables are constant with respect to time, they are described by the differential equations :
the equations are represented by the linear equation
and the
matrix and 6 vector become: -ci
0 0 0
T h e initial conditions for these equations are the respective values of the parameters. The number of state variables is increased from to 2n 2 = k . However, Equation 1 is now no longer linear with respect to the state variables. The initial conditions for all of the state variables are not known, but instead, a record of discrete input-output data is known. Therefore the parameter estimation must be formulated as a multipoint boundary value problem. The mathematical treatment and solution of linear boundary value problems are well understood. It is therefore desirable to formulate nonlinear boundary value problems in such a way that linear solution techniques may be used. Quasilinearization, which stems from Kantorovich’s extension of the wellknown Newton-Raphson procedure to function space, treats the nonlinear problem as the limit of a sequence of linear problems. T o establish notation and illustrate the method used by quasilinearization to linearize nonlinear differential equations, consider a first-order system described by the equation:
-xi 0 0 0 0
m(t
- ei) 0 0 0 0
-Gim’(t 0 0 0 0
- e,)
+
(4) where c = -
1 a1
9 G - = -
a1
The constants c, G, and 0 are to be determined from discrete observations on z ( t ) and m ( t ) . By application of the generalized Kewton-Raphson formula (subscript i denotes the
If the k state variables are defined in terms of the change of the state variables from one solution to the nest, the equations become homogeneous, making i t unnecessary to calculate a particular solution. T h e necessary equations are:
Matrix d is defined as before. During each iteration of quasilinearization it mill be necessary to determine numerically the solution to this set of linear first-order differential equations. Appearing in this set of equations is the derivative of the manipulated variable, and this is not directly available. However, the integrated form of the equation will not contain this derivative. The solution ? to a linear differential equation of the form
i=Iz+6 VOL. 9 NO. 1 FEBRUARY 1970
I&EC FUNDAMENTALS
29
1 ASSUME
) I+ - I
-
1
c---.
SOLUTION
I , NO
REGEN
CONVERGENCE
can be expressed in terms of homogeneous solutions j = 1,2, . . . , k , and the particular solution Q p as follows: IC
The coefficients ~ j , =j 1,2, . . . , k , would be the k initial conditions for the differential equation if available. But in this case, the coefficientsp 1 must be determined from the available experimental observations. Quasilinearization treats this as an overspecified multipoint boundary value problem. Thus, the coefficients are chosen to minimize the sum of squares of the deviations from the observed points, or mathematically, M
min 111,
*
Mn
N=l
[x.i+l(tN)
- xobs(tN)]
(9)
This reduces to the simultaneous solution of k linear algebraic equations, which can be applied with equal facility to either Equation 6 or 7 . The particular solution is not needed when the formulation in Equation 7 is used. Thus, the number of calculations are significantly reduced by defining the state variables in terms of changes instead of actual values. Programming Considerations
The computer program to determine the parameters in the dynamic model using the method of quasilinearization was written in Fortran IV for use on a 32K IBN-7040. U p to sixth-order differential equations of the general form in Equation 1 could be used for the model. One of the inherent problems with this method is the large demands placed on core storage for storing the fundamental solutions (Equation 7 ) . This can be circumvented by sacrificing computer running time to reduce core storage requirements, which was done in the program. For example, if the model is third-order, inclbding the dead time, and 200 data points are used with four points being calculated after each data point (giving 1000 points per solution) to improve the accuracy of the integration, the storage of the fundamental solutions would require 72K (1000 points, eight components of state vector, eight homogeneous solutions, and one particular solution). The 36
lbEC FUldDAMENTALS
VOL. 9 NO. 1 FEBRUARY 1970
program whose outline is shown in Figure 1 requires only 3072 storage locations for the same purpose, but the fundamental solution must be generated twice. Each of these blocks performs a basic operation required by quasilinearization and is written in subroutine form. T o initialize the iterations for quasilinearization requires assumptions of values for each state variable at every point in time corresponding to observed data points. The function of the subroutine Assume is to generate this initial solution for any given set of parameters and forcing function. The missing state variables, the output, and the necessary derivatives are calculated from the differential equation being used as the model. These values of the state variables must be stored for later use. The most time-consuming calculation during each iteration of quasilinearization is the numerical calculation of the fundamental solutions. The function of the subroutine Solution is to generate these fundamental solutions. T o minimize storage requirements, all the homogeneous solutions and the particular solution (if used) are generated simultaneously, with only the most recent values being retained. The calculations are kept to a minimum by taking advantage of the many zeroes that appear in the matrices. Since the complete solutions are not stored, it is necessary to calculate the fundamental solutions twice during each iteration. The first time the solutions are calculated the necessary information is generated to calculate the new parameter values via Equation 9. The subroutine RIatrix takes this information and forms k algebraic equations which are solved simultaneously in subroutine Gauss by Gauss reduction (using row interchange) for the new parameters. The second time the solution is generated the new parameters are used and subroutine Regen updates the remaining state variablesnamely, the output and its derivatives. By generating the solutions twice, considerable computer storage is saved b u t the calculation time is nearly doubled. (However, this was faster than using bulk storage.) The convergence is then checked to see if the parameters changed sufficiently during the last iteration to warrant going through the iterative process again. -411 numerical integration throughout this program was done using rectangular integration. T o increase the accuracy, the interval between observations was subdivided. The forcing function used in each subinterval was the linearly interpolated value between the observed values. Nelsa et al. (1968) report a survey of numerical integration methods on quasilinearization, including Newtonian, Gaussian quadrature, backward differences, and rectangular integration. Their conclusion was that rectangular integration is adequate. Identification of Dead Time
The termm(t - 0) which appears in the differential equations is calculated by looking back in time a period of e sampling periods. If the dead time is not an integral number of sampling periods, linear interpolation is used between the two nearest points. It is possible that on certain iterations as the solution is approached a negative dead time will appear. Although this situation is physically unrealizable for real systems, this apparently does not interfere with the convergence procedure. The convergence of qua.ilinearization is quadratic and very rapid when it does occur. Unfortunately, the parameter space for this particular problem in which convergence will occur is rather narrow, Table I shows the effect of increasing the order of the model on the convergence region. The data
Table Order
of Model
Order of Data
I.
Effect
c3
C4
of Order of M o d e l on Convergence c2
1
Correct parameters Starting point
1
Correct parameters Starting point
2
Correct parameters Starting point
3 3
Correct parameters Starting point
3
4
Correct parameters Starting point Correct parameters Starting point
4 4
4
2.5 2.0 2.5 2.4 table
Order of Model
Order of Data
3
Correct parameters Starting point
3
Correct parameters Starting point
3
3
c4
2.5 1.5 1.0 0 2.5 1.0 1.2 1.3 1.5 2.0 2.5 2.0 2.5 2.4
II.
2.5 0.9 1.0 2.5 1.5 1.0 0 2.5 1.0 1.2 1.3 1.5 2.0 2.5 2.0 2.5 2.4 Effect
Cl
G
2.5 0.6 0.8 1.0 2.5 0.9 1.0 2.5 1.5 1.0 0 2.5 1.0 1.2 1.3 1.5 2.0 2.5 2.0 2.5 2.4
2.5 0.6 0.8 1.0 2.5 0.9 1.0 2.5 1.5
0 0
1.0
0 2.5 1.0 1.2 1.3 1.5 2.0 2.5 2.0 2.5 2.4
2.5 1 .o 1.2 1.3 1.5 2.0 0 0 2.5 2.4
Significant Fig.
7 5
50 50 50
6 6 6
Yes
7
50 50
6 6
Yes Yes Yes
4 5 6
50
6 6 6
Iterations
No Yes Yes
No
No
50
50
Yes Yes Yes Yes
20 16 7 7
50 50 50 50 50
Yes
8
50
6
50
6
Significant Fig.
No
of Data on Convergence
c2
C1
G
2.5 1.2
2.5 1.2
2.5 1.2
2.5 1.2
2.5 1.5
2.5 1.5
2.5 1.5
2.5 1.5
1.0
1.0
1.0
1.0
c3
e 2.5 0.6 0.8 1 .o 2.5 0.9 1.0 0 0
No. of Data Points
Converged
for each test were generated from the appropriate order model equation using a damped sine as the forcing function. The correct parameters were therefore known and a perfect fit was possible. Table I shows that as the order of the model increases, the convergence region becomes smaller. This means that' the initial trial vector of parameters must be closer to the solution vector as the order of the model increases. Including dead time in the model, systems from first- to third-order were correctly identified. Fourth-order data using 50 data points were impossible to identify even when the trial solution was extremely close to the true solution. T h e order of the equation for which quasilinearizat,ion will be successful is limited by the fact that only one state variable-namely, s(t)-is observed and the other 12 - 1 must, be calculated. For this particular model equation includiiig dead time, t'hird-order appears to be the limit. Table I also indicates that the dead time parameter, e, is more difficult to determine than an additional dynamic parameter. The third-order system without dead time was correctly ideirtified down to a starting point of zero. The same sjstern with a dead time would not convergc a t a starting point of 1.0. On the other hand, in going from a second-
e 2.5 0
Iteration
Converged
No. of Data Points
5 5 5 11
Yes Yes Yes Yes No
50 40 30 20 10
Yes Yes Yes Yes Yes Yes
50 50 50 50 50 50
order system to a third-order system, very little difference was noted in the decrease of the convergence region. Higher than third-order systems-namely, fourth-were a1w correctly identified when the dead time was not included in the model. Effect of Data on Convergence
The convergence region for this specific problem is apparently affected very little by the number of data points used. Table I1 shows t h a t for a third-order system not including dead time 20 data points were sufficient. A fourthorder system including dead time could not be identified even when the number of data points was increased to 200. In the majority of tests 50 data points were used and if ('onvcrgerice did not occur, increasing the amount of data did not help. The accuracy of the final answer was a function of the number of data points and the accuracy of the data. Table I1 s h o w the effect of reducing the number of significwit figures. Conwrgewc \vas iiot :jtlversel\- :iffected hy the iiurir1)c.r of sigiiificaiit figures, siirce tlic saiiic Iiunilwr of iter:riioiis were required \Then 3, 2, or 1 sigiiificant figures were used. The VOL. 9 NO. 1 FEBRUARY 1970
l&EC FUNDAMENTALS
31
Table 111.
L 1 4 SOLUTION
MATRIX
Correct parameters Starting point General
GAUSS
Change ASSUME
PI = PZ = 0
Solution-Assume (PI = P2 = 0 )
‘7‘ STOP
Figure 2. zation
Solution-Assume program outline of quasilineari-
accuracy of the final answer was, however, affected by the number of significant figures. M i n i m i z i n g N u m b e r of Calculations
The method of quasilinearization has thus been demonstrated to determine the parameters in a dynamic model successfully (up to third-order including dead time) from input, output data. T h e convergence of this method, when it does occur, is quadratic and thus very rapid as far as the numher of iterations is concerned. IIowever, the number of calculations that must be performed during each iteration is extremely large. Because of the particular form of the differential equation for this problem and by utilizing knowledge of the physical sit’uation, the number of calculations can be reduced. Any reduction in the number of calculations increases the efficiency of the method as far as computer running time is concerned. The deletion of the particular solution by defining the state variables iii terms of the changes in the parameters was discussed above. For a second-order model not, including dead time this reduces the number of calculations in the solution of the differential equations by 20%. Since about 90% of the computer running time is required to generate the fundamental solutions, a 20% reduction represents a significant decrease. T h e necessity of generating all of the fundamental solutions the second time during each iteration to update the missing state variable can be circumvented by calculating them directly from the model equation-that is, by using subroutine Assume in Figure 1. This modification is shown in Figure 2. For a second-order system this means a saving of 35% in computer running time. For higher order systems the savings will even be greater. A second-order system not including dead time \vas correctly identified by the original general program in five iterations. This means that within the parameter space where convergence will occur the method of qiiasilinearization would be competitive with other techniques for solving the same problem, such as optimum seeking methods. However, tho utility of this method would be significantly enhanced if the ~‘oiivcrgeiicerc,gion could lie increased. The effect on the convergence regioii of the previously 32
l&EC FUNDAMENTALS
VOL 9 NO. 1 FEBRUARY 1970
Effect of Program Changes on Convergence c2
C1
g
2.5
2.5
2.5
5.0 12.0 20.0 5.0 12.0 20.0 5.0 12.0 5.0 12.0
5.0 12.0 20.0 5.0 12.0 20.0 5.0 12.0 5.0 12.0
5.0 12.0 20.0 5.0 12.0 20.0 5.0 12.0 5.0 12.0
Iterations
verged
5 6
Yes Yes
5 6
Yes Yes
5
Yes No Yes NO
Con-
No
No 5
described program changes was investigated using a secondorder system not including dead time. Fifty data points were generated as in the previous tests and three significant figures were used. Three initial starting points were chosen so that two (5 and 12) were within the parameter space of convergence for the original general program and the third (20) \vas outside. The results of this test are shown in Table 111. Defining the state variables in terms of the changes of the parameters had no effect on the convergence of quasilinearization. This was expected, since the same values are calculated a t each iteration for the two programs (Table 111). Although the same values are not obtained on successive iterations when subroutine Assume generates the new solutions, its effect on convergence did not appear to be significant. If it can be assumed that the state of the system is known before the forcing function IS introduced, the dimensionality of the problem can be reduced. Sormally for a second-order system, not i~icludiiigdead time, five values of p must be determined, corresponding to the output x(O), its derivative ~ ’ ( o )and , the constant parameters cl, cd, and G. If the initial state of the system is kno~vn,~ ( 0and ) ~ ’ ( 0are ) known, and only the constant parameters c,, cp, and G remain t o be determined. It was espected that a reduction in the dimensionality of the problem would enlarge the convergence region. However, Table 111 indicates that this is not the case. When parameters p1 and p2 were restricted to a value of zero, convergence was not achieved for a starting point of 12, whereas the general program was successful from this starting point. The same results were observed using subroutine Assume to generate the new solutions (Table 111). Increasing Convergence Region
Several methods have been mentioned in recent literature for circumventing the problem of convergence. Donnelly and Quon (1968) developed the niet,hod of data perturbation by which the data are altered until convergence occurs for the initially assumed solution. This solution is then used as the initial solution with the original data. If ronvergence does not occur, the data are prrturbed again and thc procedure js repeated until convergence occurs ~vithoutaltering the data. In this work, two other approaches were investigated to increase the convergence region, incorporation of Marquardt’s method (Marquardt, 1963) and reduction of step size. quasilinearizatioli as originally prol’osed uses Newton’s method csrlusively :inti the coiivergence is thiis limited to a11 area ncvir the trup solution. The mctho(1 of steepestt tlescriit, on the other hand, is kiiowi to work very \vel1 in regions far re-
Table IV. CP
Correct parameters Starting point hlarquardt’s method
2.5
Comparison of Convergence Using Marquardt’s Method 9
C1
2.5
5 20 20 50 50 50
5 20 20 50 50 50
100
100
5 20 20 50 50 50 100
(10)
= g
The values of p for the nest iteration were calculated directly from this equation in the original version of quasilinearization. If we let pt be the solution for p on the i t h iteration, then
pi+]
=
pi
+
Converged
Starting Value of X
Rate of Changing X per Iteration
2.5
moved from the true solution but very poorly in areas near the solution. XIarquardt’s method provides for an automatic transition from the method of steepest descent to Sewton’s method as the niininiuni of a function is approached. This method can be incorporated into quasilinearization if a least squares criterion is used. By the proper use of steepest descent and Sewtoil’s method, the convergence region should be increased. The solution of Equation 9 effectively reduces to the solution of the equation
zp
Iterations
fj+l
+
where is the change in p from the ith to the i lth iteration. Substituting into Equation 1 and rearranging gives:
Inserting hlarquardt’s modification into this equation gives (hlarquardt suggests that be nornialized before this is done) :
z
For X = 0 this equation reduces to the original version of quasilinearization, for which convergence is similar to that of Ken-ton’s method. For large values of A, the change f i + l will be in the direction of steepest descent. Using the test of convergence previously discussed, the effect of llarquardt ’s method was investigated. Various starting values of were tried along with various rates (including that proposed by Narquardt) of changing the value of X from iteration to iteration. As the solution is approached, the w l u e of is decreased to take advantage of both steepest descent and Sewton’s method in the regions where they work best. The efficiency of llarquardt’s method is thus greatly iiifluenced by both the starting value of X and the rate a t which it is varied. Table IV shows that a t a starting point of 20 the most efficient convergence was obtained with a starting X of 0.01 and by decreasing it by a factor of 10 during each iteration. A t a starting point of 50 convergence would not occur when X was varied between iterations. A t this starting point it was found best to hold X constant a t 0.001. At a starting point of 100 the search trajectory for steepest descent was found to be in the wroiig direction. Xewton’s method moved in the
5 7 16 25 36
Yes Yes Yes
0 0.01 10
NO
1
Yes Yes No
0.901 0.002
0.001
0 X/10 X/lO X/2
0 0 0
right direction but the step size was too large for convergence to occur. The problem of determining what value of X should be used and at, what rate it should be varied is not peculiar to the incorporation of hlarquardt’s method into quasilinearizatioIi. Even when llarquardt’s method is used as a n optimization technique by itself, the same problem exists. S o good general procedure is known for determining the value of A. However, this test demonstrated that llarquardt’s method docs increase the convergence region. Reduction of Step Size
Since the previous test indicated that Kewton’s method often gave correct search trajectories but that too large a step size prevented convergence, it would appear advantageous to take only a portion of the calculated step initially. The whole step could be utilized i n regions near the true solution. Two methods of varying the step size were tested: The step size was varied in a preassigned manner, and the step size was adjusted as part of the iteration cycle. The problem with the first method is that the poorer the assumed solution, the smaller the initial step size must be to assure coilvergelice. Any general function that could be used to adjust the step size for a wide variety of starting points would thus be rather inefficient for starting points close to the true solution. On the other hand, the m e of a function that is adjusted for each starting point’ will require several trials before the best one is found. To demonstrate that this method will work, a function was chosen which could be adjusted for various starting points and asymptotically approach 1. Portion of step size
=
1- (Iz)-~
B y adjusting the value of h the rate a t rvhich the function approached unity was changed. The second method, which adjusts the step size as part of the iteration cycle, does not have this problem and would therefore be more efficient. This method was programmed to check the sum of the squares and if no decrease was noted over the previous iteration the step size was cut in half until a decrease occurred or the step size was iiisignificarit, a t which point the program terminated. This procedure does not, guarantee the location of the opti~ i i u n :I.? i the search trajectory for S e n t o n ’ s method occasioiially is in a direction i i i which the criterion function increases. This could have been corrected, however, by moviiig i n the negative gradient direction when Sewton’s method failed to show an improvement. T h e test discussed previously was ueed to check the effect of t,hese two methods oil the colivergence region. Using the preassigned change in the step size always required several trials before the right function was found. -4s iiidicated in VOL. 9 NO. 1 FEBRUARY 1970
l&EC FUNDAMENTALS
33
Nomenclature Table V.
Comparison of Convergence Using Reduced Step Size e1
Correct parameters Starting point Portion of step size reassigned Portion of step size calculated
ct
g
Iteration
h
Converged
2 . 5 2.5 2.5 20 20 20 50 50 50 100 100 100 20 20 20 50 50 50 100 100 100 200 200 200
11 13 18 6 6
7 8
0.8 0.91 0.95
Yes Yes Yes Yes Yes Yes Yes
Table V, convergence was achieved for starting points of 20, 50, and 100. The second method, which calculates the step size internally, proved to be very efficient and was successful up to a starting point of 200. For this particular problem the search trajectory of Newton’s method was always in a direction in which the criterion function decreases. Summary
The method of quasilinearization identifies correctly the parameters in a dynamic model, including dead time, up to third-order from input-output data. The accuracy of the answer is affected by the number of significant figures in the data and can be increased b y increasing the amount of data. The number of calculations, and thus the computer running time, for this problem can be greatly reduced b y defining the state variables in terms of the change in the parameters and by using the model equation instead of the fundamental solutions to update the output and its derivatives. The convergence region for this problem was demonstrated to be increased by the incorporation of Alarquardt’s method into quasilinearisation or by using only a portion of the step size calculated in quasilinearization. The over-all efficiency of this method for this particular problem is thus probably comparable with optimization methods. Acknowledgment
The authors express their appreciation for support of this work by a Project T H E M I S contract administered for the U. S. Department of Defense by the U. S. Air Force Office of Scientific Research.
34
ILEC FUNDAMENTALS
VOL. 9 NO. 1 FEBRUARY 1970
parameters related to time constants of system coefficient matrix constant vector llal
dun steady-state system gain identity matrix iteration number 2%
+2
total number of observations level of manipulated variable a t time t order of model equation time constant vector for linear equations derived from Equation 9 vector of state variables level of system output a t time t fundamental solution of system equations particular solution of system equations coefficient matrix for linear equations derived from Equation 9
GREEKLETTERS e- = dead time in units of sampling time 6 = change in parameters = vector of weighting factors for fundamental solutions = change in ,iifrom iteration to iteration ( = scalar constant in hlarquardt’s method,, literature Cited
Bellman, R. E., Kalaba, R. E., “Quasilinearization and Nonlinear Boundary Value Problems,” American Elsevier, New York, 1965. Donnelly, J. K., Quon, D., “Computation of Best Fit Parameters in Non-Linear Mechanistic Rate Equations Using Quasilinearization and Data Perturbation,” Second Joint AIChE-IIQPR Meeting, Tampa, Fla., Preprint 19F, May 1968. Lee, E. S., “Quasilinearization, Difference Approximations, and Nonlinear Boundary Value Problems,” Sixtieth Annual A.1.Ch.E. Meeting, Kew York, N . Y., Preprint 36B, November 1967. Lee, E. S., IXD. ENG.CHEY.FUNDAMENTALS 7,152-7 (1968). Marquardt, D. W., J . SOC.Ind. A p p l . Math. 2 , No. 2, 431-41 (1963). Melsa, J. L., Pillmier, R., J.,.Bottorff, W. W., Steinway, W. J., “Research in the Ap licatlon of Modern Automatic,,Control Theory to Nuclear gocket Dynamics and Control, NASA N67-34462 (July 1968). RECEIVED for review Kovember 6, 1968 ACCEPTEDSeptember 22, 1969 Technical Report 6, LSU-T-TR-6. Research performed under Project THEMIS; AFOSR Contract F-44620-68-C-0021
