QUASILINEARIZATION AND T H E ESTIMATION OF PARAMETERS IN DIFFERENTIAL EQUATIONS E. STANLEY LEE’ Research and Deuelopment Department, Phillips Petrole urn Co., Bartlesuille, Okla.
The estimation problem is formulated as a two-point or multipoint boundary-value problem. The classical least squares criterion is used. To illustrate the approach, the Peclet group and the reaction rate group in a tubular reactor with axial mixing are obtained.
N CHEMICAL engineering applications, it is frequently necessary
I to estimate the parameters or coefficients in systems of differential equations. Since generally the parameters cannot be measured directly, the measurable variables are the dependent variables of the differential equations, and the differential equations cannot be solved analytically in most cases, it is not a simple matter to identify these parameters. Estimating parameters and coefficients in systems of differential equations represents a particular problem in a large class of mathematical problems, which have been called the “inverse problem.” There are two classes of mathematical problems in engineering and physical sciences : to predict the behavior of a process with a given mathematical representation of the process and to determine the mathematical representation of the process from various given observations. Most of the mathematical analysis in the literature is concerned with the first problem. T h e inverse problem is much more complex and very little has been done in this difficult and yet challenging area of applied mathematics. This parameter estimation is treated as a two-point or multipoint boundary-value problem in this paper (Bellman et al., 1963, 1965). T h e resulting set of differential equations of boundary-value type problem is solved by the quasilinearization technique. Because of the effectiveness of this technique, the present approach forms a promising tool for estimating parameters in differential equations. Although the examples are based upon data which do not contain any noise, the technique can be extended to experimental data with noise and measurement errors. A Simple Estimation Problem T o illustrate the approach, consider the chemical reaction
A?A+B
(1)
which is taking place in a homogeneous tubular flow chemical reactor with axial mixing. Assume that the change of volume in the reactor is negligible. T h e following equation can be established easily by the use of material balance on reactant A
-1 -d2x _ -dx
P dt2
dt
- Rx2 =
0
where P is the dimensionless Peclet group Lv/D, R is the reaction rate group k L / v , x is the concentration of reactant A, and t is the dimensionless reactor length which varies between 0 and 1. t is obtained by dividing the actual position along the axial direction of the reactor by the total reactor length, L. Present address, Kansas State University, Manhattan, Kan. 152
ICLEC F U N D A M E N T A L S
T h e boundary conditions are (Danckwerts, 1954; Wehner and Wilhelm, 1956) (34
where x c represents the concentration of reactant A before entering the reactor and is a known quantity; and x(0) represents the concentration of A just after entering the reactor. There is a discontinuity of the concentration x of the reactant A a t the entrance to the tubular reactor, caused by axial mixing. Let us assume that although Equation 2 represents the tubular model, the two constant parameters, P and R, are unknown quantities and must be determined experimentally by measuring the concentration of the reactant, x , at various positions o f t . Or in other words, by x i n g the following measured or experimental values: dexp)
(ts) = bs, S = 1, 2,
,
. ., rn
(4)
with rn 2 2 and 0 5 ts 5 t,, we wish to determine the two unknown parameters. T h e quantities bs are known values and are obtained by measuring x experimentally at various positions of ts. T h e number of the experimental values, m, must be larger than or equal to the number of the unknown constant parameters. T h e superscript (exp) denotes that the values of x are experimental values. Since these parameters form part of the differential Equation 2, their values cannot be obtained easily from experimental data unless Equation 2 can be solved analytically. This kind of mathematical limitation frequently limits the manner in which an experiment can be performed. This parameter estimation is the same mathematical problem which appears in adaptive control process. A Multipoint Boundary-Value Problem
First, the case in which m = 2 is considered. It is assumed that the experimental errors resulting from obtaining b i and bz are small and thus bl and b2 represent the true values of x a t tl and t2. Equation 2 can be rewritten as: dx - =Y dt
(5)
2 = Py + PRx2 dt I t is convenient to consider the unknown parameters, P and R, as dependent variables parallel to x and as functions of the
independent variable I ! , Thus, in addition to the above two equations, the two constant unknown parameters can be represented by the following differential equation :
dR - -- 0 dt
(7)
Computational Procedure
T h e estimation problem can now be approached by the quasilinearization technique. However, for the purpose of simplicity, first a simpler form of this problem will be considered. Instead of Boundary Condition 9a, the following boundary condition will be used :
d-p-- 0 dt T h e boundary conditions for Equations 5 to 8 are: xe = x ( 0 )
-
1 -y(O), P
at t = 0
(94
y ( t j ) = 0 , a t t = tf at t = t~
(9c)
x(t2) = bz,
a t t = tz
(94
=
where c is a known value. Now the problem is composed of Equations 5 to 8, 9b, and 14. T h e other two boundary conditions can be obtained from either Equations 9c and 9d, or Equation 10. T h e systems of differential Equations 5 to 8 can be linearized easily by Recurrence Relation 12.
(9b)
bl,
x(tl)
(14)
at t = 0
x ( 0 ) = c,
xn+ I’ = Yn+ 1 ~n+l’=
+ Pnxn2Rn+1 +
+
P n y n + ~ 2 PnRaxnxn+l (rn+
Systems 5 to 9 form ,a multipoint boundary-value problem. tl and t z are two discrete values o f t within the interval [O, t f ] .
RnxnZ)Pn+1
-
R,+I’
0
=
(1 5 4
(3 PnRnx,’
+ P~Y,)
(15b)
Pn+l’ = 0 least Squares Approach
For m
> 2, the classical least squares criterion can be used.
T h e object is to determine the constant system parameters so that the sum of the squares of the deviations is minimized. Instead of the boundary conditions Equations 9c and 9d, one can obtain these two conditions by minimizing the following expression
T h e two given boundary conditions are:
l(0) =
c
yn+l(tf) =
0
xn+
T h e other two boundary conditions can be obtained either from
m
Q
[x(ts)
=
- bslZ
S=l
where the minimization is over parameters P and R , and x ( t s ) is obtained by solving E:quations 5 and 6.
x,+t(t1) =
61
(174
x,+i(tz)
bz
(17b)
=
or by minimizing m
Generalized Newton-R!aphson Formula
Q,+I
Equations 5 to 8 are nonlinear equations, which can be linearized by the generalized Newton-Raphson formula. Consider the general system of nonlinear differential equations
!?? = f * ( x l , x z , . . ., x y , t ) , i dt
= 1, 2,
. . ., M
(11)
=
S=l
[ ~ , + l ( t s )- b s I 2
(18)
Since the systems of Equations 15 are linear, they can be solved by the use of the principle of superposition (Lee, 1966). With one initial condition 16a given, only three sets of homogenenous solutions are needed for the four equations in 15. Thus, the general solutions of Equation 1 5 are:
with appropriate boundary conditions. Equation 11 can be linearized by the vector equation (Lee, 1966, 1968), dxn+l _ _ -f ( x n , t ) dt
+ J(xn)(xn+l - xn)
(12)
where x,+l, xn, and f represent M-dimensional vectors with component Xl,n+1, x?,n+1, . . ., X,u,n+G X l , n , XZ,n, . . ., X.w,n; andfl, f2, . , . ,&, respectively. T h e Jacobi matrix J ( x n ) 3
where subscript p is used to indicate particular solutions and subscripts hl, hz, and ha denote the first, second, and third sets of homogeneous solutions, respectively. UI,,+ 1, a’,,,+ 1, and u3,,,+ 1 are integration constants to be determined from the three boundary conditions. Equation 19 can be written in matrix form :
If we assume that x, is the known value and is obtained from previous calculations a:nd x,+1 is unknown, Equation 12 will always be linear.
xn+ l ( t ) =
xp(n+ 1) ( t )
+
Xh(n+ 1) (t)an+ I
(20)
T h e state vector x,+l(t) and the particular solution vector x ~ ( 1)~( t + ) are defined as: VOL. 7
NO. 1
FEBRUARY 1968
153
and a,+ I represents the integration constant vector with components a~,,+ I, a ~ >I, ~and + ~ 3 1. ,T h e~ homogeneous ~ solution matrix is defined as:
If only two experimental values, b l and bz, are given, the integration constants, u Z , ~ + ~and 0 3 , ~ + 1 , can be obtained from Equation 17. Substituting 17 into the first equation of 19, two algebraic equations are obtained. Combining these two equations with 28, the numerical values of the three integration constants can be obtained. If more than two experimental values are given, which is usually the case, Equation 18 can be used to obtain Q Z ,,+ 1 and ~ 3 , 1. ~ + By using the particular and homogeneous solutions a t various positions of t s , S = 1, 2, . . . , m, the following m equations can be obtained from the first equation of 19
S=l,2, T h e particular and homogeneous solutions will be chosen in such a way that they satisfy the given initial condition, 16a. T h e set of particular solutions can be obtained by integrating Equation 15 with the following initial values: xp,n+l(O)
=
C,
. . . ,m
(29)
First the integration constant al,,,+ 1 is eliminated from Equation 29 by using 28. T h e results from this elimination can then be substituted into Equation 18:
0, R p , n + l ( O ) = 0,
~ p , n + l ( O )=
Pp,n+l(O) =
0
(23)
T h e homogeneous forms of Equation 15 are: =
Xn+ I' Yn+l'
=
+2
PnYnt~
(244
yn+ 1
+
PnRnXnXn+l
Pn~n'Rnfl
(Yn
+
+ Rnx,2)Pn+I
(24b)
R,+i'
=
0
(24c)
Pn+l'
=
0
(24d)
Since all the particular and homogeneous solutions are known quantities, Equation 30 can be represented by:
T h e homogeneous solutions can be obtained by integrating Equation 24 with the following initial values:
T h e initial values, Equations 23 and 25, are chosen in such a way that a t t = 0 the first equation of 19 satisfies Condition 16a. At t = 0 the following two relationships can be obtained from the third and fourth equations of 19, and the initial values 23 and 25:
T h e only unknown quaxtities on the right-hand side of Equation 32 are the two integration constants, U Z , ~ + 1 and 1. I t can be shown by Equations 19 and 27 that the minimization of Expression 32 with respect to Rn+l ( t ) and Pn+~ ( t is) equivalent to the minimization of 32 with respect to u z,n+ 1 and a3 ,n+ I. This minimum can be obtained by various analytical or numerical methods. At present, the extreme value will be obtained by partial differentiation of 32 with respect to U Z , ~ + and ~ 3 , ~ + 1 .By these differentiations and by setting the results to equal to zero, the following two equations can be obtained: rn
Rn+l(0)
=
Pn+I(o)
u2,n+1,
=
a3,n+l
(26)
Since both R,+l(t) and P n + l ( t ) are constant functions, it is evident that Relationships 26 are true not only for t = 0, but also for 0 5 t 5 tl. Rn+l(t)
= aZ,n+b
Pn+I(t) =
a3,n+l
(27)
T h e particular and homogeneous solutions will be considered known and are determined computationally by the use of the initial values (Equations 23 and 25). T h e integration constant al,,+ can be expressed as a function of U Z , ~ + Iand U S , ~ + I . At t = t f , the following equation can be obtained by substituting Equation 16b into the second equation of 19 al,n+l
= -bp,n+~(t/)
+
+
a~,n+l~h2,n+l(t/)
03,n+ lYh3,n+ l ( t / ) I / Y h l . n + l ( t / )
154
l&EC FUNDAMENTALS
m
(28)
These two expressions are the other two boundary conditions. From Equation 33, the numerical values of U Z , ~ +1 and I can be obtained. T h e value of ul,,+l can then be obtained from Equation 28. Once the integration constants are known, the general solution of Equation 15 can be obtained from 19. With x n + l , y n + l , R n + l , and P n + l known, now a n improved 1 in 15. set of values can be obtained by making n = n T h e iterative procedure can be continued until the desired results are obtained.
+
~
Numerical Results
To test the effectiveness of this approach, some numerical experiments have been performed by using the results obtained by Lee (1966). These results (see Table I) are obtained by solving Systems 2 and 3 with the following numerical values: P=6,
R=2,
x,=1
(34)
T h e problem now is to find the unknown constant parameters P and R by using .Boundary Conditions 16 and the observations of x ( e x p ) ( t S ) , S = 1, 2, . , . , m, which are listed in Table I. I n other words, we are using the actual numerical solutions as the observations. T h e numerical values used are At = 0.01
t, == 1
m = 10
c = 0.83129
x o ( t t ) = 0.83129
~ ~ ( t l ; )
= -0.5
(35)
with At = t k + l - t i and k = 0, 1, 2, . . ., .V. T h e symbol At represents the integration step size. T h e known solutions of x listed in Table I are not used as the initial approximations. Using the Runge-Kutta integration method and the numerical values listed in 35, it has been found that the solutions will not converge to the correct values even if the exact values of P and R, which are 6 and 2 , respectively, are used as their initial approximations. This is caused by the very poor initial approximations used for x and J . This difficulty could have been avoided if better initial approximations for x and y had been used. T o avoid these difficulties, the first values of a? and 03, which are equal to R and P, respectively, according to Equation 27 are changed to within the following range before they are used in Equation 20 to obtain the general solutions
T h e purpose of these restrictions is to limit the values of a2 and and hence R and P, to within reasonable ranges. If these restrictions are used, the above convergence problems are not encountered. T h e influence of the initial approximations, RO(tK) and Po(tl;) k = 0 , 1, . , , , IT, upon the convergence rate of the constant parameters are listed in Table 11. Except for the different initial approximations for R and P, all other values are the same as those given by 35 and 36. Only the convergence rates for the two constant parameters are shown. I n general, the variables x and j have approximately the same convergence rates as the two constant parameters. For experimental purposes, the known numerical solutions for x and y are not used as the initial approximations in the above computations. If these known solutions, which are given in Table I, were used, the number of iterations needed would no doubt be reduced. T h e problem of using only two observations, or m = 2, is
also solved. Since this case is only a special case of the problem solved above, no additional discussions are necessary. Nonlinear Boundary Condition
Let us consider the original problem without simplifying Boundary Condition 9a. T h e system of equations is still represented by 15. However, the two given boundary conditions are : (374
Since none of the initial conditions are given explicitly, four sets of homogeneous solutions are needed. T h e general solutions for 15 can be represented by the following matrix equation :
where xn+l(t) and xpcn+l)(t)are the state and particular solution vectors, respectively, defined previously. T h e integration constant vector is now four-dimensional. T h e homogene) defined as ous solution matrix Xh(n+~ , ( t is
To obtain the one set of particular and four sets of homogeneous solutions, the following initial values will be used :
a3,
LO1
Numerical Values Used as Experimental Data
ts 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0.8 0.9 1 .o
X(ezP)(tg)=
0.50764 0.47017 0,43862 0.41305 0.39476 0.38727
bS
o o 01
Lo
o o
iJ
where I represents a unit matrix. By using the above initial values and the third and fourth equations of 38, again it can be shown that Rn+1 = a3,,,+ 1 and Pn+ I = a4?n+ 1. With the particular and homogeneous solutions known numerically, the integration constants can be obtained in the following manner. At t = 0, substituting the initial values of 40 and 41 into the first, second, and fourth equations of 38, the following relationships can be obtained :
Table II. Table 1.
rl
Convergence Rates of the Parameters
1 .o 1 .o
Rob)
PO(tk) Iteration 1
R(tk) 2.5923
2 3 4 5 6 7 8
2,7723 2,2235 2.0077 2.0000 2.0000 2.0000
1 .oooo
5 .O 5 .O p(fk)
1 .0000
5.1481 8.5792 6,0997 5.9817 5.9998 6.0000
6.0004
VOL. 7
NO. 1
R(tk) 1.0000 1.9317 2.0855 1 .9308 1.9981 1 ,9999 1.9999 2.0000
P(td 10.000 8.8433 3.9649 6.1562 5.9945 6.0010 6.0014 5.9996
FEBRUARY 1968
155
xn+
l(0)
+
= xe
(424
al,n+ 1
Yn+l(O)
= az,n+1
(42b)
Pn+l(o)
=
(424
a4,n+l
distribution is generated in the IBM 7094 computer by a random number subroutine. T h e numerical values used are summarized as follows:
X, al,n+ 1
= az,n+
(43)
l/(14,n+ 1
At t = t l , the second equation of 38 becomes: Yn+ I ( t / ) = yp,n+
l(tf)
aj,n+ lyhj,n+ l ( t f )
(44)
j=1
= bs,
~o(tJ =
(13
5
10,
15
(14
5
10
Besides Equation 37, which has been used to obtain Equations 43 and 45, two more boundary conditions are needed. These two conditions can be obtained by minimizing Expression 18. T h e computational results of the particular and homogeneous solutions a t various ts, S = 1, 2, . . , , m can be substituted into the first equation of 38. =
Xp,n+l(ts)
+
4 aj,n+1xhj,n+l(ts),
. . .,m
a numerical experiment is carried out and the results are listed in Table 111. T h e same fast convergence rate as that listed in Table I1 is obtained. Systems of Differential Equations
Consider the general system of nonlinear differential equations
(46)
By using Equations 42 and 45, the integration constants and U Z , ~ +1 can be eliminated from 46 :
i
~ ~ (= 0 )~
where 1
A = YhZ(tl)
+
(48) yhl(tl)/a4
+
For simplicity, the second subscript (n 1) has been omitted from the particular and homogeneous solutions and the integration constants in the above two equations. By substituting 47 and 48 into Equation 18, the desired expression is obtained. Symbolically, this desired expression can be written as : rn
[fS(a3,n+1,a4,n+I)
S=l
- bsIZ
(49)
T h e two integration constants, ( 1 3 , n f l and u4,,,+1, can be obtained by minimizing Expression 49. However, since P,+ l(0) appears nonlinearly in Boundary Condition 37a, function fs is also nonlinear in (14,n+l. If the differentiation approach were used to contain the extreme of 49, we would find that a4,,+ 1 is no longer single-rooted. I n the present case, it is not too difficult to find out which root of (14,n+1 minimizes 49. However, if a large number of the unknown constant parameters appear nonlinearly in the boundary condition, the problem of finding which set of roots is the desired one can be very time-consuming. T o avoid this difficulty and also to test some other approaches, a random search procedure has been used to find the minimum of 49. This procedure is described in detail by Lee (1964). T h e random number with Gaussian l&EC FUNDAMENTALS
= 1, 2,
. . ,, M
(53)
with boundary conditions
x r ( t f ) = x j f , j = 1, 2,
=
(52)
S=
U I , ~ +1
Qn+l
. . ., N
j-1
1, 2,
156
(51)
With the following starting values
RO(tk) = Po(tJ = 5 , k = 0, 1, 2,
xn+l(tS)
-0.5
. . ., ili. T h e 10 observed data points, S = 1, 2, . . ., m = 10 and ts = 0.1, 0.2, . . .,
15
=
1.0,
~ o ( t k )=
(50)
1.O are, again, obtained from Table I. To ensure convergence, the following restrictions on the allowable values of R and P, and hence of (13 and (14, are used
Combining Equations 37b, 43, and 44, one obtains: @,n+l
= 1.0,
m = 10
with k = 0, 1, 2, xexP(ts)
4
+
tf = 1,
At = 0.01,
Combining Equations 37a and 42 one gets:
~ 0q ,= M i
. . ., M I
(544
+ 1, . . ., M
(54b)
where P represents the constant parameters to be estimated. Thus, P corresponds to the Peclet group and the reaction rate group in the previous examples. I n most practical situations, only some of the dependent variables can be measured and parameters P cannot be measured directly. We assume that this is the case and only the first ml variables of x l ( t ) can be measured, 0 < ml < M . T h e following observations or measurements for the first m l variables of x are given
Since there are K unknown parameters, m m l > K. T h e equality sign can be used only if the measurement errors in obtaining b s ( j ) are small and thus bs(j) represent the true values of x,(ts). However, in actual applications this is seldom true. Thus, only the case with m m l > K will be considered. I n vector form, Equation 53 becomes dx -
dt
Table 111.
= f[x(t),
P]
(56)
Convergence Rates with Nonlinear Boundary Condition
Iteration
R(tk)
P(tk)
1 2 3 4 5
1,0000 1,8837 1.9705 1.9997 1.9997
9.9552 7.4656 6.0440 6.0056 6.0056
T o estimate the values of the unknown parameters, the following K differential equatiions can be formed
d- P= dt
o
yp(n+l)(0)
=
-
[O
0
. . .
XM,+10
XMO
-
0
.
*
*
O]T
K elements
M I elements
(57)
1 0 o 0 1 o 0 0 1
Now the unknown parameters P are considered as functions o f t and the problem can be stated as follows: Given a total of m m l observations on the first m l variables, find the K initial conditions P(0) for Equation 57. T h e least squares criterion again can be used and the problem becomes one of finding the K initial conditions so that the following expression is minimized mi
. ..
...o ...o . .
0 0 . .
0 0 . .
. .
. . .1 . . .o
. . . .
I t has been assumed that each of the m l measurable variables has a n equal amount of measurement errors. If the measurement errors are different for different variables, a weighting factor can be assigned it0 each variable and minimization can be performed according to the reliability of the data. T h e above problem (can be solved in the manner discussed above. First, Equation 56 can be linearized by the generalized Newton-Raphson formula. After the nth iteration, these linearized equations are
+
. . .
. . .
m
where y represents the ( M K)-dimensional vector with components X I , x2, . . ., x . ~ ,2'1, P z , . . . , PK and
. o
0 0 0 . .
0 . . .
0
. .
. o
0 0 0 . .
. . . . . .
. .o
. .
0 0 1 0 0 1
o...o o...o o...o
0 0
1 . . . 0
*
-
-
0 0 0 . .
0 0 v
M 1 rows
,
.
. . . .1
K rows
T h e above initial conditions are chosen in such a way that a t t = 0 Equation 61 satisfies the given initial conditions, Equation 54b. M I of the ( M I K) integration constants, an+1,can be obtained by using the given final conditions, Equation 54a. T h e other K integration constants can be obtained in the following manner. By substituting Equation 61 into Equation 58 a t various positions of ts, S = 1, 2, . . ., m, a n equation containing only a n f l as the unknown quantities can be obtained. M I of the ( M I K) unknown quantities can be eliminated by using Equation 54a. K equations can be obtained by differentiating this equation with respect to an+l and by setting the result to equal to zero. T h e other K integration constants can be obtained by solving these K equations. Once the integration constant vector is obtained, the results for the (n 1)st iteration can be obtained by using Equation 61 and the newly obtained particular and homogeneous solutions, 2)nd iteration can be obtained in exactly the T h e next (n same way. This procedure can be continued until no more improvement can be obtained.
+
+
+
T h e Jacobi matrix J ( y n , P,) is a n M x ( M K) matrix. Equations 57 and 59 together with the original given boundary conditions, Equation 54, constitute a linear boundary value problem. T h e K missing initial conditions can be obtained by minimizing Equation 58. T h e superposition principle can be used to solve these linear equations. Since ( M - MI) initial conditions are given, only ( M I K) homogeneous solutions are needed. T h e general solution can be represented by
+
Yn+l(t)
= YP(n+l)(t)
f Y h ( n + l ) ( t ) an+i
where y p represents the particular solutions, homogeneous solution matrix:
I
Xlhl,n+l
Yh(n+I)(t)
=
X.%lhl.nfl
Plhl.n+l
PKhl,n+l
Y h ( n + 1 )( t )
XlhZ.n+l
...
Xlh(.nii+K),n+l
X.MhZ.n+l
.. .
X.MMh(Mi+K),n+l
PlhZ,n+l
...
Plh(Mi+K),n+l
.......,. ......... PKh2,n+l
*
PKh(Mi+K),n+l
(61) is the
1 1
(62)
and a n f l is the ( M I 3- K)-dimensional integration constant vector. T h e particular and homogeneous solutions, y p ( n + l()t ) and Y h ( n + l ) ( t ) , can be obtained by using the following initial conditions
+
+
Discussion
T h e quasilinearization approach is a very useful tool for estimating the unknown system parameters. T h e parameters and the solution of the differential equations representing the system are obtained simultaneously. Furthermore, any known information about the parameters can also be utilized in estimating the initial approximations or starting values for the parameters. T h e latter aspect is especially important for the updating process in on-line computer control. For simplicity, the exact solution of the concentration x has been used as the experimental data. Obviously, the above procedure can also be used when the experimental data contain noise and measurement errors, provided that enough data are available. Furthermore, instead of the least squares criterion, any other criterion can be used for this estimation. If statistical data concerning noise and measurement errors were available, a criterion better than the least squares could be obtained. VOL 7
NO.
1
FEBRUARY 1968
157
Obviously, the above procedure can be applied to many chemical engineering problems. This approach is especially suited for on-line updating of parameters of a computer-controlled process. Another area of application is the estimation of chemical reaction rate constants from raw kinetic data with certain assumed or known mechanism of reaction. Since the estimation procedure is based on the quasilinearization technique, this approach is particularly suited for systems of differential equations with mixed boundary conditions. T h e above procedure can be generalized in various ways. I t can be used to estimate parameters or coefficients whose values are not constants but are functions of the independent variable, t. Instead of Equation 10, the following weighted criterion can be used
Nomenclature = integration constant = experimental value of x a t position t s = mean mass axial dispersion coefficient = Jacobi matrix, defined by Equation 13 = specific chemical reaction rate = total length of reactor = total number of experimental values of x =
= = =
= = = = = = =
dimensionless Peclet group, or unknown constant parameter defined by Equation 10 reaction rate group dimensionless reactor length, or independent variable final value o f t integration step size flow velocity of reaction mixture concentration of reactant A, or dependent variable concentration of reactant A before entering reactor experimental or measured value of x dependent variable
SUBSCRIPTS h = homogeneous solution n = nth iteration, assumed known 1 = (n 1)st iteration, assumed unknown n p = particular solution
+
+
literature Cited
-
and [xs bs]’ A[x, - b,] represents the quadratic form associated with the matrix A . The superscript T indicates the transpose of the vector within the bracket. T h e expansion of this quadratic form leads to a weighted sum of squares of the elements x - b, with the weighting determined by the elements of A . If A is a unit matrix, this quadratic form reduces to Equation 10. This generalized criterion has been used by Kalman (Kalman and Bucy, 1961).
Bellman, R., Kagiwada, H., Kalaba, R., “Quasilinearization, System Identification, and Prediction,” RAND Corp., Santa Monica, Calif., RM-3812-PR (August 1963). Bellman, R., Kalaba, R., “Quasilinearization and Nonlinear Boundary-Value Problems,” American Elsevier, New York, 1965. Danckwerts, P. V., Chem. Eng. Sct. 2, l(1953). Kalman, R . E., Bucy, R. S., J . Basic Eng. 83, 95 (1961). Lee, E. S., A.Z.Ch.E. J . 10, 309 (1964). Lee, E. S., Chem. Eng. Scz. 21, 183 (1966). Lee, E. S., “Quasilinearization and Invariant Imbedding,” Academic Press (in preparation), 1968. LVehner, J. F., FVilhelm, R. H., Chem. Eng.Scz. 6 , 89 (1956). RECEIVED for review December 5, 1966 ACCEPTEDAugust 8, 1967
LINEAR P R O G R A M M I N G 5-Transform Design of Digital
Controllersf o r Regulator Tvstems CARL I . H U B E R ’ A N D RICHARD I . KERMODEZ Department of Chemical Engineering, Carnegie-Adellon University, Pittsburgh, Pa. 152 13
Digital controllers for linear regulator control systems can readily be designed b y linear programming coupled with the Z-transform method. This method permits designs based on performance criteria and classes of deterministic input functions which are not suitable for conventional design techniques. The method is described in detail and illustrated by an example. N THE past 10 years a considerable amount of work has been l d o n e on the analysis and design of digital or sampled-data control systems. Articles by Koppel (1966a, 1966b), Min and Williams (1961), Mosler et al. (1966), and Williams and Min (1959) indicate that significant work is being done by chemical engineers in the design of sampled-data control sys-
1 Present address, E. I . du Pont de Nemours 8r Co., Wilmington, Del. 2 Present address, University of Kentucky, Lexington, Ky.
158
l&EC FUNDAMENTALS
tems. Most of the design methods for digital systems utilize the Z-transform method. This method can be used to reduce the linear differential equations describing a sampled-data system to linear algebraic equations. Kuo (1963), Ragazzini and Franklin (1958), and Tou (1959) have shown how to apply familiar design tools such as the Kyquist plot, Bode diagram, and root locus plot to the design of digital systems. Recent papers by Fegley (1964, 1965), Porcelli and Fegley (1964), Propoi (1963), Torng (1964), Whalen (1962), and Zadeh (1962) discuss various aspects of the linear programming
