Model Bridging: Connection between Simulation Model and Neural Network

TL;DR

This paper introduces a 'model-bridging' framework that connects machine learning models with simulation models using kernel mean embeddings, enabling efficient, interpretable predictions without high computational costs, demonstrated in manufacturing simulations.

Contribution

The paper proposes a novel framework that efficiently links machine learning and simulation models through kernel mean embeddings, reducing computational costs and enhancing interpretability.

Findings

Framework successfully predicts outcomes and interpretable parameters.

Applied to manufacturing simulations like production and fluid dynamics.

Reduces computational costs compared to traditional calibration methods.

Abstract

The interpretability of machine learning, particularly for deep neural networks, is crucial for decision making in real-world applications. One approach is replacing the un-interpretable machine learning model with a surrogate model, which has a simple structure for interpretation. Another approach is understanding the target system by using a simulation modeled by human knowledge with interpretable simulation parameters. Recently, simulator calibration has been developed based on kernel mean embedding to estimate the simulation parameters as posterior distributions. Our idea is to use a simulation model as an interpretable surrogate model. However, the computational cost of simulator calibration is high owing to the complexity of the simulation model. Thus, we propose a ''model-bridging'' framework to bridge machine learning models with simulation models by a series of kernel mean embeddings to address these difficulties. The proposed framework enables us to obtain predictions and interpretable simulation parameters simultaneously without the computationally expensive calculations of the simulations. In this study, we apply the proposed framework to essential simulations in the manufacturing industry, such as production simulation and fluid dynamics simulation.

Tables

Table 1. Table 1: Comparison between machine learning models and simulation models.

	machine learning model	simulation model
interpretability of parameter	un-interpretable	interpretable
computational cost of the model	not expensive	expensive

Algorithm 1: Model bridging
1) Pre-learning for each dataset :
Input: Dataset ${\{X_l^n,Y_l^n\}}_{l=1}^L$,
machine learning model $f_{\mathrm{ml}}(x,\xi )$
and simulation model $f_{\mathrm{sim}}(x,\theta )$
Output: $\{{\hat{\mu }}_1^{\mathrm{ml}},\dots ,{\hat{\mu }}_L^{\mathrm{ml}}\}$ and $\{{\hat{\mu }}_1^{\mathrm{sim}},\dots ,{\hat{\mu }}_L^{\mathrm{sim}}\}$
for $l=1$ to $L$ do
Estimation for ${\hat{\mu }}_l^{\mathrm{ml}}$ by Eq. (7)
Estimation for ${\hat{\mu }}_l^{\mathrm{sim}}$ by Eq. (1)
end for
2) Learning for model-bridging function $\hat{T}$:
Input: $\{{\hat{\mu }}_1^{\mathrm{ml}},\dots ,{\hat{\mu }}_L^{\mathrm{ml}}\}$ and $\{{\hat{\mu }}_1^{\mathrm{sim}},\dots ,{\hat{\mu }}_L^{\mathrm{sim}}\}$
Output: Model-bridging function $\hat{T}$
Learning for $\hat{T}$ by Eq. (3)
3) Prediction:
Input: Dataset $\{X_{L+1}^n,Y_{L+1}^n\}$ and $X_{L+1,n+1}$
Output: ${\hat{Y}}_{L+1,n+1}$ and ${\hat{\theta }}_{L+1}$
Estimation for ${\xi }_{L+1}$ and ${\hat{\mu }}_{L+1}^{\mathrm{ml}}$
Prediction for ${\hat{Y}}_{L+1,n+1}$ by $f_{\mathrm{ml}}(x;{\xi }_{L+1})$
Estimation for ${\hat{\mu }}_{L+1}^{\mathrm{MB}}=\hat{T}\left({\hat{\mu }}_{L+1}^{\mathrm{ml}}\right)$ by Eq. (5)
Sampling for ${\hat{\theta }}_{L+1}^{\mathrm{MB}}$ by Eq. (6)

Table 3. Table 2: Summary of true and estimated parameters in the experiment for realistic production simulation. T BF subscript 𝑇 BF T_{\rm BF} represents the mean time between failures, and T R subscript 𝑇 𝑅 T_{R} represents the mode time of repair for each process. The estimated parameters are mean and standard deviation (in parentheses) of posterior mean for one-leave-out cross-validation.

Process	Saw		Coat		Inspection		Harden		Grind		Clean
	$T_{\mathrm{BF}}$	$T_R$	$T_{\mathrm{BF}}$	$T_R$	$T_{\mathrm{BF}}$	$T_R$	$T_{\mathrm{BF}}$	$T_R$	$T_{\mathrm{BF}}$	$T_R$	$T_{\mathrm{BF}}$	$T_R$
Param.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	${\theta }_9$	${\theta }_{10}$	${\theta }_{11}$	${\theta }_{12}$
${\theta }^{\left(0\right)}$	100	25	150	5	100	20	150	5	75	15	120	20
${\theta }^{\left(1\right)}$	100	25	150	5	80	20	150	5	75	15	120	20
${\hat{\theta }}^{\mathrm{sim}}$	100.5	25.1	153.2	5.0	104.2	17.3	146.5	5.1	73.7	14.6	96.8	20.2
$\left(x\le 30\right)$	(9.5)	(2.2)	(14.2)	(0.4)	(10.1)	(1.9)	(14.2)	(0.3)	(6.7)	(1.2)	(15.7)	(2.3)
${\hat{\theta }}^{\mathrm{sim}}$	100.5	24.7	153.6	5.0	84.3	22.8	148.2	5.0	73.8	15.2	115.4	20.4
$\left(x>30\right)$	(8.3)	(1.4)	(10.6)	(0.3)	(7.4)	(0.9)	(8.9)	(0.3)	(4.6)	(1.0)	(9.8)	(1.8)
${\hat{\theta }}^{\mathrm{MB}}$	102.5	21.2	165.1	6.0	98.6	17.0	179.0	5.5	74.0	17.8	94.6	12.7
$\left(x\le 30\right)$	(12.7)	(1.8)	(12.6)	(0.4)	(9.9)	(1.8)	(12.1)	(0.3)	(6.8)	(1.6)	(14.5)	(2.3)
${\hat{\theta }}^{\mathrm{MB}}$	104.1	21.1	165.3	6.1	86.7	19.1	180.8	5.4	74.0	18.1	99.3	13.2
$\left(x>30\right)$	(10.1)	(2.3)	(14.2)	(0.5)	(11.8)	(2.1)	(8.1)	(0.5)	(4.5)	(1.2)	(21.2)	(1.6)

