Quantification of Valve Stiction and Dead Band in Control Loops

Article pubs.acs.org/IECR

Quantification of Valve Stiction and Dead Band in Control Loops Based on the Harmonic Balance Method Alancardek P. Araujo,*,† Celso J. Munaro,*,‡ and Moacir Rosado Filho*,‡ †

Department of Mathematics and ‡Department of Electrical Engineering, Federal University of Espiŕito Santo, Brazil ABSTRACT: The presence of stiction in control valves often causes oscillations in control loops, with negative effects on quality and cost of goods. To address this issue, it is necessary to quantify this stiction to decide about maintenance or to implement compensators that can improve control loop performance until the next plant stop. The describing function (DF) method is a well-known scheme to predict the period and amplitude of limit cycles in control loops, requiring the knowledge of linear and nonlinear parameters of the system model. In the present study, a method has been proposed to estimate these nonlinear parameters using the parameters of linear model and amplitude and period of limit cycle produced by nonlinearity. A procedure is proposed to overcome the case of unknown process model. In addition, generalization of this method to the case of parametric uncertainties in the linear part of the control loop has also been presented. The result is a simple and efficient algorithm that can be easily extended to other nonlinearities. Furthermore, the conditions for existence and uniqueness of solution for dead band and stiction estimations have been obtained. The usefulness of the proposed method has been demonstrated through simulations and its applications to a pilot plant and to real industrial data.

■

INTRODUCTION Control loops containing valves often oscillate due to valve stiction. The presence of oscillations increases fluctuation of process variables, decreasing production quality. Valve stiction is well-known as one of the main factors that contributes to destabilization and poor performance of control loops. There are several methods in the literature for detection and quantification of stiction in control valves. In an earlier study,1 a method for quantifying stiction based on the identification of a Hammerstein system has been presented. The nonlinear part of this Hammerstein system consists of a simplified one-parameter model for the stiction. The linear part consists of a process modeled by a transfer function. In this methodology, the stiction model parameter is estimated together with the coefficients of the process transfer function. A method to detect and quantify stiction using only routine operation data of the control loop has also been presented1,2 The presence of nonlinearities is confirmed through the bicoherence test, and, in affirmative cases, stiction is estimated by fitting an ellipse to a filtered controller output (OP) × controlled variable (PV) plot. The most immediate limitation of this method lies in the fact that the shape of the OP × PV plot depends on the controller parameters, delay, and time constant of the process, etc. In another study,3 a method for estimating the parameters of a two-parameter stiction model has been presented. It involves searching for model parameters (S,J) in a grid of values. For each pair (S,J) of the grid, a process is identified and the output signal is estimated. The estimated parameters and process are those that produce the minimum mean squared error between measured and estimated output. Another method of estimating the parameters (S,J) of a two-parameter stiction model based on global optimization techniques has also been presented.4 First, an initial estimate for the parameter S of stiction is calculated using the ellipse method2 and then, a search for an optimum point near the initial estimate is performed, using genetic algorithms or a © 2012 American Chemical Society

standard search algorithm. If the initial estimate is not good, then the obtained optimal point may not be an adequate estimate. The search technique based on genetic algorithms requires a great computational effort. Furthermore a method for estimating dead band in stable process control loops has also been proposed.5 The dead band estimation is obtained as a function of small variations of control and output signals around the reference signal. Another study6 proposed a new method for detection and quantification of nonlinearities, emphasizing the dead band, based on the describing function (DF). The present study has extended these concepts to the case of valves with stiction which are represented by the two-parameter stiction model.7 The model considers the control loop shown in Figure 1, in which C

Figure 1. Control loop with sticky valve represented by the nonlinearity 5.

represents a PID controller, 5 is the nonlinearity that represents the valve stiction, P is a linear process, yr is the reference signal, u is the control signal and y is the output signal. Once it has been confirmed that the oscillation of u is due to stiction, it is important to quantify it. Here, stiction quantification is performed by estimating the parameters S and J of the twoparameter stiction model.7 Choudhury et al.3 and Jelali4 used this same stiction model, and estimated the parameters S and J by the methods cited above. The present study proposes a new Received: Revised: Accepted: Published: 14121

October October October October

12, 2011 5, 2012 5, 2012 6, 2012

dx.doi.org/10.1021/ie2023417 | Ind. Eng. Chem. Res. 2012, 51, 14121−14134

Industrial & Engineering Chemistry Research

Article

Describing Function Method. The DF method is an extension of the classical method of frequency domain analysis for nonlinear systems such as those represented by the block diagram in Figure 1. This method associates with the nonlinearity 5 a complex function denoted by N, which is as follows: let u = A sin (ωt) be a sinusoidal signal in the input of the nonlinearity. The nonlinearity output u5 in general is not sinusoidal, but it can be considered the fundamental component of the decomposition of uN in the Fourier series, M sin (ωt + ϕ). If Aejωt and Mej(ωt+ϕ) denote the complex representations of the sine wave input and output fundamental component of 5 , respectively, then by definition, the DF of the nonlinearity 5 is given by

methodology for estimating the parameters S and J of stiction which is based on the DF method. The methodology requires the knowledge of the linear part loop, i.e., the controller C and process P. Although only estimates of S are used for compensation purposes, when both S and J are computed best estimates are produced for S, since both are used to explain the stiction behavior.

■

STICTION, DEAD BAND, AND THE DESCRIBING FUNCTION METHOD Stiction and Dead Band in Control Valves. The definition of stiction in control valves used in this article conforms to the definition presented in an earlier study,8 which is based on the input−output behavior of a real valve with stiction. More precisely, stiction is a property of an element such that its smooth movement in response to a variable input is preceded by a static part (dead band+stickband), denoted by S, followed by an abrupt jump, called slip jump, denoted by J. The parameter S is expressed as a percentage of the input signal, while the parameter J is expressed as a percentage of the output signal. Figure 2 shows the input versus output graph of the control valve with stiction, emphasizing the parameters S and J.

N (A , ω) =

Me j(ωt + ϕ) Ae jωt

(1)

The DF method is primarily used in limit-cycles prediction in control systems with nonlinearities, as in the block diagram of Figure 1. To this end, the control system of Figure 1 is considered to have zero reference, and the output y is a periodic function. Hence, the control signal u is also periodic with amplitude A and frequency ω. If N(A,ω) denotes the DF of the nonlinearity 5 , by the classical frequency domain analysis, y = C(jω) N(A,ω) P(jω)(0 − y) so that if G = C P and y ≠ 0, then

G(jω) = −

1 N (A , ω)

(2)

