Tuning a Soft Sensor's Bias Update Term. 1. The Open-Loop Case

Article pubs.acs.org/IECR

Tuning a Soft Sensor’s Bias Update Term. 1. The Open-Loop Case Yuri A. W. Shardt* and Biao Huang* Chemical and Materials Engineering Department, University of Alberta, Edmonton, Alberta, Canada T6G 2 V4 ABSTRACT: The difficulty in measuring certain types of process variables rapidly has encouraged the use of soft sensors, which can determine the values of difficult to measure process variables based on easily available secondary process variables. A bias update term that allows the system to take into consideration disturbances in the system is often included in such soft sensor systems. The first part of this two-part series considers the bias update for the open-loop situation, including the ideal case, the case where there is measurement delay, the case with multirate sampling, and the case where there is a combination of measurement delay and multirate sampling. Proposed tuning rules are provided for all cases in order to obtain optimal open-loop tracking of the quality variable, especially in the presence of slow or drifting disturbances. Simulation and experimental validation of the results is presented.

1. INTRODUCTION In many industrial processes, certain key quality variables, such as concentration, cannot be measured sufficiently quickly or accurately enough to allow for making economic decisions or implementing closed-loop control. In both cases, a commonly used approach is to design a soft, or inferential, sensor to determine by using easier to obtain process data the value of the unknown quality variables.1 This basic soft sensor can be made more precise by incorporating the rarely available quality variable measurements from laboratory analysis through a bias update term.1,2 Three different types of soft sensors can be distinguished on the basis of the type of model used for forecasting the quality variable, namely, model-based soft sensors, inferential-based soft sensors, and identification-based soft sensors.3−5 Model-based soft sensors are developed on the basis of either first-principles or black-box modeling of the system.5 Inferential-based soft sensors incorporate into the model-based soft sensor information about secondary, or auxiliary, variables that are strongly correlated with the quality variable. Finally, identification-based soft sensors use various types of statistical methods, such as partial least-squares (PLS),6−8 principal component analysis (PCA),6 independent component analysis (ICA),9 and support vector machines (SVM),7 to determine appropriate relationships between the measured/known variables and the quality variables. More advanced models can be created by combining multiple models to deal with changes in operating points1,7 or to incorporate uncertainties in the system. Adding a bias update term to the soft sensor allows the incorporation of measured information about the quality variable. In many cases, it has been assumed that the bias update term is a simple constant correcting term, for which it has been observed that, in open-loop calibration, good results are obtained. Soft sensors can be run in either open loop or closed loop. Open-loop soft sensors are used in forecasting quality variables for economic analysis that is not used for feedback control,5,7,9 fault detection,5 or calibrating the soft sensor for closed-loop commissioning, while closed-loop soft sensors are used to forecast quality variables for control purposes.7 It should be noted that an open-loop soft sensor does not necessarily imply © 2012 American Chemical Society

that there is no feedback in the system. Rather, it implies that the soft sensor is used solely to provide the inferred value of the signal and that this value is not fed back to a controller. In industrial implementation, as well as during the simulation of soft sensors, when there are slow or drifting disturbances, it has been observed that by closing the loop, the performance of the soft sensor with a constant bias update term may deteriorate, even if its open-loop performance is satisfactory. This is especially the case if the bias update term is constructed so that it simply takes the difference between forecast and measured values and uses this difference as the new bias update. Thus, a detailed examination of this problem will be considered over this paper and the following paper in this issue using both simulations and experimental results with the objective of proposing tuning rules for the bias update term that can deal with slow or drifting disturbances. The first paper will consider the open-loop soft sensor, while the second paper will consider the closed-loop soft sensor. Therefore, the objectives of this first paper are (1) to theoretically analyze a soft sensor system with bias update in the open-loop case, (2) to theoretically analyze the system for different sampling rates and measurement delays, and (3) to propose appropriate tunings for the bias update terms in all cases. As well, validation of the proposed method using simulation and experimental results will be presented.

2. SOFT SENSORS WITH BIAS UPDATE 2.1. Problem Statement. The assumption is made that the open-loop soft sensor system can be described by Figure 1, where Gp is the “true” process transfer function, Ĝ p the assumed process transfer function, Gl the disturbance transfer function, and GB the bias update transfer function. Furthermore, the assumption is made that the measured input into the process can be given as ut, the disturbance by et, which is modeled as a white noise sequence with mean μ and variance Received: Revised: Accepted: Published: 4958

