A New Approach to Mixed - American Chemical Society

Ind. Eng. Chem. Res. 2002, 41, 6107-6119

6107

A New Approach to Mixed H2/H∞ Optimal PI/PID Controller Design Chyi Hwang* and Chun-Yen Hsiao Department of Chemical Engineering, National Chung Cheng University, Chia-Yi 621, Taiwan

This paper proposes a new approach to the problem of designing optimal proportional-integralderivative (PID) controllers that minimize an H2 norm associated with the set-point response subjected to a constraint on the H∞ norm of disturbance rejection. The proposed design approach consists of constructing the feasible domain in the controller gain space and searching over the domain for the optimal gain values of minimizing the H2-norm objective function. The construction of the feasible domain that satisfies the requirements of closed-loop stability and H∞-norm constraint is achieved through analytic characterization of the domain boundary with the notion of principal points associated with the value set of a differentiable mapping. The constructed feasible domain saves greatly on the computational effort requied to search for the H2-optimal controller gain values that satisfy both the stability and disturbance rejection requirements. 1. Introduction The proportional-integral-derivative (PID) control algorithm has played a dominant role in industrial process control systems for over 60 years.1 Consequently, many different methods have been proposed for determining the three controller parameters, i.e., proportional, integral, and derivative gains, to meet different requirements of various control applications. The proposed PID controller design methods can be broadly grouped into three categories. The first category consists of methods that make use of frequency-response information. The early Ziegler-Nichols (Z-N) tuning method2 is based on the knowledge of only one point (critical point) on the process Nyquist curve. Recently developed Z-N-like methods3-11 try to automate and improve the PID controller tuning using a slight increased amount of process frequency response data. The second category includes model-based controller design methods.12-14 These methods require a process transfer function model, and they derive PID controller settings with goals such as pole placement and internal model control. The third category consists of design methods15,16 that determine the PID controller parameters by minimizing some integral of feedback error performance criteria. This class of methods can be applied to general transfer function models and typically requires a numerical optimization method to search for the optimal PID controller parameters. Because the PID control algorithm has widespread use in industrial reality, there is a significant incentive to apply modern insights to the improvement of PID controller design methods. Several authors have recently attempted to apply modern H∞ control theory to the design of PID controllers. Two approaches have been adopted to achieve this goal. One is to first design a controller using H∞ control theory and then try to reduce it to a PID structure. Grimble17 showed that H∞observation-weighted control might result in controllers * To whom all correspondence should be addressed. Phone: 886-5-272-1208. Fax: 886-5-272-1206. E-mail: chmch@ ccu.edu.tw.

with a PID controller structure. Also, the use of loop shaping suggested by McFarlane and Glover18,19 might give rise to controllers with a PID controller structure.20-22 The other approach is the direct search for optimal PID controller parameters that minimize some performance criterion such as the H2 norm of the feedback error to step setpoint input while satisfying H∞ constraints on robustness and/or disturbance rejection.23-25 In particular, A° stro¨m et al.26 considered the design of PI controllers by minimizing the integral of error due to a step disturbance input subject to constraints on maximum sensitivity and/or complementary sensitivity. Recently, multivariable and multiloop PID controllers have been designed in the H∞ control framework (see Bao et al.27 and Zheng28 and references therein). In these cited works, the LMI approach was used to solve the design problem. The purpose of this paper is to present a new approach to the design of optimal PID controllers that minimizes an H2-norm performance index subject to an H∞-norm constraint on external disturbance attenuation. This mixed H2/H∞ optimal PID controller design problem has been previously studied by Chen et al.23 and Krohling.24 Chen et al.23 proposed a hybrid solution approach consisting of a genetic algorithm with binary codification for minimization of the integral of squared error (ISE) performance index, together with a numerical algorithm for evaluating a disturbance rejection constraint. Their method requires the evaluation of the H∞ norm and stability for each candidate solution and is only applicable to delay-free systems. Instead of using ISE as a nominal H2-norm performance index, Krohling24 solved with two real-coded genetic algorithms the mixed H2/H∞ optimal PID control design problem for the minimization of the integral of time-weighted squared error (ITSE) in conjunction with an H∞ disturbance attenuation constraint. Our improvements on this previous attempt lie in three aspects. First, we consider the general case where the process can be described by a rational transfer function with a pure delay or by an irrational transfer function. Hence, our result is

10.1021/ie010963m CCC: $22.00 © 2002 American Chemical Society Published on Web 10/30/2002

6108

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002

ment for the controller gain settings is that the feedback control system in Figure 1 be internally stable, i.e., the following three transfer functions

E(s) 1 ) R(s) 1 + Gc(s) Gp(s) Figure 1. Block diagram of the PID feedback control system.

applicable to a larger class of plants, especially to most chemical processes in which a transportation delay or a distributed delay always exists. Second, we employ a generalized quadratic error functional associated with the step response of setpoint input as the performance index to be minimized. The generalized quadratic cost functional includes the commonly used ISE and ITSE performance criteria as special cases. Last and most importantly, we propose in this paper to identify explicitly the feasible domain in the controller gain space with the notion of principal points29-31 of characterizing the value set boundary of a differentiable mapping. With the controller parameters selected from the feasible domain, the constraints on the internal stability and the disturbance attenuation H∞ norm are automatically guaranteed. The explicit construction of a feasible gain domain can greatly reduce the computational effort required to search for the optimal PID controller gain values. The remainder of the paper is organized as follows: Section 2 describes the mixed H2/H∞ optimal PI/PID control design problem. In section 3, the procedure for constructing a feasible gain domain is developed. It is based on the analytical characterization of the value set boundary via the notion of the principal points associated with a differentiable mapping. In section 4, parametric expressions are derived for computation of the H2-norm performance indices for plants with a pure time delay. In section 5, the optimization method for searching for the minimum H2-norm PID controller parameters from the feasible gain domain is described. Two design examples using the proposed approach are illustrated in section 6. Finally, conclusions are given in section 7. 2. Problem Description Consider the block diagram of the feedback control system shown in Figure 1. In this block diagram, Gc(s) and Gp(s) are transfer functions of the controller and the plant to be controlled, respectively, and R(s), D(s), and Y(s) denote the Laplace transforms of the setpoint input variable r(t), disturbance input variable d(t), and process output variable y(t), respectively. For PID feedback control, the controller transfer function is given by

Gc(s) ) kp +

ki + kds s

(1)

where kp, ki, and kd are the gains of the proportional, integral, and derivative control actions, respectively. The tuning of PID controller Gc(s) is to determine the values of its three gains. Usually, the primary require-

Gp(s) Y(s) ) ≡ Gd(s) D(s) 1 + Gc(s) Gp(s)

(2)

Gc(s) U(s) ) R(s) 1 + Gc(s) Gp(s) must be bounded-input-bounded-output (BIBO) stable. Note that a transfer function G(s) is BIBO stable if and only if it is proper and analytic on the closed right-half plane of the s plane. Hence, the three parameters of the PID controller Gc(s) must satisfy the following stability constraint

1 + Gc(σ+jω) Gp(σ+jω) * 0

for σ g 0

(3)

where j ) x-1. Let L2 denote the space of signals with finite energy and define the L2 norm as

||f||2 )

x∫ | f(t)| ∞

2

0

dt

(4)

