Optimization of temperature programmed sensing for gas identification using micro-hotplate sensors

Sensors and Actuators B 53 (1998) 24 – 43

Optimization of temperature programmed sensing for gas identification using micro-hotplate sensors Tekin A. Kunt 1,a, Thomas J. McAvoy a,*, Richard E. Cavicchi b, Steve Semancik b a

Institute for Systems Research, Department of Chemical Engineering, Uni6ersity of Maryland at College Park, Chemical Engineering & The Institute for Systems Research, A.V. Williams Building 115, College Park, MD 20742, USA b National Institute of Science and Technology, Chemical Science and Technology Laboratory, Gaithersburg, MD 20899, USA Received 23 December 1997; received in revised form 14 July 1998; accepted 15 July 1998

Abstract Micro-hotplate chemical gas sensors, such as those being developed at the National Institute of Standards and Technology (NIST) by micromachining Si, can be operated in a temperature-pulsed mode, due to their small size and mass. In the temperature-pulsed mode of operation, different gases give different dynamic responses (i.e. signatures) depending on the temperature program used. In this paper a new methodology is presented to optimize the operation of micro-hotplate gas sensors for discriminating volatile organic compounds at a fixed concentration, while minimizing the detection time. The extension of the methodology to cases where concentrations vary is currently under investigation. The Wavelet Network method is applied to accurately predict the sensor’s response for a given temperature profile. Once a dynamic model is obtained, it is used for off-line optimization of the temperature profile, i.e. the maximization of the difference between two gas signatures. The difference between two response curves was initially measured by a metric based on the Euclidean distance. This metric was then modified using the Haar wavelet transformation. The methodology was implemented in a case study in which either methanol or ethanol had to be detected in air, but the methodology is generic, and it can be applied to any two gases. © 1998 Elsevier Science S.A. All rights reserved. Keywords: Optimization; Programmed sensing; Micro-hotplate sensors

1. Introduction Although Taguchi sensors have been available since the 1970s, the idea of using temperature modulation on a single sensor to selectively detect gases is relatively new and unexplored. Since the early l990s, a few investigators have proposed the thermal cycling of a chemical sensor as a promising technique for selective detection of flammable and toxic gases in air [1]. There exist, however, different paradigms attempting to explain conductance changes in SnO2 based chemical sensors when the temperature is cycled continuously. There are also different methods to process the resulting conductance versus time signatures.

* Corresponding author. Tel:. + 1 301 405 6632; fax: +1 301 314 9920; e-mail [email protected],[email protected] 1 Aspen Technology Inc., 14701 St. Marys Lane, Houston, Tx 77079-2995, USA

Clifford [2,3] appears to be the first researcher to investigate the dynamic change of sensor conductance for individual temperature steps at several different oxygen partial pressures. He used the barrier layer theory [4] to explain the effect of temperature on the sensor response. According to the barrier layer theory, two factors, bulk conductance and the surface concentration of ionosorbed oxygen, change the macroscopic conductance of an SnO2 film. Both factors can be modeled as functions of temperature and time at a constant oxygen pressure using dynamic conductance measurements. Clifford concluded that selective gas detection could be accomplished using either many sensors at different but fixed temperatures, or by sequential operation of a single sensor at several temperatures. Long time constants, however, would complicate measurements made with Taguchi gas sensors in the latter case. Clifford also provided a detection model for which much empirical testing and measurements were required to compute the related parameters.

0925-4005/98/$ - see front matter © 1998 Elsevier Science S.A. All rights reserved. PII S0925-4005(98)00244-5

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

An investigation by Wlodek and co-workers [5] assumed that the chemical reaction between surface oxygen and a reducing gas was the primary mechanism of tin oxide sensing. This formulation required using many empirical parameters. Although the approach increased the general understanding of gas-surface reaction kinetics, it did not provide predictive models from which one could infer conductance versus time signatures. Nakata and co-workers [6] realized the potential of using multidimensional information obtained by applying a sinusoidal temperature change to a semiconducting metal oxide (SnO2) sensor. They analyzed the resulting dynamic nonlinear conductance response using fast Fourier transform (FFT) [7] methods. According to their analysis, higher harmonics in the FFT signal had a physicochemical significance for discriminating between hydrocarbons. The emphasis, however, was again on processing of the resulting signals rather than finding an accurate mapping between the temperature and conductance response. Moreover, they did not investigate the effect of the modulation frequency and the temperature range used in their experiments. Previous gas identification experiments with a single sensor had assumed a particular temperature modulation, either sinusoidal with a fixed frequency [6,8], or cyclic pulsed heating such as linear ramps [9]. In those examples, the speed of operation was about 20 – 25 s per cycle, with rates limited by the sensing devices employed. Advanced microsensor design, however, allows for faster implementation. In this paper, we propose a methodology to optimize the operational modes of microsensors for performing selective and fast detection of similar chemicals. Dynamic modeling and optimization of micro-hotplate sensors are presented after describing temperature programmed sensing and its advantages. This methodology is illustrated for the detection of methanol or ethanol in air, but the methodology is generic and can be applied to any two gases at a fixed concentration. Although this experimental test can be criticized as being somewhat narrow, two points can be noted. First, the approach discussed can be extended to cases where gas concentration varies, and this extension is currently under study. Second, the primary objective of the present paper is to show how data based dynamic modeling and optimization techniques can be used to take advantage of the temperature programmability of the microhotplate sensors studied. The results presented show that these techniques result in significantly improved sensor responses for the problem discussed.

25

Researchers at the National Institute of Standards and Technology (NIST) designed and fabricated a gas sensor consisting of an array of micro-hotplates on a silicon wafer, using commercially available complementary metal oxide semiconductor (CMOS) technology [10,11]. Each individual sensor contains three electrically conducting, functional layers that are electrically insulated from other layers: a polysilicon resistor for heating, an aluminum hotplate for distributing heat and measuring the temperature, and four aluminum contact pads for connection to the sensing film (see Fig. 1). SnO2 is used as the sensing film and is deposited onto the heated surface of the micro-hotplate in an ultrahigh vacuum chamber. The area of each sensor is around 104 mm2. This miniaturized device has a very fast thermal time constant, which is estimated to be 1 ms [12], and a thermal efficiency of about 8°C mW − 1 [10]. The fast thermal response time allows rapid temperature programmed sensing. Rapid pulsing of the sensor temperature is also possible, and is employed in the programming to separate thermally controlled chemical

