Efficient real driving emissions calibration of automotive powertrains under operating uncertainties

The steady-state calibration of automotive powertrains is typically based on the assumption of one specific drive cycle and perfectly controllable operating conditions. During real operation, however, these assumptions are violated, which implies that the calibration might in fact not be optimal. Therefore, in order to achieve reliable performance in a real-world setting, these uncertainties have to have been considered already during the calibration process. In this article, a stochastic optimization approach that takes the mentioned operating uncertainties into account by including probability distributions of the disturbances, is suggested. Furthermore, an approximation is derived of the distribution of the optimization performance criterion that greatly reduces the computational load during optimization compared with Monte Carlo sampling. Simulation results show that the proposed probabilistic approach leads to lower expected values of emissions and consumption when compared with deterministic optimization approaches ignoring the stochastic influences.


Introduction
Car manufacturers face the challenge of combining both customer demand for high driving comfort and performance on the one hand, and increasingly rigorous environmental legislation on the other-see e.g. Langouët et al. (2008). Due to elaborate testing procedures on given drive cycles in laboratories under perfectly controllable and repeatable conditions, vehicles may seemingly meet legislative restrictions. However, when driven in the real world with therefore significantly varying operating conditions, fuel consumption and cumulative emissions differ notably from test bed results and manufacturers' specifications-see e.g. Mock et al. (2012). Yet, unforeseen occurrences during the operation and unknown or unmodelled behaviour of the drivetrain itself are generally not considered in the calibration and control of automotive powertrains.
Apart from that, the European Union (EU) has recently implemented new regulations that require more elaborate real-world testing for assessing real driving emissions (RDE). Tests are now to take place on the road under real driving and environmental conditions. This implies that basing, and possibly overfitting, the calibration on one laboratory drive cycle might lead to results that do not reflect the performance during real-world tests, jeopardizing the type-approval process. Additionally, supposedly identical vehicles might show different emission and consumption values because of deviations in series production and material ageing, which further calls into question laboratory values. Therefore, in order to take realistic cycle variability and operating uncertainties into account, this work suggests a probabilistic approach based on the stochastic optimization of fuel consumption and emissions in order to obtain more reliable powertrain performance data under a broad range of operating conditions. The calibration of powertrains is generally a broad field. Owing to growing technical complexity, model-based approaches are typically the preferred methods in order to circumvent massive manual tuning efforts-see e.g. Wong et al. (2018), Isermann and Sequenz (2016) and Kruse, Kurz, and Lang (2010). The most common application of model-based calibration is the determination of optimal engine control unit (ECU) parameters for gasoline and diesel engines. A basic engine calibration task typically consists of multiple steps-see Röpke and von Essen (2008): (i) defining efficient measurement experiments using design of experiments (DoE) methods; (ii) conducting measurements on a test bench; (iii) using the data to estimate a model of the engine; (iv) mathematical optimization of the engine parameters; (v) filling and smoothing of the ECU maps.
This basic calibration is usually based on steady-state measurements and models, as discussed in Roepke (2014). After basic calibration, steady-state maps can be further adjusted based on the transient behaviour of the engine-see Röpke and von Essen (2008). A full hardware-in-the-loop-based calibration workflow is presented in Lee et al. (2018). This article, however, is focused on steady-state based calibration.
An overview of various optimization formulations discussing local, global and mixed approaches is given in Castagné et al. (2008). In terms of modelling, DoE is a standard approach for measurement experiments. The goal is to design the experiments in such a way that maximal information is obtained for modelling while keeping the time-consuming measurement effort low. Common model types are polynomial models, neural networks, local model networks and Gaussian process models-see Roepke (2014), Kruse, Kurz, and Lang (2010) and Kianifar, Campean, andbreak Richardson (2013).
As mentioned above, one shortcoming of the existing solutions is the fact that uncertainties in the real-world operation of powertrains are usually not considered in the calibration process. There are, however, innumerable different uncontrollable, unknown and random influences that affect the performance of a powertrain-see e.g. Kokkolaras et al. (2005) and Wu (2012). Some external factors are e.g. the operating conditions like air temperature, air pressure and surface topography. Additionally, the driver himself exerts great influence on fuel consumption and emissions, which is especially relevant for drive cycle based engine calibration. The operation of the powertrain is also inherently uncertain, because the sensors and actuators that are deployed in the control of different parts of the power unit do not operate without errors, but only within certain tolerances-see e.g. Criens et al. (2015) and Tschanz et al. (2013).
Another source of uncertainty stems from the use of a specific drive cycle in the calibration process. One example, which is also part of the currently implemented EU regulations, is the 'worldwide harmonized light vehicle test cycle' (WLTC). In contrast to the past, EU regulations now additionally require real-world RDE-testing, which means that a vehicle is not only tested on test benches but on real roads. There are numerous criteria that such a test drive has to fulfil in order to qualify for valid RDE-testing-see European Parliament and Council (2016). For example, the requirements include specific acceleration characteristics, minimal and maximal stop times and a certain partition of the cycle into different speed segments. In order to maximize the probability of passing the typeapproval tests, it is advisable to take possible RDE cycles into consideration in the calibration process. Using only one specific cycle introduces the risk of overfitting the calibration to that cycle-see e.g. Wasserburger, Hametner, and Didcock (2020) and Eriksson and Nielsen (2014). Moreover, the RDEspecifications allow for quite a few variations between two admissible RDE cycles, which means that it is not reasonable to base the entire calibration on just one single cycle. Figure 1. The engine is subjected to various random disturbances. Therefore, the distribution of the engine's output has to be taken into consideration for optimization.
In order to deal with these issues, this article suggests including the involved uncertainties as random variables in the model used for optimization. As a result, the optimization objective becomes a random variable itself. The new optimization objective can then be defined, for instance, as the expected value. Alternatively, other risk measures like 'value-at-risk' can be considered as in Wasserburger, Hametner, and Didcock (2020) and Mourat, Eckstein, and Koch (2020). A different approach is presented in Cook and Jarrett (2018), where the optimization goal is to shape the cumulative distribution function of the considered quantity so that it matches a desired target.
In this work, as an example, the proposed expected value optimization is conducted for a basic steady-state engine map calibration of a diesel engine considering input and cycle uncertainties as illustrated in Figure 1. None the less, the suggested approach can also readily be applied to any other technical system faced with stochastic properties or disturbances.
In Wasserburger, Hametner, and Didcock (2020), stochastic optimization for engine calibration was discussed using a set of varied drive cycles and Monte Carlo sampling for the evaluation of the risk measures to be optimized. The present work further develops and extends that approach and introduces a few critical new developments: using a fully distribution based problem formulation, an efficient approximation method for the optimization objective is derived. By means of this approximation, the relevant statistics of the objective distribution can be estimated without the need for time-consuming Monte Carlo simulation, where all disturbances needed to be simulated individually, which makes the stochastic approach computationally feasible for practitioners. For that purpose, a multivariate probability distribution over the engine operating space needs to be modelled, based on a given set of drive cycles representing the car usage for which the engine is to be optimized. The distribution estimation is achieved by means of Gaussian mixture modelling (GMM)-see e.g. Bishop (2006) and Hastie, Tibshirani, and Friedman (2001). Another very common non-parametric method for estimating distributions is kernel density estimation. The joint distribution is then described by a superposition of the kernel functions, one for each data point-see Sheather (2004).
In this article, a set of RDE-compliant drive cycles is utilized in the modelling, which implies that the resulting calibration is optimal for this class of cycles. Furthermore, disturbances in the ECU parameters, which were ignored in Wasserburger, Hametner, and Didcock (2020), are considered as well. Overall, the article offers a holistic approach for the risk-conscious and robust calibration of technical systems that are exposed to different types of operational uncertainties.
This article is organized as follows: in Section 2, a more specific and formal problem definition as well as an outline of the proposed solution approach is given. In Section 3, the disturbance modelling, the approximation of the output distribution and its optimization are presented. The results for different optimization approaches are shown, compared and discussed in Section 4. Finally, some conclusions are drawn in Section 5.