Equation 2 is called the harmonic balance equation. It follows that a necessary condition for the existence of a limit-cycle in the system presented in Figure 1 is the existence of a pair (A0,ω0) satisfying eq 2. The prediction of the frequency and amplitude of the oscillation of u is achieved through the resolution for (A,ω) of eq 2. For a large class of nonlinearities which includes the twoparameter stiction model considered in this study, the DF does not depend on the frequency of the input signal u, that is, N(A,ω) = N(A). In the DF definition, the input u is considered to be sinusoidal. In the general case, the oscillation of u due to the nonlinearity is not sinusoidal. Thus, the reliability of the DF method is directly related to how well the oscillation of u can be approximated by a sinusoidal oscillation. One way to check this condition is by calculating the total harmonic distortion of the signal u (THD(u)). If u is sinusoidal, then THD(u) = 0, and the closer this measure is to zero, the better is the approximation of u by a sinusoid. A condition for u to be well approximated by a sinusoid is that the linear part G has good low-pass filtering properties. This is a condition met by many linear processes controlled by PID controllers. There are certain conditions of nontransversality,9 which may lead to erroneous predictions in the DF method; for example, if one considers the negative inverse curves of the describing function (A → −1/N(A)) parametrized by the amplitude and the Nyquist graph of G. Equation 2 admits a solution for A and ω if, and only if, these curves intersect (Figure 3). However, if this intersection is not transverse (i.e., tangent) then the system may or may not reveal oscillations due to the approximate nature of the DF method. In an earlier study,10 the following DF for the two-parameter stiction model was derived:

Figure 2. Typical input−output behavior of a sticky valve.

The set of physically relevant values for the stiction parameters S and J is 3 = {(S , J ); 0 < J ≤ S} When S = J = 0, there is no stiction and the case J = 0 < S is the dead band case, which may also occur in control valves.5 A model for control valve stiction is a nonlinear operator that reproduces, at least roughly, the main stiction characteristics shown in Figure 2. Stiction models can be divided into two categories: physics-based models and data-driven models. Physical models involve several parameters that are normally difficult to obtain in practice and usually require a great computational effort in simulations, making them almost useless for routine industrial applications. In comparison, data-driven models typically depend on at most two parameters and are easy to implement and demand less computation effort. This study employs the two-parameter data-driven model of valve stiction.7 The input/output graph of this model agrees with the graph shown in Figure 2. Therefore, it can be concluded that it reproduces the expected characteristics of the control valve’s stiction in accordance with the above-mentioned definition. This model also has the advantage of reproducing the dead band behavior when J = 0, allowing the approach of these two nonlinearities in a unified manner. 14122



■

Article

STICTION AND DEAD BAND QUANTIFICATION BY THE DESCRIBING FUNCTION METHOD Let us consider the nonlinear control loop presented in Figure 1 and assume that the control signal u is oscillating with amplitude A0 and frequency ω0. If the oscillation of u is due to stiction, then by the DF method, the harmonic balance equation G(jω) = −

1 NSJ (A)

(5)

or, equivalently, the equation NSJ (A) = −

with

⎛π ⎞ 1 ⎛A ⎜ sin 2ϕ − 2A cos ϕ − A⎜ + ϕ⎟ ⎝2 ⎠ πA ⎝ 2 ⎞ + 2(S − J )cos ϕ⎟ ⎠

Nd(A) = −

X (A ) = −

1 ⎛⎜ A A −3 + cos 2ϕ + 2A sin ϕ 2 πA ⎝ 2 ⎞ − 2(S − J )sin ϕ⎟ ⎠

⎛ A − S⎞ ⎟ ϕ = arcsin⎜ ⎝ A ⎠

(7)

■

CONSISTENCY OF THE ESTIMATING EQUATIONS BASED ON DESCRIBING FUNCTIONS Stiction. In the DF expression of the stiction model (eq 3), the variables S and J always appear in the form S/A and J/A. In the estimation procedure by the DFSQ method, the amplitude A is fixed and equal to A0, where A0 is the measured amplitude of the oscillating control signal. Thus, to simplify the analysis, the new variables x = S/A and y = J/A and the map N : D ⊂ 2 → 2 , N(x,y) = X(x,y) + jY(x,y), are considered, which are defined by

Using some trigonometric identities one can simplify the previous expression to obtain the simplest formula (3) for the DF of the two-parameter stiction model: NSJ(A) = X(A) + j Y(A) with ⎛ J⎞ S ⎞⎛ S 1 ⎛π ⎜ + arcsin⎜1 − ⎟⎜1 − +2 ⎟ ⎝ A ⎠⎝ A A⎠ π⎝2 S ⎛⎜ S ⎞⎞ 2 − ⎟⎟ A⎝ A⎠⎠

X (A ) =

X (x , y ) = (3)

The DF for the dead band with parameter d, Nd(A) = X(A) + jY(A) is obtained from previous DF for the stiction model by making J = 0 and S = d; that is,

Y (x , y ) =

1 ⎛⎜ π + arcsin(1 − x) + (1 − x + 2y) π⎝2 ⎞ x(2 − x) ⎟ ⎠ 1 ( −x(2 − x) + 2y(1 − x)) π

(8)

From eqs 8 it can be easily observed that 0 ≤ x ≤ 2. The physically relevant case occurs when 0 ≤ J ≤ S. Thus, 0 ≤ y ≤ x. Therefore, the domain of map N(x,y) is the closed triangular region of the plane 2 , defined by D = {(x,y); 0 ≤ y ≤ x, 0 ≤ x ≤ 2}. As observed in the previous section, the parameters S and J of the stiction model are determined by solving eq 6 for (S,J), with A = A0 and ω = ω0. Using the new variables x = S/A0 and y = J/A0,

⎛ d ⎞⎛ d⎞ d⎛ d ⎞⎞ 1 ⎛π ⎜2 − ⎟⎟ X(A) = ⎜⎜ + arcsin⎜1 − ⎟⎜1 − ⎟ ⎟ ⎝ A ⎠⎝ A⎠ A⎝ A⎠⎠ π⎝2 Y (A ) = −

1 G(jω)

is satisfied for some d > 0, at a frequency ω ≈ ω0 and an amplitude A ≈ A0. Similar to the stiction case, quantifying the dead band is done by solving eq 7 for d with A = A0 and ω = ω0. The methodology for stiction quantification based on resolution for S and J of eq 5 will from now on be called the describing function stiction quantification (DFSQ) method, and the same methodology applied to dead band quantification will be called the describing function dead band quantification (DFDBQ) method.

Y (A ) =

J⎛ S⎞ S ⎞⎞ 1⎛ S⎛ Y (A ) = ⎜ − ⎜ 2 − ⎟ + 2 ⎜ 1 − ⎟⎟ ⎝ ⎠ ⎝ ⎝ A A A A ⎠⎠ π

(6)

is satisfied for some pair of parameters (S,J) at a frequency ω ≈ ω0 and an amplitude A ≈ A0. The stiction quantification is performed by estimating the parameters S and J in eq 6. To this end, eq 6 is solved for (S,J) with A = A0 and ω = ω0. The solution (Ŝ,J)̂ , if it exists, produces the estimated parameters of the stiction model. If the nonlinearity in the control loop presented in Figure 1 is dead band and the control signal oscillates with amplitude A0 and frequency ω0, then by the DF method, the equation