Fig. 1. A micro-hotplate sensor (see text for details).

2. Micro-hotplate sensors Recent advances in micromachining have enabled development of a new generation of chemical sensors.

Fig. 2. Temperature programmed sensing: (a) temperature pulse amplitude (20 – 450°C); (b) pulse duration (10 – 300 ms); (c) delay (5 ms); (d) number of pulses in a cycle (2 – 64).

26

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 3. (a) The conductance response in methanol and ethanol for sensor screening tests; (b) the corresponding temperature pulses.

effects from temperature dependent changes in the electrical characteristics of the sensing films. In the pulse-programmed mode of operation, the sensing surface is pulsed from a base temperature to variable elevated temperatures, instead of keeping the temperature constant. The sensor’s conductance is always measured at a reference temperature (usually room temperature) between the temperature pulses so that only the effects of chemical factors are recorded. Early studies using micro-hotplate chemical gas sensors showed that different gases provided characteristic signatures [9,13]. Although chemical species from various subgroups such as formaldehyde and methanol were observed to have quite different signatures, two alcohols such as methanol and ethanol produced similar traces for a linear ramp temperature program. The results demonstrated that temperature programmed sensing was a very powerful method for generating useful information using a single sensor. While the response patterns could be altered by changing the temperature pulse sequence, a systematic approach to optimize the sequence was not presented [9]. While micro-hotplate sensor devices are typically fabricated in arrays, the work presented here involves only experiments on single micro-hotplate gas sensors. The concepts described in this paper could be applied to multiple elements. Note that the information content for a single, temperature-programmed microsensor might be higher than that for multiple sensors operating at constant temperatures.

3. Methodology This paper aims to use systems engineering tools to answer the following question: ‘armed with a micro-hotplate sensor that can be quickly thermally-switched, how does one select the best temperature scheduling for fast and selective identification of chemical species?’ Temperature scheduling of a microsensor is defined by the temperature profile (pulse amplitudes), the number of pulses in this sequence, the duration of each temperature pulse, and the delay before measuring the conductance at room temperature (see Fig. 2). Temperature transient measurements have shown that waiting 5 ms before the conductance is measured allows thermal equilibrium to be attained [12]. While we have examined variation of the delay time to larger times in our overall studies, we do not report those results here. Conductance is always measured at room temperature to observe only conductance changes associated with interfacial chemical effects, and eliminate any temperature-dependent change in the inherent conductivity of the sensing film. In order to solve the optimization problem, a two step procedure is introduced: (1) dynamic modeling, and (2) off-line optimization. The first step maps the sensor’s operating condition (amplitudes of the temperature pulses) to the sensor’s gas-specific dynamic response. Off-line optimization finds the ‘best’ and the ‘fastest’ mode of operation using empirical models when detecting similar gases in air (one should note that a trade-off may exist between the ‘best’ and the

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

‘fastest’ mode of operation). The optimization problem should be formulated and solved using an appropriate distance metric, set of constraints, and initial parameters.

27

oped to assess the reliability of the empirical models. These tests can be applied as constraints when optimizing the sensor’s operation.

4.1. Sensor characteristics 4. Dynamic modeling of micro-hotplate gas sensors To date, very little use has been made of the thermal transient information in a sensor’s response [14]. This has been partly due to the slow time constant inherent in thick film and relatively large metal oxide semiconductor gas sensors. On the other hand, the micro-hotplate conductometric gas sensors developed at NIST can be heated/cooled from room temperature to 500°C in a few milliseconds [9]. Dynamic modeling techniques are required to fully take advantage of the fast response characteristics of these micro-hotpate gas sensors. The methodology developed for empirical modeling of microsensors consists of six basic steps that are arranged into two main groups: experimental and computational aspects. Experimental considerations involve: (1) the transient response measurements, (2) the sensor characteristics, and (3) the selection of the training sequence, followed by computational steps: (4) the model order selection, (5) the functional mapping, and (6) the model validation. Transient response measurements help define the important variables in operating the micro-hotplate gas sensors. Once these variables are fixed, step-like sensor screening experiments are performed to observe the short and long term responses of the sensing film to a specific gas. Finally, the temperature sequence is generated using a continuity requirement which limits the maximum temperature difference between two successive temperature pulses. The resulting temperature sequence is applied to obtain the necessary input – output database for synthesizing empirical models. Since the sensor’s conductance yi at time i will depend on the previous conductance values as well as the temperature history, ui, the next conductance value can be described as yi + 1 =F[yi, yi − 1,..., yi − (ny − 1), ui + 1, ui,..., ui − (nu − 2)]

(1)

where nu and ny denote maximum time lags (model orders) in the input and output, respectively. The model order involves the number of past u and y values, i.e. ui, ui − 1, yi, yi − 1,etc. that are needed to predict yy + 1 accurately. Model orders can be computed using the false nearest neighbors approach [15,16]. Several methods have been compared for obtaining the function F. The Wavelet Network technique [17] has thus far provided the best results in accurately predicting both short and long term response of the sensor [18]. Model validation and robustness tests have also been devel-