Problem definition and solution outline
The goal in the discussed calibration task is to determine engine parameters u that minimize a weighted sum of integrated fuel consumption and emission. The integration is approximated by a weighted sum of representative operating points using steady-state fuel and emission models. The engine's operating space, denoted by O ⊂ R 2 , consists of all physically possible combinations of engine speed and load. For each such operating point x ∈ O the optimal parameter vector u(x) ∈ R n has to be determined. For each of the n parameters, there is a mapping j , stored in the ECU, relating the operating point with the optimal parameter value: The aim of this work is to determine these engine maps j , so that fuel consumption and emissions are minimized. One common approach is to minimize the objective evaluated on one drive cycle. For that purpose, a finite set of discrete operating points representing the cycle are selected, and the optimization is performed for each operating point individually-see e.g. Langouët et al. (2008), Castagné et al. (2008) and Isermann and Sequenz (2016). Afterwards, the engine maps are constructed based on the results of the individual optimizations. However, owing to various smoothness constraints on the engine maps, the results cannot just be interpolated, which means that the final maps generally deviate from the optimal values at the representative operating points. Therefore, in this work, the smoothness constraints are already considered in the optimization by adding a penalty term to the objective, which penalizes the second-order derivatives of the maps. This implies that the operating points cannot be minimized individually any more, but the advantage is that no post-processing of the engine maps is necessary. For a given drive cycle in the operating space O, a discrete representation is given by the pair D = (X , W), where X = {x 1 , . . . , x N } ⊂ O is the set of representative operating points and W = [w 1 , . . . , w N ] ∈ R N is the vector of corresponding weights, where each weight roughly corresponds to the time spent in the operating point. Let f f and f e denote fuel consumption and emissions. Both depend on the engine's operating point x i ∈ R 2 and the respective selection of the parameters Additionally, a penalty function g is used for enforcing the smoothness of the resulting engine maps by penalizing high second derivatives in the engine maps. The second derivatives are approximated by means of a difference approximation. For that purpose, the operating points are arranged in a graph where, for each pair of connected edges, the second derivative can be approximated by a linear combination of the engine parameters. By denoting the vector of all engine parameter vectors as such a linear combination can be written as a kū , for k = 1, . . . , K, where K is the number of edgepairs in the graph and a k is a row vector with the same number of elements asū. The penalty function can then be defined as For a specific drive cycle, the optimization objective is then defined bỹ where p f , p e and p s are the respective weights.
Minimizing the objective function (3) would mean that only one drive cycle and no disturbances were considered. In Wasserburger, Hametner, and Didcock (2020), a sample consisting of numerous drive cycles was taken into account simultaneously by evaluating (3) for each cycle and thereby approximating a distribution of the performance criterion (3). Based on that Monte Carlo simulation and the obtained estimated distribution, either the corresponding expected value, or other risk measures such as the conditional value-at-risk, were optimized.
As mentioned before, another important source of uncertainty related to the engine maps arises during the operation of the vehicle: the demanded ECU parameters u i defined by the engine maps are set by various control units. However, the involved sensors measuring the relevant states for the controllers and the actuators performing the control actions are not exact and operate only within certain tolerance bands. Furthermore, external effects can cause deviations from the desired parameter setpoints. Therefore, the parameter values that are realized during operation vary around the target value defined in the engine map. An unfavourable disturbance combination could thus unexpectedly lead to higher fuel consumption and emissions if this uncertainty is not accounted for in the calibration of the ECU, as demonstrated in e.g. Herrmann et al. (2012) and Tschanz et al. (2013).
One obvious drawback of a Monte Carlo simulation based approach is the considerably increased calculation time stemming from the numerous model evaluations. Taking into account even more uncertainties, such as those mentioned above, will further increase the computational burden. However, this may render the optimization intractable for manufacturers. Therefore, in order to address these issues, the present article introduces the following contributions.
(1) In addition to drive cycle variation, disturbances in the engine maps during operation are considered as well.
(2) A probabilistic representation of a set of drive cycles is designed using Gaussian mixture modelling.
(3) An efficient approximation of the distribution statistics is developed, eliminating the timeconsuming re-evaluation of the cost function for different disturbance realizations.
In order to tackle these tasks, the objective (3) is extended with additive disturbances in the operating points and in the engine parameters: The disturbances ε x i and ε u i are multidimensional random variables that follow distributions that will be specified later. The disturbance vector ε x i represents the uncertainty in the operating points stemming from the variability in the considered drive cycles, while ε u i is the disturbance of the ECU parameters. The vector εū is defined equivalently to (1). Since the cost function C(u 1 , . . . , u N ) in (4) is random itself, it cannot be optimized directly. Instead, the expected value can be minimized, leading to the following minimization task: Each u i ∈ R n consists of n = 6 parameters: air mass flow, boost pressure, exhaust gas recirculation, main timing, rail pressure and swirl position. Thus, the number of unknowns is n × N. The constraints for the optimization, symbolically represented by the set U, contain upper and lower bounds for the resulting engine maps. The upper and lower bound for each parameter also depend on the operating point, so each u i has different limit-vectors.
The question now is how to calculate the expected value in the objective of (5). The naive approach of simply simulating many realizations of different disturbances and drive cycles and taking the average as an approximation of the expected value is possible but computationally intensive: for N operating points, M 1 drive cycles and an ECU disturbance sample size M 2 , the number of evaluations of the models f f and f e equals N × M 1 × M 2 for one approximation of the expected value. This value is further increased drastically during the optimization, when the expected value, which is the minimization objective, has to be evaluated many times. In cases where f f and f e are somewhat complicated and time-consuming to evaluate, this leads to very long, potentially unacceptable computation times.
In order to remedy this issue, an approximation of the expected value using a second-order Taylor approximation is suggested instead. It is based on the following observation: for a twice continuously differentiable function f : R → R and a zero-mean random variable Y with variance σ 2 > 0, the expected value of a disturbed function evaluation at an operating point x can be approximated as follows: Since the first derivative disappears owing to E [Y] = 0, a second-order approximation is necessary. In order to approximate the objective function in (5), a multidimensional equivalent of (6) will be used. Instead of performing a Monte Carlo simulation for each iteration of the optimization solver, only a few function evaluations in order to approximate the derivatives are carried out, which greatly reduces the computational load. Note that, in higher dimensions, the relevant approximation formulas are much more elaborate than the illustrative example (6). In order to apply this approximation of the expected value, the distributions of the uncertainties ε x i and ε u i need to be known or estimated. In this article, zero-mean normal distributions for the engine parameter disturbances ε u i are assumed. The variances are based on expert knowledge and depend primarily on measurement and control accuracy, which are more or less known quantities. The uncertainties in the operating space ε x i are estimated based on a multitude of driving cycles, upon which the calibration is to be based. Using GMM, a probability distribution over the engine operating space is derived. The operating points are chosen to be the derived means of the contained normal distributions. Movement around those operating points are interpreted as the disturbances ε x i , which are consequently also normally distributed random variables with zero mean and a covariance structure stemming from the GMM. Using these distributions and the respective approximations, the optimization can be performed much faster while also considering simultaneous disturbances in the operating points and the ECU parameters.