Figure 3. Intersection between negative inverse curves of DF and Nyquist graph of G.

NSJ (A) = X(A) + jY (A),

1 G(jω)

d ⎛ d⎞ ⎜2 − ⎟ ⎝ A⎠ πA (4) 14123



Article

f |∂U is injective. Then f homeomorphically maps U onto f(U). In particular, f |U is injective. In what follows, it is shown that the map N(x,y) satisfies the hypotheses of the Massey Theorem. Lemma 1. N : int(D) → 2 is a local diffeomorphism. Proof: The derivative of the map N(x,y) at a point (x,y) ∈ int (D) is given by

the problem of solving eq 6 is equivalent to that of solving the eq 9 for (x,y) 1 N (x , y ) = − G(jω0) (9) If (x̂,ŷ) is a solution for eq 9, then Ŝ = A0x̂ and J ̂ = A0ŷ are the estimated parameters of the stiction. Therefore, the consistency analysis of the estimator equation for the parameters S and J depends on the behavior of the application (x,y) → N(x,y). Existence of Solutions. The problem of the existence of solutions for eq 9 can be described in terms of the surjectivity of the map N(x,y); that is, there is (x,y) satisfying eq 9 if, and only if, −1/G(jω0) belongs to the image N(D) or, equivalently if, and only if, G(jω0) belongs to the set −1/N(D) = {−1/N(x,y);(x,y) ∈ D}. The set −1/N(D) = {−1/N(x,y);(x,y) ∈ D} is the unbounded region of the complex plane depicted in Figure 4. The curve Γ1 is

⎡ 2 ⎢− 2(2x − y + xy − x ) 2 x(2 − x) ⎢ π π 2x − x 2 N ′(x , y) = ⎢ ⎢ 2(y − x + 1) 2(x − 1) − ⎢⎣− π π

⎤ ⎥ ⎥ ⎥ ⎥ ⎥⎦

The Jacobian determinant of N(x,y) at a point (x,y) ∈ int (D) is det N ′(x , y) =

4y π

2

x(2 − x)

(10)

It follows from eq 10 that det N′(x,y) = 0 in int (D) if, and only, if y = 0. Therefore, N′(x,y) is nonsingular at every point (x,y) ∈ int (D). By the inverse function theorem, N|int(D) is a local diffeomorphism and, therefore, a local homeomorphism. ■ Theorem 1. The map N:D → N(D), defined in eq 8, satisfies the hypotheses of Massey Theorem and, therefore, is a homeomorphism. Proof: Using the same notations as in the Massey Theorem, with f(x,y) = N(x,y) defined by the eqs 8 and 9, U = int (D) (the interior of the triangular domain D) and ∂U the triangle with vertices (0.0) (2.0) and (2.2), it follows that (i) n = 2, U ⊂ 2 is open, with U̅ = U∪∂U compact. Moreover, f = N is continuous on U, because its coordinate functions X(x,y) and Y(x,y) are continuous on U. (ii) f = N is a local homeomorphism on U = int (D), due to Lemma 1. It only remains be proven that f |∂U is injective. Let us represent the boundary ∂U by three rectilinear parametrized paths: (i) γ1(t) = (t,0), 0 ≤ t ≤ 2 (horizontal segment connecting (0,0) to (2.0)); (ii) γ2(t) = (2,t), 0 ≤ t ≤ 2 (vertical segment connecting (2,0) to (2.2)); (iii) γ3(t) = (t,t), 0 ≤ t ≤ 2 (straight line segment connecting (0,0) to (2.2)). It can now be shown that f(γi(t)), where i = 1,2,3 is injective on the interval [0,2]: (i) f(γ1(t)) = N(t,0) = X(t,0) + iY(t,0) is injective on the interval [0,2]. In fact, it is sufficient to prove that the component X(t,0) is injective on the interval [0,2]. Since X′(t,0) = −2(2t − t2)1/2/π < 0, ∀ t ∈ (0,2), X(0,0) = 1 e X(2,0) = 0, it follows that X(t,0) is decreasing in [0,2], in particular X(t,0) is injective. (ii) f(γ2(t)) = N(2,t) = X(2,t) + iY(2,t) is injective on the interval [0, 2]. In fact, it is sufficient to prove that the component Y(2,t) is injective on the interval [0,2]. Y(2,t) = −2t/π is clearly injective on the interval [0,2]. (iii) f(γ3(t)) = N(t,t) = X(t,t) + iY(t,t) is injective on the interval [0,2]. In fact, it is sufficient to prove that the component Y(t,t) is injective on the interval [0,2]. However, Y(t,t) = −t2/π, and, therefore Y(t,t) is clearly injective on the interval [0,2]. It follows from (i), (ii), and (iii) that f |∂U is injective. ∂U is not a one-dimensional differentiable manifold because it is not a closed smooth curve. However, ∂U can be approximated by closed smooth curves (Figure 5), along with the property of f restricted to these curves being injective.

Figure 4. Stiction region.

the image of the parametrized curve x → −1/N(x,0), x ∈ [0,2], which coincides with the negative inverse curve of the dead band DF. The curve Γ2 is the image of the parametrized curve x → −1/ N(x,x), x ∈ [0,2]. If G(jω0) ∈ Γ1 then the estimations are of the form (x̂,0) or (Ŝ,0), that is, dead band only and when G(jω0) ∈ Γ1 the estimations are of the form (x̂,x̂) or (Ŝ,J)̂ with Ŝ = J.̂ The curves Γ1 and Γ2 and the negative y-axis form the boundary of a region called stiction region (Figure 4). Therefore, eq 9 admits solution if, and only if, the frequency response G(jω0) belongs to the stiction region. Uniqueness of Solutions. The problem of the uniqueness of solutions of eq 9 can be equivalently described in terms of the injectivity of the map N(x,y) as defined in eq 8. In fact, if N:D → N(D) is injective, then given any z0 ∈ N(D), there exists a unique z = (x,y) ∈ D such that N(x,y) = z0. The problem of determining when a map of the plane is injective is quite delicate and involves certain assumptions about the behavior of the map in its domain and the domain boundary. W. S. Massey’s theorem from Topology11 identify the conditions sufficient for a local homeomorphism to be injective. Theorem (Massey). Let U ⊂ n be an open subset with compact closure U̅ and f : U̅ → n be a continuous map such that f |U is a local homeomorphism. Suppose that ∂U is a closed, connected, and orientable manifold of dimension n − 1 such that 14124



Article

ϕ(d)̂ = min ϕ(d) = min |Nd(A 0) + 1/G(jω0)|2 0 ≤ d ≤ 2A 0

0 ≤ d ≤ 2A 0

(13)