Quick initial assessment of the sensing capabilities for a new sensor, and periodic and rapid checks of an already operational sensor, can be important to ensure proper performance. For these purposes, sensor characteristics can be obtained using some basic, fast and exploratory experiments. A linear temperature ramp, i.e. temperature pulses starting from room temperature and rising up to the maximum allowable temperature with fixed increments, is one approach that has been applied for scanning device performance [9]. Another possible temperature modulation is a step-like sensor screening experiment developed in this work. In this mode, both the slow and fast dynamics are revealed by repeating the same temperature pulse amplitude many times and cleaning the surface with high temperature pulses between different temperatures. An example of such a step-like temperature profile with cleaning pulses in between is shown in Fig. 3. In this example, there are 200 pulses at each temperature. After 200 pulses at a particular temperature, 200 high temperature pulses at 425°C (maximum allowable temperature for that particular sensor) are applied to initialize the surface for the next temperature pulses. In this scheme, seven different temperatures are investigated starting from 50 to 350°C with 50°C increments. Two characteristics are striking in Fig. 3: (1) the transients for conductance responses are a function of temperature, (2) the response behavior for methanol and ethanol is dynamically different for pulses around 300°C, although it is similar at other temperatures. While the conductance in methanol gas goes down with temperature pulses, the response in ethanol goes up for 300°C pulsing. The different dynamics are an important factor when the temperature profile is being optimized for the separation of similar compounds, such as methanol and ethanol. Indeed, as will be seen below, the optimization algorithm tends to choose a temperature region around 300°C for the best separation of methanol and ethanol. The aim of these experiments was to see how the film conductance responds when the sensor is subject to the same temperature pulses many times, so that the response should go to some steady state value in the limit. In the analysis, however, it is not the actual steady state value, but the way that the conductance converges to that value, that is more important, since a dynamic temperature program will eventually be used in sensing applications. A similar approach for understanding the effect of temperature on the sensing mechanisms was developed by Clifford [2]. He investigated

28

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 4. Two successive cycles (solid and dashed lines) for the sensor response in 80 ppm hydrogen gas.

the temperature stimulated kinetic response of Taguchi sensors by applying a special sequence of heater voltages and measuring the conductance for those larger, slower thermally-controlled devices. Using dynamic measurements, he observed that the initial rapid change in conductance was followed by a slow relaxation towards its final steady state value. He also deduced that the temperature dependent exponential decay times, which are on the order of 104 s, were responsible for the hysteresis and long stabilization times often associated with sensor responses. After observing different response characteristics during the step-like sensor screening experiments and trying to develop models from the resulting data, a continuity requirement in the temperature profile was adopted to prevent sudden changes in the surface condition. When large step changes in temperature are used, these changes disrupt the sensor surface in a manner that results in data that cannot be modeled accurately. As a result when training sequences are designed, a constraint is included to assure that the difference between two successive temperature pulses is not greater than a fixed value such as 40°C. The generation of the training sequence is further regulated by the requirement that a random temperature is selected and the slope is calculated so that the selected temperature was reached in ten pulses. The selection of ten pulses was done by dividing the available range of temperatures (20–425°C) by 40. In other words, starting from 20°C, the temperature can go up to the maximum temperature of 425°C in ten steps. The particular selection of the training scheme allows the user

to explore a smaller, but more meaningful excitation region. The experimental data set used in modeling had one input (temperature) and one output (conductance). The test gas was hydrogen (80 ppm with balance nitrogen). The sampling rate was 0.4 s per measurement. The temperature profile consisted of 10000 temperatures, and the same temperature profile was repeated several times to confirm that the small features observed in the conductance response were not noise. The entire experiment required 2 h and 5 min. The data set was partitioned into two parts: training data consisting of 7000 samples, and model evaluation data consisting of 3000 samples. Fig. 4 shows a section of the data set for two successive cycles. Note that the y scale in Fig. 4 is on the order of 10 − 7 V − 1, which is extremely small. Thus, both general features and small details are very reproducible from cycle to cycle.

4.2. Modeling results To model the sensor response, the first step is to select the model orders, nu and ny using the training data set. Then, the problem of finding the appropriate functional mapping, F, reduces to searching for optimal parameters for a specified model structure. Well-known regression algorithms such as linear least squares, Gauss–Newton methods, etc. exist to compute the parameters for the selected model structure. In this paper, four different modeling approaches are considered as the functional mapping, F. These four model structures are: (1) linear autoregressive exogenous

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

29

Fig. 5. Comparison of different models on the test set (Part 1).

(ARX) [19], (2) neural network partial least squaresnonlinear autoregressive with exogenous inputs (NNPLS-NARX) [20], (3) continuous time models, such as Runge Kutta neural net (NNRUNGE) [21], and (4) wavelet networks (WNET) [17].

Initially, simple model structures are considered. The complexity of the model structure is only increased when the predictive model obtained has failed to capture the systems dynamics. The predictive capability of each model is compared by computing the normalized

Fig. 6. Comparison of different models on the test set (Part 2).

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

30

Fig. 7. Comparison of WNET models with different size of training data: 3000 samples versus 7000 samples (Sensor-1).

mean squared error (NMSE) for the multi-step ahead prediction over the test data set: MSE +





1 N % (y − yˆi )2/ % (y − y¯ )2/(n − 1) N i=1 i

(2)

where yˆ denotes the model prediction, and the term in the denominator is an estimate of the variance. Akaike’s information criterion (AIC) [22] was used to propose a model order for the linear ARX model. The

AIC gave the values ny = 10, nu = 20, and zero order dead time. The zero order dead time for the input was expected since the current conductance was highly correlated with the current temperature pulse. The coefficients for the linear ARX model were computed using the linear least squares method. The NMSE for this linear ARX model was 0.805. Multi-step ahead predictions using the linear ARX model are shown in Fig. 5. In the multi-step ahead

Fig. 8. The T 2 statistic comparing the WNET model with the actual data (3000 samples).

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

31

Fig. 9. The T 2 statistic comparing the WNET model with the actual data (7000 samples).