July 7, 2011 January 24, 2012 February 26, 2012 February 27, 2012 dx.doi.org/10.1021/ie201456z | Ind. Eng. Chem. Res. 2012, 51, 4958−4967

Industrial & Engineering Chemistry Research

Article

Figure 2. Schematic representation of an actual process.

in order to obtain a more appropriate comparison with the sampled process data, yt. It should be noted that the zero-order hold has a hold time of N samples. It is possible to rearrange this version into the one given by Figure 2. First, one should note that Ĝ B is a discrete filter. Next, it should be noted that both streams contain the same sampler and delay, so these two components can be placed after the subtraction. Finally, one should note that although the actual bias update transfer function that will be designed will change its form, placing the zero-order hold after Ĝ B makes simulations easier, since a single update value will be available for the given interval (Figure 3).

Figure 1. Schematic of an open-loop soft sensor with bias update.

σ2, the true output yt, and the soft sensor output ym,t. Finally, the assumption is made that yt is not directly measurable or cannot be measured fast enough, and, thus, there is a need to design a soft sensor to provide an estimate for the value of yt with sufficient sampling rate. However, the soft sensor output yα,t is based on the model Ĝ p, which unavoidably has error. This error is often corrected by adding the bias update term, yβ,t, which is based on the infrequently sampled yt. It should be noted that while initially the case where the soft sensor has the same inputs as the process is considered, the results obtained will be extended to the more general case, where the soft sensor can have different inputs as well as a different model from the process. From Figure 1, it can be seen that the following relationships hold

Figure 3. Detailed examination of the bias update.

yt = Gpu t + G le t

2.2. Open-Loop Investigation. Before considering the closed-loop case in the second part (the following paper in this issue), it is instructive to first consider the open-loop case and determine what conditions can be satisfied in order to obtain good tracking of the actual process output. Initially, it will be assumed unrealistically that GB = Ĝ B; that is, the data can be obtained instantaneously fast and whenever required. Strictly speaking, in this case, there is no need for a soft sensor. However, an analysis of this situation will provide interesting results that can be used in further cases. With comparison of yt and ym,t given in eq 2, it can be seen that, in order to match the two terms, it is necessary to have that at least

yα,t = Gp̂ u t yβ,t = GB(ym,t − yt ) ym,t = yα,t + yβ,t

(1)

Using eq 1, the open-loop transfer functions for ym,t and yt are written as yt = Gpu t + G le t

(2)

⎛ Ĝ − G G ⎞ −GBG l p B p⎟ ym,t = ⎜⎜ ⎟u t + 1 − G e t 1 G − B ⎠ B ⎝

(3)

GB = −1 1 − GB

2.1.1. Some Observations on GB, the Bias Update Term. In Figure 1, no assumptions have been made about the sampling rate at which the process is measured. The assumption is made that the controller update is every n samples. In a real process, there are two issues to consider: (1) yt is sampled every N samples, which is often much larger than the controller update period. This behavior can be modeled by adding a zero-order hold to the bias update term that has a hold of N samples. (2) A time delay in obtaining the data is present, which implies that although yt is sampled at time a, its value due to processing is not available until time a + d, where d is the time delay in processing the data. This behavior can be modeled by adding a d-sample time delay into the bias update term. Given the above discussion, the actual bias update system can be represented as shown in Figure 2. Note that a zero-order hold, which can be considered to be equivalent to a sampler, and a time delay term are added to the soft sensor output, ym,t,

(4)

Now, it can be noted that this equation is unsolvable in general. If we assume that Ĝ B = KB, that is, a constant gain, then eq 4 can only be achieved as KB → ±∞; that is, ⎛ KB ⎞ ⎜ ⎟ = −1 KB→±∞⎝ 1 − KB ⎠ lim

(5)

Under this situation, ym,t will approach yt, since the input term will become ⎛ Ĝ − K G ⎞ B p⎟ ⎜ p ⎜ 1−K ⎟ = Gp KB→±∞⎝ B ⎠ lim

(6)