for f(t) ∈ L2. Let yd(t) denote the output response due to the presence of external disturbance d(t) ∈ L2. Then, under the assumption that the feedback control system is internally stable, the desired disturbance attenuation level γ < 1 is specified by

Gp(s) ||yd||2 ) || || e γ d(t)∈L2 ||d||2 1 + Gc(s) Gp(s) ∞ sup

(5)

where ||G(s)||∞ denotes the H∞ norm of a BIBO-stable transfer function and is defined as

||G(s)||∞ ) sup |G(jω)|

(6)

ω∈[0,∞)

Hence, to achieve the specified disturbance attenuation level γ, the gain parameters of the PID controller Gc(s) must be chosen to stabilize the closed-loop system and satisfy the following constraint on the H∞ norm of the disturbance transfer function Gd(s)

||Gd(s)||∞ ) sup

|Gp(jω)|

ω∈[0,∞)|1

+ Gc(jω) Gp(jω)|

eγ

(7)

Now, we consider the set-point tracking ability of the feedback control system in Figure 1. From eq 2, we know that the error e(t) due to a unit step setpoint input R(s) ) 1/s has its Laplace transform as follows

E(s) ) L {e(t)} )

R(s) 1 + Gc(s) Gp(s)

)

1 s+(kds +kps+ki)Gp(s) 2

(8)

Ind. Eng. Chem. Res., Vol. 41, No. 24, 2002 6109

where eq 1 has been used. An analytical measure of the tracking ability for the feedback control system in Figure 1 can be defined as follows

∫0∞tme2(t) dt

Jm )

(9)

where m is a nonnegative integer. It is noted that, if the feedback control system satisfies the internal stability requirement (eq 3), the Laplace transform of the time-weighted error variable tm/2e(t), denoted by Em/2(s), belongs to the H2 Hardy space consisting of all functions G(s) that are analytical in the closed right-half plane of the s plane and have a finite H2 norm, defined as

||G||2 )

∫-∞∞G(jω)G(jω) dω

1 j2π

(10)

where G h denotes the complex conjugate of G and G(jω) ) G(-jω) if G(s) is a function with real coefficients. Because Em/2(s) is a real-coefficient function of s and belongs to H2, the use of Parseval’s identity32 leads to the equality

Jm )

∫0∞tme2(t) dt ) 1 ∞ ∫ E (jω) Em/2(-jω) dω ) ||Em/2||2 j2π -∞ m/2

k

(11)

(12)

subject to the constraints in eqs 3 and 5 where k ) (kp, ki, kd). To solve this constrained optimization problem, we suggest first to construct the feasible gain domain

Kf ) {k: kp g 0, ki g 0, kd g 0, 1 + Gc(σ+jω) Gp(σ+jω) * 0, ∀ω ∈ [0,∞) and σ g 0; ||Gd(s)||∞ e γ}

K∞(γ) ) {k: kp g 0, ki g 0, kd g 0; ||Gd(s)||∞ e γ} (14) Hence, effort is first devoted to construct the H∞ constraint domain K∞(γ). 3.1. Construction of H∞ Constraint Domain K∞(γ). It follows from eq 6 that a lower bound of the H∞ norm ||Gd(s)||∞ is given by |Gd(jω/)| for any ω/ ∈ [0,∞). Let γ > |Gd(jω/)|. Then, γ e ||Gd(s)||∞ if and only if the function [Gd(s) Gd(-s)]/γ2 - 1 has zeros on the imaginary axis, i.e., [Gd(jω) Gd(-jω)]/γ2 - 1 ) 0 for some ω ∈ [0,∞). Hence, for a specified disturbance attenuation level γ, it follows from eqs 1 and 7 that, if the controller parameters kp, ki, and kd satisfy the equation

Φ(kp,ki,kd,ω,γ) ≡

[ (

)]

ki 1 |Gp(jω)|2 - |1 + kp - j - kdω Gp(jω)|2 ) 0 (15) 2 ω γ for some ω ∈ [0,∞), then γ e ||Gd(s)||∞. In other words, for given kd and γ, ||Gd(s)||∞ < γ if and only if

Hence, the performance index Jm is an H2 norm. Also noted is that the performance indices J0 and J1 are usually referred to as the ISE and ITSE performance indices, respectively. For the generalized quadratic performance index Jm, the time weighting factor tm with m > 0 provides the measure with a flexibility of placing less emphasis on the error during the initial transient period, which can consequently lead to a step response that is less oscillatory than those based on the minimization of an ISE performance index. Because the setpoint tracking error of the feedback control in Figure 1 can be reduced only by adjusting the three PID gain parameters, the aim of PID controller tuning for achieving a good tracking ability is to minimize the performance index Jm. Hence, the optimal PID control design problem, taking both tracking and disturbance rejection into account, can be formulated as the following mixed H2/H∞ optimization problem

min Jm(k)

later, the construction of the feasible gain domain Kf amounts to the construction of the H∞ constraint domain

(13)

and then to search the optimal gain values from the constructed feasible domain Kf that minimizes H2 norm Jm. The next section is devoted to the construction of the feasible gain domain Kf. 3. Construction of Feasible Gain Domain The purpose of this section is to construct the feasible gain domain Kf described by eq 13. As will be shown

(kp, ki) ∉ E(kd,γ) ≡ {(kp, ki): φ(kp,ki,ω;kd,γ) ) 0, ∀ω ∈ [0,∞)} (16) where φ(kp,ki,ω;kd,γ) ) Φ(kp,ki,kd,ω,γ) with kd and γ being regarded as parameters. Moreover, it is obvious that maxω∈[0,∞) |Gd(jω)| ) γ if and only if (ki, kd) ∈ ∂E(kd,γ), where ∂E(kd,γ) denotes the boundary of the H∞ constraint set E(kd,γ). In the following discussion, E(kd,γ) will be referred to as an H∞ constraint set. Let Gp(jω) be represented as

Gp(jω) ) Pr(ω) + jPi(ω) ) Pa(ω)ejθ(ω)

(17)

where Pr(ω), Pi(ω), Pa(ω), and θ(ω) are real functions of ω. Then, the function φ can be written as

φ(kp,ki,ω;kd,γ) )

[

2

(ki-ω kd) 1 2 P (ω) - 1 + 2Pr(ω)kp + 2Pi(ω) + 2 a ω γ 2 ki (18) Pa2(ω)kp2 + Pa2(ω) - kd ω

(

)]

For a fixed frequency ω and a given pair (kd, γ), the above quadratic equation defines an ellipsoidal curve in the kp-ki plane. It is obvious that the union of the curves corresponding to every ω ∈[0,∞) constitutes a band in the kp-ki plane, which is exactly the H∞ constraint set E(kd,γ) defined by eq 16. It is noted that, for a pair (kp, ki) that lies outside the band E(kd,γ) in the kp-ki plane, the corresponding disturbance transfer function Gd(s) has the property that maxω∈[0,∞)|Gd(jω)| < γ. Hence, in the kp-ki plane, the domain K∞(γ) is the region outside the band E(kd,γ) within which the disturbance transfer function Gd(s) is stable. In the following, we apply the notion of principal points29-31 to characterize the boundary of the H∞ constraint set E(kd,γ), i.e., to characterize ∂E(kd,γ). For notational simplicity, we temporarily omit the parameters kd and γ from φ and E. From a geo-

6110


metrical viewpoint, the equation