prediction mode, the network predicts the next output using previous predictions for output instead of actual output values in Eq. 1. In other words, the temperature profile alone is enough for predicting the conductance response of gas i (i.e. y i =Fi (u)). As can be seen in Fig. 5, the fit is not satisfactory for a linear model. The linear model cannot predict large changes in the actual conductance measurements. Since a linear ARX model was not able to capture variations in the actual conductance values, a nonlinear NNPLS-NARX model was applied. The NNPLS uses partial least squares (PLS) to remove the linear correlation among the input variables and employs a nonlinear neural network model in the score space [20]. The NNPLS model structure was selected as 3×6× 1 (three inputs with ny =1 and nu =2, six hidden neurons, and one output) using the validation data set method discussed in [20]. The results for long-term NNPLSNARX predictions are also shown in Fig. 5. While performing slightly better than the linear ARX model, it is not accurate enough to follow the actual conductance changes. The NMSE of this model was 0.537. The NNRUNGE method was previously used to accurately model bifurcations and long-term dynamics for a complex autonomous system involving surface reactions [21]. Since the actual conductance on a microhotplate chemical gas sensor depends greatly on the surface reactions, it was anticipated that this method might be useful in modeling a sensor’s transduction mechanism. The original method, however, had to be modified to handle the input – output data structure. In this approach, a neural network is identified to model

the first derivative of conductance using the actual temperature and conductance values. The neural network is integrated using the Runge–Kutta method instead of differentiating the experimental data numerically, which would amplify the noise in the data. For the NNRUNGE method, the model structure was selected as 2× 5× 1 by trial and error. This structure has two inputs since the neural network predicts the time derivative of the conductance at time t using the current conductance and temperature. For clarity, the result for this model is plotted separately in Fig. 6, with the same time frame enabling an easier comparison between previous model’s (see Fig. 5). The NNRUNGE model does a reasonable job in capturing slow dynamics in the data. It clearly reconstructs each peak up to a certain accuracy. Fast dynamics, however, cannot be represented by this model. Since NNRUNGE is composed of an inherent integrator, it smooths the actual response curves. Small features, which are not noise, are not reproduced in the predictions. A comparison of the models frequency response plot with the actual data also showed that this structure captures low frequencies at the expense of neglecting the high frequency behavior present in the experimental values. The NMSE was 0.180. The results with the NNRUNGE method show that a modeling approach with multi-resolution capability is required to accurately model the conductance responses measured with micro-hotplate chemical gas sensors. Over the past few years, wavelets have been extensively used to accurately represent both time and frequency information of signals, hence overcoming the limita-

32

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 10. (a) Normalized conductance response in methanol (solid line) and ethanol (dashed line); (b) the corresponding linear temperature pulse.

tions of Fourier analysis [23]. The recently proposed wavelet networks approach (see Appendix A for details) combines the multi-resolution feature of the wavelet transformation with neural networks, and it is capable of approximating any function to an arbitrary accuracy. In this example, Mexican hat wavelets (see Appendix A) are used with up to four scale levels. The initial model structure is synthesized using a stepwise selection by the orthogonalization procedure [24], which selected 88 wavelets out of 247 possible wavelets by the procedure outlined in Appendix A. This initial model is further trained for four epochs using the Gauss –Newton method. The whole training took about 4 h on a SPARC system 12 workstation. The WNET model produces an excellent fit which follows both slow and fast dynamics as shown in Fig. 6. The NMSE for WNET predictions was 0.147.

4.3. Model 6alidation When the data-based models are used for future predictions, some error prediction mechanism is needed to indicate how accurate the model is. Since the accuracy of nonlinear black box models depends greatly on the training data, it is dangerous to extrapolate these models outside of the training region. One way of 2 The equipment/software manufacturer is identified only to specify research procedure. Its mention in no way implies recommendation or endorsement of that product by the National Institute of Standards and Technology.

avoiding extrapolation is to check whether the current operating conditions are consistent with the data set from which the original model is synthesized. Multivariate statistical analysis tools such as principal component analysis (PCA) have been used for this purpose [25]. The correlation structure of a new set of data can be compared to the training set by using hypothesis tests such as the squared prediction error (SPE) and Hotelling’s T 2 test as explained in Appendix B. The SPE measures the change in the correlation structure, whereas the T 2 test detects the change in the operating conditions. Another difficulty arises when the NARX models are used for multi-step ahead predictions. The stability of the NARX models remains an open question. If these models are iterated using the previous predictions instead of the past values, the error between the true and predicted values might increase continuously, up to the point where the model becomes unstable. This phenomenon was observed by changing the number of samples in the training data set. A second WNET model was developed using 3000 training data samples instead of 7000. As seen in Fig. 7, the WNET model using 3000 data samples for training cannot follow the conductance excursions from 520 to 565 s. Further Gauss–Newton iterations for this WNET do not improve the predictions. This observed mismatch between the actual values and predictions occurs for no apparent reason, although increasing the number of samples in the training data set seems to remove this temporary mismatch.

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

33

Fig. 11. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum (SSD) temperature profile.

Since including more and more samples in the training set makes the model more robust, one explanation for this mismatch is that there are not enough examples in the training data set to prevent the model from extrapolating. So a metric is needed to warn the user when the model is not providing reliable predictions. Hotelling’s T 2 statistic was considered as such a test. This statistic measures a new sample’s distance from the multivariate mean in the reduced dimension space. The calculation of the T 2 test is presented in Appendix B. First, an input space was constructed using the training data (3000 samples). Since the model order was ny = 1, nu,=2, the data matrix, X, consisted of three columns and 2999 rows: u2 u1  Á y1 à y2 u3 u2 à à y3 u4 u3 à X=à — — — à à y i u i+1 u i à à — — — à Äy N − 1 u N u N − 1Å

using the experimental data (see Fig. 8). The T 2 statistic went out of the 95% confidence limits for the samples between 540 and 560 s. Similarly, the WNET model trained with 3000 data points failed to predict the sensors response for the same range as shown in Fig. 7. The same procedure was repeated using 7000 data points instead of 3000, and the T 2 values for the WNET model using 7000 training samples were compared with the experimental data. As can be seen in Fig. 9, not only did the WNET model stay below the 95% confidence limits, but the experimental data was also below the limit. This example shows that the T 2 statistic is a useful metric for detecting when accuracy of the predictive model is low. This statistic is used as a constraint in the off-line optimization of micro-hotplate sensors. Moreover, the same statistic can be used to decide size requirements for the training data set.

(3)

