- Email: [email protected]
Closed loop Reaction Curve method for Identification of TITO systems
Third International Conference on Advances in Control and Optimization of Dynamical Systems March 13-15, 2014. Kanpur, India
Closed loop Reaction Curve method for Identification of TITO systems V. Dhanya Ram*. M.Chidambaram**
*Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai -600036, India (e-mail: [email protected]) **Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai -600036, India (e-mail: chidam@iitm .ac.in)
Abstract: A closed loop reaction curve method is presented to identify FOPTD (First Order plus Time Delay) model parameters of a Two Input Two Output (TITO) multivariable stable system which is controlled by decentralized PI/PID (Proportional Integral/ Proportional Integral Derivative) controllers. First, the responses and the interactions that are obtained for a step change in the set points separately are modeled by a SOPTD (Second Order plus Time Delay) models and integrated plus SOPTD models respectively. The parameters of these models are obtained in time domain by the extension of Yuwana and Seborg method originally proposed for SISO (Single Input Single Output) systems. From the identified closed loop transfer function matrix and the given controller transfer function matrix, the open loop transfer function matrix is obtained from the relation between the open loop and the closed loop transfer functions model. From the derived expressions for the open loop transfer functions, simplified FOPTD models are fitted in the s domain. The closed loop responses and interactions using the decentralized PI /PID controllers and the identified model are found to be matching satisfactorily with that of the original system. Further improved values of model parameters are obtained by using these identified model parameters as initial guess in the standard least square optimization method. For arbitrary values of the initial guess, the optimization method does not convergence. Keywords: model identification, closed loop, decentralized controllers, reaction curve method, optimization or transient response data. Yuwana and Seborg (1982) proposed an analytical method to calculate the model parameters of a FOPTD transfer function using closed loop response data. Kavdia and Chidambaram (1996) extended the method to an unstable FOPTD (First Order plus Time Delay) model by using the step response data of the closed loop system with a proportional controller. Ananth and Chidambaram (1999) extended this closed loop method for identifying the unstable FOPTD parameters using PI (Proportional Integral) or PID (Proportional Integral Derivative) controllers. Pramod and Chidambaram (2000, 2001) presented closed loop identification by an optimization method for stable and unstable SISO systems using the step response of PID-controlled systems. They reported a simple method for obtaining the initial guesses of the parameter values of the model transfer function. Padma Sree and Chidambaram (2002) extended this method to unstable transfer function model with a zero to estimate the model parameters by matching the closed loop step response of the actual process with that of the model. Rajapandiyan and Chidambaram (2012a) have proposed an optimization method for the closed loop identification of MIMO systems using decentralized controllers. A standard least square optimization method is used to obtain the parameters of the FOPTD models by matching the closed loop step response of the model with the actual process. Rajapandiyan and
1. INTRODUCTION Chemical industries often encounter multi-input multi-output (MIMO) systems. The presence of interactions among the loops in MIMO systems causes the control of these systems difficult. Also, it makes the identification of MIMO systems difficult. Transfer function model identification is important for the design of controllers. Open loop identification and closed loop identification are the widely available methods for identification. In open loop identification method, an excitation is given in each input variable one at a time and the corresponding output response is obtained. From the output responses, the transfer function is identified. Even though open loop identification is easier and simple, they are sensitive to disturbances and noise and also cannot be applied to unstable systems. Closed loop identification is insensitive to disturbances and measurement noises. An excitation of a known magnitude is applied to each set point one at a time and the output responses of the process are obtained. In closed loop identification, the main and the interaction responses are required for the identification of transfer function model. A better model fitting is obtained in the closed loop identification than in the open loop identification. Sundaresan and Krishnaswamy (1978) proposed four methods for estimating the dominant time constant and dead time of a given process from its moments, s-plane frequency 978-3-902823-60-1 © 2014 IFAC
989
10.3182/20140313-3-IN-3024.00135
2014 ACODS March 13-15, 2014. Kanpur, India
Chidambaram (2012b) have extended the method to the closed loop identification of SOPTD multivariable systems by using a combined step-up and step down responses. Melo and Friedly (1992) developed a frequency response characteristic using a closed loop procedure to identify MIMO systems. They converted the time domain responses and interactions of the closed loop system (controlled by decentralized Proportional controllers) by frequency response. Basically they have extended the frequency response technique for a SISO system proposed by Rajkumar and Krishnaswamy (1975) to a MIMO system. The frequency response characteristics (as Bode Plots) are presented. However, the method using P controllers introduces offset. They matched the frequency response of the identified model and the actual model. However no process model was explicitly obtained from their technique. The method was applied to identify stable systems. Ham and Kim (1998) proposed a closed loop identification method for multivariable systems having decentralized controller using a rectangular pulse change of set point. They replaced the Laplace parameter in the equation for the process transfer function with jω and for a given value of frequency the process transfer function was calculated. The model parameters are found out by using an optimization technique. The method was applied to identify only stable systems.
Fig.1. (b) TITO process with set point given to yr2 2.1 Identification of individual responses Consider a 2x2 transfer function matrix (G). kpij e-θ s G Gp12 where Gpi j= Gp= p11 τijs+1 Gp21 Gp22 ij
(1)
Let the controller settings be defined by the following matrix. G c=
Gc11 0
0 Gc22
where Gcij=k c ij(1+
1 +τDij s) τIij
(2)
The PI/PID controllers are selected so as to get a closed loop under damped response. Initial controller settings are taken based on the knowledge of the gain alone, using Davison’s method. Consider Fig.1 (a). A step change of known magnitude is given to the set point yr1 and the other set point is kept unchanged. The main response obtained is y11 and the interaction response is y21. Similarly refer Fig. 1 (b). A step change of the same magnitude is given as yr2 and the other set point is kept unchanged. The main response is denoted as y22 and the interaction as y21. yr1 0 Let the set point matrix be given as Yr = (3) 0 yr2 y11 y12 (4) and the output matrix be Y = y y22 21
In the present work, decentralized PI / PID controllers are used. The responses and interactions are modelled by the extension of Yuwana and Seborg (1982) method and the transfer function matrix for the closed loop system is obtained. Using the relation between the open loop and the closed loop transfer function matrices, the open loop transfer function is obtained in Laplace (s) domain. The responses of the obtained and the actual transfer functions matrix with the original controller settings are compared. Better results are obtained if the parameters of the identified transfer function models are used as initial guess for optimisation and the transfer function model is identified. In this study, the higher order models are approximated to a FOPTD model. Since the method is proposed for stable systems, many of the higher order stable systems can be approximated to a FOPTD system.
The transfer function of the closed loop responses is assumed to be an under-damped SOPTD system. From the closed loop response curve for each case, the value yp1, yp2, ym1 and T are noted (Refer Fig (9)). Using the formulas given by Yuwana and Seborg (1982), the value of τe and ζ are noted for each case. The closed loop time delay is noted from the responses. Due to the presence of the integral actions, the closed loop process gain kpc for the main responses (y11 and y22) are taken as 1. For the interaction responses (y21 and y12), kpc is found out by solving the Laplace inverse of y21 and y12 respectively.
2. PROPOSED METHOD
Now the responses can be approximated in the form Main response: yii (s)= Fig.1. (a) TITO process with set point given to yr1
kpcii e-θiis
s(τeii 2 s 2 +2ζii τeii s+1)
(5)
Interaction response: y (s)=
990
kpcij e-θijs (τeij 2 s2 +2ζij τeij s+1)
(6)
2014 ACODS March 13-15, 2014. Kanpur, India
The simulation is done for 15 sec with a sample time of 0.01. The transfer function of the main response is assumed as a SOPTD model.
2.2 To identify the transfer function for each case Friedly and Melo (1992) solved the equation for the transfer function of a closed loop response as G = YY (I − YY )
G#
Gpii (s) = (7)
$y12y21-y11y22%Gc22 s2 +$y11Gc22-y12Gc21%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
Gp21=
(8)
$y11y22-y12y21%Gc21 s2 +$y21Gc22-y22Gc21%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
Gp12 = Gp22 =
(9)
$y11y22 - y12y21%Gc12 s2 +$y12Gc11-y11Gc12%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12] $y12y12 - y11y22%Gc11 s2 +$y22Gc11-y21Gc12%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
(10) (11)
Gpij (s) =
Substituting the equations for y11, y12, y21, y22, Gc11, Gc21, Gc12, Gc22, the values of Gp11, Gp21, Gp12 and Gp22 are obtained as a function of s. A graph is plotted between Gp vs s for each case. By equating Gp11 to a FOPTD model, ie., Gp11=
kp11 e-θ11s τ11s+1
Gp12 = ( τ11s+1 Gp22 kp21 e-θ
kp11 e-θ11s 21s
τ21s+1
kp12 e-θ12s τ12s+1 kp22 e
-θ22s
τ22s+1
)
,
y11 (s)=
yii (s) =
1.807s+1
4.689 e-0.2 s
5.8 e-0.4 s
2.174s+1
1.801s+1
)
(13)
y22 (s)=
0.263+0.1852/s 0 0 0.163+0.09209/s
(17)
kpii (ps+1) e*θ++ s
-(τeii 2 s2 +2ζii τeii s+1)
(18)
(1+0.4448s)
s(0.4969s, +0.5057s+1)
1 s(0.4464s2 +0.3276 s+1)
y21 (s) =
The controller settings are designed using Davison’s (1976) method so as to get an under damped response Gc(s) =
e-0.2s
e-0.2s
(19)
For the response y22, the first peak is at 1.1162, but the first valley is at 0.895. This posed a problem in response identification as the obtained damping coefficient is much less. Hence, an approximate model is obtained by drawing a smooth curve from the first peak and taking around the midpoint of first valley and the steady state value 1, and reaching the second peak. Thus the responses are fitted as:
Consider the transfer function given by Chien et al (1999):
Gp(s)=(
1
Now y11 (s)=
3.1 Example 1
-11.64 e-0.4 s
s(0.4969s, +0.5057s+1)
The value of p can be obtained from the Laplace inverse of (18) which is given in the Appendix B. Substituting the values of kp11, τe, ζ, y, t, the value of p is obtained as 0.4448.
3. SIMULATION EXAMPLE
4.572s+1
(16)
A better matching of the response is obtained when the y11 is of the form of a SOPTD model with numerator dynamics as
(12)
The main and the interaction responses of the obtained transfer function matrix are compared with the actual process using the same controller settings. A better result can be obtained if these parameters are used as the initial guess and the transfer function model is identified through an optimization method.
22.89 e-0.2 s
kpcij s e-θijs 2 2 (τeij s +2ζijτeijs+1)
To obtain the value of kp21 and kp12, the Laplace inverse of (6) is taken. Appendix B gives the Laplace inverse of (6). Substituting the values of τe, ζ, y, t, in each case of interaction response, the values of kp21 and kp12 are obtained. The obtained responses can be written as:
suitable curve fitting method is applied to solve for kp11, τ11 and ϴ11. The procedure is repeated to Gp21,Gp12 and Gp22 in order to obtain kp21, τ21 and ϴ21, kp12, τ12 and ϴ12 , kp22, τ22 and ϴ22.Thus the transfer function matrix is identified as Gp11 Gp (s) = Gp21
(15)
The analytical method proposed by Yuwana and Seborg is used to identify the model (refer appendix A). For this the first and the second peak values, the first minimum valley and the period of oscillation are noted. From this the effective time constant (τe) and zeta (ζ) are obtained. Closed loop delay is taken same as the open loop delay (ϴ). Process Gain (kp) is found from the final steady state value. In order to identify y22, the above procedure is repeated. The value of kp11 and kp22 are taken as 1. Table 1 shows the observed peak values, the minimum valley, calculated period of oscillation, time constant, zeta and time delay for main and interaction responses. The response can be written in the form of (5). The transfer function of the interaction is assumed as‘s’ times of a SOPTD model.
Substituting Gp,y,yr,I,Gc in the matrix form and solving, the equations for Gp11,Gp21,Gp12 and Gp22 are obtained as below. Gp11 =
kpcii e*θ++ s
(τeii 2 s2 +2ζii τeii s+1)
y12 (s) =
(14)
991
e-0.4s
(20)
e-0.2s
(21)
(0.5087s2 +0.9896s+1) 0.3658
-0.5372 (0.4814s2 +0.3304s+1)
e-0.4s
(22)
2014 ACODS March 13-15, 2014. Kanpur, India
Gc(s) = 7
The fitted response and the actual responses are shown in Fig (2).The above four equations from (19) to (22) are substituted in (8) to (11) and a graph is plotted between each Gp vs s in each case. Fig.3 shows the Gp vs s graph of the identified responses. Each of the Gp is fitted to a FOPTD model, and the values of kp, τ and ϴ are found out in each case. For this 10 values of Gp for 10 values of s are taken in each case and are solved by using fmincon in MatLab. The obtained values of kp, τ and ϴ in each transfer function element are shown in Table 2. Fig. 4 shows the response curve for the identified model and the actual model using the same controller settings. From Fig.4, it can be seen that the responses are satisfactorily matching. The slight mismatching may be due to the approximation of the main and the interaction responses to a second order delay system. Even in SISO systems an accurate matching of the actual and the identified model is not reported. MIMO system identification encounters more problems because of the interactions from all the transfer functions. However, the model fitting can be further improved by using an optimization method such that the model parameters are selected to minimize the sum of the squared errors between the model and the actual process responses.
y11 (s)=
n- Number of the data points in the process and the model responses. The optimisation problem is solved by using the MATLAB routine lsqnonlin with the trust-region reflective algorithm. For the case of limits for optimization, the lower bound of kp (for positive values) is taken as 1/20th of the initial guess and the upper bound as 20 times the initial guess. For negative values of kp the lower bound of kp is fixed as 20 times the initial guess and the upper bound as 1/20th the initial guess value. For the case of time constant the lower limit is kept as 1/20th of the initial guess and upper limits as 5 times the initial guess. The lower bound for time delay is taken as 1/4th of the initial guess value and the upper bound as 4 times the initial guess value. The number of iterations taken is 31and the time taken for convergence is 36.97 sec. The sum of the ISE value obtained is 1.2865 X 10-8. The computational work is done on Intel Core i5 / 3.10 GHz personal computer. The converged parameters are shown in Table 3. Fig. 5 shows the identified and the actual responses.
6.6 e
10.9s+1
0.03 2s
-0.09s
8
(25)
1
e-s
(26)
s(22.7252s, +2.8556s+1) (1+1.0616s)
e-s
(27)
Similarly for main response y22, closed loop response obtained by Y-S method was found to be not matching with the actual response. Hence, by using fractional overshoot method, the obtained values are ζ =0.4331and τe =4.5475. The response can be written as y22 (s)=
1 s(20.6798s2 +3.939 s+1)
y21 (s) =
-18.9 e-3 s -19.4 e
s(26.639s, +1.441s+1)
Now y11 (s) =
Consider the transfer function given by Wood and Berry (1973): ) -3 s
−0.03-
0
0
A better matching of the response is obtained when the y11 is of the form of a SOPTD model with numerator dynamics as (18). p is obtained from the Laplace inverse of (18) given in Appendix B. Substituting the values of kp11, τe, ζ, y, t, the value of p is obtained as 1.0616.
3.2 Example 2
21s+1
+0.6s
A good matching of the identified model using Yuwana and Seborg method was not found with the actual response. It has been reported that Yuwana and Seborg method provides poor parameters for processes with large time delays. This is due to the use of a first order Pade approximation for time delay in the closed loop transfer function denominator (Jutan and Rodriguez, 1984). Hence effective time constant (τe) and zeta (ζ) are calculated by fractional over shoot method (Appendix C). The obtained values are ζ =0.4044 and τe =4.7671.
i= 1, 2; j=1, 2
Gp(s)=( 16.7s+1 -7 s
s
The simulation is done for 150 sec with a sample time of 0.01. The transfer function of the main responses is assumed as a SOPTD model of the form of (15) and the response model is assumed in the form of (5). The analytical method proposed by Yuwana and Seborg is used to identify the model response (refer appendix A). Table 4 shows the observed peak values, the minimum valley, calculated period of oscillation, time constant, zeta and time delay for main and interaction responses. The transfer function of the interaction is assumed as ‘s’ times a SOPTD model of the form of (16) and the interaction response model of the form (6). The value of kp11 and kp22 are taken as 1. To obtain the value of kp21 and kp12, the Laplace inverse of (6) is taken. Appendix B gives the Laplace inverse of (6). Substituting the values of τe, ζ, y, t, in each case of interaction response, the values of kp21 and kp12 are obtained. The obtained responses by using Yuwana and Seborg method can be written as:
Minimize ϕ = ∑526 ∑536 ∑4[0123 (tn) – yij(tn)]2 ∆t (23) (kpmij, τmij, ϴmij) ϕ = objective function
12.8 e- s
0.05
0.2+
(24)
y12 (s)=
14.4s+1
The controller settings are taken from Viswanathan et al (2001)
e-3s
4.4832 (32.0763s2 +2.5973s+1) 3.1235
(30.848s2 +3.6135s+1)
e-7s
e-3s
(28) (29) (30)
The fitted response and the actual responses are shown in Fig. (6).
992
2014 ACODS March 13-15, 2014. Kanpur, India
As mentioned earlier, equations from (27) to (30) are substituted in (8) to (11). A graph is plotted between each Gp vs s in each case. Fig.7 shows the Gp vs s graph of the identified responses. Each of the Gp is fitted to a FOPTD model, and the values of kp, τ and ϴ are found out in each case by using fmincon in MATLAB. The obtained values of kp, τ and ϴ in each transfer function element are shown in Table 5. When these obtained values are taken as initial guess for optimization, a good matching of the identified response with the actual response is obtained. The number of iterations taken is 14 and the time taken for convergence is 68.1011 sec. The sum of the ISE value obtained is 3.4131 X 10-7. The converged parameters are shown in Table 6. Fig. 8 shows the identified and the actual responses. The method can be extended to 3x3 systems or 4 x 4 systems. In these cases, the response models are identified for step changes in each set point one at a time.
Fig.2. Comparison plot of fitted response ( __ ) and actual response(----) (Example 1) (Closed loop reaction curve response)
Table 1. Main and Interaction responses (Example 1): Parameters
y11
y21
y12
y22
yp1 yp2 ym1 T τe ζ ϴ
1.369 1.034 0.891 4.745 0.704 0.358 0.2
0.331 0.0261 -0.124 4.7570 0.7132 0.3354 0.2
-0.5486 0.254 2.244 0.6938 0.2381 0.4
1.1162 0.9475 4.330 0.6681 0.2452 0.4
Table 2. Transfer function parameters obtained from Gp vs s curve (Example 1) Gp11 Gp21 Gp12 Gp22
kp 10.2890 5.1600 -15.1857 7.3514
τ 1.9820 2.6449 2.7003 1.2367
ϴ 0.2635 0.1773 0.3978 0.6447
Fig.3. Plot of Gp vs s curve (Example 1)
Table 3. Converged parameters from optimization using the parameters in Table 2 as initial guess (Example 1)
Gp11 Gp21 Gp12 Gp22
kp 22.8890 4.6888 -11.6401 5.8003
τ 4.5719 2.1749 1.8075 1.7998
ϴ 0.2002 0.2003 0.4001 0.4002
Fig.4. Comparison plot of the closed loop response of identified vs actual transfer function model for the same controller settings (Example 1)
993
2014 ACODS March 13-15, 2014. Kanpur, India
Fig.5. Response plot of the identified and actual transfer model using optimization ( _____ identified response; ------actual response ( both the responses over lap)) (Example 1)
Fig. 6. Comparison plot of fitted response ( __ ) and actual response(----)(Example 2) (Closed loop reaction curve response)
Table 4. Main and Interaction responses (Example 2): Paramete rs yp1 yp2 ym1
y11 1.2493 1.1095 0.8449
y21 0.5687 0.0916 -0.3205
y12 0.2826 0.0297 -0.1003
y22 1.212 1.08 0.9335
T τe ζ ϴ
32.75 5.1613 0.1396 1
36.56 5.6636 0.3393 7
36.9050 5.5541 0.3253 3
31.7 4.8853 0.2497 3
Table 5. Transfer function parameters obtained from Gp vs s curve (Example 2):
Gp11 Gp21 Gp12 Gp22
kp 11.2633 7.2743 -16.7834 -20.4995
τ 24.2644 24.0863 19.8797 12.7444
Fig.7. Plot of Gp vs s curve (Example 2)
ϴ 2.8368 7.4754 2.9102 4.5955
Table 6. Converged parameters from optimization using the parameters in Table 2 as initial guess (Example 2): kp
τ
ϴ
Gp11
12.7565
16.6364
1.0011
Gp21
6. 5759
10.8514
7.0034
Gp12
-18.8176
20.9070
2.9975
Gp22
-19.3505
14.3516
3.0013
Fig.8. Closed loop response plot of the identified and actual transfer matrix using optimization ( _____ identified response; ------- actual response ( both the responses over lap)) (Example 2)
994
2014 ACODS March 13-15, 2014. Kanpur, India
Chien, I. L., Huang, H. P., Yang, J. C.(1999).A Simple Multiloop Tuning method for PID Controllers with No Proportional Kick. Ind. Eng. Chem. Res., 38, 1456-1468. Chidambaram, M. (1998), Applied Process Control, Allied Publishers, New Delhi, India. Davison, E.J. (1976). Multivariable tuning regulators: The feed-forward and robust control of general servomechanism problem. IEEE Trans Auto. Contr, AC-21, 35-41. Ham,T.W. and Kim,Y.H.(1998). Process Identification and PID Controller Tuning in Multivariable Systems. Journal of Chemical Engineering of Japan, 31(6), 941949. Jutan, A., and Rodriguez, E.S. (1984). Extension of a new method of on-line controller tuning. Can. J. Chem. Eng., 62, 802-807. Kavdia, M. and Chidambaram, M. (1996). On-Line Controller tuning for unstable systems. Computers chem.Engng, 20(3), 301-305. Melo,D.L. and Friedly,J.C, (1992).On-Line, Closed-Loop identification of Multivariable Systems. Ind.Eng.Chem.Res., 31, 274-281. Padma Sree,R .and Chidambaram,M. (2002). Identification of unstable transfer model with a zero by optimization method. J.Indian Inst.Sci., 82, 219-225. Pramod,S. and Chidambaram, M. (2000). Closed loop identification of transfer function model for unstable bioreactors for tuning PID controllers. Bioprocess Engineering, 22,185-188 Pramod,S.and Chidambaram,M. (2001).Closed loop identification by Optimization method for tuning PID controllers. Indian Chem.Engr., 43 (2),90-94. Rajakumar,A. and Krishnaswamy, P, R. (1975). Time to Frequency Domain Conversion of Step Response Data. Ind. Eng. Chem., Process Des. Dev., 14 (3), 1975. Rajapandiyan, C.and Chidambaram, M. (2012a). Closed loop Identification of Multivariable systems by Optimization method. Ind. Eng. Chem. Res., 51, 1324- 1336. Rajapandiyan,C. and Chidambaram,M. (2012b). ClosedLoop Identification of Second-Order Plus Time Delay (SOPTD) Model of Multivariable Systems by optimization Method. Ind. Eng. Chem. Res., 51, 96209633. Sundaresan,K.R. and Krishnaswamy,P.R. (1978), Estimation of Time Delay Time Constant Parameters in Time, Frequency, and Laplace Domains. The Canadian Journal of Chemical Engineering, 56, 257-262. Pradeep K. Viswanathan, Wei Khiang Toh, and G. P. Rangaiah (2001), Closed-Loop Identification of TITO Processes Using Time-Domain Curve Fitting and Genetic Algorithms. Ind. Eng. Chem. Res., 40, 28182826. Wood, R. K.; Berry, M. W. (1973), Terminal Composition Control of a Binary Distillation Column. Chem. Eng. Sci., 28, 1707-1717. Yuwana, M. and Seborg, D.E. (1982). A new method for online controller tuning. AIChE J., 28,434-440.
4. CONCLUSIONS The closed loop reaction curve method for SISO systems proposed by Yuwana and Seborg (1982) is extended to identify multivariable systems. The problem associated with this method is brought out. Simulation application is given for transfer function matrix of 2x2 multivariable systems. Two examples are given. The method involves first by fitting the main actions and the interactions of the closed loop system for step changes in the set points by a SOPTD transfer functions with (1/s) term and the interaction by SOPTD model. From the analytical expressions, the open loop transfer function matrices are obtained using the closed loop transfer function matrix and the controller transfer function matrix by plotting Gp vs s curve. The obtained model is used with the same PID controllers and the responses/interactions are matched with that of the actual systems. A method to improve the matching involves by taking the approximate model parameters as the guess values of the least square optimization method. Any arbitrary guess values do not give convergence of the optimization problem.
5. NOMENCLATURE Gp(s) Gc(s) gp, gc, gm kpm, kp, kpc τm, τ θm, θ yr,ymr Y y,ym yp1 yp2 ym1 T τe ζ t kc τI τD ∆t FOPTD SOPTD ISE
process transfer function matrix controller transfer function matrix process, controller and model transfer functions model and process steady state gains closed loop process gain model and process time constant model and process time delay process and model set point closed loop response matrix closed loop process and model response first peak in the closed loop response (refer Fig.9) second peak in the closed loop response (refer Fig.9) first valley in the closed loop response(refer Fig.9) period of oscillation (refer Fig.9) effective time constant damping coefficient time controller gain Integral time constant Derivative time constant sampling time first order plus time delay system second order plus time delay system Integral of squared error REFERENCES
Ananth, I. and Chidambaram, M. (1999). Closed-loop identification of transfer function model for unstable systems. Journal of the Franklin Institute, 336, 10551061. 995
2014 ACODS March 13-15, 2014. Kanpur, India
Appendix A. FIRST APPENDIX
Fig. 9 shows the closed loop response of an under damped process The expression for the effective time constant (τ) and the damping coefficient (ζ) can be derived as: 9∞
9:;
9<;
9∞
9<;
v 1=
9∞
9<,
v 2=
ζ =
ζ5 =
9∞
(
KL =
(A-2)
=4 (>; )
[?, @A=4(>; )B, ]D.F
(A-3)
=4 (>, )
[G?, @A=4(>, )B, ]D.F
ζ = ζ1 + ζ2 T=
(A-1)
5?IJ
M(
,
(A-4) (A-5)
ζ )D.F
(A-6)
,
ζ )D.F
5?
(A-7)
Appendix B. SECOND APPENDIX y(s)=
(1+ps)e-θs
s(τe , s2 @2ζτe s+1
(B-1)
Laplace inverse of (B-1) is
y(t)=A [sinbexp( )]B+{1- expN O [qζ sin(b) + cos(b)]B pq
-ζt’
-ζt’
τe
τe
τe
(B-2)
y(s)=
kp e-θs 2 τe s +2ζτe s+1
(B-3)
Laplace inverse of (B-3) is y(t)=
kp q τe
[exp( )sinb] -ζt’ τe
where t’=t-ϴ ; b=(1-ζ2)0.5
(B-4) RS
IJ
; q=(1-ζ2) -0.5
Appendix C. THIRD APPENDIX Max Overshoot final value
τeff =
=e
TT1-ζ2 2π
πζ
T1-ζ2
(C-1) (C-2) 996
Closed loop Reaction Curve method for Identification of TITO systems V. Dhanya Ram*. M.Chidambaram**
*Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai -600036, India (e-mail: [email protected]) **Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai -600036, India (e-mail: chidam@iitm .ac.in)
Abstract: A closed loop reaction curve method is presented to identify FOPTD (First Order plus Time Delay) model parameters of a Two Input Two Output (TITO) multivariable stable system which is controlled by decentralized PI/PID (Proportional Integral/ Proportional Integral Derivative) controllers. First, the responses and the interactions that are obtained for a step change in the set points separately are modeled by a SOPTD (Second Order plus Time Delay) models and integrated plus SOPTD models respectively. The parameters of these models are obtained in time domain by the extension of Yuwana and Seborg method originally proposed for SISO (Single Input Single Output) systems. From the identified closed loop transfer function matrix and the given controller transfer function matrix, the open loop transfer function matrix is obtained from the relation between the open loop and the closed loop transfer functions model. From the derived expressions for the open loop transfer functions, simplified FOPTD models are fitted in the s domain. The closed loop responses and interactions using the decentralized PI /PID controllers and the identified model are found to be matching satisfactorily with that of the original system. Further improved values of model parameters are obtained by using these identified model parameters as initial guess in the standard least square optimization method. For arbitrary values of the initial guess, the optimization method does not convergence. Keywords: model identification, closed loop, decentralized controllers, reaction curve method, optimization or transient response data. Yuwana and Seborg (1982) proposed an analytical method to calculate the model parameters of a FOPTD transfer function using closed loop response data. Kavdia and Chidambaram (1996) extended the method to an unstable FOPTD (First Order plus Time Delay) model by using the step response data of the closed loop system with a proportional controller. Ananth and Chidambaram (1999) extended this closed loop method for identifying the unstable FOPTD parameters using PI (Proportional Integral) or PID (Proportional Integral Derivative) controllers. Pramod and Chidambaram (2000, 2001) presented closed loop identification by an optimization method for stable and unstable SISO systems using the step response of PID-controlled systems. They reported a simple method for obtaining the initial guesses of the parameter values of the model transfer function. Padma Sree and Chidambaram (2002) extended this method to unstable transfer function model with a zero to estimate the model parameters by matching the closed loop step response of the actual process with that of the model. Rajapandiyan and Chidambaram (2012a) have proposed an optimization method for the closed loop identification of MIMO systems using decentralized controllers. A standard least square optimization method is used to obtain the parameters of the FOPTD models by matching the closed loop step response of the model with the actual process. Rajapandiyan and
1. INTRODUCTION Chemical industries often encounter multi-input multi-output (MIMO) systems. The presence of interactions among the loops in MIMO systems causes the control of these systems difficult. Also, it makes the identification of MIMO systems difficult. Transfer function model identification is important for the design of controllers. Open loop identification and closed loop identification are the widely available methods for identification. In open loop identification method, an excitation is given in each input variable one at a time and the corresponding output response is obtained. From the output responses, the transfer function is identified. Even though open loop identification is easier and simple, they are sensitive to disturbances and noise and also cannot be applied to unstable systems. Closed loop identification is insensitive to disturbances and measurement noises. An excitation of a known magnitude is applied to each set point one at a time and the output responses of the process are obtained. In closed loop identification, the main and the interaction responses are required for the identification of transfer function model. A better model fitting is obtained in the closed loop identification than in the open loop identification. Sundaresan and Krishnaswamy (1978) proposed four methods for estimating the dominant time constant and dead time of a given process from its moments, s-plane frequency 978-3-902823-60-1 © 2014 IFAC
989
10.3182/20140313-3-IN-3024.00135
2014 ACODS March 13-15, 2014. Kanpur, India
Chidambaram (2012b) have extended the method to the closed loop identification of SOPTD multivariable systems by using a combined step-up and step down responses. Melo and Friedly (1992) developed a frequency response characteristic using a closed loop procedure to identify MIMO systems. They converted the time domain responses and interactions of the closed loop system (controlled by decentralized Proportional controllers) by frequency response. Basically they have extended the frequency response technique for a SISO system proposed by Rajkumar and Krishnaswamy (1975) to a MIMO system. The frequency response characteristics (as Bode Plots) are presented. However, the method using P controllers introduces offset. They matched the frequency response of the identified model and the actual model. However no process model was explicitly obtained from their technique. The method was applied to identify stable systems. Ham and Kim (1998) proposed a closed loop identification method for multivariable systems having decentralized controller using a rectangular pulse change of set point. They replaced the Laplace parameter in the equation for the process transfer function with jω and for a given value of frequency the process transfer function was calculated. The model parameters are found out by using an optimization technique. The method was applied to identify only stable systems.
Fig.1. (b) TITO process with set point given to yr2 2.1 Identification of individual responses Consider a 2x2 transfer function matrix (G). kpij e-θ s G Gp12 where Gpi j= Gp= p11 τijs+1 Gp21 Gp22 ij
(1)
Let the controller settings be defined by the following matrix. G c=
Gc11 0
0 Gc22
where Gcij=k c ij(1+
1 +τDij s) τIij
(2)
The PI/PID controllers are selected so as to get a closed loop under damped response. Initial controller settings are taken based on the knowledge of the gain alone, using Davison’s method. Consider Fig.1 (a). A step change of known magnitude is given to the set point yr1 and the other set point is kept unchanged. The main response obtained is y11 and the interaction response is y21. Similarly refer Fig. 1 (b). A step change of the same magnitude is given as yr2 and the other set point is kept unchanged. The main response is denoted as y22 and the interaction as y21. yr1 0 Let the set point matrix be given as Yr = (3) 0 yr2 y11 y12 (4) and the output matrix be Y = y y22 21
In the present work, decentralized PI / PID controllers are used. The responses and interactions are modelled by the extension of Yuwana and Seborg (1982) method and the transfer function matrix for the closed loop system is obtained. Using the relation between the open loop and the closed loop transfer function matrices, the open loop transfer function is obtained in Laplace (s) domain. The responses of the obtained and the actual transfer functions matrix with the original controller settings are compared. Better results are obtained if the parameters of the identified transfer function models are used as initial guess for optimisation and the transfer function model is identified. In this study, the higher order models are approximated to a FOPTD model. Since the method is proposed for stable systems, many of the higher order stable systems can be approximated to a FOPTD system.
The transfer function of the closed loop responses is assumed to be an under-damped SOPTD system. From the closed loop response curve for each case, the value yp1, yp2, ym1 and T are noted (Refer Fig (9)). Using the formulas given by Yuwana and Seborg (1982), the value of τe and ζ are noted for each case. The closed loop time delay is noted from the responses. Due to the presence of the integral actions, the closed loop process gain kpc for the main responses (y11 and y22) are taken as 1. For the interaction responses (y21 and y12), kpc is found out by solving the Laplace inverse of y21 and y12 respectively.
2. PROPOSED METHOD
Now the responses can be approximated in the form Main response: yii (s)= Fig.1. (a) TITO process with set point given to yr1
kpcii e-θiis
s(τeii 2 s 2 +2ζii τeii s+1)
(5)
Interaction response: y (s)=
990
kpcij e-θijs (τeij 2 s2 +2ζij τeij s+1)
(6)
2014 ACODS March 13-15, 2014. Kanpur, India
The simulation is done for 15 sec with a sample time of 0.01. The transfer function of the main response is assumed as a SOPTD model.
2.2 To identify the transfer function for each case Friedly and Melo (1992) solved the equation for the transfer function of a closed loop response as G = YY (I − YY )
G#
Gpii (s) = (7)
$y12y21-y11y22%Gc22 s2 +$y11Gc22-y12Gc21%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
Gp21=
(8)
$y11y22-y12y21%Gc21 s2 +$y21Gc22-y22Gc21%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
Gp12 = Gp22 =
(9)
$y11y22 - y12y21%Gc12 s2 +$y12Gc11-y11Gc12%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12] $y12y12 - y11y22%Gc11 s2 +$y22Gc11-y21Gc12%s
&$y11y22-y12y21%s2 -(y11+y22)s+1'[Gc11Gc22-Gc21Gc12]
(10) (11)
Gpij (s) =
Substituting the equations for y11, y12, y21, y22, Gc11, Gc21, Gc12, Gc22, the values of Gp11, Gp21, Gp12 and Gp22 are obtained as a function of s. A graph is plotted between Gp vs s for each case. By equating Gp11 to a FOPTD model, ie., Gp11=
kp11 e-θ11s τ11s+1
Gp12 = ( τ11s+1 Gp22 kp21 e-θ
kp11 e-θ11s 21s
τ21s+1
kp12 e-θ12s τ12s+1 kp22 e
-θ22s
τ22s+1
)
,
y11 (s)=
yii (s) =
1.807s+1
4.689 e-0.2 s
5.8 e-0.4 s
2.174s+1
1.801s+1
)
(13)
y22 (s)=
0.263+0.1852/s 0 0 0.163+0.09209/s
(17)
kpii (ps+1) e*θ++ s
-(τeii 2 s2 +2ζii τeii s+1)
(18)
(1+0.4448s)
s(0.4969s, +0.5057s+1)
1 s(0.4464s2 +0.3276 s+1)
y21 (s) =
The controller settings are designed using Davison’s (1976) method so as to get an under damped response Gc(s) =
e-0.2s
e-0.2s
(19)
For the response y22, the first peak is at 1.1162, but the first valley is at 0.895. This posed a problem in response identification as the obtained damping coefficient is much less. Hence, an approximate model is obtained by drawing a smooth curve from the first peak and taking around the midpoint of first valley and the steady state value 1, and reaching the second peak. Thus the responses are fitted as:
Consider the transfer function given by Chien et al (1999):
Gp(s)=(
1
Now y11 (s)=
3.1 Example 1
-11.64 e-0.4 s
s(0.4969s, +0.5057s+1)
The value of p can be obtained from the Laplace inverse of (18) which is given in the Appendix B. Substituting the values of kp11, τe, ζ, y, t, the value of p is obtained as 0.4448.
3. SIMULATION EXAMPLE
4.572s+1
(16)
A better matching of the response is obtained when the y11 is of the form of a SOPTD model with numerator dynamics as
(12)
The main and the interaction responses of the obtained transfer function matrix are compared with the actual process using the same controller settings. A better result can be obtained if these parameters are used as the initial guess and the transfer function model is identified through an optimization method.
22.89 e-0.2 s
kpcij s e-θijs 2 2 (τeij s +2ζijτeijs+1)
To obtain the value of kp21 and kp12, the Laplace inverse of (6) is taken. Appendix B gives the Laplace inverse of (6). Substituting the values of τe, ζ, y, t, in each case of interaction response, the values of kp21 and kp12 are obtained. The obtained responses can be written as:
suitable curve fitting method is applied to solve for kp11, τ11 and ϴ11. The procedure is repeated to Gp21,Gp12 and Gp22 in order to obtain kp21, τ21 and ϴ21, kp12, τ12 and ϴ12 , kp22, τ22 and ϴ22.Thus the transfer function matrix is identified as Gp11 Gp (s) = Gp21
(15)
The analytical method proposed by Yuwana and Seborg is used to identify the model (refer appendix A). For this the first and the second peak values, the first minimum valley and the period of oscillation are noted. From this the effective time constant (τe) and zeta (ζ) are obtained. Closed loop delay is taken same as the open loop delay (ϴ). Process Gain (kp) is found from the final steady state value. In order to identify y22, the above procedure is repeated. The value of kp11 and kp22 are taken as 1. Table 1 shows the observed peak values, the minimum valley, calculated period of oscillation, time constant, zeta and time delay for main and interaction responses. The response can be written in the form of (5). The transfer function of the interaction is assumed as‘s’ times of a SOPTD model.
Substituting Gp,y,yr,I,Gc in the matrix form and solving, the equations for Gp11,Gp21,Gp12 and Gp22 are obtained as below. Gp11 =
kpcii e*θ++ s
(τeii 2 s2 +2ζii τeii s+1)
y12 (s) =
(14)
991
e-0.4s
(20)
e-0.2s
(21)
(0.5087s2 +0.9896s+1) 0.3658
-0.5372 (0.4814s2 +0.3304s+1)
e-0.4s
(22)
2014 ACODS March 13-15, 2014. Kanpur, India
Gc(s) = 7
The fitted response and the actual responses are shown in Fig (2).The above four equations from (19) to (22) are substituted in (8) to (11) and a graph is plotted between each Gp vs s in each case. Fig.3 shows the Gp vs s graph of the identified responses. Each of the Gp is fitted to a FOPTD model, and the values of kp, τ and ϴ are found out in each case. For this 10 values of Gp for 10 values of s are taken in each case and are solved by using fmincon in MatLab. The obtained values of kp, τ and ϴ in each transfer function element are shown in Table 2. Fig. 4 shows the response curve for the identified model and the actual model using the same controller settings. From Fig.4, it can be seen that the responses are satisfactorily matching. The slight mismatching may be due to the approximation of the main and the interaction responses to a second order delay system. Even in SISO systems an accurate matching of the actual and the identified model is not reported. MIMO system identification encounters more problems because of the interactions from all the transfer functions. However, the model fitting can be further improved by using an optimization method such that the model parameters are selected to minimize the sum of the squared errors between the model and the actual process responses.
y11 (s)=
n- Number of the data points in the process and the model responses. The optimisation problem is solved by using the MATLAB routine lsqnonlin with the trust-region reflective algorithm. For the case of limits for optimization, the lower bound of kp (for positive values) is taken as 1/20th of the initial guess and the upper bound as 20 times the initial guess. For negative values of kp the lower bound of kp is fixed as 20 times the initial guess and the upper bound as 1/20th the initial guess value. For the case of time constant the lower limit is kept as 1/20th of the initial guess and upper limits as 5 times the initial guess. The lower bound for time delay is taken as 1/4th of the initial guess value and the upper bound as 4 times the initial guess value. The number of iterations taken is 31and the time taken for convergence is 36.97 sec. The sum of the ISE value obtained is 1.2865 X 10-8. The computational work is done on Intel Core i5 / 3.10 GHz personal computer. The converged parameters are shown in Table 3. Fig. 5 shows the identified and the actual responses.
6.6 e
10.9s+1
0.03 2s
-0.09s
8
(25)
1
e-s
(26)
s(22.7252s, +2.8556s+1) (1+1.0616s)
e-s
(27)
Similarly for main response y22, closed loop response obtained by Y-S method was found to be not matching with the actual response. Hence, by using fractional overshoot method, the obtained values are ζ =0.4331and τe =4.5475. The response can be written as y22 (s)=
1 s(20.6798s2 +3.939 s+1)
y21 (s) =
-18.9 e-3 s -19.4 e
s(26.639s, +1.441s+1)
Now y11 (s) =
Consider the transfer function given by Wood and Berry (1973): ) -3 s
−0.03-
0
0
A better matching of the response is obtained when the y11 is of the form of a SOPTD model with numerator dynamics as (18). p is obtained from the Laplace inverse of (18) given in Appendix B. Substituting the values of kp11, τe, ζ, y, t, the value of p is obtained as 1.0616.
3.2 Example 2
21s+1
+0.6s
A good matching of the identified model using Yuwana and Seborg method was not found with the actual response. It has been reported that Yuwana and Seborg method provides poor parameters for processes with large time delays. This is due to the use of a first order Pade approximation for time delay in the closed loop transfer function denominator (Jutan and Rodriguez, 1984). Hence effective time constant (τe) and zeta (ζ) are calculated by fractional over shoot method (Appendix C). The obtained values are ζ =0.4044 and τe =4.7671.
i= 1, 2; j=1, 2
Gp(s)=( 16.7s+1 -7 s
s
The simulation is done for 150 sec with a sample time of 0.01. The transfer function of the main responses is assumed as a SOPTD model of the form of (15) and the response model is assumed in the form of (5). The analytical method proposed by Yuwana and Seborg is used to identify the model response (refer appendix A). Table 4 shows the observed peak values, the minimum valley, calculated period of oscillation, time constant, zeta and time delay for main and interaction responses. The transfer function of the interaction is assumed as ‘s’ times a SOPTD model of the form of (16) and the interaction response model of the form (6). The value of kp11 and kp22 are taken as 1. To obtain the value of kp21 and kp12, the Laplace inverse of (6) is taken. Appendix B gives the Laplace inverse of (6). Substituting the values of τe, ζ, y, t, in each case of interaction response, the values of kp21 and kp12 are obtained. The obtained responses by using Yuwana and Seborg method can be written as:
Minimize ϕ = ∑526 ∑536 ∑4[0123 (tn) – yij(tn)]2 ∆t (23) (kpmij, τmij, ϴmij) ϕ = objective function
12.8 e- s
0.05
0.2+
(24)
y12 (s)=
14.4s+1
The controller settings are taken from Viswanathan et al (2001)
e-3s
4.4832 (32.0763s2 +2.5973s+1) 3.1235
(30.848s2 +3.6135s+1)
e-7s
e-3s
(28) (29) (30)
The fitted response and the actual responses are shown in Fig. (6).
992
2014 ACODS March 13-15, 2014. Kanpur, India
As mentioned earlier, equations from (27) to (30) are substituted in (8) to (11). A graph is plotted between each Gp vs s in each case. Fig.7 shows the Gp vs s graph of the identified responses. Each of the Gp is fitted to a FOPTD model, and the values of kp, τ and ϴ are found out in each case by using fmincon in MATLAB. The obtained values of kp, τ and ϴ in each transfer function element are shown in Table 5. When these obtained values are taken as initial guess for optimization, a good matching of the identified response with the actual response is obtained. The number of iterations taken is 14 and the time taken for convergence is 68.1011 sec. The sum of the ISE value obtained is 3.4131 X 10-7. The converged parameters are shown in Table 6. Fig. 8 shows the identified and the actual responses. The method can be extended to 3x3 systems or 4 x 4 systems. In these cases, the response models are identified for step changes in each set point one at a time.
Fig.2. Comparison plot of fitted response ( __ ) and actual response(----) (Example 1) (Closed loop reaction curve response)
Table 1. Main and Interaction responses (Example 1): Parameters
y11
y21
y12
y22
yp1 yp2 ym1 T τe ζ ϴ
1.369 1.034 0.891 4.745 0.704 0.358 0.2
0.331 0.0261 -0.124 4.7570 0.7132 0.3354 0.2
-0.5486 0.254 2.244 0.6938 0.2381 0.4
1.1162 0.9475 4.330 0.6681 0.2452 0.4
Table 2. Transfer function parameters obtained from Gp vs s curve (Example 1) Gp11 Gp21 Gp12 Gp22
kp 10.2890 5.1600 -15.1857 7.3514
τ 1.9820 2.6449 2.7003 1.2367
ϴ 0.2635 0.1773 0.3978 0.6447
Fig.3. Plot of Gp vs s curve (Example 1)
Table 3. Converged parameters from optimization using the parameters in Table 2 as initial guess (Example 1)
Gp11 Gp21 Gp12 Gp22
kp 22.8890 4.6888 -11.6401 5.8003
τ 4.5719 2.1749 1.8075 1.7998
ϴ 0.2002 0.2003 0.4001 0.4002
Fig.4. Comparison plot of the closed loop response of identified vs actual transfer function model for the same controller settings (Example 1)
993
2014 ACODS March 13-15, 2014. Kanpur, India
Fig.5. Response plot of the identified and actual transfer model using optimization ( _____ identified response; ------actual response ( both the responses over lap)) (Example 1)
Fig. 6. Comparison plot of fitted response ( __ ) and actual response(----)(Example 2) (Closed loop reaction curve response)
Table 4. Main and Interaction responses (Example 2): Paramete rs yp1 yp2 ym1
y11 1.2493 1.1095 0.8449
y21 0.5687 0.0916 -0.3205
y12 0.2826 0.0297 -0.1003
y22 1.212 1.08 0.9335
T τe ζ ϴ
32.75 5.1613 0.1396 1
36.56 5.6636 0.3393 7
36.9050 5.5541 0.3253 3
31.7 4.8853 0.2497 3
Table 5. Transfer function parameters obtained from Gp vs s curve (Example 2):
Gp11 Gp21 Gp12 Gp22
kp 11.2633 7.2743 -16.7834 -20.4995
τ 24.2644 24.0863 19.8797 12.7444
Fig.7. Plot of Gp vs s curve (Example 2)
ϴ 2.8368 7.4754 2.9102 4.5955
Table 6. Converged parameters from optimization using the parameters in Table 2 as initial guess (Example 2): kp
τ
ϴ
Gp11
12.7565
16.6364
1.0011
Gp21
6. 5759
10.8514
7.0034
Gp12
-18.8176
20.9070
2.9975
Gp22
-19.3505
14.3516
3.0013
Fig.8. Closed loop response plot of the identified and actual transfer matrix using optimization ( _____ identified response; ------- actual response ( both the responses over lap)) (Example 2)
994
2014 ACODS March 13-15, 2014. Kanpur, India
Chien, I. L., Huang, H. P., Yang, J. C.(1999).A Simple Multiloop Tuning method for PID Controllers with No Proportional Kick. Ind. Eng. Chem. Res., 38, 1456-1468. Chidambaram, M. (1998), Applied Process Control, Allied Publishers, New Delhi, India. Davison, E.J. (1976). Multivariable tuning regulators: The feed-forward and robust control of general servomechanism problem. IEEE Trans Auto. Contr, AC-21, 35-41. Ham,T.W. and Kim,Y.H.(1998). Process Identification and PID Controller Tuning in Multivariable Systems. Journal of Chemical Engineering of Japan, 31(6), 941949. Jutan, A., and Rodriguez, E.S. (1984). Extension of a new method of on-line controller tuning. Can. J. Chem. Eng., 62, 802-807. Kavdia, M. and Chidambaram, M. (1996). On-Line Controller tuning for unstable systems. Computers chem.Engng, 20(3), 301-305. Melo,D.L. and Friedly,J.C, (1992).On-Line, Closed-Loop identification of Multivariable Systems. Ind.Eng.Chem.Res., 31, 274-281. Padma Sree,R .and Chidambaram,M. (2002). Identification of unstable transfer model with a zero by optimization method. J.Indian Inst.Sci., 82, 219-225. Pramod,S. and Chidambaram, M. (2000). Closed loop identification of transfer function model for unstable bioreactors for tuning PID controllers. Bioprocess Engineering, 22,185-188 Pramod,S.and Chidambaram,M. (2001).Closed loop identification by Optimization method for tuning PID controllers. Indian Chem.Engr., 43 (2),90-94. Rajakumar,A. and Krishnaswamy, P, R. (1975). Time to Frequency Domain Conversion of Step Response Data. Ind. Eng. Chem., Process Des. Dev., 14 (3), 1975. Rajapandiyan, C.and Chidambaram, M. (2012a). Closed loop Identification of Multivariable systems by Optimization method. Ind. Eng. Chem. Res., 51, 1324- 1336. Rajapandiyan,C. and Chidambaram,M. (2012b). ClosedLoop Identification of Second-Order Plus Time Delay (SOPTD) Model of Multivariable Systems by optimization Method. Ind. Eng. Chem. Res., 51, 96209633. Sundaresan,K.R. and Krishnaswamy,P.R. (1978), Estimation of Time Delay Time Constant Parameters in Time, Frequency, and Laplace Domains. The Canadian Journal of Chemical Engineering, 56, 257-262. Pradeep K. Viswanathan, Wei Khiang Toh, and G. P. Rangaiah (2001), Closed-Loop Identification of TITO Processes Using Time-Domain Curve Fitting and Genetic Algorithms. Ind. Eng. Chem. Res., 40, 28182826. Wood, R. K.; Berry, M. W. (1973), Terminal Composition Control of a Binary Distillation Column. Chem. Eng. Sci., 28, 1707-1717. Yuwana, M. and Seborg, D.E. (1982). A new method for online controller tuning. AIChE J., 28,434-440.
4. CONCLUSIONS The closed loop reaction curve method for SISO systems proposed by Yuwana and Seborg (1982) is extended to identify multivariable systems. The problem associated with this method is brought out. Simulation application is given for transfer function matrix of 2x2 multivariable systems. Two examples are given. The method involves first by fitting the main actions and the interactions of the closed loop system for step changes in the set points by a SOPTD transfer functions with (1/s) term and the interaction by SOPTD model. From the analytical expressions, the open loop transfer function matrices are obtained using the closed loop transfer function matrix and the controller transfer function matrix by plotting Gp vs s curve. The obtained model is used with the same PID controllers and the responses/interactions are matched with that of the actual systems. A method to improve the matching involves by taking the approximate model parameters as the guess values of the least square optimization method. Any arbitrary guess values do not give convergence of the optimization problem.
5. NOMENCLATURE Gp(s) Gc(s) gp, gc, gm kpm, kp, kpc τm, τ θm, θ yr,ymr Y y,ym yp1 yp2 ym1 T τe ζ t kc τI τD ∆t FOPTD SOPTD ISE
process transfer function matrix controller transfer function matrix process, controller and model transfer functions model and process steady state gains closed loop process gain model and process time constant model and process time delay process and model set point closed loop response matrix closed loop process and model response first peak in the closed loop response (refer Fig.9) second peak in the closed loop response (refer Fig.9) first valley in the closed loop response(refer Fig.9) period of oscillation (refer Fig.9) effective time constant damping coefficient time controller gain Integral time constant Derivative time constant sampling time first order plus time delay system second order plus time delay system Integral of squared error REFERENCES
Ananth, I. and Chidambaram, M. (1999). Closed-loop identification of transfer function model for unstable systems. Journal of the Franklin Institute, 336, 10551061. 995
2014 ACODS March 13-15, 2014. Kanpur, India
Appendix A. FIRST APPENDIX
Fig. 9 shows the closed loop response of an under damped process The expression for the effective time constant (τ) and the damping coefficient (ζ) can be derived as: 9∞
9:;
9<;
9∞
9<;
v 1=
9∞
9<,
v 2=
ζ =
ζ5 =
9∞
(
KL =
(A-2)
=4 (>; )
[?, @A=4(>; )B, ]D.F
(A-3)
=4 (>, )
[G?, @A=4(>, )B, ]D.F
ζ = ζ1 + ζ2 T=
(A-1)
5?IJ
M(
,
(A-4) (A-5)
ζ )D.F
(A-6)
,
ζ )D.F
5?
(A-7)
Appendix B. SECOND APPENDIX y(s)=
(1+ps)e-θs
s(τe , s2 @2ζτe s+1
(B-1)
Laplace inverse of (B-1) is
y(t)=A [sinbexp( )]B+{1- expN O [qζ sin(b) + cos(b)]B pq
-ζt’
-ζt’
τe
τe
τe
(B-2)
y(s)=
kp e-θs 2 τe s +2ζτe s+1
(B-3)
Laplace inverse of (B-3) is y(t)=
kp q τe
[exp( )sinb] -ζt’ τe
where t’=t-ϴ ; b=(1-ζ2)0.5
(B-4) RS
IJ
; q=(1-ζ2) -0.5
Appendix C. THIRD APPENDIX Max Overshoot final value
τeff =
=e
TT1-ζ2 2π
πζ
T1-ζ2
(C-1) (C-2) 996