Equations

$$\text{vspace}-0.1cmp(y|x,\theta )=\frac{1}{\sqrt{2\pi {\sigma }_0^2}}\exp \left\{-\frac{1}{2{\sigma }_0^2}{\Vert y-f_{\mathrm{sim}}(x;\theta )\Vert }^2\right\},\text{vspace}-0.1cm$$

$${\hat{\mu }}_{\theta |YX}=\sum _{j=1}^m{w}_j{k}_{\theta }(\cdot ,{\theta }_j)$$

$$(w_1,...,w_m{)}^T$$

$$G_y$$

$${\mathbf{k}}_y(Y^n)$$

$$\hat{T}=\arg \underset{T\in \mathcal{F}}{\max }\frac{1}{L}\sum _{l=1}^L\Vert {\hat{\mu }}_l^{\mathrm{sim}}-T({\hat{\mu }}_l^{\mathrm{ml}}){\Vert }_{\mathcal{F}}^2+\lambda \Vert T{\Vert }_{\mathcal{F}}^2,$$

$$\kappa (\mu ,{\mu }^{\prime })=\exp \left\{-\frac{1}{2{\sigma }_{\mu }^2}{\Vert \mu -{\mu }^{\prime }\Vert }_{\mathcal{H}}^2\right\}\in \mathcal{F},$$

$${\hat{\mu }}_{L+1}^{\mathrm{MB}}=\hat{T}({\hat{\mu }}_{L+1}^{\mathrm{ml}})=\sum _{l=1}^L{v}_l{\hat{\mu }}_l^{\mathrm{sim}}\in \mathcal{F},$$

$$G_{\mu }$$

$${\mathbf{k}}_{\mu }({\hat{\mu }}_{L+1}^{\mathrm{ml}})$$

$${\hat{\theta }}_{L+1,j}=\arg \underset{\theta }{\max }\sum _{l=1}^L\sum _{j^{\prime }=1}^m{v}_l{w}_{l,j^{\prime }}{k}_{\theta }(\theta ,{\theta }_{l,j^{\prime }})-\frac{1}{j}\sum _{j^{\prime }=1}^{j-1}{k}_{\theta }(\theta ,{\theta }_{j^{\prime }})\in {\mathbb{R}}^{d_{\theta }},$$

$$\kappa ({\hat{\mu }}_l^{\mathrm{ml}},{\hat{\mu }}_{l^{\prime }}^{\mathrm{ml}})=\exp \left\{-\frac{1}{2{\sigma }_{\mu }^2}{\Vert {\hat{\mu }}_l^{\mathrm{ml}}-{\hat{\mu }}_{l^{\prime }}^{\mathrm{ml}}\Vert }_{\mathcal{H}}^2\right\}\in \mathcal{F}$$

$$=\exp \left\{-\frac{1}{{\sigma }_{\mu }^2}\left(1-\sum _{j=1}^m\sum _{j^{\prime }=1}^{m^{\prime }}{k}_{\xi }({\xi }_{l,j},{\xi }_{l^{\prime },j^{\prime }})\right)\right\}.$$

$$y={\hat{\mu }}_{Y|X,l}(x)=\sum _{i=1}^n{u}_{l,i}(x)k_y(\cdot ,Y_{l,i}),$$

$$\kappa ({\hat{\mu }}_l^{\mathrm{ml}},{\hat{\mu }}_{l^{\prime }}^{\mathrm{ml}})=\exp \left\{-\frac{1}{{\sigma }_{\mu }^2}\left(1-\sum _{i=1}^n\sum _{i^{\prime }=1}^{n^{\prime }}{u}_{l,i}{u}_{l^{\prime },i^{\prime }}{k}_y(Y_{l,i},Y_{l^{\prime },i^{\prime }})\right)\right\}.$$

$$\Vert {\hat{\mu }}_{L+1}^{\mathrm{MB}}-{\hat{\mu }}_{L+1}^{\mathrm{sim}}{\Vert }_{\mathcal{H}}^2$$

Full text

1 Introduction

The interpretability of machine learning, especially for deep neural networks, is crucial for decision making in real-world applications. In recent years, many studies have addressed the interpretability of neural networks [6, 4, 15]. One of the approaches is replacing the un-interpretable machine learning model with a surrogate model, which has a simple structure for interpretation. This approach is a type of model compression. For instance, the “distillation” of a neural network model [7] is one of the representative methods for model compression for replacing a complex model with a simplified model; however, there is no interpretability for a small surrogate neural network model. There are some methods to obtain an interpretable model, such as LIME [20], SHAP [14], and a method combined with a rule-based model [22]. These methods do not provide a clear pathway toward obtaining the interpretability of a neural network, as there are limitations to obtain local interpretability regarding the decision boundary of the prediction result [6, 4].

Another approach for understanding the target system is by employing a simulation that might be outside the scope of conventional machine learning. In some application fields, simulations such as multi-agent simulation, traffic simulation, production simulation, or simulation of the dynamics of the physical system have already been used to understand the target system and to predict future behavior. Simulation modeling is implemented to describe the fundamental law of the objective system, using human knowledge with interpretable simulation parameters. The recently developed “simulator calibration” [9, 10, 3] is a method in which the simulation parameters are estimated as posterior distributions in the context of machine learning. Simulator calibration can provide a predictive result with interpretable simulation parameters. Our idea is to use a simulation model as an interpretable surrogate model. However, the difficulty of simulator calibration is attributed to a substantial computational cost; it typically takes more than one hour owing to the complexity of the simulation model (Table 1). In real-world applications, a predictive result and its reason often should be required to obtain within a minute.