Stochastic engine calibration
In this section, the stochastic set-up is presented and the objective function approximation is derived.

Drive cycles and Gaussian mixture modelling
The goal is to consider the variability of different drive cycles in the optimization without resorting to repeated evaluation of the complete cycles in order to estimate the expected value. For that purpose, a description of the distribution of the given sample of drive cycles in the operating space O that allows the evaluation of the entirety of the drive cycles using just a few operating points and deviations around them is needed. An advantageous formulation therefore consists of a number of centres, which can serve as operating points, and corresponding covariance structures around those centres. A possible model structure that meets these requirements is a GMM, which is basically a superposition of normal distributions. Let f i be the probability density functions of multivariate normal distributions with mean μ i ∈ R d and covariance matrix x,i ∈ R d×d . Then a GMM is defined by the density function with N i=1 π i = 1. A GMM can be used for modelling distributions but also for soft clustering, where each individual distribution describes one cluster and each data point is associated with a probability of belonging to a specific cluster. The a-priori probability is given by the positive, constant weights π i , which are also often called responsibilities in this context. For a more detailed introduction to GMMs, the reader is referred to Hastie, Tibshirani, and Friedman (2001) and Bishop (2006).
Here, a distribution over the two-dimensional engine operating space O made up by engine speed and load is to be obtained. Therefore, d = 2. The engine calibration is to be performed for RDEcompliant drive cycles, as discussed in Section 1. A large sample of generated RDE cycles totalling approximately 180 h are given in the time-velocity domain with a sampling time of 1 s. These cycles were generated using the algorithm suggested in , where models for the conditional probabilities of acceleration and deceleration are simulated in an iterative fashion leading to RDE-compliant velocity trajectories. Note that of course the suggested GMM and also the entire calibration process can be based on arbitrary sets of cycles. However, the example in this work is focused on engine calibration with the obligatory RDE testing in mind.
Since the cycles are given in the time-velocity domain, they first have to be transformed to the engine operating space, because the engine maps are defined on that space. This can be done either by means of complete vehicle simulations of various complexities or by simply assuming certain vehicle parameters and a gear selection strategy. The velocity and the gear choice then determine the engine speed. The engine torque can be calculated from the engine speed and the power consumption due to all tractive forces. Finally, the torque is transformed into the mean effective pressure, which is here the second component of the operating space, besides engine speed. An example of an RDE-compliant drive cycle in the velocity domain and its transformation into the operating space are depicted in Figure 2. There is one data point per second. However, as the focus of this work is the estimation and usage of the distribution and not the question of how to obtain drive cycles, the precise process of the cycle transformation and its analysis are considered to be beyond the present scope.
Transforming all cycles of the sample to the operating space leads to a large two-dimensional data set. This data set is now used for the estimation of a GMM. The calculation is based on maximum likelihood estimation using the expectation-maximization (EM) algorithm-see e.g. Bishop (2006). The EM algorithm, given initial values for component means, covariances and responsibilities, iterates over expectation and likelihood maximization steps until convergence. However, the algorithm might converge to a local minimum. That is why the algorithm is run multiple times from different sets of initial conditions in order to find the global optimum. The solution with the largest likelihood is then selected as the final solution. Alternative estimation methods that make use of entropy matching techniques and particle swarm optimization, which are also suitable for generalized Gaussian mixture distributions, can be found in Fan, Lin, and Wu (2008) and Fan and Lin (2009).
One degree of freedom in the estimation is the number of components, which has to be specified before the estimation. A review of available selection techniques is given in McLachlan and Rathnayake (2014). Among the most common methods are Akaike's and the Bayesian information criterion (AIC, BIC). BIC has been observed to be a reliable model selection tool for mixture models, while AIC sometimes overestimates the number of components-see McLachlan and Rathnayake (2014) and Keribin (2000). Estimating GMMs for different numbers of components and calculating the respective AIC and BIC values results in the curves depicted in Figure 3(a). BIC assumes its absolute minimum at N = 15 components. AIC keeps decreasing for larger models, but the improvements slow down a lot after N = 15, which also indicates that N = 15 is a sensible choice.
The resulting GMM with 15 components and its distribution function f are illustrated in the contour plot in Figure 3(b). The squares indicate the means μ i of the components and the sizes of the squares reflect their respective responsibilities π i . Also, for each component the covariance matrix x 2,i is estimated. As visible in Figure 3(b), the utilized sample of RDE cycles and the transformation to the operating space primarily lead to engine speeds between 2000 and 2500 revolutions per minute (rpm) and low loads. There is only a low probability that high loads will be contained in the data. These observations depend of course on the transformation algorithm, and especially the gear shifting strategy. Furthermore, the covariance matrices x,i are generally not diagonal, which means that engine speed and load are correlated within the components.
Regarding the evaluation of the objective in (5), the estimated GMM will be used in the following way. • The component centres μ i will serve as the operating points x i in (4).
• The component proportions π i will serve as weights w i . • The component covariance matrices x,i will be used in the distribution of ε This also nicely sums up the advantages of the suggested approach using GMM. Everything that is needed for the approximation of the optimization objective is immediately obtained from the modelling. This would not be the case if, for example, other clustering techniques or kernel density estimation were used.
Lastly, note that the ε x i are not strictly speaking disturbances. They just describe the probability of scattering around the calculated centres μ i of the drive cycle during a typical RDE-compliant drive. However, in order to keep the nomenclature simple, ε x i will be referred to as disturbances in the operating space.