φ(kp,ki,ω) ) 0

(19)

defines a two-dimensional (2-D) surface S in the threedimensional gain-frequency space, i.e., in the space with kp, ki, and ω as coordinates. The projection of the 2-D surface S onto the kp-ki plane is exactly the H∞ constraint set E defined by eq 16. By regarding the kpki plane as the complex plane, the H∞ constraint set E can be viewed as the image of the surface S under the complex-valued mapping

In a similar manner, we can show that eq 25 also holds if ∂φ/∂ki * 0. Moreover, it follows from eq 23 that the relation in eq 25 holds even if (∂φ/∂kp, ∂φ/∂ki) ) (0, 0). Hence, let the set of points on the surface φ(kp,ki,ω) ) 0 while satisfying the relation ψ(kp,ki,ω) ) 0 for ω ∈ [0,∞). Then, it is clear that the image of P under the mapping g includes the boundary of the value set g(S), that is

∂g(S) ⊂ g(P)

(26)

(23)

It is noted that P is a one-dimensional manifold in the (kp, ki, ω) space. Hence, by tracing this manifold with a path-following algorithm, we can construct the boundary of the H∞ constraint set E ) g(S). Before leaving this subsection, it should be noted that the characterization of the domain boundary ∂K∞(γ) is based on regarding kd as parameter. Alternative characterizations can also be obtained in a similar way by regarding kp or ki as a parameter. Because the boundary equations φ(kp,ki,ω;kd,γ) ) 0 and ψ(kp,ki,ω;kd,γ) ) 0 are both differentiable, the nonempty domain K∞(γ) enclosed by these boundaries is compact but not necessarily convex.44,45 Also noted is that the size of the domain K∞(γ) depends on the value of γ specified. If the specified disturbance attenuation level γ is too small, K∞(γ) could become an empty set. This would mean that no PID controller exists that can reduce the H∞ norm of the disturbance transfer function Gd(s) to the specified level γ. 3.2. Stability Domain Ks. The primary concern in tuning the PID controller Gc(s) is that Gc(s) must be a stabilizing controller. The stability domain Ks is defined by

Assuming that ∂φ/∂kp * 0, we can eliminate dkp from the above two equations to obtain

Ks ) {k: kp g 0, ki g 0, kd g 0; 1 + Gc(σ + jω) Gp(σ + jω) * 0 ∀σ, ω ∈ [0,∞)} (27)

g(kp,ki,ω) ) kp + jki

(20)

E ) g(S) ) {kp + jki: ∀(kp, ki, ω) ∈ S}

(21)

That is

To characterize the boundary of the value set g(S), we first note that the differential of the mapping g is given by

dg )

∂g ∂g ∂g dω dk + dk + ∂kp p ∂ki i ∂ω

) dkp + j dki

(22)

Because (kp, ki, ω) ∈ S, the differential changes dkp, dki, and dω are not independent but satisfy the following relation

dφ(kp,ki,ω) )

∂φ ∂φ ∂φ dk + dk + dω ) 0 ∂kp p ∂ki i ∂ω

( )

∂φ ∂φ ∂ki ∂ω + j dki dω dg ) ∂φ ∂φ ∂kp ∂kp ≡ gidki + gω dω

The critical stability boundaries of Ks satisfy the equation

1 + Gc(jω) Gp(jω) ) 0 (24)

where the signs of dki and dω can be chosen arbitrarily. If the plane vectors gi and gω are not collinear, then for any point v in the neighborhood of g(kp, ki, ω) with |g v| < 0+, we can choose dki and dω such that g + dg ) v. In this case, the image g is not on the boundary of the value set g(S). Hence, the necessary condition for a point (kp, ki, ω) on the surface φ(kp, ki, ω) ) 0 with its image under the mapping g lying on the boundary of the value set g(S) is that the plane vectors gi and gω be collinear. Mathematically, this condition is represented by

[ ]

∂φ ∂φ / ∂ki ∂kp 1 ∂φ ∂φ ) / ∂φ ∂φ ∂ω ∂kp - / ∂ω ∂kp 0

which is equivalent to

[ (

|1 + kp - j

)]

ki - kdω Gp(jω)|2 ) 0 ω

(29)

Because the equation φ(kp,ki,ω;kd,γ) ) 0 defined in eq 15 reduces to the above equation as γ f ∞, it can easily be shown that the H∞ constraint set E(kd,γ) includes the critical boundaries of the stability domain Ks described by eq 29. As a result, we have the following inclusion relation

K∞(γ) ⊂ Ks

-

(28)

(30)

Hence, a PID controller with its three gains lying in the H∞ constraint domain K∞(γ) is guaranteed to be a stabilizing controller.

)0

or

4. Evaluation of the H2 Norm Jm

ψ(kp,ki,ω) ≡

∂φ )0 ∂ω

(25)

In this section, we derive parametric formulas for the evaluation of the H2 norm Jm with m ) 0 and


m ) 1. To this end, let the process to be controlled be described by the transfer function

Gp(s) ) ≡

b0 + b1s + ‚‚‚ + bn-1sn-1 -hs e a0 + a1s + ‚‚‚ + ansn B(s)e-hs A(s)

E(s) )

(31)

B0(s) A0(s) + A1(s)e-hs

L {te(t)} ) - E′(s) A0(s) + A1(s)e-hs

A′0(s) + [A′1(s) - hA1(s)]e-hs

A0(s) + A1(s)e-hs

A0(s) + A1(s)e-hs

≡ -[E0(s) - E(s) E1(s)]

J2 )

A0(s) + A1(s)e-hs A0(-s) + A1(-s)ehs

(34)

+

A0(-s) + A1(-s)ehs

A0(-s) + A1(-s)ehs

≡ Z(s) + Z(-s)

] (36)

where

1 Z0(s) ) A0(s) B0(s) B0(-s) 2

With the product-to-sum expansion in eq 36, the integral J0 becomes

J0 )

)

Y0(-s) + Y1(-s)ehs

A0(s) + A1(s)e

B0(-s) -hs

-hs Z0(-s) + Z1(-s)ehs 1 Z0(s) + Z1(s)e + ∆(s) A0(s) + A1(s)e-hs A0(-s) + A1(-s)ehs

1 j∞ ∫0∞te(t) te(t) dt ) j2π ∫-jωE′(s) E′(-s) ds

B0(s) + B1(s)e-hs C0(-s) + C1(-s)ehs

[

B0(s)

1 Z1(s) ) - A1(s) B0(s) B0(-s) 2

The derivation of expressions for the parametric evaluation of Jm is based on the following product-to-sum expansion33

where

E(s) E(-s) )

(33)

-1 j∞ ∫0∞te(t) e(t) dt ) j2π ∫-jωE(s) E′(-s) ds

A0(s) + A1(s)e-hs

1 A (-s)[B0(s) C0(-s) + 2∆(s) 1 1 A (-s) B0(s) C1(-s) B1(s) C1(-s)] + ∆(s) 0

It can be shown that, if C0(s) ) B0(s) and C1(s) ) B1(s), then Y0(-s) ) X0(-s) and Y1(-s) ) X1(-s). For the purpose of easy reference, the product-to-sum expansion formula in eq 35 is derived in Appendix A. Now, with E(s) given in eq 32, we have the following product-to-sum expansion for the integrand of J0

