Ind. Eng. Chem. Res. 1998, 37, 4725-4733
4725
Multiloop Version of the Modified Ziegler-Nichols Method for Two Input Two Output Processes Qing-Guo Wang,* Tong-Heng Lee, and Yu Zhang Department of Electrical Engineering, National University of Singapore, Singapore 119260, Singapore
In this paper, a multiloop version of the modified Ziegler-Nichols method1 is proposed for two input two output (TITO) processes. The key idea is to regard each loop together with its corresponding interactions from all other loops as an equivalent single input, single output (SISO) plant and design for it a SISO controller by shifting a given point on the Nyquist curve of the equivalent SISO plant to a desired position. This is, however, a nonlinear problem because the equivalent plants depend on unknown controllers. A novel approach is thus presented to obtain the controller parameters which achieve the above design objective. Unlike many other methods that emphasize the suppression of interactions, our method is developed to channel the effect of interactions to individual loops to increase the speed of loop responses. Simulation examples are provided to show the effectiveness of the proposed method and comparisons are made with the biggest log modulus tuning (BLT) method. 1. Introduction Processes with inherently more than one variable at the output to be controlled are frequently encountered in the industries and are known as multivariable or multi-input, multioutput (MIMO) processes. Interactions usually exist between control loops, which account for the renowned difficulty in their control compared with single input, single output (SISO) processes. The goal of designing a controller to achieve satisfactory loop performance has hence posed a great challenge. Depending on the application and requirement, either fullcross coupled or multi-loop controllers can be adopted for MIMO processes. Although multivariable controllers can provide explicit suppression of interactions, their designs are usually more complex and their implementation is inevitably more costly. Multi-loop controllers, sometimes known as decentralized controllers, have much simpler structures and fewer tuning parameters to handle than the full-cross coupled controllers. In addition, in the event of actuator or sensor failure, multiloop controllers are relatively easy to stabilize manually, because only one loop is directly affected by the failure.2,3 Hence for processes with modest interactions, multi-loop controllers are often more favored than multivariable controllers. Many multi-loop controller design methods have been reported in the literature. In the biggest log modulus tuning (BLT) method4 presented by Luyben, the familiar Ziegler-Nichols rule is modified with the inclusion of a detuning factor, which determines the tradeoff between stability and performance of the system. There, individual controllers are designed for their respective loops by first ignoring all interactions. The calculated controller gains are then scaled by the detuning factor to guarantee stability. Despite simple computations involved, the design regards interactions as elements obstructing system stability and attempts to dispose of them rather than control them to increase the speed of * To whom all correspondence should be addressed. E-mail:
[email protected]. Telephone: (+65)-874 2282. Fax: (+65)-779 1103.
individual loops. It is hence too conservative to exploit process structures and characteristics for best achievable performance. In the sequential loop closing method, the loops are closed one after another, and those previously closed are fitted with appropriate controllers. The main drawback of this method is that the design must proceed in a very ad hoc manner. Design decisions made when the first one or two loops are closed may have deleterious effects on the behavior of the remaining loops. The interactions are well taken care of only if the loops are of considerably different bandwidths and the closing sequence starts from the fastest loop. This assumption of high gains in the loops that have already been closed can rarely be justified, except at low frequencies.5 Generally, the closing of a new loop will bring its interaction into all the previously closed loops. In this paper, we present a new method for the design of multi-loop controllers, which has none of the above drawbacks and avoids the conservativeness of the BLT method because all possible interactions from both loops are handled simultaneously. It is a multi-loop extension of the original Astrom and Hagglund’s modified Ziegler-Nichols method, the idea of which is to move a given point on the process Nyquist curve to a desired position.1 In face of loop interactions, a direct application of the method to the diagonal loop elements in a SISO way may not be appropriate. To take into account the multivariable interactions, we view each loop as an independent equivalent process with all possible interactions lumped into it and design controllers to meet the specifications that a given point on the Nyquist curve be moved to a desired position for each equivalent process. A novel approach is presented to obtain the controller parameters which achieve the above design criterion. Simulation examples are given to show the effectiveness of our method and comparisons are made with the BLT method. The paper is organized as follows. In Section 2, the SISO modified Ziegler-Nichols method1 is reviewed and the selection of the “optimal” desired position is discussed on the basis of ref 6. Section 3 sets the design equations for the multi-loop case of the modified Zie-
10.1021/ie980313r CCC: $15.00 © 1998 American Chemical Society Published on Web 11/19/1998
4726 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
gler-Nichols method. A graphic procedure to solve these equations is developed in Section 4. Simulation examples are given in Section 5. Conclusions are drawn in Section 6. 2. The Modified Ziegler-Nichols Method Two classical tuning methods for proportional and integral/proportional-integral-derivative (PI/PID) controllers were proposed by Ziegler and Nichols.7 These methods are still commonly used, either in their original forms or in modifiied forms. Astrom and Hagglund1 interpreted the Ziegler-Nichols frequency domain method as finding controller parameters so that the critical point, defined as the point where the Nyquist curve intersects the negative real axis, is moved to the point of -(0.6 + j0.28). A general framework which can incorporate such a method and others is formulated and named in ref 1 as the modified Ziegler-Nichols method (MZN). It is described briefly as follows. Given a point A on the Nyquist curve of the process g(s):
A } g(jω) ) raej(-π+φa) find a controller k(s) to compensate for g(s) such that this point is moved to
B } g(jω) k(jω) ) rbej(-π+φb)
(1)
When the controller is chosen as
(
k(s) ) kc 1 +
)
1 + TDs TIs
a solution to eq 1 is
kc ) ωTD -
rb cos(φb - φa) ra
1 ) tan(φb - φa) ωTI
(2) (3)
Note that the controller gain kc is uniquely given by eq 2, but TI and TD are not. An additional condition must thus be introduced to determine these two latter parameters uniquely. A common method is to specify a relation between TI and TD as
TD ) RTI
(4)
where R is recommended8 to be equal to 0.25. Another method to specify R is given below by eq 10, which can generally yield better performance. In many practical cases such as high noise, it is desirable not to use derivative action and a PI controller:
(
k(s) ) kc 1 +
)
β ω
{
tan(φb - φa) + x4R + tan2(φb - φa) R > 0 2R β) 1 R)0 tan(φb - φa) and TD is then given by eq 4. To illustrate the generality of the above framework, note that the Ziegler-Nichols method corresponds to φa ) 0, rb ) 0.66, and φb ) 0.44 for PID and φa ) 0, rb ) 0.4078, and φb ) -0.1964 for PI. The refined Ziegler-Nichols method9 implies φa ) 0, and the values for rb and φb are related to the process normalized gain or normalized dead time through a certain heuristic rule. Design methods in which the specifications are given in terms of gain margin or phase margin can also be interpreted in this general framework. The AstromHagglund phase margin method10 is obtained by setting rb ) 1 and φb ) φm, where φm is the specified phase margin. The choice of φb ) 0 and rb ) 1/Am, where Am is the gain margin, yields their gain margin method. Astrom and Hagglund8 use the relay feedback to estimate the critical point and move it to a desired position on the complex plane, 0.5e-j3π/4, which corresponds to φa ) 0, rb ) 0.5, and φb ) π/4. This design method is considered to be a combination of gain and phase margin specification. Therefore, all of these methods fall into the MZN framework. However, note that the performances of these tuning methods could vary significantly. This is not surprising because a critical point is moved to different desired positions in these methods. In addition, some of these tuning methods may not produce satisfactory closed-loop responses in certain circumstances. For instance, the Ziegler-Nichols tuning laws usually result in rather oscillatory set-point responses and the Astrom-Hagglund phase margin method may not be appropriate to tune PID controllers for processes with a large time delay.6 It is thus important for one to know where to move the critical point for best performance with the help of some knowledge of the process. Zhuang and Atherton6 discuss tuning of the PI/PID controller to achieve optimum integral time-weighted square error (ITSE) performance criterion. Their optimum tuning can be interpreted as the movement of the critical point to a desired position specified by
(
φb ) -arctan
)
1 0.166π(1.935κ + 1)
rb ) 0.361/cos φb
(6) (7)
for the PI controller, or
1 TIs
is preferred. The PI controller can be viewed as a special case of a PID controller with TD ) 0, or R ) 0. It is clear that if a PI controller is adopted, φb - φa should be chosen negative to make the specification achievable, because a PI controller always gives phase lag. With eq 4, TI is obtained from eq 3 as
TI )
where
(5)
φb ) 0.59(1 - 0.97e-0.45κ)
(8)
rb ) 0.614(1 - 0.233e-0.347κ)
(9)
for the PID controller, where κ is the process-normalized gain:
κ)
|
g(j0) g(jωc)
|
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4727
Figure 1. Equivalent g1(s) for a 2 × 2 plant.
and ωc is the process critical frequency. The parameter R for PI and PID controllers is given by R ) 0 and
R)
0.413 3.302κ + 1
(10)
respectively. According to our experience, the modified ZieglerNichols method plus the above optimum choice of the desired point B gives best performance among all the existing methods using information only on static gains and critical points of the process. In the subsequent two sections, the MZN will be extended to the multi-loop case for the control of multivariable plant, where the use of eqs 6 and 10 is optional, but recommended.
Consider a stable 2 × 2 process described via the transfer function matrix G(s) as Y(s) ) G(s)U(s) or
][ ]
y1(s) g (s) g12(s) u1(s) ) 11 y2(s) g21(s) g22(s) u2(s)
Assume that proper input-output pairing has been made to the process. If the process is inherently poorly paired, the relative gain array (RGA) method5 can be employed to make the necessary arrangement. The process is to be controlled in a negative feedback configuration by the multi-loop controller:
K(s) )
[
k1(s) 0 k2(s) 0
]
The resultant control system is shown in Figure 1. The controller design objective is to find K(s) such that both loops achieve satisfactory performance. The boxed portion in Figure 1 can be viewed as an individual SISO process with an equivalent transfer function g1(s) between input u1 and output y1. It follows that g1(s) can be obtained5 as
g1 ) g11 -
g12g21 k2-1 + g22
(11)
Similarly, the equivalent plant between u2 and y2 is given by
g2 ) g22 -
Ai ) gi(jωi) ) raiej(-π+φai)
i ) 1, 2
g21g12 k1-1 + g11
(12)
(13)
be the two given points on the Nyquist curves of gi(s), i ) 1,2. They are to be moved respectively to the desired points
Bi ) gi(jωi)ki(jωi) ) rbiej(-π+φbi)
3. Setup of Design Equations
[ ][
It is thus clear that in order to achieve good loop performance, k1(s) and k2(s) should be designed for g1(s) and g2(s) instead of the plant diagonal element g11(s) and g22(s), as is done in many other methods.4 In the sequential loop closing method,5 k1(s) is first designed for g11(s) (assuming that loop closing starts from loop one). After that, k2(s) is designed for the equivalent transfer function g2(s). While the performance of loop two will be consistent with what is expected from the design, the performance of loop one will not, because the closing of the second loop will bring its interactions into the first loop, which has not been taken into account when k1(s) is designed. In this paper, the modified Ziegler-Nichols method described in the preceding section is applied to the equivalent transfer function g1(s) and g2(s); i.e., the controller ki(s), i ) 1,2 is designed such that a given point on the Nyquist curve of gi(s), i ) 1,2 is moved to a desired point. It is recommended that the selection of these two points follow the same rules as those in the preceding section. In this way, the optimum settings for the multi-loop controller can be obtained with all the benefits derived from the well-developed SISO tuning. Let
i ) 1, 2
(14)
Note that, unlike the SISO case, because of the dependence of g1 (or g2) on k2 (or k1), ω1 and thus k1 (or ω2 and k2) cannot be determined until k2 (or k1) has been fixed. This forms a circle and causes a major design difficulty. Let k1(s) and k2(s) be of PID type, i.e.,
(
ki(s) ) kci 1 +
)
1 + TDis TIis
i ) 1, 2
which reduces to PI type when TDi ) 0. Now the modified Ziegler-Nichols method specified in eq 2, 4, and 5 is applied to each ki(s). Invoking first eq 4 and eq 5 to each ki(s) yields
(
ki(s) ) kci 1 +
)
ωi βi + Ri s βis ωi
i ) 1, 2
(15)
{
where
tan(φbi - φai) + x4Ri + tan2(φbi - φai) R > 0 i 2Ri βi ) 1 Ri ) 0 tan(φbi - φai) (16) and Ri > 0 is for PID and Ri ) 0 is for PI. However, eq 2 cannot be used directly because kc1 is related to ra1 which depends on the unknown kc2. It follows from the phase part of eq 14 that
ki(jωi) ) kci(1 + j tan(φbi - φai))
i ) 1, 2
4728 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
Then, eq 14 can be written as
4. Solving the Equations Consider eq 17 first. Because the equation is complex, it can be broken down to two real equations:
gi(jωi)kci ) cos(φbi - φai)rbiej(-π+φai) Bringing eqs 11, 12, and 15 in gives
{ {
g11(jω1) -
} }
Im[∆1]kc1kc2 + Im[a1]kc1 + Im[b1]kc2 + Im[c1] ) 0
g12(jω1)g21(jω1) kc1 ) 1 + g22(jω1) ω2 β2ω1 1 + j R2 ω2 β2ω1
[ (
kc2
Re[∆1]kc1kc2 + Re[a1]kc1 + Re[b1]kc2 + Re[c1] ) 0
)]
cos(φb1 - φa1)rb1ej(-π+φa1)
g22(jω2) -
g12(jω2)g21(jω1) kc2 ) 1 + g11(jω2) ω1 β1ω2 1 + j R1 ω1 β1ω2
[ (
kc1
)]
cos(φb2 - φa2)rb2ej(-π+φa2)
where
a1 )
g11(jω1) ω2 β2ω1 1 + j R2 ω2 β2ω1
[ (
)]
b1 ) -cos(φb1 - φa1)rb1ej(-π+φa1)g22(jω1) cos(φb1 - φa1)rb1ej(-π+φa1) c1 ) ω2 β2ω1 1 + j R2 ω2 β2ω1
[ (
)]
kc1 can be expressed in terms of kc2 using the first equation as
or
kc1 ) -
g11(jω1) kc1 + β2ω1 ω2 1 + j R2 ω2 β2ω1
[ (
)]
Re[∆1]kc2 + Re[a1]
(19)
The expression is next substituted into the second equation to obtain
cos(φb1 - φa1)rb1ej(-π+φa1)g22(jω1)kc2 + ∆1kc1kc2 ) j(-π+φa1)
cos(φb1 - φa1)rb1e β2ω1 ω2 1 + j R2 ω2 β2ω1
[ (
)]
(17)
g22(jω2) kc2 + ω1 β1ω2 1 + j R1 ω1 β1ω2
[ (
Re[b1]kc2 + Re[c1]
The above equation is quadratic on kc2 and has two solutions. First, only real solutions should be taken for physical realizability. This requires
)]
cos(φb2 - φa2)rb2ej(-π+φa2)g11(jω2)kc1 + ∆2kc2kc1 ) cos(φb2 - φa2)rb2ej(-π+φa2) (18) ω1 β1ω2 1 + j R1 ω1 β1ω2
[ (
(Re[∆1] Im[b1] - Im[∆1] Re[b1])kc22 + (Re[∆1] Im[c1] - Im[∆1] Re[c1] + Re[a1] Im[b1] Im[a1] Re[b1]) kc2 + (Re[a1] Im[c1] Im[a1] Re[c1]) ) 0 (20)
)]
where
∆1 ) g11(jω1)g22(jω1) - g12(jω1)g21(jω1) ∆2 ) g11(jω2)g22(jω2) - g12(jω2)g21(jω2) When φai, φbi, rbi, and Ri are specified, there are altogether four unknowns, namely kc1, kc2, ω1, and ω2, in eqs 17 and 18. These two equations are complex and they can be broken down to four real equations. Because the number of unknowns equals the number of real equations, intuitively, we may expect that the problem can be easily solved. Unfortunately this is not the case. The main difficulty lies in the nonlinearities of the equations which are further complicated by the coupling among the equations. In what follows, a novel procedure for solving the above equations is presented with a graphic approach.
(Re[∆1] Im[c1] - Im[∆1] Re[c1] + Re[a1] Im[b1] Im[a1] Re[b1])2 - 4(Re[∆1] Im[b1] - Im[∆1] Re[b1]) × (Re[a1] Im[c1] - Im[a1] Re[c1]) g 0 (21) Because a1, b1, c1, and ∆1 are functions of ω1 and ω2, eq 21 restricts (ω1, ω2) to some admissible region, and its search is not a simple task. Our strategy here is to make an initial guess on ω2, say, ω2(0), and then find all ω1 to meet eq 21. Such ω1 will be in a union of intervals. We find that all of the intervals are not needed, but only one is sufficient for our design. Notice that when kc2 ) 0, i.e., the second loop is kept open, g1(s) is just g11(s) as seen from eq 11. According to eq 13, ω1 should be the ω11, satisfying
arg g11(jω11) ) -π + φa1
(22)
By changing the controller gain kc2, a different amount of interactions is brought into the first loop and thus g1(s) also changes. When kc2 varies from -∞ to +∞, the corresponding ω1 values in eq 21 form a set. As from eq 11, g1(s) is continuous on kc2, and the set can be expected to be an interval containing ω11. Thus, this j 1], is sufficient particular interval, denoted by [ω1,ω
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4729
for our design because it already covers all possible gains. This interval is here called the principal interval of eq 21. Apart from this principal interval, there are possibly many other intervals which satisfy eq 21. However, the corresponding g1(jω1) will correspond to a point with the phase of, say, -3π + φa1. If such a point is used in design, it may lead to the violation of the Nyquist criterion, although it seems that a given point has been moved to a “safe position”. In view of the above analysis, only the principal j 1] is needed, though eq 21 yields more. The interval [ω1,ω interval is found by starting with ω1 ) ω11 and increasing and decreasing ω1 until the equality in eq 21 is reached. Once the principal interval is determined, eq 20 is solved at each ω1 in the interval to obtain two real solutions for kc2. Note that some solutions may cause an unstable closed loop and thus must be discarded. It follows from the Nyquist stability criterion that for closed-loop stability, kc2 must be positive if g2(0) > 0, or kc2 has to be negative if g2(0) < 0. Because of the presence of integrators in k1 and k2, one sees from eq 12 that
g2(0) ) g22(0) -
g21(0)g12(0) g11(0)
(23)
Thus, the correct solutions for kc2 are those which make g2(0)kc2 > 0 and are used for our design. The rest are discarded. With correct kc2, the corresponding kc1 is readily obtained from eq 19. When ω1 varies in the principal interval for a given ω2(0), the resultant solutions are denoted by kc1(ω1,ω2(0)) and kc2(ω1,ω2(0)). They are plotted in a kc1 - kc2 plane as a curve C1:
j 1]} C1 ) {(kc1(ω1,ω2(0)), kc2(ω1,ω2(0)))|ω1 ∈ [ω1,ω
(24)
So far, only eq 17 is considered. A similar discussion can be made regarding eq 18 and the point of (kc1(ω1(0),ω2), kc2(ω1(0),ω2)) varying with ω2 in the principal j 2], defined similarly to [ω1,ω j 1] will form interval [ω2,ω the second curve C2:
C2 ) {(kc1(ω1(0),ω2), kc2(ω1(0),ω2))|ω2 ∈ [ω2,ω j 2]} (25) The intersection of the two curves, C1 and C2, can be found as
(kc1(ω1*,ω2(0)), kc2(ω1*,ω2(0))) ) (kc1(ω1(0),ω2*), kc2(ω1(0),ω2*)) } Z (26) At the intersection Z, the gains are equal, i.e., kc1(ω1*,ω2(0)) ) kc1(ω1(0),ω2*) and kc2(ω1*,ω2(0)) ) kc2(ω1(0),ω2*). If their arguments are also equal, that is, if ω1* ) ω1(0) and ω2* ) ω2(0), then kci and ωi* would be a solution of the simultaneous eqs 17 and 18 and would fulfill our design. Unfortunately, this is usually not the case; that is, ω1* * ω1(0) and ω2* * ω2(0), in general. To search for a solution, we may iterate on ω1 and ω2 by substituting ω1* for ω1(0) and ω2* for ω2(0) and repeating the above procedure one after another. However, extensive simulations indicate that as long as the initial ω1(0) and ω2(0) are chosen in their principal intervals, the approximate solution in eq 26 with ω1 ) ω1* and ω2 ) ω2* can yield satisfactory control performance. Any iterations will generate very little improvement and thus are not worth carrying out. Therefore, we recom-
mend that the algorithm start with ω1(0) ) ω11 and ω2(0) ) ω22, and stop at Z to produce an approximate solution
ω1 ) ω1*, ω2 ) ω2*
(27a)
kc1 ) kc1(ω1*,ω2(0))
(27b)
kc2 ) kc2(ω1*,ω2(0))
(27c)
The above development can be summarized as the design procedure below, where g1(0) and ω22 are defined respectively by
g1(0) ) g11(0) -
g12(0)g21(0) g22(0)
(28)
and
arg g22(jω22) ) -π + φa2
(29)
Multi-loop Modified Ziegler-Nichols Tuning Procedure. Given data: G(s), type of ki (PI or PID), and the specifications (φai, φbi, and rbi). (1) For PI controller, set Ri ) 0; for PID controller, set Ri ) 0.25 or the value from eq 31 or the user’s own choice. Calculate βi from eq 16. (2) Determine ω11 from eq 22. Take ω2(0) ) ω22 from j ] of eq 21. Solve eq 22, find the principal interval [ω1,ω eq 20 for kc2 for ω1 in the principal interval with kc2 satisfying g2(0)kc2 > 0. Obtain the corresponding kc1 from eq 19. Plot the resultant solutions (kc1(ω1,ω2(0)), kc2(ω1,ω2(0))) in the kc1 - kc2 plane to form the curve C1. (3) Draw C2 similarly to (2). (4) Determine the intersection Z of C1 and C2 and find ωi and kci from eq 27. (5) Set TIi ) βi/ωi and TDi ) RiTIi. In addition to G(s), the above procedure must be supplied with other information on given and desired points and type of controllers for the run. They should be specified correctly to generate a reasonable solution. The detailed instructions and guidelines are presented. Consider first the specification of the given points Ai. For an SISO process, it is generally accepted that the process critical point provides vital information for controller design and thus it is recommended that the points Ai be chosen as the critical points of the equivalent processes gi(s). Note that for gi(0) > 0, the critical point is the place where the Nyquist curve intersects the negative real axis, so that we have φai ) 0. If gi(0) < 0, it is quite clear that φai should equal -π, instead. Though the above are generally recommended, other points can also be chosen as Ai if desired. Importantly, the sign of gi(0) should still be taken into account in this matter to position points Ai correctly, as just done with regard to the critical point. Next, we consider the selection of the points Bi or the choice of φbi and rbi. We recommend that the users follow the rules of eqs 6-9 in Section 2. These parameters are given there in terms of the process normalized gain, which is unknown for gi(s) prior to our multi-loop controller design. We have to approximate by
κi ≈
|
gi(0)
gii(jωcii)
|
(30)
4730 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
Figure 2. Curves C1, C2, C1New, and C2New in Example 1 for PID design.
where ωcii is the critical frequency of gii(s) and gi(0) is readily evaluated from eq 23 and eq 28. In practice, a perfect plant model is never obtainable and hence the robustness of the proposed tuning method should be addressed. Note that the specification of the desired point (or equivalently φbi and rbi) can be viewed as a combination of gain and phase margin specifications. Different desired points will yield different closed-loop robustness. In the presence of large modeling errors or severe nonlinearity, it is advisable to choose the desired point which corresponds to large gain and phase margins for the purpose of stable and robust closed-loop performance. The last question to answer is whether we should launch a PI or a PID controller for a given process. Generally, a PID controller may give more tight control of the loop than does a PI controller. For PID, Ri values are recommended to be
Ri )
0.413 3.302κi + 1
(31)
However, for processes which can be approximated with a first-order model, a PI controller is adequate8 and derivative action is of little use for these processes. Usually, the derivative term is also turned off when the process has a large noise or a long dead time. In the context of multi-loop control, a fairly useful analysis tool for the assessment of the process controllability is the Niederlinski index4. The positiveness of the Niederlinski index is necessary for the integral stability of the system, and a process with negative Niederlinski index will be either wrongly paired or inherently difficult to control by a multi-loop controller. For such processes, our method may not generate satisfactory performance and is thus not recommended for use. It can be easily shown that if the Niederlinski index is positive then the static gain of the equivalent plant gi(s) and that of the loop diagonal plant gii(s) will be of the same sign. As a result, φai can be chosen on the basis of the sign of gii(0), instead of the sign of gi(0). 5. Simulations We now apply the method to some typical processes and compare its performance with that of the BLT
Figure 3. Step responses in Example 1 for PID design (s proposed; - - - BLT).
method4. In all of the examples below, we follow the recommendations made in the preceding section, and show our design with solid lines and the BLT design with dashed lines. The PID controller is implemented in the form where the derivative action is on the measurement signal to avoid “derivative kick”.11 Example 1. Consider the well-known Wood/Berry binary distillation column plant:12
[ ]
y1(s) ) y2(s)
[
12.8e-s 16.7s + 1 6.60e-7s 10.9s + 1
][ ] [ ]
-18.9e-3s 3.8e-8s 21.0s + 1 u1(s) + 14.9s + 1 d(s) 4.9e-3s -19.4e-3s u2(s) 14.4s + 1 13.2s + 1
Assume that a multi-loop PID controller is adopted. Because g1(0) > 0 and g2(0) < 0, we have φa1 ) 0 and φa2 ) -π. The normalized gains for the equivalent processes, κi, are evaluated from eq 30 as κ1 ) 12.8 and κ2 ) 4.08. Thus, it follows that φb1 ) 0.588, φb2 ) 0.499, rb1 ) 0.612, and rb2 ) 0.579. With the above data, the multi-loop MZN tuning procedure can be invoked: Step (1) Because PID controllers are used, Ri values are set from eq 31 as R1 ) 0.0558 and R2 ) 0.167. βi values are calculated from eq 16 as β1 ) 13.3 and β2 ) 4.56. Step (2) ω11 is determined from eq 22 as ω11 ) 1.61. Take ω2(0) ) ω22 ) 0.565 from eq 29 and the principal j 1] ) [1.48, 1.65]. Equations interval is found as [ω1,ω 20 and 19 are then solved for ω1 in [1.48,1.65] with kc2 satisfying g2(0)kc2 > 0 and C1 in eq 24 is shown in Figure 2.
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4731
Figure 4. Step responses in Example 1 for PI design (s proposed; - - - BLT).
Figure 5. Control system robustness in Example 1 for PI design (s ideal; - - - gain decrease; -‚-‚ gain increase; ‚‚‚ saturation).
Step (3) Draw C2 in eq 25 in a manner similar to that in step (2) in Figure 2. Step (4) The intersection Z of C1 and C2 in Figure 2 is determined and the solution is found from eq 27 as ω1 ) ω1* ) 1.643, ω2 ) ω2* ) 0.505, kc1 ) kc1(ω1*,ω2(0)) ) 0.949, and kc2 ) kc2(ω1*,ω2(0)) ) -0.131. To verify that further iterations will generate little improvement, we substitute ω1* for ω1(0) and ω2* for ω2(0) and run steps (2) and (3) again. The resultant curves, CiNew, are also plotted in Figure 2. From the figure, we see that CiNew and Ci are very close to each other, and the new solution found from the intersection ZNew of CiNew and C2New is ω1 ) 1.645, ω2 ) 0.504, kc1 ) 0.940, and kc2 ) -0.132, which does not differ significantly from that obtained without iteration. Therefore, a similar performance can be expected from the new solution, and any iterations are not worth trying. Step (5) Calculate TIi and TDi as TI1 ) β1/ω1 ) 8.09, TI2 ) β2/ω2 ) 9.03, TD1 ) R1TI1 ) 0.452, and TD2 ) R1TI1 ) 1.51. The multi-loop PID controller is thus formed as
gives better loop performance and shorter settling time for the decouplings, especially for the slow (second) loop. In practice, PI controllers are usually more commonly used than PID controllers. If a PI controller is adopted, one has φa1 ) 0, φa2 ) -π, φb1 ) -0.0741, φb2 ) -0.213, rb1 ) 0.362, and rb2 ) 0.369. The multi-loop MZN method is then activated to yield the following: Step (1) Because PI controller is used, Ri values are set as 0. βi values are then calculated from eq 16 as β1 ) 13.5 and β2 ) 4.64. Step (2) ω11 ) 1.61 and ω22 ) 0.565, the same as in the PID case. Take ω2(0) ) ω22 ) 1.61 from eq 29, and j 1] ) [1.58, 1.79]. the principal interval is found as [ω1,ω Equations 20 and 19 are then solved for ω1 in [1.58, 1.79] and the resulting C1 in eq 24 is drawn. Step (3) Draw C2 in eq 25 in a manner similar to that in Step (2). Step (4) The intersection Z of C1 and C2 is determined and the solution is found from eq 27 as ω1 ) ω1* ) 1.60, ω2 ) ω2* ) 0.487, kc1 ) kc1(ω1*,ω2(0)) ) 0.699, and kc2 ) kc2(ω1*,ω2(0)) ) -0.0895. Step (5) Set TI1 ) β1/ω1 ) 8.42, TI2 ) β2/ω2 ) 9.52, and TD1 ) TD2 ) 0. The multi-loop PI controller is thus obtained as
{
(
K(s) ) diag 0.949 1 +
1 + 0.452s , 8.09s
)
(
-0.131 1 +
1 + 1.51s 9.03s
)}
The step responses of the resultant feedback system to unit set-point changes followed by load disturbance change of -0.5 are shown in Figure 3 with solid lines. The corresponding step responses using the BLT method4 with kc1 ) 0.375, TI1 ) 8.29, kc2 ) -0.075, and TI2 ) 23.6 are shown with dashed lines. The proposed method
{
(
K(s) ) diag 0.699 1 +
1 1 , -0.0895 1 + 8.42s 9.52s
)
(
)}
The step responses of the resultant feedback system to the same set-point and load disturbance changes as in the PID case are shown in Figure 4 (solid lines) together with the BLT method (dashed lines). Again, it is observed that the proposed method gives better loop performance and shorter settling time for the decou-
4732 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
Figure 6. Step responses in Example 2 (s proposed; - - - BLT).
Figure 7. Step responses in Example 3 (s proposed; - - - BLT).
plings, especially for the slow (second) loop. Through the comparison of PI and PID performance with our design, one notes that for this example, PI and PID controllers give similar performance, with the PID controller offering slightly faster loop responses. To illustrate the robustness of the closed-loop system, the process is rewritten as kG(s), where G(s) is the nominal model and k is the scaler gain now perturbed. The corresponding step responses to the same set-point and load disturbance changes are shown in Figure 5, in which the responses with 20% gain decrease and 20% gain increase are shown with dashed line and dasheddotted line, respectively. Note that the control performance from the perturbed process (dashed lines) is very similar to the ideal control performance (solid lines). A real process more or less has some nonlinearities. When the nonlinearity is modest, our method can be applied directly. For instance, a saturation nonlinearity is introduced into the process G(s) in Figure 1 such that ui, i ) 1,2 are limited by ui ∈ [-0.4, 0.4]. The corresponding responses are shown in Figure 5. It is observed that u1 does run into saturation; however, the control performance from the nonlinear process (dotted lines) is very similar to that of the linear process (solid lines). Example 2. The process studied by Luyben4 has the following transfer function matrix:
Our procedure yields
[
-0.101e-12s 0.126e-6s 1 (45s + 1)(48s + 1) G(s) ) 60s +-8s -0.12e-8s 0.094e 38s + 1 35s + 1
]
{ (
K(s) ) diag 68.4 1 +
1 + 2.85s , 24.8s
)
(
-29.5 1 +
1 + 3.73s 15.8s
)}
The step responses of the resultant feedback system to unit set-point changes followed by load disturbance changes are shown in Figure 6 with solid lines. Because there is no load disturbance model for this process, load disturbance changes of -30 and 30 are applied directly on the two process inputs, u1 and u2, respectively. Here, the magnitude of the load disturbance changes is chosen as 30 to make the load response peak comparable with that caused by the unit set-point changes, to exhibit both types of responses clearly in the figure. The corresponding step responses for the BLT method are shown with dashed lines in the figure. Again, the proposed controller exhibits better loop and decoupling performance than does the BLT method. Example 3. The 24 tray tower separating methanol and water examined by Vinante and Luyben4 has the following transfer function matrix:
[
-1.3e-0.3s -2.2e-s 7s + 1 G(s) ) 7s + 1 -2.8e-1.8s 4.3e-0.35s 9.5s + 1 9.2s + 1
]
Because all of its elements are first-order in nature, we use a simple PI controller. It is designed with our
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4733
method as
Literature Cited
{
(
K(s) ) diag -1.46 1 +
1 1 , 3.40 1 + 4.30s 5.84s
)
(
)}
The step responses of the resultant feedback system to unit set-point changes followed by load disturbance changes are shown in Figure 7 with solid lines, in which load disturbance changes of 0.5 and 1 are applied directly on the two process inputs, u1 and u2, respectively. The proposed PID controller also gives better loop and decoupling performance than does the BLT method (dashed lines). One can see from the above examples that with our method, the performance enhancement over the BLT method is more noticeable in the slow loops than in the fast loops. Because slow loops dominate the system performance, their improvement is more desirable and meaningful, and our method meets this purpose. 6. Conclusions In this paper, a multi-loop version of the modified Ziegler-Nichols method was presented. The idea is to shift the critical point of each equivalent SISO plant, obtained by regarding each loop together with its corresponding interactions from all other loops, to a desired position. Such a design problem for multi-loop controllers will naturally lead to coupled nonlinear equations and a novel algorithm is presented to solve them. Simulation examples are provided to show the effectiveness of the proposed method and comparisons are made with the BLT method. The performance improvement in set-point and load disturbance responses as well as decoupling is significant with the proposed method.
(1) Astrom, K. J.; Hagglund, T. Automatic Tuning of PID Controllers; Instrument Society of America: Research Triangle Park, NC, 1988. (2) Palmor, Z. J.; Halevi, Y.; Krasney, N. Automatic tuning of decentralized PID controllers for MIMO processes. J. Process Control 1996, 42, 1174-1180. (3) Skogestad, S.; Morari, M. Robust performance of decentralized control systems by independent designs. Automatica 1989, 25, 119-125. (4) Luyben, W. L. Simple method for tuning SISO controllers in multivariable systems. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 654-660. (5) Maciehowski, J. M. Multivariable feedback design; AddisonWesley: Workingham, U.K., 1986. (6) Zhuang, M.; Atherton, D. P. Automatic tuning of optimum PID controllers. Proc. IEE, Part D 1993, 140, 216-224. (7) Ziegler, J. G.; Nichols, N. B. Optimum settings for automatic controllers. Trans. ASME 1943, 65, 433-444. (8) Hagglund, T.; Astrom, K. J. Industrial adaptive controllers based on frequency response techniques. Automatica 1991, 27, 599-609. (9) Hang, C. C.; Astrom, K. J.; Ho, W. K. Refinements on the Ziegler-Nichols tuning formula. Proc. IEE, Part D 1991, 138, 111118. (10) Astrom, K. J.; Hagglund, T. Automatic tuning of simple regulators with specifications on phase and amplitude margins. Automatica 1984, 20, 645-651. (11) Astrom, K. J.; Hugglund, T.; Hang, C. C.; Ho, W. K. Automatic tuning and adaptation for PID controllers - A survey. Control Eng. Practice 1993, 1, 699-714. (12) Wood, R. K.; Berry, M. W. Terminal composition control of a binary distillation column. Chem. Eng. Sci. 1973, 28, 17071717.
Received for review May 20, 1998 Revised manuscript received September 14, 1998 Accepted September 15, 1998 IE980313R