If G(jω0) is close enough to the negative inverse curve of the dead band DF Γ1 then the global minimum d̂ is close to the parameter d corresponding to the intersection between the Nyquist plot of G and Γ1. By the DF method this d parameter corresponding to the intersection is an approximation to the real dead band. From expressions of the dead band DF (eq 4) the domain of the function ϕ is the interval [0, 2A0]. The existence of global minimum in eq 13 is guaranteed by the continuity of ϕ and the compactness of the interval [0,2A0]. Thus, the estimation of d parameter is guaranteed. Numerical Methodology. Dead Band Estimation. In the dead band DF (eq 4), the change of variable x = d/A is considered. Then, the expression of the DF becomes N(x) = X(x) + jN(x), with

Figure 5. Smoothing a vertex in ∂U.

Therefore, N(x,y) satisfies the hypotheses of the Massey Theorem, and is, therefore, a homeomorphism. ■ Corollary 1: If G(jω0) belongs to a stiction region then eq 9 admits a unique solution (x̂,ŷ) and also, (Ŝ = A0x̂) and J ̂ = A0ŷ are the estimated stiction parameters. Proof: If G(jω0) belongs to the stiction region then −1/G(jω0) belongs to image N(D). From Theorem 1 N:D → N(D) is a homeomorphism, therefore there exists a unique (x̂,ŷ) ∈ D such that N(x̂,ŷ) = −1/G(jω0). Since x = S/A and y = J/A one has to consider that Ŝ = A0x̂ and J ̂ = A0ŷ. Because Theorem 1 is true on the boundary of the stiction region, one has the following. Corollary 2: If G(jω0) belongs to the curve in Γ1 then there exists a unique x̂ ∈ [0,2] such that G(jω0) = −1/N(x̂,0) and d̂ = A0x̂ (dead band only). Observation 1. If G(jω0) belongs to the negative imaginary axis then this case corresponds to x = 2 in eq 8. This fact significantly simplifies the solution of eq 9; in fact, in this case the estimated parameters are given by

S ̂ = 2A 0 Ĵ = −

πA 0 2 imag(G(jω0))

X (x ) =

⎞ 1 ⎛⎜ π + arcsin(1 − x) + (1 − x) x(2 − x) ⎟ ⎝ ⎠ π 2

1 Y (x) = − x(2 − x) π (14)

where 0 ≤ x ≤ 2. If G(jω0) is the frequency response of the linear part of the control loop corresponding to the oscillation frequency ω0, then, as presented in section about stiction and dead band quantification by the DF method, to quantify the dead band by the DF method, the harmonic balance eq 7 is solved for d using the measured amplitude A = A0 and measured frequency ω = ω0. This problem is equivalent to solving eq 15 for x with −1/ G(jω0) = x0 + jy0.

(11)

N (x) = X(x) + jY (x) = x0 + jy0

(15)

In this case, the estimated dead band is given by d̂ = A0x̂, with x̂ being the solution of eq 15. The problem of solving the nonlinear eq 15 for the variable x according to eq 14, is transformed into a minimization problem in one variable. Let us considering the function F defined by

(12)

Dead Band. It is assumed that the nonlinearity in the control loop in Figure 1 is dead band and that the loop is oscillating with the control signal u having amplitude A0 and frequency ω0. By the DF method, the Nyquist graph of G intersects the negative inverse curve of the DF of the dead band. This intersection occurs at a point corresponding to an amplitude Ai and a frequency ωi, so that G(jωi) = −1/Nd(Ai), for some d > 0. If the intersection is transversal and G has good low-pass filtering properties, then Ai ≈ A0 and ωi ≈ ω0; that is, the amplitude and frequency corresponding to the intersection of the curves are very close to the amplitude and frequency measurements of the control signal. There are two cases to be considered: (1) G(jω0) ∈ Γ1, that is, the frequency response G(jω0) belongs to the negative inverse curve of the DF of the dead band. In this case, as a consequence of Corollary 2, there is a unique d̂ > 0 satisfying eq 7, with A = A0 and ω = ω0. (2) G(jω0) ∉ Γ1. This situation is much more common, because any regular curve has zero measure in the plane, and measurement errors in ω0 rarely occur in G(jω0) ∈ Γ1. Let us consider the function ϕ(d) = |Nd(A0) + 1/G(jω0)|2. In the previous case (G(jω0) ∈ Γ1), there exists a unique d̂ > 0 satisfying eq 7 or equivalently ϕ(d̂) = 0 where d̂ is the global minimum of the function ϕ. Now if G(jω0) ∉ Γ1 then there does not exist a parameter d̂ such that ϕ(d̂) = 0. Therefore ϕ(d) > 0 for all d. In this case the proposed dead band estimation d̂ is the global minimum of the function ϕ, that is,

F(x) = (X(x) − x0)2 + (Y (x) − y0 )2

Then, x̂ is a solution of eq 15 if, and only if, F(x̂) = 0. Since F(x) ≥ 0 for all x, it follows that x̂ can be estimated as global minimum of the function F, that is, x̂ is such that F(x)̂ = min F(x) 0≤x≤2

(16)

As shown in the previous subsection, if G(jω0) ∈ Γ1, that is, the frequency response G(jω0) belongs to the negative inverse curve of the DF of the dead band Γ1 shown in Figure 4, then the solution x̂ in eq 15 (d̂ = x̂A0) satisfies F(x̂) = 0. On the other hand, if G(jω0) ∉ Γ1, then F(x) > 0 ∀x and the global minimum in eq 16 satisfies F(x̂) > 0. As seen in case 2 of the previous section, if G(jω0) is close enough to Γ1, then the estimation d̂ = x̂A0 is close to the d parameter corresponding to the intersection between the Nyquist curve of G(s) and the Γ1. This is an approximation to the real dead band by the DF method. Therefore, in any case, the d parameter of the dead band always can be estimated. Estimation of Parameters S and J of the Two-Parameter Stiction Model. This is similar to the dead band case. If G(jω0) is the frequency response of the linear part corresponding to the 14125



Article

oscillation frequency ω0, and A0 is the measured amplitude of the control signal, then to estimate the parameters S and J of the stiction model by the DF method, the harmonic balance eq 6 is solved for (S,J), using A = A0 and ω = ω0. Considering the DF expression of the stiction model in the variables x = S/A and y = J/A in eq 8 and x0 + jy0 = −1/G(jω0), the estimation procedure consists in solving eq 17 for (x,y) in the triangular domain D = {(x,y);0 ≤ y ≤ x, 0 ≤ x ≤ 2}. N (x , y) = X(x , y) + jY (x , y) = x0 + jy0

where KN(s) represents the four Kharitonov polynomials associated with N(s), KD(s) represents the four Kharitonov polynomials associated with D(s), SN(s) and SD(s) are segment plants defined by SN(s) = λKiN(s) + (1 − λ)KkN(s), (i,k) ∈ {(1,2), (1,3),(2,3),(3,4)} and S D (s) = λK Di (s) + (1−λ)K Dk (s), (i,k)∈{(1,2),(1,3),(2,3),(3,4)}, in both cases 0 ≤ λ ≤ 1. In this section, the frequency domain template set G(jω) plays an important role because it substitutes the frequency response point in the estimation procedure. The following Theorem correlates the border of the frequency domain template set with the extremal transfer function (eq 18) and justifies using the extremal transfer function GE(s) to construct the border of the frequency domain template ∂G(jω). Theorem: Using the notations above

