Inverse Parametric Uncertain Identification using Polynomial Chaos and high-order Moment Matching benchmarked on a Wet Friction Clutch

TL;DR

This paper introduces an efficient inverse method combining Polynomial Chaos and moment matching for identifying probabilistic models of system parameters, demonstrated on a wet clutch system, improving accuracy and reducing computational effort.

Contribution

It presents a novel inverse uncertainty identification approach using high-order moments and polynomial chaos, significantly reducing simulation costs and increasing model accuracy.

Findings

Achieved a tenfold reduction in required simulations.

Improved log-likelihood by approximately 4%.

Enhanced accuracy of output probability density estimation by up to 47%.

Abstract

A numerically efficient inverse method for parametric model uncertainty identification using maximum likelihood estimation is presented. The goal is to identify a probability model for a fixed number of model parameters based on a set of experiments. To perform maximum likelihood estimation, the output probability density function is required. Forward propagation of input uncertainty is established combining Polynomial Chaos and moment matching. High-order moments of the output distribution are estimated using the generalized Polynomial Chaos framework. Next, a maximum entropy parametric distribution is matched with the estimated moments. This method is numerically very attractive due to reduced forward sampling and deterministic nature of the propagation strategy. The methodology is applied on a wet clutch system for which certain model variables are considered as stochastic. The number of required model simulations to achieve the same accuracy as the brute force methodologies is decreased by one order of magnitude. The probability model identified with the high order estimates resulted into a true log-likelihood increase of about 4% since the accuracy of the estimated output probability density function could be improved up to 47%.

Tables

Table 1. Table 1 : Wiener-Askey polynomial chaos.

distribution, $f_X$	polynomials, ${\varphi }_i$	support, $𝒳$
Gaussian	Hermite	$[-\infty ,\infty ]$
Gamma	Laguerre	$[0,\infty ]$
Beta	Jacobi	$[-1,1]$
Uniform	Legendre	$[-1,1]$

Table 2. Table 2 : Summary of log-Likelihood estimation.

	$\log L$	PDF est.	MAE	$S_{gPC}$	$S_{MC}$	$S_{qMC}$
${\alpha }_G^{\ast }$	176.7	Gauss	0.183	25	458	577
		ME	0.155	25	629	1171
${\alpha }_{ME}^{\ast }$	183.7	Gauss	0.117	25	265	370
		ME	0.088	25	383	1009

Equations

$${\alpha }^{\ast }$$

$$=\underset{\alpha }{\max }\sum _{l\in \mathcal{N}}\log f_{Y_l}\left(y_l|\alpha \right)$$

$$y(\mathrm{x})={\sum }_{i=1}^{\infty }{c}_i{\psi }_i(\mathrm{x})$$

$$y^{(d)}(\mathrm{x})={\sum }_{i\in \mathcal{I}}{c}_i{\psi }_i(\mathrm{x})$$

$$⟨{\psi }_i,{\psi }_j⟩\equiv \mathrm{E}\left\{{\psi }_i{\psi }_j\right\}={\int }_{\mathcal{X}}{\psi }_i(\mathrm{x}){\psi }_j(\mathrm{x})f_{\mathrm{X}}(\mathrm{x})\text{d}\mathrm{x}$$

$$⟨{\psi }_i,{\psi }_j⟩$$

$$=\prod _{k=1}^n⟨{\varphi }_{i_k}^{(k)},{\varphi }_{j_k}^{(k)}⟩={\delta }_{ij}$$

$${\mu }_m=\mathrm{E}\left\{Y^m\right\}=\mathrm{E}\left\{y(\mathrm{X}{)}^m\right\}={\int }_{\mathcal{X}}y(\mathrm{x}{)}^m{f}_{\mathrm{X}}(\mathrm{x})\text{d}\mathrm{x}$$

$${\mu }_m^{(d)}$$

$$=\sum _{i_1\in \mathcal{I}}\cdots \sum _{i_m\in \mathcal{I}}{c}_{i_1}\cdots c_{i_m}⟨{\psi }_{i_1}\cdots {\psi }_{i_m}⟩$$

$${\mu }_1^{(d)}$$

$${\mu }_2^{(d)}$$

$$\mathrm{X}=\mathrm{x}(\mathrm{\Theta }|\alpha )\to Y=y\left(\mathrm{x}\left(\mathrm{\Theta }|\alpha \right)\right)\equiv \eta \left(\mathrm{\Theta }|\alpha \right)$$

$$c_i(l|\alpha )=⟨\eta (l,\mathrm{\Theta }|\alpha ),{\psi }_i⟩={\int }_{\vartheta }\eta (l,\theta |\alpha ){\psi }_i(\theta )f_{\Theta }(\theta )\text{d}\theta$$

$$c_i(l|\alpha )\approx \sum _{j\in \mathcal{Q}}\eta \left(l,{\theta }_j|\alpha \right){\psi }_i\left({\theta }_j\right)w_j$$

$$f_Y(y|\alpha )=\phi \left(\frac{y-{\mu }^{(d)}(\alpha )}{{\sigma }^{(d)}(\alpha )}\right)$$

$$f_Y(y|\lambda )=\exp \left(-{\sum }_{m=0}^M{\lambda }_m{y}^m\right)$$

$$\underset{\lambda }{\max }\log {\int }_{\mathcal{X}}\exp \left(-{\sum }_{m=1}^M{\lambda }_m{x}^m\right)\text{d}x+{\sum }_{k=1}^M{\lambda }_m{\mu }_m^{(d)}(\alpha )$$

$${\mu }_m^{(d)}=\sum _{\text{uline}i\in {\mathcal{I}}_{mp}}\left(\genfrac{}{}{0ex}{}{m}{i_1\cdots i_p}\right)\cdot \prod _{j=1}^n⟨\prod _{k=1}^p{{\varphi }_{\text{uline}k(j)}^{(j)}}^{i_k}⟩\cdot \prod _{k=1}^p{c}_k^{i_k}$$

$${\mu }_m^{(d)}=\sum _{\text{uline}i\in {\mathcal{I}}_{m,p}^{\ast }}{a}_{\text{uline}i}^{n,d,m}\prod _{k=1}^p{c}_k^{i_k}$$

$$y(x)$$

$$T_1({\omega }_m,{\omega }_1)={\omega }_m^2\frac{f_t(\frac{{\omega }_1}{{\omega }_m})}{f_c(\frac{{\omega }_1}{{\omega }_m}{)}^2}.$$

$$T_b({\omega }_2)=T_{b0}+b{\omega }_o$$

$$\dot{\mathrm{s}}$$

$$\mathrm{A}$$

$$p=\max (0,\frac{z-z_p}{z_M-z_p})p_{hc}$$

$$T_c({\omega }_1,{\omega }_2,p,z)={\delta }_t(z)\cdot \alpha p+(1-{\delta }_t(z))\cdot \frac{\gamma }{z_M-z}\cdot ({\omega }_1-R{\omega }_2)$$

$$J_1{\dot{\omega }}_1$$

$$J_2{\dot{\omega }}_2$$

