8270
Ind. Eng. Chem. Res. 2007, 46, 8270-8272
RESEARCH NOTES Identification of Integrating Processes with Deadtime and Inverse Response Anupam Srivastava† and Ashok K. Verma* Department of Chemical Engineering Technology, Institute of Technology, Banaras Hindu UniVersity, Varanasi-221005, India
A simple method for identifying parameters for integrating processes with deadtime and inverse response has been presented. Explicit expressions for the identification parameters have been derived. These expressions involve numerical integration. The present method is compared to another method (Luyben’s method). The present method gives more-reliable results and is free from the problems that are related to convergence. 1. Introduction Integrating processes with deadtime and inverse response have been observed in the case of adiabatic tubular reactors1 and boiler steam drums.2 The open-loop transfer function for such a process is given by2
GM(s) )
Kp(-τ1s + 1)e-Ds s(τ2s + 1)
(1)
The response for the process may appear as shown in Figure 1. The figure shows that the deadtime can easily be determined as the time until which the value of y does not change. Determination of the other parameters from the response (Kp, τ1, and τ2) involves the solution of highly nonlinear equations; hence, numerical procedures must be used. Luyben2 used MATLAB software to obtain these parameters from the knowledge of times t1 and t2, at which the response y from step test is zero, and ymin, which is the minimum value of y at t1. It involves numerical solution of the nonlinear equations. In this paper, a new explicit procedure to identify τ1, τ2, and Kp is obtained. The present method is compared to Luyben’s method.2
[ ( ) ]
τ1 + τ2 -t1/τ2 dy ) ∆uKp 1 e )0 dt τ2
(3)
τ2 ) (τ1 + τ2)e-t1/τ2
(4)
or
Equation 4 is used later to eliminate the nonlinearity in the identification procedure proposed in this work. At t ) t1, the response ymin is obtained. Substituting eq 4 into eq 2, we obtain
2. Development of the Identification Procedure The response for D ) 0 is given by
y ) ∆uKp[t - (τ1 + τ2) + (τ1 + τ2)e-t/τ2]
Figure 1. Identification parameters from the step test.
(2)
For D * 0, the time may be shifted by D, so that this condition is obtained. This observation implies that D must be determined first from the step response y, and the other three parameters should be obtained later. After shifting the time by D, the value of y has a minimum at t ) t1, which is obtained by equating the first derivative at t1 to zero:
ymin ) t1 - (τ1 + τ2) + (τ1 + τ2)e-t1/τ2 ∆uKp ) t1 - (τ1 + τ2) + τ2 ) t1 - τ1 or
τ1 ) t1 * To whom correspondence should be addressed. Tel.: 915422575993. E-mail address:
[email protected]. † Present address: Institut fu ¨ r Technische Thermodynamik und Thermische Verfahrenstechnik, Universita¨t Stuttgart, Germany.
ymin ∆uKp
(5)
Equation 5 may be used to estimate τ1 if Kp is known. One may note that, because ymin is negative (τ1 > t1).
10.1021/ie0614778 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/26/2007
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8271 Table 1. Comparison of the Present Method to Luyben’s Method for Estimating Kp, τ1, and τ2. Actual Values τ1
Luyben’s Method2
Present Method
Kp
t1
t2
ymin
τ1
τ2
Kp
τ1
1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 2 2 5 5 5 10 10 10
0.5 1 2 5 10 0.5 1 2 0.5 1 2 0.5 1 2
0.694 0.694 0.694 0.694 0.694 0.812 0.812 0.812 0.914 0.914 0.914 0.956 0.956 0.956
1.598 1.598 1.598 1.598 1.598 1.754 1.754 1.754 1.888 1.888 1.888 1.944 1.944 1.944
-0.3069 -0.1534 -0.6137 -1.5343 -3.0685 -0.0945 -0.1891 -0.3781 -0.0442 -0.0884 -0.1768 -0.0234 -0.0469 -0.0938
1.0037 1.0037 1.0037 1.0037 1.0037 1.0054 1.0054 1.0054 1.0036 1.0036 1.0036 1.0038 1.0038 1.0038
0.9965 0.9965 0.9965 0.9965 0.9965 1.9671 1.9671 1.9671 4.9615 4.9615 4.9615 9.8757 9.8757 9.8757
0.4954 0.9909 1.9817 4.9543 9.9086 0.4887 0.9774 1.9549 0.4931 0.9862 1.9724 0.4905 0.981 1.962
0.9986 1.0007 1.0302 1.0999 1.1612
2 2 2
1 1 1
0.5 1 2
1.105 1.105 1.105
2.835 2.835 2.835
-0.4507 -0.9014 -1.8028
2.0023 2.0023 2.0023
1.0156 1.0156 1.0156
0.5023 1.0045 2.0091
5
1
10
1.8
6.01
-32.0824
5.0063
1.0144
10.0062
5.1479
0.9332
9.4434
0.5
3.2
5
0.47
0.98
-0.1771
0.5285
2.0518
3.0282
0.6055
0.9637
1.3052
a
τ2
Values Estimated from eqs 12-14
τ2
Kp
1.007 0.5044 1.0033 1.1011 0.9559 1.8544 0.8332 3.8487 0.7258 6.5594 optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka optimizer is stucka
The optimizer is stuck at a minimum that is not a root.
∫0t
1
τ2 )
y dt
∆uKp
- (t12/2) + τ1t1 (7)
τ1 - t1
The values of t1 and t2 are known. The area under the curve (Figure 1) may be estimated numerically. To obtain Kp, let us consider the area under the curve from t ) 0 to t ) t2, which is given by
∫0t
2
y dt
∆uKp
)
t22 - (τ1 + τ2)t2 - (τ1 + τ2)τ2e-t2/τ2 + (τ1 + τ2)τ2 2 (8)
Figure 2. MATLAB code used to estimate Kp, τ1, and τ2.
At t ) t2, y ) 0; therefore, eq 2 gives the following condition: Let us integrate eq 2 to obtain the area under the curve (represented by the shaded portion in Figure 1) from t ) 0 to t ) t1.
∫0
t1
[
y dt ) ∆uKp
t2 - (τ1 + τ2)t - (τ1 + τ2)τ2e-t/τ2 2
]
t1
0
or Again substituting eq 4 into eq 6, the nonlinear terms can be
∫0t
1
y dt
∆uKp
)
t12 - (τ1 + τ2)t1 - (τ1 + τ2)τ2e-t1/τ2 + (τ1 + τ2)τ2 2 (6)
removed:
∫0t
1
t12 ) - (τ1 + τ2)t1 - τ22 + (τ1 + τ2)τ2 ∆uKp 2 y dt
t12 ) - τ1t1 + (τ1 - t1)τ2 2 or
(τ1 + τ2)e-t1/τ2 ) (τ1 + τ2) - t2
(9)
Substituting eq 9 into eq 8, we get
∫t
2
t2 1 0 y dt τ1 ) 2 t2 ∆uKp
(10)
Eliminating τ1 from eqs 5 and 10 and rearranging the formula yields the following expression for Kp:
[ ] ∫0t
2
Kp )
1 ∆u
ymin -
y dt
t2
t1 - (t2/2)
(11)
The integration ∫t02 y dt can be estimated numerically. Thus, all three parameters (Kp, τ1, and τ2) may be estimated, using eqs 11, 5, and 7, respectively. 3. Identification Procedure To obtain the values of parameters D, τ1, τ2, and Kp, the step response y, as a function of time, must be known. The procedure described in the following discussion considers the discrete data of y.
8272
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
(1) If the signal is noisy, use suitable filters. For scattered data, use interpolation methods to obtain a sufficient number of data. (2) From the step response, D is determined as the time until which y ) 0. The time is shifted by D. (3) Determine the values of ymin, t1, and t2 from the shifted response, as given in MATLAB code (see Figure 2). (4) Obtain the area under the curve from t ) 0 to t ) t2 numerically. Use eq 11 to estimate Kp. (5) Use eq 5 to estimate τ1. (6) Estimate the area under the curve from t ) 0 to t ) t1. Use eq 7 to estimate τ2. A smaller sampling frequency and a bad signal-to-noise ratio near the minimum of the inverse response will adversely affect the accuracy of the estimates of D, t1, and t2. As a result, the accuracy of the estimated values of Kp, τ1, and τ2 will also become worse. 4. Results and Discussion To evaluate the present identification method and compare it to Luyben’s method,2 the values of t1, t2, and ymin were determined for the process with known values of τ1, τ2, and Kp, using following equations:
( )
(12)
ymin ) (t1 - τ1)∆uKp
(13)
t1 ) -τ2 ln
τ2 τ1 + τ2
Equations 12 and 13 were obtained by rearrangement of eqs 3 and 5, respectively. Interestingly, note that t1 is independent of Kp. The value of ymin is proportional to Kp. Putting y ) 0 at t ) t2, eq 2 gives
t2 - (τ1 + τ2) + (τ1 + τ2)e
-t1/τ2
)0
(14)
After knowing the value of t1, eq 14 was solved using EXCEL’s SOLVER add-in to obtain the value of t2. It also may be estimated using MATLAB’s “fzero” function. Equation 14 has a root at t ) 0. It should be avoided by specifying the region of search excluding t ) 0. Equation 14 indicates that t2 also is independent of Kp. The inverse response to unit step input (∆u ) 1) was generated using eq 2 for a process with no delay time (D ) 0). The sampling rate of 100 Hz was determined to be sufficient. It is equivalent to errors of ∼1% and ∼9% in the estimation of times t1 or t2, if these were ∼1 and ∼10 s, respectively. The estimated values of y, and interval of samples, were stored in a file. The values of y were loaded from the file. From these values of y, the values of ymin and t1 were determined using MATLAB’s “min” function. The values of t2 were determined by noticing the change in the sign of y. Both the integrals were evaluated using a trapezoidal numerical integration method (MATLAB’s
“trapz” function). Luyben did not give the method that was used to determine the values of ymin, t1, and t2;2 therefore, the same values also were used to estimate the parameters Kp, τ1, and τ2, using Luyben’s method, by applying the code that he presented.2 To avoid the convergence problem, the number of iterations was increased from 1000 to 10 000 in some cases. The results are presented in Table 1. It was observed that Luyben’s method2 did not converge in many cases, whereas the present method is free from this problem, because the latter involved explicit equations. The numerical integration does not introduce any such problem. In the cases when it converged, the error, in some cases, was greater than that in the present case. For example, for values of τ1 ) 1, τ2 ) 1, and Kp ) 5, the calculated values using Luyben’s method were τ1 ) 1.1612, τ2 ) 0.7258, and Kp ) 6.5594. The values using present methods were similar to the actual values. Thus, the present method gives more-reliable results. The results may be affected in the presence of noise or in the case of a low sampling rate. Identification of the process for noisy data is under investigation. 5. Conclusions A simple identification procedure was proposed for integrating processes with reverse response and deadtime. Explicit expressions for Kp, τ1, and τ2 were derived. The expressions involved numerical integration. The present method gave morereliable results. As compared to Luyben’s method, in the present procedure there is no problem related to convergence. To obtain good estimates of the parameters, accurate values of t1, t2, and ymin were necessary. Nomenclature D ) deadtime (s) GM(s) ) process open-loop transfer function Kp ) process steady-state gain s ) Laplace transform variable t1 ) time at which the response reaches a minimum value (s) t2 ) time at which the response passes through the original steady-state value (s) y ) response of the step change in the process ymin ) minimum value of the response at t1 τ1, τ2 ) time constants (s) Literature Cited (1) Luyben, W. L. Tuning Proportional-Integral Controllers for Processes with Both Inverse Response and Deadtime. Ind. Eng. Chem. Res. 2000, 39, 973. (2) Luyben, W. L. Identification and Tuning of Integrating Processes with Deadtime and Inverse Response. Ind. Eng. Chem. Res. 2003, 42, 3030.
ReceiVed for reView November 20, 2006 ReVised manuscript receiVed October 10, 2007 Accepted October 15, 2007 IE0614778