Then, the PCA algorithm (see Appendix B) was applied to the X matrix. Two principal components were kept in the PCA model, which was able to explain 87% of the total variance. The 95% confidence limit on T 2 was calculated to be 6.05 based on the first two principal components. Next, multi-step ahead predictions for the test data set were projected on the first two principal components and compared with the T 2 values obtained

5. Optimization of micro-hotplate gas sensors The aim of off-line optimization is to increase selectivity by choosing the ‘best’ temperature profile for microsensor operation. For a given temperature profile u, let the predicted conductance response for methanol and ethanol be y MeOH and y EtOH, respectively. Since a functional mapping, Fi, can be developed for each gas i, following the approach of the previous section, an off-line optimization problem can be formulated as [26] max

u d(y MeOH, y EtOH)

(4)

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

34

Fig. 12. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum (Haar: [2 3]) temperature profile.

where y MeOH = FMeOH (u) and y EtOH =FEtOH (u) and u is the temperature profile. A metric is needed to measure and to quantify the distance, d, between two response curves. One simple and straightforward metric is the sum of squared differences (SSD) between two response curves (conductance measurements in two different gases for the same cyclic temperature program). The normalized SSD measures the area between two response curves, and is defined as NSSD=

n y MeOH −y EtOH 22 = % (y MeOH −y EtOH )2/n i i n i=1

(5) where n is the number of temperature pulses in a cycle. Other parameters such as the pulse duration are chosen based on previous experience. Later, these assumptions are relaxed and the problem formulation is generalized to accommodate different parameters. The search space for the optimal temperature profile is over a limited subset of realizable actual temperature pulses, u= (u1 u2…un ). ul 5 ui 5uu

(6)

where ul and uu are the lower and upper limits, respectively. The lower limit, ul, is considered as 20°C, while the upper limit, uu, is chosen to be 425°C, and depends on the sensor structure. The lower and upper limits (ul and uu, respectively) for the magnitude of the temperature pulses are the most basic constraints used in the formulation of the optimization problem. Another set of constraints, however, arises

from the way that the random training sequences are generated. The continuity requirement, which was adopted when generating the random training pulses, dictates that the difference between two successive temperature pulses cannot be greater than 40°C in order to prevent sudden drastic changes in the surface condition and to interrogate the sensor response in a consistent manner. Without this constraint it was difficult to develop accurate predictive models for the sensors studied. This constraint is included in the optimization by the following equations: ui + 1 − ui 5 40°C

(7)

i= 1,…., n− 1 un − u1 5 40°C

(8)

Eq. (8) arises from the cyclic nature of the temperature input. Furthermore, a constraint is needed to prevent the models from extrapolating or producing results that are not consistent with the previous observations. Hotelling’s T 2 test is proposed as such a constraint. The optimization algorithm searches for the optimal temperature profile which maximizes the distance, d, without violating the 95% limit on the T 2 values. T 2(y MeOH)5 T 2limit 2

T (y

EtOH

)5 T

(9)

2 limit

(10) 2 limit

The same 95% limit, T , is used for both response curves, provided that both WNET models are trained with the same number of data points, and the same number of latent variables is kept when computing the PCA model for these gases.

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

35

Fig. 13. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum (Haar: [0]) temperature profile.

5.1. Optimization of temperature profile using NSSD d metric for a fixed number of pulses Before presenting the optimization results, the normalized conductance responses in methanol and ethanol for the case of a linear temperature profile are shown in Fig. 10. The reasons for normalization and its implementation are as follows. Micro-hotplate chemical gas sensors possess an advantage in dealing with drift problems. The conductance response, which is mea-

sured on a microsensor using a cyclic temperature program, can be decomposed into two basic parts: (1) its shape, and (2) the mean and the variance of the shape. The conductance response for each cycle of temperature profile is normalized to zero mean and unit variance. The mean and the variance are calculated for each cycle of the response rather than using the mean and the variance of a signal with multiple cycles, since these parameters are slowly changing during an experiment. If the drift is modeled as a linear change in the

Fig. 14. (a) Normalized conductance response in methanol (solid) and ethanol (dashed); (b) for temperature pulses: 275 and 380°C.

36

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 15. (a) Normalized conductance response in methanol (solid) and ethanol (dashed); (b) for ten temperature pulses: five at 275 and five at 380°C.*

sensors response, it can be shown that the normalization for each cycle preserves the original shape. As a result, using temperature programmed sensing mode with a proper signal processing approach eliminates the adverse effects of sensor drift in detecting chemical compounds. The off-line optimization problem was solved using

the sequential quadratic programming (SQP) algorithm provided in the MATLAB1 optimization toolbox [27]. Fig. 11 shows the resulting response curves (both predictions and experimental results) after the optimization. Note that the predictions for both gases fit very well with the experimental results. The computed optimal profile consists of temperature oscillations around

Fig. 16. (a) Conductance response in ethanol; (b) conductance response in methanol; (c) the corresponding temperature pulses; temperature pulse durations are 30 (solid), 40 (dashed), 50 (dot dashed), and 100 ms (dotted).

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

37

Fig. 17. (a) Normalized conductance response in methanol (solid line) and ethanol (dashed line); (b) the linear temperature pulse (D =100 ms).

300°C with a small amplitude. The frequency of these oscillations depends on the particular distance metric used in the off-line optimization. In contrast to the linear results of Fig. 10, the NSSD for the actual experiments is calculated to be 3.89 using this optimal temperature profile. When the gas is interrogated with the optimal temperature program, the methanol response follows the temperature profile, while ethanol gives an out of phase response. A different distance metric based on a modified response vector, 6, can be obtained with the Haar transform [28]. A signal, y, with n= 2m + 1 data points can be decomposed into m +1 scales, starting from scale zero up to scale m (0 1… m). Lower scales (the first few components of 6) contain information about the general shape of the actual signal y, while higher scales describe the details of the signal. The distance metric can be modified by keeping some of the scales and discarding others. Initially, n is set to 64, hence m is five. When the scales two and three were kept in the objective function, the numerical optimization algorithm converged to the result shown in Fig. 12 with NSSD =3.90. Again, the response curves in methanol and ethanol are completely out of phase with each other. The shape of the optimal temperature profile is of interest. Although both optimization algorithms started with the same initial temperature profile, the Haar transform with scales [2,3] converged to slower oscillations in the temperature profile. This phenomenon occurs because removing higher scales in optimization