We propose a “model-bridging (MB) ” framework to predict using a machine learning model as well as obtain interpretable simulation parameters simultaneously without expensive calculation of a simulation model. The idea of this framework is to map the un-interpretable parameters of the machine learning model and the interpretable parameters of the simulation model (Fig. 1). The algorithm has to learn the relation of the posterior distribution estimated from each dataset between the machine learning model and the simulation model in advance; this framework can be considered as a meta-learning for each dataset.

Let us consider the example of production simulation for predicting the efficiency of manufacturing production, implementing a series of processes for production (example in Fig. 4). The production simulation aims to obtain a production efficiency and the reason for it simultaneously within a minute to improve the production efficiency. We formulate this problem setting. Assume that we obtain the dataset $\{X^{n},Y^{n}\}=\{X_{1},...,X_{n},Y_{1},...,Y_{n}\}$, where input $X_{i}\in\mathbb{R}^{d_{x}}$ is the number of products to be manufactured in unit time and output $Y_{i}\in\mathbb{R}^{d_{y}}$ is the efficiency of production. The simulation parameter $\theta\in\mathbb{R}^{d_{\theta}}$ is the elapsed time for each process, which undergoes a probabilistic behavior. The parameter $\theta$ is interpretable and helpful in understanding the system and decision making. Thus, we need to obtain the prediction $\hat{Y}_{n+1}$ for new data $X_{n+1}$ as well as obtain the interpretable simulation parameters $\theta$ representing the elapsed time of each process, which provides information regarding the occurrence of “bottleneck processes.” Here, the observed data and its generation process are considered to drift gradually, for example, the daily production of the factory due to the load of labors and machine environment factors such as temperature. The detailed assumption is described in a later section.

Note that this study considered a different problem setting from the conventional methods with simplified surrogate models, such as LIME, SHAP, and rule-based model; the interpretable model of the proposed method, i.e., simulation model, is complex and computationally expensive. There is no existing method for solving this new problem setting, where it is difficult to show the baseline for the evaluation. Experimentally, we confirm that the estimation of model bridging is reasonable in comparison with simulator calibration as a baseline with a significantly fast process owing to no execution of the simulation.

The main contribution of this paper is to propose a novel framework for bridging machine learning and simulation, which has never been discussed before from the context of machine learning and to demonstrate its effectiveness in real-world applications. The technical contribution is to expend the distribution-to-distribution regression on reproducing kernel Hilbert space (RKHS), as a suitable method for bridging function. The rest of this paper is organized as follows. We briefly review a series of applications of kernel mean embedding as the building blocks for the proposed framework. Subsequently, we propose the model-bridging framework. Finally, we confirm the accuracy of the proposed method for three cases of simulation.

2 Related Works

We briefly introduce simulator calibration and distribution regression based on kernel mean embedding as a building block of the proposed framework.

2.1 Simulator Calibration

“Simulator calibration” [10] is a method for estimating the simulation parameter as the posterior distribution to reproduce real data. Simulator calibration is an example of data assimilation. The simulation model is treated as a regression function $f_{\rm sim}(x;\theta)$ by combining a series of kernel mean embedding methods. The conventional statistical methods of parameter estimation are not applicable to simulator calibration owing to the properties of the likelihood function: intractable or nondifferentiable. When Gaussian noise is employed with regression function $f_{\rm sim}(x;\theta)$, the likelihood is expressed as

[TABLE]

where $\sigma_{0}>0$ is a constant of observation noise. This likelihood function is nondifferentiable owing to the simulation model $f_{\rm sim}(x;\theta)$. The posterior mean to be obtained is formulated as $p(\theta|X^{n},Y^{n})=p(Y^{n}|X^{n},\theta)\pi(\theta)/Z(X^{n},Y^{n})$, where $\pi(\theta)$ is the prior distribution and $Z(X^{n},Y^{n})$ is the regularization constant. In this application, simulator calibration estimated the simulation parameter $\theta$ as a kernel mean of the posterior distribution by using kernel approximated Bayesian computation (kernel ABC) [17, 5]. After obtaining the kernel mean of the posterior distribution, a posterior sample is obtained using kernel herding [1].

2.2 Application of Kernel Mean Embedding

As an application of kernel mean embedding [16], we briefly review the kernel ABC and kernel herding. The kernel mean embedding is a framework for mapping distributions into a RKHS $\mathcal{H}$ as a feature space. Kernel herding is a sampling method from the embedded distribution in RKHS that has the opposite operation of kernel mean embedding. Figure 2 shows a schematic of the relation of kernel ABC and kernel herding.

Kernel ABC: Kernel ABC [17, 5] is a method for computing the kernel mean of the posterior distribution from a sample of parameter $\theta$, generated by the prior distribution $\pi(\theta)$. The assumption is that the explicit form of the likelihood function is unavailable, while the sample from the likelihood function is available. The kernel ABC allows us to calculate the kernel mean of the posterior distribution as follows. First, the sample $\{\theta_{1},...,\theta_{m}\}$ is generated from prior distribution $\pi(\theta)$ and pseudo-data $\{\bar{Y}_{1}^{n},...,\bar{Y}_{m}^{n}\}$, as a sample from $p(y|x,\theta_{j})$ for $j=1,...,m$. Next, the empirical kernel mean of the posterior distribution

[TABLE]

is calculated, where $k_{\theta}$ is a kernel of $\theta$. Weight $w_{j}$ is calculated by

[TABLE]

The $\delta\geq 0$ is a regularization constant, $I$ is an identity matrix, and $k_{y}$ is a kernel of $y$. The kernel $k_{y}(\bar{Y}^{n}_{j},Y^{n})$ indicates the “similarity” between pseudo-data $\bar{Y}_{j}^{n}$ and real data $Y^{n}$. The calculations of the kernel mean corresponds to the estimation of the posterior distribution as an element in $\mathcal{H}$.