Disturbance in the ECU parameters
The ECU parameters are disturbed by the random vectors ε u i , which stem from measurement uncertainty and control errors. Their distributions are not estimated here but assumed based on expert knowledge and tolerances of actuators and sensors specified by suppliers. The ε u i , i = 1, . . . , N, are hence assumed to be a family of independent, identically distributed n-dimensional normal distributions: with u being a diagonal matrix with entries σ 2 u,1 , . . . , σ 2 u,n . This means that uncertainties in one engine parameter do not influence the uncertainties in the other parameters and the disturbances arising at one operating point do not affect the performance at other operating points and they are equally distributed for all operating points. Additionally, it is assumed that all uncertainties are zeromean, which means that there is no systematic error assumed in the engine operation.

Approximation of the objective function
In this subsection an approximation of the expected value of (4) based on the relation (6) is derived. Note that the weights w i are replaced by the estimated GMM component proportions π i and the operating points x i are replaced with the GMM component centres μ i : It suffices to show the derivation of the expected value of one of the summands in the sum and of the penalty term. Therefore, a function f is considered as a placeholder for either f f or f e . Furthermore, the index i will be suppressed for clarity. As mentioned above, u i is a function of the operating point μ i . Therefore, the uncertainty in that relation has to be considered as well by virtue of u i = u(μ i + ε x i ): where J (μ) is the Jacobian matrix of the engine map = ( 1 , . . . , n ) T . The disturbance vector ε has dimension d + n = 2 + 6 = 8. Applying a multivariate second-order Taylor approximation yields where H f is the Hessian matrix of f.
The partial derivatives of f can either be calculated analytically or, if that is not feasible, approximated numerically using adequate difference approximations. The final step is the calculation of the variances and covariances of the disturbance vector ε. Since ε, defined in (10), consists of two distinct components, a case discrimination is necessary: • Case B: j > 2: Define an index i = j−2. This is done in order to ease the notation because of the fact that, for example, the first component of ε u shows up in the third component of ε, and so on-see (10).
Here, A i denotes the ith row, and A i,j the jth entry in the ith row of a matrix A. The last term in (15) has a zero-expectation because ε x and ε u are uncorrelated. For the calculation of E ε j ε k for j = k, a case discrimination is necessary as well: where ¬j indicates the index in {1, 2} which is not j.
Plugging the variances E[ε 2 j ] and covariances E ε j ε k into (12), the approximate expected value can be evaluated. Note that the derived expected value is only one of N summands in the sum in (9). As mentioned before, all these summands can be evaluated equivalently. However, the calculation of the expected value of the penalty term g is different. Using (2) yields Since the distribution of εū was assumed to be normal with a zero-mean and known covariances, Y k := a k (ū + εū) is also normally distributed with known expected value and variance. Therefore, the expected penalty can be calculated as where h k is the known probability density function of Y k . Plugging everything into (12) and (9), the expected value of the cost function can be approximated. Note that, since the sampling time in the drive cycle data is 1 s, the approximated quantity is the expected cost per second. Dividing by the average velocity of the data, the expected cost per kilometre can be obtained.
In order to show the efficiency of the proposed approximation, the accuracy and calculation effort of the suggested method are compared with those of a basic Monte Carlo approach, as used in Wasserburger, Hametner, and Didcock (2020). For that approximation, the sample of drive cycles, which was used to estimate the GMM in Section 3.1, is evaluated on a regular grid over the operating space with as many nodes as there are components in the GMM, using a specific set of engine maps. The weight of each node is determined by the number of data points that are closest to each node. For each drive cycle and each node, a sample of input disturbances is generated and the respective functions are evaluated. The sum over the nodes is calculated, transformed into cost per kilometre and averaged over the cycles. By adjusting the number of evaluated drive cycles and the input disturbance sample size, more or less accurate approximations can be achieved and compared. Additionally, the accuracy is also compared with the naive single cycle approach where all uncertainties are completely ignored, as would be the case in standard, deterministic engine calibration.
As a benchmark for assessing the accuracy of these approximations, the drive cycle costs are calculated without the use of a grid. This means that each data point in the drive cycle data set is evaluated with a large, individual set of input disturbances. Then, as before, the average of the costs per kilometre of the different cycles is determined.
In addition, the number of function evaluations (evaluations of the functions f f and f e ) are compared. These functions can be complex and computationally expensive to evaluate, so the number of function evaluations matters greatly for any optimization algorithm. For the Monte Carlo simulations, the number of evaluations is simply determined by the number of operating points, drive cycles and simulated input disturbances. For the suggested approximation approach, the number of function evaluations is defined by the number of operating points and the approximation formula used for the second derivatives in (12). In this work, a central difference approximation of the derivatives is deployed. The results of the comparison are illustrated in Figure 4.
The proposed approximation performs very well with an error of only 2.6%. While a similar accuracy could, of course, also be achieved using Monte Carlo simulation, the computational load is multiple orders of magnitude larger in that case. For a simulation using the same number of function evaluations (MC 4 ) as the proposed approximation, the error is increased almost fivefold. The naive approach, which disregards any uncertainty, is by far the computationally least costly, but an approximation error of 43.9% is unacceptably large for the naive method to be used in a reasonable, robust engine calibration.