(17)

The problem of solving eq 17 is converted to the following minimization problem: Consider the function F(x , y) = (X(x , y) − x0)2 + (Y (x , y) − y0 )2

∂G(jω) ⊂ G E(jω)

where ∂C denotes the boundary of set C. Estimation Procedure. Consider the control loop presented in Figure 1. It is assumed here that the nonlinearity is the dead band or the two-parameter stiction model and in both cases Nλ(A) denotes its DF with λ = d in the dead band case and λ = (S,J) in the stiction case. Let us suppose that the control and process transfer functions are not exactly known but are defined by interval transfer functions C(s) and P(s), respectively. In this case the linear part transfer function is also an interval transfer function G(s) = C(s)/P(s). If the control signal u is oscillating due to nonlinearity with measured amplitude A0 and measured frequency ω0, then the harmonic balance equation

Then, (x′,y′) is a solution of eq 17 if, and only if, F(x′,y′) = (0,0). In this case (x′,y′) is a global minimum for F. On the basis of this fact, the estimation procedure adopted here consists of minimizing F(x,y) for (x,y) in the triangular domain D. The estimated parameters of the stiction model are given by Ŝ = A0x̂ and J ̂ = A0ŷ. Unlike the dead band case, where the DF curve Γ1 is fixed and rarely contains the frequency response point G(jω0), in the stiction case the negative inverse curve of the DF of the stiction with the estimated parameters Ŝ and J,̂ always contains the frequency response point G(jω0). Therefore, in both cases the problem of estimating the nonlinearity parameters is converted into an optimization (minimization) problem, in one variable in the case of dead band, and in two variables in the case of stiction.

G(jω0) = −

■

KN (s) S (s ) ∪ N SD(s) KD(s)

1 Nλ(A 0)

(20)

is satisfied for some parameter value λ, where G(s) denotes the unknown nominal linear part transfer function. In the uncertain case the nominal frequency response G(jω0) is an interior point of the frequency domain template G(jω0) and if the uncertainty intervals that define the interval transfer function G(s) are small enough, then the set G(jω0) is well approximated by a polygon whose vertices belong to the set formed by the 16 frequency responses of the Kharitonov plants Gik(jω0) associated with the interval transfer function G(s). The proposal here is to solve the harmonic balance equation for λ for each frequency response Gik(jω0):

THE CASE OF AN UNCERTAIN LINEAR PART In the previous section, the case of the estimation of the dead band and stiction under exact knowledge of the linear part of the loop was considered. In practical situations, the controller is usually known but the process transfer function is not precisely known. In this section, the methodologies for the estimation of the dead band and stiction developed in the previous section will be generalized to deal with the case of uncertainties in the process transfer function. A well-known technique from robust control theory12 will be used here. Kharitonov Theory. Uncertain plants can be represented by interval transfer functions, these in turn are defined by quotients of interval polynomials G(s) = N(s)/D(s). For each interval polynomial N(s) and D(s) there exists four extremal polynomials called Kharitonov polynomials, denoted by KiN(s), i = 1,2,3,4 for N(s) and KkD(s), k = 1,2,3,4 for D(s). Several properties of the interval polynomials, which are infinite families of polynomials, can be deduced from the properties associated with the four Kharitonov polynomials. If G(s) = N(s)/D(s) is an interval transfer function then G(jω) is a compact region in the complex plane called the frequency domain template, and the boundary of this region is formed by segments, arcs, and vertices. The set formed by vertices of the boundary of the frequency domain template is contained in the set formed by the 16 frequency responses of the Kharitonov plants Gik(jω) = KiN(jω)/KkD(jω), i, k = 1,2,3,4. Associated with the interval transfer function G(s) = N(s)/ D(s), one has the extremal transfer function GE(s), defined by GE(s) =

(19)

Gik (jω0) = −

1 Nλ(A 0)

, i , k = 1, 2, 3, 4

(21)