Kernel Herding: Kernel herding [1] is a method used for sampling data from the kernel mean representation of a distribution, which is an element of the RKHS. Kernel herding can be considered as an opposite operation to that of kernel ABC. Kernel herding greedily obtains samples of $\theta$ by updating Eqs.(1) and (2) as given in Chen et al. [1].

2.3 Distribution Regression

Distribution regression is a regression for $d_{x}$-dimensional “distributions” represented by samples. In contrast, normal regression is regression for $d_{x}$-dimensional “point.” There are several studies of distribution regression, including distribution-to-distribution regression [19] and distribution-to-point regression [21, 12]. Oliva et al. [19] employed the idea of approximating a density function by kernel density estimation, rather than using RKHS. Szabó et al. [21] proposed the distribution-to-point with the kernel ridge regression method on RKHS; however, no methods are available for distribution-to-distribution regression on RKHS.

3 Proposed Framework: Model Bridging

We propose a novel framework to bridge the un-interpretable machine learning model and the interpretable simulation model. In this study, we assume a machine learning model, such as a Bayesian neural network (BNN) [18], as a parametric model and a Gaussian process as a non-parametric model. This proposed framework is applicable to any model. In this section, first, we confirm the problem setting and framework of model bridging. Second, we propose the algorithm of distribution-to-distribution regression, which is suitable for the proposed framework. Thereafter, we propose the formulation of the input of distribution-to-distribution regression for the parametric model, assuming BNN, and the non-parametric model, assuming the Gaussian process. Figure 3 and Alg. 1 shows an overview of the framework.

3.1 Problem Setting, Assumption, and Usage of Model Bridging

We define the problem setting of the model-bridging framework. Let $L$ be dataset $\{X_{l}^{n},Y_{l}^{n}\}_{l=1}^{L}$ $(X^{n}_{l}\in\mathbb{R}^{n\times d_{x}},Y^{n}_{l}\in\mathbb{R}^{n\times d_{y}})$, given in the pre-learning phase. For simplicity of explanation, we use the unique number of the data $n$ and sample size $m$ for all datasets. However, it can be different numbers generally, such as $n_{l}$ and $m_{l}$. The purpose is to predict $\hat{Y}_{L+1,n+1}$ and simultaneously obtain interpretable simulation parameter $\hat{\theta}^{\rm MB}_{L+1}$ to reproduce $Y_{L+1,n+1}=f_{\rm sim}(X_{L+1,n+1};\hat{\theta}^{\rm MB}_{L+1})$ without the expensive calculation of simulation model $f_{\rm sim}(x;\theta)$, when we obtain new dataset $\{X^{n}_{L+1},Y^{n}_{L+1}\}$. The assumptions of the problem setting are as follows. These assumptions are prevalent for many applications of a simulation.

• •

The existing simulation model $f_{\rm sim}(x;\theta)$ with interpretable simulation parameter $\theta\in\mathbb{R}^{d_{\theta}}$ and a machine learning model $f_{\rm ml}(x;\xi)$ that is sufficiently accurate to predict a typical regression problem while having un-interpretable parameter $\xi\in\mathbb{R}^{d_{\xi}}$.

• •

The cost of simulator calibration is much higher than that of learning for the machine learning model. For instance, it takes more than one hour for simulator calibration of one dataset $\{X^{n}_{l},Y^{n}_{l}\}$, while learning of BNN takes less than a minute.

• •

Dataset $\{X_{l}^{n},Y_{l}^{n}\}$ has dependency of parameter $\theta_{l}$ for each $l=1,...,L$. Let us assume the following situation. $\{X^{n}_{l},Y^{n}_{l}\}$ is obtained in one day with the same conditions, described as parameter $\theta_{l}$, while the conditions are changed for the following day, described as $\theta_{l+1}$.

• •

The time for offline calculation of simulator calibration is sufficient, while the time for prediction is restricted.

Once we obtain the model-bridging function as a mapping from the machine learning model to the simulation model, we can obtain an accurate prediction for $\hat{Y}_{L+1,n+1}$ by both the machine learning model and interpretable $\hat{\theta}^{\rm MB}_{L+1}$ by the simulation model for new dataset $\{X^{n}_{L+1},Y^{n}_{L+1}\}$ without an expensive calculation from the simulation model.

3.2 Distribution-to-Distribution Regression

We present the regression algorithm between the conditional kernel mean of the machine learning model $\mu^{\rm ml}\in\mathcal{H}$ and that of the simulation model $\mu^{\rm sim}\in\mathcal{H}$, as a model-bridging function $\mu^{\rm sim}=T(\mu^{\rm ml})$. We develop the algorithm based on kernel ridge regression, which is suitable for kernel mean input and output on RKHS. This is the extension of the distribution-to-point regression method proposed by Szabó et al. [21] for the distribution output.

3.2.1 Kernel Ridge Regression for Kernel Mean

The formulation to be solved as an analogy of normal kernel ridge regression is as follows:

[TABLE]

where $\lambda\geq 0$ is a regularization constant. $\mathcal{F}$ is a function space of kernel mean embeddings following Christmann et al. [2] and $\|\cdot\|_{\mathcal{F}}$ is its norm. The difference from ordinary kernel ridge regression is that the inputs and outputs are kernel means. Therefore, we define kernel $\kappa\in\mathcal{F}$, as a function of kernel mean $\mu\ \in\mathcal{H}$. We employ a Gaussian-like kernel as

[TABLE]

where constant $\sigma_{\mu}>0$ is the width of kernel $\kappa$ and $\|\cdot\|_{\mathcal{H}}$ is RKHS norm. The kernel $\kappa$ is also a positive definite kernel [2].

Following the representor theorem of kernel ridge regression [11], the estimated model-bridging function $\hat{T}$ for new $\hat{\mu}^{\rm ml}_{L+1}$ is described as

[TABLE]

where ${\bf v}=(v_{1},...,v_{L})^{T}=(G_{\mu}+\lambda LI)^{-1}{\bf k}_{\mu}(\hat{\mu}^{\rm ml}_{L+1})\ \ \in\mathbb{R}^{L}$. Gram matrix $G_{\mu}$ and the vector ${\bf k}_{\mu}(\hat{\mu}^{\rm ml}_{L+1})$ are described as follows:

[TABLE]

3.2.2 Kernel Herding from Kernel Mean $\hat{\mu}^{\rm MB}_{l}$

After obtaining the kernel mean of $\hat{\mu}^{\rm MB}_{L+1}$, kernel herding can be applied to sample $\hat{\theta}^{\rm MB}_{L+1}=\{\hat{\theta}_{L+1,1},...,\hat{\theta}_{L+1,m}\}$, where $\hat{\theta}_{L+1,j}\in\mathbb{R}^{d_{\theta}}$. The explicit form of the update equation for sample $j=1,...,m$ iteration of kernel herding with kernel mean $\hat{\mu}^{\rm MB}_{L+1}$ is as follows:

[TABLE]

for $j=2,...,m$. For initial state $j=1$, the update equation constitutes only the first term of Eq. (6). The weight of $w_{l,j}$ is calculated by kernel ABC for dataset $\{X^{n}_{l},Y^{n}_{l}\}$ in Eq.(2).

3.3 Input of Distribution-to-Distribution Regression

We present the explicit formulation for calculating the kernel means of the machine learning model $\hat{\mu}^{\rm ml}_{l}$, as an input of the distribution-to-distribution regression. First, we present the formulation of BNN, as BNN is a useful model for many applications as a parametric Bayesian model. Second, we present the formulation for Gaussian process regression as a non-parametric Bayesian model. We consider the Gaussian process regression as a non-parametric alternative to BNN. The equivalence between the Gaussian process and BNN with one hidden layer with infinite nodes is well known [18]. Furthermore, a recent study reveals the kernel formulation that is equivalent to multi-layered BNN, as an extension of the Gaussian process [13]. We directly obtain empirical kernel mean without calculation of kernel from parameters in the parametric model for using Gaussian process regression.

3.3.1 Parametric Model: Bayesian Neural Network

We assume the BNN model $f_{\rm ml}(x;\xi)$ with a few hidden layers, where $\xi$ is a parameter, such as weights for each node and bias terms of each layer. We can obtain the posterior distribution of $\xi_{l}$ for $l=1,...,L$ by the Markov Chain Monte Carlo method or variational approximation. Then, the empirical kernel mean of the posterior distribution is represented as $\hat{\mu}^{\rm ml}_{l}=\sum_{j=1}^{m}k_{\xi}(\cdot,\xi_{l,j})\ \in\mathcal{H}$ for $l=1,...,L$ dataset, where $k_{\xi}$ is kernel of $\xi$.

We employ Gaussian-like kernel $\kappa$ as an function of $\hat{\mu}^{\rm ml}_{l}\in\mathcal{H}$ as

[TABLE]

The relation $\langle\hat{\mu}^{\rm ml}_{l},\hat{\mu}^{\rm ml}_{l^{\prime}}\rangle=\sum_{j=1}^{m}\sum_{j^{\prime}=1}^{m^{\prime}}k_{\xi}(\xi_{l,j},\xi_{l^{\prime},j^{\prime}})$ is used, where $\langle\cdot,\cdot\rangle$ represents the inner product.

3.3.2 Non-Parametric Model: Gaussian Process Regression

We use the Gaussian process regression as a non-parametric model. In this case, we can directly obtain empirical kernel mean without the calculation of kernel from parameters, such as $\xi$ in the parametric model. We present that the prediction with Gaussian process regression can be considered as a conditional kernel mean. As a result of Gaussian process regression, we can express the mean of the predictive distribution in general for $l$-th dataset $\{X^{n},Y^{n}\}$ as

[TABLE]

where $u_{l,i}(x)=\{(G_{x}+n\lambda^{\prime}I)^{-1}{\bf k}_{x}(x)\}_{i}$. The $G_{x}$ is the Gramm matrix, $\lambda^{\prime}\geq 0$ is regularization constant, and $k_{x}$ is kernel of $x$. This formulation is clear if we remember the equivalence between Gaussian process regression and kernel ridge regression [8]. As a predictor of $\hat{Y}_{l,n+1}$ for new $X_{l,n+1}$, we can calculate $\hat{Y}_{l,n+1}=\hat{\mu}_{Y|X,l}(X_{l,n+1})$. Note that this $\hat{\mu}_{Y|X,l}$ is interpreted as the kernel mean of the $Y^{n}_{l}$ conditioned by $X^{n}_{l}$. Thus, we can use $\hat{\mu}_{Y|X,l}$ as the input of the distribution-to-distribution regression, represented as $\hat{\mu}^{\rm ml}_{l}$. We employ Gaussian-like kernel $\kappa$ as an function of $\hat{\mu}^{\rm ml}_{l}$ as

[TABLE]

There is a difference between the proposed non-parametric method and parametric method. In the parametric method, the distribution of parameter $\xi_{l}$ is the input of the distribution-to-distribution regression, while in the non-parametric method, the distribution of data $Y^{n}_{l}$ conditioned by $X^{n}_{l}$ is the input.

4 Experiment

This framework is widely applicable to various domains for industries that include multi-agent simulation, traffic simulation, and simulation of dynamics of physical systems, such as thermomechanics, structural mechanics, and electromagnetic mechanics. We present the applications of the model-bridging framework for three simulations: 1) a simple production simulation to show and explain the framework effectiveness in detail, 2) a realistic production simulation to show the capability for realistic application, and 3) a simulation of fluid dynamics to show of the applicability to a wide variety of simulation fields. The detailed information on the three experiments is provided in the supplemental material owing to page limitations.

Through the three experiments, we confirm that model bridging enables us to predict $\hat{Y}_{L+1,n+1}$ for new $X_{L+1,n+1}$, using a machine learning model, and obtain the interpretable simulation parameter $\hat{\theta}^{\rm MB}_{L+1}$ without an expensive calculation of the simulation when we obtain the $L+1$-th dataset $\{X^{n}_{L+1},Y^{n}_{L+1}\}$. We also investigate the accuracy of the estimation of the parameter $\hat{\theta}^{\rm MB}_{L+1}$ and compared the execution time with simulator calibration as a baseline. Note that these experiments cannot be compared with other state-of-the-art surrogate approaches, such as LIME, SHAP, and other methods with model compression, because no other methods exist that simultaneously obtain the prediction result and interpretable simulation parameters.