1 j∞ ∫0∞e(t) e(t) dt ) j2π ∫-jωE(s) E(-s) ds

X0(s) + X1(s)e-hs

1 A (s)[B0(-s) C0(-s) + B1(s) C1(-s)] 2∆(s) 0 1 A (-s) B1(s) C0(-s) ∆(s) 1

Y1(-s) ) -

)

where the prime denotes differentiation with respect to s. If the control system in Figure 1 is internally stable, then the H2 norms Jm for m ) 0, 1, 2 can be expressed as follows

J1 )

Y0(-s) )

1 A (s)[B0(s) C0(-s) + B1(s) C1(-s)] + 2∆(s) 1 1 A (s) B1(s) C0(-s) ∆(s) 0

+

B0(s)

J0 )

X1(s) ) -

(32)

where B0(s) ) A(s), A0(s) ) sA(s), and A1(s) ) (kds2 + kps + ki)B(s). Moreover, the time-weighted error variable te(t) has the following Laplace transform

B′0(s)

1 A (s)[B0(s) C0(-s) + B1(s) C1(-s)] 2∆(s) 0 1 A (s) B0(s) C1(-s) ∆(s) 1 ∆(s) ) A0(s) A0(-s) - A1(s) A1(-s)

where h is the time delay. Here, it is noted that the delay-free part of Gp(s) is assumed to be strictly proper because all chemical processes have the capacity of storing energy and/or materials. With this process transfer function, the error E(s) in the PID control loop due to a unit step set-point change, R(s) ) 1/s, can be written from eq 8 as

)-

X0(s) )

(35)

j∞ 1 Z(-s) ds Z(s) ds + ∫-j∞ j2π

1 j2π

(37)

Because E(s) is assumed to be analytical on the closed half-plane of the s plane, we can deform the integration contour (from -j∞ to j∞) in the right-hand side of eq 37 so that all of the poles of E(s) and all of the zeros of ∆(s) on the negative imaginary axis are to the left of the contour. Hence, with the closed integral contours ΓR and ΓL shown in Figure 2, eq 37 can be rewritten as

J0 ) where

1 1 I Z(s) ds + I Z(-s) ds + SR - SL (38) j2π ΓR j2π ΓL

6112


SR ) lim Γf∞

Γ 2π

Γ 2π

SL ) lim Γf∞

π/2 Z(Γejθ)ejθ dθ ∫-π/2

∫π/23π/2Z(-Γejθ)ejθ dθ

Applying the Cauchy residue theorem32 to the two complex integrals above, we obtain

J0 ) -Res[Z(s),0] -

∑1 Res[Z(s),sk] + ∑2 Res[Z(-s),sk] + SR - SL

(39)

where Res[Z(s),sk] denotes the residue of Z(s) at sk, the first summation is over all zeros sk of ∆(s) that either have positive real parts or lie on the positive imaginary axis (excluding the origin), and the second summation is over all other zeros -sk of ∆(s) different from the origin. Because ∆(-s) ) ∆(s), we have SL ) -SR and

Res[Z(-s),-sk] ) - Res[Z(s),sk)]

(40)

Figure 2. Integration contour.

into eq 42, we can write the product E(s) E′(-s) as

for sk lying inside the integration contour ΓR. Hence, we finally obtain the expression for evaluating J0 as follows

J0 ) -Res[Z(s),0] - 2

∑1 Res[Z(s),sk] + 2SR

(41)

E(s) E′(-s) ) X(s) + Y(-s) where

X(s) ) U(s) - P(s)

Now, we consider the case for evaluation of J1. With the expression in eq 33 for E′(s) and the product-to-sum expansion in eq 36 for E(s) E(-s), the product E(s) E′(-s) can be expanded as

Y(-s) ) V(-s) - Q(-s) - Z(-s) E1(-s) With this expansion, the integral J1 in eq 34 can be evaluated as

E(s) E′(-s) ) E(s) E0(s) - E(s) E(-s) E1(-s) ) E(s) E0(s) - [Z(s) + Z(-s)]E1(-s) (42) ) E(s) E0(s) - Z(s) E1(-s) Z(-s) E1(-s)

J1 ) -

1 j2π

j∞ E(s) E′(-s) ds ∫-j∞

)-

1 j2π

j∞ [X(s) + Y(-s)] ds ∫-j∞

) Res[X(s),0] +

Substituting the product-to-sum expansions

E(s) E0(-s) ) )

A0(s) + A1(s)e-hs A0(-s) + A1(-s)ehs U0(s) + U1(s)e

∆(s)[A0(s) + A1(s)e

where

-hs

SX ) lim

]

∆(s)[A0(-s) + A1(-s)e ] ≡ U(s) + V(-s)

(43)

and

Z(s) E1(-s) ) Z0(s) + Z1(s)e-hs

[A′1(-s) - hA1(-s)]ehs

∆(s)[A0(s) + A1(s)e-hs]

A0(-s) + A1(-s)ehs

P0(s) + P1(s)e-hs ∆2(s)[A0(s) + A1(s)e-hs]

π/2 X(Γejθ)ejθ dθ ∫-π/2

Γ Γf∞ 2π

+

hs

+

Q0(-s) + Q1(-s)ehs ∆2(s)[A0(-s) + A1(-s)ehs] ≡ P(s) + Q(-s)

(46)

-hs

V0(-s) + V1(-s)ehs

)

∑1 Res[X(s),sk] -

∑2 Res[Y(-s),sk] - SX + SY

B0′(-s)

B0(s)

(45)

(44)

SY ) lim Γf∞

Γ 2π

∫π/23π/2Y(-Γejθ)ejθ dθ

Before leaving this section, note that, although expressions for the evaluation of the H2 norms Jm are given here only for m ) 0 and m ) 1, the expressions for the evaluation of J2 or Jm with m > 2 can be obtained by performing sequential product-to-sum expansions similar to those given in eqs 42-46 for E(s) E′(-s). 5. Optimization Method In this paper, the problem of mixed H2/H∞ optimization (eq 12) is solved via the differential evolution algorithm developed by Price and Storn.34 Like other evolutionary algorithms, such as genetic algorithms,35 evolutionary strategy,36 evolutionary programming,36 and simulated annealing,37 the differential evolution algorithm (DEA) belongs to the class of direct, probabilistic search techniques inspired by the process of


Differential evolution (DE) is a parallel direct search method that searches for the n-dimensional parameter vector x ) (x/1, x/2, ..., x/n) over the finite parameter + + + space X ≡ [x1 , x1 ] × [x2 , x2 ] × ‚‚‚ × [xn , xn ] such that the cost functional I(x) evaluated at x/ is a minimum. It utilizes Np n-dimensional parameter vectors

x(i) j , j ) 0, 1, ..., Np - 1

(47)

as a population for each iteration of the minimization. The initial population is chosen randomly and should try to cover the entire parameter space X uniformly. Just as with most evolutionary algorithms, DE generates a better next generation of the population by using the operations of mutation, recombination, and selection. The mutation operation is used to explore the unsearched domain. DE mutates a parameter vector xk to a mutated vector y ) (y1, y2, ..., yn) by adding to it the weighted difference between two population vectors xl and xm, i.e.

k, l, m ∈ {0, 1, ..., Np - 1} (48)

y ) xk + F(xl - xm),

