- Email: [email protected]
Identification of Static Boundary Model Based on Gaussian Process Classification
Preprints, 8th IFAC International Symposium on Preprints, IFAC Preprints, 8th 8th IFAC International International Symposium on on Advances in Automotive Control Symposium Preprints, 8th IFAC International Symposium on Advances in Automotive Control Available online at www.sciencedirect.com Advances in Automotive Control June 19-23, 2016. Norrköping, Sweden Advances in2016. Automotive Control June 19-23, Norrköping, Sweden June June 19-23, 19-23, 2016. 2016. Norrköping, Norrköping, Sweden Sweden
ScienceDirect
IFAC-PapersOnLine 49-11 (2016) 787–792
Identification of Static Boundary Model Identification Identification of of Static Static Boundary Boundary Model Model Based on Gaussian Process Classification Based on Gaussian Process Classification Based on Gaussian Process Classification ∗∗∗ H. Oyama ∗∗ M. Yamakita ∗∗ K. Sata ∗∗ ∗∗ A. Ohata ∗∗∗ ∗ M. Yamakita ∗ ∗ K. Sata ∗∗ ∗∗ A. Ohata ∗∗∗ ∗∗∗ H. Oyama ∗ H. Oyama M. Yamakita K. Sata A. Ohata H. Oyama M. Yamakita K. Sata A. Ohata ∗ ∗ Dept. of Mechanical and Control Engineering, Tokyo Institute of ∗ of Mechanical and Control Engineering, Tokyo Institute of ∗ Dept. Dept. and Control Engineering, Technology, Tokyo, Japan Dept. of of Mechanical Mechanical and(e-mail: [email protected]) Engineering, Tokyo Tokyo Institute Institute of of Technology, Tokyo, Japan (e-mail: [email protected]) ∗∗ Technology, Tokyo, Japan (e-mail: [email protected]) Unit Management System Development Div., TOYOTA Technology, Japan (e-mail: [email protected]) ∗∗ Advanced Tokyo, ∗∗ Advanced Unit Management System Development Div., TOYOTA ∗∗ Advanced MOTOR Unit Management Management System Development Development Div., TOYOTA TOYOTA CORPORATION, Shizuoka, Japan Advanced Unit System Div., MOTOR CORPORATION, Shizuoka, Japan MOTOR CORPORATION, Shizuoka, Japan (e-mail: kota [email protected]) MOTOR CORPORATION, Shizuoka, Japan (e-mail: kota [email protected]) ∗∗∗ (e-mail: [email protected]) INC., kota Higashifuji Technical Center, Toyota Motor (e-mail: kota [email protected]) ∗∗∗ TECHNOVA ∗∗∗ TECHNOVA INC., Higashifuji Technical Center, Toyota Motor ∗∗∗ Corporation, TECHNOVA INC., Higashifuji Technical Center, Shizuoka, Japan (e-mail: [email protected]) TECHNOVA INC., Higashifuji Technical Center, Toyota Toyota Motor Motor Corporation, Shizuoka, Japan (e-mail: [email protected]) Corporation, Shizuoka, Shizuoka, Japan Japan (e-mail: (e-mail: [email protected]) [email protected]) Corporation,
Abstract: Due to the needs for CO2 emission reduction, the control region for automotive Abstract: Due to for CO reduction, the control region for 2 emission Abstract: Due to the the needs needs for CO emission reduction, theand control region for automotive automotive 2 engine system is pushed towardfor theCO boundary between normal abnormal engine operations Abstract: Due to the needs reduction, the control region for automotive 2 emission engine system is pushed toward the boundary between normal and abnormal engine operations engine system is pushed toward the boundary between normal and abnormal engine operations such assystem knocking and misfiring. When an engine is operated outside of admissible region, the engine is pushed toward the boundary between normal and abnormal engine operations such as knocking and misfiring. When an engine is operated outside of admissible region, the such as knocking and misfiring. When an engine is operated outside of admissible region, the engine or the catalyst is seriously damaged or passengers feel very uncomfortable. Therefore such as knocking and misfiring. When an engine is operated outside of admissible region, the engine or the catalyst is seriously damaged or passengers feel very uncomfortable. Therefore engine or the the catalyst catalyst is boundary seriously damaged damaged or passengers passengers feel very uncomfortable. Therefore the identification of theis is highly demanded. In feel this very paperuncomfortable. we propose a Therefore design of engine or seriously or the identification of the boundary is highly demanded. In this paper we propose aa design of the identification the highly demanded. this we propose experiment (DoE) of strategy based onis Gaussian processIn classifier (GPC) the expectation the identification of the boundary boundary isthe highly demanded. In this paper paper we with propose a design design of of experiment (DoE) strategy based on the Gaussian process classifier (GPC) with the expectation experiment (DoE) strategy based on the Gaussian process classifier (GPC) with the expectation propagation(DoE) (EP) strategy algorithm. In the experimental result, classifier the proposed algorithm can identify experiment based on the Gaussian process (GPC) with the expectation propagation (EP) algorithm. In the experimental result, the proposed algorithm can identify propagation (EP) algorithm. In result, can the boundaries of an engine benchmark problem, provided byproposed the joint algorithm research committee of propagation (EP) algorithm. In the the experimental experimental result, the the proposed algorithm can identify identify the boundaries of an engine benchmark problem, provided by the joint research committee of the Society boundaries of an an engine engine benchmark problem, provided by the joint joint and research committee of of Automotive Engineers (JSAE) and Society of by Instrument Control Engineers the boundaries of benchmark problem, provided the research committee of the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers (SICE). the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers (SICE). (SICE). (SICE). © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Gaussian processes; Classification; Identification algorithms; Engine modelling; Keywords: Gaussian Gaussian processes; processes; Classification; Classification; Identification Identification algorithms; algorithms; Engine Engine modelling; Keywords: Engine control; Model based control; Statistical design. Keywords: Gaussian processes; Classification; Identification algorithms; Engine modelling; modelling; Engine control; control; Model based control; control; Statistical design. Engine Model based Statistical design. Engine control; Model based control; Statistical design. 1. INTRODUCTION region are the convex hull based method (e.g. Rapid Hull 1. INTRODUCTION INTRODUCTION region hull method (e.g. Rapid Hull 1. region are are the the convex convex hull based based method (e.g. Rapidand Hulla Determination Algorithm, Renniger et al.(e.g. (2005)) 1. INTRODUCTION region are the convex hull based method Rapid Hull Determination Algorithm, Renniger et al. (2005)) and a Determination Algorithm, Renniger et al. (2005)) and a way that minimize the variance of the parameters (e.g. Determination Algorithm, Renniger et al. (2005)) and a To solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. DDTo solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. Doptimal design, Pukelsheim (2006)). However, the region To solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. Dproblem, automakers have developed high functional enTo solve the energy issue the global design, Pukelsheim (2006)). However, the region problem, automakers have and developed high environmental functional enen- optimal optimal design, the of normal enginePukelsheim operations (2006)). seems toHowever, be non-convex and problem, automakers have developed high functional optimal design, Pukelsheim (2006)). However, the region region gines. Although the engine control system has been comproblem, automakers have developed high functional enof normal engine operations seems to be non-convex and gines. Although Although the the engine engine control control system system has has been been comcom- of of normal engine operations seems to be non-convex and in the latter case, we usually deal with quite simple models normal engine operations seems to be non-convex and gines. plicated, the development becoming shorter to in the latter case, we usually deal with quite simple models gines. Although the engine period controlis system has been complicated, the development period is becoming shorter to in the latter case, we usually deal with quite simple models such as polynomial models. Therefore the basic identified plicated, the development period is becoming shorter to in the latter case, we usually deal with quite simple models meet market It is difficult to isdevelop products satisplicated, the needs. development period becoming shorter to such as polynomial models. Therefore the basic identified meet market market needs. It is is difficult difficult to develop develop products satissuch as polynomial Therefore the identified models never fit models. the actual boundaries. To address this meet needs. It to products satissuch as will polynomial models. Therefore the basic basic identified fying requirements in the short term through conventional meet market needs. It is difficult to develop products satismodels will never fit the actual boundaries. To address this fying requirements in the short term through conventional models will never fit the actual boundaries. To address this issue, a DoE strategy based on support vector machine will never fit the actual boundaries. To address this fying requirements in the short through conventional development environments. Forterm such a background the models fying requirements in the short term through conventional issue, a DoE strategy based on support vector machine development environments. For such a background the issue, a DoE strategy based on support vector machine (SVM) is proposed in Kampmann et al. (2012). This development environments. For such a background the issue, a DoE strategy based on support vector machine model based development (MBD), which isbackground a development development environments. For such a the (SVM) is proposed in Kampmann et al. (2012). This model based based development development (MBD), (MBD), which which is is aa development development (SVM) Kampmann et al. method is canproposed find andin not only but This also model is proposed indescribe Kampmann et convex, al. (2012). (2012). This approach using virtual models and simulators, attracts (SVM) model based development (MBD), which is a development method can find and describe not only convex, but also approach using virtual models and simulators, attracts method can find and describe not only convex, but also non-convex boundaries. On the other hand, this strategy method can find and describe not only convex, but also approach using virtual models and simulators, attracts much attention as one solution to the problem. approach using virtual models and simulators, attracts non-convex boundaries. On the other hand, this strategy much attention as one solution to the problem. non-convex boundaries. On the other hand, this strategy can determine only search directions for additional inputs. much attention as one solution to the problem. non-convex boundaries. On the other hand, this strategy much attention as one solution to the problem. determine only directions for additional inputs. To reduce CO2 emission, the control region for automotive can canfind determine only search search directions forslowly additional inputs. To boundaries, inputs directions are changed and usually can determine only search for additional inputs. To reduce reduce CO CO2 emission, the control region for automotive To find boundaries, inputs are changed slowly and usually To emission, the control region for automotive 2 engine system is pushed toward the boundary between To find boundaries, inputs are changed slowly and usually To reduce CO emission, the control region for automotive in small discrete steps. 2 To find boundaries, inputs are changed slowly and usually engine is toward the between small discrete steps. engine system system is pushed pushed toward the boundary boundary between in normal and abnormal engine operations such as knocking in small discrete steps. engine system is pushed toward the boundary between in small discrete steps. normal and abnormal engine operations such as knocking In this paper we a DoE strategy based on the normal and abnormal engine operations such as knocking and misfiring. The abnormal operation may damage the In this paper we propose normal and abnormal engine operations such as knocking propose aa DoE strategy based on the and misfiring. The abnormal operation may damage the In this paper we propose DoE on Gaussian process classifier with thebased expectation and misfiring. The abnormal operation may damage the In this paper we propose a(GPC) DoE strategy strategy based on the the engine seriously. Therefore the control method satisfying and misfiring. The abnormal operation may damage the Gaussian process classifier (GPC) with the expectation engine Therefore the method Gaussian process with the (EP) classifier algorithm(GPC) in Rasmussen etexpectation al. (2006) engine seriously. seriously. Therefore theis control control method satisfying satisfying Gaussian process classifier (GPC) with the expectation boundary constrains robustly highly demanded. To de- propagation engine seriously. Therefore the control method satisfying (EP) in et (2006) boundary constrains constrains robustly robustly is is highly highly demanded. demanded. To To dede- propagation propagation (EP) algorithm algorithm in Rasmussen Rasmussen et al. al.accurate (2006) and Minka (2001). This strategy can describe (EP) algorithm in Rasmussen et al. (2006) boundary sign such aconstrains controllerrobustly throughisMBD, course, accurate boundary highlyof demanded. To de- propagation and Minka (2001). This strategy can describe accurate sign such a controller through MBD, of course, accurate and Minka (2001). This strategy can describe accurate boundary models and provide an efficient input sign such a controller through MBD, of course, accurate and Minka (2001). This strategy can describe accurate numerical of these boundaries (Watan- boundary models and provide an efficient input design sign such amodels controller through MBD,are of required course, accurate design numerical models of these boundaries are required (Watanboundary models provide an design algorithm using theand generated That is,input this method numerical models of these boundaries are required (Watanboundary models and providemodel. an efficient efficient input design abe et al. (2014) and Kieft (2014)). To identify boundaries numerical models of these boundaries are required (Watanalgorithm using the generated model. That is, this method abe et al. (2014) and Kieft (2014)). To identify boundaries algorithm using the generated model. That is, this method can describe non-convex region and can determine the using the generated model. That is, this method abean etengine al. (2014) (2014) and Kieftthe (2014)). To method identify needs boundaries at testand bench, common a lot algorithm abe et al. Kieft (2014)). To identify boundaries describe non-convex region and can determine the at an an engine engine test test bench, bench, the the common common method method needs needs aa lot lot can can describe non-convex region and can determine the next measurement point near the boundaries. By this at can describe non-convex region and can determine the of time and cost. The main objective method of this paper to next measurement point near the boundaries. By this at an engine test bench, the common needs ais lot of time and cost. The main objective of this paper is to next measurement point near the boundaries. By this direct determination of additional inputs as points, the of time and cost. The main objective of this paper is to next measurement point near the boundaries. By this obtain models withofanthis efficient of timeaccurate and cost.boundary The main objective paper input is to direct determination of additional inputs as points, the obtain accurate boundary models with an efficient input direct determination of points, the calibration time is expected to beinputs much as direct determination of additional additional inputs asshorter. points, The the obtain accurate boundary models with an efficient input design algorithm. obtain accurate boundary models with an efficient input calibration time is expected to be much shorter. The design algorithm. calibration time is expected to be much shorter. The validity of the proposed method is demonstrated by using design algorithm. algorithm. calibration time is expected to be much shorter. The design of the method is by using In the field of automotive engine design, we usually con- validity validity of benchmark the proposed proposedproblem method(Watanabe is demonstrated demonstrated by using an engine et al. by (2014)), validity of the proposed method is demonstrated using In the field of automotive engine design, we usually conan engine benchmark problem (Watanabe et al. (2014)), In the field of automotive engine design, we usually consider the static case. This means that we obtain an experian engine benchmark problem (Watanabe et al. (2014)), In the field of automotive engine design, we usually conprovided by the joint research committee of the Society of an engine benchmark problem (Watanabe et al. (2014)), sider static case. This that we an experiby the joint research committee of the Society of sider the the static case. This means means that we obtain obtain anThe experimental data after an engine reaches a steady state. ap- provided provided by the joint research of the Society of sider the static case. This means that we obtain an experiAutomotive Engineers (JSAE)committee and Society of Instrument provided by the joint research committee of the Society of mental data after an engine reaches a steady state. The apAutomotive Engineers (JSAE) and Society of Instrument mental data after an engine reaches a steady state. The approach data to determine inputs reaches is calleda design of experiment Automotive Engineers (JSAE) and Society of Instrument mental after an engine steady state. The apand Control (SICE). The empirical combustion Automotive Engineers (JSAE) and Society of Instrument proach to determine inputs is called design of experiment proach to inputs is of (DoE). The most common to design the control and and Control Control Engineers Engineers (SICE). (SICE). The The empirical empirical combustion combustion proach to determine determine inputsapproaches is called called design design of experiment experiment and Control Engineers (SICE). The empirical combustion (DoE). The most common approaches to design the control (DoE). The most common approaches to design the control (DoE). The most common approaches to design the control Copyright 2016 IFAC 801 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2016, IFAC (International Federation of Automatic Control) Copyright 2016 IFAC 801 Copyright © 2016 IFAC 801 Peer review© of International Federation of Automatic Copyright ©under 2016 responsibility IFAC 801Control. 10.1016/j.ifacol.2016.08.115
IFAC AAC 2016 788 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
profile with engine operating conditions and the knock model are implemented in this benchmark. The organization of this paper is as follows. In Section 2, we define the boundary modeling. We present the GPC and EP approximation algorithm in Section 3 and Section 4. In Section 5, we derive an input design algorithm based on GPC. Then in Section 6, we demonstrate the effectiveness of the proposed method by numerical simulations. Conclusions are given in Section 7. 2. BOUNDARY MODELING First of all, we explain about some definitions of the boundary modeling. We consider a discrete time nonlinear system represented by xk+1 = f (xk , uk , wk ), (1) yk = g(xk , uk , vk ),
Fig. 1. Example of boundary modeling strategy. To avoid engine damage, the inputs are changed slowly and usually in small discrete steps.
where x ∈ RN x , u ∈ RN u , and y ∈ RN y are the state vector, the input vector, and the output vector respectively. The subscript k denotes k-step value, w ∈ RN w and v ∈ RN v denote parameters of system uncertainty. We define the boundary model h(z) : RN z → R as follows: � > 0, outside of admissible domain on boundary (2) h(zk ) = 0, < 0, inside of admissible domain,
In this section we introduce Gaussian process classification in Rasmussen et al. (2006) and Bishop (2006). GPC is a Bayesian kernel classifier derived from Gaussian process priors. We assume that we have a training data set D = {(zi , ai )|i = 1, . . . , n}, where zi ∈ RD , ai ∈ {−1, 1}, and n denote input data, output data, and the number of samples respectively. The training set can be expressed as D = (Z, a) by using Z = [z1 z2 · · · zn ] and a = [a1 a2 · · · an ]. Given this training data set, the classification problem is to output the correct class label for a new input data z∗ . In GPC, we predict the class label of z∗ through the class posterior probability p(a∗ |D, z∗ ).
where zk ∈ RN z is the decision vector and we assume that H := {zk |h(zk ) < 0} becomes a connected set. 2.1 Static boundary modeling The goal of boundary modeling is to represent the true boundary numerically as h(z) in (2). However, in cases where systems are very complicated, it is difficult to determine the composition of the decision vector. The way how to acquire efficient data sets for identification is also an important problem. Therefore, in practice, we usually identify the boundary model using steady state data points to simplify the problems. In this paper, we define the modeling with steady state data set as static boundary modeling. In this � static�Tcase we can derive the decision vector as z = xT , uT . Here we define the steady state of the system x∞ ∈ RN x and constant inputs u∞ ∈ RN u , which satisfy x∞ = f (x∞ , u∞ , 0). Our objective of static boundary modeling is to identify the model as follows: �� �� � > 0, outside of admissible domain x∞ = 0, on boundary hS (3) u∞ < 0, inside of admissible domain. When elements of the state x∞ are determined only by constant input u∞ , we can omit the elements from decision vector z. The common approaches to identify the static boundary model are the convex hull based method (e.g. Rapid Hull Determination Algorithm, Renniger et al. (2005)) and Doptimal design in Pukelsheim (2006). These strategy can determine only search directions for additional inputs. To avoid engine damage, inputs are changed slowly and usually in small discrete steps (see Fig. 1). This procedure needs a lot of time and cost. 802
3. GAUSSIAN PROCESS CLASSIFICATION
In GPC, we assume that the probability of a class label depends on the value of some latent real value function f (z) ∈ R. That is, for binary classification, p(a = +1|z, f (z), D) = p(a = +1|f (z)). The probability of belonging to a class label ai = +1 for an input sample zi is monotonically increasing function related to the value fi (z) ∈ R of f (z). This monotonic relationship is defined according to a squashing function such as logistic and probit functions. For example, fi = f (zi ) can be derived as follows: 1 , logistic function p(ai = +1|fi ) = 1 + exp(−fi ) (4) Φ(fi ), probit function, where Φ normal function Φ(z) = √ is the cumulative �z 2 (1/ 2π) exp(−(t /2))dt. −∞
The prediction of the output probability for the input sample z∗ can be obtained by integration over the latent function f∗ as follows: � (5) p(a∗ = +1|D, z∗ ) = p(a∗ |f∗ )p(f∗ |D, z∗ )df∗ . The second term of the integral represents the distribution of the latent variable corresponding to z. It can be obtained by marginalization of f = [f1 f2 · · · fn ]. p(f∗ |D, z∗ ) =
�
p(f∗ |Z, z∗ , f )p(f |D)df .
(6)
where p(f |D) is the posterior over the latent variables. This can be reformulated by the Bayes’ rule as
IFAC AAC 2016 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
p(f |D) = p(a|f )p(f |Z)/p(a|Z) n p(ai |fi ) p(f |Z)/p(a|Z), =
(7)
i=1
where p(a|f ) is the likelihood function and it can be expressed by using one of the squashing functions given in (4). p(a|Z) is the marginal likelihood and p(f |Z) is the GP prior of the latent variables. The GP prior over function can be written as 1 1 T −1 p(f |Z) = (f − µ) K (f − µ) , exp − 1 2 (2π)n/2 |K| 2 (8) where the mean µ is usually assumed to be zero and each term of a covariance matrix Kij is a function of zi and zj . This covariance function is the kernel k(zi , zj ). The covariance function is parameterized by the hyperparameters θ, which we can learn from the training data D. 4. EXPECTATION PROPAGATION
The EP approximation aims at minimizing locally the Kullback-Leibler divergence measure KL(p(f |D)�q(f |D)) between the true posterior p(f |D) and its Gaussian approximation q(f |D). By using the Bayes’ rule, the posterior p(f |D) an be rewritten as follows: n 1 p(f |D) = p(f |Z) p(ai |fi ), (9) L i=1
where the normalizing term L is the marginal likelihood p(a|Z). The idea of the EP approximation is to approximate the likelihood by a product of local likelihoods ti estimated through unnormalized functions in the latent variables fi ’s ˜ i, µ ˜ i N (fi |˜ p(ai |fi ) ∼ ˜i , σ ˜i2 ) = L µi , σ ˜i2 ), (10) = ti (fi |L
˜ i, µ ˜i , σ ˜i2 are site parameters. The approximated where L likelihood is equal to n n 2 ˜ ˜ ˜ i, (11) ti (fi |Li , µ ˜i , σ ˜i ) = N (˜ µ, Σ) L i=1
2 2 T ˜ = diag σ where µ ˜ = [µ1 µ2 · · · µn ] and Σ ˜1 σ ˜2 · · · σ ˜n2 The Gaussian approximation q(f |D) to the posterior becomes n 1 ˜i, µ q(f |D) = ti (fi |L ˜i , σ ˜i2 ) = N (˜ µ, Σ) p(f |Z) LEP i=1 ˜ −1 µ ˜ −1 )−1 , with µ = ΣΣ ˜ and Σ = (K−1 + Σ
The EP approximation estimates the each ti approximation iteratively using the cavity distribution as ˜j, µ q−i (fi ) ∝ p(f |Z) tj (fi |L ˜j , σ ˜j2 )dfj , (13) j�=i
where the subscript −i denotes the case −i is removed. The estimation process is performed by repeating four steps on each training sample until convergence. (i) The cavity distribution is computed as 2 q−i (fi |zi , ai ) = N (fi |µ−i , σ−i )
2 (σi−2 µi − σ ˜i−2 µ ˜i ) with µ−i = σ−i
˜i−2 )−1 . (14) and σ−i = (σi−2 − σ (ii) The cavity distribution is combined with the true likelihood r(fi ) = q−i (fi |zi , ai )p(ai |fi ). (15) (iii) The Gaussian approximation to the non-Gaussian marginal is computed ˆ i N (ˆ qˆi (fi ) = L µi , σ ˆ2) ∼ = r(fi |zi , ai ), i
Since the integrals in (5) and (6) are not analytically tractable due to the non-Gaussian nature of the likehood terms, analytical approximation or sampling method (e.g. Monte Carlo sampling) is applied. The common analytical approximation methods are Laplace and EP approximations. In this section we summarize the EP algorithm in Rasmussen et al. (2006) and Minka (2001). In general, GPC with EP approximation (GPC-EP) show better results than GPC with Laplace approximation (GPC-LP) in our DoE strategy.
i=1
789
(12)
and LEP = p(a|Z) is the EP approximation to the normalizing term L. 803
2 ai σ ˆ−i N (zi ) , 2 Φ(zi ) 1 + σ−i 4 σ ˆ−i N (zi ) N (zi ) 2 z , σ ˆi2 = σ−i + − i 2 )Φ(z ) (1 + σ−i Φ(zi ) i 2 . with zi = ai µ−i / 1 + σ−i (16)
ˆ = Φ(zi ), µ L ˆi = µ−1 +
(iv) The parameters of the approximation ti are updated to make the posterior have the desired marginal as follows: −2 −2 −1 µ ˜i = σ ˜i2 (ˆ σi−2 µ ˆi − σ−i µ ˜−i ), σ ˜i = (ˆ σi−2 − σ−i ) , 2 √ ˜i ) (µ−i − µ 2 exp ˜i = L ˆ i 2π σ 2 + σ L ˜−i −i 2 +σ 2 ) . (17) 2(σ−i ˜−i 4.1 Predictions
The prediction of the sample z∗ is evaluated by using the Gaussian approximation p(a∗ = +1|D, z∗ ) ∼ (18) = q(a∗ = +1|D, z∗ ) = p(a∗ |f∗ )q(f∗ |D, z∗ )df∗ ,
where q(f∗ |D, z∗ ) is a Gaussian with mean and variance given as follows: ˜ −1 µ ˜ (19) µ∗ = k(z∗ )T (K + Σ) ˜ −1 k(z∗ ), (20) σ∗2 = k(z∗ , z∗ ) − k(z∗ )T (K + Σ)
and k(z∗ ) = [k(z1 , z∗ ) k(z2 , z∗ ) · · · k(zn , z∗ )]T is a vector of kernel covariances between z∗ and all the training samples. If the probit form is taken for the squashing function, then prediction can be evaluated analytically as µ∗ . (21) q(a∗ = +1|D, z∗ ) = Φ 1 + σ∗2 4.2 Hyperparameters
The choice of the hyperparameters θ is important, since it has a direct impact on the accuracy of the GPC. The hyperparameters can be determined by the maximization
IFAC AAC 2016 790 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
problem of the marginal likelihood with respect to θ. The EP approximation to the marginal likelihood can be found from the normalization LEP = q(a|Z, θ) � n � ˜i, µ ti (fi |L ˜i , σ ˜i2 )df . = p(a|Z, θ)
(22)
i=1
After some algebraic manipulations, we can obtain the marginal log likelihood as 1 ˜ − 1µ ˜ −1 µ ˜T (K + Σ) log LEP = − log |K + Σ| ˜ 2 2 n n � ai mu−i 1 � 2 2 + log Φ � log(σ−i +σ ˜−i ) + 2 2 1 + σ−i i=1 i=1 +
n � (µ−i − µ ˜i )2 2 2 ). 2(σ−i + σ−i i=1
(23)
Fig. 2. Input design based on GPC. This approach can determine the next measurement ”point” near the boundaries.
5. METHOD OF INPUT DESIGN The Gaussian process classifier can generate not only a predictive probability model but also an approximate predictive variance model. In this section we propose an efficient method of input design using this approximate predictive variance model generated by GPC. Our proposed method is based on the idea that the predictive variance will have large value near the boundaries. In this paper we obtain output data of outside of the admissible domain as ai = −1 and output data of inside of the admissible domain as ai = +1. In GPC, we consider the probability of belonging to the class label a∗ = +1, p(a∗ = +1|D, z∗ ). We summarize the input design method in Algorithm 1, where l is the number of additional inputs. This algorithm is also illustrated by Fig. 2. The parameter P ∈ R is the probability threshold value. We determine new measurement point only in the domain where predictive probabilities are lower than P . If large P is given, a violation of the admissible domain likely to be caused although an aggressive search may be performed. The parameter Ph ∈ R is a probability threshold for the boundary model.
6. NUMERICAL RESULT In this section we demonstrate the effectiveness of the input design based on GPC-EP. First we evaluate the performance of two common GPC approximation methods using a simple boundary. Then we apply the proposed method to engine benchmark problem and generate the static boundary models. To introduce the framework of automatic relevance determination (ARD), we consider the kernel function as � � Nz 1 � k(zi , zj ) = θ0 exp − (24) θm (zim − zjm )2 , 2 m=1
where zim is the mth element of zi . In this framework, we deal with hyperparameters for each input variables. These numerical simulation are performed on Matlab R2015b using GPML Matlab Code available at Rasmussen et al. (2006). 6.1 Identification of the unit circle boundary
Algorithm 1 Input Design Based on GPC 1: Step 0: Set some classified points. 2: Loop: For k = 0, 1, . . . , l − 1, perform: Step 1 Generate a predictive probability model by GPC algorithm (Section 3 and Section 4) and estimate the admissible domain Hk := {z∗ |p(a∗ = +1|D, z∗ ) < P }. Step 2 Determine a next measurement point zk+1 = arg max σ∗ (e.g. in (20)). z∗ ∈Hk
Step 3 Apply the point zk+1 to the system and obtain an additional output data ak . 3: Output: Generate a predictive probability model by GPC and estimate the boundary model h(z∗ ) =: p(a∗ = +1|D, z∗ ) − Ph .
804
Fig. 5. Initial condition of unit circle modeling. We set some classified points. Since GPC-EP commonly shows good performance for binary classification problems, GPC-EP is also expected to be suitable for static boundary modeling. Here we
IFAC AAC 2016 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
791
Fig. 3. (a) GPC with Laplace approximations, and (b) GPC with EP approximations (10 additional inputs). 2σ of GPC-EP provides large score near the boundary.
Fig. 4. Example of boundary identification with Algorithm 1. In this simulation, the EGR valve was closed, the phase of VVT is fixed at 20 degree and the parameters P and Ph are 0.30. (a) Input data after 200 iterations. These points are placed near the boundaries. (b) Identified static boundary models after 100, 150 and 200 iterations. This figure also shows the 900 test data for this validation.
805
IFAC AAC 2016 792 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
compare the performance between GPC-EP and GPCLP (Rasmussen et al. (2006) and Bishop (2006)). Both approximation methods are shown with same value of the hyperparameters. The initial condition is shown in Fig. 5. We set some classified points on input limitaion and near the boundary. Fig. 3 illustrates the results of GPCs performance after 10 iterations in Algorithm 1. Both methods can estimate the boundary by probability model, but GPC-EP shows more effective models. The most marked difference between two methods is the predictive variance model. The predictive 2σ∗ model provides large score only near the true boundary. This feature of GPC-EP is suitable for our DoE strategy since additional inputs are computed by predictive variance model. It can be expected that we place additional inputs near the true boundaries efficiently and generate an accurate boundary model with small amount of data. 6.2 Engine benchmark problem We apply the proposed method with GPC-EP to the engine benchmark problem in Watanabe et al. (2014). The detail of engine benchmark problem is not discussed here, but we consider the decision vector of static boundary model hS (z) consists of 5 variables, that is the engine speed, the opening angle of the throttle valve, the phase shift of VVT, the spark advance, and the lift of EGR valve. For visualization, we show the case of 3 inputs and 1 output (see Fig. 4). In this example, we generate the static boundary model corresponding to the engine speed, the opening angle, and spark advance. The engine operating conditions are divided into 100 × 100 × 100 grids. It can be seen that we determine the inputs near the boundaries and generate non-convex boundary models. Please note that the additional inputs are not placed far from the admissible domain. We also verify the accuracy of
Fig. 6. Classification accuracy between the identified boundary model and the test data. The accuracy is validated with True Positive (TP) and True Negative (TN). The False Positive (FP) is 0.3-1.0 % of the test data set. identified boundary models using the test data (see Fig. 6). The 900 test data on the engine operating conditions are acquired in advance. As the additional inputs increase, the number of correct decision (TP or FP) tends to increase between 50 and 200 additional inputs. On the other hand, the increase seems reach the limit in more than 200 inputs case. In this example, the limit of classification accuracy is about 97%. The incorrect decision rate in abnormal operation domain (FN) is very low (between 0.3% and 1.0%). As the additional inputs increase, the number of incorrect decision decreases. 806
The proposed algorithm can be also directly applied to high dimensional cases. On the otherhand, large data sets may cause computational problems, which depend on the number of data points. For a situation like that, we need to apply an approximation method such as SPGPs (Snelson et al. (2006)). 7. CONCLUSION In this paper, we proposed a DoE strategy based on GPCEP. This strategy describes accurate boundary models and also provided an efficient input design algorithm using the generated model. This method obtaines additional inputs as points, therefore the calibration time is expected to be much shorter than the common methods and the approach based on SVM. The effectiveness of the proposed method was demonstrated by a simple model and the JSAE and SICE engine benchmark problem. This benchmark has abnormal engine operation models such as knocking and misfiring. In the experimental result GPC-EP provided better performance than GPC-LP and input design algorithm based on GPC-EP can provide accurate boundary models of the engine simulator. Future work will be focused on the identification of dynamical boundary models. Our goal is to control an engine robustly near the boundaries by using the boundary models as inequality constraints. REFERENCES S. Watanabe, and A. Ohata. Benchmark Problem for Near Boundary Operation Control for Automotive Engine. In Proceeding of 6th Conference on Simulation and Testing for Automotive Electronics, pp. 101-110, 2014 N. Kieft. Evaluation of different design space description methods for analysing combustion engine operation limits. Leiden Institute of Advanced Computer Science (LIACS), Faculty of Science, Leiden University, 2014. P. Renniger, and M. Aleksandrov. Rapid hull determination: a new method to determine the design space for model based approaches. Design of Experiments (DoE) in Engine Development II. Expert-Verlag, Renningen, 2005. F. Pukelsheim. Optimal Design of Experiments. Classics in Applied Mathematics. SIAM (Society for Industrial and Applied Mathematics), 2006. G. Kampmann, N. Kieft, and O. Nelles. Support Vector Machines for Design Space Exploration. In Proceeding of the World Congress on Engineering and Computer Science 2012, volume 2, pp. 1116-1121, 2012. C. E. Rasmussen, and C. K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006. C. M. Bishop. Pattern recognition and machine learning. Springer-Verlag New York, LLC, 2006. T. P. Minka. A family of algorithms for approximate Bayesian inference. Ph. D. thesis, MIT, 2001. C. E. Rasmussen, and C. K. I. Williams. Gaussian Processes for Machine Learning Software. Online, available at http://www.gaussianpro cess.org/gpml/code/matlab/doc/. E. Snelson, Z. Ghahramani. Sparse Gaussian Processes using Pseudo-inputs. Advances in Neural Information Processing Systems 18, 2006.
ScienceDirect
IFAC-PapersOnLine 49-11 (2016) 787–792
Identification of Static Boundary Model Identification Identification of of Static Static Boundary Boundary Model Model Based on Gaussian Process Classification Based on Gaussian Process Classification Based on Gaussian Process Classification ∗∗∗ H. Oyama ∗∗ M. Yamakita ∗∗ K. Sata ∗∗ ∗∗ A. Ohata ∗∗∗ ∗ M. Yamakita ∗ ∗ K. Sata ∗∗ ∗∗ A. Ohata ∗∗∗ ∗∗∗ H. Oyama ∗ H. Oyama M. Yamakita K. Sata A. Ohata H. Oyama M. Yamakita K. Sata A. Ohata ∗ ∗ Dept. of Mechanical and Control Engineering, Tokyo Institute of ∗ of Mechanical and Control Engineering, Tokyo Institute of ∗ Dept. Dept. and Control Engineering, Technology, Tokyo, Japan Dept. of of Mechanical Mechanical and(e-mail: [email protected]) Engineering, Tokyo Tokyo Institute Institute of of Technology, Tokyo, Japan (e-mail: [email protected]) ∗∗ Technology, Tokyo, Japan (e-mail: [email protected]) Unit Management System Development Div., TOYOTA Technology, Japan (e-mail: [email protected]) ∗∗ Advanced Tokyo, ∗∗ Advanced Unit Management System Development Div., TOYOTA ∗∗ Advanced MOTOR Unit Management Management System Development Development Div., TOYOTA TOYOTA CORPORATION, Shizuoka, Japan Advanced Unit System Div., MOTOR CORPORATION, Shizuoka, Japan MOTOR CORPORATION, Shizuoka, Japan (e-mail: kota [email protected]) MOTOR CORPORATION, Shizuoka, Japan (e-mail: kota [email protected]) ∗∗∗ (e-mail: [email protected]) INC., kota Higashifuji Technical Center, Toyota Motor (e-mail: kota [email protected]) ∗∗∗ TECHNOVA ∗∗∗ TECHNOVA INC., Higashifuji Technical Center, Toyota Motor ∗∗∗ Corporation, TECHNOVA INC., Higashifuji Technical Center, Shizuoka, Japan (e-mail: [email protected]) TECHNOVA INC., Higashifuji Technical Center, Toyota Toyota Motor Motor Corporation, Shizuoka, Japan (e-mail: [email protected]) Corporation, Shizuoka, Shizuoka, Japan Japan (e-mail: (e-mail: [email protected]) [email protected]) Corporation,
Abstract: Due to the needs for CO2 emission reduction, the control region for automotive Abstract: Due to for CO reduction, the control region for 2 emission Abstract: Due to the the needs needs for CO emission reduction, theand control region for automotive automotive 2 engine system is pushed towardfor theCO boundary between normal abnormal engine operations Abstract: Due to the needs reduction, the control region for automotive 2 emission engine system is pushed toward the boundary between normal and abnormal engine operations engine system is pushed toward the boundary between normal and abnormal engine operations such assystem knocking and misfiring. When an engine is operated outside of admissible region, the engine is pushed toward the boundary between normal and abnormal engine operations such as knocking and misfiring. When an engine is operated outside of admissible region, the such as knocking and misfiring. When an engine is operated outside of admissible region, the engine or the catalyst is seriously damaged or passengers feel very uncomfortable. Therefore such as knocking and misfiring. When an engine is operated outside of admissible region, the engine or the catalyst is seriously damaged or passengers feel very uncomfortable. Therefore engine or the the catalyst catalyst is boundary seriously damaged damaged or passengers passengers feel very uncomfortable. Therefore the identification of theis is highly demanded. In feel this very paperuncomfortable. we propose a Therefore design of engine or seriously or the identification of the boundary is highly demanded. In this paper we propose aa design of the identification the highly demanded. this we propose experiment (DoE) of strategy based onis Gaussian processIn classifier (GPC) the expectation the identification of the boundary boundary isthe highly demanded. In this paper paper we with propose a design design of of experiment (DoE) strategy based on the Gaussian process classifier (GPC) with the expectation experiment (DoE) strategy based on the Gaussian process classifier (GPC) with the expectation propagation(DoE) (EP) strategy algorithm. In the experimental result, classifier the proposed algorithm can identify experiment based on the Gaussian process (GPC) with the expectation propagation (EP) algorithm. In the experimental result, the proposed algorithm can identify propagation (EP) algorithm. In result, can the boundaries of an engine benchmark problem, provided byproposed the joint algorithm research committee of propagation (EP) algorithm. In the the experimental experimental result, the the proposed algorithm can identify identify the boundaries of an engine benchmark problem, provided by the joint research committee of the Society boundaries of an an engine engine benchmark problem, provided by the joint joint and research committee of of Automotive Engineers (JSAE) and Society of by Instrument Control Engineers the boundaries of benchmark problem, provided the research committee of the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers (SICE). the Society of Automotive Engineers (JSAE) and Society of Instrument and Control Engineers (SICE). (SICE). (SICE). © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Gaussian processes; Classification; Identification algorithms; Engine modelling; Keywords: Gaussian Gaussian processes; processes; Classification; Classification; Identification Identification algorithms; algorithms; Engine Engine modelling; Keywords: Engine control; Model based control; Statistical design. Keywords: Gaussian processes; Classification; Identification algorithms; Engine modelling; modelling; Engine control; control; Model based control; control; Statistical design. Engine Model based Statistical design. Engine control; Model based control; Statistical design. 1. INTRODUCTION region are the convex hull based method (e.g. Rapid Hull 1. INTRODUCTION INTRODUCTION region hull method (e.g. Rapid Hull 1. region are are the the convex convex hull based based method (e.g. Rapidand Hulla Determination Algorithm, Renniger et al.(e.g. (2005)) 1. INTRODUCTION region are the convex hull based method Rapid Hull Determination Algorithm, Renniger et al. (2005)) and a Determination Algorithm, Renniger et al. (2005)) and a way that minimize the variance of the parameters (e.g. Determination Algorithm, Renniger et al. (2005)) and a To solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. DDTo solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. Doptimal design, Pukelsheim (2006)). However, the region To solve the energy issue and the global environmental way that minimize the variance of the parameters (e.g. Dproblem, automakers have developed high functional enTo solve the energy issue the global design, Pukelsheim (2006)). However, the region problem, automakers have and developed high environmental functional enen- optimal optimal design, the of normal enginePukelsheim operations (2006)). seems toHowever, be non-convex and problem, automakers have developed high functional optimal design, Pukelsheim (2006)). However, the region region gines. Although the engine control system has been comproblem, automakers have developed high functional enof normal engine operations seems to be non-convex and gines. Although Although the the engine engine control control system system has has been been comcom- of of normal engine operations seems to be non-convex and in the latter case, we usually deal with quite simple models normal engine operations seems to be non-convex and gines. plicated, the development becoming shorter to in the latter case, we usually deal with quite simple models gines. Although the engine period controlis system has been complicated, the development period is becoming shorter to in the latter case, we usually deal with quite simple models such as polynomial models. Therefore the basic identified plicated, the development period is becoming shorter to in the latter case, we usually deal with quite simple models meet market It is difficult to isdevelop products satisplicated, the needs. development period becoming shorter to such as polynomial models. Therefore the basic identified meet market market needs. It is is difficult difficult to develop develop products satissuch as polynomial Therefore the identified models never fit models. the actual boundaries. To address this meet needs. It to products satissuch as will polynomial models. Therefore the basic basic identified fying requirements in the short term through conventional meet market needs. It is difficult to develop products satismodels will never fit the actual boundaries. To address this fying requirements in the short term through conventional models will never fit the actual boundaries. To address this issue, a DoE strategy based on support vector machine will never fit the actual boundaries. To address this fying requirements in the short through conventional development environments. Forterm such a background the models fying requirements in the short term through conventional issue, a DoE strategy based on support vector machine development environments. For such a background the issue, a DoE strategy based on support vector machine (SVM) is proposed in Kampmann et al. (2012). This development environments. For such a background the issue, a DoE strategy based on support vector machine model based development (MBD), which isbackground a development development environments. For such a the (SVM) is proposed in Kampmann et al. (2012). This model based based development development (MBD), (MBD), which which is is aa development development (SVM) Kampmann et al. method is canproposed find andin not only but This also model is proposed indescribe Kampmann et convex, al. (2012). (2012). This approach using virtual models and simulators, attracts (SVM) model based development (MBD), which is a development method can find and describe not only convex, but also approach using virtual models and simulators, attracts method can find and describe not only convex, but also non-convex boundaries. On the other hand, this strategy method can find and describe not only convex, but also approach using virtual models and simulators, attracts much attention as one solution to the problem. approach using virtual models and simulators, attracts non-convex boundaries. On the other hand, this strategy much attention as one solution to the problem. non-convex boundaries. On the other hand, this strategy can determine only search directions for additional inputs. much attention as one solution to the problem. non-convex boundaries. On the other hand, this strategy much attention as one solution to the problem. determine only directions for additional inputs. To reduce CO2 emission, the control region for automotive can canfind determine only search search directions forslowly additional inputs. To boundaries, inputs directions are changed and usually can determine only search for additional inputs. To reduce reduce CO CO2 emission, the control region for automotive To find boundaries, inputs are changed slowly and usually To emission, the control region for automotive 2 engine system is pushed toward the boundary between To find boundaries, inputs are changed slowly and usually To reduce CO emission, the control region for automotive in small discrete steps. 2 To find boundaries, inputs are changed slowly and usually engine is toward the between small discrete steps. engine system system is pushed pushed toward the boundary boundary between in normal and abnormal engine operations such as knocking in small discrete steps. engine system is pushed toward the boundary between in small discrete steps. normal and abnormal engine operations such as knocking In this paper we a DoE strategy based on the normal and abnormal engine operations such as knocking and misfiring. The abnormal operation may damage the In this paper we propose normal and abnormal engine operations such as knocking propose aa DoE strategy based on the and misfiring. The abnormal operation may damage the In this paper we propose DoE on Gaussian process classifier with thebased expectation and misfiring. The abnormal operation may damage the In this paper we propose a(GPC) DoE strategy strategy based on the the engine seriously. Therefore the control method satisfying and misfiring. The abnormal operation may damage the Gaussian process classifier (GPC) with the expectation engine Therefore the method Gaussian process with the (EP) classifier algorithm(GPC) in Rasmussen etexpectation al. (2006) engine seriously. seriously. Therefore theis control control method satisfying satisfying Gaussian process classifier (GPC) with the expectation boundary constrains robustly highly demanded. To de- propagation engine seriously. Therefore the control method satisfying (EP) in et (2006) boundary constrains constrains robustly robustly is is highly highly demanded. demanded. To To dede- propagation propagation (EP) algorithm algorithm in Rasmussen Rasmussen et al. al.accurate (2006) and Minka (2001). This strategy can describe (EP) algorithm in Rasmussen et al. (2006) boundary sign such aconstrains controllerrobustly throughisMBD, course, accurate boundary highlyof demanded. To de- propagation and Minka (2001). This strategy can describe accurate sign such a controller through MBD, of course, accurate and Minka (2001). This strategy can describe accurate boundary models and provide an efficient input sign such a controller through MBD, of course, accurate and Minka (2001). This strategy can describe accurate numerical of these boundaries (Watan- boundary models and provide an efficient input design sign such amodels controller through MBD,are of required course, accurate design numerical models of these boundaries are required (Watanboundary models provide an design algorithm using theand generated That is,input this method numerical models of these boundaries are required (Watanboundary models and providemodel. an efficient efficient input design abe et al. (2014) and Kieft (2014)). To identify boundaries numerical models of these boundaries are required (Watanalgorithm using the generated model. That is, this method abe et al. (2014) and Kieft (2014)). To identify boundaries algorithm using the generated model. That is, this method can describe non-convex region and can determine the using the generated model. That is, this method abean etengine al. (2014) (2014) and Kieftthe (2014)). To method identify needs boundaries at testand bench, common a lot algorithm abe et al. Kieft (2014)). To identify boundaries describe non-convex region and can determine the at an an engine engine test test bench, bench, the the common common method method needs needs aa lot lot can can describe non-convex region and can determine the next measurement point near the boundaries. By this at can describe non-convex region and can determine the of time and cost. The main objective method of this paper to next measurement point near the boundaries. By this at an engine test bench, the common needs ais lot of time and cost. The main objective of this paper is to next measurement point near the boundaries. By this direct determination of additional inputs as points, the of time and cost. The main objective of this paper is to next measurement point near the boundaries. By this obtain models withofanthis efficient of timeaccurate and cost.boundary The main objective paper input is to direct determination of additional inputs as points, the obtain accurate boundary models with an efficient input direct determination of points, the calibration time is expected to beinputs much as direct determination of additional additional inputs asshorter. points, The the obtain accurate boundary models with an efficient input design algorithm. obtain accurate boundary models with an efficient input calibration time is expected to be much shorter. The design algorithm. calibration time is expected to be much shorter. The validity of the proposed method is demonstrated by using design algorithm. algorithm. calibration time is expected to be much shorter. The design of the method is by using In the field of automotive engine design, we usually con- validity validity of benchmark the proposed proposedproblem method(Watanabe is demonstrated demonstrated by using an engine et al. by (2014)), validity of the proposed method is demonstrated using In the field of automotive engine design, we usually conan engine benchmark problem (Watanabe et al. (2014)), In the field of automotive engine design, we usually consider the static case. This means that we obtain an experian engine benchmark problem (Watanabe et al. (2014)), In the field of automotive engine design, we usually conprovided by the joint research committee of the Society of an engine benchmark problem (Watanabe et al. (2014)), sider static case. This that we an experiby the joint research committee of the Society of sider the the static case. This means means that we obtain obtain anThe experimental data after an engine reaches a steady state. ap- provided provided by the joint research of the Society of sider the static case. This means that we obtain an experiAutomotive Engineers (JSAE)committee and Society of Instrument provided by the joint research committee of the Society of mental data after an engine reaches a steady state. The apAutomotive Engineers (JSAE) and Society of Instrument mental data after an engine reaches a steady state. The approach data to determine inputs reaches is calleda design of experiment Automotive Engineers (JSAE) and Society of Instrument mental after an engine steady state. The apand Control (SICE). The empirical combustion Automotive Engineers (JSAE) and Society of Instrument proach to determine inputs is called design of experiment proach to inputs is of (DoE). The most common to design the control and and Control Control Engineers Engineers (SICE). (SICE). The The empirical empirical combustion combustion proach to determine determine inputsapproaches is called called design design of experiment experiment and Control Engineers (SICE). The empirical combustion (DoE). The most common approaches to design the control (DoE). The most common approaches to design the control (DoE). The most common approaches to design the control Copyright 2016 IFAC 801 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2016, IFAC (International Federation of Automatic Control) Copyright 2016 IFAC 801 Copyright © 2016 IFAC 801 Peer review© of International Federation of Automatic Copyright ©under 2016 responsibility IFAC 801Control. 10.1016/j.ifacol.2016.08.115
IFAC AAC 2016 788 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
profile with engine operating conditions and the knock model are implemented in this benchmark. The organization of this paper is as follows. In Section 2, we define the boundary modeling. We present the GPC and EP approximation algorithm in Section 3 and Section 4. In Section 5, we derive an input design algorithm based on GPC. Then in Section 6, we demonstrate the effectiveness of the proposed method by numerical simulations. Conclusions are given in Section 7. 2. BOUNDARY MODELING First of all, we explain about some definitions of the boundary modeling. We consider a discrete time nonlinear system represented by xk+1 = f (xk , uk , wk ), (1) yk = g(xk , uk , vk ),
Fig. 1. Example of boundary modeling strategy. To avoid engine damage, the inputs are changed slowly and usually in small discrete steps.
where x ∈ RN x , u ∈ RN u , and y ∈ RN y are the state vector, the input vector, and the output vector respectively. The subscript k denotes k-step value, w ∈ RN w and v ∈ RN v denote parameters of system uncertainty. We define the boundary model h(z) : RN z → R as follows: � > 0, outside of admissible domain on boundary (2) h(zk ) = 0, < 0, inside of admissible domain,
In this section we introduce Gaussian process classification in Rasmussen et al. (2006) and Bishop (2006). GPC is a Bayesian kernel classifier derived from Gaussian process priors. We assume that we have a training data set D = {(zi , ai )|i = 1, . . . , n}, where zi ∈ RD , ai ∈ {−1, 1}, and n denote input data, output data, and the number of samples respectively. The training set can be expressed as D = (Z, a) by using Z = [z1 z2 · · · zn ] and a = [a1 a2 · · · an ]. Given this training data set, the classification problem is to output the correct class label for a new input data z∗ . In GPC, we predict the class label of z∗ through the class posterior probability p(a∗ |D, z∗ ).
where zk ∈ RN z is the decision vector and we assume that H := {zk |h(zk ) < 0} becomes a connected set. 2.1 Static boundary modeling The goal of boundary modeling is to represent the true boundary numerically as h(z) in (2). However, in cases where systems are very complicated, it is difficult to determine the composition of the decision vector. The way how to acquire efficient data sets for identification is also an important problem. Therefore, in practice, we usually identify the boundary model using steady state data points to simplify the problems. In this paper, we define the modeling with steady state data set as static boundary modeling. In this � static�Tcase we can derive the decision vector as z = xT , uT . Here we define the steady state of the system x∞ ∈ RN x and constant inputs u∞ ∈ RN u , which satisfy x∞ = f (x∞ , u∞ , 0). Our objective of static boundary modeling is to identify the model as follows: �� �� � > 0, outside of admissible domain x∞ = 0, on boundary hS (3) u∞ < 0, inside of admissible domain. When elements of the state x∞ are determined only by constant input u∞ , we can omit the elements from decision vector z. The common approaches to identify the static boundary model are the convex hull based method (e.g. Rapid Hull Determination Algorithm, Renniger et al. (2005)) and Doptimal design in Pukelsheim (2006). These strategy can determine only search directions for additional inputs. To avoid engine damage, inputs are changed slowly and usually in small discrete steps (see Fig. 1). This procedure needs a lot of time and cost. 802
3. GAUSSIAN PROCESS CLASSIFICATION
In GPC, we assume that the probability of a class label depends on the value of some latent real value function f (z) ∈ R. That is, for binary classification, p(a = +1|z, f (z), D) = p(a = +1|f (z)). The probability of belonging to a class label ai = +1 for an input sample zi is monotonically increasing function related to the value fi (z) ∈ R of f (z). This monotonic relationship is defined according to a squashing function such as logistic and probit functions. For example, fi = f (zi ) can be derived as follows: 1 , logistic function p(ai = +1|fi ) = 1 + exp(−fi ) (4) Φ(fi ), probit function, where Φ normal function Φ(z) = √ is the cumulative �z 2 (1/ 2π) exp(−(t /2))dt. −∞
The prediction of the output probability for the input sample z∗ can be obtained by integration over the latent function f∗ as follows: � (5) p(a∗ = +1|D, z∗ ) = p(a∗ |f∗ )p(f∗ |D, z∗ )df∗ . The second term of the integral represents the distribution of the latent variable corresponding to z. It can be obtained by marginalization of f = [f1 f2 · · · fn ]. p(f∗ |D, z∗ ) =
�
p(f∗ |Z, z∗ , f )p(f |D)df .
(6)
where p(f |D) is the posterior over the latent variables. This can be reformulated by the Bayes’ rule as
IFAC AAC 2016 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
p(f |D) = p(a|f )p(f |Z)/p(a|Z) n p(ai |fi ) p(f |Z)/p(a|Z), =
(7)
i=1
where p(a|f ) is the likelihood function and it can be expressed by using one of the squashing functions given in (4). p(a|Z) is the marginal likelihood and p(f |Z) is the GP prior of the latent variables. The GP prior over function can be written as 1 1 T −1 p(f |Z) = (f − µ) K (f − µ) , exp − 1 2 (2π)n/2 |K| 2 (8) where the mean µ is usually assumed to be zero and each term of a covariance matrix Kij is a function of zi and zj . This covariance function is the kernel k(zi , zj ). The covariance function is parameterized by the hyperparameters θ, which we can learn from the training data D. 4. EXPECTATION PROPAGATION
The EP approximation aims at minimizing locally the Kullback-Leibler divergence measure KL(p(f |D)�q(f |D)) between the true posterior p(f |D) and its Gaussian approximation q(f |D). By using the Bayes’ rule, the posterior p(f |D) an be rewritten as follows: n 1 p(f |D) = p(f |Z) p(ai |fi ), (9) L i=1
where the normalizing term L is the marginal likelihood p(a|Z). The idea of the EP approximation is to approximate the likelihood by a product of local likelihoods ti estimated through unnormalized functions in the latent variables fi ’s ˜ i, µ ˜ i N (fi |˜ p(ai |fi ) ∼ ˜i , σ ˜i2 ) = L µi , σ ˜i2 ), (10) = ti (fi |L
˜ i, µ ˜i , σ ˜i2 are site parameters. The approximated where L likelihood is equal to n n 2 ˜ ˜ ˜ i, (11) ti (fi |Li , µ ˜i , σ ˜i ) = N (˜ µ, Σ) L i=1
2 2 T ˜ = diag σ where µ ˜ = [µ1 µ2 · · · µn ] and Σ ˜1 σ ˜2 · · · σ ˜n2 The Gaussian approximation q(f |D) to the posterior becomes n 1 ˜i, µ q(f |D) = ti (fi |L ˜i , σ ˜i2 ) = N (˜ µ, Σ) p(f |Z) LEP i=1 ˜ −1 µ ˜ −1 )−1 , with µ = ΣΣ ˜ and Σ = (K−1 + Σ
The EP approximation estimates the each ti approximation iteratively using the cavity distribution as ˜j, µ q−i (fi ) ∝ p(f |Z) tj (fi |L ˜j , σ ˜j2 )dfj , (13) j�=i
where the subscript −i denotes the case −i is removed. The estimation process is performed by repeating four steps on each training sample until convergence. (i) The cavity distribution is computed as 2 q−i (fi |zi , ai ) = N (fi |µ−i , σ−i )
2 (σi−2 µi − σ ˜i−2 µ ˜i ) with µ−i = σ−i
˜i−2 )−1 . (14) and σ−i = (σi−2 − σ (ii) The cavity distribution is combined with the true likelihood r(fi ) = q−i (fi |zi , ai )p(ai |fi ). (15) (iii) The Gaussian approximation to the non-Gaussian marginal is computed ˆ i N (ˆ qˆi (fi ) = L µi , σ ˆ2) ∼ = r(fi |zi , ai ), i
Since the integrals in (5) and (6) are not analytically tractable due to the non-Gaussian nature of the likehood terms, analytical approximation or sampling method (e.g. Monte Carlo sampling) is applied. The common analytical approximation methods are Laplace and EP approximations. In this section we summarize the EP algorithm in Rasmussen et al. (2006) and Minka (2001). In general, GPC with EP approximation (GPC-EP) show better results than GPC with Laplace approximation (GPC-LP) in our DoE strategy.
i=1
789
(12)
and LEP = p(a|Z) is the EP approximation to the normalizing term L. 803
2 ai σ ˆ−i N (zi ) , 2 Φ(zi ) 1 + σ−i 4 σ ˆ−i N (zi ) N (zi ) 2 z , σ ˆi2 = σ−i + − i 2 )Φ(z ) (1 + σ−i Φ(zi ) i 2 . with zi = ai µ−i / 1 + σ−i (16)
ˆ = Φ(zi ), µ L ˆi = µ−1 +
(iv) The parameters of the approximation ti are updated to make the posterior have the desired marginal as follows: −2 −2 −1 µ ˜i = σ ˜i2 (ˆ σi−2 µ ˆi − σ−i µ ˜−i ), σ ˜i = (ˆ σi−2 − σ−i ) , 2 √ ˜i ) (µ−i − µ 2 exp ˜i = L ˆ i 2π σ 2 + σ L ˜−i −i 2 +σ 2 ) . (17) 2(σ−i ˜−i 4.1 Predictions
The prediction of the sample z∗ is evaluated by using the Gaussian approximation p(a∗ = +1|D, z∗ ) ∼ (18) = q(a∗ = +1|D, z∗ ) = p(a∗ |f∗ )q(f∗ |D, z∗ )df∗ ,
where q(f∗ |D, z∗ ) is a Gaussian with mean and variance given as follows: ˜ −1 µ ˜ (19) µ∗ = k(z∗ )T (K + Σ) ˜ −1 k(z∗ ), (20) σ∗2 = k(z∗ , z∗ ) − k(z∗ )T (K + Σ)
and k(z∗ ) = [k(z1 , z∗ ) k(z2 , z∗ ) · · · k(zn , z∗ )]T is a vector of kernel covariances between z∗ and all the training samples. If the probit form is taken for the squashing function, then prediction can be evaluated analytically as µ∗ . (21) q(a∗ = +1|D, z∗ ) = Φ 1 + σ∗2 4.2 Hyperparameters
The choice of the hyperparameters θ is important, since it has a direct impact on the accuracy of the GPC. The hyperparameters can be determined by the maximization
IFAC AAC 2016 790 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
problem of the marginal likelihood with respect to θ. The EP approximation to the marginal likelihood can be found from the normalization LEP = q(a|Z, θ) � n � ˜i, µ ti (fi |L ˜i , σ ˜i2 )df . = p(a|Z, θ)
(22)
i=1
After some algebraic manipulations, we can obtain the marginal log likelihood as 1 ˜ − 1µ ˜ −1 µ ˜T (K + Σ) log LEP = − log |K + Σ| ˜ 2 2 n n � ai mu−i 1 � 2 2 + log Φ � log(σ−i +σ ˜−i ) + 2 2 1 + σ−i i=1 i=1 +
n � (µ−i − µ ˜i )2 2 2 ). 2(σ−i + σ−i i=1
(23)
Fig. 2. Input design based on GPC. This approach can determine the next measurement ”point” near the boundaries.
5. METHOD OF INPUT DESIGN The Gaussian process classifier can generate not only a predictive probability model but also an approximate predictive variance model. In this section we propose an efficient method of input design using this approximate predictive variance model generated by GPC. Our proposed method is based on the idea that the predictive variance will have large value near the boundaries. In this paper we obtain output data of outside of the admissible domain as ai = −1 and output data of inside of the admissible domain as ai = +1. In GPC, we consider the probability of belonging to the class label a∗ = +1, p(a∗ = +1|D, z∗ ). We summarize the input design method in Algorithm 1, where l is the number of additional inputs. This algorithm is also illustrated by Fig. 2. The parameter P ∈ R is the probability threshold value. We determine new measurement point only in the domain where predictive probabilities are lower than P . If large P is given, a violation of the admissible domain likely to be caused although an aggressive search may be performed. The parameter Ph ∈ R is a probability threshold for the boundary model.
6. NUMERICAL RESULT In this section we demonstrate the effectiveness of the input design based on GPC-EP. First we evaluate the performance of two common GPC approximation methods using a simple boundary. Then we apply the proposed method to engine benchmark problem and generate the static boundary models. To introduce the framework of automatic relevance determination (ARD), we consider the kernel function as � � Nz 1 � k(zi , zj ) = θ0 exp − (24) θm (zim − zjm )2 , 2 m=1
where zim is the mth element of zi . In this framework, we deal with hyperparameters for each input variables. These numerical simulation are performed on Matlab R2015b using GPML Matlab Code available at Rasmussen et al. (2006). 6.1 Identification of the unit circle boundary
Algorithm 1 Input Design Based on GPC 1: Step 0: Set some classified points. 2: Loop: For k = 0, 1, . . . , l − 1, perform: Step 1 Generate a predictive probability model by GPC algorithm (Section 3 and Section 4) and estimate the admissible domain Hk := {z∗ |p(a∗ = +1|D, z∗ ) < P }. Step 2 Determine a next measurement point zk+1 = arg max σ∗ (e.g. in (20)). z∗ ∈Hk
Step 3 Apply the point zk+1 to the system and obtain an additional output data ak . 3: Output: Generate a predictive probability model by GPC and estimate the boundary model h(z∗ ) =: p(a∗ = +1|D, z∗ ) − Ph .
804
Fig. 5. Initial condition of unit circle modeling. We set some classified points. Since GPC-EP commonly shows good performance for binary classification problems, GPC-EP is also expected to be suitable for static boundary modeling. Here we
IFAC AAC 2016 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
791
Fig. 3. (a) GPC with Laplace approximations, and (b) GPC with EP approximations (10 additional inputs). 2σ of GPC-EP provides large score near the boundary.
Fig. 4. Example of boundary identification with Algorithm 1. In this simulation, the EGR valve was closed, the phase of VVT is fixed at 20 degree and the parameters P and Ph are 0.30. (a) Input data after 200 iterations. These points are placed near the boundaries. (b) Identified static boundary models after 100, 150 and 200 iterations. This figure also shows the 900 test data for this validation.
805
IFAC AAC 2016 792 June 19-23, 2016. Norrköping, Sweden
H. Oyama et al. / IFAC-PapersOnLine 49-11 (2016) 787–792
compare the performance between GPC-EP and GPCLP (Rasmussen et al. (2006) and Bishop (2006)). Both approximation methods are shown with same value of the hyperparameters. The initial condition is shown in Fig. 5. We set some classified points on input limitaion and near the boundary. Fig. 3 illustrates the results of GPCs performance after 10 iterations in Algorithm 1. Both methods can estimate the boundary by probability model, but GPC-EP shows more effective models. The most marked difference between two methods is the predictive variance model. The predictive 2σ∗ model provides large score only near the true boundary. This feature of GPC-EP is suitable for our DoE strategy since additional inputs are computed by predictive variance model. It can be expected that we place additional inputs near the true boundaries efficiently and generate an accurate boundary model with small amount of data. 6.2 Engine benchmark problem We apply the proposed method with GPC-EP to the engine benchmark problem in Watanabe et al. (2014). The detail of engine benchmark problem is not discussed here, but we consider the decision vector of static boundary model hS (z) consists of 5 variables, that is the engine speed, the opening angle of the throttle valve, the phase shift of VVT, the spark advance, and the lift of EGR valve. For visualization, we show the case of 3 inputs and 1 output (see Fig. 4). In this example, we generate the static boundary model corresponding to the engine speed, the opening angle, and spark advance. The engine operating conditions are divided into 100 × 100 × 100 grids. It can be seen that we determine the inputs near the boundaries and generate non-convex boundary models. Please note that the additional inputs are not placed far from the admissible domain. We also verify the accuracy of
Fig. 6. Classification accuracy between the identified boundary model and the test data. The accuracy is validated with True Positive (TP) and True Negative (TN). The False Positive (FP) is 0.3-1.0 % of the test data set. identified boundary models using the test data (see Fig. 6). The 900 test data on the engine operating conditions are acquired in advance. As the additional inputs increase, the number of correct decision (TP or FP) tends to increase between 50 and 200 additional inputs. On the other hand, the increase seems reach the limit in more than 200 inputs case. In this example, the limit of classification accuracy is about 97%. The incorrect decision rate in abnormal operation domain (FN) is very low (between 0.3% and 1.0%). As the additional inputs increase, the number of incorrect decision decreases. 806
The proposed algorithm can be also directly applied to high dimensional cases. On the otherhand, large data sets may cause computational problems, which depend on the number of data points. For a situation like that, we need to apply an approximation method such as SPGPs (Snelson et al. (2006)). 7. CONCLUSION In this paper, we proposed a DoE strategy based on GPCEP. This strategy describes accurate boundary models and also provided an efficient input design algorithm using the generated model. This method obtaines additional inputs as points, therefore the calibration time is expected to be much shorter than the common methods and the approach based on SVM. The effectiveness of the proposed method was demonstrated by a simple model and the JSAE and SICE engine benchmark problem. This benchmark has abnormal engine operation models such as knocking and misfiring. In the experimental result GPC-EP provided better performance than GPC-LP and input design algorithm based on GPC-EP can provide accurate boundary models of the engine simulator. Future work will be focused on the identification of dynamical boundary models. Our goal is to control an engine robustly near the boundaries by using the boundary models as inequality constraints. REFERENCES S. Watanabe, and A. Ohata. Benchmark Problem for Near Boundary Operation Control for Automotive Engine. In Proceeding of 6th Conference on Simulation and Testing for Automotive Electronics, pp. 101-110, 2014 N. Kieft. Evaluation of different design space description methods for analysing combustion engine operation limits. Leiden Institute of Advanced Computer Science (LIACS), Faculty of Science, Leiden University, 2014. P. Renniger, and M. Aleksandrov. Rapid hull determination: a new method to determine the design space for model based approaches. Design of Experiments (DoE) in Engine Development II. Expert-Verlag, Renningen, 2005. F. Pukelsheim. Optimal Design of Experiments. Classics in Applied Mathematics. SIAM (Society for Industrial and Applied Mathematics), 2006. G. Kampmann, N. Kieft, and O. Nelles. Support Vector Machines for Design Space Exploration. In Proceeding of the World Congress on Engineering and Computer Science 2012, volume 2, pp. 1116-1121, 2012. C. E. Rasmussen, and C. K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006. C. M. Bishop. Pattern recognition and machine learning. Springer-Verlag New York, LLC, 2006. T. P. Minka. A family of algorithms for approximate Bayesian inference. Ph. D. thesis, MIT, 2001. C. E. Rasmussen, and C. K. I. Williams. Gaussian Processes for Machine Learning Software. Online, available at http://www.gaussianpro cess.org/gpml/code/matlab/doc/. E. Snelson, Z. Ghahramani. Sparse Gaussian Processes using Pseudo-inputs. Advances in Neural Information Processing Systems 18, 2006.