Optimization
The final step is the actual optimization of the engine maps. Using the derivations of the above sections, problem (5) can be solved numerically. Since the smoothness of the engine maps is handled by means of the penalty function (2), only upper and lower bounds for the engine parameters remain as constraints. In this work, the optimization is performed by a modified particle swarm optimization (PSO) that uses a neighbourhood topology where not all particles can communicate directly with each other, so that the risk of getting stuck in a local minimum is mitigated-see e.g. Engelbrecht (2007). Particle swarm optimization has the advantage that it does not depend on a lot of metaparameters and is therefore easy to tune. Furthermore, it is simple to parallelize the evaluations of the objective function for each particle in each iteration and can therefore be sped up considerably.
The optimization goal considered in this work, as defined in (5), is the minimization of the expected value of the cost associated with fuel consumption, emissions and map smoothness. Alternatively, if one is not only interested in average performance but in limiting the spread of the different cost realizations, one can also consider the variance of the cost distribution in addition and minimize a weighted sum of expected value and variance as a combined criterion. Other options are the risk measures considered in Wasserburger, Hametner, and Didcock (2020), like value-at-risk (VaR) or conditional value-at-risk (CVaR).
The functions f f and f e here represent local quadratic model networks, as proposed in Hametner and Jakubek (2013) and Hametner and Jakubek (2007). The local models are chosen to be fully quadratic models and the global model output is a weighted sum of the local model outputs. The modelling is based on the steady-state measurement data of a diesel engine.
Note that the GMM-centres μ i are not necessarily the same as the engine map grid points x i . The decision variables are the engine parameters at the grid points u(x i ). In every iteration step, the engine maps are interpolated in order to calculate u(μ i ) so that the approximation (9) can be evaluated for obtaining the particle fitness in the PSO. The same regular grid is also used for the approximation of the second derivatives of the engine maps for (2).