where xl and xm are randomly chosen from the current population and F, which is guided41 to lie in the interval [0.5, 1], is a constant weighting factor that controls the amplification of the differential variation (xl - xm). To ensure that the mutated vector y lies in the parameter space X, the following modification of mutated vector components is needed

{

xj if yj < xj yj ) yj if xj e yj < x+ j + y x+ > x if j j j

Figure 3. Boundaries of feasible domains and stability domains for kd ) 0.1-2.2 with γ ) 0.8.

natural evolution. It simulates the principle of evolution, i.e, survival of the fittest, and deals with a population of potential individuals through repeated applications of “evolutionary” operators such as mutation, crossover, and selection. The DEA also produces individuals with successively improved fitness and converges, hopefully, to the fittest individuals representing optimal solutions. As such, the DEA can be viewed as a general-purpose optimization technique for fitness functions. It has been successfully applied to the solution of global optimization problems related to some concrete multimodal, nondifferentiable, and combinatorial functions.38-40 The adaptation of the differential evolution algorithm used here is motivated by its following advantages: (1) ability to handle nondifferentiable and multimodal objective functions; (2) parallelizability to cope with computation-intensive objective functions; (3) ease of use, i.e., few control variables to steer the minimization, and variables that are robust and easy to choose; (4) good convergence properties, i.e., consistent convergence to the global minimum in consecutive independent trials; and (5) ease of implementation and operation, i.e., use of floating-point numbers to represent points in continuous space and addition to search the continuum.

(49)

To increase the potential diversity of the perturbed parameter vectors, the recombination or crossover operation is used. The crossover of two parameter vectors x ) (x1, x2, ..., xn) and y ) (y1, y2, ..., yn) yields a new parameter vector t ) (t1, t2, ..., tn) as follows

{

x r eC tj ) yi rj > Cr j j r

(50)

where {r1, r2, ..., rn} is a sequence of random numbers rj ∈ (0, 1) and Cr ∈ (0, 1) is the crossover probability. As guided by Storn,41 the crossover probability Cr must usually be less than 1 (e.g., 0.3). If no convergence can be achieved, however, Cr ∈ (0.8, 1) then helps. DE uses a simple selection process for producing better offspring. Each vector xj in the current population is targeted once for recombination with a mutated vector y to produce a trial vector t. The cost or object function of the trial vector t is compared with that of the target vector xj. If the cost of the target vector xj is lower than that of the trial vector, the target vector xj is allowed to advance to the next generation of the population. Otherwise, the trial vector t becomes a member of the next generation of the population. With the above description of differential evolution, we can outline in the following subsection a search algorithm for finding x/ ∈ X that minimizes the cost functional I(x).

6114


Figure 4. Boundaries of the H∞ constraint sets E(kd,0.8) for kd ) 1.0, 1.5, 2.0, and 2.5.

Figure 5. Comparison of the magnitude of the frequency responses of Gd(s) for example 6.1.

Differential Evolution Algorithm. (0) Choose population size Np, weighting factor F, crossover constant Cr, and maximum number of iterations M. (1) Generate an initial population of Np parameter vectors x(0) j ∈ X, j ) 0, 1, ..., Np-1 and evaluate the associated costs I(0) j ) I(x(0) ). (2) Set the iteration index i ) 1. (3) Set the j target index j ) 0. (4) Choose three random integer numbers k, l, m ∈ {0, 1, ..., Np - 1}. (5) Generate a (i) (i) mutated vector y ) x(i) k + F(xl - xm ) and modify y according to eq 49 such that y ∈ X. (6) Generate a sequence of n random parameters rl ∈ (0, 1), l ) 1, 2, ..., n. (7) Obtain a trial vector t from the target vector x ) x(i) j and the mutated vector y by the crossover operation (eq 50). (8) Evaluate the cost of the trial vector (i+1) t and denote it by It. If It > I(i) ) x(i) j , then xj j and

Figure 6. Comparison of unit step responses for example 6.1. (i+1) ) I(i) ) t and I(i+1) ) It. (9) I(i+1) j j ; otherwise, xj j Increment the target index j by 1. If j < Np, then go to step 4. (10) Increment the iteration index i by 1. If i < M, then go to step 3. (11) Let k ) arg min I(M) j . Then, the optimal solution is given by x/ ) x(M) . (12) Stop. k

6. Illustrative Examples This section presents two examples to demonstrate the proposed approach to solving the mixed H2/H∞optimal PI/PID control design problem stated in eq 12. The solution procedure includes the following steps: (i) application of the boundary generation method described in section 3 to obtain the boundary of the sets E(kd,γ) for various values of kd so as to construct the feasible gain domain K∞(γ); (ii) definition of the cost


Figure 7. Feasible domains K∞(γ)0.6) in the kp-kd plane for ki ) 0, 0.01, ..., 0.07.

functional to be minimized as

Im(k) )

{

Jm(k) if k ∈ K∞(γ) 2 Jm(k) + λ[G(k)] if k ∈ K∞(γ)

(51)

where λ is an arbitrarily large number, say, 106 and G(k) is a penalty function, which is usually set to unity; and (iii) application of the differential evolution algorithm described in section 5 to search for the optimal controller gain k that minimizes the cost functional Im(k). It is noted that, by setting the amplification factor F to be less than unity, say, F ) 0.7, in the differential evolution search, the searching domain can be made to shrink and move into the feasible domain K∞(γ) as the search iteration proceeds. As a result, the final solution is independent of the penalty parameter λ. In working out the following example, the onedimensional manifold M described by φ(kp,ki,ω;kd,γ) ) 0 and ψ(kp,ki,ω;kd,γ) ) 0 is traced by the spherical (0) method.42,43 Starting from a given point q(0) ) (k(0) p , ki , (0) ω ) ∈ M, this method traces the manifold M by solving the following equations recursively

φ(q

(k+1)

))0

Figure 8. Comparison of feasible domains K∞(0.6) and K∞(0.8) for ki ) 0.02 and 0.05.

6.1. Example 1. Consider the feedback control system shown in Figure 1, in which the plant transfer function is given by25

Gp(s) )

The aim is to find the settings of a PI/PID controller Gc(s) such that the H2 performance J0 is minimized subject to the following H∞ constraint

Gp(s) || || e γ ) 0.8 1 + Gc(s) Gp(s) ∞ To solve this mixed H2/H∞-optimal PI/PID controller design problem, we first construct the gain feasible domain K∞(γ). For the specified transfer function, Gp(jω) can be written as

Gp(jω) ) Pr(ω) + jPi(ω) where

Pr(ω) )

ψ(q(k+1)) ) 0

(52)

2 2 2 (k+1) (k+1) - q(k) - q(k) - q(k) [q(k+1) 0 0 ] + [q1 1 ] + [q2 2 ] -

2 ) 0 (k) (k) to obtain a sequence of points q(k) ) (q(k) 0 , q1 , q2 ), k ) 1, 2, ..., M, with the distance between two successive points being less than a prespecified value , which is set to 0.01 in the actual computation. In general, the boundaries of the feasible gain domain K∞(γ) cannot be described by simple mathematical expressions. In practice, an approximate gain domain K ˜ ∞(γ) with a piecewise linear boundary is used. To obtain K ˜ ∞(γ), we first grid the parameter kd and then apply the method described in subsection 3.1 to construct the boundary of the H∞ constraint domain K∞(γ) in the kp-ki plane for each gridded kd point. Finally, the approximate domain K ˜ ∞(γ) is reconstructed from the sliced K∞(γ) using the algorithm described in Appendix B.

