Modeling and Optimization of a Robust Gas Sensor

! ! ! ! CIplus Band 3/2016 Modeling and Optimization of a Robust Gas Sensor Margarita A. Rebolledo C., Sebastian Krey, Thomas Bartz-Beielstein, Oliver Flasch, Andreas Fischbach, Jörg Stork  ! ! ! ! ! ! ! ! MODELING AND OPTIMIZATION OF A ROBUST GAS SENSOR Margarita A. Rebolledo C., Sebastian Krey, Thomas Bartz-Beielstein, Oliver Flasch, Andreas Fischbach, J¨org Stork SPOTSeven Lab, TH K¨oln University of Applied Sciences [email protected] Abstract In this paper we present a comparison of different data driven modeling methods. The first instance of a data driven linear Bayesian model is compared with several linear regression models, a Kriging model and a genetic programming model. The models are build on industrial data for the development of a robust gas sensor. The data contain limited amount of samples and a high variance. The mean square error of the models implemented in a test dataset is used as the comparison strategy. The results indicate that standard linear regression approaches as well as Kriging and GP show good results, whereas the Bayesian approach, despite the fact that it requires additional resources, does not lead to improved results. Keywords: Bayesian modeling, BMA, Design of experiments, Genetic programming, Linear regression, Lasso, Kriging. 1. Introduction Theoretically, there are many advantages for the implementation of Bayesian analysis [5]. The use of Bayesian models might represent a good alternative for industrial applications as they produce more informative results. The generation of a data-driven model to optimize the development of a carbon-monoxide sensor provides an opportunity to test these assertions on limited and sparse data. As a first approach, Bayesian robust linear regression is implemented and compared to standard regression methods and a genetic programming approach. Our goal is to learn the difference in performance from the tested methods when applied to this kind of data and to set future considerations for working with Bayesian models in a more demanding fashion. In recent years the need to reduce air pollution levels has gained more importance in the automotive industry. The efficiency increase of the 1 2 motor combustion process plays an important role for the reduction of pollution levels. This efficiency can be indirectly measured by monitoring the concentrations of carbon monoxide and other harmful gases released into the atmosphere. This paper focuses on the modeling and optimization of a carbon monoxide in-situ sensor. The sensor should be able to discern the carbon-monoxide concentration apart from the other exhaust gases. This is a difficult goal, because the sensor is exposed to and influenced by the other gases. Thus, the sensor output is not expected to be a direct result of the concentration of the gas of interest. Instead it will be the result of an underneath process influenced by all the other gases. At the end of the analysis we hope to obtain models from different methods with an improved sensitivity to carbon-monoxide concentrations. The models will be compared in order to check the performances differences and possible improvement opportunities. This paper is structured as follows: Section 2 describes the research configuration, i.e., data and experimental designs. Key features of the algorithms are introduced in Sec. 3. Section 4 presents results from the experiments. Finally, a discussion of the results in given in Sec. 5. 2. Problem 2.1 Data Description The data was collected following a response surface design of experiments (RS-DoE). This design constraints itself to the maximum and minimum expected concentration values of each gas under normal working conditions. Given the cost and time consumption required for the experiments, only a limited amount of samples could be measured. The minimum number of samples required to have a good system description and the real limit of possible realizable samples in the industrial testing station was balanced. Finally, a sample size of 80 was chosen. A summary of the data is shown in Table 1. This application example is anonymized due to confidentiality reasons. The data were standardized, meaning that every sample had its mean subtracted and was then divided by the standard deviation. The different gases are denominated as the variables X1 to X7. The values of interest correspond to the columns denominated Y1 and Y2, which are the sensor measurements. All the models will use this dataset as the training set. A general idea of the system behavior can be obtained by examining the correlation between the system output and inputs as shown in Table2. Some assumptions can be made about the influence each variable has on the sensor output: not all parameters seem to have the same influence on the sensor output. Also, the sensors do not behave Modeling and Optimization of a Robust Gas Sensor 3 Table 1: Overview of the standardized dataset used to generate the models of the sensors. Here every input of the model is denoted by an Xand every sensor output is denoted by an Y. X1 X2 X3 X4 X5 X6 X7 Y1 Y2 Minimum −1.13 −1.21 −1.16 −1.13 −1.15 −1.17 −1.00 −1.94 −2.06 1st Quartal −1.13 −1.21 −1.16 −1.13 −1.15 −1.17 −0.82 −0.63 −0.58 Median 0.09 0.03 0.12 0.08 0.08 0.05 −0.39 0.06 0.09 Mean 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3rd Quartal 1.30 1.26 1.40 1.29 1.28 1.28 0.59 0.66 0.67 Maximum 1.30 1.26 1.40 1.29 1.31 1.28 3.79 2.32 2.28 Table 2: Correlation between the system output and inputs for the training dataset X1 X2 X3 X4 X5 X6 X7 Y1 0.34 −0.19 −0.27 0.73 0.01 −0.00 −0.21 Y2 0.31 −0.16 −0.18 0.78 0.00 −0.03 −0.22 identically. Figure 1 shows the effect the two most strongly correlated parameters, X1 and X4, respectively, have on the sensors signal. A second dataset, denominated test set, which follow the characteristic of the previously described training set was made available to validate the results of the obtained models. 2.2 Experimental Design Considerations The data described in Sec. 2.1 was retrieved during tests based on an experimental design suitable for fitting models using the response surface methodology (RSM) [13]. First experimental results taken from screening design experiments indicated that a first order polynomial model is not sufficient, due to cross-sensitivity of the sensors. Therefore it was decided to use a RSM with two-factor interactions quadratic effects. In this case central composite designs would have been a logical choice. They are a combination of a box design, typically a full factorial or fractional factorial design and additional star or center points. Afull factorial design (FFD) with three levels for each of the six factors to estimate main effects and all quadratic terms would lead to 37= 2187 experiment runs at minimum. Choosing to divide factors 4 −1.0 −0.5 0.0 0.5 1.0 −2 −1 0 1 2 Concentration X1 Measured response Y1 (a) −1.0 −0.5 0.0 0.5 1.0 −2 −1 0 1 2 Concentration X4 Measured response Y1 (b) Figure 1: Scatter plots showing the general behavior of Y1 with respect to: a) the influence of X1 and b) the influence of X4 in two levels, a two level FFD would still need at least 27= 128 runs, without any repetitions or center points. A valid system description would need even more runs. To reduce the number of runs and be capable of fitting second order polynomial models a Box-Behnken design comes into consideration. But as a linear constraint on the input variables is limiting the sum of their values, almost all standard designs methods does not meet the requirements. Therefore, a more flexible design is needed. Using the statistical software JMP the RS-DoE was generated following the given constraints and applying the I-optimality criterion [12, 11]. The design was run and optimized a total of 80 times with 80 data points in order to obtain the best possible inference accuracy. 3. Algorithms 3.1 OLS A linear model estimated by ordinary least squares is the natural first modeling attempt for data generated by an experimental design. Our baseline model for the comparison of the different modeling methods is Modeling and Optimization of a Robust Gas Sensor 5 the linear main effects model f1: ˆy=β0+ 7 X i=1 βixi. In this work we used a RSM design, so beside the main effects the parameters of all two-way interactions and quadratic terms of the input variables can be estimated. This results in the full linear model f7: ˆy=β0+ 7 X i=1 βixi 7 X i=1 7 X j=i βijxixj. Based on the full linear model f7we applied variable selection based on an analysis of variance to get a more sparse model, which can be better interpreted. With a F-Test p-value of α= 0.01 as the decision boundary for the inclusion into the final model, we obtained the model f2: ˆy=β0+ 4 X i=1 βixi+β14x1x4+β34x3x4. The mean squared error (MSE) as defined in Eq. 4 on page 11, was used for our comparisons. While the full linear model f7has a lower training error, i.e., MSE of 0.11 for sensor Y1 and 0.10 for sensor Y2, compared to the baseline model (MSE of 0.24 for sensor Y1 and 0.23 for sensor Y2), the prediction performance on the test dataset is very weak (MSE of 7.76 and 9.08). This is a strong indicator of overfitting. The model f2has a MSE of 0.79 for sensor Y1 and 0.80 for sensor Y2, which is comparable (for sensor Y1) and little higher (for sensor Y2) than the baseline model. The residual standard error for the training set is lower (0.43 compared to 0.52 for sensor Y1 and 0.41 compared to 0.51 for sensor Y2) resulting in narrower confidence intervals for the parameters. The adjusted coefficient of determination (adjusted R2) is 0.81 compared to 0.73 for sensor Y1 and 0.83 compared to 0.74 for sensor Y2. This means the inclusion of the two two-way interactions X1:X4 and X3:X4 has a large contribution for the explanation of the variance in the dataset, while the input variables X5, X6 and X7 have very little or no contribution and can be left out of the model. 3.2 Lasso The Least Absolute Shrinkage and Selection Operator (Lasso) implements a selection method for linear models [9]. It selects solutions 6 with fewer parameter values, effectively reducing the number of variables upon which the given solution is dependent. The Lasso trains a linear model with a L1prior as regularizer. Give a set of input measurements X={xi}n i=1 and an outcome measurement y, the lasso fits a linear model ˆy=β0+ p X i=1 βixi. Let α≥0 be a constant. The Lasso uses the following optimization criterion: min β 1 2n||Xβ −y||2 2under the constraint ||β||1≤α, (1) where ||·||1and ||·||2denote the L1and L2-norm, respectively. The positive constant αis a tuning parameter. For large αvalues, the constraint ||β||1≤αin Eq. 1 has no effect and the usual linear least squares regression is performed. For smaller values of α, the solutions are shrunken versions of the least squares estimates. Decreasing the values of αforces the coefficients βi’s to become zero, i.e., choosing αresults in selecting the number of predictors to use in a regression model. The Lasso can recover the exact set of non-zero weights (under certain conditions). Coordinate descent is used to fit the coefficients. 3.3 Kriging Kriging or Gaussian process regression is a method of interpolation [10]. The nobservations in an arbitrary data set, Y={yi}n i=1 can be associated as a single point sampled from some multivariate (n-variate) Gaussian distribution. The observations and the Gaussian process are related to each other by the covariance or kernel function k(xi, xj). Kernel functions compute the distance between two samples in an arbitrary metric and apply a radial function to this distance. The squared exponential kernel, also known as the Gaussian radial basis function (RBF) kernel, is used in our study. This kernel is given by k1(xi, xj) = σ2exp(−θkxi−xjk2 2) with θ=1 2l2.(2) The RBF kernel can be interpreted as a similarity measure, because values of this kernel decrease with distance. They range between zero (in the limit) and one. The length parameter lin Eq. 2 determines the effect of other observations during interpolation at new xvalues. The RBF kernel was selected in our study, because Gaussian processes with Modeling and Optimization of a Robust Gas Sensor 7 this kernel generate smooth functions. Since noisy data were analyzed in our study, the white noise kernel k2(xi, xj) = σ2δ(xi, xj), where δ(xi, xj) denotes the Kronecker delta function, was added to the RBF kernel. Hence, we used the kernel function k=k1+k2. 3.4 Robust Bayesian Modeling From previous knowledge of the combustion process it is expected that not only the main predictors but also the interactions between predictors have an effect on the sensor reading. As a general rule, if the data contains Kvariables then the expected number of possible models will be 2K. The total number of variables in the dataset with all the interactions included accounted to 22, that is 4.19 ×103possible model combinations to describe the sensor reading. To reduce the dimensionality of the problem Bayesian model averaging (BMA) is implemented. This provides a way to account for the uncertainty in model selection and provide in average a better predictive ability [3]. BMA is implemented in the statistical programing language R using the Bayesian Model Sampling (BMS) package [8]. The results show that the predictors X1, X3, X4, X1:X4, X3:X4 and X2 seem to be the most important for a good model. Bayesian modeling is the mathematical relocation of credibility of parameters values for a model according to what can be inferred from the data. As the first taken approach the reduced model containing only 6 out of the 22 variables is defined using a linear relationship. The model was implemented using Just Another Gibbs Sampler (JAGS), which is a program for Bayesian modeling using Markov Chain Monte Carlo (MCMC) [6] and rjags [7] as a link between R and JAGS. The sensor responses Y1 and Y2 are modeled following a non standardized Student’s t-distribution. This distribution was selected assuming that the variance present in the sensor output, illustrated in Figure 1, served as an indicator of variance in the model response. The mean of the distribution is defined by the canonical linear formula of the linear regression. The spread of the data was set to have a wide range of probable values defined by an uniform distribution. The normality factor, expressed as an exponential distribution, have preference for values close to one. Given the limited prior information available for the experiment, weakly informative priors are assigned to the parameters. The coefficients priors, βi, are defined to have a normal distribution centered around zero and a large variance. The t-distribution normality factor υfavors values smaller than 30 and σallows for a wide enough 8 distribution. The prior distributions were chosen as follows (Eq. 3): α∼N(0,4), βi∼N(0,4), υ ∼Exp(30), σ ∼U(−1−4,10) (3) The MCMC simulation are executed on the defined model to sample the posterior distribution of the parameters of interest, α, βi, σ, and υ. The chains were specified to run 500 adaptive iterations, followed by 1,500 burn-in iterations. Afterwards, 15,000 samples were taken from the posterior distribution with a thinning factor of 20 steps. Investigating their trace plots and diagnostic statistics of the resulting MCMC object reveals that the chains have converged. A value for the GelmanRubin diagnostic statistic [2] of under 1.1 suggest a good convergence. The effective sampling size (ESS) backups this assumption. The visual and numeric diagnostics allows us to think that the resulting MCMC sampling is representative and accurate of the posterior distribution of the different parameters. The posterior distributions obtained from the MCMC sampling for each parameter coefficient can be seen in Table 3 together with the high density intervals (HDI) of 95%. The MSE for the fitted models is 0.16 and 0.15 for the sensor Y1 and Y2, respectively. Table 3: Posterior mean for the coefficients βifor i= 1, ..7 and the parameters σand υfor the models of Y1 and Y2. The lower HDI (LHDI) and upper HDI (U-HDI) limits are indicated for each entry. Y1 B0 B1 B2 B3 B4 B5 B6 σ υ Mean −0.01 0.32 −0.14 −0.28 0.72 −0.23 −0.14 0.41 37.67 L-HDI −0.09 0.23 −0.24 −0.38 0.63 −0.33 −0.42 0.34 5.13 U-HDI 0.09 0.42 −0.04 −0.18 0.82 −0.14 −0.05 0.49 117.66 Y2 B0 B1 B2 B3 B4 B5 B6 σ υ Mean −0.00 0.29 −0.10 −0.20 0.78 −0.26 −0.13 0.39 34.20 L-HDI −0.09 0.20 −0.19 −0.29 0.69 −0.35 −0.22 0.31 4.44 U-HDI 0.09 0.39 −0.10 −0.11 0.88 −0.17 −0.04 0.47 113.50 3.5 Genetic Programming Genetic programming is an evolutionary algorithm that searches the set of symbolic expressions defined by a set of basis expressions (building blocks) for expressions that minimize one or multiple loss (fitness) functions. Symbolic regression is the application of genetic programming to