$$(J_1+\frac{J_2}{R^2}){\dot{\omega }}_1=T_1({\omega }_m,{\omega }_1)-\frac{1}{R}{T}_b\left(\frac{{\omega }_1}{R}\right)$$

$$u_i$$

$${\tilde{\theta }}_i$$

$$\mathrm{x}$$

$$\mathcal{L}[p,\lambda ]=\mathcal{W}[p]+{\sum }_{k=0}^K{\int }_{\mathcal{X}}{\lambda }_k(x^k-{\mu }_k)p(x)\text{d}x$$

Full text

††journal: Mechatronics

1 Introduction

The increasing performance demands in design and control of mechatronic applications lead to an upward trend in the need for accurate and robust models [1]. These mathematical relations allow to make predictions on the behavior and performance of the system in a virtual computational environment. It enables efficient design processes, leading to faster successive optimization iterations, so that more reliable products can be made at lower production cost.

In this research, we consider the modeling of a wet clutch system. Wet friction clutches are hydraulic-mechanical devices used to transmit torque from an input shaft to an output shaft by means of friction. These systems are used in various types of automatic transmissions to selectively engage gear elements. During engagement, multiple plates make contact which enables power transfer from motor to load. Contrary to dry clutches, the plates are bathed in oil to assure better heat conduction of the friction losses. Wet friction clutches thus result in higher torques which make them suitable and primordial in off-road vehicles and agricultural machines [2]. Models of clutch systems have been widely used for fault diagnosis purposes [3] and numerical optimization of both design and control [4]. Clutch models identified for the wet clutch considered in this paper have been previously implemented in condition monitoring [5] and control [6]. Additionally, model based control techniques applied on the wet clutch setup, have proven shorter convergence time and can avoid excessive control inputs that lead to unsafe operations [7].

Although the value of having accurate and robust models for mechatronic applications, as illustrated with the wet friction clutch, is straightforward, the construction of accurate models remains a cumbersome process. Since these systems are plagued by their intrinsic complexity and nonlinear behavior [2], variations in the physical phenomena over time due to e.g. wear of plates or degradation and centrifugal effects of the oil in the wet friction clutches; model discrepancies are inevitable [8]. Due to these uncertaintes, deterministic models become less useful. Alternatively, one can consider a system model of stochastic nature that predicts an output distribution for given control input. Models with ingrained parametric uncertainty have proven increased robustness in the field of fault diagnosis [9] and control [10, 11]. In this research we adopt the available system model and focus on identifying a probabilistic model to several lumped model parameters. Mostly (mechatronic) models are derived from first principles. Lumped parameter values are then determined empirically by fitting the model to the experiments. Considering that this strategy generates a model that can explain the experiment set, it is no direct proof of the validity of the physical model. It may occur that the parameter value fitted on this experiment set compensates for an inherent model shortcoming. Rather than increasing the model complexity, we propose to couple a parameterized input probability distribution to some of the lumped model parameters [12]. We assume that the uncertainty structure (i.e. the parameterized input probability distribution) is fixed, and only its parameters, such as mean and standard deviation (e.g. for a normal distribution), need to be identified based on what parameters best explain a number of experiments.

Stochastic system identification can be performed by means of prevailing methods as (quasi) Monte Carlo (MC) sampling techniques for assessing the propagation of the uncertainty parameter to the output variables (response) of the system model [13, 14]. These techniques are however curtailed by the computational inconvenience that comes with the numerous forward simulations required to achieve an acceptable degree of accuracy. Furthermore, in an inverse uncertainty identification setting, the uncertain parameters need to be optimized with respect to the correspondence of model responses to measurements, leading to even higher computational costs. More recent work has revived the exploitation of generalized Polynomial Chaos (gPC) expansions that can lead to significant gains in terms of computational feasibility [15, 16, 17, 18, 19]. In the gPC framework, the propagation of input uncertainty to output variables is realized by developing the stochastic subspace through a polynomial series expansion. We propose in this paper to choose the polynomial basis so that it is orthogonal with respect to the joint probability density function of the random input variables. Hence, an efficient mathematical context emerges that allows to express the statistical moments as a function of the polynomial coefficients. Once these moments are available, the conditional Probability Density Function (PDF) of the outcome of the experiment can be retrieved through the well-known method of moments [20, 21].

To our knowledge, application of the gPC framework to estimate probabilistic moments of a stochastic output model has remained limited to mean and variance estimates. In this work we devised an efficient algorithm to calculate high-order moments as well, so to increase the information content that can be passed to the moment matching algorithm. The major computational bottleneck is thereby isolated and can be executed as an offline step. Figure 1 illustrates the incorporation of the identified output distributions for each experiment within a Maximum Likelihood Estimation (MLE) framework [22]. This methodology enables to find sufficient parametric uncertainty that a model needs to include so that different experiments can be explained [23]. A Genetic Algorithm (GA) is used here to optimize the parameterized distributions of the model parameters, because of its exploratory characteristics [24]. First, we elaborate our inverse uncertainty identification method to define different methods and concepts that are included within our theoretical framework. Next we address the wet clutch system and discuss in detail our practical set-up. Lastly, the developed methods are applied on the clutch application to identify the input probability model characterized by high-order moments.

2 Inverse uncertainty identification method

2.1 Introduction

We consider physics based nonlinear forward models, $\mathcal{Y}:\mathcal{N}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$, used to generate predictions, $Y_{l}=\mathcal{Y}(l;\bm{\mathrm{x}})$, for the univariate outcome of experiments, $y_{l}$, that are characterized by some index $l\in\mathcal{N}$. The index $l$ is assumed to contain sufficient information to render the experiments deterministic from the perspective of the model, i.e. the information contains the operational settings that can be adjusted by an experimenter to the mechatronic system. The variable $\bm{\mathrm{x}}\in\mathbb{R}^{n}$ contains any (lumped) circumstantial parameters that, once identified, are assumed to remain constant for all experiments. We consider a number ($N$) of experiments within the index set $\mathcal{N}=\{1,\cdots,N\}$. When we assume that the measurements are deterministic, the mismatch stems either from a lack of circumstantial knowledge or suggests a genuine shortcoming of the model. In the former case, a latent variable may be present while in the latter case, the model does not, or not correctly, take into account all underlying physical phenomena.

In recent literature issues are addressed by associating a random variable, $\bm{\mathrm{X}}\in\mathcal{X}\sim\mathcal{P}_{X}$, to the model parameter set $\bm{\mathrm{x}}$ [12]. This approach renders the forward model stochastic and rather than exact predictions the model generates a distribution of possible outcomes. Although this strategy may deny an underlying physical reality, it grants the model practical use as it hands the user a measure for how far off a prediction may be. With respect to this concept, our contribution is twofold.

First, we propose an efficient numerical method to obtain an approximation of the nonlinear transformed output distribution of the random variable $Y_{l}=\mathcal{Y}(l;\bm{\mathrm{X}})$, by moment matching with a parametric distribution. The more moments are taken into account, the better the approximation will be. To facilitate cheap yet reliable moment estimates we engage the gPC framework and extend it so that high-order moments can be extracted rigorously. Such benefits both the MLE, as nonlinear effects can be taken into account, as well as the post identification usage of the stochastic model.

