ARTICLE pubs.acs.org/IECR
Extending Ideal PID Tuning Rules to the ISA Real Structure: Two Procedures and a Benchmark Campaign Martina Maggio and Alberto Leva* Dipartimento di Elettronica e Informazione, Politecnico di Milano, Via Ponzio 34/5, 20133 Milano, Italy ABSTRACT: Although virtually any installed PID controller encompasses a derivative filter, and very frequently takes the standard ISA form, many tuning rules refer to the ideal controller structure, which is subsequently made realizable by introducing the so-called “derivative filter” a posteriori. Doing so can sometimes lead to undesired results. This paper investigates the matter and discusses two very simple but rigorous procedures to extend ideal PID tuning rules to the ISA controller. A benchmark simulation campaign, referring to several well-assessed rules, demonstrates the opportunity and efficacy of the proposed extensions. A quite realistic test referring to a typical process control case is also presented.
’ INTRODUCTION Among the various implementations of the proportional integralderivative (PID) controller, the International Society of Automation (ISA) one1 is considered particularly flexible and practical. In the first place it is a real PID; i.e., it comprises the socalled “derivative filter”. Then, it permits for example complex zeros, contrary to the so-called “series” or “interacting” forms, and facilitates switching on and off one or more of the three actions, even at runtime, since the state of the two dynamic ones is confined to each of them. Moreover, when the regulator is switched from manual (or tracking) to automatic mode, the ISA forms allows for an easier (albeit necessarily arbitrary) recomputation of the derivative part state. For the above reasons, to the best of the authors’ knowledge, said form is among the most widely used in industrial codes. Nonetheless, most of the tuning rules proposed so far do not refer to that structure, or do not exploit it completely. According to an excellent review,2 a significant part of said rules (approximately 54%) sticks to the ideal PID case without any derivative filter, only about 8% refer to the ISA form or analogous ones, and the rest use either series-like forms or some specialized expression. Tuning methods for the real PID—ISA or not—are overall still a minority (less than 10% according to the quoted reference), and in most of their industrial applications the derivative filter is merely introduced a posteriori,3 or the derivative action is even just switched off.4 Indeed, although tuning the derivative filter is important,57 only a very few rules do; see, e.g., refs 810. Discussing the reasons of the scenario just sketched would be too long here. It is sufficient to suggest that tuning rules for the ideal PID, given its simpler structure, tend to be more flexible as for the tuning objectives (set point tracking, load disturbance rejection, stability margins, and so on) than rules directly conceived for the real PID. In fact, relevant issues such as the performancerobustness trade-off are seldom addressed when the real structure is considered; see again refs 9 and 10 and the papers quoted therein. In any case, rules originally conceived for the ideal PID form the backbone of most industrial autotuners, and are well understood and accepted. It is therefore interesting to extend said rules to the real PID, provided that such an r 2011 American Chemical Society
extension still allows each particular rule to achieve the specific goal it was conceived for in the ideal PID case. This paper aims at the envisaged extension and focuses on the ISA one. In fact, as will be shown, besides the mentioned advantages the ISA PID exhibits a notable keenness to produce undesired effects if the derivative filter is not tuned, backing up the “industrial myth that derivative action does not work”.11 Of course, “tuning rules” do not exhaust the variety of tools available to the skilled engineer for an efficient control design, and more refined or methodologically grounded techniques can be used by experts. However, easier and safer tuning rules are indeed a means to have also nonspecialists achieve satisfactory control results. With a view to that, removing structural flaws (such as negligence or incorrect treatment of the PID derivative term) is inherently a benefit. The paper is organized as follows. First, the relevance of the problem is shown in general, explaining what are the advantages of extending ideal PID rules instead of simply abandoning the ideal structure in the rule derivation. Then, two extension procedures are presented. Said procedures are subsequently applied to 25 tuning rules for the ideal PID, pursuing different goals, with a small but significant batch of processes as a benchmark. Aggregate results from the simulation campaign are shown and analyzed, as the main purpose of this paper, to testify to the generality and usefulness of the proposed extensions; a particular case is also commented on, to evidence that sometimes a rigorous extension methodology is mandatory to preserve the tuning quality ensured by the (ideal) tuning rules. Finally, a temperature control case is addressed involving a very detailed simulation model and noise, to show the practical applicability of the proposal. The scope is here limited to model-based rules, but of course the presented ideas are general, and could be applied, for example, to relay-based tuning. Received: October 11, 2010 Accepted: June 29, 2011 Revised: June 20, 2011 Published: June 29, 2011 9657
dx.doi.org/10.1021/ie102065y | Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
’ RELEVANCE OF THE PROBLEM The ideal or “textbook” continuous-time PID in the onedegree-of-freedom form considered herein is described by the error to control transfer function: 1 þ sτd ð1Þ RID ðsÞ ¼ k 1 þ sτi while the real ISA form is given by 1 sTd RISA ðsÞ ¼ K 1 þ þ sTi 1 þ sTd =N
ð2Þ
or equivalently, for some of the following computations, by 1 sTd RISA ðsÞ ¼ K 1 þ þ , ν ¼ 1=N sTi 1 þ sνTd ð3Þ Notice that the ideal regulator is obtained from eqs 2 and 3 with N = ∞ and ν = 0, respectively. In eqs 2 and 3, the derivative filter has two effects: one (well studied) is a radical modification of the PID high-frequency behavior; the other (less studied) is a displacement of the regulator zeros. While the former effect is essentially a benefit, as it results in a reduction of the control sensitivity to measurement noise at the cost of a generally moderate phase margin decrease, the latter may conversely alter the controller behavior in the control band, thus being potentially detrimental if not accounted for properly. To prove how relevant the zeros’ displacement problem can be, at least two reasoning paths can be followed, which are presented below. Sensitivity Analysis. The first path is based on a simple sensitivity analysis. If eq 1 is made proper by merely introducing a fixed value of ν, the numerator of eq 3, normalized for convenience with respect to the gain K, is nðs, νÞ ¼ 1 þ sðτi þ τd νÞ þ s2 ðτi τd þ τi τd νÞ
ð4Þ
while the numerator of eq 1, correspondingly normalized by a factor of k, is obviously ni ðsÞ ¼ nðs, 0Þ ¼ 1 þ sτi þ s2 τi τd
ð5Þ
Denoting by zi a zero of eq 5, the displacement of that zero induced by ν is quantified in terms of sensitivity by 0 11 ∂nðs, νÞ B∂nðs, νÞ A Sν, i ¼ @ ∂ν ∂s s ¼ zi , ν ¼ 0
sτd þ s2 τi τd ¼ ::: ¼ τi þ 2sτi τd
s ¼ zi , ν ¼ 0
ð6Þ
s ¼ zi
Apparently, eq 6 is totally general with respect to the ideal tuning rule considered, as it contains the PID parameters only, and can therefore be applied downstream the rule itself. If only the amplitude of the zeros’ displacement produced by a ν selected a priori is of concern, a synthetic sensitivity index is provided by the quantity S : ¼ maxðjSν, 1 j, jSν, 2 jÞ
first example is worth reporting right now, to see how important the addressed matter can be at least in a particular case. Example 1. Consider the first order plus dead time (FOPDT) process pðsÞ ¼
esθ 1 þ s
ð8Þ
normalized for convenience to unity gain (which is however irrelevant for the zeros’ analysis) and unity time constant, thereby allowing interpretation of θ as the normalized delay. The normalized delay is a relevant quantity in PID tuning, and in the relative literature it also termed sometimes the “controllability index”, as in several works it is taken as a measure of “how difficult a process is to control”. There is not the space here to explain the reasons for that, except for noticing that the reason basically resides in the different behavior of the process (thus the open loop) transfer function magnitude response phase, depending on whether the delay or the rational dynamics dominate. In any case, in the PID (auto)tuning domain, the fact is well established: the interested reader can refer to ref 13 for some more words on the matter, and for some detailed references. If an ideal PID is tuned for eq 8 by cancellation after replacing the term esθ with its (1,1) Pade approximation (1 sθ/2)/(1 + sθ/2), which is standard practice in the vast population of IMC-like rules, the result is τi = 1 and τd = θ/2, which substituted into eq 6 together with the expression of the zeros of eq 5 yields Sν1, 2
pffiffiffiffiffiffiffiffiffiffiffiffiffi ðθ 2Þ 1 2θ - ð3θ 2Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2θ 1 2θ
ð9Þ
As θ approaches 0.5, Sν1,2 tend to infinity, and small modifications of N can cause so large perturbations of the real PID zeros to make the corresponding closed-loop transients very different from those forecast with the ideal PID. 0 As a final remark after quitting the particular example above, by means of eq 6 one can study not only the amplitude but also the direction of the zeros’ displacement. This is surely of interest, because—quite intuitively—there will be directions that are “worse” than others as for the resulting deterioration of the closed-loop behavior. However, such a study is not within the scope of this work, since the rationale of the extensions presented later on is to have the zeros not move at all, and is thereby deferred to future research. Perturbation Analysis. In the second reasoning path, the zeros’ displacement is viewed as a perturbation of the open-loop transfer function, hence casting the problem into (quite simply) a robust control one. To this end, define for simplicity a normalized complex variable σ := sτi, and notice that denoting by ξ the (clearly positive) damping factor of the ideal PID’s zeros results in τi/τd = 4ξ2. The regulator Rz obtained from eq 1 by merely adding the derivative filter pole without accounting for the zeros’ displacement is then σ2 1 þ σ þ 2 4ξ Rz ðσÞ ¼ k σ σ 1 þ 4Nξ2
ð7Þ
that is analyzed later on, when dealing with the benchmark chosen to assess the proposed extensions’ efficacy. However, a 9658
ð10Þ
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
and the mentioned (multiplicative) perturbation is Wz ðσÞ : ¼ ¼
respectively. In addition, if qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi ̅ ξ ¼ ξ ¼ 0:5 1 þ 2≈0:7768
Rz ðσÞ RID ðσÞ RID ðσÞ σð1 þ σÞ 4Nξ2 1 þ σ þ
σ2 4ξ2
!
ð11Þ
Defining a normalized frequency ψ := ωτi, the frequency response W(jψ) has limits 0 and 1/N for ψ going to 0 and ∞,
Wz, max
the magnitude of that frequency response has a maximum at qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2jξ 8ξ2 þ 1 þ 4ξ2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ψmax ¼ ð12Þ 16ξ4 8ξ2 1 and of value
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8ξ2 þ 1ð64ξ6 þ 32ξ4 þ 4ξ2 Þ þ 256ξ6 þ 32ξ4 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8ξ2 þ 1ð1024ξ10 þ 1024ξ8 64ξ6 32ξ4 ÞN 2 þ ð256ξ6 þ 32ξ4 ÞN 2
Notice that eq 13 depends only on the PID parameters; thus also this second reasoning path is totally general with respect to the considered ideal tuning rule. Figure 1 reports the inverse (in decibels) of Wz,max (that notoriously provides a magnitude overbound for the nominal complementary sensitivity frequency response in order to preserve stability) as a function of ξ and N. As can be seen, there are cases where “not unreasonable” values of N and ξ (having a tuning rule produce complex PID zeros is not the most frequent situation but happens) make 1/Wz,max dangerously approach the 0 dB axis. Notice also (details are omitted for brevity) that ψmax ranges from 0.5 to 2 for ξ approximately in the interval (0.25,0.65). Since 1/τi is typically around the desired cutoff, which therefore occurs around ψ = 1, it is not infrequent that the minimum of 1/Wz,max falls near the nominal cutoff frequency. Most likely the effect of the zeros’ displacement cannot jeopardize stability, but according to the presented analysis, the shape of the obtained closed-loop transients can be affected to a relevant extent (an example follows later on). It is also worth observing that, by viewing the zeros’ displacement as a perturbation of the open-loop transfer function, just induced by the controller instead of the process, interesting considerations can be done on controller fragility, for example analyzing the variations of the peak sensitivity. This is however deferred to future works. Concluding Remarks. Besides showing the problem relevance, the simple analysis above gives suggestions on how to face it. In fact, any ideal-to-real tuning rule extension must not only shape the high-frequency PID behavior “correctly”, but also preserve its aspect in the control band as closely as possible
ð13Þ
(which is not guaranteed by any prespecified N). A viable and simple way to do so is to maintain the zeros’ position, which is therefore the basis of the presented methodology. From a more practical standpoint, however, one may wonder why extend tuning rules for the ideal PID, instead of just refraining from the use of that structure in deriving a rule. The main reason is that conceiving rules for eq 1 is tendentiously easier than for eq 2, particularly in the model-based case, and at least in the frequent situation where explicit tuning formulas are sought (e.g., when autotuners need realizing on process controllers that are often based on low-end hardware, with limited computational power). In such situations, if the latter regulator is used, owing to the denominator phase contribution depending on N, prescribing many relevant loop properties often requires solving some equations numerically. Along such reasoning, a viable approach is to start from the ideal PID, derive a rule by pursuing a desired objective, and then have the derivative filter “not perturb too much” that objective—an idea is backed up by the literature, according to the quoted references.
’ TWO EXTENSION PROCEDURES Assume that some tuning rule has yielded the three parameters k, τi, and τd of eq 1. The problem is to determine the four parameters K, Ti, Td, and N of eq 2. To stay abstracted from any particular tuning formula, according to the scope of this work, in the following it is only assumed that some estimated (or nominal) values for the cutoff frequency and the phase margin, termed respectively ω^c and j^m, are known. Obtaining said quantities is possible with any type of ideal PID tuning rule, generally at a moderate computational cost. The zeros of eq 1 are pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi τi - τi 2 4τi τd ideal ð14Þ s̅ 1, 2 ¼ 2τi τd while those of eq 2 are s̅ 1,ISA 2 ¼
NTi Td -
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N 2 ðTi 2 4Ti Td Þ 2NTi Td þ Td 2 2Ti Td ð1 þ NÞ
ð15Þ
Figure 1. Values of 1/Wz,max in decibels (dB) as a function of ξ and N.
Equating the two couples, based on the considerations above, yields two of the four (real) equations needed. A third equation is found by equating also the generalized gains of eqs 1 and 2, which 9659
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
provides
provides
k K ¼ τi Ti
ð16Þ
K ¼k
2
As for the fourth equation, two reasoning paths are here followed. The first one is relative to the high-frequency control sensitivity, defined as RðsÞ CðsÞ ¼ 1 þ RðsÞ PðsÞ
ð17Þ
where P(s) is the transfer function of the process under control. When the frequency ω goes to infinity, the magnitude of the frequency response of the control sensitivity function approaches the “high-frequency regulator gain” R∞ ¼ lim jRðjωÞj
ð18Þ
ωf∞
If eq 1 is taken as R(s), then R∞ is apparently ∞, while in the ° for R∞, case of 2 its value is K(1 + N). Given a desired value R∞ eq 2 can be parametrized by adding ° ¼ Kð1 þ NÞ R∞
ð19Þ
to the three equations coming from 15 and 16. Solving the soobtained system provides K ¼ k Td ¼
k2 τd ; ° τ R∞ i
Ti ¼ τi
° ° τi τd R∞ ðR∞
k τd þ kÞ ; ° 2ðτ kτ Þ R∞ i d 2
2
kτd ° R∞ N ¼
Td ¼
k
R>1
ð21Þ
The real PID eq 2 is thus parametrized by adding eq 21 instead of eq 19 to the three real equations given by 1416 above. This
2
N ¼
b c Rτi ω bc þ 1 R2 τ i τ d ω bc 1 Rτi ω
Relating the extension rule 22 to the phase margin decrease accepted for the derivative filter introduction provides a simple clue to selecting R, since that reduction (estimate) is c ¼ argðRISA ðjb ωc ÞÞ argðRideal ðjb ωc ÞÞ Δj m ¼ arctanð1=RÞ
ð23Þ
c in the range 25° are reasonable in any where values of Δj m c has an influence on the highcase of practical interest. Also Δj m frequency control sensitivity, with lower values resulting in a looser bound for that quantity.
’ A BENCHMARK SIMULATION CAMPAIGN To test the proposed extension procedures, a benchmark set of four processes is here employed, each one characterized by a single parameter p and exhibiting different dynamics. The set is
° ° τ d þ τ i R∞ ðR∞ kÞ ° τ kτ Þ kðR∞ i d
To apply the extension rule 20, it is necessary to select a ° . If some information is available on the suitable value for R∞ measurement noise level and the admissible actuator stress, said choice can be made from technological facts. More interesting for ° needs to this work is the opposite case, where the selection of R∞ come from the tuning model, and possibly the design parameters used in the rule that synthesized eq 1. In this case, a reasonable idea is to constrain the PID high-frequency gain below a chosen ° = kR∞kmf. multiple of the midfrequency value kmf, i.e., setting R∞ A suitable value for kR∞ can be chosen in the range 330, while for simplicity kmf can be estimated as k, with the approximation being very good for ideal PIDs with real distinct zeros, and still decent for PIDs with complex zeros as long as the damping factor remains in the range provided by the considered rules in nonpathological cases (i.e., when the process is “suitable for PID control” according, e.g., to works like ref 1). Of course, lower values of kR∞ result in a tighter bound on the high-frequency control sensitivity. The second reasoning path followed herein refers to the phase margin reduction caused by the derivative filter’s pole. If the tuning rule for eq 1 provides some estimate ω^c of the closed-loop cutoff frequency, it is possible to fix the position of the angular frequency N/Td of the second pole of eq 2 with respect to ω^c, for example by imposing that
bc 1 Rτi ω ; Rb ωc
b c Rτi ω bc þ 1 R2 τ i τ d ω ; b c 1Þ Rb ωc ðRτi ω
P1 ðsÞ ¼ 2
Ti ¼
ð22Þ
1 ð1 þ sÞð1 þ psÞð1 þ p2 sÞð1 þ p3 sÞ
p ∈ ½0:1; 1
1 ps p ∈ ½0:1; 2:5 ð1 þ sÞ3 1 þ ps p ∈ ½0:05; 1:8 P3 ðsÞ ¼ ð1 þ sÞð1 þ 0:1sÞ 1 P4 ðsÞ ¼ p ∈ ½0:3; 0:9 1 þ 2ps þ s2 P2 ðsÞ ¼
ð20Þ
N ¼ Rb ωc , Td
bc 1 Rτi ω ; bc Rτi ω
ð24Þ
and comprises high-order, undershooting, fastslow or overshooting, and oscillatory dynamics. For the sake of synthesis a smaller benchmark set is used with respect to other works such as ref 14; nonetheless the variety of dynamics contained in eqs 24 is sufficient to support the claims of this paper. The absence of models with delay is motivated in the first place by the intention of evidencing that the problems addressed here are relevant even in the case of rational dynamics only. In addition, especially in process control, systems that really contain a delay are a minority, and are well identifiable. More frequent, thus more interesting here, is the case of high-order (and possibly nonminimumphase) rational dynamics. An example with a delay process is here reported to avoid leaving the matter totally untouched, although a complete treatise of it needs deferring to future works. Based on the equation set 24, the extension procedures just presented are here applied to 25 ideal PID tuning methods, listed in Table 1. Said methods pursue different goals, privileging, e.g., set point tracking versus disturbance rejection, and share the FOPDT model structure PðsÞ ¼ μ
esD 1 þ sT
ð25Þ
that in all the tests is parametrized based on process step responses with the well-known method of areas. The choice of the used methods in the enormously vast panorama of available ones2 was made by selecting those most frequently 9660
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
Table 1. Ideal PID Tuning Methods Considered in This Work no.
authors
τi
k
ref
τd
notes/remarks
1
Ziegler and Nichols
15
aT/μD
2D
0.4D
reduced τd, a = 1.2
2
Ziegler and Nichols
15
aT/μD
2D
0.4D
reduced τd, a = 2
3
Åstr€om and H€agglund
1
0.94T/μD
2D
0.4D
reduced τd
4
Chien et al.
16
0.95T/μD
2.38D
0.42D
regulatory, 0% overshoot
5
Chien et al.
16
1.2T/μD
2D
0.42D
regulatory, 20% overshoot
6
Chien et al.
16
0.6T/μD
D
D/2
servo, 0% overshoot
7
Chien et al.
16
0.65T/μD
1.36D
0.47D
servo, 20% overshoot
8
Cohen and Coon
17
1:35T=D þ 0:25 μ
Tð2:5D=T þ 0:46ðD=TÞ2 Þ 1 þ 0:61D=T 0.75
0:37D 1 þ 0:2D=T
9
Murrill
18
0.92
(1.44/μ)(T/D)
0.95
0.48T(D/T)1.14
regulatory, min IAE
(T/1.1)(T/D)
0.56T(D/T)1.01
regulatory, min ISE
(T/0.88)(T/D)
0.77
10
Murrill
18
(1.50/mu)(T/D)
11
Murrill
18
(1.36/mu)(T/D)0.95
(T/0.84)(T/D)0.74
0.38D
regulatory, min ITAE
12
Zhuang and Atherton
19
(1.47/μ)(T/D)0.97
(T/1.12)(T/D)0.75
0.55T(D/T)0.95
regulatory, min ISE, 0.1 e D/T < 1.05
19
(1.52/μ)(T/D)0.74 (1.47/μ)(T/D)0.97
(T/1.13)(T/D)0.64 (T/0.94)(T/D)0.73
0.55T(D/T)0.85 0.44T(D/T)0.94
regulatory, min ISE, 1.05 e D/T e 2 regulatory, min ISTSE, 0.1 eD/T e 1.05
(1.52/μ)(T/D)0.73
(T/0.96)(T/D)0.6
13 14
Zhuang and Atherton Zhuang and Atherton
19
0.44T(D/T)0.85
regulatory, min ISTSE, 1.05 e D/T e 2
(1.53/μ)(T/D)
(T/0.97)(T/D)
0.75
0.41T(D/T)0.93
regulatory, min ISTES, 0.1 e D/T < 1.05
(1.59/μ)(T/D)0.71
(T/0.96)(T/D)0.0.6
0.41T(D/T)0.85
regulatory, min ISTES, 1.05 e D/T e 2
T 0:74 0:13D=T T 0:8 0:15D=T T 1:20 0:37D=T T 1:05 0:22D=T T 0:99 0:24D=T T 0:92 0:17D=T T 0:98 0:25D=T T 0:89 0:17D=T
0.35T(D/T)0.91
servo, min IAE
0.31T(D/T)0.93
servo, min ITAE
0.49T(D/T)0.99
servo, min ISE, 0.1eD/T < 1.05
0.49T(D/T)0.71 0.39T(D/T)0.99
servo, min ISE, 1.05 e D/T e 2 servo, min ISTSE, 0.1 e D/T < 1.05
0.38T(D/T)0.71
servo, min ISTSE, 1.05 e D/T e 2
0.32T(D/T)0.89
servo, min ISTES, 0.01 e D/T < 1.05
0.32T(D/T)0.83
servo, min ISTES, 1.05 e D/T e 2
T + D/2
TD=2 T þ D=2
λ = max(0.25D,0.2T)
T+D/2
TD=2 T þ D=2
λ=0
0.96
15
Rovira et al.
20
(1.09/μ)(T/D)0.87
16
Rovira et al.
20
(0.97/μ)(T/D)0.85
17
Zhuang and Atherton
19
(1.05/μ)(T/D)0.9
19
(1.15/μ)(T/D)0.57 (1.04/μ)(T/D)0.91
18
Zhuang and Atherton
(1.14/μ)(T/D)0.84 19
Zhuang and Atherton
19
0.9
(0.97/μ)(T/D)
(1.06/μ)(T/D)0.58 20
Dahlin
21
T þ D=2 μðD þ λÞ
21
Dahlin
21
T þ D=2 μðD þ λÞ
22
ðT þ D=2Þð0:7303 þ 0:5307T=DÞ μðT þ DÞ
22 23 24 25
Wang, Juang, and Chan Lee et al. Lee et al. Åstr€om and H€agglund
23
T þ
D2 2ðD þ λÞ D2 2ðD þ λÞ
23
T þ
12
0:2 þ 0:45T=D μ
0:5TD T þ D=2
T+D/2
D2 þ
τi μðD þ λÞ
τi D3 þ D þ λ 3ðD þ λÞ 2λi
IMC, λ = max(0.25D,0.2T)
τi μðD þ λÞ
τi D3 D2 þ þ D þ λ 3ðD þ λÞ 2λi
IMC, λ = 0
þ 0:8T D0:4D D þ 0:1T
TD=2 0:3D þ T
AMIGO
encountered in the applications, of course to the best of the authors’ experience; extending the set to other rules is apparently straightforward. For the sake of clarity, it is worth noticing that the aim of the presented benchmark campaign is not to compare the tuning methods, that are in fact so heterogeneous as to make such a comparison hard to formalize. On the contrary, the goal is to show that the proposed extensions work satisfactorily (and better than “fixing N”) on so wide and various a range of methods as to allow considering them a general tool. With reference to Table 1, just two notes are worth adding. For methods 13, with respect to the original rules, the derivative time τd is here reduced by a factor of 0.8. Such reductions are frequently encountered in industrial versions of those and similar rules, generally resulting in less nervous (thus preferable) control behavior. For methods 21 and 24, where λ is interpreted as the desired closed-loop dominant time constant, setting λ = 0 is not generally advised, although in several cases low values of that parameter, in comparison to the widely used, standard choice λ = max(0.25D,0.2T), can be introduced to avoid sluggish loop behaviors. In this work λ = 0 has the meaning of a “very
aggressive” tuning, opposite to the definitely “cautious” one yielded by the standard choice mentioned above. Also, only one-degree-of-freedom (1-dof) tuning rules are considered. There are two reasons for that. First, but less important, they are the majority. Second, it is easy to see that a two-degree-of-freedom (2-dof) PID can be represented as a feedback 1-dof PID part plus a unity-gain set point prefilter, both affected by the zeros’ displacement problem: thus, dealing with the former dictates also how to treat the latter. Finally, to evaluate the effects of the proposed extensions, both set-point- and loaddisturbance-related indices are used for uniformity with the way similar results are presented in the literature. Notice, however, that the conclusions drawn later on are backed up also by the sole load disturbance responses, where the 1- or 2-dof nature of the PID is irrelevant, and that is often the preferred means to evaluate an autotuner—see, e.g., ref 24—avoiding the “overemphasis on the set point response” evidenced in works such as ref 25. Sensitivity of the PID Zeroes to N. In Figure 2, the quantity S defined by eq 7 is plotted against parameter p for the four processes 24 and the 25 methods considered; see Table 1. 9661
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
Figure 2. Quantification of the PID zeros’ sensitivity to a ν selected a priori for the benchmark processes of eqs 24, and the tuning methods of Table 1.
To better evidence the potentially more critical cases given by ideal PIDs with complex zeros, the curves evidence where the damping factor of√those zeros is less than 1 (thick blue line segments) and of 2 (double thickness red segments), excluding on the other hand the cases where the phase margin computed with the tuning model 25 is smaller than 10°. As can be seen, the values of S span more than 1 order of magnitude and can be definitely high, whence the general relevance of the displacement phenomenon. Note also how different the behavior of S is with the various methods, and in some cases the presence of small p intervals provides sharp peaks. Specific Case Study. One may question the importance of preserving the PID zeros’ position, for example objecting that the really significant quantities are tuning quality indicators such as the cutoff frequency and the phase margin, on which the zeros’ position seems to have a limited impact. In principle such an objection can make sense, but the problem here is slightly different. The question is not what is the effect of displacing some zeros in a regulator, but what is the effect of displacing said zeros while at the same time adding a pole. In fact, the figures of Table 2 later on do provide an answer, but to make things more evident, before showing aggregate results, in this example an interesting case is taken from the benchmark simulation campaign of the next example, and briefly commented on. The case is relative to method 25 (the AMIGO), and process type 2. Figure 3 shows the closed-loop response of the controlled variable to a load disturbance unit step, with parameter p taking the values 0.1, 0.5, 1.0, and 1.2. The red plot is obtained with the ideal PID, the blue one is obtained with the real PID yielded by ^ = 2°, the magenta one the second proposed procedure and Δj m is obtained with the value of N computed by the same procedure but with the ideal PID parameters, and the green one is obtained with the ideal PID made real with N fixed at 10. As can be seen, the computed N applied to the ideal PID, i.e., not preserving the zeros’ position, best approximates the time domain behavior of
the ideal one. However, with that PID, as p goes from 0.1 to 1.2 the high-frequency control sensitivity increases from 48.18 to 160.63, while with the second procedure the same quantity decreases from 36.52 to 17.94. The extension with the fixed (and quite standard) N gives even lower high-frequency control sensitivities (from 19.69 down to 11.15), but the transients’ degradation is apparent. In synthesis, the proposed procedure slightly degrades the time domain transient, but the effect only appears for high values of p and is comparable to that introduced by the computed N applied to the real PID. On the other hand, the proposed extension—thanks to the zeros’ position preservation—gives far lower high-frequency control sensitivities, and is definitely safer than the fixed-N one. It is worth noticing that several other cases similar to the one just shown were found in the simulation campaign results, with no evident distribution pattern. This corroborates the idea that the zeros’ sensitivity irregularities of Figure 2 may unexpectedly result in significant differences between the control results as forecast with the tuned ideal PID and those obtained with the real ISA one, unless a rigorous extension procedure (like those proposed here) is adopted. Aggregate Results. This section reports some very synthetic data from the full simulation campaign conducted with each process in eqs 24, for parameter p spanning the indicated ranges, and with each of the 25 tuning rules of Table 1. Table 2 c used in summarizes the results. The values of kR∞ and Δj m the campaign were respectively 20 and 4°—some further comments on the matter are given later on. As anticipated, each rule has a particular purpose. Nonetheless, in coherence with the variety of situations that an autotuner has to withstand irrespectively of which rule is used, all rules were applied to all processes and evaluated with the same indices (the ratio of Integral of the Squared Error (ISE), Integral of the Squared Time multiplied Error (ISTE), and the Integral of Time multiplied Absolute Error (ITAE) obtained with the real PID over the same quantities computed with the ideal one) for both the set point tracking (sp) 9662
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
Table 2. Aggregate Results from the Benchmark Simulation Campaign P1
R∞
c Δj m
c , id Δj m
NLO
NMI
NHI
N
P2
P3
P4
min
max
min
max
min
max
min
max
19.00829
19.66745
19.14573
19.51282
19.00522
19.20202
19.00288
19.51282
ISE sp
0.99997
1.00127
1.00002
1.00113
0.99993
1.00000
0.99989
1.00002
ISTE sp
0.99250
0.99993
0.99401
1.01194
0.98808
0.99985
0.99173
0.99991
ITAE sp
0.98091
1.01171
0.98603
1.02503
0.98435
1.00223
0.98408
1.00010
ISE ld
0.97947
0.99998
0.98642
1.01283
0.99467
0.99999
0.97978
0.99999
ISTE ld
0.93790
1.01714
0.95864
1.03240
0.99481
1.00018
0.97076
1.00045
ITAE ld
0.95513
1.00194
0.97830
1.01277
0.99824
1.00016
0.97756
1.00000
N ISE sp
0.12669 0.99992
9.23741 1.00279
0.14724 0.99838
9.47323 1.02240
0.16427 0.99966
8.14407 1.00000
0.06092 0.99967
14.74584 1.00215
ISTE sp
0.98560
0.99966
0.98795
1.27584
0.95415
0.99944
0.95131
1.02917
ITAE sp
0.95754
1.03440
0.96585
1.94914
0.95325
1.00602
0.91919
1.24490
ISE ld
0.95847
0.99995
0.97170
1.44732
0.97974
0.99998
0.88215
1.07225
ISTE ld
0.87256
1.04472
0.90290
2.78911
0.97823
1.00065
0.89123
1.22025
ITAE ld
0.90559
1.01107
0.95404
1.91413
0.99318
1.00076
0.87021
1.13297
N
0.05838
60.87213
7.03237
840.03044
0.03806
0.33867
0.15080
3.92611
ISE sp ISTE sp
1.00003 0.97497
1.00257 1.03582
1.00001 0.99157
1.00061 1.00691
0.99989 0.70318
1.00068 0.99321
0.99986 0.97539
1.00839 1.13481
ITAE sp
0.95540
2.10662
0.98160
1.01428
0.86265
1.00000
0.93079
6.76373
ISE ld
0.92682
0.99962
0.98075
1.00734
0.73648
0.99802
0.94939
1.19759
ISTE ld
0.88776
2.48608
0.93767
1.01826
0.71795
1.00397
0.78191
4.74196
ITAE ld
0.81537
1.00117
0.96753
1.00684
0.79180
1.00682
0.90135
1.01772
N
4.00000
4.00000
4.00000
4.00000
4.00000
4.00000
4.00000
4.00000
ISE sp
0.99992
1.00702
1.00009
1.00525
0.99966
1.00000
0.99966
1.00014
ISTE sp ITAE sp
0.97771 0.93528
0.99968 1.13722
0.97809 0.96751
1.06245 1.14867
0.94782 0.92385
0.99931 1.00729
0.96631 0.93953
0.99963 1.00780
ISE ld
0.92397
0.99992
0.94358
1.07444
0.97562
0.99997
0.91827
0.99996
ISTE ld
0.74727
1.25660
0.84191
1.22252
0.97520
1.00079
0.86927
1.00297
ITAE ld
0.80818
1.00993
0.90331
1.09576
0.99215
1.00061
0.90996
1.00000
10.00000
10.00000
10.00000
10.00000
10.00000
10.00000
10.00000
10.00000
ISE sp
0.99995
1.00256
1.00003
1.00202
0.99986
1.00000
0.99981
1.00003
ISTE sp
0.98675
0.99987
0.98906
1.02342
0.97763
0.99972
0.98478
0.99983
ITAE sp ISE ld
0.96547 0.96207
1.03044 0.99997
0.97769 0.97448
1.05063 1.02581
0.96993 0.98989
1.00383 0.99999
0.97116 0.96287
1.00035 0.99998
ISTE ld
0.88314
1.03826
0.92450
1.06935
0.99001
1.00035
0.94434
1.00093
ITAE ld
0.91460
1.00381
0.95877
1.02766
0.99688
1.00029
0.95869
1.00000
50.00000
50.00000
50.00000
50.00000
50.00000
50.00000
50.00000
50.00000
ISE sp
0.99999
1.00048
1.00001
1.00052
0.99997
1.00000
0.99996
1.00001
ISTE sp
0.99686
0.99997
0.99757
1.00451
0.99536
0.99994
0.99676
0.99996
ITAE sp
0.99230
1.00357
0.99393
1.00926
0.99406
1.00092
0.99370
1.00001
ISE ld ISTE ld
0.99162 0.97497
0.99999 1.00611
0.99459 0.98276
1.00477 1.01173
0.99794 0.99801
1.00000 1.00008
0.99206 0.98868
1.00000 1.00016
ITAE ld
0.98190
1.00073
0.99140
1.00431
0.99933
1.00006
0.99130
1.00000
N
N
and the load disturbance rejection (ld) results. Also, the best and worst performances are reported for each process (P14) and extension type: from top to bottom the two proposed procedures, the extension obtained by computing N according to the second procedure but not accounting for the zeros’ displacement (i.e., c -based extensions but leaving K, Ti, computing N with the Δj m and Td equal respectively to k, τi and τd), and finally the use of three (low, midrange, and high) fixed values of N. Clearly a wide distance between the best and the worst indices may denote an unsatisfactory extension on almost the totality of the cases, but also just one case
yielding particularly bad results: the rationale here is to flag both situations as undesired, since in real-life application autotuners typically do not allow selectiion of the tuning rule, irrespective of the aptitude of said rule to attain the particular tuning goal at hand. Of course, all the cases were excluded in which a rule would be applied outside the applicability limits specified by its authors, and also in which the phase margin computed with the FOPDT model was lower than a threshold set to 35° owing to the necessity of autotuners to maintain a relatively high stability degree as evaluated with the model, which may easily be an underestimation of the really achieved one. 9663
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
Figure 3. Particular case from the benchmark processes of eqs 24, and the tuning methods of Table 1.
Figure 4. Delay case drawn from ref 5.
On the basis of Table 2, the following can be stated. First, with the two proposed extensions, the indices’ degradation is comparable to that obtained with medium or high values of N, although said extensions produce a high N “only when needed”, and in general yields lower values of that parameter. Second, and quite surprising, without the previous considerations of the problem relevance, c -based procedure but disregarding computing N with the Δj m the zeros’ displacement produces on average the worst results. More in general, when a rule is well suited for the tuning problem at hand, the proposed extensions gives comparable or better real PIDs, with typically lower values of the highfrequency control sensitivity. On the other hand, when the considered rule is not well suited for the tuning problem, the proposed extensions tend to mitigate (or at least not to enhance) the possible additional degradation caused by the derivative filter, while the same is not true for the other extensions (except again the “high N” one, but recall that with the proposed extensions lower values of N tend to arise). Notice that in the campaign kR∞ was selected to give a moderately tight bound on the control sensitivity, while the c is looser. Other values corresponding bound given by Δj m were tested, allowing the conclusion that the high-frequency control sensitivity can be governed within a wide enough range, invariantly yielding advantages over nonrigorous extension approaches.
Figure 5. Temperature control plant.
In synthesis, applying the proposed procedures with possibly an upper limit for N (say from 30 to 50, the value taken for NHI in this example) appears to improve the tuning results when said rules need using with the ISA PID. Example with a Delay Process. The case analyzed in this section is drawn from ref 5 and has two purposes. The first is to show the advantages of the proposed extensions in a delaydominated case. Of course this does not exhaust the matter, as the Pade-based design approach used in that paper is just one of many alternatives, but it nonetheless refers to a wellfunctioning method and is thereby a representative case. The second (and principal) purpose is to further evidence the zeros’ displacement relevance. 9664
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research
ARTICLE
Table 3. Regulator Parameters in the Temperature Control Test
AMIGO Wang, Juang, and Chan Chien, 0% overshoot, regulatory Rovira, minimum IAE, servo
k
τi
τd
K
Ti
Td
N
kR∞ = 3 c = 4° Δj
0.76
311.73
34.10
0.73
300.36
24.03
2.11
0.93
835.00
33.53
0.90
807.03
6.72
0.24
kR∞ = 3 c = 4° Δj
1.55
166.60
29.40
1.46
156.80
21.44
2.19
1.29
1097.96
30.51
1.26
1069.98
3.33
0.12
m
m
Figure 6. Temperature control results with the Chien (left) and the Rovira (right) rules. The top row shows set point and controlled variable and the bottom row reports the control signal; see Table 3.
In ref 5, an ISA PID with N fixed at 10 (PID2) is compared to an output-filtered PID (PID3) having the derivative filter time constant as an additional parameter, noticing that in some cases the latter produces a smoother control behavior (less step-like changes in the control variable, for example). The test reported in ref 26 is here repeated. The process is an FOPDT with unity gain and time constant, and a delay of 5. The recommended PID setting is implemented as PID2 and PID3, and also as an ISA PID obtained with the R∞-based extension proposed here, and kR∞ = 3. Figure 4 reports the results, showing that the differences between PID2 and PID3 are basically due to the zeros’ displacement negligence: the ISA PID extended as proposed here better approximates the ideal one than PID3 does, also without incurring high-frequency sensitivity problems.
second with the proposed extensions. Both were used: with the kR∞-based one that parameter was set to 3 as an operator would be trained to do in the presence of a noisy measurement, while c -based one the parameter was set to 4°, thus with the Δj m allowing for a certain phase margin reduction in a view to smooth the high-frequency control behavior. Parameters are summarized in Table 3. Figure 6 shows some closed-loop transients in response to a 5 °C set point step at t = 1000, and a disturbance consisting in a drop of pHin from 2 to 1.3 bar at t = 6000, while pHout is the atmospheric pressure. Only the (regulatory) Chien and the (servo) Rovira rules are shown for brevity. Observe again the advantages in terms of control sensitivity, here obtained in a realistic situation, with a complex process dynamics, and definitely a coarse identification of the used model.
’ A REALISTIC TEST To show how the proposed extensions would be used in practice, the case is briefly addressed of temperature control in a vessel. The plant scheme is illustrated in Figure 5: a heating fluid is taken from a line at pressure pHin, with its flow being modulated by a valve operated by the PID TC, and heats a vessel through a coil, from which it is discharged at pressure pHout. A simulator was realized with the object-oriented language Modelica,27,28 using accurate, nonlinear, first-principles differential and algebraic equations along the approach explained in ref 29. The plant model alone is composed of 1754 equations, with 255 state variables. To further improve realism, noise is introduced in the transmitter TT that measures the controlled temperature T. In accordance with the “practical” character of this example, a step test was simulated on the valve command x, and the recorded data were used to parametrize an FOPDT model by hand, with a ruler and the tangent method, obtaining μ = 7, T = 800 s, and D = 70 s. Four rules from Table 1 were then applied, and for each rule two real PIDs were generated, with the first keeping the ideal parameters and taking the standard choice N = 10 and the
’ CONCLUSIONS AND FUTURE WORK A methodology was proposed to extend tuning rules for the ideal PID to the real ISA form, to take advantage of the greater flexibility of ideal rules while synthesizing a realistic controller structure. The proposed methodologies allows the derivation of extremely simple extension procedures, in favor of industrial realizability. Two extension procedures were tested on a benchmark set of processes and 25 well-established tuning rules. Based on such a wide-spectrum campaign, the yielded advantages are quite evident. Also, some tests on realistic cases (one of which was here reported) indicate the practical applicability of the proposed procedures. Future work will concentrate on several aspects: widening the set of rules considered (including, e.g., nonmodel-based ones), and the set of processes (e.g., with delay dominated ones), further investigating the role of the model parametrization procedure, or of the tuning experiment at large, refining the parameter selection mechanism, that albeit not strictly necessary is surely helpful for the operator, and finally 9665
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
Industrial & Engineering Chemistry Research implementing the proposed procedures in real industrial regulators.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
ARTICLE
(27) Mattsson, S.; Elmqvist, H.; Otter, M. Physical system modeling with Modelica. Control Eng. Pract. 1998, 6, 501–510. (28) The Modelica Association, Modelica home page. http://www. modelica.org, 19972010. (29) Casella, F.; Leva, A. Modelling of thermo-hydraulic power generation processes using Modelica. Math. Comput. Modell. Dyn. Syst. 2006, 12, 19–37.
’ REFERENCES (1) Åstr€om, K.; H€agglund, T. PID Controllers: Theory, Design and Tuning, 2nd ed.; Instrument Society of America: Research Triangle Park, NC, USA, 1995. (2) O’Dwyer, A. Handbook of PI and PID Controller Tuning Rules; Imperial College Press: London, 2006. (3) Li, Y.; Ang, K.; Chong, C. Patents, software, and hardware for PID control—an overview and analysis of the current art. IEEE Control Syst. Mag. 2006, No. Feb, 42–54. (4) Ang, K.; Chong, G.; Li, Y. PID control system analysis, design, and technology. IEEE Trans. Control Syst. Technol. 2005, 13, 559–576. (5) Luyben, W. Effect of derivative algorithm and tuning selection on the PID control of dead-time processes. Ind. Eng. Chem. Res. 2001, 40, 3605–3611. (6) Isaksson, A.; Graebe, S. Derivative filter is an integral part of PID design. IEE Proc.: Control Theory Appl. 2002, 149, 41–45. (7) Cominos, P.; Munro, N. PID controllers: recent tuning methods and design to specification. IEE Proc.: Control Theory Appl. 2002, 149, 46–53. (8) Leva, A.; Colombo, A. On the IMC-based synthesis of the feedback block of ISA-PID regulators. Trans. Inst. Meas. Control 2004, 26, 417–440. (9) Vilanova, R. IMC based Robust PID design: tuning guidelines and automatic tuning. J. Process Control 2008, 18, 61–70. (10) Alcantara, S.; Pedret, C.; Vilanova, R. On the model matching approach to PID design: Analytical perspective for robust Servo/ Regulator tradeoff tuning. J. Process Control 2010, 20, 596–608. (11) Reference 10, p 602. (12) Åstr€om, K.; H€agglund, T. Revisiting the Ziegler-Nichols step response method for PID control. J. Process Control 2004, 74, 635–650. (13) Reference 12, p 638. (14) Åstr€om, K.; H€agglund, T. Benchmark systems for PID control. Presented at the IFAC Workshop on Digital Control—Past, Present, and Future of PID Control, Terrassa, Spain, 2000. (15) Ziegler, J.; Nichols, N. Optimum Settings for Automatic Controllers. Trans. ASME 1942, 64, 759–768. (16) Chien, K.; Hrones, J.; Reswick, J. On the automatic control of generalised passive systems. Trans. ASME 1952, 74, 175–185. (17) Cohen, G.; Coon, G. Theoretical consideration of retarded control. Trans. ASME 1953, 75, 827–834. (18) Murrill, P. Automatic Control of Processes; International Textbook Co.: Scranton, PA, USA, 1967. (19) Zhuang, M.; Atherton, D. Automatic tuning of optimum PID controllers. IEE Proc.: Control Theory Appl. 1993, 140, 216–224. (20) Rovira, A.; Murrill, P.; Smith, C. Tuning controllers for setpoint changes. Instrum. Control Syst. 1969, 42, 67–69. (21) Dahlin, E. Designing and tuning digital controllers. Instrum. Control Syst. 1968, 41, 77–81. (22) Wang, F.; Juang, W.; Chan, C. Optimal tuning of cascade PID control systems. Proceedings of the Second IEEE Conference on Control Applications, Vancouver, Canada; IEEE: New York, 1993; pp 825828. (23) Lee, Y.; Park, S.; Lee, M.; Brosilow, C. PID controller tuning for desired closed-loop responses for SI/SO Systems. AIChE J. 1998, 44, 106–115. (24) Leva, A. Autotuning process controller with improved load disturbance rejection. J. Process Control 2005, 55, 223–234. (25) Shinskey, F. Process control: as taught versus as practiced. Ind. Eng. Chem. Res. 2002, 41, 3745–3750. (26) Reference 5, Figure 5. 9666
dx.doi.org/10.1021/ie102065y |Ind. Eng. Chem. Res. 2011, 50, 9657–9666