4.1 Common Setting of Experiments

In practice, the effective hyperparameter for model-bridging function to be tuned is the regularization constant $\lambda$ for distribution-to-distribution regression. The hyperparameter $\lambda$ can stabilize the calculation of the inverse Gram matrix. This hyperparameter should be determined by cross-validation. Further, as a common hyperparameter of the kernel method, the width of the kernel must be selected to measure the similarity between the data. We employed a Gaussian kernel for $k_{y}$, $k_{x}$, $k_{\theta}$, and $k_{\xi}$ for all experiments. The typical setting of the width of the kernel, in practice, is the median of Euclid distance of the input data of a kernel. In all the experiments performed in this study, we apply this setting and confirmed that all kernels perform adequately. We used a PC equipped with a 3.4-GHz Intel core i7 quad-core processor and 16GB memory. The main computational cost is for a simulation in the pre-learning phase of this framework.

4.2 Experiment with Simple Production Simulator

4.2.1 Setting

Production simulators are widely used simulation software for discrete and interconnection systems to model various processes, such as production, logistics, transportation, and office works. We used a WITNESS, a popular software package of production simulation111https://www.lanner.com/en-us/technology/witness-simulation-software.html. We examined the regression problem using a simulation model that has a simple four-dimensional simulation parameter $\theta\in\mathbb{R}^{4}$ (Fig. 4). We defined the simulation input $x=X_{i}\in\mathbb{R}$ as the number of products to be manufactured, output $Y_{i}=f_{\rm sim}(X_{i},\theta)\in\mathbb{R}$ as the total time to manufacture all the $X_{i}$-th products, and parameter $\theta$ as the time required for each procedure on the production line. Moreover, the time required for “ASSEMBLY” is $\mathcal{N}(\theta_{1},\theta_{2})$, and that for “INSPECTION” is $\mathcal{N}(\theta_{3},\theta_{4})$, where $\mathcal{N}(\mu_{\rm ND},\sigma_{\rm ND})$ is the normal distribution with mean $\mu_{\rm ND}$ and standard deviation $\sigma_{\rm ND}$. We assumed that the elapsed time of each process would increase considerably, owing to an increasing load, if the number of products to be manufactured also increases. To create this situation artificially, we set different true parameters between the observed data region $\theta^{(0)}$ and the predictive region $\theta^{(1)}$. We set $\theta^{(0)}=(2,0.5,5,1)^{T}$ if $x\leq 110$ and $\theta^{(1)}=(3.5,0.5,7,1)^{T}$ if $x>110$. The shift in parameters $\theta_{1}$ and $\theta_{3}$ between $\theta^{(0)}$ and $\theta^{(1)}$ is the sigmoid function. For each $l$-th dataset, the observed data of size $n=50$ and sample size $m=100$ are generated by $q_{l}(x)=\mathcal{N}(\chi_{l},5)$, where $\chi_{l}$ is uniform distribution in $[70,130]$ for $l=1,...,L$. The number of training datasets, $L$, is 100. We defined the prior distribution as the uniform distribution over $[0,5]\times[0,2]\times[0,10]\times[0,2]$. We used a BNN having two fully connected hidden layers with three nodes and bias nodes for each layer, as a machine learning model. The activation function is ReLU. The regularization constant is $\lambda=1.0\times 10^{-6}$ for this experiment.

4.2.2 Result

The execution time of model bridging is 9.6 [s] in the presented computational environment, while simulator calibration requires about 3.1 [h] for $L+1$-th dataset. Simulator calibration requires $m\times n$ execution of simulation and each simulation takes 2 [s] in this case. As representatives of all test datasets, the results of the model-bridging framework for two different datasets ($l$-th and $l^{\prime}$-th dataset), which are randomly selected, are shown in Fig. 5. Figure 5 (A) shows the observed data for $l$-th dataset as red squares and the $l^{\prime}$-th dataset as red triangles. The solid line and dashed line are the fitted results by BNN with variational approximation. Figure 5 (B) shows the estimated posterior distributions of simulation parameters by model bridging $\hat{\theta}^{\rm MB}_{l}$ each dataset. The red markers show the true parameter, green markers show the estimated result of simulator calibration, and blue markers show the mean of the distribution. Each square denotes the $l$-th dataset and each triangle denotes the $l^{\prime}$-th dataset. We can see a reasonably accurate estimation of $\theta_{l}^{\rm MB}$ and $\theta_{l^{\prime}}^{\rm MB}$ model-bridging framework in comparison with simulator calibration for the two different datasets with two different true parameters $\theta$.

Note that from the perspective of interpretability, we can clearly see the practical effectiveness of simultaneously obtaining the prediction result with interpretable parameters, such as “elapsed time of a process.” From these interpretable parameters, we can understand that the production efficiency is decreased ($l^{\prime}$-th dataset in Fig.5) mainly because of the increased elapsed time of “INSPECTION” ($=\theta_{1}$).

We also investigated $\|\hat{\mu}^{\rm MB}_{L+1}-\hat{\mu}^{\rm sim}_{L+1}\|^{2}_{\mathcal{H}}$ to confirm the convergence of the proposed distribution-to-distribution regression in the model-bridging framework. The detailed formulation for the numerical calculation is presented in the supplementary material. Figure 6 shows the mean and standard deviation of one-leave-out cross-validation of the test dataset. The horizontal axis shows the number of training datasets. We can see the convergence for bias that originates from simulator calibration.

4.3 Experiment with Realistic Production Simulator

4.3.1 Setting