Secondly, we propose a fast computational method to associate a probability model to the variable, $\bm{\mathrm{X}}$, based on a number of experimental observations, $\{y_{l}\}_{l\in\mathcal{N}}$. To that end we put forward a parametric model, $\mathcal{P}_{X}(\bm{\mathrm{\alpha}})$, for the probability of the random variable, $\bm{\mathrm{X}}$, parameterized by the variable $\bm{\mathrm{\alpha}}$, and we identify an optimal value for $\bm{\mathrm{\alpha}}$ by means of Maximum Likelihood Estimation (MLE). Hereinafter, for notational convenience we make no distinction between model $\mathcal{P}_{X}(\bm{\mathrm{\alpha}})$ and parameter $\bm{\mathrm{\alpha}}$, and refer to it as the (input) probability model.

2.2 Maximum Log-Likelihood Estimate

The proposed identification procedure is formulated as an inverse problem (1). That is, we want to find the probability model, $\bm{\mathrm{\alpha}}^{*}$, that maximizes the natural $\log$ of the likelihood, $L(\bm{\mathrm{\alpha}})$. The likelihood of a given probability model, $\bm{\mathrm{\alpha}}$, is calculated as the product (or summation of the natural logarithm) of the probability density functions, $f_{Y_{l}}(y_{l}|\bm{\mathrm{\alpha}})$ $l\in\mathcal{N}$. Probability $f_{Y_{l}}(y_{l}|\bm{\mathrm{\alpha}})$ is the conditional relative likelihood that the value of the stochastic model output, $Y_{l}$, of experiments $l$, as obtained through simulation given probability model $\bm{\mathrm{\alpha}}$, would equal the experimentally observed value, $y_{l}$. The higher the value of $L$, the more likely the forward model.

[TABLE]

With regard to the construction of the MLE cost function, $\log L$, it is clear that the computational bottleneck is predominated by the numerical evaluation of the conditional probabilities, $f_{Y_{l}}(y_{l}|\bm{\mathrm{\alpha}}),l\in\mathcal{N}$. Hence, we will require an efficient numerical procedure to quantify the probabilities, $f_{Y_{l}}$.

A novel procedure of such kind is described here next.

2.3 Moment Estimation

To evaluate the conditional output probability density functions, $f_{Y_{l}}\left(y_{l}|\bm{\mathrm{\alpha}}\right)$ we will engage the generalized Polynomial Chaos (gPC) expansion framework. The gPC framework accounts for a number of advantageous mathematical conditions to propagate uncertainty. Such are accomplished by approximating the forward nonlinear model with a polynomial expansion. By subsequently choosing the polynomial basis so that it satisfies orthogonality conditions with respect to the probability density function of the random input variables, an efficient mathematical context emerges that allows to express the statistical moments in function of the coefficients associated to the expansion. Conventional application is limited to the mean and variance. In this work we propose a novel method to incorporate high-order moments as well.

2.3.1 Generalized polynomial chaos

According to the polynomial approximation theorem, any smooth function, $y:\mathbb{R}^{n}\rightarrow\mathbb{R}$, is equivalent to an infinite polynomial series [25]

[TABLE]

From a computational perspective our interest is however reserved to the $d$-th order approximation. That is, let $\mathcal{P}^{d}_{n}$ be the $n$-variate polynomial space of at most degree $d$ and let $\bm{\mathrm{\psi}}=\{\psi_{i}\}_{i=1}^{p}$ serve as basis for $\mathcal{P}^{d}_{n}$ with $p=\frac{(n+d)!}{n!d!}$. The $d$-th order approximation is then given by (2) for given polynomial coefficients $\{c_{i}\}_{i=1}^{p}$. An $n$-variate basis $\bm{\mathrm{\psi}}$ can be constructed from $n$ univariate bases $\bm{\mathrm{\phi}}^{(k)}=\{\phi^{(k)}_{j}\}_{j=0}^{d},k\in\{1,\dots,n\}$. Consider the basis vector elements, $\psi_{|\uline{i}|\leq d}=\prod_{k=1}^{n}\phi^{(k)}_{\uline{i}(k)}$, with multi-index $\uline{i}=(i_{1},\dots,i_{n})$ and where $|\uline{i}|$ is defined as $\sum_{k=1}^{n}i_{k}$. For notational convenience we exploit the bijection between index $\uline{i}$ and $i$ taking values in $\mathcal{I}=\{1,\dots,p\}$.

[TABLE]

Within the context of uncertainty propagation this representation allows to establish advantageous computational conditions by a distinct choice of basis, $\bm{\mathrm{\psi}}$. Assume that $\bm{\mathrm{x}}$ is composed of $n$ independently distributed random variables, $X_{k}$, with known supports and PDF, $f_{X_{k}}:\mathcal{X}_{k}\subseteq\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$. The joint support, $\mathcal{X}$, and PDF, $f_{\bm{\mathrm{X}}}$, are given by $\bigotimes_{k}\mathcal{X}_{k}$ and $\prod_{k}f_{X_{k}}$, respectively. Now recall that our goal is to propagate the input uncertainty on $\bm{\mathrm{x}}$ to output $y$. By choosing the univariate bases, $\bm{\mathrm{\phi}}^{(k)}$, so that they satisfy an orthogonality condition w.r.t. the PDFs, $f_{X_{k}}$, associated to the respective variables, $X_{k}$, the statistical moments can be calculated in function of the polynomial coefficients.

Orthogonality of a basis is established in function of an inner product definition. We take interest in the inner product defined in (3). A polynomial basis is orthogonal w.r.t. $f_{X}$ if it satisfies the orthogonality condition $\langle\psi_{i},\psi_{j}\rangle=\delta_{ij}\langle\psi_{i}^{2}\rangle$.

[TABLE]

Note that if we construct the multivariate basis as described above and so that the generating univariate bases, $\bm{\mathrm{\phi}}^{(k)}$, satisfy the orthogonality condition w.r.t. to the PDFs, $f_{X_{k}}$, also the multivariate basis, $\bm{\mathrm{\psi}}$, will satisfy the orthogonality condition w.r.t. the joint PDF, $f_{\bm{\mathrm{X}}}$, considering that

[TABLE]

2.3.2 Stochastic relation

Now recall that we desire to characterize the stochastic properties of the output by quantifying its stochastic moments, $\mu_{m}$

[TABLE]

When we substitute the $d$-th order polynomial approximation of the output model in (4), we retrieve an approximate expression for the $m$-th moment in function of the polynomial expansion coefficients.

[TABLE]

From this expression one may easily verify that the first two moments can be estimated as shown below. Efficient estimation of the high-order moments is discussed in 2.5.

[TABLE]

We emphasize that this result only holds when the series is expanded over a basis that satisfies the orthogonality condition w.r.t. $f_{\bm{\mathrm{X}}}$. Several standard probability distributions are associated to known polynomial families by the Wiener-Askey scheme [15]. An overview is presented in Table 1. If the input stochasticity does not correspond with a standard distribution, a variable transformation can be used.