forces the solution to a smoother form. When scales [2,3] were used, the distance is calculated for only 22 + 23 = 12 elements, instead of 64 as in the case of directly using the SSD metric. The same procedure can be repeated by using scale zero alone in the objective function. For scale zero, there is only one element (20 = 1). The resulting optimal temperature profile, and conductance responses are shown in Fig. 13, with NSSD equal to 3.85 for the actual measurements. Furthermore, the final optimal profile (Fig. 13b), can be approximated as a step change from a lower temperature ( 300°C) to a higher one ( 400°C). One should note that the optimized profile cannot be exactly a step change, since that solution is infeasible for the optimization problem with the continuity constraints (Eqs. 7 and 8). However, there may be a practical advantage to simpler profiles. Another set of experiments is performed to show the importance of the dynamics for the temperature profile optimization. Two temperatures (275 and 380°C) are selected based on (Fig. 13b). The microsensor is pulsed back and forth between these two temperatures. The measured conductance responses in methanol and ethanol are shown in Fig. 14. Both conductances go up and down in the same way, which is out of phase with the temperature excitation. In this figure, there is very little difference between the two conductance responses. However, if the same temperatures are repeated five times each, the measured conductance responses starts to show some differences, especially for the low temperature pulses (see Fig. 15). When the temperature pulse

38

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 18. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum temperature profile (D= 100 ms) for n =16.

is switched from 380 to 275°C, both conductances increase first, but then the conductance in methanol decreases, while the conductance in ethanol keeps increasing. This pair of experiments showed that more than two temperature pulses are required to observe the differences between conductance responses in similar reducing gases. The minimization of the cycle length for analyte response separation is considered next.

5.2. Minimum detection time The next objective is to minimize the cycle length, while still differentiating similar gases. In order to observe the effect of different pulse durations, D, on the sensor’s dynamics, characterization experiments are applied for various pulse durations with both methanol and ethanol. The top plot in Fig. 16 shows a portion of conductance responses in ethanol for different pulse durations (D =30, 40, 50, 100 ms), while the middle plot is for methanol. The corresponding temperature pulses are shown in the bottom plot of the same figure. When temperature pulses are switched from 450 to either 250 or 300°C, the conductance response in methanol (middle plot) increases monotonically, independent of the pulse duration. In ethanol, however, changing the pulse duration affects the observed response characteristics. For shorter pulse durations (e.g. D =30 ms), the conductance response in ethanol behaves similar to the response in methanol: it increases

with repeated temperature pulses. This behavior changes as the temperature duration is increased, so that the response in ethanol becomes completely different than that in methanol for a pulse duration of 100 ms. Hence, the duration is set to a minimum viable value of D= 100 ms, and WNET models are obtained for both methanol and ethanol gases. Before discussing the time-minimized results, the sensor’s response for a linear ramp temperature profile is presented in Fig. 17 for purposes of comparison. The linear ramp profile consists of 82 temperature pulses, i.e. u =(20, 25, 30,..., 420, 425). It takes about 10 s to complete one cycle of the temperature profile, since each pulse duration is 100 ms. The normalized conductance response curves (experimental data) for methanol and ethanol (shown in Fig. 17) are very similar to each other. The NSSD between these two curves is 0.07. The optimization of the temperature profile is carried out for varying number of temperature pulses in a cycle (n = 16, 20, 40, 58 and 64), and the results for profile optimization are shown in Figs. 18–22, respectively. The NSSD is used as the objective function. NSSD for the actual measurements and the cycle duration are reported for each case in Table 1, along with the linear ramp result. The largest NSSD is obtained for 58 temperature pulses, which takes 6.67 s to complete one cycle. The best result (both small cycle length and large NSSD) is obtained, however, for n= 20, since the cycle duration is decreased 3-fold while NSSD increases by two. It takes 2.3 s to complete the optimal temperature

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

39

Fig. 19. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum temperature profile (D = 100 ms) for n =20.*

profile presented in the lower half of the figure. The NSSD for the actual measurements is 1.28. Here, the phase shift between the two response curves was much more apparent compared to the linear ramp results. The optimization algorithm proposed a triangular shape for the pulse sequences starting at 155 going up

to 420°C with approximately 25°C increments. The above result was obtained with pulse durations of 100 ms. If one wanted to decrease the cycle length by decreasing the pulse duration, the difference between the two response curves would decrease. In order to illustrate this fact, the same optimal temperature profile

Fig. 20. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum temperature profile (D= 100 ms) for n =40.

40

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

Fig. 21. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum temperature profile (D= 100 ms) for n =58.*

was performed with 30 ms pulse durations. The results for these experiments are shown in the top plot of Fig. 23. The response in methanol did not change much when the pulse duration was decreased, while the response in ethanol deteriorated, and became similar to the one in methanol. The NSSD for this case (D= 30

ms) is 0.30, almost five times lower than the optimal case (D= 100 ms). The reason for this deterioration is that the response in ethanol is highly affected by the pulse duration, as already illustrated in Fig. 16. If the pulse duration is not long enough, then the dynamic response is not different from that for methanol. Vari-

Fig. 22. (a) Actual experiments in methanol (solid line) and ethanol (dashed line) gases, model predictions are shown by circles for methanol and plus for ethanol models; (b) the optimum temperature profile (D= 100 ms) for n =64.

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43 Table 1 A summary of the profile optimization results for different number of pulses (n) in a cycle n

Cycle duration

NSSD

16 20 40 58 64 82a

1.84 2.30 4.60 6.67 7.36 9.43

0.66 1.28 0.36 2.66 1.29 0.07

a

The linear ramp result.