e-1.5s (s + 1)3

Pi(ω) )

ω(ω2 - 3) sin 1.5ω - (3ω2 - 1) cos 1.5ω (ω2 + 1)3 ω(ω2 - 3) cos 1.5ω - (3ω2 - 1) sin 1.5ω (ω2 + 1)3

With the above expressions for Pr(ω) and Pi(ω), we have the defining equation φ ) 0 in eq 18 for the H∞ constraint set E(kd,γ) with γ ) 0.8. Following the derivations given in section 3, the additional equation required to construct the boundary of the H∞ constraint set E(kd,0.8) is given by

ψ(kp,ki,ω;kd,0.8) )

∂φ )0 ∂ω

Using the spherical algorithm,42,43 we traced the boundary curves of the feasible domains in the kp-ki plane for various values of kd. The traced boundary curves are shown in Figure 3. Figure 4a-d shows the boundaries of the sets E(kd,0.8) for kd ) 1.0, 1.5, 2.0, and 2.5, along with the curves defined by φ ) 0 for various values of

6116


Figure 9. Boundaries of the H∞ constraint sets E(ki,0.8) for ki ) 0.01, 0.05, 0.1, and 0.15. Solid lines, boundaries of E(ki,0.8); dashed lines, boundaries of the stability domain.

ω. Having constructed the feasible domain K∞(0.8), we applied the differential evolution algorithm with control variables F ) 0.6 and Cr ) 0.67 to search for the optimal controller gain in the feasible domain K∞(0.8) that minimizes the H2 norm J0. The obtained settings for the PID controller are given by (k/p, k/i , k/d) ) (0.988, 0.297, 1.5). The minimum H2 norm attained by this optimal PID control is J0 ) 3.306. The magnitude of the frequency response of the disturbance transfer function Gd(s) with the designed controller is shown in Figure 5. It can be seen that the maximum magnitude of the frequency response of the disturbance transfer function Gd(s) is indeed less than the specified level γ ) 0.8. The time response of the closed-loop system subject to a unit step reference input is shown in Figure 6. For comparison, we relaxed the H∞-norm constraint on the disturbance transfer function Gd(s) and performed the minimization of H2 norm J0. The PID controller parameters obtained by minimizing J0 only are given by (kp, ki, kd) ) (1.333, 0.313, 1.802), which renders a minimum H2 norm J0 ) 3.142 and an H∞ norm ||Gd(s)||∞ ) 1.144. The magnitude of the frequency responses and the set-point responses for the designed PID controllers are compared in Figures 5 and 6, respectively. It can be seen that there is a tradeoff between the H2 norm J0 and the H∞ norm ||Gd(s)||∞. Also, it can be seen from Figure 6 that the PID control design with specification of a constraint on the disturbance attenuation level leads to a less oscillatory response than that without specification of an H∞-norm constraint.

6.2. Example 2. Consider a high-vacuum distillation column in which the transfer function between the viscosity and reflux flow is given by46

Gp(s) )

0.57 e-18.7s 2 (8.6s + 1)

This process has a long time delay. In constructing the feasible domain K∞(γ), we grid the parameter ki rather than the parameter kd and find the boundary of the H∞ constraint set ∂E(ki,γ) in the kp-kd plane. Figure 7 shows the feasible domains K∞(γ)0.6) for ki varying from 0.00 to 0.07 in a step of 0.01. It can be seen from this figure that the feasible domain K∞(γ)0.6) is convex for each ki and its size decreases as ki is increased. Figure 8 compares the sizes of the feasible domains K∞(γ)0.6) and K∞(γ)0.8) for ki ) 0.02 and 0.05. As expected, the feasible domain becomes smaller as the disturbance attenuation level decreases (becomes more stringent). Figure 9 shows the boundaries of the H∞ constraint set E(ki,0.8) for ki ) 0.01, 0.05, 0.1, and 0.15. It can be seen in Figure 9d that the feasible domain disappears for ki ) 0.15. With the H∞ constraint level γ ) 0.6, the minimum-ISE PID controller parameters searched from the feasible domain K∞(0.6) are given by (kp, ki, kd) ) (1.3847, 0.0576, 14.007). The corresponding minimum H2 norm is J0 ) 29.345. By relaxing the H∞ constraint, the optimal PID controller parameters are obtained as (kp, ki, kd) ) (1.7565, 0.0597, 16.8508). The H2 performance index for this controller is J0 ) 28.441, which is smaller than that associated with the H∞constraint PID controller. For comparison, we applied


Figure 10. Comparison of the magnitude of the frequency responses of Gd(s) for example 6.2.

evaluating the quadratic cost functional greatly facilitates the serach for the optimal controller parameters that minimize the H2 norm associated with the set-point response. To obtain the global optimal solution, the search was done by the simple heuristic differential evolution algorithm. The effectiveness of the proposed PID design method was demonstrated with two numerical examples. Finally, it should be mentioned that the proposed approach to the design of H2/H∞-optimal PID controllers is viable even if the constraint domain K∞(γ) is nonconvex in the gain space. In such a case, a more sophisticated description of the approximate domain K ˜ ∞(γ) would be required, and an interactive graphical design procedure would be helpful. The latter approach is often adapted for parameter space design methods44,45 because the shape of the feasible domain K∞(γ) is difficult to predict via simple mathematical formulas, particularly for the case of higher-order systems. Acknowledgment This work was supported by the National Science Councils of the Republic of China under Grants NSC89-2214-E-194-014 and NSC-89-2214-E-194-018. Appendix A This appendix contains the derivation of the productto-sum expansion given in eq 35. To begin, we let z ) e-hs and write polynomials A(s) and A(-s) briefly as A and A h , respectively. Thus, we can write eq 35 as

Figure 11. Comparison of unit step responses for example 6.2.

the Ziegler-Nichols tuning method2 to obtain the PID controller settings (kp, ki, kd) ) (1.7527, 0.0529, 14.51658). The associated H2 performance index was computed to be J0 ) 28.584. Figure 10 compares the magnitude of the frequency responses of the disturbance transfer function Gd(s) obtained using the above three PID controllers. The unit step responses of control systems using these three PID controllers are compared in Figure 11. It can be seen in these figures that the disturbance attenuation capability of the designed H2/ H∞ PID controller is superior to that of the controller tuned by the Ziegler-Nichols method. Also noted is that the set-point response of the PID control system tuned by the proposed method exhibits less oscillatory behavior compared to those tuned by the Ziegler-Nichols method and by simple minimization of the H2 norm. 7. Conclusions A parameter space design method has been proposed for PID controllers that minimize the H2 norm associated with the set-point response while satisfying an H∞ norm on the load disturbance rejection. The method is based on an effective approach to the constructino of the feasible controller parameter domain within which the specified H∞-norm constraint is automatically satisfied. The construction of the feasible domain is achieved bu using the notion of principal points to derive an analytic expression to describe the boundaries of an H∞ constraint set defined by the maximum H∞ norm. By tracing the H∞ constraint set boundary using a pathfollowing algorithm, the feasible domain boundaries in the kp-ki parameter plane can be easily obtained for each value of kd. The availability of the feasible controller parameter domain and the parametric formula for