2.3.3 Variable transformation

A set of models with parameterized uncertainty parameters $\bm{\mathrm{X}}$ can now be characterized. To that end we introduce the parametric transformation, $\bm{\mathrm{x}}(\cdot|\bm{\mathrm{\alpha}})$, that maps standard random variable, $\bm{\mathrm{\Theta}}\in\vartheta$, with entries that are distributed according to one of the distributions in Table 1, to the random variable, $\bm{\mathrm{X}}$. The multivariate expansion basis, $\bm{\mathrm{\psi}}$, can then be generated from the corresponding standard polynomials. We emphasize that the forward model becomes a function of the variable, $\bm{\mathrm{\Theta}}$. Consequently, the conditional PDF of the outcome is now fully determined by parameter $\bm{\mathrm{\alpha}}$, and the distribution of $\bm{\mathrm{\Theta}}$.

[TABLE]

2.3.4 Coefficient determination

The polynomial coefficients from (2) can be quantified numerically [18] by projection of the forward model on the polynomial space exploiting the properties of the inner product. All terms but one will vanish resulting in

[TABLE]

for all $i\in\mathcal{I}$. Since the scope of this work is on general nonlinear forward models, we can not evaluate the associated integral explicitly. We therefore approximate the integrals using Gaussian-quadrature (8). For details we refer to [17].

An univariate Gaussian-quadrature of order $q$ is defined as a signature of the same size. A signature is defined as a set $\{(w_{j},\bm{\mathrm{\theta}}_{j})\}_{j\in\mathcal{Q}},\mathcal{Q}=\{1,\dots,q\}$, where $\bm{\mathrm{\theta}}_{j}$ and $w_{j}$ are associated to the position and the weight attributed to each element, respectively. The signature will depend on the weighting function, $f_{\Theta}$, and is as such related to the polynomial basis. A quadrature of order $q$ is exact for polynomials up to degree $2q-1$. The order hence affects the number of coefficients that can be retrieved correctly and therefore the accuracy of the polynomial approximation. In the multivariate case, a full tensor product of univariate quadrature rules can be considered. We note that the number of collocation points will thus scale exponentially with the number of dimensions.

[TABLE]

Note that, similar to MC techniques, the gPC approach requires a number of evaluations of the forward model, $\eta(l,\bm{\mathrm{\theta}}_{j}|\bm{\mathrm{\alpha}})$. The difference is that the estimation of the statistical moments is realized through a mathematical detour. The focus of approximation in the gPC framework is on modeling the polynomial response function whilst that of the MC approach is on the direct estimation of the output distribution. The accuracy of gPC depends on the capacity of the basis to capture the nonlinearity of the forward model rather than on the capacity of the sampling method to properly represent the input uncertainty by the spatial distribution of the sample points. It is the prevailing consensus that an equivalent level of accuracy can be achieved with only a fraction of the input points of any MC approach.

2.4 PDF fitting with the method of moments

Once the approximate stochastic moments, $\mu_{l,m}^{(d)}(\bm{\mathrm{\alpha}})$, of the $l$-the experiment for given input stochastic model, $\bm{\mathrm{\alpha}}$, are obtained a model for the conditional PDF of the random outcome, $Y_{l}$, can be extracted by matching these moments to a known parametric PDF, the so called methods of moments. The complexity of the parametric PDF should match the stochasticity of the random input, the expected symmetry of the stochasticity of the outcome and the nonlinearity of the forward model. Two possibilities are discussed for extracting the conditional PDF.

Gaussian distribution: The most straightforward and prevailing approach is to fit a normal distribution to the output PDF based on estimates of the mean and variance. The conditional PDF is then approximated by

[TABLE]

• 2.

Maximum entropy distribution: To cope with an increased number of moments, one can fit a maximum entropy distribution. Here, we seek a distribution $f_{Y}$ that maximizes the entropy $\mathcal{W}=-\int_{\mathcal{Y}}f_{Y}(y)\log f_{Y}(y)\text{d}y$ subject to $M+1$ moment matching constraints. The solution to this optimization problem is given by the following parameterized PDF (see Appendix A for derivation)

[TABLE]

Substituting the latter into the original constrained optimization problem, yields the following dual unconstrained optimization problem that can be solved to find values for the parameters $\lambda_{m}$

[TABLE]

A solution exists for any $M$. However, in order to capture fourth order asymmetric effects (i.e. bimodal densities), the number of moments should $\geq 4$ so that the polynomial in the exponent can have two local maxima.

Worth mentioning is the common conception in the literature [26, 27] to obtain the moments, or even directly the PDF, by applying MC techniques on the polynomial expansion (2) which is relatively cheap to evaluate. In this approach the expansion is used as a response function and the benefit of gPC is only exploited through the supposed superior convergence rate of the expansion for polynomials corresponding the input distribution (2). When compared to its peers, gPC is usually applied in this sense [17, 26, 28].

2.5 Estimation of high-order moments

The expression presented in (5) becomes inefficient for $m>2$. Therefore we expand the power of the sum by applying the multinomial theorem to obtain an expression which is numerically far more efficient. Remark that there is only one summation and each term in the series is rendered unique. Moreover, we decompose the high-order multivariate inner products into a product of univariate inner products. As a result we can compute the high-order inner products numerically using a univariate instead of a multivariate quadrature. The multi-index set $\mathcal{I}_{m,p}$ is defined as $\{\uline{i}\in\mathcal{I}_{m,p}:|\uline{i}|=m\}$.

[TABLE]

Details on deriving equation (9) can be found in Appendix C. The algorithmic execution time still demands an upper-limit for $m$. That is because $\dim(\mathcal{I}_{m,p})=\binom{m+p-1}{m}$ where $p$ itself is a combinatorial number. Note that the orthogonality property extends to the high-order inner products, as they occur in (5) or (9) for $m>2$, however in a less predictable fashion.

Hence in practice we only require the multinomial powers $\uline{i}$ whose corresponding inner product is nonzero. Such are stored in the set $\mathcal{I}^{*}_{m,p}$. The corresponding multinomial coefficient and inner product are compressed into coefficients $a_{\uline{i}}^{n,d,m}$ and stored for later use. One may verify that a significant fraction of $\mathcal{I}_{m,p}$ result into zero valued inner products, see Fig. 2.

[TABLE]

To compute the set $\mathcal{I}_{m,p}^{*}$, we devised an algorithm that calculates the next multi-index, satisfying $|\uline{i}|=m$, from the previous multi-index. Because our algorithm generates the admissible multi-indices one by one, it admits further increase of the computational efficiency. The process of generating the multi-indices and the corresponding inner product can be executed in parallel avoiding storage of the entire set $\mathcal{I}_{m,p}$.