We used a model to reproduce a real metal-processing factory that manufactures valves from metal pipes, with six primary processes: “saw,” “coat,” “inspection,” “harden,” “grind,” and “clean,” in the order shown in the supplementary material. Each process is composed of complex procedures, such as the preparation rule, waiting, and machine repair during trouble. The purpose of this production simulation is also to predict the total production time $Y_{i}\in\mathbb{R}$ when the number of units $X_{i}\in\mathbb{R}^{3}$ for three types of products to be manufactured is set. Each of the six processes contains two parameters of machine downtime owing to failure: mean time between failures ($T_{\rm BF}$) and mode time required for repair ($T_{\rm R}$). We defined these parameters as twelve-dimensional parameter $\theta\in\mathbb{R}^{12}$ (see Table 2). The distribution of the mean time between failures is represented as a negative exponential distribution. The distribution of the time required for repair is represented as an Erlang distribution with the mode time and shape parameter set at three.

Similar to the simple experiment discussed in the previous section, we set the true parameter as $\theta^{(0)}$ if $x\leq 30$ and $\theta^{(1)}$ if $x>30$. The summary of the true parameter is shown in Table 2. The shift in parameter $\theta_{5}$ between $\theta^{(0)}$ and $\theta^{(1)}$ is a sigmoid function. The observed data of size is $n=30$ by $\mathcal{N}(\chi_{l},3)$ where $\chi_{l}$ is generated by uniform distribution in $[20,40]$. The number of parameter samples is $m=50$ and the number of datasets is $L=40$. We defined the prior distribution as the uniform distribution over $[60,140]\times[15,35]\times[100,200]\times[3,10]\times[60,140]\times[15,35]\times[100,200]\times[3,10]\times[50,100]\times[10,20]\times[100,200]\times[15,35]$. We use Gaussian process regression as a machine learning model. The hyperparameter of the regularization constant $\lambda$ is $0.1$.

4.3.2 Result

The execution time of model bridging is 1.1 [s] in the presented computational environment, while simulator calibration requires about 1.3 [h] for $L+1$-th dataset. The simulator calibration requires $m\times n$ execution of simulation, and each simulation takes 3 [s] in this case. The results of the mean and standard deviation of the estimated parameters for one-leave-out cross-validation are shown in the bottom rows in Table 2. All parameters of estimation by the model-bridging framework are accurate within the standard deviation for $\theta^{(0)}$ and $\theta^{(1)}$, respectively. We can see the effectiveness of a high-dimensional parameter space with a realistic experiment. From the estimated result of simulation parameters, we can understand that the difference in $\theta_{5}$ results in different predictions for each situation, while other parameters are constant. This insight obtained from the interpretable parameters leads to improvements in the production process.

4.4 Experiment with Simulator for Fluid Dynamics

4.4.1 Setting

Through computer-aided engineering (CAE) simulations, we confirmed that our model-bridging algorithm is applicable to the simulation of fluid-dynamics systems. We employed the typical benchmark in this field, called “cavity flow experiment,” with OpenFOAM®222https://www.openfoam.com/ 333https://www.openfoam.com/documentation/tutorial-guide (Fig. 7 (A)). We considered a two-dimensional squared space called “cavity” fulfilled with a fluid having an unknown Reynolds number. The Reynolds number is used to help predict flow patterns and velocities in fluid dynamics. Turbulent flow is somewhat challenging to predict, even though it is ubiquitous in real-world situations. In this experiment, input $X_{i}\in\mathbb{R}$ is the velocity of the material on top of the cavity; the output $Y_{i}\in\mathbb{R}$ is the velocity at the particular point (see Fig. 7 (A)); and parameter $\theta\in\mathbb{R}$ is the Reynolds number (see supplementary material for details). The number of data $n=50$, the number of samples $m=41$, and the number of dataset $L=41$ are generated by different true $\theta_{l}(=\theta^{\rm true}_{l})$. The prior distribution is defined as the uniform distribution over $[20000,65000]$. We used Gaussian process regression as a machine learning model. The hyperparameter of regularization is $\lambda=1.0^{-5}$.

4.4.2 Result

The execution time of model bridging is 2.6 [s] in the presented computational environment, while the simulator calibration requires about 9 [h]. Each simulation takes about 17 [s] in this case. Figure 7 (B) shows the estimated result of $\hat{\theta}^{\rm sim}$ by simulator calibration and $\hat{\theta}^{\rm MB}$ by model bridging as a function of true $\theta$ for $L=41$ dataset with one-leave-out cross-validation. The dashed line in Fig. 7 (B) shows $\theta^{\rm true}=\hat{\theta}^{\rm MB}(=\hat{\theta}^{\rm sim})$ to ensure that the estimation is accurate if the result is on the dashed line. We can see a reasonable estimation of $\hat{\theta}^{\rm MB}$. The result of the velocity prediction of velocity $Y_{i}$ is also reasonably accurate (see the supplementary material). Human experts can understand why the Reynolds number causes such flow of fluid.

5 Discussion

There are many possible options to be discussed in the proposed framework for the individual-use case. In this study, we assume the given observed dataset as the problem setting. Further, there are two other possible ways for problem setting with the assumption of the data generation process: 1) generate data from $f_{\rm ml}(x;\xi)$ and 2) generate data from $f_{\rm sim}(x;\theta)$. Considering another case with these assumptions of data generation might be meaningful, e.g., when the real observed data are limited or when the simulation has high confidence. Another option to be discussed is the parametric or non-parametric regression model for the model-bridging function $\hat{T}$. In this study, we present the practical effectiveness of the model-bridging framework, while a theoretical analysis of the asymptotic behavior of this framework is still desired.

6 Conclusion

We propose a novel framework named “model bridging” to bridge from the un-interpretable machine learning model to the simulation model with interpretable parameters. The model-bridging framework enables us to obtain precise predictions from the machine learning model as well as obtain the interpretable simulation parameter simultaneously without the expensive calculations of a simulation. We confirmed the effectiveness of the model-bridging framework and accuracy of the estimated simulation parameter using production simulation and simulation of fluid dynamics, which are widely used in the real-world manufacturing industry.

A. Explicit Formulation of Norm of Empirical Kernel Mean