Results and discussion
In this section, the results of the stochastic optimization based on the efficient objective approximation are presented and discussed. Several different optimization approaches are compared: (a) stochastic optimization using the approximation approach of Section 3; (b) stochastic optimization using Monte Carlo simulation based on the regular grid for cost evaluation; (c) deterministic optimization based on a single cycle without input disturbance evaluated on the regular grid; (d) deterministic optimization based on a single cycle without input disturbance evaluated on the GMM-nodes.
The grid of operating points for the evaluation of (b) and (c) is the same as that used for the engine map interpolation mentioned before in Section 3.4. It consists of N = 15 operating points, which is the same as the number of components of the GMM in order to enable a fair comparison. The weights are chosen as p f = 1, p e = 1 and p s = 10 6 . Note: Each row represents one stochastic or deterministic engine calibration. For each optimization approach, an engine calibration consisting of the six engine maps is obtained. The calibrations can be compared by evaluating the objective function for a large number of disturbance realizations and drive cycles, as in the benchmark evaluation in Section 3.3. Note that the RDE drive cycles used in the evaluation are not the same as those used for the optimization. The results of this evaluation can be found in Table 1 and Figure 5.
The figure shows histograms of the distributions of the cost function based on simulated input disturbances and 1000 RDE-compliant drive cycles. In addition, in Table 1 the sample-based estimates of the expected value, variance, VaR and CVaR are listed for each calibration, the latter two both with a confidence level of 0.9. First of all, it is evident that the results from the two stochastic optimizations are virtually the same. The only difference between the two optimization approaches is the way the objective function is approximated. As demonstrated in Figure 4, the Monte Carlo simulation-based approach (b) involves a much larger number of function evaluations and therefore a greatly increased computation time during the optimization process. The fact that the results are so similar proves the usefulness of the proposed optimization approach, since the computational load is greatly reduced without impairing the quality of the solution.
Furthermore, the calibrations from the stochastic optimizations (a) and (b) are superior to the classic deterministic approaches (c) and (d). For the deterministic formulations, the variance is larger and also the average cost is much higher. Disregarding the various uncertainties and disturbances in the optimization process evidently leads to calibrations that perform poorly under realistic, uncertain and varying conditions. However, in terms of the undisturbed cost, which is the objective value without input disturbances and only for the one drive cycle upon which the deterministic optimization was based, the two calibrations (c) and (d) perform best. This is as expected, because these calibrations were specifically designed for that single cycle and without uncertainties in mind.
This reveals the issue of overfitting a calibration to a specific drive cycle and ignoring related uncertainties. Even though the calibrations (c) and (d) perform better on the single cycle, they are far inferior calibrations overall, which can be seen in all the other performance indicators of Table 1. Lastly, the results of calibration (d) are better than those of (c) in all respects. This shows that, even if the uncertainty is not incorporated in the optimization, the performance can be improved by considering operating points that better reflect the car usage defined by the drive cycles. Using the operating points obtained through the GMM provides a viable alternative for the deterministic optimization approach (c).
Some intersection plots of the resulting engine maps are shown in the left and middle column of Figure 6. For different levels of the mean effective pressure, two of the engine parameters, exhaust gas recirculation (EGR) and swirl position, are depicted over the engine speed range for all four calibrations. Here, 100% EGR means that the EGR valve is completely open. The maps are smooth and quite similar for low loads but vary slightly for higher loads. In the right-hand column, the cost distributions for trajectories along these engine map intersections are illustrated. The trajectories are simply defined by increasing engine speed and constant mean effective pressures spread around the set values. The histograms show that not only fuel consumption and emissions, but also the uncertainty in the results, increase with higher loads. However, the cost variability is much lower for the calibrations (a) and (b). This means that the stochastic optimization takes the engine maps to those areas in the parameter space where the uncertainties do not influence the outcome as much.