The general idea is presented in Fig. 2. The multinomial power problem can be addressed as reorganizing a stack of $m$ items over $p$ possible locations. We use two logical operators, push ($\mathcal{P}$) and fork ($\mathcal{F}$). A row is ‘pushed’ if the last location is nonzero. Otherwise the trailing nonempty location is ‘forked’ which involves pushing $1$ item from this location to the next and placing the remaining items back to the last location. The algorithm is initiated with all items in the last position. The calculation time of the parallel offline step is shown in Fig. 3 and could be limited, e.g. $(n,d,m)=(4,3,5)$, to $18$\mathrm{s} nonetheless that $575,757$ inner products are validated. The generation of the multi-indices itself took only about $1.2$\mathrm{s}. Calculations were performed on a 2.1GHz, 4 GB RAM laptop.

Above procedure allowed us to compute up to the 5th moment quite efficiently. As discussed in section 2.4, this allows us to fit bimodal output distributions rigorously.

2.6 Forward uncertainty propagation

In order to illustrate the forward propagation of uncertainty concatenating the techniques described in the sections 2.3 to 2.5, we apply the methodology on an illustrative univariate nonlinear forward model, $y(x)$.

[TABLE]

For convenience we further assume that the random variable, $X$, is distributed according to the standard normal distribution, $\mathcal{N}(0,1)$. The forward propagation of this random variable passed through the function $y(x)$, results into a camelback output distribution as is depicted in Fig. 4. The lower right plot shows the density of the stochastic input, $X$. The density of the output random variable, $Y$, is plotted in the upper left graph. Additionally, an MC sample set of $N=50$ is visualized by the black dots with corresponding $10$ bin histogram in gray.

Figure 5(a) depicts the polynomial coefficients for varying quadrature order, $q$, and using Hermite chaos. One may observe that the higher order coefficients are poorly approximated with low order quadratures. Recall that a quadrature of order $q$ is exact for polynomials up to degree $2q-1$. Now, assuming that model $y(x)$ contains polynomials up to degree $d$ and that we desire to approximate the $d$-th coefficient, then the integrand in (7) will contain a polynomial of degree $2d$. In order to be exact, the quadrature should therefore have order $d+1$ at least. Remark that the coefficient estimate becomes increasingly precise for $q$ surpassing $d+1$, as can be seen for e.g. $i=3$. That is because, $y(x)$ is truly a polynomial of degree $\gg d$.

In Fig. 5(b) we compare the stochastic moments numerically obtained from the true PDF111For a monotonic function, $y(x)$, with inverse, $x(y)$, this is equal to, $f_{Y}(y)=\lvert x^{\prime}(y)\rvert\cdot f_{X}(x(y))=\lvert y^{\prime}(x(y))\rvert^{-1}\cdot f_{X}(x(y))$. with those obtained from gPC approximation (2) for varying $d$ and $q$. Note that we have chosen the quadrature order so that $q\geq d+1$. As a result, the poorly approximated coefficients from Fig. 5(a) are never considered in the moment computation. This explains why even for low polynomial order, the higher order moment approximations remain very precise.

In conclusion, the moments are employed to fit parametric distributions to the output distribution, $f_{Y}$. Results are depicted in Fig. 6 when using the Gaussian and maximum entropy distribution. The numerical correspondence between two PDFs is quantified by the Earth Mover’s Distance (EMD) (Appendix B). Each fit is compared to the actual PDF. The PDF fit obtained with the exact moments serves as an upper bound for the amount of information that is contained within the moments.

3 Wet clutch shifting time

3.1 System description

The proposed methodology is applied to identify the parametric uncertainty of a wet clutch system. A clutch is a mechanical device that connects a motoring unit to a load that transforms power into useful work. The objective is to ensure smooth and on demand coupling of the driving shaft to the driven shaft. A wet clutch system realizes this by means of contact friction between two series of friction plates, see Fig. 7 for a schematic cross section [6]. By construction these friction plates are fully submerged in oil so that slip induced heat can be transferred adequately. This avoids thermal overheat and increases the clutch’s life time. A wet clutch uses the same hydraulic system to activate the engagement process. By pressurizing the hydraulic chamber, a freewheeling piston is pressed against the series of interlocking friction plates realizing a power transfer. The pressure in the hydraulic chamber is controlled via a current signal that activates a servo-valve linking the hydraulic chamber to a pressurized reservoir. The pressure difference over the piston determines its position and therewith the clutch engagement.

Since the wet-clutch set-up should be low cost, additional sensors are avoided and a feedforward control strategy using the parameterized current profile shown in Fig. 8 is applied. The signal is characterized by three parameters. An initial maximum current pulse, parameterized by its duration $\Delta t$, gets the piston as close as possible to the friction plates without making contact. During this phase, the hydraulic chamber fills up and initial torque transfer commences due to Couette flow phenomena. The pressure drop invoked by reducing the current to $u_{0}$ avoids a hard collision. Once the plates are in contact a second pulse, defined by $\Delta u$, is generated to have a proper pressure build-up for synchronizing the load. At synchronous speed the torque transfer caused by slip speeds drops away. Therefore, hazardous jerks should be avoided by having slow synchronization. This explains the downward flank of the current pulse. After syncing the pressure is increased to guarantee the plates cannot release again by accident and start slipping. The shifting time, defined as the time interval between initialization of the feedforward control signal and the moment that both shafts obtain the same angular speed, is therefore considered as key characteristic of the engagement process.

3.2 Shifting time model

The shifting time is derived from a dynamic simulation of the clutch engagement. Prior research was on the building of a dynamic model of the engagement process based on first principle physics. Model parameter values were then identified via experimental system identification (from isolated subsystem testing….). The exact model identification is not in the scope of this paper and a general overview is given. We refer to [29, 30, 5, 2] for further reading. Figure 9 illustrates the simplified dynamical scheme of the wet clutch system.

The motor is controlled towards a constant speed $\omega_{m}$ delivering a torque $T_{m}$ to the system. The motor is connected with a torque converter that drives the input shaft of the clutch with torque $T_{1}$ at rotational velocity $\omega_{1}$. The system can be modeled using a nonlinear torque ratio function $\frac{T_{1}}{T_{m}}=f_{t}(\frac{\omega_{1}}{\omega_{m}})$ and a capacity factor function $\frac{\omega_{m}}{\sqrt{T_{m}}}=f_{c}(\frac{\omega_{1}}{\omega_{m}})$ provided by the manufacturer. Consequently, the driving torque $T_{1}$ of the input shaft can be written as

[TABLE]

The clutch is connected to a ratio selector with gear ratio $R=\frac{\omega_{in}}{\omega_{out}}$. Thereafter, the output shaft operates at $\omega_{2}$ as result of the driving torque $T_{2}$. Furthermore, the rotational movement of the output shaft and flywheel, characterized by inertia $J_{2}$, are counteracted by an additional brake torque $T_{b}$.

[TABLE]

The system is fully engaged when both shafts have an equivalent speed $\omega_{1}=R\omega_{2}$. The dynamics of this engagement process can be described by a succession of two distinct phases.

3.2.1 Asynchronous phase