ous factors might cause the observed dynamic differences: diffusion at the surface and in the near subsurface region, variations in the heats of reactions, different reaction paths, etc. Further experiments have to be carried out in order to decide which of these factors are predominantly responsible for producing different dynamics for similar compounds.

6. Conclusions This paper has discussed optimal tuning of the operational temperature profile for micro-hotplate chemical gas sensors. The approach involves off-line

41

optimization and dynamic data-based modeling techniques. Temperature modulation was employed in this study to increase response selectivity between methanol and ethanol gases injected into air. Databased dynamic models have been developed for predicting conductance responses of each gas for a given temperature profile. Among the different dynamic modeling techniques studied, the wavelet networks method gave very accurate predictive models. Using these predictive models in an off-line optimization scheme, an optimal temperature profile was computed and validated through actual experiments. This optimal temperature program produced methanol and ethanol responses that were out of phase as temperature was controlled, as a function of time. Using optimization the cycle length for analyte separation was reduced from 10 to 2 s, while increasing the distance metric 20-fold for the identification of these two gases. Although methanol and ethanol gases have been chosen to illustrate the optimal tuning of microhotplate gas sensors, the methodology described in this paper should be extensible to a range of analytes for various applications. This study has demonstrated that the specific operation of a micro-hotplate sensor can be optimized for a particular application using a structured design approach in the temperature programmed sensing mode.

Fig. 23. Optimum temperature profile performed with 30 ms pulse durations.

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

42

2 error, N i = 1 (fn (xi )− yi ) , procedure.

Appendix A. Wavelet networks In this paper, wavelet networks (WNET) [17,29] are used for the dynamic modeling of a micro-hotplate gas sensor’s response. The approach allows one to exploit the similarity between the discrete inverse wavelet transform and neural networks with one hidden layer.Given a wavelet function C:RdR, its associated wavelet network is written as p

