1. Introduction
Linear optimal control theory has been extensively studied since Bellman’s dynamic programming [1] and Pontryagin et al.’s minimum principle [2] were published early in the 1960s. In linear optimal control, the linear quadratic regulator (LQR) problem is to find the optimal feedback control law to minimize the quadratic performance criterion with respect to the states and inputs via weighting matrices Q and R [3-5]. Thus the central issue remains of how to relate the weighting matrices in the quadratic performance criterion to classical specifications in the time and frequency domains for the control system design [6-9].
Various analytic methods have been developed for frequency domain specifications such as disturbance attenuation, noise rejection, and good command following for the LQR [10, 11] and linear quadratic Gaussian / linear transfer recovery (LQG/LTR) [12, 13], based on the loopshaping method. These existing methods have difficulty directly matching the quadratic performance indices with the filter dynamics to the frequency loop shape [3, 4].
On the other hand, eigenvalue positioning methods were proposed for the time domain approach to select the appropriate weighting matrices, Q and R, to locate the closed loop poles inside the desired region or match them with the specific eigenvalues through the ARE [4]. Existing pole positioning approaches are limited in their analysis of time domain performances such as overshoot, rising and settling time [4]. LQ-proportional integral (PID) control methods have been also developed separately in time and frequency domains [14, 15].
We propose a new optimal control method for LQ-PID design to simultaneously satisfy the design specifications in the time- and frequency-domain without any filter or pole-zero cancelation [16]. For frequency-domain performances, we improved the existing limiting behaviour method of the loop transfer function to more exactly interpret the loop shape with respect to all design factors, based on Kalman’s equality [17, 18]. This loop shaping makes it possible to consider the magnitudes of the LQ loop system independently in low and high frequency ranges for good command following, disturbance rejection, and attenuation of noise effect [19]. Time domain specifications such as percent overshoot, rise time, and settling time are related to the design factor by target function matching regardless of the pole-zero cancellation [16]. Overshoot and the Routh- Hurwitz stability criterion are formulated by the coefficient analysis with respect to the normalized third-order closed-loop system [20]. The other time-related specifications are addressed by the timeregulation property of the target function [21]. In the proposed method, the LQ-PID control parameters are determined by optimizing to minimize the cost function, which is subject to constraints to meet the combined time - and frequency-domain specifications for classical performances. We choose a combination of integral absolute error (IAE) and 2-norm of the control gains for the cost function of the optimal problem [22].
The primary contributions of the proposed method include as follows:
(i) The proposed method provides the direct matching manner to relate the weighting factor, Q, to the loop shape of the loop transfer function in frequency domain. (ii) The proposed target function matching method improves the existing time domain approach of LQPID control by considering the specific time domain performances of transient response such as overshoot, rising time, settling time. (iii) Both of time and frequency domain specifications are combined by the proposed optimization problem in order to be simultaneously satisfied.
The paper is organized as follows. In Section 2, we introduce the LQ-PID formulation to design the optimal PID controller by transforming the PID control into the LQ approach. The design methodologies for the proposed controller are discussed as the optimization problem with constraints subject to the time- and frequency-domain specifications in Section 3. Numerical examples are given in Section 4 to show the effectiveness of the proposed design method, compared to other design methods. Finally, Section 6 contains some concluding remarks and remaining challenges.
2. LQ-PID Control
Briefly, we introduce the LQ-PID control method proposed by Suh and Yang [14], in which the optimal linear feedback control law is matched to the conventional PID control for the second order system by augmenting the integral of the output variable as a new state.
Consider the following second order model:
where y(t) , u(t) , ζ and ωn are the output variable, the control variable, the damping ratio and the natural frequency. The initial conditions of y(t) and dy(t) / dt are specified.
The new augmented state variables are:
The state-space representation of the augmented system becomes:
where
If we apply the different input gain according to and as
where Ac is the closed loop system matrix.
The time delay system can be represented as the system without the time delay as where [23].
The quadratic performance criterion is considered for the LQR formulation:
with the assumption that the weighting factor, Q , is positive semi-definite and symmetric Q = QT ≥ 0 , and ρ is a positive value.
The linear feedback control law to regulate the state, x(t) , is obtained as follows:
where the optimal gain matrix, K = KT , is obtained by solving the ARE:
Let the components of K and G be:
which we will hereafter refer to as the optimal Kalman gains of the LQ-PID.
The optimal control law of Eq. (9) is represented as the following PID control formula:
The state feedback system with the control input, u(t) , of Eq. (8) is obtained by substituting Eq. (12) into Eq. (3):
So, the closed loop system matrix, Ac , is:
and its characteristic equation becomes:
The closed loop transfer function of the second order system with the conventional PID controller can be expressed as follows:
where and the transfer function of Eq. (1) is
The characteristic equation of Eq. (19) is
Fig. 1 shows the block diagram of LQR of the state feedback system with the augmented state variable, forming PID structure where the optimal tuning of the conventional PID controller is achieved by LQR.
Fig. 1.Structure of LQ-PID control
By corresponding Eq. (18) to Eq. (16), the parameters of the PID controller are matched with the optimal Kalman gains of the LQ-PID structure from Eq. (13).
It has the following relationship
such that the PID controller is tuned by selecting the weighting matrices, Q and ρ , because the parameters of the PID controller are related to the solution, K , of the ARE through Eqs. (19) and (12), and ωn is given by the system plant.
Let the weighting matrix, Q , be as Eq. (20):
where N = [n0 n1 n2] is assumed as a partitioned matrix of the weighting matrix, Q .
3. Design Methodology
In this section we describe the proposed design procedure of the LQ-PID controller to meet both the frequency and time domain specifications.
Given a specific set of design requirements in the time domain such as overshoot, and rise and settling times, the feasible area of design parameters will be determined by pole matching to the normalized simple third order system [15, 18]. The parameters are then connected to weighting factors, Q and ρ , by the ARE. For the frequency specifications, we consider the limiting behaviour for the loop shaping method in which Kalman’s inequality is applied [16, 21]. The time- and frequency-domain approaches are combined by formulating the optimal problem with constraints in terms of Q and ρ to satisfy both time- and frequency-domain specifications. We consider the cost function of the proposed optimal problem as the weighted summation with the 2-norm of control gains and the integral of the absolute error in order to avoid the saturation of input energy and reduce the overall effect of the error.
3.1 Combined constrains for time and frequency domain
3.1.1 Frequency domain performances
In the frequency domain, there are many performance indices such as gain margin, phase margin, infinite norm of the sensitivity function, and the complementary sensitivity function. The LQ-PID control guarantees the robust stability of the optimal system, which is verified by Kalman’s inequality [15]. The inequality verifies that the sensitivity performance is improved at all frequencies,
|S( jω)| ≤ 1 , for all ω [3, 4, 17, 19].
The sensitivity function is defined in the usual way,
where the loop transfer function gLQ(s) of LQR is obtained as follows:
In order to associate the frequency loop shaping procedure with the weighting factors, Q and ρ , we use Kalman’s identity for the single input single output case:
where gOL(jω) = N( jωI − A)−1B with a partitioned matrix, N, of the weighting matrix, Q,. Eq. (23) can be represented as:
Deriving the magnitude of the complex function in Eq. (24), we have:
Eq. (25) can be approximated as Eq. (26) only for low and high frequencies since |gLQ( jω)| >> 1 and |gLQ( jω)| << 1 , respectively [15].
The magnitudes of the loop transfer function, gLQ(s) , are approximated via the limiting behaviours of gLQ(jω) at low and high frequencies [19], respectively, as:
and
Let each asymptotic line of |gLQ( jω)| at the low and high frequency be lL(ω) and lH(ω) , respectively.
We consider a design procedure that involves explicitly shaping the magnitude of the loop transfer function, gLQ(s) , for frequency domain performances.
In classical loop shaping, better performances such as good command following and disturbance rejection are obtained when the magnitude of the loop transfer function, gLQ(s) , is larger for low frequencies at |gLQ( jω)| >> 1 [19]. Fortunately, robust stability, noise attenuation, and control energy reduction are valid primarily in high frequencies at |gLQ( jω)| << 1 . The shape of |gLQ( jω)| then should not invade either the low and high frequency barriers, Q and β (ω) , respectively, as described in Fig. 2.
Fig. 2.Typical shape of |gLQ( jω)| with barriers for the requirements
For this adequate loop shaping, we can develop a loop shaping technique in which the design parameter, N , plays a key role in maintaining |gLQ(jω)| beyond the required boundaries. Specifically, this means that we match the weighting matrix,Q , of the LQ problem directly with the performances of the frequency domain to satisfy the boundary conditions via the approximation of limiting behavior.
Let Ωr be a boundary frequency of the command following barrier,α (ω) , and scalar r m be a magnitude of α (ω) , as shown in Fig. 2. Set the barrier,α (ω) , as:
Since lL(ω) should exist on the upper side of the barrier,α (ω) , the asymptotic line of |gLQ( jω)| at the low frequency must satisfy:
Solving Eq. (30) for the range of the parameters, n0 and n1 , gives:
Let Ωr be a boundary frequency of the senor noise, and be the inverse of the maximum value of the modeling error in Fig. 2.
The barrier,β (ω) , can then be expressed as:
The inequality constraint of lH(ω) for high frequencies is derived from Eq. (28) in the same manner:
Hence this constraint can be rearranged with the respect to the element n2 of the partial matrix, N of Q as
Finally, we can yield the constraints of Eqs. (31) and (34) with respect to the weighting factors,Q and ρ , of the LQ problem to satisfy frequency domain performances, adjusting the loop shaping of |gLQ( jω)| .
3.1.2 Time domain performances
In the time domain, the desired transient performances are typically characterized by overshoot, and rise, and settling times. Designs for controllers that improve these transient performances of a second order closed loop system in the time domain are well explained in current literatures; however higher order systems are generally analyzed via the dominant pole concept as the reduced second order system is hard to formulate [24].
In this paper, we consider the special formulation of the third order closed loop system as the target function which is referred to [14] as:
The shape of the unit step response of Eq. (35) is determined by parameters p, r, and ωt. Furthermore, it is found that the overshoot is maintained constantly, regardless of the variations in natural frequency, ωt, when parameters p and r are fixed. This is verified from the formulation of the output, y(t) [25].
where
Since the time, t , and the natural frequency, ωt , are united, we can modify the time, t , to a new time, τ , using ωt as the time ratio, as follows:
where τ = * tωt
The output y(τ ) is represented as the response of the normalized third-order closed-loop system (NTCS), described as Eq. (35) with ωt = 1 . Hence the time ratio of the same step response y(τ) can be regulated independently on variations of the natural frequency, ωt as shown in Fig. 3.
Fig. 3.Step response of the third order system according to variations in ωt
We consider two steps to formulate the constraints of time domain specifications. First, we determine the constraints with respect to p , r , and ωt to satisfy the performance requirement in the time domain. For the overshoot requirement, the constraint function is illustrated as follows:
where fOS(r, p) represents the maximum step response function consistent with the variables, p and r .
Fig. 4 shows that the contours of maximum percentage overshoot with respect to the step reference can be plotted in the coefficient plane of p and r . If the Routh-Hurwitz stability criterion is applied to the normalized system transfer function of Eq. (35), Eq. (40) is required for simple absolute stability with h = 1 :
Fig. 4.Coefficient plane for overshoot in NTCS
A gain margin of 0.33, which is allowable in practice, corresponds to the line of Eq. (40) for h = 1.5 [18]. The region between the two contours at h = 1 and h = 1.5 is the unallowable design area, which is distinct by oblique in the coefficient plane of Fig. 4.
Hence, except in terms of area, the overshoot coefficient is settled according to the boundary of the required overshoot percentage. The natural frequency, ωt, can be addressed as the design parameter, which plays a role in regulating the rise and settling time under the overshoot determined by p and r .
where ωST and ωRT depend on the values of the natural frequency, ωt , to satisfy the required settling and rising time, respectively. In other words, the system parameters p and r are determined to meet the given overshoot requirement, and ωt is considered as a factor used to separately satisfy the given settling time and rising time requirements.
Secondly, the system parameters p , r , and ωt should be formulated as a function with respect to the optimal Kalman gains of the LQ-PID feedback system via pole matching to equalize the characteristic equation of the target function Eq. (16) of NTCS.
Given the above Eq. (42), the following three equations are derived.
Using Eq. (43), weighting matrices,Q , and ρ , can be expressed in terms of the variables of the target transfer function via the ARE to satisfy the time domain performances, as follows:
By substituting Eq. (43) to Eq. (44), the components [n0 n1 n2 ] of weighting matrix, Q are expressed as Eq. (45) in terms with the variables, p , r , and ωt for time domain performances, given input weighting, ρ , and system parameters, ωn , ζ and c of the plant
where f1and f2 are the functions to meet time domain specifications in terms of the system variables. Eq. (39) and (41) then can be represented by the function of the optimal Kalman gains through Eqs. (42), (43), and (44) as follows:
Therefore we can determine the feasible constraints Eqs. (46) and (47) in terms of positive weighting matrices, Q and ρ , of the LQ problem to meet transient step response requirements such as overshoot, rise, and settling time in time domain.
3.2 Formulation of an optimization problem
In the above sections, we derived that the constraint conditions Eqs. (31), (34), (46), and (47) to satisfy the combined time- and frequency-domain requirements. The conditions are represented as inequality functions with respect to weighting matrices,Q and ρ , of LQ problem, obtained via the loop shaping method and the pole matching method by the ARE. In this subsection, we formulate the optimal problem to determine optimal values of the design parameter N to minimize the object function subject to the constraints combined by time- and frequencydomain performances. We choose the weighted combination of two objectives, the 2-norm of PID gain and IAE. This objective of the optimal problem can vary depending on the control purpose of the designer in practice.
The first objective term is represented as Eq. (48) in order to avoid a high gain problem when we take ρ << 1for the cheap control as shown by Eq. (27) [26].
The second objective therm is the IAE index as Eq. (50) for the time response to consider the error arising from the difference between a nominal plant and a closed loop target function Eq. (35):
These two objectives are combined by weighting vector, γ , as follows:
where the weighting vector can be chosen according to the importance of each objective.
In order to combine the time- and frequency-domain performances, the optimization problem is formulated with respect to the object functions in Eq. (50) and the constraints in Eqs. (31), (34), (46), and (47), as follows:
minimize (n0, n1, n2) subject to
fos(ρ,n0,n1,n2)≤Γ, ∀τ fWT(ρ,n0,n1,n2)≥max(ωST,ωRT), fixed p, r n0,n1,n2>0 for positive definite
where ωn andζ are given by the process, and Ωr , mr , Ωn, , Γ , ωST , ωRT , and ρ are given as the design specifications.
The design procedure is summarized as
Summary of Design Procedure
Step 1: Set parameter values of a given plant : (ζ , ωn ) .
Step 2: Specify the frequency domain constraints with respect to the required design specifications (mr , Ωr , Ωn , via the loop shaping method.
Step 3: Specify the time domain constraints with respect to the required design specifications (Γ, ωST, ωRT) via target function matching.
Step 4: Perform the optimization problem of Eq. (51) with the specific value of ρ to determine the weighting factor, Q = NT N , subject to the constraints obtained by Step 2 and 3.
Step 5: Determine the optimal PID control gain, G , by solving the ARE with the system parameters and the weighting matrix Q and R = ρ I determined in Step 4.
Step 6: Finish the process to determine the adequate controller if all of required design specifications are satisfied, if not reset the constraints, beginning with Step 2.
4. Simulation Examples
Example 1: Consider the second order system with a small delay, referred in [27] as
For the optimization, the constraints and weighting vector γ are given as follows
fOS (ρ, n0, n1, n2) ≤ 9% fWT (ρ, n0, n1, n2) ≥ωST (8 sec) γ = [1 0.5], ρ = 0.0001
The frequency and time domain constraints are chosen to avoid the loop shape from the frequency barriers Q and β (ω) , and to satisfy that the overshoot and settling time of the target system are less than 9% and 8 sec, respectively.
According to the design procedures, design parameter N is determined as N = [0.0099 0.0151 0.005] , and PID controller parameters are obtained as Kp = 1.5027, Ti = 0.9967 and 0.5005 Td = via the ARE.
In Table 1, it is shown that the proposed method only yields the adequate loop shape required to avoid invading the frequency barrier, a good step response to meet the time response requirements of less than 9% overshoot and 8 sec settling time, as well as simultaneously minimizing the cost function. Figs. 5 and 6 present the comparisons of the simulation results for frequency and time responses, respectively.
Table 1.System performance comparisons of time and frequency responses for example 1.
Fig. 5.Comparison of the frequency shape of the loop transfer function with various other methods
Fig. 6.Comparison of step response with various other methods for example 1
The results of the simulation are compared with other PID control methods by Yang [27], Ho [28], Grassi [29], and Olof [30]. As shown in the Table 1 and Fig. 5, the responses of the proposed and Olof methods satisfy the frequency-domain specifications. Among them, only the proposed method meets the required time-domain specifications described as settling time (8 sec) and overshoot (9%) in Fig. 6.
Example 2: A fired heater problem is taken from [31] as an example of practical process; this is one of the most common examples with long time delay for the practical process. The open loop dynamic of the system is described by the transfer function, Gl(s) , at low fuel gas feed rates:
For the optimization, the constraints and weighting vector γ are given as follows.
fOS (ρ, n0, n1, n2) ≤ 10% fWT (ρ, n0, n1, n2) ≥ωST (100 sec) γ = [1 0.5], ρ = 0.0001
The design parameter N and PID controller are obtained as
N = [0.196 3.08 5.33]×10-3 Kp = , 0.221 Ti = 0.02, and Td = 0.4567 respectively.
Table 2 shows the results of the proposed controller and those by Grassi [29], Olof [30], Ho [28], and Huang [32]. With the exception of the Ho method, they satisfy the time domain specifications with respect to the percent overshoot (10%) and settling time (100 sec). On the other hand, the loop shapes at the low and high frequency are obtained only by the proposed method, and Ho’s. Only by the proposed method presents the allowable good step response for the time domain performance simultaneously satisfied in addition to the frequency response to guarantee a magnitude higher than 15 db at the low frequency and lower than -40 db at high frequency. Figs. 7 and 8 present a comparison of loop shapes in the frequency domain and step responses in the time domain, respectively.
Table 2.System performance comparisons of time and frequency responses for example 2.
Fig. 7.Comparison of the frequency shape of the loop transfer function with other methods for example 2
Fig. 8.Comparison of the step response with other methods for example 2
5. Conclusion
This paper proposes a new tuning method to combine time- and frequency-domain requirements for an LQ-PID controller. The weighting matrices, Q and R, are formulated to improve limiting behaviour of the loop transfer function by Kalman’s equality and target function matching is used for frequency and time domain performances, respectively. Both of them are simultaneously considered as the constraints of the optimal problem with the cost function in which the 2-norm of optimal control law and IAE index are applied. The PID control coefficients are obtained directly through the ARE, including the weighting factors Q and R. Simulation results have been presented to show the effectiveness of the proposed method in comparison to others. In conclusion, the proposed method guarantees not only the robust stability and design simplicity inherent in LQ-PID optimal control, but also considers the time- and frequency-domain performances for the second order system.
참고문헌
- R. Bellman and S. Dreyfus, “Applied dynamic programming”, NJ: Princeton University Press, 1962.
- L. Pontryagin, V. Boltyanskii, R. Gamkrelidze, and E. Mishchenko, “The mathematical theory of optimal processes”, Wiley Press, 1962.
- B. Anderson and J. Moore, “Optimal control: Linear quadratic methods”, Prentice Hall Press, 1989.
- P. Dorato and V. Cerone, “Linear-quadratic control: An introduction”, Prentice Hall Press, 1995.
- Y. Shih and C. Chen, “On the weighting factors of the quadratic criterion in optimal control,” Int. J. Control, vol. 19, pp. 947-955, 1974. https://doi.org/10.1080/00207177408932688
- M. Nekoui and H. Bozorgi, “Weighting matrix selection method for lqr design based on a multiobjective evolutionary algorithm,” Advanced Materials Research, vol. 383, pp. 1047-1054, 2012.
- S. Ghoreishi and M. Nekoui, “Optimal weighting matrices design for LQR controller based on genetic algorithm and PSO,” Advanced Materials Research, vol. 433, pp. 7546-7553, 2012.
- E. V. Kumar, J. Jerome, and K. Srikanth, “Algebraic approach for selecting the weighting matrices of linear quadratic regulator,” Green Computing Communication and Electrical Engineering (ICGCCEE), 2014, pp. 1-6.
- S. Das, I. Pan, K. Halder, S. Das, and A. Gupta, “LQR based improved discrete PID controller design via optimum selection of weighting matrices using fractional order integral performance index,” Applied Mathematical Modelling, vol. 37, pp. 4253-4268, 2013. https://doi.org/10.1016/j.apm.2012.09.022
- J. Moore and D. Mingori, “Robust Frequency-shaped LQ Control,” Automatica, vol. 23, pp. 641-646, 1987. https://doi.org/10.1016/0005-1098(87)90060-4
- Y. Shi, T. C. Becker, S. Furukawa, E. Sato, and M. Nakashima, “LQR control with frequency-dependent scheduled gain for a semi-active floor isolation system,” Earthquake Engineering & Structural Dynamics, vol. 43, pp. 1265-1284, 2014. https://doi.org/10.1002/eqe.2352
- D. L. Silva, C. F. Paula, and L. H. Ferreira, “On the target feedback loop parameterization for H∞/LTR control,” In Industrial Electronics (ISIE), 2013 IEEE International Symposium on, 2013, pp. 1-5.
- G. Stein and M. Athans, “The LQG / LTR procedure for multivariable feedback control design,” IEEE Trans. Automat. Contr., vol. 32, pp. 105-114, 1987. https://doi.org/10.1109/TAC.1987.1104550
- B. Suh and J. Yang, “A Tuning of PID Regulators via LQR Approach,” J. Chem. Eng. Japan, vol. 38, pp. 344-356, 2005. https://doi.org/10.1252/jcej.38.344
- J. Yang and B. Suh “A new loop-shaping procedure for the tuning LQ-PID regulator,” J. Chem. Eng. Japan, vol. 40, pp. 575-589, 2007. https://doi.org/10.1252/jcej.40.575
- M. Thomson, P. Cassidy, and D. Sandoz, “Automatic tuning of PID controllers using a combined time-and frequency-domain method,” Trans Inst. MC, vol. 11, pp. 40-47, 1989. https://doi.org/10.1177/014233128901100105
- J. Burl, “Linear optimal control: H2 and H∞ methods,” Addison Wesley Press, 1998.
- R. Kalman, “When is a linear control system optimal?,” J. Basic Eng., vol. 86, pp. 51-60, 1964. https://doi.org/10.1115/1.3653115
- M. Athans, “Lecture note on multivariable control system,” M. I. T., Lecture Note, Ref. No. 860224/6234, 1986.
- D. Towill, “Analysis and synthesis of feedback compensatedthird-order control systems via the coefficient plane,” The Radio and Electronic Eng., vol. 32, pp. 119-131, 1966. https://doi.org/10.1049/ree.1966.0064
- A. Marques, Y. Geerts, M. Steyaert, and W. Sansen, “Settling time analysis of third order systems,” Electronics, Circuits and Systems, IEEE International Conference, 1998, pp. 505-508.
- R. C. Dorf and R. H. Bishop, “Modern control systems,” 12th ed., Prentice Hall, 2010.
- J. He, Q. Wang, and T. Lee, “PI/PID controller tuning via LQR approach,” Chem. Eng. Sci., vol. 55, pp. 2429-2439, 2000. https://doi.org/10.1016/S0009-2509(99)00512-6
- N. Nise, “Control systems engineering,” John Wiley & Sons Press, 5th ed., 2008.
- K. Ogata, “Modern control engineering,” Prentice Hall Press, 4th ed., 2002.
- H. Kwakernaak and R. Sivan, “Linear optimal control systems,” Wiley Interscience Press, 1967.
- Y. Yang, Q. Wang, and W. Dai, “Robust PID controller design for gain and phase margins,” J. Chem. Eng. Japan, vol. 35, pp. 874-879, 2004.
- W. Ho, C. Hang, and L. Cao, “Tuning of PID controllers based on gain and phase margin specifications’, Automatica, vol. 31, pp. 497-502, 1995. https://doi.org/10.1016/0005-1098(94)00130-B
- E. Grassi and K. Tsakalis, “PID controller tuning by frequency loop-shaping’ Proc. IEEE CDC, Kobe, Japan, Dec. 1996, pp. 4776-4781.
- G. Olof, “Design of robust pid controllers with constrained control signal activity,” Licentiate Thesis, Lund University, 2009.
- Y. Wang and J. Rawlings, “A new robust model predictive control method. II: examples,” J. Proc. Cont., vol. 14, pp. 249-262, 2004. https://doi.org/10.1016/S0959-1524(02)00132-4
- H. Huang, M. Lee, and C. Chen, “Inverse-based design for a modified PID controller,” J. Chin. Inst. Chem. Eng., vol. 31, pp. 225-236, 2000.