We present the explicit formulation of $\|\hat{\mu}^{\rm MB}_{L+1}-\hat{\mu}^{\rm sim}_{L+1}\|^{2}_{\mathcal{H}}$. The key for the calculation is the relation $\langle\hat{\mu},\hat{\mu}^{\prime}\rangle=\sum\sum k_{\theta}(\theta,\theta^{\prime})$ for $\theta$ kernel. The RKHS norm between estimated $\hat{\mu}^{\rm MB}_{L+1}$ and target $\hat{\mu}^{\rm sim}_{L+1}$ is described as follows:

[TABLE]

Then, the norm of the empirical kernel mean is obtained.

B. Detailed Setting of Experiment for Simple Production Simulation

A production simulation is a general-purpose simulation software package for discrete and interconnected systems, which is used to model various processes such as production, logistics, transportation, and office work. All processes are implemented in WITNESS, which is a general-purpose simulation software package for discrete and interconnected systems. The purpose of the production simulation in this experiment is to predict the total production time when the number of products to be manufactured is set. Figure 4 shows a typical assembly process for one product with four parts used in this experiment. The product consists of a “TOPS” part, a “BOTTOMS” part, and two “SCREWS.” The products assembled in the “ASSEMBLY” machine are inspected by the “INSPECTION” machine before shipping. The “INSPECTION” machine starts when four assembled products arrive, and it can inspect four assembled products simultaneously. The parameters $\theta_{1}$ and $\theta_{2}$ represent the mean and standard deviation in a normal distribution of elapsed time for the “ASSEMBLY” machine, respectively. The parameters $\theta_{3}$ and $\theta_{4}$ represent the mean and standard deviation in a normal distribution of elapsed time for the “INSPECTION” machine, respectively.

C. Details of the Experiment with Realistic Production Simulator

As a realistic experimental setting for a factory manufacturing valves, the process details are described below. All processes are implemented in WITNESS, which is a general-purpose simulation software package for discrete and interconnection systems. Figure 8 is an illustration of the simulation model for a realistic experiment. We defined the total production time as $Y_{i}\in\mathbb{R}$ when the number of units $X_{i}\in\mathbb{R}^{3}$ for three types of products to be manufactured. Each type of product has different elapsed times at the “SAW” process. Figure 9 shows the two datasets and regression results by BNN as representatives for all datasets. The discrete data is originated from a series of batch processing in the simulation model.

**Cutting process: ** The first phase of the manufacturing process begins with the arrival of a pipe having the same diameter and length, 30 cm. The pipes arrive at a fixed time interval based on the vendor’s supply schedule. Subsequently, each pipe is cut to 10-cm sections along the length. Thus, three parts can be obtained from one pipe. For the cutting process, a worker who performs changeover, repair, and disconnection is assigned. The worker goes for a lunch break once every eight hours. Thereafter, the parts are transferred from the cutting process to the coating process on a conveyor belt.

**Coating process: ** The cut parts are coated for protection. In the coating machine, six parts are batch-processed at once. The coating material must be prepared in the coating machine prior to the part. Otherwise, the parts will be degraded by the heat. When the parts ride on the belt conveyor, the sensor detects them, and the coating material is prepared.

Inspection process: After being coated, each part is placed in the inspection waiting buffer before the inspection step. The inspector will remove the parts individually from the waiting buffer and inspect the coating quality. If the part fails the quality inspection, the inspector places the part in the recoating waiting buffer. The coating machine must process the parts of the recoating buffer preferentially. When the part passes the quality inspection, the inspector sends the part to the curing step.

**Harden process: ** In the harden (quenching) process, up to 10 parts are processed simultaneously on a first-come first-out basis, and each part is quenched for at least one hour.

**Grind process: ** The quenched parts are polished for satisfying the customer’s specifications. Two polishing machines with the same priority are available. Each machine uses special jigs to process four parts simultaneously. Each of the two polishing machines produces two different types of valves. Further, 10 jigs exist in the system, and when not in use, they are placed in the jig storage buffer. A loader fixes the four parts with a jig and sends it to the polishing machine. The polishing machine sends the jig and four parts to the unloader after the polishing is finished. The unloader sends the finished parts to the valve storage area and the jig to the jig return area. The two types of valves are separated and placed in a dedicated valve storage buffer. As the jig needs to be used again, it is returned from the jig return conveyor to the jig storage buffer.

**Cleaning process: ** The valves removed from the valve storage area are cleaned before shipment. In the washing machine, five stations are available where the valves can be placed one at a time, and the valves are cleaned in these stations. Up to 10 valves of each type can be washed simultaneously. When the valve type is changed, the cleaning head must be replaced.

D. Details of the Experiment with Fluid-Dynamics Simulator

These experiments are performed using the general-purpose open CAE simulators using the Finite Element Method (FEM) solver; OpenFOAM®. OpenFOAM ® includes some realistic problems as a tutorial. The example of “cavity flow” is one of the typical benchmarks. Figure 10 shows illustrations of the experiment of “cavity flow” that are simulated in this study. Most of the initial experimental settings, such as the number of FEM meshes, cavity size, and so on, were the same as in the tutorial. The range of the Reynolds number (Re) used was $10000<{\rm Re}<70000$. These settings cause the typical turbulent flow, as shown in Fig. 10.

Time-series data are observed where the center of flow is moving, as shown in Fig. 10 (a) to (d). Part (a) of Fig. 10 shows the state at $t=0$ [s] (initial state), (b) at $t=1.0$ [s], (c) at $t=2.5$ [s], and (d) at $t=5.5$ [s]. In this experiment, $X_{i}\in\mathbb{R}$ is velocity of the material on top of the cavity; $Y_{i}\in\mathbb{R}$ is the velocity at the particular x-marked point at $t=5.5$ [s] (Fig. 10 (d)); and parameter $\theta\in\mathbb{R}$ is the Reynolds number. Figure 11 shows the relation between $X_{i}$ and $Y_{i}$ for five different dataset as representatives. It is difficult to select the parametric statistical model as a regression function owing to the non-trivial relation.