h0 + C h0 + Y B0 + B1z C h 1z-1 X0 + X1z Y h 1z-1 + ) (a1) A0 + A1z A h +A h z-1 A0 + A1z A h +A h z-1 0

1

0

1

Combining the two terms on the right-hand side and then equating the numerators of both sides of the resulting equation, we have

h0 + C h 1z-1) ) (B0 + B1z)(C h 1z-1)(X0 + X1z) + (A0 + A1z)(Y h0 + Y h 1z-1) (a2) (A h0 + A Expanding the above equation and then equating the coefficients of like powers zk, k ) -1, 0, 1, we obtain the following three equations for the four unknowns X0, X1, Y h 0, and Y h1

h1 + A h 1X0 ) B0C h1 A0Y h 1 + A0Y h0 + A h 0X0 + A h 1X1 ) B1C h 1 + B0 C h0 A 1Y A1Y h0 + A h 0X1 ) B1C h0

(a3) (a4) (a5)

Solving eqs a3 and a4 for Y h 1 and X1, respectively, yields

Y h1 ) -

A h1 B0C h1 X0 + A0 A0

(a6)

X1 ) -

A1 B1C h0 Y h0 + A h0 A h0

(a7)

Upon substitution of eqs a6 and a7 for Y h 1 and X1, respectively, eq a4 becomes

6118


A1B0C h 1B1C h1 A h0 ∆ ∆ X0 + Y h0 + + ) B0 C h 0 + B1C h1 A0 A h0 A0 A h0 (a8) where

h 0 - A1A h1 ∆ ) A0A

(a9)

Equation a8 can be split into the following two equations for the unknowns X0 and Y h 0, respectively

A1B0C h1 1 ∆ X0 + ) (B0C h 0 + B1C h 1) A0 A0 2

(a10)

h0 1 A h 1B1C ∆ Y h0 + ) (B0C h 0 + B1C h 1) A h0 A h0 2

(a11)

subintervals. (3) At the lth grid point of kd, the range of + ki is determined to be [ki (l), ki (l)], where l ) 0, 1, ..., + md. (4) the interval [ki (l), ki (l)] is gridded into mi(l) equally spaced subintervals. (5) At the grid point (l, j) of the kd-ki plane, the range of kp is determined to be + [kp (l,j), kp (l,j)], where l ) 0, 1, ..., md and j ) 0, 1, ..., mi(l). On the basis of the above description of the approximate gain domain K ˜ ∞(γ), we determine whether a point (k/p, k/i , k/d) belongs to the approximate gain domain by the following algorithm: Algorithm for Testing Whether (k/p, k/i , k/d) ∈ K ˜ ∞(γ).

begin procedure / + if kd e kd e kd , then ∆d ) (k+ d - kd )/md; l ) integer part of

From eqs a10 and a11, the unknowns X0 and Y h 0 can be readily obtained as

1 A (B C h + B1 C h 1) - A1B0C h1 2 0 0 0 X0 ) ∆

(a12)

1 h + B1C h 1) - A h 1B1C h0 A h (B C 2 0 0 0 Y h0 ) ∆

(k/d - kd )/∆d; / k˜ i ) ki (l) + [ki (l+1) - ki (l)](kd - l∆d)/∆d; + + + / k˜ + i ) ki (l ) + [ki (l+1) - ki (l)](kd - l∆d)/∆d; / ˜+ if k˜ i e ki e k i , then

˜∆i ) (k˜ + i - k i )/mi; r ) integer part of (k/i - k˜ i )/∆i;

(a13) ) kk˜ p,r p (l,r) +

Substituting eqs a12 and a13 into eqs a6 and a7, respectively, we obtain Y h 1 and X1 as follows

1 - A h + B1C h 1) + A h 0B0C h1 h (B C 2 1 0 0 Y h1 ) ∆ 1 - A1(B0C h 0 + B1C h 1) + A0B1C h0 2 X1 ) ∆

/ [kp (l+1,r) - kp (l,r)][kd - (kd + l∆d)]/∆d; ) kk˜ p,r+1 p (l,r+1) + / [kp (l+1,r+1) - kp (l,r+1)][kd - (kd + l∆d)]/∆d;

(a14)

+ ) k+ k˜ p,r p (l,r) + + / [k+ p (l+1,r) - kp (l,r)][kd - (kd + l∆d)]/∆d;

(a15)

+ ) k+ k˜ p,r+1 p (l,r+1) + + / [k+ p (l+1,r+1) - kp (l,r+1)][kd - (kd + l∆d)]/∆d;

Here, it is noted that the split of eq a8 is by no means unique. The two split equations given in eqs 69 and 70 have the advantage that, when C h0 ) B h 0 and C h1 ) B h 1, we have

h0 Y h0 ) X

(a16)

h1 Y h 1) X

(a17)

˜ p,r + (k˜ p,r+1 - k˜ p,r )[k/i - (k˜ k˜ p ) k i + r∆i)]/∆i; + + + ˜ p,r + (k˜ p,r+1 - k˜ p,r )[k/i - (k˜ k˜ + p ) k i + r∆i)]/∆i; / ˜+ if k˜ p e kp e k p , then

˜ ∞(γ) (k/p, k/i , k/d) ∈ K

As a result, we have the following product-to-sum expansion -1

1

(k/p, k/i , k/d) ∉ K ˜ ∞(γ); G(k) )

-1

X0 + X1z X h 1z h 1z B0 + B1z B h0 + B h0 + X ) + A0 + A1z A + A z-1 A0 + A1z A h +A h z-1 0

else

0

/ ˜+ min(|k/p - k˜ p |, |kp - k p |)

(a18)

1

else (k/p, k/i , k/d) ∉ K ˜ ∞(γ); G(k) )

Appendix B This appendix outlines the construction of the approximate H∞ constraint domain K ˜ ∞(γ) and provides an algorithm to test whether a given point (kp, ki, kd) belongs to K ˜ ∞(γ). The construction of K ˜ ∞(γ) is summarized as follows: + (1) The range of kd is determined to be [kd , kd ]. (2) The + interval [kd , kd ] is gridded into md equally spaced

/ ˜+ min(|k/i - k˜ i |, |ki - k i |)

else (k/p, k/i , k/d) ∉ K ˜ ∞(γ); G(k) ) / ˜+ min(|k/d - k˜ d |, |kd - k d |)

end procedure