ti )) + g¯ fp (x) = % wi C(a*(x− i

(11)

i=1

where d is the input dimension, p is the number of wavelet bases, ai is the dilation parameter, ti is the translation parameter, wi  R, ai  Rd, ti  Rd, and * denotes a component-wise product of two vectors. The parameter g¯ is introduced to help in dealing with nonzero mean functions, since the wavelet C has zero mean. Here, C belongs to the set of square-integrable functions, i.e. C (x) L 2 (R), and satisfies an admissibility condition, which implies that the mean value of C (x) be zero [30],

&

+ −

C(x) dx = 0

(12)

The wavelet function C is chosen as the so-called ‘Mexican hat’, C(x) = ( x 2 − d) exp



− x 2 2



(13)

where x 2 = x Tx. All the parameters ai, ti, wi, g¯ and p of the network in Eq. 11 are to be adapted using training data. An outline of the training (learning) procedure is given below whereas the details are in [31]: Given: A set of training data (xi, yi ) where xi is an input (1× d), and yi is an output (1× 1) with i =1,…, N, then (1) set g¯ to be the arithmetic mean of the training data. (2) Construct a regressor library W of discretely dilated and translated versions of the wavelet C. The training data is scanned and a subset of some regular (dyadic) lattice of dilated (ai) and translated (ti ) versions of C, whose supports contain the training data points, is selected. (3) Find wi by the least squares method. (4) Select the ‘best’ p regressors from W which minimize Akaike’s final prediction error criterion (FPE) using stepwise selection by orthogonalization. FPE is defined as FPE(fp )=

1+ np /N 1 N % (f (x ) − yi )2 1=np /N 2N i = 1 p i

(14)

and np = p(d+ 2)+ d+1

where np is the total number of parameters in a wavelet network. (5) Adjust the parameters (ai, ti, wi ) of the resulting network to further decrease the prediction

a

quasi–Newton

Appendix B. Calculation of T 2 values Principal component analysis [32] is used to reduce the dimensionality of variables, that is to find an orthogonal basis which will explain most of the variance in the observation space. The (n× m) data matrix X with n observations and m variables, is decomposed using singular value decomposition (SVD) X= USP T

(16)

where U is a (n× n) orthogonal matrix and has the eigenvectors of XX T as its columns. Similarly, the (m× m) matrix P is composed of the eigenvectors of X TX, and finally S is a diagonal matrix consisting of the square roots of the ordered eigenvalues of X TX. In PCA, US is defined as the score matrix T, and the columns of P are called the loadings. Thus, X= TP T

(17)

X can be written as the bilinear product X= t1p T1 + t2p T2 + …+ tmp Tm

(18)

where subscripts denote the column number, and the score matrix T can be obtained from T= XP

(19)

By choosing a smaller number of terms r, where rm, the original data matrix X can be approximated as X. : X. = t1p T1 + t2p T2 + …+ trp Tr

(20)

The PCA model is thus obtained by keeping the first r terms in the projection. In the PCA model, the original data matrix X is accurately represented using a smaller number of fictitious variables (principal components) which are linear combinations of the original variables. The difference between the original data set and its approximation is expressed as the squared prediction error (SPE) SPE= (X− X. )2

(21)

The SPE, also called the Q statistic, is a widely used statistic associated with PCA models. It describes the lack of fit between the original data set and its projection to a smaller number of principal components. Hotelling’s T 2 statistic, however, measures the variation in each sample within the PCA model. T 2 is defined as the sum of the normalized squared scores [33] T 2 = TL − 1T T

(15)

by

(22)

where T is now the matrix of r score vectors from the PCA model and L − 1 is the diagonal matrix containing the inverse of the eigenvalues associated with the r eigenvectors (principal components) kept in the model.

T.A. Kunt et al. / Sensors and Actuators B 53 (1998) 24–43

The F distribution is used to calculate the statistical confidence limits for the T 2 values: T 2r, n =

r(n− 1) F n− r r, n − r

(23)

where n is the number of observations (number of rows in the data matrix X) and r is the number of principal components retained in the PCA model. A geometrical interpretation of the Q and T 2 statistics is as follows. Q is the measure of the distance off the hyperplane described by the PCA model, and in fact, the Euclidean distance between a data point and its projection to the principal component axis is equal to the square root of Q. Hence, the confidence limit on Q defines a region around the hyperplane based on the PCA model. If the covariance structure captured in the PCA model is not enough to explain a sample, then the SPE will exceed the limit. Hotelling’s T 2 statistic, however, measures the distance from the multivariate mean (which is the intersection of the principal components) to the sample’s projection onto the principal component basis.

References [1] Y. Hiranaka, T. Abe, H. Murata, Gas-dependent response in the temperature transient of SnO2 gas sensors, Sensors and Actuators B 9 (3) (1992) 177. [2] P.K. Clifford, D.T. Tuma, Characteristics of semiconductor gas sensors II. transient response to temperature change, Sensors and Actuators B, 3: 233–254, 1982–1983. [3] P.K. Clifford, Mechanisms of gas detection by metal oxide surfaces, PhD thesis, Carnegie-Mellon University, 1981. [4] S.R. Morrison, Chemical Physics of Surfaces, Plenum Press, New York, 1990. [5] S. Wlodek, K. Colbow, F. Consadori, Kinetic model of thermally cycled tin oxide gas sensor, Sensors and Actuators B 3 (1991) 123 – 127. [6] S. Nakata, S. Akakabe, et al., Gas sensing based on a nonlinear response: discrimination between hydrocarbons and quantification of individual components in a gas mixture, Anal. Chem. 68 (1996) 2067 – 2072. [7] A. Oppenheim, R. Schafer, Discrete-Time Signal Processing, Prentice-Hall, New York, 1989. [8] A. Heilig, N. Baˆrsan, et al., Gas identification by modulating temperatures of SnO2-based thich film sensors, Sensors and Actuators B 43 (1 – 3) (1997) 45–51. [9] R.E. Cavicchi, J.S. Suehle, et al., Fast temperature programmed sensing for micro-hotplate gas sensors, IEEE Electron Device Lett. 16 (6) (1995) 1–3. [10] J.S. Suchle, R.E. Cavicchi, M. Gaitan, S. Semancik, Tin oxide gas sensor fabricated using CMOS micro-hotplates and in-situ processing, IEEE Electron Device Lett. 14 (3) (1993) 118 – 120.

.

43

[11] R.E. Cavicchi, J.S. Suehle, et al, Micro-hotplate temperature control for sensor fabrication, study and operation in: Proceedings of the Fifth International Meeting on Chemical Sensors, 1994, pp. 1136 – 1139. [12] T.A. Kunt, Dynamic modeling and optimization of micro-hotplate gas sensors, PhD thesis, University of Maryland, College Park, MD, 1997. [13] S. Semancik, R.E. Cavicchi, M. Gaitan, J.S. Suehle, Temperature-controlled, micromachined arrays for chemical sensor fabrication and operation. [14] J.W. Gardner, P.N. Bartlett, A brief history of electronic noses, Sensors and Actuators B 18 – 19 (1994) 211 – 220. [15] C. Rhodes, M. Morari, Determining the model order of nonlinear input/output systems directly from data, in: Proceedings of the American Control Conference, Seattle, WA, 1995, pp. 2190– 2194. [16] J.D. Bomberger, D.E. Seborg, Estimation of model order from input – output data applied to radial basis function network identification in: Proceedings of the IFAC Symposium Adchem ‘97, Banff, Canada, 1997, pp. 31 – 36. [17] Q. Zhang, A. Benveniste, Wavelet networks, IEEE Trans. Neural Netw. 3 (1992) 88 – 898. [18] T. Kunt, T.J. McAvoy, R.E. Cavicchi, S. Semancik, Dynamic modeling and optimization of micro-hotplate chemical gas sensors, in: Proceedings of the IFAC Symposium Adchem ‘97, Banff, Canada, 1997, pp. 91 – 95. [19] L. Ljung, System Identification: theory for the User, PrenticeHall, Engelwood Cliffs, NJ, 1987. [20] S.J. Qin, T.J. McAvoy, Nonlinear FIR modeling via a neural net PLS approach, Comput. Chem. Eng. 20 (2) (1996) 147–159. [21] R. Rico-Martinez, K. Krischer, Discrete versus continuous time nonlinear signal processing of Cu electrodissolution data, Chem. Eng. Commun. 118 (1992) 25 – 48. [22] L. Ljung, System Identification: theory for the User, PrenticeHall, Engelwood Cliffs, NJ, 1987, ch. 16. [23] O. Rioul, M. Vetterli, Wavelets and signal processing, IEEE Signal Process. Mag., October 1991. [24] S. Chen, S. Billings, W. Luo, Orthogonal least squares methods and their application to non-linear system identification, Int. J. Control 50 (5) (1989) 1873 – 1896. [25] A. Negiz, A. Cinar, Monitoring of multivariable dynamic processes and sensor auditing, in: Proceedings of the IFAC Symposium Adchem ‘97, Banff, Canada, 1997, pp. 55 – 60. [26] T. Kunt, L. Ratton, et al., Towards the development of an artificial nose for chemical process applications, Comput. Chem. Eng. (Suppl.) 20 (1996) S1437 – S1442. [27] MathWorks, Natrick, MA, Optimization Toolbox User’s Guide, December 1992. [28] T.J. McAvoy, et al., A comparative study of signal processing techniques for clustering microsensor data, Sensors and Actuators B 41 (1 – 3) (1997) 105 – 120. [29] Q. Zhang, Wavelet Network Version 2.1. anonymous FTP: ftp.irisa.fr: /local/wavenet, May 1994. Public domain MATLAB toolbox. [30] I. Daubechies, The wavelet transform, time-frequency localization and signal analysis, IEEE Trans. Inf. Theory 36 (1990). [31] Q. Zhang, Using wavelet network in nonparametric estimation, Technical Report 833, IRISA, 1994. [32] T.W. Anderson, An Introduction to Multivariate Statistical Analysis, Wiley, New York, 1984. [33] Eigenvector Research, Manson, WA, BLS Toolbox for Use with MATLAB, version 1.5 edition, 1995.