This implies that, in order to obtain perfect tracking, the bias updating gain should be selected as large as possible. In addition to the above solution to the soft sensor tracking problem, it is also possible to consider the design of the soft 4959

dx.doi.org/10.1021/ie201456z | Ind. Eng. Chem. Res. 2012, 51, 4958−4967


Article

obtain a simple solution that both is stable and satisfies eq 10, we consider the case where

sensor on the basis of its performance in the long run, that is, as z → 1 or asymptotic convergence of the performance. The assumption is made that the bias update term can be written as nA ∑i = α i z −i G B = n 0 −j B βz ∑j= 0 j

α 0 = α1 = ... = αd − 1 = 0 αd = − 1 β0 = 1

(7)

β1 = β2 = ... = βd − 1 = 0 βd = −1

where α and β are arbitrary real numbers, such that β0 is not equal to zero. Then the following relationship should hold in order that the condition stated in eq 4 holds asymptotically, ⎛ G (z−1) ⎞ B ⎟ = −1 lim ⎜⎜ −1 ⎟ − 1 z → 1⎝ 1 − G B(z ) ⎠

that is, the bias update term, GB, which includes the effects due to the time delay, is written as

(8)

GB =

Substituting eq 7 into eq 8 and simplifying give n

n

B β − ∑A α ∑j= 0 j i=0 i

= −1 (9)

In order for this transfer function to be stable, the roots of the denominator (in terms of z) must lie inside the unit circle. Assuming that this has been satisfied, then eq 9 reduces to

β0 = 1

∑ βj = 0

(15)

1 − GB = 1 + KBz−d

(16)

will be roots = −KB1/ d

(17)

However, eq 4 requires that KB → ±∞, which implies that the roots as given by eq 17 will also go to infinity. Since they are outside the unit circle, this implies that the resulting forecast of the output value, given by ym, will be unstable. Despite the common industrial practice of a constant bias updating term, these results show that using a constant bias updating gain will not provide satisfactory tracking in the presence of time delay. Instead, the proposed tuning, given by eq 14, should be considered. 2.2.2. Multirate System in GB. For the multirate system, where the sampling rate for yt is much lower, there is an actual need for a soft sensor. This situation will now be analyzed. As was previously noted, adding an N-sample zero-order hold to GB is equivalent to saying that only information at a certain time period is available to update the bias term. Further the assumption is made that the only time at which the bias update term is changed is at kN, k ∈ . As well, this implies that the bias update term can only have data every N samples. Therefore, one way to represent Ĝ B, which excludes the effects due to the zero-order hold, at the sample instances kN is as Ĝ B(z−N), where Ĝ B is the original unmodified bias update transfer function; that is,

(11)

which will give a denominator that contains a single term of the form z−j, for which it is rather easy to verify that the denominator will not have any unstable roots. It can be noted that k represents an arbitrarily selected value, while l is all integers between 0 and max(nA,nB), not equal to k. 2.2.1. Time Delay in GB. As was previously noted, adding a time delay term to GB is equivalent to saying that the information about the true process takes time to measure. Once again, it is not strictly speaking necessary to use a soft sensor in this case as the obtained data can be simply shifted to obtain the process output. However, as for the ideal case, the results will be useful in analyzing more complex configurations of soft sensors later on. For simplicity the asymptotic analysis developed for the perfect case will be adopted here to simplify the analysis. A time delay of d samples will cause the first d terms of the numerator to be zero; that is, α 0 = α1 = ... = αd − 1 = 0

(14)

with all other values being equal to zero. The roots of the denominator given as

(10)

Equation 10 provides the primary constraint on attaining good asymptotic tracking between the measured and forecast variables. To obtain a solution that requires a minimal amount of computational effort in determining the roots, the following two additional conditions can be imposed on α and β ⎧ ⎪ αk − βk ≠ 0 ∃ k such that 0< k < max(nA ,nB) ⎨ ⎪ ⎩ αl − βl = 0 ∀ l ≠ k

1 − z −d

αd = −KB

nB

j=0

− z −d

It is easy to verify that eq 10 is satisfied and that the auxiliary condition for stability given by eq 11 is also verified, since only the first term (j = 0) is nonzero, with all other values being zero. Furthermore, for a constant bias updating gain, the following situation will be obtained,