The first phase is characterized by a speed difference $\omega_{1}>R\omega_{2}$ of the clutch shafts. During this phase the hydraulic piston chamber starts filling up. The resulting overpressure $p_{hc}$ is measured and depicted in Figure 10(a). The dynamics of the pressure build-up can be captured by defining a state $\bm{\mathrm{s}}=\begin{bmatrix}p_{hc}&\dot{p}_{hc}&\tilde{z}\end{bmatrix}^{T}$. The state variable $\tilde{z}$ can easily be scaled and truncated within a feasible area to obtain the piston position $z\in[0,z_{M}]$, for which $z_{M}$ is the position that assures full contact between the friction plates. The dynamics are modeled by an affine system with parameters identified to optimally match the pressure $p_{hc}$. Mark that the control input $u$ is the feedforward current signal of the valve as was illustrated in Fig. 8.

[TABLE]

The detailed pressure dynamics are out of scope and we refer to prior research for more information [2, 29]. During the filling phase the plates do not make contact. However, the oil between the friction plates will behave as a planar Couette flow. The resulting shear stresses within the fluid will cause an initial torque transmission and acceleration of the load speed $\omega_{2}$. This phenomena is characterized by the constant $\gamma$, capturing the geometry of the plates and the fluid viscosity. A second influence on the torque transmission is the pressure on the plates $p$. This pressure is a fraction of the total hydraulic overpressure $p_{hc}$ as is illustrated in Fig. 10(a). There is assumed that the pressure build-up on the plates initiates when the piston reaches a constant threshold $z_{p}$.

[TABLE]

The torque transfer due to plate pressure is modeled with parameter $\alpha$, based on the geometry of the plates and fluid characteristics. A naive approach to model the transferred clutch torque $T_{c}$ would be to switch between Couette flow and plate friction once $z$ equals $z_{p}$ in a discrete way. However, the sudden torque shift would cause unrealistic behavior. Therefore, a transient function $\delta_{t}\in[0,1]$ is defined (Fig. 10(b)) to have a smooth transition between the aforementioned sources of torque transfer.

[TABLE]

The expression for the counteracting clutch torque $T_{c}$ on the input shaft holds when $\omega_{1}>R\omega_{2}$. The driving torque of the output shaft equals $T_{2}=RT_{c}$ and is illustrated in Fig. 11(a). Hence, one can formulate the nonlinear system dynamics during the asynchronous phase by considering the torque equilibrium of both shafts.

[TABLE]

.

3.2.2 Synchronous phase

The system shifting process arrives in the synchronous phase once an equivalent speed for both shafts is obtained ($\omega_{1}=R\omega_{2}$). Henceforth, the shafts can be considered as fully engaged and the system becomes independent of the control signal $u$. The driving torque of the output shaft equals $T_{2}=RT_{1}$ since perfect torque transmission is assumed (Fig. 11(a)).

[TABLE]

3.2.3 Shifting time

The purpose of this dynamical model is to determine the shifting time, which is key for various applications of the wet clutch. This is determined by measuring the time interval between initialization of the feedforward control signal and the moment that both shafts obtain the same angular velocity (note that the mutual speed does not need to equal the initial input velocity). Fig. 11(b) illustrates the acceleration of the output shaft $\omega_{2}$ towards a synchronous speed with the input shaft. The figure clearly illustrates the discrepancy between the measured shifting time $y_{l}$ and corresponding simulation $Y_{l}$.

3.3 Practical set-up

The test set-up considered for this paper is depicted in Fig. 12. An AC electric motor ($30$ $\mathrm{kW}$) is controlled to a constant speed via a high bandwidth motor drive. The motor is connected to a controlled transmission via a torque converter. Within this transmission, the clutch can be controlled to engage with the proper gear. Furthermore, the output shaft of the clutch is connected to a load transmission. The load consists of a flywheel ($2,5$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{2}$) and a brake. This system was previously used to study smooth gearbox controllers [6].

3.4 Design of experiments

We experimentally tested the influence of several control and operating conditions on the shifting time $y_{l}$ on the wet clutch test set-up. We performed an up-shift from neutral to first gear for various control parameters and conditions, spanning a test scenario space. More specifically, 18 distinct feedforward control signals were tested. We compared 3 different values for $\Delta t$ and $u_{0}$ and 2 values for $\Delta u$. Further we verified 3 different motor speeds $\omega_{m}$ (1200, 1350 and 1500 $\mathrm{r}\mathrm{p}\mathrm{m}$) and 2 different friction load conditions $T_{b}$ (low, high). For these, a full factorial design of $l\in\{1,\cdots,108\}$ experiments has been tested. As an illustration of the modeling deficiency, we compare the actual measurements, $y_{l}$, with the simulated shifting times, $Y_{l}$, in Fig. 13. The global trend is captured by the deterministic model, however a distinct discrepancy between measurements and simulation is present.

3.5 Parametric uncertainty model

The objective of this research is to associate and identify a probabilistic model for several of the lumped model parameters. As mentioned in section 3.2, several of these parameters were identified empirically based on isolated subsystem measurements. The inertia of the shafts are derived from materials properties and geometrical considerations. Experiments in steady state regime are used to determine values for the viscous and Coulomb friction coefficients. On the available set-up no measurements of the oil temperature are possible and the temperature effect on the dynamics is therefore eliminated by performing all experiments after warm-up of the machine. However, our model contains various other variables that could not be determined or measured directly. Initially, we determined the parameter values using a least-squared procedure on the test scenario set. The comparison in Fig. 13 already illustrated the disadvantages of this strategy. Using a sensitivity analysis, we could identify two key parameters: $c_{2}$, the oil pressure bias in piston position computation (11) and the initial piston position $z_{p}$ (12). Therefore, we define the variables $x_{1}$ and $x_{2}$ that serve as scaling factor for respectively $c_{2}$ and $z_{p}$. We aim at identifying a probabilistic model for the random variable $\bm{\mathrm{X}}=(X_{1};X_{2})$. We propose a normal distribution, $\mathcal{P}_{\bm{\mathrm{X}}}(\bm{\mathrm{\alpha}})=\mathcal{N}(\bm{\mathrm{\mu}},\mathrm{\Sigma})$, with mean $\bm{\mathrm{\mu}}=(\mu_{1};\mu_{2})$ and covariance $\mathrm{\Sigma}$, corresponding the variable transformation, $\bm{\mathrm{X}}=\bm{\mathrm{\mu}}+\mathrm{\Sigma}\bm{\mathrm{\Theta}},\bm{\mathrm{\Theta}}\sim\mathcal{N}(0,\mathrm{I})$. Further we assume that both variables are independent so that $\mathrm{\Sigma}=\mathrm{diag}(\sigma_{1},\sigma_{2})$. Therefore, the probabilistic input model can be defined by $\bm{\mathrm{\alpha}}=(\mu_{1},\mu_{2},\sigma_{1},\sigma_{2})$. Note that a normal transformation implicitly assumes that forward model, $\mathcal{Y}$, can be evaluated for any value $\bm{\mathrm{x}}\in\mathbb{R}^{2}$. This can be avoided by clipping the distribution, constricting the stochastic domain of $\bm{\mathrm{X}}$ to a feasible area $\bm{\mathrm{x}}\in[\bm{\mathrm{x}}^{-},\bm{\mathrm{x}}^{+}]$, i.e. we are only interested in the part of the normal distribution in between these boundaries. This is facilitated by the transformation, $\mathrm{clip}(\bm{\mathrm{\Theta}})$. The random variables $\tilde{\bm{\mathrm{\Theta}}}_{i}=\mathrm{clip}(\bm{\mathrm{\Theta}}_{i})$ should be distributed according to $\frac{1}{\Phi(\theta_{i}^{+})-\Phi(\theta_{i}^{-})}\varphi(\tilde{\theta}_{i}),\tilde{\theta}_{i}\in[\theta_{i}^{-},\theta_{i}^{+}]$ where $\theta^{-}_{i}=\frac{1}{\sigma_{i}}(x^{-}_{i}-\mu_{i})$ and $\theta^{+}_{i}=\frac{1}{\sigma_{i}}(x^{+}_{i}-\mu_{i})$222Here $\varphi(\cdot)$ and $\Phi(\cdot)$ are standard notation for respectively the PDF and the cumulative density function (CDF) of the standard normal distribution.. First we map $\bm{\mathrm{\Theta}}_{i}$ to a uniform distribution in between the values $\Phi(\theta^{-}_{i})$ and $\Phi(\theta^{+}_{i})$. Then we use the inverse transform sampling method to map the uniform variable back onto the standard normal distribution but now limited to the interval $[\theta_{i}^{-},\theta_{i}^{+}]$. In conclusion the parametric probability model is defined through the transformations