Conclusion and outlook
In this work a novel, probabilistic approach for the optimization of technical systems under random disturbances arising during operation was presented. By taking disturbances in the decision variables and the operating cycles directly into account in the objective function, a stochastic optimization formulation was attained. For that purpose, the related disturbances were modelled as random variables. A multitude of drive cycles was considered by estimating the multivariate distribution over the operating space. The stochastic optimization problem was solved by minimizing the expected value of the distribution of the random objective function. An additive disturbance formulation based on the modelled distributions was the basis for an efficient approximation of the expected cost. As an application example for the proposed optimization approach, a steady-state based calibration of a diesel engine was performed. It was shown that the incorporation of drive cycle and input related uncertainties is vital for achieving a reliable and robust engine calibration. Disregarding the variability can lead to over-fitted calibrations that perform poorly under a broader range of operating conditions. In addition, the proposed approximation method allows for greatly reduced computation time without sacrificing the accuracy and quality of the solution when compared to a more time-consuming Monte Carlo sampling approach for the evaluation of the objective distribution.
Future research may be focused on the consideration of transient effects and the development of a dynamic optimization algorithm dealing with the discussed stochastic influence. The input related uncertainties might have to be modelled as stochastic processes rather than uncorrelated normal distributions. Evaluating and dealing with the effect of random disturbances on dynamic models might pose an interesting challenge. Furthermore, a dynamic approach would raise the question how to identify the most representative transients through the operating space for a dynamic optimization.
In that context, accurate modelling of the nonlinear, dynamic engine behaviour is also paramount. It is a well-known fact that the modelling of engines and powertrains is a very difficult task, especially when dynamic operation is considered. Since model-based calibration relies heavily on the quality of the models utilized, future research must also take this into account. Knowing and handling the limitations of model accuracy, be it static or dynamic models, directly in the optimization might therefore be a promising research direction as well. If model uncertainty can be quantified depending on the parameter and operating space, model uncertainty could be included in the stochastic optimization, for example, through the consideration of model confidence intervals.