n

A α ∑i = 0 i

(13)

n

Ĝ B =

(12)

A α z−iN ∑i = 0 i

n

∑i =B 0 βiz−iN

(18)

This will force all of the values of α and β between kN and (k + 1)N to be zero. Since, over this interval, where there will be no update from the system, the last available value will be

This set of zero terms will introduce issues with determining the stability of the system. However, the key result given as eq 10 from the perfect case will still apply to this situation. To 4960



Article

used, the difference between the measured and forecast values is denoted as yd; that is, yd = ym − y (19)

Finally, eq 28 can be rewritten to give ym̂ = (IN − B)̂ −1Gp̂ IN u ̂ − (IN − B)̂ −1Bŷ ̂

It should be noted that since IN − B̂ is lower triangular with nonzero diagonal terms, its determinant cannot be zero, and hence an inverse can always be found. Furthermore, the inverse can be obtained as follows:

where ym is the forecast process output and y is the actual process output. At any time interval, t, the following equation describes the behavior of the updating term, ym,t = Gp̂ u t + GB(z−N )yd ,t

(29)

(20)

where ym,t represents the value of ym at time instance t and yd,t represents the value of yd at time instance t. At the time point t = N, eq 20 can be rewritten as ym, N = Gp̂ uN + GB(z−N )yd , N

(21)

This can be obtained by applying the standard manual inversion method by forming the matrix [A|Im] and then row reducing it so that the result [Im|A−1] is obtained.10 This allows eq 27 to simplify to

while at t = N + 1, eq 20 can be rewritten as ym, N + 1 = Gp̂ uN + 1 + GB(z−N )yd , N + 1

(22)

However, it should be noted that, due to the fact that the data are only available every N samples, at this time instance, yd,N+1 is not available. Instead, it should be replaced by the last available value, yd,N; that is, eq 22 can be rewritten as ym, N + 1 = Gp̂ uN + 1 + GB(z−N )yd , N

(23)

Furthermore, using eq 19, eq 23 can be rewritten to give ym, N + 1 = Gp̂ uN + 1 + GB(z−N )(ym, N − yN )

(24)