[TABLE]

The clipping boundaries, $\bm{\mathrm{x}}^{-}$ and $\bm{\mathrm{x}}^{+}$, of the stochastic variables are chosen so that the space for which a feasible solution of the shifting time can be obtained for all 108 experiments, is maximized. Figures 14(a) and 14(b) illustrate respectively the parametric model output and associated clipped normal input distribution. Figure 14(c) demonstrates the corresponding output distribution for each experiment. Visual comparison with Fig. 13 demonstrates that the least-squared approach treats the problem in a non appropriate manner333We may note that the least-squares approach maximizes the objective $-\sum\frac{1}{2}(y_{l}-\mathcal{Y}_{l})^{2}$. Comparison with the $\log$-Likelihood approach reveals that this is equivalent with assuming that each experiment has an output probability density $\sim\exp(-\frac{1}{2}(y_{l}-\mathcal{Y}_{l})^{2})$. In other words the least-squares implicitly assumes that each experiment is normally distributed with variance $1$..

4 Results and discussion

In this section we apply the methodologies described in section 2 on the wet clutch system as was described in section 3. First we demonstrate the forward propagation of uncertainty by concatenating the techniques described in the sections 2.3 and 2.4, simply to determine the corresponding output distribution by moment matching for a given input distribution. Secondly we use the MLE method to identify the optimal parameter $\bm{\mathrm{\alpha}}^{*}$ that best describes observed experiments. We compare results obtained with both the Gaussian and max entropy parametric PDFs. To obtain a benchmark PDF we estimated the output distributions empirically using Monte Carlo (MC) uncertainty propagation of 200 000 simulation samples and fitted a spline curve on the associated histograms. The associated $\log L$ values are referred to as the true values.

4.1 PDF fitting by the method of moments

First we demonstrate the PDF estimation strategy combining high-order gPC moments with the method of moments. Corresponding the variable transformation and associated parametric probability model that were motivated in section 3.5, and in concordance with Table 1, we make use of Hermite polynomials. We used polynomial order $d=4$ and corresponding univariate quadrature order $q=5$, i.e. $25$ function simulations were required to estimate the output distribution for a single experiment. Due to numerical restrictions (see section 2.5), we are limited to the first $m=4$ moments.

Estimates for the output distribution of the 16th experiment are visualized in Fig. 15(a). We compare the MC reference distribution with Gaussian and max entropy fits. We also included PDF estimates that are based on the stochastic moments extracted from the MC sample set. Based on these results we can compare how much information could be extracted by the gPC approach with respect to the amount of information that is contained within the first four moments about the true distribution. The minor difference between the gPC PDF estimates and the MC equivalent PDF estimates suggests that gPC is a viable approach with sparse function evaluations.

For a fair comparison besides the visual verification, we calculate the Earth Mover’s Distance (EMD) (see Appendix B) for every estimate w.r.t. the estimated true output distribution. Different values are presented in Fig. 15(b). These results confirm that high-order estimates obtained through the extended gPC framework are capable of extracting statistical information from significantly fewer function evaluations when compared to the brute force methodologies such as MC.

4.2 Maximum log-likelihood estimation

We demonstrated that the high-order gPC moment matching strategy provides cheap yet accurate output PDF estimates for given input uncertainty. Consequently it can be engaged as an efficient alternative to evaluate the $\log$-Likelihood of a number of experiments at a mere fraction of the computational cost associated to brute force strategies such as MC. Next we apply the method to identify the optimal parameter set $\bm{\mathrm{\alpha}}^{*}$ associated to the probability model described in section 3.5. To solve problem (1), we used a genetic algorithm (GA) with population size 50 using the $\mathrm{Matlab}$ optimization toolbox. We compared optimization results for both Gaussian and max entropy PDF abstractions. The performance of the optimal parameter sets are verified using MC estimates of the $\log$-Likelihood corresponding the identified probability models.

4.2.1 Gaussian output distribution

Using the Gaussian output distribution and gPC moment estimates, the following optimal parameter set was identified, $\bm{\mathrm{\alpha}}^{*}_{G}=(0.91,0.55,0.077,0.014)$. Figure 16 compares measurement with the corresponding output distributions. The solution has an estimated $\log$-Likelihood of $\log\hat{L}^{G}(\bm{\mathrm{\alpha}}^{*}_{G})=188.5$ using the moment matching approach. When we evaluate the $\log$-Likelihood using the MC approach, we obtain a smaller value $\log L(\bm{\mathrm{\alpha}}^{*}_{G})=176.7$.

4.2.2 Maximum entropy output distribution

Using the max entropy distribution we obtain the following optimal parameter, $\bm{\mathrm{\alpha}}^{*}_{ME}=(0.91,0.55,0.081,0.046)$. Note that we retrieve the same mean as with the Gaussian output distributions, $\bm{\mathrm{\mu}}$, but that both covariance values have increased. Figure 17 visualizes the corresponding output distributions and measurements. Using the max entropy moment matching approach, the $\log$-likelihood is estimated on $\log\hat{L}^{ME}(\bm{\mathrm{\alpha}}^{*}_{ME})=189$. Verification using MC delivers a slightly smaller $\log$-Likelihood $\log L(\bm{\mathrm{\alpha}}^{*}_{ME})=183.7$. Note that the difference is significantly smaller compared to the $\log$-Likelihood mismatch obtained with the Gaussian PDF estimates.

4.3 Discussion

With respect to the global $\log L$ estimates, we remark the following. In both the max entropy and the Gaussian case the $\log L$ is overestimated using the proposed gPC moment matching method. However, for the max entropy distribution, the mismatch between true $\log L$ and estimated $\log L$ is significantly smaller (approximately 55%). Moreover the true $\log L$ obtained with $\bm{\mathrm{\alpha}}^{*}_{ME}$ exceeds that obtained with $\bm{\mathrm{\alpha}}^{*}_{G}$ by 4%. Both of these observations confirm that the use of high-order stochastic moments benefits the identification of an optimal input probability model and that the gPC framework can be used to obtain accurate estimates of these moments.