If λik̂ , i, k = 1,2,3,4 denotes the 16 solutions of eq 21, then considering the extremal values min{λik̂ } and max{λik̂ } it is possible to determine the estimation intervals of the parameters: Id = [min{λik̂ },max{λik̂ }] for λ = d in the dead band case and IS = ̂ }], IJ = [min{λ2ik ̂ },max{λ2ik ̂ }] for λ = (S,J) in [min{λ̂1ik},max{λ1ik the stiction case. The justification for the validity of this estimation procedure is based on convexity arguments applied to the frequency domain template G(jω 0 ) which is well approximated by a polygon for small uncertainty intervals. Figure 6 illustrates a frequency domain template G(jω0) associated with an interval transfer function G(s) contained in the stiction region, with one of its vertices Gik(jω0) displayed.

■

THE CASE OF UNKNOWN PLANT The estimation of stiction by DF method requires the controller and process models and parameters. The controller parameters are available in the digital controller. When only routine data are

(18) 14126



Article

One can see that tuning rules in general put a zero between the origin and the pole of the transfer function, usually close to the pole. When the Nyquist curve of G(s) = C(s)P(s) is plotted, the proximity of pole and zero results in a plot that can be approximated by the Nyquist plot of G(s) = C(s)P(s) ≈

KK p Tsi

(29)

This approximation becomes even better if the oscillation frequency due to stiction is smaller than 1/Ti, which is common as the examples from industry will show. This fact simplify enough the estimation of the stiction parameters, because in this case S and J can be calculated directly using the eqs 30 and 31, which follows directly from eqs 11, 12, and 29. Figure 6. A frequency domain template.

S ̂ = 2A 0

available, the controller parameters can be obtained via closedloop identification,13 using SP, OP, and PV. The model for the process operating around the set point can be obtained via an open-loop step response. This is a usual procedure for controller tuning. When only input−output data are available, the joint closed loop identification of linear process and control valve should be performed.1 If the procedure succeeds, stiction parameters are estimated. However, this methodology is computationally cumbersome. The procedure for modeling here proposed is quite simpler and is strongly based on common practices in industrial processes. PID controllers are used in more than 95% of the industrial control loops.14 Additionally, each control loop uses a well-known controller that assures the desired performance. Therefore, the controller structure is normally known, which makes easy the identification of the controller parameters. Flow control loops, e.g., use in general a PI controller. The transfer function of a flow process is given by K P(s) = (22) τs + 1

Ĵ =

(23)

T Kp = i Kλ

(24)

The use of the frequency response method relates the module and angle of y(t) and u(t), respectively, by |Y (jω)| = |P(jω)||U (jω)| =

θ 2

(28)

Ti = τ +

(32) (33)

The time constant τ can be calculated by tan(φ(jω)) ω and its substitution in eq 32yields τ=

|K | =

|Y (jω)| tan(φ(jω))2 + 1 |U (jω)|

(34)

(35)

Actually, eqs 34 and 35 can be used to calculate the parameters of a transfer function P(s) if the signal y(t) is known for a given signal u(t) = A sin(ωt). The application of the signal u(t) = 6 sin(2t) with P(s) = (2/(s + 1)) and stiction parameters S = 5 and J = 1 in the system of Figure 7, results in the signals v(t) and y(t) shown in Figure 8. The phase of v(t) is delayed and the module of v(t) is decreased by the nonlinearity, with the same effect in y(t). The effect of the nonlinearity on module and angle can be predicted by the DF method. The module and phase of the stiction DF (eq 3) with A = 8, S = 7, and S/J = 1, ..., 14 is shown in Figure 9. One can see that S/J reduces (increases) the module and lags (leads) the angle of the stiction DF. If greater values of A are used, this relation is preserved, but the effect of S/J on module and angle is reduced.

and for PI controller the IMC method gives (27)

|U (jω)|

φ(jω) = arg(Y (jω)) − arg(U (jω)) = arctan(ωτ )

(25)

2τ + θ 2Kλ

|K | 2 2

ωτ +1

When the time delay is not negligible, the process is normally well represented by K −θs P(s) = e (26) τs + 1

Kp =

(31)

Figure 7. Open loop stiction + linear process used to approximate calculation of gain K.

where λ is the desired closed loop constant, and the controller C(s) is given by ⎛ 1 ⎞ C(s) = K p⎜1 + ⎟ Tsi ⎠ ⎝

πTA i 0 ω0 2KK p

To use eq 30 and eq 31, the process gain K is required, which is done in the sequence. Let us consider in the block diagram shown in Figure 7, the process P(s) = (K/(τs + 1)) and the stiction parameters S = J = 0 (linear).

The time delay is negligible, and if exists, is due to stiction. The time constant τ comes from the control valve actuator. A suitable tuning for this control loop is given by the IMC method,15 with

Ti = τ

(30)

14127



Article

Therefore, this is the method that will be used in this paper to calculate the process gain K to be used in eq 31.

■

APPLICATIONS Dead Band Estimation in a Control Loop with a PI Controller and an Integrating Process. A nonlinear control loop, like the one presented in Figure 1 is considered. The system was simulated with the following elements and parameters: Process: 1 P(s) = s(s + 1)

Figure 8. Signals of the system of Figure 7 with stiction for sinusoidal input.

PI Controller: 0.5s + 0.1 C(s) = s The nonlinearity 5 is the deadband. Table 2 shows the estimation results. d̂ is the estimated dead band, ed(%) is the estimation error percentage and A0 is the Table 2. Results of Dead Band Estimations

Figure 9. Angle and module of the stiction DF.

This observation together with eq 35, allows one to conclude that the calculation of K using u(t) is approximate to the calculation of K using v(t), since when S/J increases, |Y(jω)|/| V(jω)| is reduced while arg(Y(jω)) − arg(V(jω)) is increased. This analysis was performed also via simulation, using P(s) = e−0.2/(s + 1) and C(s) = (s + 2)/s from the section ″Stiction Estimation in a Control Loop with PI Controller and First Order Process with Time Delay″, with Ts = 0.1 s. The transfer function P(s) was estimated as an ARX model using signals u from controller output and y from process output. The gain of this transfer function is shown in Table 1 for S = 7 and S/J = 1, ..., 14. The gain K calculated in all situations was close to the real value 1.

d (%)

d̂

ed (%)

A0

5 10 20 30 40 50 60 70 80 90

4.91 9.81 19.62 29.44 39.25 49.06 58.87 68.68 78.49 88.31

1.8830 1.8833 1.8829 1.8825 1.8823 1.8821 1.8820 1.8819 1.8818 1.8817

0.0395 0.0790 0.1580 0.2370 0.3160 0.3950 0.4740 0.5531 0.6321 0.7111

amplitude of the control signal. The shape of the negative inverse curve of the DF of the dead band (Figure 10) does not depend on

Table 1. Calculation of Gain K for Different Values of S/J S/J

K

14.0000 7.0000 4.6667 3.5000 2.8000 2.3333 2.0000 1.7500 1.5556 1.4000 1.2727 1.1667 1.0769 1.0000

0.9997 1.0028 0.9981 1.0022 1.0011 1.0030 1.0022 1.0055 1.0038 1.0041 1.0044 1.0087 1.0079 1.0074

Figure 10. Dead band negative inverse curve and Nyquist plot of the linear part G.

the parameter d, since in the expression of the DF in eq 4 the parameter always appears in the form d/A. The most interesting consequence of this fact is that the parameter d influences the oscillation amplitude but, since the shape of the curve does not change, the point of intersection with the Nyquist plot of G does not change, and hence the frequency does not change either. In 14128



Article

Table 3. Results of Stiction Parameters Estimations S

Ŝ

eS(%)

J

Ĵ

eJ(%)

A

Ai

eA(%)

ω

ωi

eω(%)

10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0

9.3271 8.8338 9.4838 9.7684 9.8797 9.9337 9.9617 9.9763 9.9836 9.9877

6.7290 11.6624 5.1622 2.3159 1.2028 0.6632 0.3833 0.2374 0.1645 0.1231

0.5 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

0.0001 0.0001 1.6698 3.0266 4.2008 5.3017 6.3632 7.4011 8.4239 9.4376

99.9743 99.9910 16.5081 0.8858 5.0192 6.0338 6.0539 5.7303 5.2987 4.8619

5.9847 5.7671 5.5649 5.5179 5.5200 5.5406 5.5696 5.6029 5.6383 5.6753

6.3839 6.1342 5.8111 5.6661 5.5937 5.5824 5.5849 5.6136 5.6474 5.6661

6.2532 5.9841 4.2374 2.6159 1.3158 0.7481 0.2732 0.1920 0.1612 0.1617

0.1555 0.1643 0.2035 0.2553 0.3084 0.3598 0.4086 0.4547 0.4983 0.5396

0.1785 0.1849 0.2112 0.2519 0.2972 0.3448 0.3920 0.4377 0.4815 0.5233

12.8917 11.1478 3.6634 1.3498 3.7824 4.3345 4.2484 3.8938 3.4848 3.1058

these simulations the measured frequency of the control signal was ω0 = 0.2701 rad/s. As the shape of the DF curve of dead band does not change, it is possible to determine the frequency ωi corresponding to the intersection between the Nyquist curve of G and the negative inverse curve of the DF of the dead band. If the DF is a good approximation then ωi must be close to ω0. In these simulations ωi = 0.2854 rad/s and the percentage variation from ω0 relative to ωi was eω = 5.36%. This small value confirms that the DF method gave a good approximation and that the dead band estimates are reliable. Stiction Estimation in a Control Loop with PI Controller and an Integrating Process. The system was simulated with the following elements and parameters: Process:

P(s) =

1 s +s 2

PI Controller:

C(s) =

Figure 11. High sensitivity region for computation of J.

0.5s + 0.05 s

The system was simulated with the following elements and parameters: Process: 1 −0.2s P(s) = e s+1

The nonlinearity 5 is the two-parameter data-driven model of valve stiction.7 Table 3 shows the estimation results. Ŝ e J ̂ are the stiction parameters used in each simulation, Ŝ and J ̂ are the estimated parameters, eS and eJ are the respective percentage errors. A and ω are respectively the amplitude and frequency of the control signal. Using S, J, and G(s) = C(s)P(s) it was possible to predict the amplitude Ai and frequency ωi by the DF method. These values were compared with A and ω to produce the percentage variations eA(%) and eω(%). As shown in Table 3, the algorithm estimates the parameter S with error eS less than 11.7% in all cases. For J = 0.5 and J = 1.0 the algorithm was not able to estimate the parameter J. This behavior is due to the proximity between the frequency response G(jω) and the negative inverse curve of the DF of the stiction with J = 0 (or dead band curve) (Figure 11). A small variation Δω in the measured frequency induces a variation Δz in the frequency response, which induces variations ΔS and ΔJ in the estimations. To assess this sensitivity, an increase of five percent was done in each frequency value ω, as presented in Table 2. Then, the percentage variations ΔS and ΔJ were calculated and plotted as a function of J used in the simulations (Figure 12). It follows from this analysis that, if the frequency response corresponding to oscillation is close to the dead band curve, then the estimated values for J are not reliable. Stiction Estimation in a Control Loop with a PI Controller and a First Order Process with Time Delay.

PI Controller: s+2 C(s) = s The nonlinearity 5 is the two-parameter data-driven model of valve stiction.7 The stiction parameters used in simulations are S = 10, J = 0.5, 1.0, ..., 9.0. Figure 13 shows the Nyquist plot of the linear part G(s) = C(s)P(s) and the DF curves corresponding to the parameters S = 10 and J = 0.5, 1.0, ···, 9.0. As shown in Figure 13, the Nyquist plot of the linear part is not close to the dead band curve, so more accurate estimates are expected, which is confirmed by the results shown in Table 4. As shown in Table 4 the stiction parameters S and J were estimated with a maximum error of less than 13.5%. Stiction Estimation in a Pilot Plant with Uncertainty in the Linear Part. The pilot plant diagram and the control valve used in the experiments are shown in Figure 14. The control valve, its stem and the packing are shown in details in Figure 14b,c. The friction was introduced into the valve by tightening the packing nut. The flow transmitter and the I/P transducer 14129



Article

Figure 12. Error on S and J for different values of S/J (S = 10 and J varies).

The parameters of the process transfer function were identified through step tests performed in the plant. The identification procedure has introduced a parametric uncertainty of 10% so that the process transfer function is represented by the following uncertain transfer function P(s) =

[1.3878, 1.6962] [0.3369, 0.4117]s + 1

The linear part transfer function G = C(s)P(s) is represented by the uncertain transfer function G(s) =

[0.8998, 1.0998]s + [2.4051, 2.9395] [0.3369, 0.4117]s 2 + s

(36)

The Nyquist template G(jω0) associated with the uncertain transfer function (eq 36) contained in the stiction region is shown in Figure 15. Γ1 denotes the negative inverse curve of the DF of the stiction with J = 0. Γ2 denotes the negative inverse curve of the DF in case S = J. Two DF curves corresponding to the minimum and maximum values of the parameter J are also shown. The vertices of the Nyquist template G(jω0) were calculated using RPC ToolBox,12 a MATLAB toolbox used to deal with problems involving the robust stability of uncertain systems. The algorithm was applied to each vertex of G(jω0) resulting in a set of estimates {(Ŝi,Jî )}. The resulting estimation intervals IS = [min Si,max Si] and IJ = [min Ji,max Ji] are IS = [12.502, 12.506] IJ = [1.01, 1.23].

Figure 13. Nyquist plot and the DF curves for the example.

works in the range 4−20 mA. The industrial controller CompactRIO from National Instruments was used to control and monitor all signals. The controller transfer function is given by C(s) = ((0.6484s + 1.733)/s). Table 4. Results of Stiction Parameters Estimations S

Ŝ

eS (%)

J

Ĵ

eJ (%)

A

ω

10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0

9.9747 9.9058 9.7123 9.5203 9.3577 9.2241 9.1120 9.0165 8.9310 8.8536

0.2531 0.9416 2.8767 4.7969 6.4228 7.7589 8.8796 9.8349 10.6899 11.4643

0.5 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

0.5670 1.0598 1.9580 2.8471 3.7722 4.7386 5.7434 6.7837 7.8436 8.9319

13.4025 5.9799 2.0991 5.0964 5.6953 5.2285 4.2764 3.0897 1.9553 0.7564

5.0013 5.0026 5.0052 5.0077 5.0098 5.0120 5.0138 5.0162 5.0178 5.0202

0.1502 0.2865 0.5248 0.7292 0.9123 1.0819 1.2431 1.3990 1.5504 1.7001

14130



Article

Figure 14. Pilot plant.

Figure 15. Frequency domain template of the linear part of pilot plant.

The stiction parameters obtained from off line tests on the valve were S = 12 and J = 1, confirming the good estimations provided by the algorithm when used in the presence of uncertainty. Application to Real Data from Industry. In this section, we show the performance of the proposed method on real industrial data. This data set was used previously as a benchmark for stiction detection methods.16 The methodology developed for unknown plant was applied to the data from the three loops. In each case, the results are compared with the estimations given by the well-known iterative optimization procedure based on a two-dimensional grid search of values of S and J and identification of a Hammerstein model.8 The choice of this method for comparison is due to its capacity to estimate parameters S and J of the stiction, since various methods in literature estimates only the parameter S. Flow Loop I. Figure 16 shows the OP, PV, and set point (SP) signals of the loop; the measured amplitude of OP was 2.5 and the measured frequency17 was 0.0175 rad/s. Closed loop identification was used for identifying the PI controller. The first two-thirds of the data were used for identification and the one-third remaining was used for validation.

Figure 16. Flow loop I.

The ARX model used for controller identification was C(z) =

b1z + b0 z + a0

(37)

and the estimated controller transfer function was 14131


Industrial & Engineering Chemistry Research C(z) =

0.1479z − 0.1398 z−1

Article

(38)

The validation of the identified controller transfer function is shown in Figure 17. Clearly this PI controller is a good model. Other controller structures should be tried, if a0 ≠ −1 or the validation test does not succeed.

Figure 17. Validation of the PI controller.

The continuous time controller transfer function is given by ⎛ 1 ⎞ ⎟ C(s) = 0.148⎜1 + ⎝ 18.43s ⎠

(39)

The estimated gain of the process was K = 11.3. Using K and the controller parameters Kp = 0.148 and Ti = 18.43 in eq 29, yields 0.09 G (s ) = (40) s

Figure 18. Flow loop II.

The stiction parameters estimated by DFSQ method were SDF = 5.0 and JDF = 0.8. The grid search method using the search space {(S,J);0 ≤ J ≤ S, 0 ≤ S ≤ 10} with step 0.1 for both variables estimated the stiction parameters as Sgrid = 4.1 and Jgrid = 1.2. Flow Loop II. Figure 18 top shows the OP, PV, and SP signals of the flow loop II and Figure 18 bottom shows the OP × PV plot. The measurements of the amplitude and frequency of the OP signal were respectively A = 1.6 and ω = 0.15 rad/s. The transfer function of the identified PI controller was

that is not a regular oscillation. The algorithm to detect oscillations can easily detect the frequency. A high pass filter is then used to generate a regular oscillation so that the amplitude of the signal can be measured. The measurements of the amplitude and frequency of the OP signal were respectively A = 1.2 and ω = 0.006 rad/s. The identified controller transfer function was ⎛ 1 ⎞ ⎟ C(s) = 0.01271⎜1 + ⎝ 48.8s ⎠

⎛ 1 ⎞ ⎟ C(s) = 0.05146⎜1 + ⎝ 1.4118s ⎠

thus, the controller parameters are Kp = 0.01271 and Ti = 48.8. The estimated process gain was K = 30.8. The substitution of these parameters in eq 29 yields 0.008 G (s ) = s

The identification of the process gain yielded K = 7.4. Thus, using this gain and the controller parameters, the linear part of the flow loop II is given by G (s ) =

0.2695 s

The stiction parameters estimated by DFSQ method were SDF = 2.4 and JDF = 1.4. The grid method using the search space {(S,J);0 ≤ J ≤ S, 0 ≤ S ≤ 4} with step 0.1 for both variables estimated the stiction parameters as Sgrid = 2.2 and Jgrid = 1.3. Table 5 summarizes the results obtained using DFSQ and grid search methods. The calculation of S using eq 30 tends to produce overestimates of its values, since the amplitude of the OP signal tend to be greater than S/2 to overcame stiction.

The stiction parameters estimated by the DFSQ method were SDF = 3.2 and JDF = 1.4. The grid search method using the search space {(S,J);0 ≤ J ≤ S, 0 ≤ S ≤ 6.8} with step 0.1 for both variables estimated the stiction parameters as Sgrid = 2.7 and Jgrid = 1.6. Flow Loop III. As shown in Figure 19, this loop exhibits a variable set point. The variable set point produces an OP signal 14132



Article

the stiction. A procedure was proposed for situations where the process model is unknown. The conditions for existence and uniqueness of the solutions were provided. A sensitivity analysis was carried out to indicate situations where the error in the estimates can increase. A simple and efficient numerical algorithm was proposed to estimate the parameters even when the process model contains parametric uncertainties. The procedure could be easily extended to other nonlinearities that causes limit cycles, with the requirement that its describing function be available. The usefulness of this method was demonstrated through simulations, and its application to a pilot plant and to real data from industry. In the latter case, a comparison to a well-known method from literature was accomplished.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; [email protected]; [email protected]. Notes

The authors declare no competing financial interest.

■

(1) Srinivasan, R.; Rengaswamy, R.; Narasimhan, S.; Miller, R. Control loop performance assessment. 2. Hammerstein model approach diagnosis. Ind. Eng. Chem. Res. 2005, 44, 6719−6728. (2) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F.; Shook, D. S. Automatic detection and quantification of stiction in control valves. Control Eng. Pract. 2006, 14, 1395−1412. (3) Choudhury, M. A. A. S.; Jain, M.; Shah, S. L. Detection and quantification of valve stiction. Proc. Am. Control Conf. 2006, 2097− 2106. (4) Jelali, M. Estimation of valve stiction in control loops using separable least-squares and global search algorithms. J. Process Control 2007, 18, 632−642. (5) Hägglund, T. Automatic on-line estimation of backlash in control loops. J. Process Control 2007, 17, 489−499. (6) Araujo, A. P.; Munaro, C. J.; Filho, M. R. Quantification of valve stiction and dead band in control loops based on harmonic balance method. Proc. 9th IEEE ICCA 2011, 1066−1073. (7) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F.. A data-driven model for valve stiction. Proceedings IFAC Symposium on Advanced Control of Chemical Processes (ADCHEM). January 11−14, 2004, Hong Kong. (8) Choudhury, M. A. A. S.; Jain, M.; Shah, S. L. StictionDefinition, modelling, detection and quantification. J. Process Control 2008, 18, 232−243. (9) Engelberg, S. Limitations of the Describing Function for Limit Cycle Prediction. IEEE Trans. Autom. Control 2002, 47, 1887−1890. (10) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F. Diagnosis of Process Nonlinearities and Valve Stiction; Springer: New York, 2008. (11) Massey, W. S. Sufficient conditions for a local homeomorphism to be injective. Topol. Its Appl. 1992, 47, 133−148. (12) Bhattacharyya, S. P.; Chapellat, H.; Keel, L. H. Robust Control. The parametric approach; Prentice-Hall Inc.: Upper Saddle River, NJ, 1995. (13) Forssell, U.; Ljung, L. Closed-loop identification revisited. Automatica 1999, 35, 1215−1241. (14) Aström, K. J.; Hägglund, T. Advanced PID control; ISAThe Instrumentation, Systems, and Automation Society: Research Triangle Park, NC, 2006. (15) Rivera, D. E.; Morari, M.; Skogestad, S. Internal Model Control. 4. PID Controller Design. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 252− 265. (16) Ordys, A.; Uduehi, D.; Johnson, M. A. Process Control Performance Assessment (From Theory to Implementation); Springer, New York, 2007.

Figure 19. Flow loop III.

Table 5. Results of the Estimation of the Stiction Parameters loop

SDF

JDF

Sgrid

Jgrid

I II III

5 3.2 2.4

0.8 1.4 1.4

4.1 2.7 2.2

1.2 1.6 1.3

REFERENCES

The estimates of J in both methods are close. Since there is no information about the real values of S and J, it is not possible to conclude about the values that are the better estimates. The values obtained via grid search method produce the smaller mean squared error. However, the author’s experience have shown that the minimum error is quite sensitive to prefilters that must be used so that the two-parameter stiction model can work properly. The minimum error is also dependent on the order of ARMAX model used to identify the process. However, the computational effort for the DFSQ method is negligible when compared to the grid search method.

■

CONCLUSIONS A new method for the quantification of dead band and stiction in control loops is presented, based on the describing function method. The proposed method requires the amplitude and the period of the control signal and the transfer function of the controller and process to estimate the parameters of the DF of 14133



Article

(17) Karra, S.; Karim, M. N. Comprehensive methodology for detection and diagnosis. Control Eng. Practice 2009, 17, 939−956.

14134


Quantification of Valve Stiction and Dead Band in Control Loops

Recommend Documents