A similar argument can be made for all of the remaining time intervals between N and 2N, for which no new update values are available. Therefore, the equations for the system can be written as for each point in the time interval [N, 2N[ as Using eq 2, ŷ can be written in terms of the input and the disturbance, to give for the N-sample time interval,

At the time point 2N, a new bias update term is available and the cycle restarts for the new time interval. Thus, we can write the following system of equations to describe the sampled data system

Combining eqs 31 and 32 gives ym̂ = (Ĝ − Ĉ(GpIN ))u ̂ − Ĉ(G lIN )e ̂

(33)

Taking the first row of the matrices in eq 33 gives the equation for ym,N at the sample interval N at which information about the system is available; that is, ⎛ Ĝ − Ĝ (z−N )G ⎞ −Ĝ B(z−N )G l p B p⎟ eN u ym, N = ⎜ + ⎜ 1 − Ĝ (z−N ) ⎟ N 1 − Ĝ B(z−N ) B ⎝ ⎠

where IN is the N × N identity matrix, B̂ is an N × N matrix, and the sampling time of this new system is N. Using the resulting matrix equation and rewriting it gives ym̂ = Gp̂ IN u ̂ + B(̂ ym̂ − y )̂

(34)

(27)

Since eq 34 has the same form as eq 3, the results obtained for the perfect case hold, except that GB will have the form given by eq 18. Specifically, this shows that a constant bias update term will be sufficient for the open-loop case with a zero-order hold

Rearranging eq 27 gives (IN − B)̂ ym̂ = Gp̂ IN u ̂ − Bŷ ̂

(28) 4961



Article

Figure 4. Forecast process values plotted against the actual process values for (left) case 1 and (right) case 2. The dashed line represents the y = x line.

variables have a transfer function given as Ga with values given by dt. Thus, the true process output can be written as

in the system. As well, solutions similar to that of the time delay case will also hold asymptotically; that is,

yt = Gput + Gad t + G le t

α0 = −1 β0 = 1

The soft sensor forecast values can then be written as

βN = −1

(35)

ym,t =

with all other terms equal to zero; that is,

Ĝ B =

(37)

(38)

Comparing eqs 37 and 38 against the original open-loop case given by eq 2, it is easy to see that the previously developed results will hold in this case, since the term due to the secondary variables has the same general form as the input variable and hence will behave the same way as before. This implies that although adding auxiliary variables may improve the accuracy of the soft sensor and, hence, reduce the magnitude of the bias update term, the general form of the bias update term will remain the same as before.

−1 1 − z −N

Gp̂ − GBGp Ĝ − GBGa − GBG l ut + a dt + et 1 − GB 1 − GB 1 − GB

(36)

It can be noted that the only nonzero term is at j = N, which satisfies the constraints given in eq 11. It should be noted that the best tracking is obtained with a large constant bias updating term. 2.2.3. Time Delay and Multirate Sampling in GB. The most realistic situation for a soft sensor is to combine a time delay with a zero-order hold; that is, it takes time to obtain the data which are sampled at a rate that is often much slower than required for forecasting. It should be noted that the time delay will not be affected by the zero-order hold, since the time delay is not a component of Ĝ B which is affected by the zero-order hold. If the time delay, d, is equal to the sampling interval, N, then the results will be very similar to that obtained for a time delay case, since GB will have a form similar to that previously obtained; that is, eqs 10 and 11 can be used to obtain the desired solution. Specifically, one should note that a constant bias update term will not suffice due to the presence of time delay. If d < N, then experimentally, it has been verified that using Ĝ B defined on the basis of the zero-order hold case and ignoring the time delay will provide the correct results. On the other hand, the solution for the case where d > N is still an open problem. However, if the time delay is an integer multiple of the zero-order hold, then the maximum of either d or N is taken and used when constructing the bias update function. This result can be derived by combining the results from the multirate derivation with the time delay component. 2.2.4. Augmenting the Estimated Plant Model. In most soft sensor applications, it is possible that additional, secondary variables may be measurable, whose impact on the process could be quantified. The assumption is made that these secondary

3. SIMULATIONS The above open-loop results will be tested on two different types of simulations: a pure discrete time system and a sampled data system. 3.1. Discrete Time Simulations. The discrete time simulations will be performed in Simulink using the following models and parameters: (1) Gp = 0.5z−1/(1 − 0.85z−1); Ĝ p will take on different values depending on the simulation. (2) Gl = (1 + 0.5z−1)/(1 − z−1). (3) The bias update term will be designed differently depending on the circumstances. (4) White noise with a mean of 0 and a variance of 1 will be used. (5) The Monte Carlo simulation was run for 10 000 time samples. 3.1.1. Ideal Case. The assumption is made in this section that the bias update function is given as GB = KB and that Ĝ p = Gp; that is, there is no plant−model mismatch. Two cases are considered: (1) KB = −1 and (2) KB = −20. Negative values of KB were selected since these provide a stable closed-loop system, and it is desired to compare the open-loop results with the closed-loop ones later on. Changing the sign of KB flips the graph about the y = x line. Figure 4 shows the forecast process values plotted against the actual process values for both cases. Ideally, the data should 4962



Article

cluster about the y = x line. From Figure 4, it can be seen that for case 1 the data are quite far from the desired location and that consequently there is offset bias. On the other hand, for case 2, where KB = −20, the data are much closer to the desired line and the offset bias is correspondingly smaller. This suggests that, although theoretically the bias correction gain should be infinite, the actual bias gain can be quite small and still obtain reasonable tracking of the true process values. 3.1.2. Time Delay Case. The simulation results are discussed for the same case as before except that there is a time delay of d = {1, 2, 4, 6, 8, 10, 12} and GB =

− z −d 1 − z −d

(39)

which satisfies the constraints for asymptotic tracking previously developed. No plant−model mismatch was assumed; that is, Ĝ p = Gp is now considered. The accuracy of the forecast was determined using the average forecast error (AFE) defined as11 AFE =

1 15000

Figure 6. Average forecast error (AFE) as a function of the zero-order hold, N.

15000

∑

|yi − ym, i |

Figure 5 shows the AFE as a function of d. As expected, as the time delay increases, the AFE increases. This can be attrib-

is a greater discrepancy between the forecast and actual values. As well, it can be noted that the effect of the zero-order hold on the accuracy of the obtained results is smaller than for the case of the time delay. Figure 7 shows the forecast values as a function of the actual process values in the presence of a drifting disturbance and a

Figure 5. Average forecast error (AFE) as a function of the measurement time delay, d.

Figure 7. Simulation results with 5-sample zero-order hold. The dashed line represents the y = x line.

uted to the fact that as the time delay increases, the relevancy of the bias update value to the current value decreases and there is a greater discrepancy between the forecast and actual values. 3.1.3. Multirate System. The simulation results are considered for the same case as before except that there is a zeroorder hold with sampling time of N = {1, 2, 4, 6, 8, 10, 12}, and

five-sample zero-order hold. This confirms that the forecast values are similar to the actual values. Increasing the value of N will not affect the tracking behavior of the soft sensor; for example, for N = 480 samples, the AFE is 16.5. 3.1.4. Both Multirate Sampling and Time Delays. Figure 8 shows simulation results for the same case as before except that there is both a time delay of d = 3, a zero-order hold with a sampling time of N = 5, and the proposed solution was used; that is,

GB =

i=1

(40)

−1 1 − z −N

(41)

which also satisfies all of the constraints for good tracking. Figure 6 shows the AFE as a function of the zero-order hold, N. As expected, as the zero-order hold period increases, which is equivalent to decreasing the sampling rate of the variable of interest, the AFE increases. This can be attributed to the fact that as the zero-order hold increases, the amount of information obtained from the process decreases in relevancy and there

GB =

−1 1 − z −5

(42)

From Figure 8, it can be seen that the proposed solution provides tracking of the output. The AFE for this case is 1.80. 3.2. Soft Sensor for a Continuous, Stirred Tank Reactor: A Sampled Data System. To investigate the results of the 4963



Article

concentration of component A, CA, whose sampling may take some time. A secondary variable was assumed to be the temperature of the reactor, which could be sampled much faster. The soft sensor was constructed by linearizing the differential equations about the nominal values. Two different models were obtained:

proposed system in a continuous time environment, the standard continuous, stirred tank reactor (CSTR) model for an irreversible, exothermic reaction given as A→B, as originally proposed by Morningred et al.,12 was used. The governing differential equations can be written as dCA ⎛ −E ⎞ V̇ ⎟ = (C A 0 − CA ) − k 0CA exp⎜ ⎝ RT ⎠ dt V

(44)

CA(s) =

⎛ −0.00266 ⎞ ⎜ ⎟ T (s ) ⎝ 0.979s + 1 ⎠

(45)

Gl =

( −ΔH )k 0CA ⎛ −E ⎞ dT V̇ ⎟ = (T0 − T ) − exp⎜ ⎝ RT ⎠ ρcp dt V ⎞⎞ ⎛ ρcc p ⎛ −hA ⎟⎟ c ̇⎜ + Vc 1 − exp⎜⎜ (Tc − T ) ρcpV ⎜ V̇ ρcc p ⎟⎠⎟ ⎝ ⎝ c ⎠

(43)

Ĝ B =

nominal values

CA T V̇ c V̇ CA0

0.1 mol/L 438.54 K 103.41 L/min 100 L/min 1 mol/L

feed temperature inlet coolant temperature tank volume heat-transfer term reaction rate constant activation energy term enthalpy of reaction liquid density constant pressure heat capacity

T0 Tc V hA k0 E/R ΔH ρ, ρc cp, cpc

350 K 350 K 100 L 7 × 103 cal/min/K 7.2 × 1018 min−1 1 × 104 K −2 × 103 cal/min 1 × 103 g/L 1 cal/g/K

(46)

−1 1 − z −1

(47)

It was assumed that there were no set point changes. The set point case will be considered in the experimental results later. The results are shown in Figure 10. The few points that are located far from the main cluster result from the initialization of the continuous and discrete time components of the model. Using eq 44 as the model, the AFE is 2.9 × 10−5, while, with eq 45 as the model, the AFE is 2.0 × 10−5. Discretizing eq 44 does not change the AFE. For comparison, using a constant bias update term of −0.5 and eq 44 as the model gives an AFE of 0.003, which is about 100 times larger than for the integratorbased solution. This suggests that, irrespective of the model selected, using the proposed solution of an integrator in the bias update term provides good tracking in the presence of a slow disturbance. It should be noted that in the open-loop case the above proposed solution provides good tracking even if the disturbance variable contains an integrator. A sample result is shown in Figure 11, where the disturbance transfer function was

Table 1. Nominal Parameter Values outlet product concentration reactor temperature coolant flow rate process flow rate feed concentration

0.1 100s + 1

The white noise input was assumed to have zero mean and a variance of 0.005. It was assumed that there was a 1 min time delay between sampling the data and obtaining the value. In all cases, the following form was assumed for the bias update term

where CA is the concentration of component A, V̇ the volumetric flow rate, V the volume of the reactor, k0 the reaction rate constant, E the activation energy, R the universal gas constant, T the temperature of the reactor, ΔH the enthalpy of reaction, ρ the density of the fluid in the tank, cp the constant pressure heat capacity, h the heat-transfer coefficient, and A the crosssectional area of the heat exchangers; the subscript 0 refers to the inlet conditions, and the subscript c refers to the coolant. The nominal parameter values are given in Table 1.

symbol

⎛ −0.00266 ⎞⎛ 0.8775 ⎞ ⎟ V ̇ (s ) ⎜ ⎟⎜ ⎝ 0.979s + 1 ⎠⎝ 7.324s + 1 ⎠ c

The soft sensor defined by eq 44 is a model based on the inputs to the process, while the soft sensor defined by eq 45 is a model based on the fast sampled output (secondary variable) from the process, the temperature. As well, eq 44 was discretized using the exact discretization formulas13 and assuming a sampling time of 1 min. These three models were used to test the effects on the performance of the bias update tuning for different types of soft sensor designs. The Simulink system used to simulate the process is shown in Figure 9. The disturbance transfer function was defined as

Figure 8. Simulation results with three-sample time delay, five-sample zero-order hold. The dashed line represents the y = x line.

parameter name

CA(s) =

Gl =

0.01 s

(48)

and the soft sensor was defined by eq 44. As expected, the tracking is good.

4. EXPERIMENTAL RESULTS 4.1. Setup. To test the above results, a heated tank system with a schematic given in Figure 12 was used.

The coolant flow rate was chosen as the manipulated variable, while the inlet concentration was chosen as the disturbance variable. The quality or controlled variable was chosen as the 4964



Article

Figure 9. Simulink diagram of the system.

Figure 10. (top left) Using the soft sensor defined by eq 44; (top right) using the soft sensor defined by the discretization of eq 44; and (bottom) using the soft sensor defined by eq 45. The dashed line is the y = x line. 4965



Article

Figure 11. Tracking behavior in the presence of an integrating disturbance and a soft sensor defined by eq 44. The dashed line is the y = x line.

Figure 13. Experimental data.

Figure 12. Process schematic.

For the heated tank system, the following parameters were used: (1) The level in tank was set to 0.2 m and controlled using the cold water flow rate cascade controller. An integrating disturbance in the level driven by white noise with mean zero and variance 0.000 01 was introduced into the system. (2) The following nominal steady-state values were used: the cold water flow rate was 5.8 kg/min, the temperature was 32.2 °C, and the steam flow rate was 15 kg/h. (3) The measured input for the soft sensor is the steam flow rate. The temperature from the system can only be obtained every 3 s. Controller update is every second. (4) On the basis of the previous experiments with this system,14 the continuous-time model between the steam and the temperature can be written as

Gp =

1.2 e−30s 66s + 1

Figure 14. Scatter plot of the measured and forecast temperatures. The dashed line is the y = x line.

steam changes and when the steam stays constant. As well, it can handle the transition between the two set points.

5. CONCLUSIONS In this first part of the two-part series, the tuning relationships for the bias update term have been developed for the ideal case, a measurement time delay, a multirate system, and a combination of both. It was shown that the type of the model used has an insignificant influence on the tracking behavior of the system. Simulations and experimental results confirmed the proposed tuning rules for the bias update term. The proposed tuning rules for the bias update term can be summarized as follows: (1) In the absence of time delay or zero-order hold, then a large constant term can be used to obtain optimal tracking of the process. (2) In the presence of time delay, then the optimal bias update term is given as

(49)

which was discretized using the exact discretization method13 and a sampling time of 1 s. This model was determined using a simple step test, so that it is an approximate model of the system, which may not be too accurate. (5) The data are sampled every second. The experiment was run for 1800 s. A step change of 5 kg/h in the steam was introduced at 900 s. (6) The bias update term was designed on the basis of the proposed methods. 4.2. Results. The data obtained are shown in Figure 13, while a scatter plot of the forecast and actual temperatures is shown in Figure 14. It can be seen that the forecast temperature is accurately predicted by the soft sensor for both when the

Ĝ B =

−1 1 − z −d

(50)

where d is the time delay. 4966



Article

(12) Monringred, J. D.; Paden, B. E.; Seborg, D. E.; Mellichamp, D. A. An adaptive nonlinear predictive controller. Chem. Eng. Sci. 1992, 47 (4), 755−762. (13) Huang, B.; Kadali, R. Dynamic Modeling, Predictive Control and Performance Monitoring: A Data-Driven Subspace Approach; SpringerVerlag: London, England, U.K., 2008. (14) Shardt, Y. A. W.; Huang, B. Closed-loop identification with routine operating data: effect of time delay and sampling time. J. Process Control 2011, 21 (7), 997−1010.

(3) In the presence of a zero-order hold, that is, the process data is only sampled every N samples, then the optimal bias update term is given as ĜB =

−1 1 − z −N

(51)

(4) In the presence of both time delay and a zero-order hold, the following three cases need to be distinguished: (a) If d ≤ N, then the optimal bias update term is given as ĜB =

−1 1 − z −N

(52)

(b) If d is an integer multiple of N, then the optimal bias update term is given as Ĝ B =

−1 1 − z −d

(53)

(c) The general result for the case where d ≥ N is still an open question. Future work will focus on solving the time delay and zeroorder hold for the case when both are present in the system, especially for the case where the time delay is larger than the sampling time, to determine appropriate tuning rules.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: 780-492-9016. E-mail: [email protected] (Y.A.W.S); biao.huang@ ualberta.ca (B.H.). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Domlan, E.; Huang, B.; Xu, F.; Espejo, A. A decoupled multiple model approach for soft sensor design. Control Eng. Pract. 2011, 19, 126−134. (2) Shao, X.; Huang, B.; Lee, J. M.; Xu, F.; Espejo, A. Bayesian method for multirate data synthesis and model calibration. AIChE J. 2011, in press. (3) Sahebsara, M.; Chen, T.; Shah, S. L. Optimal fast-rate soft-sensor design for multi-rate processes. Proceedings of the 2006 American Control Conference, Minneapolis, MN, USA; IEEE: Washington, DC, 2006; pp 976−981. (4) Fortuna, L.; Graziani, S.; Rizzo, A.; Xibilia, M. G. Soft Sensors for Monitoring and Control of Industrial Processes; Springer-Verlag: London, England, U.K., 2007. (5) Kadlec, P.; Gabrys, B.; Strandt, S. Data-driven soft sensors in the process industry. Comput. Chem. Eng. 2009, 33, 798−814. (6) Flynn, D.; Ritchie, J.; Cregan, M. Data mining techniques applied to power plant performance monitoring. DYCOPS 2005, Prague, Czech Republic; Elsevier: 2005. (7) Kadlec, P.; Grbić, R.; Gabrys, B. Review of adaptation mechanism for data-driven soft sensors. Comput. Chem. Eng. 2011, 35, 1−24. (8) Aumi, S.; Mhaskar, P. Integrating data-based modeling and nonlinear control tools for batch process control. AIChE J. 2011, in press. (9) Lee, J.-M.; Yoo, C.; Lee, I.-B. Statistical process monitoring with independent component analysis. J. Process Control 2004, 14, 467− 485. (10) Anton, H. Elementary Linear Algebra, 8th ed.; John Wiley & Sons: Hoboken, NJ, USA, 2000. (11) Liu, J.; Chen, D.-S.; Shen, J.-F. Development of self-validating soft sensors using fast moving windows least squares. Ind. Eng. Chem. Res. 2010, 49, 11530−11546. 4967


Tuning a Soft Sensor's Bias Update Term. 1. The Open-Loop Case

Recommend Documents