Figure 18 compares the estimates $\log\hat{L}_{l}^{G}(\bm{\mathrm{\alpha}}^{*}_{G})$ and $\log\hat{L}_{l}^{ME}(\bm{\mathrm{\alpha}}^{*}_{ME})$, for each individual experiment $l\in\{1,\cdots,108\}$, w.r.t. to their true value. One may note that for the majority of the experiments the $\log L$ is estimated with great precision for either output distribution. This indicates that for the majority of the experiments, the true output distribution is about Gaussian and no additional information is contained within the high-order moments. However for a small subset of experiments, the Gaussian overestimates the true value significantly due to its incapacity of modeling asymmetric distributions. As a result, measurements that are highly improbable and would strongly affect the global $\log L$ are penalized less by the Gaussian distribution. This effect was already visualized in Fig. 15.

Following the argument made above, we compare the EMD metric averaged over all experiments in Fig. 19, as it proofs to be valuable to quantify the estimated output distribution’s accuracy. We found that the use of the maximum entropy distribution with 4 moments could reduce the EMD by approximately 42% for $\alpha_{G}^{*}$ and 47% for $\alpha_{ME}^{*}$ when compared with the averaged EMD value obtained with Gaussian output distributions. The minor difference compared to the EMD of the PDF estimates based on the stochastic moments directly extracted from the MC sample set indicates that the gPC can correctly estimate the statistical moments. The improved accuracy can therefore be fully attributed to the higher order moments gPC framework.

Figure 20 illustrates the mean absolute error (MAE) of the $\log$-Likelihood estimates of all individual experiments $l\in\{1,\cdots,108\}$ compared to the MC reference obtained by 200 000 simulation samples. The accuracy of the brute force propagation techniques increases with respect to the number of samples used for uncertainty propagation. The quasi Monte Carlo (qMC) technique considers an equally distributed grid of $S$ samples. The gPC results require a magnitude less expensive function evaluations to achieve the same accuracy as the brute force techniques. The interpolated values $S_{MQ}$ and $S_{qMQ}$ are depicted in Table 2. A comparison between the gPC techniques indicates that the MAE could be reduced by 15.3% for $\bm{\mathrm{\alpha}}_{G}^{*}$ and 24.8% for $\bm{\mathrm{\alpha}}_{ME}^{*}$ due to the incorporation of higher order estimates in the novel gPC framework.

5 Conclusion

In this paper we proposed and discussed an efficient numerical method that is tailored to parametric model uncertainty identification using maximum likelihood estimation. Here the objective is to find a parametric input distribution that best explains a series of observations. A novel method is proposed to determine the output distribution for given input probability model making use of the generalized Polynomial Chaos (gPC) expansion framework. We extended the gPC framework so that high-order stochastic moments can be obtained efficiently. These high-order estimates can then be fed to a Gaussian or maximum entropy distribution. This strategy enables a reduction in the required number of function evaluations to obtain an adequate estimation of the statistical moments of the output distribution. The method is applied to calibrate a model for the shifting time for wet clutch engagement based on a series of measurements verified in the lab. It is shown that output distribution estimations via high-order moment matching excels traditional identification methods that are based on mean and variance estimates whilst the computational cost remains the same due to efficient techniques that were developed. The high-order moment matching resulted into a $\log$-likelihood increase of about 4% since the accuracy of the estimated output probability density function could be improved up to 47% compared to Gaussian distributions. This methodology reveals high potential when scaling up the number of uncertain parameters that need to be identified through inverse uncertainty identification. In addition to the inclusion of parametric uncertainty, future work will investigate the incorporation of model structure uncertainty.

Acknowledgments

Wannes De Groote holds a doctoral grant strategic basis research (3S07219) of the Fund for Scientific research Flanders (FWO). Some of the computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI. This study was performed in the frameworks of the EVIT ICON project. This research was supported by Flanders Make, the strategic research center for the manufacturing industry.

Appendix A Maximum entropy distribution

The maximum entropy distribution, $\hat{p}$ provides a solution to the classical problem where an estimate of univariate density $p$ is sought from knowledge of its moments [31]. The extent to which a density may be determined from its moments is a topic of discussion in mathematical literature. In practice only a finite number of moments, say $K+1$, are usually available so that there exists an infinite variety of functions whose first $K+1$ moments coincide and a unique reconstruction of $p$ is simply impossible. To remedy this issue, $\hat{p}$ is sought to maximize $\mathcal{W}[p]=-\int_{\mathcal{X}}p(x)\log p(x)\text{d}x$ under the condition that the first $K+1$ moments are equal to the true moments $\mu_{k}$. The Lagrangian $L$ corresponding to this problem is

[TABLE]

Nullifying the functional variation to $p$ yields

[TABLE]

Without loss of generality we may further assume that $\hat{p}$ is normalized so that $\mu_{0}=1$. Accordingly, we can express $\lambda_{0}$ in function of the remaining Lagrangian multipliers

[TABLE]

One may substitute this solution into the original Lagrangian (15) so to obtain the dual unconstrained problem which can be solved accordingly.

[TABLE]

Appendix B Earth mover’s distance

The Earth Mover’s distance (EMD) is a metric for the dissimilarity between two signatures that are each a compact representation of a distribution. In formal mathematics the metric is known as the Wasserstein distance between two probability distributions. The metric was coined in computer science due to an apparent analogy to its definition that can be modeled as the solution to a transportation problem [32].

Given are two signatures $P=\{(u_{i},\bm{\mathrm{x}}_{i})\}_{i=1}^{n}$ and $Q=\{(v_{j},\bm{\mathrm{y}}_{j})\}_{j=1}^{m}$. An insightful interpretation is that $u_{i}$ and $v_{j}$ represent the amount of dirt at the respective positions $\bm{\mathrm{x}}_{i}$ and $\bm{\mathrm{y}}_{j}$. The EMD between the signatures $P$ and $Q$ is defined as the minimal (normalized) work required to reconfigure $P$ into $Q$ moving around the dirt. Formally that is

[TABLE]

where $f_{ij}$ represents the dirt transferred from position $\bm{\mathrm{x}}_{i}$ to $\bm{\mathrm{y}}_{j}$, and, $d(\bm{\mathrm{x}}_{i},\bm{\mathrm{y}}_{j})$ the distance to be covered for that transport.

The transportation problem above is a special linear programming problem.

Appendix C Derivation of equation (9)

Details are given about the derivation of formula (9). We initiate the derivation from equation (5). First we expand the power of the sum using the multinomial theorem. Then we regroup the terms so that the high-order inner product is isolated from the coefficients. Finally we substitute expression $\prod_{j=1}^{n}\phi^{(j)}_{\uline{k}(j)}$ for $\psi_{\uline{k}}$ so to obtain equation (9).

[TABLE]