Literature Cited (1) A° stro¨m, K. J.; Ha¨gglund, T. PID Controllers: Theory, Design, and Tuning, 2nd ed.; ISAsThe Instrumentation, Systems, and Automation Society: Research Triangle Park, NC, 1995. (2) Ziegler, J. G.; Nichols, N. B. Optimum Settings for Automatic Controllers. Trans. ASME 1942, 64 (11), 759. (3) Blickley, G. Modern Control Started with Ziegler-Nichols Tuning. Control Eng. 1990, 37, 11. (4) Bobal, V. Self-tuning Ziegler-Nichols PID Controller. Int. J. Adaptive Control Signal Process. 1995, 9 (2), 213. (5) De Paor, A. M.; O’Malley, M. Controllers of Ziegler-Nichols Type for Unstable Process with Time Delay. Int. J. Control 1989, 49 (4), 1273. (6) De Paor, A. M. Fiftieth Anniversary Celebration of the Ziegler-Nichols PID Controller. Int. J. Electr. Eng. Educ. 1993, 30 (4), 303. (7) Hang, C. C.; A° stro¨m, K. J. Refinements of the ZieglerNichols Tuning Formula for PID Autotuners. Adv. Instrum. 1988, 43, 1021. (8) Hang, C. C.; A° stro¨m, K. J.; Ho, W. K. Refinements of the Ziegler-Nichols Tuning Formula. IEE Proc. D: Control Theory Appl. 1991, 138 (2), l11. (9) Mantz, R. J.; Tacconi, E. J. Complementary Rules to Ziegler and Nichols’s Rule for a Regulating and Tracking Controller. Int. J. Control 1989, 49 (5), 1465. (10) Ott, M. G.; Wojsznis, W. K. Autotuning: From ZieglerNichols to Model Based Rules. Adv. Instrum. Control 1995, 50, 323. (11) Wang, Q. G.; Lee, T. H.; Zhang, Y. Multiloop Version of the Modified Ziegler-Nichols Method for Two Input-Two Output Processes. Ind. Eng. Chem. Res. 1998, 37 (12), 4725. (12) Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. (13) Rivera, E. E.; Morari, M.; Skogestad, S. Internal Model Control. 4 PID Controller Design. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 252-265. (14) Chien, I. L.; Fruehauf, P. S. Consider IMC Tuning to Improve Controller Performance. Chem. Eng. Prog. 1990, 86 (10), 33-41. (15) Graham, D.; Lathrop, R. C. The Synthesis of Optimum Transient Response: Criteria and Standard Forms. AIEE Trans. II: Appl. Ind. 1953, 72, 273. (16) Zhuang, M.; Atherton, D. P. Automatic tuning of optimum PID controllers. IEE Proc. D: Control Theory Appl. 1993, 140 (3), 216-224. (17) Grimble, M. J. H∞ Controllers with A PID Structure. Trans. ASME 1990, 112 (3), 325. (18) McFarlane, D. C.; Glover, K. Robust Controller Design Using Normalized Coprime Factor Plant Descriptions; SpringerVerlag: Berlin, 1990. (19) McFarlane, D. C.; Glover, K. A Loop-Shaping Design Procedure Using H∞ Synthesis. IEEE Trans. Autom. Control 1992, 37 (6), 759. (20) Grassi, E.; Tsakalis, K. S. PID Controller Tuning by Frequency Loop-shaping: Application to Diffusion Furnace Temperature Control. IEEE Trans. Control Syst. Technol. 2000, 8 (5), 842. (21) Panagopoulos, H.; A° stro¨m, K. J. PID Control Design and H∞ Loop Shaping. Int. J. Robust Nonlinear Control 2000, 10 (15), 1249. (22) Tan, W.; Liu, J.; Tam, P. K. S. PID Tuning Based on Loopshaping H∞ Control. IEE Proc. Control Theory Appl. 1998, 145 (6), 485. (23) Chen, B. S.; Cheng, Y. M.; Lee, C. H. A Genetic Approach to Mixed H2/H∞ Optimal PID Control. IEEE Control Syst. Mag. 1995, 15 (5), 51. (24) Krohling, R. Design of PID Controller for Disturbance Rejection: A Genetic Optimization Approach. In Second International Conference on Genetic Algorithms in Engineering Systems: Innovations and Applications; IEE/INSPEC: Edison, NJ, 1997; p 498.

(25) Panagopoulos, H.; A° stro¨m, K. J.; Ha¨gglund, T. Design of PID Controllers Based on Constrained Optimization. In Proceedings of the 1999 American Control Conference; IEEE Press: Piscataway, NJ, 1999; Vol. 6, p 3858. (26) A° stro¨m, K. J.; Panagopoulos, H.; Ha¨gglund, T. Design of PI Controllers Based on Non-Convex Optimization. Automatica 1998, 34 (3), 585. (27) Bao, J.; Forbes, J. F.; McLellan, P. J. Robust Multiloop PID Controller Design: A Successive Semidefinite Programming Approach. Ind. Eng. Chem. Res. 1999, 38 (9), 3407. (28) Zheng, F.; Wang, Q. G.; Lee, T. H. On the Design of Multivariable PID Controllers via LMI Approach. Automatica 2002, 38 (3), 517. (29) Polyak, B. T.; Kogan, J. Necessary and Sufficient Conditions for Robust Stability of Linear Systems with Multiaffine Uncertainty Structure. IEEE Trans. Autom. Control 1995, 40 (7), 1255-1260. (30) Chen, J. J.; Hwang, C. Value Set of Polynomial Families with Coefficients Depending Nonlinearly on Perturbed Parameters. IEE Proc. D: Control Theory Appl. 1998, 145 (1), 304-307. (31) Hwang, C.; Hsiao, C. H. Solution of a Non-Convex Optimization Arising in PI/PID Control Design. Automatica 2002, 38 (11), 1895-1904. (32) Conway, J. B. Functions of One Complex Variable, 2nd ed; Springer-Verlag: New York, 1978. (33) Marshall, J. E.; Gorecki, H.; Korytowski, A.; Walton, K. Time-Delay Systems Stability and Performance Criteria with Applications; Ellis Horwood: New York, 1992. (34) Price, K.; Storn, R. Differential Evolution. Dr. Dobb’s J. Software Tools Prof. Programmer 1997, 22 (4), 18. (35) Golderg, D. E. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1899. (36) Back, T. Evolutionary Algorithms in Theory and Practice; Oxford University Press: New York, 1996. (37) Azencott, R. Simulated Annealing; John Wiley & Sons: New York, 1992. (38) Price, K.; Storn, R. Differential EvolutionsA Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Global Optim. 1997, 11, 341. (39) Wang, F. S.; Chiou, J. P. Optimal Control and Optimal Time Location Problems of Differential-Algebraic Systems by Differential Evolution. Ind. Eng. Chem. Res. 1997, 36 (12), 5348. (40) Lee, M. H.; Han, C. H.; Chang, K. S. Dynamic Optimization of a Continuous Polymer Reactor Using a Modified Differential Evolution Algorithm. Ind. Eng. Chem. Res. 1999, 38 (12), 4825. (41) Storn, R. On the usage of differential evolution for function optimization. In Proceedings of the 1996 Biennial Conference of the North American Fuzzy Information Processing Society; IEEE Press: Piscataway, NJ, 1996; p 519. (42) Yamamura, K. Simple Algorithms for Tracing Solution Curves. IEEE Trans. Circuits Syst. 1993, 40 (8), 537. (43) Yamamura, K.; Sekiguchi, T. A Modified Spherical Method for Tracing Solution Curves. IEICE Trans. Fundam. 1995, E78-A (9), 1233. (44) Datta, A.; Ho, M. T.; Bhattacharyya, S. P. Structure and Synthesis of PID Controllers; Springer-Verlag: New York, 2000. (45) Besson, V.; Shenton, A. T. Interactive control system design by a mixed H∞-parameter space method. IEEE Trans. Autom. Control 1997, 42 (7), 946. (46) He, J. B.; Wang, Q. G.; Lee, T. H. PI/PID controller tuning via LQR approach. Chem. Eng. Sci. 2000, 55, 2429.

Received for review November 29, 2001 Revised manuscript received August 16, 2002 Accepted September 19, 2002 IE010963M

A New Approach to Mixed - American Chemical Society

Recommend Documents